高效忠实地实现错误函数erff（）
误差函数与标准正态分布密切相关，并且经常出现在自然科学以及其他领域。 例如，它在定价选项时用于财务。 虽然它首先添加到ISO C99，然后以函数erf() ， erff()的forms添加到C ++，但直到最近才出现了至少一个流行的C / C +工具链。 许多项目仍然使用他们自己的错误函数实现，通常基于旧文献的近似，例如Abramowitz和Stegun ，后者又回到
Cecil Hastings Jr，“数字计算机的近似值”。 普林斯顿大学出版社，1955年
在现代计算中，超越函数的忠实圆形实现通常被视为数学库的最低准确度标准; 这样的标准仍然允许高性能实现。 当函数返回最大误差小于1 ulp的结果与整个输入域中的数学值相比时，函数被称为忠实舍入。 当使用IEEE-754单精度操作实现时，较早发布的算法不能提供忠实的圆形结果。
现代计算机硬件提供称为融合乘法 – 加法（或简称FMA）的浮点运算，它计算浮点乘法，然后进行相关的浮点加法，以便在加法中使用完整的未接地产品，并且只有单个舍入发生在操作结束时。 IBM于1990年推出的这种融合操作在许多计算中提供了更高的准确性和更高的性能。 它可用于当今最流行的两种CPU架构（ARM和x86）以及GPU。 它已通过fmaf()和fmaf()函数在C和C ++中fmaf() 。
假设FMA本身是由硬件支持的，那么如何构建单精度错误函数erff() ，它既忠实圆润又高效？ 优选地，代码应该是可矢量化的，可能在次要代码修改之后。
观察误差函数的图形 ，可以观察到函数关于原点是对称的; 因此，近似可以简单地限于正半平面。 此外，图形可以分成两个部分，这些部分之间的边界接近x = 1.在靠近原点的部分中，误差函数是相当线性的并且垂直占优势，而在远离原点的部分中它是水平支配的，渐近地以指数衰减的方式接近统一。
一个合理的结论是简单的多项式逼近x * p（x）适用于接近零的段，而另一个段很好地近似为1 – exp（x * q（x）），其中q（x）是第二个多项式近似。 基于误差函数的泰勒级数展开，原点附近的段的近似实际上应该是x * p（x 2 ）的forms。
然后，第一个任务是找到两个段之间的切换点。 我使用了一种实验方法，从切换点0.875开始，逐步向1.0行进。 对于切换点的每个值，我使用Remez算法生成零和切换点之间的误差函数的初始极小极大近似值。 然后，在准确性方面，使用搜索试探法对系数值进行进一步改进，基于多项式的评估，通过使用FMA操作的Horner方案。 重复增加切换点，直到得到的近似值的最大误差超过1 ulp。 通过这个过程，我确定两个近似段之间的最佳边界为0.921875。 这导致近似x * p（x 2 ），其实现的最大误差仅小于1 ulp。
Remez算法还用于为多项式q（x）提供初始极小极大值计算。 早期很清楚，q（x）和内部使用的近似值（）在他们的误差相互加强或相互补偿的方式之间存在相当大的相互作用。 这意味着q（x）系数的最佳选择将与exp()的实现紧密相关，并且必须作为初始系数集的启发式细化的一部分加以考虑。 因此，我决定使用自己的expf()实现，以便将自己与任何特定的库实现隔离开来。 作为最低要求， expf()本身需要忠实地舍入，并且可能必须符合这种方法工作的稍微更严格的错误限制，尽管我没有尝试确切地确定有多紧。 在这种情况下，我自己的expf()提供了0.87161 ulps的错误限制，结果certificate是足够的。
由于需要使用exp（）的段是慢路径，因此我选择将Estrin方案用于多项式q（x）的低阶项，以便提高指令级并行性，从而提高性能。 这样做对准确性的影响可以忽略不计。 出于准确性原因，必须将Horner方案用于多项式的高阶项。 观察两个多项式的最不重要系数，可以观察到两者都是1.128 …，因此我们可以通过将系数分成（1 + 0.128 …）来略微提高精度，这有助于使用FMA执行与x的最终乘法。
最后，我能够实现erff() ，其中两个代码路径中的每一个都实现了不到1 ulp的最大误差，这是通过针对更高精度参考的详尽测试建立的。 因此，该function忠实地圆润。 使用FMA是这一成功的关键组成部分。 根据工具链，下面显示的C99代码可以按原样进行矢量化，或者可以手动修改它，以便两个代码路径同时与最终选择的所需结果一起计算。 高性能数学库包括expf()可向量化版本，应该使用它来代替我的自定义函数my_expf() 。 并非所有expf()矢量化实现都提供了足够的精度，而对于其他实现，多项式q（x）中的系数的调整将是必要的。
如果使用自定义版本的expf() ，就像我在这里所做的那样，出于性能原因，人们可能希望用更快的机器特定代码替换对ldexpf()的调用。
/* 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; } /* compute error function. max ulp error = 0.99993 */ float my_erff (float a) { float r, s, t, u; t = fabsf (a); s = a * a; if (t > 0.921875f) { // 0.99527 ulp r = fmaf (0x1.222900p-16f, t, -0x1.91d2ccp-12f); // 1.72948930e-5, -3.83208680e-4 u = fmaf (0x1.fd1336p-09f, t, -0x1.8d6300p-06f); // 3.88393435e-3, -2.42545605e-2 r = fmaf (r, s, u); r = fmaf (r, t, 0x1.b55cb0p-4f); // 1.06777847e-1 r = fmaf (r, t, 0x1.450aa0p-1f); // 6.34846687e-1 r = fmaf (r, t, 0x1.079d0cp-3f); // 1.28717512e-1 r = fmaf (r, t, t); r = my_expf (-r); r = 1.0f - r; r = copysignf (r, a); } else { // 0.99993 ulp r = -0x1.3a1a82p-11f; // -5.99104969e-4 r = fmaf (r, s, 0x1.473f48p-08f); // 4.99339588e-3 r = fmaf (r, s, -0x1.b68bd2p-06f); // -2.67667342e-2 r = fmaf (r, s, 0x1.ce1a46p-04f); // 1.12818025e-1 r = fmaf (r, s, -0x1.8126e0p-02f); // -3.76124859e-1 r = fmaf (r, s, 0x1.06eba6p-03f); // 1.28379151e-1 r = fmaf (r, a, a); } return r; }