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

Question

标准正态分布的概率密度函数定义为e^-x²/2 / √(2π)。这可以直接呈现为 C 代码。示例单精度实现可能是：

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

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

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

Answer 1

问题中提出的准确性问题可以通过使用有限数量的双float或双double计算来有效且高效地解决，使用融合乘加 (FMA) 运算。

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

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

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π) 的精度有限。这可以通过对提供两倍本机精度的常量使用 head:tail 表示来改进，并计算：

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

下面的文章指出，这种技术甚至可以为一些任意精度的常量产生正确的舍入乘法：

Nicolas Brisebarre 和 Jean-Michel Muller，"Correctly rounded multiplication by arbitrary precision constants"，IEEE 计算机汇刊，卷。 57，第 2 期，2008 年 2 月，第 165-174 页

结合这两种技术，并处理一些涉及 NaN 的极端情况，我们得出以下 float 基于 IEEE-754 binary32 的实现：

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 实现看起来几乎相同，除了使用的常量值不同：

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);
}

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

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

Accurate computation of PDF of standard normal distribution using C standard math library

c

math

floating-point

floating-accuracy