c - 错误函数 erff() 的高效忠实圆整实现

Question

误差函数与标准正态分布密切相关，在自然科学和其他领域中经常出现。例如，在为期权定价时，它用于金融领域。虽然对它的支持首先添加到 ISO C99，随后以函数的形式添加到 C++ erf() , erff() ，直到最近，它才出现在至少一个流行的 C/C+ 工具链中。许多项目仍然使用自己的误差函数实现，通常基于旧文献中的近似值，例如 Abramowitz and Stegun ,
反过来又回到

Cecil Hastings Jr，“数字计算机的近似值”。普林斯顿大学出版社，1955 年

在现代计算中，超越函数的忠实圆整实现通常被视为数学库精度的最低标准；这样的标准仍然允许高性能实现。当函数返回最大误差小于 1 ulp 的结果时，它被称为忠实舍入。与整个输入域的数学值进行比较。当仅使用 IEEE-754 单精度运算来实现时，较旧的已发布算法无法提供忠实舍入的结果。

现代计算机硬件提供了一种称为融合乘加(或简称 FMA)的浮点运算，它计算一个浮点乘法，然后是一个相关的浮点加法，这样在加法中使用完整的未舍入乘积，并且只有一个单次舍入发生在操作结束时。这种融合运算由 IBM 于 1990 年推出，在许多计算中提供了更高的准确性和更高的性能。它可用于当今两种最流行的 CPU 架构(ARM 和 x86)以及 GPU。它已通过 fmaf() 在 C 和 C++ 中公开。和 fmaf()职能。

假设硬件本身支持 FMA，那么如何构建单精度误差函数 erff()哪个既忠实又高效？最好是，代码应该是可向量化的，可能是在对代码进行较小的修改之后。

Answer 1

查看 graph of the error function ，可以观察到该函数关于原点对称；因此，近似值可以简单地限制在正半平面上。此外，该图可以分成两部分，它们之间的边界在 x = 1 附近。在靠近原点的部分中，误差函数是相当线性和垂直主导的，而在远离原点的部分中，它是水平支配的，并以指数衰减的方式渐近接近统一。

一个合理的结论是，简单的多项式近似 x * p (x) 适用于接近零的段，而另一段可以很好地近似为 1 - exp (x * q (x))，其中 q (x)是二次多项式近似。基于泰勒级数展开的误差
函数，原点附近线段的近似值实际上应该是 x * p (x2) 的形式。

第一个任务是找到两个段之间的切换点。为此，我使用了一种实验方法，从 0.875 处的切换点开始，并逐步将其推向 1.0。对于切换点的每个值，我使用 Remez algorithm 生成了零和切换点之间误差函数的初始极小极大近似。然后在准确性方面进一步改进，使用搜索
基于使用 FMA 操作通过霍纳方案评估多项式的​​系数值的启发式方法。重复增加切换点，直到产生的近似值的最大误差超过 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's scheme 用于多项式 q (x) 的低阶项，以提高指令级并行性和性能。这样做对准确性的影响可以忽略不计。出于精度原因，多项式的高阶项必须使用霍纳方案。查看两个多项式的最低有效系数，可以观察到两者都是 1.128...，因此我们可以通过将系数拆分为 (1 + 0.128...) 来稍微提高精度，这有助于使用 FMA 执行与 x 的最终乘法。

最后，我能够实现 erff() 的实现，其中两个代码路径中的每一个都实现了略低于 1 ulp 的最大误差，这是通过针对更高精度引用的详尽测试确定的。因此，该函数被忠实地四舍五入。 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;
}