这里给出的普通CDF的实现是单精度近似，其已经用double替换了float，因此仅精确到7或8个有效(十进制)数字。

对于Hart的双精度近似的VB实现，请参见West对累积正态函数的更好近似的图2。

编辑：我将West的实现转换为C ++：

double

phi(double x)

{

static const double RT2PI = sqrt(4.0*acos(0.0));

static const double SPLIT = 7.07106781186547;

static const double N0 = 220.206867912376;

static const double N1 = 221.213596169931;

static const double N2 = 112.079291497871;

static const double N3 = 33.912866078383;

static const double N4 = 6.37396220353165;

static const double N5 = 0.700383064443688;

static const double N6 = 3.52624965998911e-02;

static const double M0 = 440.413735824752;

static const double M1 = 793.826512519948;

static const double M2 = 637.333633378831;

static const double M3 = 296.564248779674;

static const double M4 = 86.7807322029461;

static const double M5 = 16.064177579207;

static const double M6 = 1.75566716318264;

static const double M7 = 8.83883476483184e-02;

const double z = fabs(x);

double c = 0.0;

if(z<=37.0)

{

const double e = exp(-z*z/2.0);

if(z

{

const double n = (((((N6*z + N5)*z + N4)*z + N3)*z + N2)*z + N1)*z + N0;

const double d = ((((((M7*z + M6)*z + M5)*z + M4)*z + M3)*z + M2)*z + M1)*z + M0;

c = e*n/d;

}

else

{

const double f = z + 1.0/(z + 2.0/(z + 3.0/(z + 4.0/(z + 13.0/20.0))));

c = e/(RT2PI*f);

}

}

return x<=0.0 ? c : 1-c;

}

请注意，我已将表达式重新排列为更熟悉的系列形式和连续分数近似值。 West代码中的最后一个神奇数字是2＆pi;的平方根，我通过利用身份acos(0)=＆frac12在第一行推迟到编译器; ＆PI ;.

我已经对魔术数字进行了三次检查，但我总是错误地输入了一些东西。如果您发现拼写错误，请评论！

John Cook在他的回答中使用的测试数据的结果是

x phi Mathematica

-3 1.3498980316301150e-003 0.00134989803163

-1 1.5865525393145702e-001 0.158655253931

0 5.0000000000000000e-001 0.5

0.5 6.9146246127401301e-001 0.691462461274

2.1 9.8213557943718344e-001 0.982135579437

我从他们同意Mathematica结果的所有数字这一事实中得到一些小小的安慰。