互补错误函数erfcf（）的可矢量化实现
互补误差函数erfc是与标准正态分布密切相关的特殊函数。 它经常用于统计学和自然科学（例如扩散问题），其中需要考虑这种分布的“尾部”，因此使用误差函数erf是不合适的。
补充误差函数在ISO C99标准数学库中可用作函数erfcf ， erfc和erfcl ; 这些随后也被采用到ISO C ++中。 因此，源代码可以很容易地在该库的开源实现中找到，例如在glibc中 。
然而，许多现有的实现本质上是标量的，而现代处理器硬件是面向SIMD的（显式地，如在x86 CPU中，或隐含地，如在GPU中）。 出于性能原因，非常需要可矢量化的实现。 这意味着需要避免分支，除非作为选择分配的一部分。 同样，未指示广泛使用表，因为并行查找通常是低效的。
如何构建单精度函数erfcf()的高效矢量化实现？ 以ulp为单位测量的准确度应与glibc的标量实现大致相同，其最大误差为3.12575 ulps（通过详尽测试确定）。 可以假设融合乘法加法（FMA）的可用性，因为此时所有主要的处理器架构（CPU和GPU）都提供它。 虽然可以忽略浮点状态标志和errno处理，但应根据ISO C的IEEE 754绑定处理非正规数，无穷大和NaN。
在研究了各种方法后，最合适的方法是下面提出的算法：
MM Shepherd和JG Laframboise，“Chebyshev逼近（1 + 2 x ）exp（ x 2 ）erfc x在0≤x <∞。” 计算数学 ，第36卷，第153期，1981年1月，第249-253页（ 在线副本 ）
本文的基本思想是创建（1 + 2 x ）exp（ x 2 ）erfc（ x ）的近似值，我们可以通过简单地除以（1 + 2 x ）来计算erfcx（ x ），并且erfc （x）然后乘以exp（ – x 2 ）。 函数的紧密有界范围，函数值大致在[1,1.3]，其一般“平坦度”适用于多项式近似。 通过缩小近似区间进一步改善该方法的数值性质：原始参数x由q =（ x -K）/（ x + K）变换，其中K是适当选择的常数，然后计算p （ q ） ，其中p是多项式。
由于erfc – x = 2 – erfc x ，我们只需要考虑通过该变换映射到区间[-1,1]的区间[0，∞]。 对于IEEE-754单精度， erfcf()在x > 10.0546875时消失（变为零），因此需要仅考虑x∈[0,10.0546875]。 对于这个范围，K的“最佳”值是多少？我知道没有数学分析可以提供答案，本文根据实验建议K = 3.75。
可以容易地确定，对于单精度计算，9度的极小极大多项式近似对于该一般附近的K的各种值是足够的。 使用Remez算法系统地生成这样的近似，其中K以1/16的步长在1.5和4之间变化，对于K = {2,2.625,3.3125}观察到最低近似误差。 其中，K = 2是最有利的选择，因为它有助于非常准确地计算（ x -K）/（ x + K），如该问题所示。
值K = 2和x的输入域表明有必要使用我的答案中的变体4，但是一旦可以通过实validation明较便宜的变体5在这里达到相同的准确度，这可能是由于非常浅q > -0.5的近似函数的斜率，这会导致参数q中的任何误差减少大约十倍。
由于erfc()计算除了初始近似之外还需要后处理步骤，因此很明显这两种计算的精度必须很高才能获得足够精确的最终结果。 必须使用纠错技术。
人们观察到（1 + 2 x ）exp（ x 2 ）erfc（ x ）的多项式近似中的最重要系数是（1 + s）forms，其中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 – x 2的计算。 FMA的可用性通常允许计算-x 2作为双float 高 ：s 低 。 e x是它自己的导数，因此可以计算e s high ：s s low as e s high + e s high * s low 。 该计算可以与先前中间结果r的乘法组合以产生r = r * e s high + r * e s high * s low 。 通过使用FMA，可以确保尽可能准确地计算最高有效项r * e s high 。
将上面的步骤与一些简单的选择结合起来处理exception情况和负面参数，可以得到以下C代码：
float my_expf (float); /* * Based on: MM Shepherd and JG 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 ，用调用expf()替换调用my_expf() expf() ，然后使用这些命令行开关编译：
/Qstd=c99 /O3 /QxCORE-AVX2 /fp:precise /Qfma /Qimf-precision:high:expf /Qvec_report=2
英特尔编译器报告说该循环已被矢量化，我通过检查反汇编的二进制代码进行了双重检查。
由于my_erfcf()仅使用倒数而不是完全除法，因此只要它们提供几乎正确的舍入结果，就可以使用快速倒数实现。 对于在硬件中提供快速单精度倒数近似的处理器，可以通过将其与具有立方收敛的Halley迭代耦合来轻松实现。 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; }