c - 使用C标准数学库准确计算标准正态分布的CDF

标准C数学库不提供计算标准正态分布CDF的函数，normcdf()但是，它确实提供了密切相关的函数：错误函数，erf()和互补错误函数，erfc()。计算CDF的最快方法通常是通过误差函数，使用预定义的常数m_sqrt1_2来表示√1/2：

double normcdf (double a)
{
    return 0.5 + 0.5 * erf (M_SQRT1_2 * a);
}

显然，这种方法在负半平面上存在大量的减法消除，不适合大多数应用。由于使用erfc()可以很容易地避免取消问题，但其性能通常比erf()低一些，因此最常用的建议计算是：

double normcdf (double a)
{
    return 0.5 * erfc (-M_SQRT1_2 * a);
}

一些测试表明，在负半平面中发生的最大ULP误差仍然很大。使用精度erfc()达到0.51ulps的双精度实现，可以观察到
normcdf()中高达1705.44 ulps。这里的问题是erfc()输入中的计算误差被erfc()固有的指数标度放大（参见answer
解释指数运算引起的误差放大率）。
下面的文章展示了如何在将浮点操作数与任意精度常数（如√1/2）相乘时获得（几乎）正确的四舍五入乘积：
Nicolas Brisebarre和Jean-Michel Muller，“任意精度常数的正确四舍五入乘法”，《计算机上的IEEE交易》，第57卷，第2期，2008年2月，第165-174页
本文所提出的方法依赖于融合乘法加法运算，该运算可用于所有通用处理器体系结构的最新实现，并通过标准的数学函数fma()在c中公开。这将导致以下版本：

double normcdf (double a)
{
    double SQRT_HALF_HI =  0x1.6a09e667f3bcd0p-01; //  7.0710678118654757e-01
    double SQRT_HALF_LO = -0x1.bdd3413b264560p-55; // -4.8336466567264567e-17

    return 0.5 * erfc (fma (-SQRT_HALF_HI, a, -SQRT_HALF_LO * a));
}

测试表明，这比以前的版本减少了大约一半的最大误差。使用与以前相同的高度精确的erfc()实现，最大观测误差为842.71 ULPS。这与提供基本数学函数的通常目标相去甚远，其误差至多为几个ulps。
是否有一种高效的方法允许精确计算normcdf()，并且只使用标准c数学库中可用的函数？

最佳答案

解决问题中概述的方法的精度限制的一种方法是有限地使用双重计算。这包括以头/尾方式将-sqrt (0.5) * a作为一对double变量h和l进行计算。产品的高阶部分h传递到erfc()，而低阶部分l则用于根据erfc()处互补误差函数的局部斜率插值h结果。
erfc（x）的导数是-2*exp（-x*x）/√π。但是，我们希望避免计算exp（-x*x）的开销。当x>0时，erfc（x）~=2*exp（-x*x）/（√π*（x+sqrt（x*x+4/π））因此，渐近地，
erfc’（x）~=-2*x*erfc（x），其结果是对于l h，erfc（h+l）~=erfc（h）-2*h*l*erfc（h）。后一项的否定可以很容易地引入到l的计算中。其中一个实现了以下双精度（使用ieee-754binary64）：

double my_normcdf (double a)
{
    double h, l, r;
    const double SQRT_HALF_HI =  0x1.6a09e667f3bcd0p-01; //  7.0710678118654757e-01
    const double SQRT_HALF_LO = -0x1.bdd3413b264560p-55; // -4.8336466567264567e-17

    /* clamp input as normcdf(x) is either 0 or 1 asymptotically */
    if (fabs (a) > 38.625) a = (a < 0.0) ? -38.625 : 38.625;

    h = fma (-SQRT_HALF_HI, a, -SQRT_HALF_LO * a);
    l = fma (SQRT_HALF_LO, a, fma (SQRT_HALF_HI, a, h));
    r = erfc (h);
    if (h > 0.0) r = fma (2.0 * h * l, r, r);
    return 0.5 * r;
}

使用以前使用的相同的erfc()实现的最大观测误差是1.96个ULPS。相应的单精度实现（使用ieee-754binary32）是：

float my_normcdff (float a)
{
    float h, l, r;
    const float SQRT_HALF_HI = 0x1.6a09e6p-01f; // 7.07106769e-1
    const float SQRT_HALF_LO = 0x1.9fcef4p-27f; // 1.21016175e-8

    /* clamp input as normcdf(x) is either 0 or 1 asymptotically */
    if (fabsf (a) > 14.171875f) a = (a < 0.0f) ? -14.171875f : 14.171875f;

    h = fmaf (-SQRT_HALF_HI, a, -SQRT_HALF_LO * a);
    l = fmaf (SQRT_HALF_LO, a, fmaf (SQRT_HALF_HI, a, h));
    r = erfcf (h);
    if (h > 0.0f) r = fmaf (2.0f * h * l, r, r);
    return 0.5f * r;
}