互补错误函数erfcf（）的矢量化实现

2018-06-30 16:29:28

互补误差函数erfc是与标准正态分布密切相关的特殊函数。它经常用于统计和自然科学（例如扩散问题），因为这种分布的“尾巴”需要考虑，因此错误函数erf的使用是不合适的。

互补错误函数在ISO C99标准数学库中作为函数erfcf ， erfc和erfcl ; 这些随后也被采用到ISO C ++中。因此源代码很容易在该库的开源实现中找到，例如在glibc中。

然而，许多现有的实现本质上是标量的，而现代处理器硬件是面向SIMD的（明确地，如在x86 CPU中，或者隐式地，如在GPU中）。出于性能原因，可矢量化的实现因此是非常需要的。这意味着需要避免分支，除非作为选择分配的一部分。同样，没有指出表的广泛使用，因为并行查找往往是低效的。

如何构建单精度函数erfcf()的高效矢量化实现？以ulp衡量的准确度应该与glibc的标量实现大致相同，最大误差为3.12575 ulps（通过穷举测试确定）。由于所有主要处理器架构（CPU和GPU）目前均提供融合乘法 - 加法（FMA），因此可以假设融合乘法 - 加法（FMA）的可用性。虽然浮点状态标志和errno可以忽略，但是应该根据ISO C的IEEE 754绑定来处理非规范，无穷和NaN。

在研究了各种方法之后，似乎最适合的方法是以下文章中提出的算法：

MM Shepherd和JG Laframboise，“（1 + 2 x）exp（x2）erfc x在0≤x <∞时的Chebyshev逼近”。计算数学，第36卷，第153号，1981年1月，第249-253页（在线副本）

本文的基本思想是创建（1 + 2 x）exp（x2）erfc（x）的近似值，从中可以通过简单除以（1 + 2 x）来计算erfcx（x），并且erfc x）然后乘以exp（-x2）。函数紧密有界的范围，函数值大致在[1,1.3]，而且它的一般“平坦性”很好地适用于多项式逼近。这种方法的数值特性通过缩小近似间隔进一步改善：原始参数x通过q =（x-K）/（x + K）变换，其中K是适当选择的常数，接着计算p（q），其中p是一个多项式。

由于erfc -x = 2-erfc x，我们只需要考虑通过这个变换映射到区间[-1，1]的区间[0，∞]。对于IEEE-754单精度，对于x> 10.0546875， erfcf()消失（变为零），因此只需考虑x∈[0，10.0546875）。在这个范围内，K的“最优”值是多少？我知道没有数学分析可以提供答案，本文根据实验提出K = 3.75。

人们可以容易地确定，对于单精度计算，9次方的最小多项式近似对于在该总体邻域中的各种K值是足够的。对于K = {2,2.625,3.3125}，观察到用Remez算法系统地产生这样的近似值，K以1/16的步长变化在1.5和4之间，最低的近似误差。其中，K = 2是最有利的选择，因为它适用于非常精确的（x - K）/（x + K）计算，如此问题所示。

值K = 2和x的输入域表明有必要使用我的答案中的变体4，但是一旦能够以实验方式证明较便宜的变体5在这里达到相同的精度，这可能是由于非常浅近似函数的斜率为q> -0.5，这会导致参数q中的任何误差减少大约十倍。

由于erfc()计算除了初始逼近之外还需要后处理步骤，很显然，为了获得足够准确的最终结果，这两种计算的准确性必须很高。必须使用纠错技术。

我们观察到（1 + 2 x）exp（x2）erfc（x）的多项式近似中的最高有效系数是（1 + s）的形式，其中s <0.5。这意味着我们可以通过分解1来更精确地表示前导系数，并且仅在多项式中使用s。因此，不是计算多项式p（q），然后乘以倒数r = 1 /（1 + 2 x），它在数学上是等价的，但在数值上有利于将核心近似计算为p（q）+1，并且使用p计算fma (p, r, r) 。

通过从倒数r计算初始商q，通过FMA计算残差e = p + 1-q *（1 + 2 x），然后使用e应用修正q = q +（e * r），再次使用FMA。

指数运算具有误差放大特性，因此必须谨慎执行e-x2的计算。 FMA的可用性允许将-x2计算为double- float shigh：slow。 ex是它自己的衍生物，所以可以计算eshigh：慢如eshigh + eshigh * slow。该计算可以与前一个中间结果r的乘积相结合，得到r = r * eshigh + r * eshigh * slow。通过使用FMA，可以确保尽可能精确地计算最重要的术语r * eshigh。

将上述步骤与几个简单的选择结合起来处理异常情况和负面的争论，可以得到以下C代码：

float my_expf (float);

/*  
 * Based on: M. M. Shepherd and J. G. Laframboise, "Chebyshev Approximation of 
 * (1+2x)exp(x^2)erfc x in 0 <= x < INF", Mathematics of Computation, Vol. 36,
 * No. 153, January 1981, pp. 249-253.  
 */  
float my_erfcf (float x)
{
    float a, d, e, m, p, q, r, s, t;

    a = fabsf (x); 

    /* Compute q = (a-2)/(a+2) accurately. [0, 10.0546875] -> [-1, 0.66818] */
    m = a - 2.0f;
    p = a + 2.0f;
    r = 1.0f / p;
    q = m * r;
    t = fmaf (q + 1.0f, -2.0f, a); 
    e = fmaf (q, -a, t); 
    q = fmaf (r, e, q); 

    /* Approximate (1+2*a)*exp(a*a)*erfc(a) as p(q)+1 for q in [-1, 0.66818] */
    p =             -0x1.a48024p-12f;  // -4.01020574e-4
    p = fmaf (p, q, -0x1.42a172p-10f); // -1.23073824e-3
    p = fmaf (p, q,  0x1.585784p-10f); //  1.31355994e-3
    p = fmaf (p, q,  0x1.1ade24p-07f); //  8.63243826e-3
    p = fmaf (p, q, -0x1.081b72p-07f); // -8.05991236e-3
    p = fmaf (p, q, -0x1.bc0b94p-05f); // -5.42047396e-2
    p = fmaf (p, q,  0x1.4ffc40p-03f); //  1.64055347e-1
    p = fmaf (p, q, -0x1.540840p-03f); // -1.66031361e-1
    p = fmaf (p, q, -0x1.7bf612p-04f); // -9.27639678e-2
    p = fmaf (p, q,  0x1.1ba03ap-02f); //  2.76978403e-1

    /* Divide (1+p) by (1+2*a) ==> exp(a*a)*erfc(a) */
    t = a + a;
    d = t + 1.0f;
    r = 1.0f / d;
    q = fmaf (p, r, r); // q = (p+1)/(1+2*a)
    e = (p - q) + fmaf (q, -t, 1.0f); // (p+1) - q*(1+2*a)
    r = fmaf (e, r, q);

    /* Multiply by exp(-a*a) ==> erfc(a) */
    s = a * a; 
    e = my_expf (-s);  
    t = fmaf (a, -a, s);
    r = fmaf (r, e, r * e * t);

    /* Handle NaN arguments to erfc() */
    if (!(a <= 0x1.fffffep127f)) r = x + x;

    /* Clamp result for large arguments */
    if (a > 10.0546875f) r = 0.0f;

    /* Handle negative arguments to erfc() */
    if (x < 0.0f) r = 2.0f - r; 

    return r;
}

/* Compute exponential base e. Maximum ulp error = 0.87161 */
float my_expf (float a)
{
    float c, f, r;
    int i;

    // exp(a) = exp(i + f); i = rint (a / log(2)) 
    c = 0x1.800000p+23f; // 1.25829120e+7
    r = fmaf (0x1.715476p+0f, a, c) - c; // 1.44269502e+0
    f = fmaf (r, -0x1.62e400p-01f, a); // -6.93145752e-1 // log_2_hi 
    f = fmaf (r, -0x1.7f7d1cp-20f, f); // -1.42860677e-6 // log_2_lo
    i = (int)r;
    // approximate r = exp(f) on interval [-log(2)/2,+log(2)/2]
    r =             0x1.6a98dap-10f;  // 1.38319808e-3
    r = fmaf (r, f, 0x1.1272cap-07f); // 8.37550033e-3
    r = fmaf (r, f, 0x1.555a20p-05f); // 4.16689515e-2
    r = fmaf (r, f, 0x1.55542ep-03f); // 1.66664466e-1
    r = fmaf (r, f, 0x1.fffff6p-02f); // 4.99999851e-1
    r = fmaf (r, f, 0x1.000000p+00f); // 1.00000000e+0
    r = fmaf (r, f, 0x1.000000p+00f); // 1.00000000e+0
    // exp(a) = 2**i * exp(f);
    r = ldexpf (r, i);
    // handle special cases
    if (!(fabsf (a) < 104.0f)) {
        r = a + a; // handle NaNs
        if (a < 0.0f) r = 0.0f;
        if (a > 0.0f) r = 1e38f * 1e38f; // + INF
    }
    return r;
}

我在上面的代码中使用了我自己的expf()实现，以将我的工作与不同计算平台上的expf()实现中的差异分开。但是，最大错误接近0.5 ulp的expf()的任何实现都应该可以正常工作。如上所示，即使用my_expf() ， my_erfcf()的最大错误为2.65712 ulps。提供了一个可矢量化的expf()可用性，上面的代码应该没有问题地矢量化。我用英特尔编译器13.1.3.198做了一个快速检查。我在循环中调用了my_erfcf() ，添加了#include <mathimf.h> ，并调用了my_expf() ，并调用了expf() ，然后使用这些命令行开关进行编译：

/Qstd=c99 /O3 /QxCORE-AVX2 /fp:precise /Qfma /Qimf-precision:high:expf /Qvec_report=2

英特尔编译器报告说，该循环已被矢量化，我通过检查反汇编的二进制代码进行了双重检查。

由于my_erfcf()只使用倒数而不是全分度，所以它可以使用快速倒数实现，只要它们提供了几乎正确的舍入结果。对于在硬件中提供快速单精度互易近似的处理器，可以通过将其与具有三次收敛的哈雷迭代耦合来轻松实现。这种x86处理器方法的（标量）例子是：

/* Compute 1.0f / a almost correctly rounded. Halley iteration with cubic convergence */
float fast_recipf (float a)
{
    __m128 t;
    float e, r;
    t = _mm_set_ss (a);
    t = _mm_rcp_ss (t);
    _mm_store_ss (&r, t);
    e = fmaf (r, -a, 1.0f);
    e = fmaf (e, e, e);
    r = fmaf (e, r, r);
    return r;
}

链接地址: http://www.djcxy.com/p/85605.html

上一篇: Vectorizable implementation of complementary error function erfcf()

下一篇: Can I use rounding to ensure determinism of atomic floating point operations?