標準正態分布函數的近似計算

#define pi (3.1415926535898)
#define a0  (0.33267)
#define a1  (0.4361836)
#define a2 (-0.1201676)
#define a3 (0.9372986)
double test(double x)
{
double t = (1.0/(1.0+a0*fabs((x))));
return exp(-x*x/2.0) * t * (a1 + a2*t + a3*t*t) / sqrt(2.0*pi);
}
double testfun(double x)
{
if(x == 0)
{
return 0.5;
}
else if(x > 0)
{
return 1.0 - test(x)/sqrt(2*pi);
}
else
{
return test(x)/sqrt(2*pi);
}
}