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

Question

标准正态分布的概率密度函数定义为 e^-x2/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⁹ ulps用于 IEEE-754 binary64 double 。

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

Answer 1

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

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

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

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

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 binary32 的 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 实现看起来几乎相同，只是使用的常量值不同:

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。