使用C标准数学库精确计算标准正态分布的CDF
标准C数学库不提供计算标准正态分布的CDF函数normcdf() 。 但它确实提供了密切相关的函数：错误函数erf()和补充错误函数erfc() 。 计算CDF的最快方法通常是通过误差函数，使用预定义的常量M_SQRT1_2来表示√½：
double normcdf (double a) { return 0.5 + 0.5 * erf (M_SQRT1_2 * a); }
显然，这在负半平面中遭受大量的减法消除，并且不适合大多数应用。 由于使用erfc()可以很容易地避免取消问题，但是erf()性能通常比erf()低一些，因此最常推荐的计算是：
double normcdf (double a) { return 0.5 * erfc (-M_SQRT1_2 * a); }
一些测试表明，在负半平面中产生的最大ulp误差仍然相当大。 使用精确到0.51 ulps的erfc()的双精度实现，可以在normcdf()观察到高达1705.44 ulps的误差。 这里的问题是erfc()输入中的计算误差被erfc()固有的指数缩放放大（有关由取幂引起的误差放大的解释，请参阅此答案 ）。
下面的文章展示了当浮点操作数与任意精度常数（如√½）相乘时，如何实现（几乎）正确舍入的产品：
Nicolas Brisebarre和Jean-Michel Muller，“正确舍入乘以任意精度常数”， IEEE Transactions on Computers ，Vol。 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则用于根据h处的互补误差函数的局部斜率插入erfc()结果。
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-754 binary64 ），实现了以下实现：
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-754 binary32 ）是：
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; }