错误函数erff（）的圆角实现
误差函数与标准正态分布密切相关，在自然科学以及其他领域中经常出现。 例如，它在定价选项中用于财务。 虽然对它的支持首先添加到ISO C99，然后以函数erf() ， erff()的形式添加到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()上是由硬件支持的，那么如何构造一个单精度误差函数erff() ，这个函数既忠实又有效？ 优选地，代码应该是可矢量化的，可能在稍微修改代码之后。
看看误差函数的图表，可以看到函数关于原点是对称的; 因此近似值可以平凡地限制在正半平面上。 此外，该图可以分成两个部分，这些部分之间的边界在x = 1附近。在更接近原点的部分，误差函数是相当线性的并且是垂直支配的，而在离原点更远的部分中以指数衰减的方式在水平方向占主导地位并逐渐接近统一。
一个合理的结论是，一个简单的多项式近似x * p（x）适用于接近零的区段，而另一个区段可以很好地近似于1 - exp（x * q（x）），其中q（x）是第二个多项式近似。 基于误差函数的泰勒级数展开式，靠近原点的线段的近似值实际上应该是x * p（x2）的形式。
接下来的第一项任务就是找到两段之间的切换点。 我使用了一个实验方法，从0.875的切换点开始，并逐渐向1.0迈进。 对于每个切换点的值，我用Remez算法生成零和切换点之间的误差函数的初始最小最大值近似值。 然后在准确性方面进一步改进，基于通过使用FMA操作的Horner方案对多项式的评估，对系数值使用搜索启发式。 重复增加切换点，直到产生的近似值的最大误差超过1 ulp。 通过这个过程，我确定了两个近似段之间的最佳边界为0.921875。 这导致近似值x * p（x2）达到最小误差几乎不到1 ulp。
Remez算法也被用来为多项式q（x）提供一个初始极大极小计算。 早期很清楚，q（x）和exp（）内部使用的近似值之间存在相当数量的相互作用，它们的错误相互补充或补偿。 这意味着q（x）的系数的最佳选择将与exp()的实现紧密相关，并且必须作为初始系数集的启发式细化的一部分来考虑。 因此，我决定使用自己的expf()实现，以便将自己与任何特定的库实现隔离开来。 作为最低要求， expf()本身需要忠实地四舍五入，并且可能必须符合这种方法工作的稍微更严格的误差，尽管我没有试图确定究竟是多么紧密。 在这种情况下，我自己实现的expf()提供了0.87161 ulps的错误限制，这足以证明这一点。
由于需要使用exp（）的段是慢路径，因此我选择使用Estrin的方案来计算多项式q（x）的低阶项，以便提高指令级并行性和性能。 这样做对准确性的影响可以忽略不计。 出于准确性的原因，Horner方案必须用于多项式的高阶项。 查看两个多项式的最小有效系数，可以发现两者都是1.128 ...，因此我们可以通过将系数分解为（1 + 0.128 ...）来稍微提高精度，这有助于使用FMA来执行最后乘以x。
最后，我能够实现erff()的实现，其中两个代码路径中的每一个通过针对较高精度参考的穷举测试而实现最大误差仅略低于1 ulp。 该功能因此被忠实地四舍五入。 使用FMA是这一成功的关键组成部分。 根据工具链的不同，下面显示的C99代码可以按原样进行矢量化处理，或者可以手动对其进行修改，以便两个代码路径可以在最后选择所需结果的同时进行计算。 高性能数学库包含一个应该用来代替我的自定义函数my_expf()的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;
}