使用 C 标准数学库精确计算标准正态分布的 PDF

Question

标准正态分布的概率密度函数定义为e < sup >-x < sup > 2 /2 /√(2π)。这可以用简单的方式转换成C代码。单精度实现的示例可能是:

float my_normpdff (float a)
{
    return 0x1.988454p-2f * my_expf (-0.5f * a * a); /* 1/sqrt(2*pi) */
}

虽然此代码没有过早的下限溢位，但存在精度问题，因为计算²/2时产生的误差会被后续的幂运算放大。人们可以通过针对更高精度引用的测试轻松证明这一点。确切的误差将根据所使用的exp（）或expf（）实现的准确性而有所不同；对于忠实四舍五入的幂函数，人们通常会观察到IEEE-754binary32单精度的最大误差约为2⁶ulps，IEEE-754binary64双精度的最大误差约为2⁹ulps。

如何以合理有效的方式解决准确性问题？一种简单的方法是采用更高精度的中间计算，例如对< code>float实现使用< code>double计算。但是，如果不容易获得更高精度的浮点运算，这种方法不适用于< code>double实现，并且如果< code>double运算比< code>float计算昂贵得多，例如在许多GPU上，这种方法对于< code>float实现可能是低效的。

Answer 1

问题中提出的精度问题可以通过使用有限数量的double-float或double--doublecomputer来有效地解决，这得益于融合乘法加法（FMA）运算的使用。

从< code>C99开始，此操作可通过标准数学函数< code>fmaf(a，b，c)和< code>fma(a，b，c)进行，它们计算a*b c，无需对中间产品进行舍入。虽然这些函数直接映射到几乎所有现代处理器上的快速硬件操作，但它们可能使用旧平台上的仿真代码，在这种情况下，它们可能非常慢。

这允许仅使用两个操作以两倍于正常精度的精度计算乘积，从而产生一对头：尾的本机精度数字：

prod_hi = a * b            // head
prod_lo = FMA (a, b, -hi)  // tail

结果的高阶位可以传递给求幂运算，而低阶位则用于通过线性插值提高结果的精度，利用e^x是其自身导数的事实：

e = exp (prod_hi) + exp (prod_hi) * prod_lo  // exp (a*b)

这使我们能够消除原始实现的大部分错误。另一个较小的计算误差来源是常数1的有限精度/√表示为（2π）。通过对提供两倍本机精度的html" target="_blank">常量使用head:tail表示法，并计算：

r = FMA (const_hi, x, const_lo * x)  // const * x

下面的论文指出，这种技术甚至可以对某些任意精度的常数产生正确舍入的乘法:

Nicolas Brisebarre和Jean Michel Muller，“任意精度常数的正确四舍五入乘法”，《IEEE计算机汇刊》，第57卷，第2期，2008年2月，第165-174页

结合这两种技术，并考虑到一些涉及NaN的角落案例，我们得出了以下基于IEEE-754二进制32的float实现：

float my_normpdff (float a)
{
    const float RCP_SQRT_2PI_HI =  0x1.988454p-02f; /* 1/sqrt(2*pi), msbs */
    const float RCP_SQRT_2PI_LO = -0x1.857936p-27f; /* 1/sqrt(2*pi), lsbs */
    float ah, sh, sl, ea;

    ah = -0.5f * a;
    sh = a * ah;
    sl = fmaf (a, ah, 0.5f * a * a); /* don't flip "sign bit" of NaN argument */
    ea = expf (sh);
    if (ea != 0.0f) ea = fmaf (sl, ea, ea); /* avoid creation of NaN */
    return fmaf (RCP_SQRT_2PI_HI, ea, RCP_SQRT_2PI_LO * ea);
}

基于IEEE-754 binary64的相应双重实现看起来几乎相同，除了使用不同的常量值：

double my_normpdf (double a)
{
    const double RCP_SQRT_2PI_HI =  0x1.9884533d436510p-02; /* 1/sqrt(2*pi), msbs */
    const double RCP_SQRT_2PI_LO = -0x1.cbc0d30ebfd150p-56; /* 1/sqrt(2*pi), lsbs */
    double ah, sh, sl, ea;

    ah = -0.5 * a;
    sh = a * ah;
    sl = fma (a, ah, 0.5 * a * a); /* don't flip "sign bit" of NaN argument */
    ea = exp (sh);
    if (ea != 0.0) ea = fma (sl, ea, ea); /* avoid creation of NaN */
    return fma (RCP_SQRT_2PI_HI, ea, RCP_SQRT_2PI_LO * ea);
}

这些实现的准确性分别取决于标准数学函数< code>expf()和< code>exp()的准确性。在C数学库提供精确舍入版本的情况下，上述两种实现的最大误差通常小于2.5 ulps。