nju*_*ffa 8 c algorithm math floating-point
互补误差函数erfc是与标准正态分布密切相关的特殊函数.它经常用于统计学和自然科学(例如扩散问题),其中需要考虑该分布的"尾部",因此使用误差函数erf是不合适的.
互补误差函数是在ISO C99标准数学库提供作为功能erfcf,erfc和erfcl; 这些随后也被采用到ISO C++中.因此,源代码可以很容易地在该库的开源实现中找到,例如在glibc中.
然而,许多现有的实现本质上是标量的,而现代处理器硬件是面向SIMD的(显式地,如在x86 CPU中,或隐式地,如在GPU中).出于性能原因,因此非常需要可矢量化的实现.这意味着需要避免分支,除非作为选择分配的一部分.同样,没有指出表的广泛使用,因为并行查找通常是低效的.
如何构建单精度函数的高效矢量化实现erfcf()？测量的准确度ulp应与glibc的标量实现大致相同,其最大误差为3.12575 ulps(通过详尽测试确定).可以假设融合乘法加法(FMA)的可用性,因为此时所有主要处理器架构(CPU和GPU)都提供它.处理浮点状态标志时errno可以忽略,非正规,无穷大和NaN应根据ISO 75的IEEE 754绑定进行处理.
在研究了各种方法之后，最合适的方法是以下论文中提出的算法：
\nMM Shepherd 和 JG Laframboise，“0 \xe2\x89\xa4 x < \xe2\x88\x9e中(1 + 2 x ) exp( x 2 ) erfc x的切比雪夫近似。” 《计算数学》，第 36 卷，第 153 期，1981 年 1 月，第 249-253 页（在线副本）
\n本文的基本思想是创建 (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是多项式。
\n由于 erfc - x = 2 - erfc x，我们只需要考虑区间 [0, \xe2\x88\x9e] 通过此变换映射到区间 [-1, 1] 。对于 IEEE-754 单精度，当xerfcf() > 10.0546875时消失（变为零），因此只需考虑x \xe2\x88\x88 [0, 10.0546875)。对于这个范围，K 的“最佳”值是多少？我知道没有任何数学分析可以提供答案，论文根据实验建议 K = 3.75。
人们可以很容易地确定，对于单精度计算，9 次极小极大多项式近似足以满足该一般附近的各种 K 值。使用 Remez 算法系统地生成此类近似值，K 在 1.5 和 4 之间以 1/16 的步长变化，观察到 K = {2, 2.625, 3.3125} 时的最低近似误差。其中，K = 2 是最有利的选择，因为它有助于非常精确地计算 ( x - K) / ( x + K)，如本问题所示。
\n值 K = 2 和x的输入域表明有必要使用我的答案中的变体 4 ，但是一旦可以通过实验证明较便宜的变体 5 在这里达到相同的精度，这可能是由于非常浅q > -0.5时近似函数的斜率，这会导致参数q中的任何误差大约减少十倍。
\n由于erfc()除了初始近似之外， 的计算还需要后处理步骤，因此显然，这两种计算的精度都必须很高，才能获得足够准确的最终结果。必须使用纠错技术。
我们观察到 (1 + 2 x ) exp( x 2 ) 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。
\n求幂具有误差放大特性，因此e - x 2的计算必须小心执行。FMA的可用性允许将 - x 2计算为双float高：低。e x是它自己的导数，因此可以将 e s high :s low计算为 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 。
将上述步骤与一些处理异常情况和负参数的简单选择相结合，得到以下 C 代码：
\nfloat my_expf (float);\n\n/* Compute complementary error function.\n *\n * Based on: M. M. Shepherd and J. G. Laframboise, "Chebyshev Approximation of \n * (1+2x)exp(x^2)erfc x in 0 <= x < INF", Mathematics of Computation, Vol. 36,\n * No. 153, January 1981, pp. 249-253.\n *\n * maximum error: 2.65184 ulps\n */ \nfloat my_erfcf (float x)\n{\n float a, d, e, p, q, r, s, t;\n\n a = fabsf (x);\n \n /* Compute q = (a-2)/(a+2) accurately. [0, 10.0546875] -> [-1, 0.66818] */\n p = a + 2.0f;\n r = 1.0f / p;\n q = fmaf (-4.0f, r, 1.0f);\n t = fmaf (q + 1.0f, -2.0f, a); \n e = fmaf (-a, q, t); \n q = fmaf (r, e, q); \n \n /* Approximate (1+2*a)*exp(a*a)*erfc(a) as p(q)+1 for q in [-1, 0.66818] */\n p = -0x1.a4a000p-12f; // -4.01139259e-4\n p = fmaf (p, q, -0x1.42a260p-10f); // -1.23075210e-3\n p = fmaf (p, q, 0x1.585714p-10f); // 1.31355342e-3\n p = fmaf (p, q, 0x1.1adcc4p-07f); // 8.63227434e-3\n p = fmaf (p, q, -0x1.081b82p-07f); // -8.05991981e-3\n p = fmaf (p, q, -0x1.bc0b6ap-05f); // -5.42046614e-2\n p = fmaf (p, q, 0x1.4ffc46p-03f); // 1.64055392e-1\n p = fmaf (p, q, -0x1.540840p-03f); // -1.66031361e-1\n p = fmaf (p, q, -0x1.7bf616p-04f); // -9.27639827e-2\n p = fmaf (p, q, 0x1.1ba03ap-02f); // 2.76978403e-1\n\n /* Divide (1+p) by (1+2*a) ==> exp(a*a)*erfc(a) */\n d = fmaf (2.0f, a, 1.0f);\n r = 1.0f / d;\n q = fmaf (p, r, r); // q = (p+1)/(1+2*a)\n e = fmaf (fmaf (q, -a, 0.5f), 2.0f, p - q); // residual: (p+1)-q*(1+2*a)\n r = fmaf (e, r, q);\n\n /* Multiply by exp(-a*a) ==> erfc(a) */\n s = a * a; \n e = my_expf (-s);\n t = fmaf (-a, a, s);\n r = fmaf (r, e, r * e * t);\n\n /* Handle NaN, Inf arguments to erfc() */\n if (!(a < INFINITY)) r = x + x;\n\n /* Clamp result for large arguments */\n if (a > 10.0546875f) r = 0.0f;\n \n /* Handle negative arguments to erfc() */\n if (x < 0.0f) r = 2.0f - r; \n\n return r;\n}\n\n/* Compute exponential base e. Maximum ulp error = 0.86565 */\nfloat my_expf (float a)\n{\n float c, f, r;\n int i;\n\n // exp(a) = exp(i + f); i = rint (a / log(2))\n c = 0x1.800000p+23f; // 1.25829120e+7\n r = fmaf (0x1.715476p+0f, a, c) - c; // 1.44269502e+0\n f = fmaf (r, -0x1.62e400p-01f, a); // -6.93145752e-1 // log_2_hi\n f = fmaf (r, -0x1.7f7d1cp-20f, f); // -1.42860677e-6 // log_2_lo\n i = (int)r;\n // approximate r = exp(f) on interval [-log(2)/2,+log(2)/2]\n r = 0x1.694000p-10f; // 1.37805939e-3\n r = fmaf (r, f, 0x1.125edcp-07f); // 8.37312452e-3\n r = fmaf (r, f, 0x1.555b5ap-05f); // 4.16695364e-2\n r = fmaf (r, f, 0x1.555450p-03f); // 1.66664720e-1\n r = fmaf (r, f, 0x1.fffff6p-02f); // 4.99999851e-1\n r = fmaf (r, f, 0x1.000000p+00f); // 1.00000000e+0\n r = fmaf (r, f, 0x1.000000p+00f); // 1.00000000e+0\n // exp(a) = 2**i * exp(f);\n r = ldexpf (r, i);\n if (!(fabsf (a) < 104.0f)) {\n r = a + a; // handle NaNs\n if (a < 0.0f) r = 0.0f;\n if (a > 0.0f) r = INFINITY;\n }\n return r;\n}\nexpf()我在上面的代码中使用了我自己的实现，以将我的工作与expf()不同计算平台上的实现差异隔离开来。但任何expf()最大误差接近 0.5 ulp 的实现都应该可以正常工作。如上图，即使用时my_expf()，my_erfcf()最大误差为 2.65184 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\n英特尔编译器报告该循环已矢量化，我通过检查反汇编的二进制代码进行了双重检查。
\n由于my_erfcf()仅使用倒数而不是完全除法，因此可以使用快速倒数实现，只要它们提供几乎正确舍入的结果即可。对于在硬件中提供快速单精度倒数近似的处理器，可以通过将其与三次收敛的 Halley 迭代相结合来轻松实现。针对 x86 处理器的此方法的（标量）示例是：
/Qstd=c99 /O3 /QxCORE-AVX2 /fp:precise /Qfma /Qimf-precision:high:expf /Qvec_report=2\n