使用 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() 实现的准确性而有所不同；对于忠实四舍五入的幂函数，通常会观察到 IEEE-754 binary32 单精度的最大误差约为 2⁶ ulps，IEEE-754 @ 约为 2⁹ ulps 987654325@双精度。

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

Answer 1

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

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

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

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

结果的高位可以传递给求幂，而低位用于通过线性插值提高结果的准确性，利用 e^{x sup> 是它自己的派生词：}

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，“任意精度常数的正确舍入乘法”，IEEE Transactions on Computers，Vol。 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);
}

对应的 double 实现基于 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);
}

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