'erf' function in C

[James Harlacher]

Hi,
Are there any C compilers which have the erf function (from
probability) as part of their math libraries? Or are there any math
libraries available to download which implement this function?
Thanks,
James Harlacher

 

[jacob navia]

Hi,
Are there any C compilers which have the erf function (from
probability) as part of their math libraries? Or are there any math
libraries available to download which implement this function?
Thanks,
James Harlacher

Yes. lcc-win32 uses Sun's implementation.

/* @(#)s_erf.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/

/* double erf(double x)
* double erfc(double x)
* x
* 2 |\
* erf(x) = --------- | exp(-t*t)dt
* sqrt(pi) \|
* 0
*
* erfc(x) = 1-erf(x)
* Note that
* erf(-x) = -erf(x)
* erfc(-x) = 2 - erfc(x)
*
* Method:
* 1. For |x| in [0, 0.84375]
* erf(x) = x + x*R(x^2)
* erfc(x) = 1 - erf(x) if x in [-.84375,0.25]
* = 0.5 + ((0.5-x)-x*R) if x in [0.25,0.84375]
* where R = P/Q where P is an odd poly of degree 8 and
* Q is an odd poly of degree 10.
* -57.90
* | R - (erf(x)-x)/x | <= 2
*
*
* Remark. The formula is derived by noting
* erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
* and that
* 2/sqrt(pi) = 1.128379167095512573896158903121545171688
* is close to one. The interval is chosen because the fix
* point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
* near 0.6174), and by some experiment, 0.84375 is chosen to
* guarantee the error is less than one ulp for erf.
*
* 2. For |x| in [0.84375,1.25], let s = |x| - 1, and
* c = 0.84506291151 rounded to single (24 bits)
* erf(x) = sign(x) * (c + P1(s)/Q1(s))
* erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0
* 1+(c+P1(s)/Q1(s)) if x < 0
* |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
* Remark: here we use the taylor series expansion at x=1.
* erf(1+s) = erf(1) + s*Poly(s)
* = 0.845.. + P1(s)/Q1(s)
* That is, we use rational approximation to approximate
* erf(1+s) - (c = (single)0.84506291151)
* Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
* where
* P1(s) = degree 6 poly in s
* Q1(s) = degree 6 poly in s
*
* 3. For x in [1.25,1/0.35(~2.857143)],
* erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
* erf(x) = 1 - erfc(x)
* where
* R1(z) = degree 7 poly in z, (z=1/x^2)
* S1(z) = degree 8 poly in z
*
* 4. For x in [1/0.35,28]
* erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
* = 2.0 - (1/x)*exp(-x*x-0.5625+R2/S2) if -6<x<0
* = 2.0 - tiny (if x <= -6)
* erf(x) = sign(x)*(1.0 - erfc(x)) if x < 6, else
* erf(x) = sign(x)*(1.0 - tiny)
* where
* R2(z) = degree 6 poly in z, (z=1/x^2)
* S2(z) = degree 7 poly in z
*
* Note1:
* To compute exp(-x*x-0.5625+R/S), let s be a single
* precision number and s := x; then
* -x*x = -s*s + (s-x)*(s+x)
* exp(-x*x-0.5626+R/S) =
* exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
* Note2:
* Here 4 and 5 make use of the asymptotic series
* exp(-x*x)
* erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
* x*sqrt(pi)
* We use rational approximation to approximate
* g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
* Here is the error bound for R1/S1 and R2/S2
* |R1/S1 - f(x)| < 2**(-62.57)
* |R2/S2 - f(x)| < 2**(-61.52)
*
* 5. For inf > x >= 28
* erf(x) = sign(x) *(1 - tiny) (raise inexact)
* erfc(x) = tiny*tiny (raise underflow) if x > 0
* = 2 - tiny if x<0
*
* 7. Special case:
* erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
* erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
* erfc/erf(NaN) is NaN
*/


static const double tiny = 1e-300,
half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
/* c = (float)0.84506291151 */
erx = 8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
/*
* Coefficients for approximation to erf on [0,0.84375]
*/
efx = 1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
efx8= 1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
pp0 = 1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
pp1 = -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
pp2 = -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
pp3 = -5.77027029648944159157e-03, /* 0xBF77A291, 0x236668E4 */
pp4 = -2.37630166566501626084e-05, /* 0xBEF8EAD6, 0x120016AC */
qq1 = 3.97917223959155352819e-01, /* 0x3FD97779, 0xCDDADC09 */
qq2 = 6.50222499887672944485e-02, /* 0x3FB0A54C, 0x5536CEBA */
qq3 = 5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */
qq4 = 1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */
qq5 = -3.96022827877536812320e-06, /* 0xBED09C43, 0x42A26120 */
/*
* Coefficients for approximation to erf in [0.84375,1.25]
*/
pa0 = -2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */
pa1 = 4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */
pa2 = -3.72207876035701323847e-01, /* 0xBFD7D240, 0xFBB8C3F1 */
pa3 = 3.18346619901161753674e-01, /* 0x3FD45FCA, 0x805120E4 */
pa4 = -1.10894694282396677476e-01, /* 0xBFBC6398, 0x3D3E28EC */
pa5 = 3.54783043256182359371e-02, /* 0x3FA22A36, 0x599795EB */
pa6 = -2.16637559486879084300e-03, /* 0xBF61BF38, 0x0A96073F */
qa1 = 1.06420880400844228286e-01, /* 0x3FBB3E66, 0x18EEE323 */
qa2 = 5.40397917702171048937e-01, /* 0x3FE14AF0, 0x92EB6F33 */
qa3 = 7.18286544141962662868e-02, /* 0x3FB2635C, 0xD99FE9A7 */
qa4 = 1.26171219808761642112e-01, /* 0x3FC02660, 0xE763351F */
qa5 = 1.36370839120290507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */
qa6 = 1.19844998467991074170e-02, /* 0x3F888B54, 0x5735151D */
/*
* Coefficients for approximation to erfc in [1.25,1/0.35]
*/
ra0 = -9.86494403484714822705e-03, /* 0xBF843412, 0x600D6435 */
ra1 = -6.93858572707181764372e-01, /* 0xBFE63416, 0xE4BA7360 */
ra2 = -1.05586262253232909814e+01, /* 0xC0251E04, 0x41B0E726 */
ra3 = -6.23753324503260060396e+01, /* 0xC04F300A, 0xE4CBA38D */
ra4 = -1.62396669462573470355e+02, /* 0xC0644CB1, 0x84282266 */
ra5 = -1.84605092906711035994e+02, /* 0xC067135C, 0xEBCCABB2 */
ra6 = -8.12874355063065934246e+01, /* 0xC0545265, 0x57E4D2F2 */
ra7 = -9.81432934416914548592e+00, /* 0xC023A0EF, 0xC69AC25C */
sa1 = 1.96512716674392571292e+01, /* 0x4033A6B9, 0xBD707687 */
sa2 = 1.37657754143519042600e+02, /* 0x4061350C, 0x526AE721 */
sa3 = 4.34565877475229228821e+02, /* 0x407B290D, 0xD58A1A71 */
sa4 = 6.45387271733267880336e+02, /* 0x40842B19, 0x21EC2868 */
sa5 = 4.29008140027567833386e+02, /* 0x407AD021, 0x57700314 */
sa6 = 1.08635005541779435134e+02, /* 0x405B28A3, 0xEE48AE2C */
sa7 = 6.57024977031928170135e+00, /* 0x401A47EF, 0x8E484A93 */
sa8 = -6.04244152148580987438e-02, /* 0xBFAEEFF2, 0xEE749A62 */
/*
* Coefficients for approximation to erfc in [1/.35,28]
*/
rb0 = -9.86494292470009928597e-03, /* 0xBF843412, 0x39E86F4A */
rb1 = -7.99283237680523006574e-01, /* 0xBFE993BA, 0x70C285DE */
rb2 = -1.77579549177547519889e+01, /* 0xC031C209, 0x555F995A */
rb3 = -1.60636384855821916062e+02, /* 0xC064145D, 0x43C5ED98 */
rb4 = -6.37566443368389627722e+02, /* 0xC083EC88, 0x1375F228 */
rb5 = -1.02509513161107724954e+03, /* 0xC0900461, 0x6A2E5992 */
rb6 = -4.83519191608651397019e+02, /* 0xC07E384E, 0x9BDC383F */
sb1 = 3.03380607434824582924e+01, /* 0x403E568B, 0x261D5190 */
sb2 = 3.25792512996573918826e+02, /* 0x40745CAE, 0x221B9F0A */
sb3 = 1.53672958608443695994e+03, /* 0x409802EB, 0x189D5118 */
sb4 = 3.19985821950859553908e+03, /* 0x40A8FFB7, 0x688C246A */
sb5 = 2.55305040643316442583e+03, /* 0x40A3F219, 0xCEDF3BE6 */
sb6 = 4.74528541206955367215e+02, /* 0x407DA874, 0xE79FE763 */
sb7 = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */

extern double exp(double);
extern double fabs(double);
double erf(double x)
{
int n0,hx,ix,i;
double R,S,P,Q,s,y,z,r;
n0 = ((*(int*)&one)>>29)^1;
hx = *(n0+(int*)&x);
ix = hx&0x7fffffff;
if(ix>=0x7ff00000) { /* erf(nan)=nan */
i = ((unsigned)hx>>31)<<1;
return (double)(1-i)+one/x; /* erf(+-inf)=+-1 */
}

if(ix < 0x3feb0000) { /* |x|<0.84375 */
if(ix < 0x3e300000) { /* |x|<2**-28 */
if (ix < 0x00800000)
return 0.125*(8.0*x+efx8*x); /*avoid underflow */
return x + efx*x;
}
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
return x + x*y;
}
if(ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
s = fabs(x)-one;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
if(hx>=0) return erx + P/Q; else return -erx - P/Q;
}
if (ix >= 0x40180000) { /* inf>|x|>=6 */
if(hx>=0) return one-tiny; else return tiny-one;
}
x = fabs(x);
s = one/(x*x);
if(ix< 0x4006DB6E) { /* |x| < 1/0.35 */
R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/0.35 */
R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7))))));
}
z = x;
*(1-n0+(int*)&z) = 0;
r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
if(hx>=0) return one-r/x; else return r/x-one;
}

double erfc(double x)
{
int n0,hx,ix;
double R,S,P,Q,s,y,z,r;
n0 = ((*(int*)&one)>>29)^1;
hx = *(n0+(int*)&x);
ix = hx&0x7fffffff;
if(ix>=0x7ff00000) { /* erfc(nan)=nan */
/* erfc(+-inf)=0,2 */
return (double)(((unsigned)hx>>31)<<1)+one/x;
}

if(ix < 0x3feb0000) { /* |x|<0.84375 */
if(ix < 0x3c700000) /* |x|<2**-56 */
return one-x;
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
if(hx < 0x3fd00000) { /* x<1/4 */
return one-(x+x*y);
} else {
r = x*y;
r += (x-half);
return half - r ;
}
}
if(ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
s = fabs(x)-one;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
if(hx>=0) {
z = one-erx; return z - P/Q;
} else {
z = erx+P/Q; return one+z;
}
}
if (ix < 0x403c0000) { /* |x|<28 */
x = fabs(x);
s = one/(x*x);
if(ix< 0x4006DB6D) { /* |x| < 1/.35 ~ 2.857143*/
R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/.35 ~ 2.857143 */
if(hx<0&&ix>=0x40180000) return two-tiny;/* x < -6 */
R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7))))));
}
z = x;
*(1-n0+(int*)&z) = 0;
r = exp(-z*z-0.5625)*
exp((z-x)*(z+x)+R/S);
if(hx>0) return r/x; else return two-r/x;
} else {
if(hx>0) return tiny*tiny; else return two-tiny;
}
}

 

[Tim Prince]

Hi,
Are there any C compilers which have the erf function (from
probability) as part of their math libraries? Or are there any math
libraries available to download which implement this function?

What happened to your web browser? Most libraries on which gcc is
implemented, including glibc and newlib, have this function.

 

[Sidney Cadot]

Hi,
Are there any C compilers which have the erf function (from
probability) as part of their math libraries? Or are there any math
libraries available to download which implement this function?

Any C99 compiler has it; gcc has it.

Best regards,

Sidney

 

[George Marsaglia]

Hi,
Are there any C compilers which have the erf function (from
probability) as part of their math libraries? Or are there any math
libraries available to download which implement this function?
Thanks,
James Harlacher

I have always considered the erf function
a damned nuisance---not for its potential use, but
for the way that it must be used, and the many complicated
methods used to get what can otherwise be expressed
by means of simple mathematics and simple programming
that lets the CPU do the work for you.

erf is often not available on particular systems,
(for example on the gcc under DOS that I use), and may
rely on obscure methods---fits by piecewise polynomials
rational functions, continued fractions, Chebychev
or asymptotic expansions for the tail---and have
less that hoped-for accuracy.
Examples are the widely used C programs based on
Algorithm AS66 Applied Statistics v22, p 424 by Hill,
from which you are lucky to get 7 digits of accuracy,
or the ponderous 300-line code from Sun Microsystems,
listed in one of the postings above.

Virtually all calls for an erf(z) evaluation
are likely to be for the normal probability distribution:

Phi(x) = integral exp(-t^2/2), t = -infinity to x,

which must be obtained via erf by means of
Phi(x) = .5 + .5*erf(x/sqrt(2)),
another source of irritation for erf (or rather Phi) users.

A procedure that can determine Phi(x) nearly to the limit
available in double precision can be based on the following
method, which I developed the 1960's:

Let phi(t) be the normal density, phi(t) = exp(-t^2/2)/sqrt(2*Pi),
so that
Phi(x) = integral phi(t), t = -infinity to x.
and
cPhi(x) = 1-Phi(x) = integral phi(t), t = x to infinity.

If you define

R(x) = cPhi(x)/phi(x)

then the function R is a well behaved function that has
a rapidly convergent Taylor expansion, an expansion
with terms easily determined by a two-step recursion.
(Try to derive the recursion yourself:
r[k+1](x)=x*r[k](x)+k*r[k-1](x),
where r[k](x) is the kth derivative of R(x).)

Furthermore, the convergence is fast ---or at least
accurate---enough that the Taylor expansion about
x = 0 can be used to determine Phi(x) to some 15-16
decimal places even for x as large as 6 or 7---farther
than most applications are likely to require.
So, a few lines of code and no tables are all one needs
for a very accurate value of Phi(x), the normal probability distribution:
------------------------------------------------------
#include <stdio.h>
#include <math.h>

double Phi(double x)
{long double s,t=0,a=1.2533141373155L,z=0,b=-1,pwr=1;
int i;
s=a+b*x;
for(i=2;s!=t;i+=2)
{ a=(a+z*b)/i;
b=(b+z*a)/(i+1);
pwr=pwr*x*x;
t=s;
s+=pwr*(a+x*b);
}
return 1.-s*exp(-.5*x*x-.91893853320467274178L);
}

int main(){
double x;
while(1){
printf("Enter your x value:");
scanf("%lf",&x);
printf("Normal Prob(X<%f)=%20.16f\n",x,Phi(x));
}
}
/* An indication of the accuracy of this Phi function
is given by the following examples, with the true
value (to 20 places) given under the resulting Phi(x):
x Phi(x)
0.123 0.5489464510164369
.54894645101643675909
1.2 0.8849303297782918
.88493032977829173198
2.4 0.9918024640754040
.99180246407540387055
6.1 0.9999999994696578
.99999999946965767370
-6.1 0.0000000005303427
.00000000053034232638
-1.1 0.1356660609463828
.13566606094638267517
7.2 0.9999999999996990
.99999999999969893721
See what your system gives and compare with
an erf function if it has one.
*/
----------------------------------------------------------
Cut, paste and try. You may like it.

George Marsaglia

(Note: With "long double"s I hope to invoke the 80-bit
floating point processor, providing a little more
accuracy that the IEEE 64-bit floating point standard.
Does your system provide similar accuracy? I use
gcc under DOS on Intel machines with W98, Me or XP.)

 

[Matthias Klaey]

Hi,
Are there any C compilers which have the erf function (from
probability) as part of their math libraries? Or are there any math
libraries available to download which implement this function?
Thanks,
James Harlacher

[...]

A procedure that can determine Phi(x) nearly to the limit
available in double precision can be based on the following
method, which I developed the 1960's:
[...]

/* An indication of the accuracy of this Phi function
is given by the following examples, with the true
value (to 20 places) given under the resulting Phi(x):
x Phi(x)
0.123 0.5489464510164369
.54894645101643675909
1.2 0.8849303297782918
.88493032977829173198
2.4 0.9918024640754040
.99180246407540387055
6.1 0.9999999994696578
.99999999946965767370
-6.1 0.0000000005303427
.00000000053034232638
-1.1 0.1356660609463828
.13566606094638267517
7.2 0.9999999999996990
.99999999999969893721
See what your system gives and compare with
an erf function if it has one.
*/

[...]

Is there someplace a table available of the "true" values of the
normal distribution up to 20 places?

With kind regards
Matthias Kläy

 

[P.J. Plauger]

I have always considered the erf function
a damned nuisance---not for its potential use, but
for the way that it must be used, and the many complicated
methods used to get what can otherwise be expressed
by means of simple mathematics and simple programming
that lets the CPU do the work for you.

erf is often not available on particular systems,
(for example on the gcc under DOS that I use), and may
rely on obscure methods---fits by piecewise polynomials
rational functions, continued fractions,

Obscure? It's hard to find a math function that doesn't use
at least one of these techniques. And why should you care,
if it gets the right answer with reasonable speed?

Chebychev
or asymptotic expansions for the tail---

Ah, now you're thinking about erfc. erf saturates quickly
for values away from zero, so it has no tail problems.

and have
less that hoped-for accuracy.

Now *that's* a legitimate gripe. We find accuracies all over
the map for even the commonest math functions. For erf and
its exotic brethren, the errors are often astonishing.

Examples are the widely used C programs based on
Algorithm AS66 Applied Statistics v22, p 424 by Hill,
from which you are lucky to get 7 digits of accuracy,

Indeed. We've found a rich literature of mediocre algorithms
for math functions. Sadly, these algorithms are widely used.

or the ponderous 300-line code from Sun Microsystems,
listed in one of the postings above.

That ponderousness is mostly about getting the erfc tail
right for positive x, which is very difficult to do. The
erf part doesn't use it. In fairness to Sun, the function
is pretty damned accurate.

P.J. Plauger
Dinkumware, Ltd.
http://www.dinkumware.com

 

[Axel Vogt]

I have always considered the erf function
a damned nuisance---not for its potential use, but
for the way that it must be used, and the many complicated
methods used to get what can otherwise be expressed
by means of simple mathematics and simple programming
that lets the CPU do the work for you. ....
A procedure that can determine Phi(x) nearly to the limit
available in double precision can be based on the following
method, which I developed the 1960's:

....

I once did this (with a small lookup table as you had it
in earlier versions) for Excel. Even in that poor setting
it worked well (and with some stupid language tricks one
could go there for 18 digits).

 

[Sidney Cadot]

Is there someplace a table available of the "true" values of the
normal distribution up to 20 places?

<wildly offtopic for c.l.c>

Let me share one of my favorite hacks with you! Go to
integrals.wolfram.com, and submit something like:

Erf[0.3`200]

(0.3`200 means 0.3, with a precision of 200 digits)

You will get back the integral of your input (as a picture,
unfortunately), which is just x multiplied by the constant you're after
of course; Mathematica is pretty good at getting the precision right.

Best regards,

Sidney

 

[Martin Dickopp]

Is there someplace a table available of the "true" values of the
normal distribution up to 20 places?

<wildly offtopic for c.l.c>

Let me share one of my favorite hacks with you! Go to
integrals.wolfram.com, and submit something like:

Erf[0.3`200]

(0.3`200 means 0.3, with a precision of 200 digits)

You will get back the integral of your input (as a picture,
unfortunately),

Let me share an even better hack with you. Open the picture in a new
browser window and change the "format=StandardForm" part of the URL to
"format=OutputForm". You'll get the result as text.

which is just x multiplied by the constant you're after of course;
Mathematica is pretty good at getting the precision right.

Actually, unless my understanding of the term "precision" is completely
wrong, Mathematica seems to be pretty /bad/ at getting the precision
right: It returns a number with only 199 digits after the decimal point!

OTOH, the Emacs calculator gets the number /and/ the precision right.

Martin

 

[Sidney Cadot]

Let me share an even better hack with you. Open the picture in a new
browser window and change the "format=StandardForm" part of the URL to
"format=OutputForm". You'll get the result as text.

Nice one...! To the OP: assuming you're using a Unix-like OS, this means
something like

curl "http://integrals.wolfram.com/
graphics.cgi?format=OutputForm&expr=Erf[0.3\`200]"

....can be used in a script to generate your table, with a bit of
scripting magic.

Actually, unless my understanding of the term "precision" is completely
wrong, Mathematica seems to be pretty /bad/ at getting the precision
right: It returns a number with only 199 digits after the decimal point!

This can be because of two things. Mathematica distinguishes "precision"
(total number of significant digits) from "accuracy" (number of
significant digits after the decimal point), so the Erf argument is
already given with an accuracy of 199 digits. Alternatively, the result
of an unary operation doesn't necessarily have the same precision as the
argument, e.g. what is sin(x), with x=1.00000E20 ?

OTOH, the Emacs calculator gets the number /and/ the precision right.

There's people in the world that use web-browsers but don't use emacs...
I kid you not.

Best regards, Sidney

 

[bv]

Are there any C compilers which have the erf function (from
probability) as part of their math libraries? Or are there any math
libraries available to download which implement this function?

I'm aware of two compilers, one C (Sun) and one Fortran (Watcom). I'm
guessing that Watcom C compiler might have it too. There are, of course,
several "erf" codes available on netlib, including one from Sun in C.
It's located in "fdlibm" lib, at http://netlib.org.

 

[George Marsaglia]

Hi,
Are there any C compilers which have the erf function (from
probability) as part of their math libraries? Or are there any math
libraries available to download which implement this function?
Thanks,
James Harlacher

[...]

A procedure that can determine Phi(x) nearly to the limit
available in double precision can be based on the following
method, which I developed the 1960's:
[...]

/* An indication of the accuracy of this Phi function
is given by the following examples, with the true
value (to 20 places) given under the resulting Phi(x):
x Phi(x)
0.123 0.5489464510164369
.54894645101643675909
1.2 0.8849303297782918
.88493032977829173198
2.4 0.9918024640754040
.99180246407540387055
6.1 0.9999999994696578
.99999999946965767370
-6.1 0.0000000005303427
.00000000053034232638
-1.1 0.1356660609463828
.13566606094638267517
7.2 0.9999999999996990
.99999999999969893721
See what your system gives and compare with
an erf function if it has one.
*/

[...]

Is there someplace a table available of the "true" values of the
normal distribution up to 20 places?

With kind regards
Matthias Kläy

Sorry, I should have mentioned that the values I used
were based on setting Digits:=20 or 30 in Maple, then
calling .5+.5*erf(x/sqrt(2.)), (Yes, Maple is still stuck
with erf for lack of a standard for the normal distribution,
although it may be avoided by the rather awkward
stats[statevalf,cdf,normald](x)

Mathematika and other extended precision math packages
would also serve. Of course the old standby of
looking up tabled values, for example the Handbook of
Mathematical Functions edited by Abramowitz and Stegun,
will still serve, but they are likely to list evenly spaced
x's and perhaps less than 20 digits for Phi(x).

George Marsaglia

 

[George Russell]

George Marsaglia wrote (snipped):

Let phi(t) be the normal density, phi(t) = exp(-t^2/2)/sqrt(2*Pi),
so that
Phi(x) = integral phi(t), t = -infinity to x.
and
cPhi(x) = 1-Phi(x) = integral phi(t), t = x to infinity.

If you define

R(x) = cPhi(x)/phi(x)

then the function R is a well behaved function that has
a rapidly convergent Taylor expansion, an expansion
with terms easily determined by a two-step recursion.

Very nice. But can someone explain why you need
R = cPhi / phi
rather than just
S = Phi/phi
? The maths seems to work out just as well in the second case, and you don't
need to keep on writing (1-) during it. Thus

Phi' = phi
phi' = -x phi

therefore
S' = (phi Phi' - Phi phi') / phi^2
= (phi^2 + x Phi phi) / phi^2
= 1 + x S

giving you the recurrence in the terms of the power series for S.

 

[Bill Daly]

If you define

R(x) = cPhi(x)/phi(x)

then the function R is a well behaved function that has
a rapidly convergent Taylor expansion, an expansion
with terms easily determined by a two-step recursion.
(Try to derive the recursion yourself:
r[k+1](x)=x*r[k](x)+k*r[k-1](x),
where r[k](x) is the kth derivative of R(x).)
...
double Phi(double x)
{long double s,t=0,a=1.2533141373155L,z=0,b=-1,pwr=1;
int i;
s=a+b*x;
for(i=2;s!=t;i+=2)
{ a=(a+z*b)/i;
b=(b+z*a)/(i+1);
pwr=pwr*x*x;
t=s;
s+=pwr*(a+x*b);
}
return 1.-s*exp(-.5*x*x-.91893853320467274178L);
}

It is not clear what role z is supposed to play here, since it is
initialized to 0 and never modified. Removing z, and the terms
involving z, corresponds to your definition of r[k](x):

double Phi(double x)
{long double s,t=0,a=1.2533141373155L,b=-1,pwr=1;
int i;
s=a+b*x;
for(i=2;s!=t;i+=2)
{ a/=i;
b/=(i+1);
pwr=pwr*x*x;
t=s;
s+=pwr*(a+x*b);
}
return 1.-s*exp(-.5*x*x-.91893853320467274178L);
}

A possible improvement is to note that the even part of R(x),
(R(x)+R(-x))/2, is equal to 1/phi(x), thus only the odd part of R(x),
(R(x)-R(-x))/2, needs to be computed:

double Phi(double x)
{long double s,t=0,b=1,pwr=x;
int i;
s=x;
for(i=2;s!=t;i+=2)
{ b/=(i+1);
pwr=pwr*x*x;
t=s;
s+=pwr*b;
}
return .5+s*exp(-.5*x*x-.91893853320467274178L);
}

I've been using a similar algorithm, using R(x) = 1/(Phi(x)*Phi(-x)),
which also has a nice Taylor series expansion with a couple of other
desirable properties: it is an even function, thus only the terms in
x^2i need to be computed; and all the coefficients are positive, thus
eliminating the loss of precision resulting from subtracting two large
numbers to get a small answer. This leads to (2*Phi(x)-1)^2 = 1 -
4/R(x), thus Phi(x) = 1/2 +/- sqrt(1/4-1/R(x)), where the sign of the
sqrt should be the same as the sign of x. This has the advantage that
when |x| is large, R(x) is also large, thus 1/R(x) is quite close to 0
and the result is fairly close to the correct value in absolute terms.

 

[bv]

A possible improvement is to note that the even part of R(x),
(R(x)+R(-x))/2, is equal to 1/phi(x), thus only the odd part of R(x),
(R(x)-R(-x))/2, needs to be computed:

double Phi(double x)
{long double s,t=0,b=1,pwr=x;
int i;
s=x;
for(i=2;s!=t;i+=2)
{ b/=(i+1);
pwr=pwr*x*x;
t=s;
s+=pwr*b;
}
return .5+s*exp(-.5*x*x-.91893853320467274178L);
}

Very nice, and a definite improvement*! It speeds up the convergence
considerably which, quite possibly, accounts for more accurate results.
i.e. w/CVF on Wintel,

x phi(x)
0.123 0.5489464510164368
1.200 0.8849303297782917
2.400 0.9918024640754040
6.100 0.9999999994696567
-6.100 0.0000000005303433
-1.100 0.1356660609463827
7.200 0.9999999999996979

* Actually, it's an understatement considering the Syziphian effort at
LANL some odd thirty years ago. See, netlib\fn lib.

 

[Jerry W. Lewis]

Matthias Klaey wrote:
...

Is there someplace a table available of the "true" values of the
normal distribution up to 20 places?

A 20 decimal place table of normal percentiles was published in JASA
65(330):635-638,June 1970. The article is "Tables of Normal
Percentiles" by John S. White. It is downloadable from JSTOR if you
have JSTOR access.

Jerry

 

[Deane Yang]

And what about the inverse cumulative normal distribution?
I used to implement an algorithm you published,
but do not like it because it is machine-dependent.
What in your opinion is the best algorithm?

 

[George Marsaglia]

And what about the inverse cumulative normal distribution?
I used to implement an algorithm you published,
but do not like it because it is machine-dependent.
What in your opinion is the best algorithm?

With the improvement on my method for
evaluating cPhi(x) as R(x)*phi(x),
when the Taylor series for R(x) is about zero,
(suggested by Bill Daly), the simple C function

double Phi(double x)
{ long double t=0,b=1,s=x,pwr=x;
int i;
for(i=2;s!=t;i+=2)
{ b/=(i+1); pwr=pwr*x*x; t=s; s+=pwr*b;}
return .5+s*exp(-.5*x*x-.91893853320467274178L);
}

should serve very well for solving the equation
Phi(X)=U
for positive X, given U>1/2, by Newton's method:
Start with an initial estimate

x=sqrt(-1.6*ln(1-U^2));

then repeat

x=x-(Phi(x)-U)/phi(x);

until you get no further changes in x.


The paper that you refer to must be where I want
to generate a normal X by means of solving
Phi(X)=U, given a random U, uniform in [0,1),
by getting an integer j, formed from certain
bits of the exponent part of the floating point
representation of U. Then use j to access tabled
values A[j] and B[j] to initalize the Taylor series
for Phi inverse, using the fractional part of U
as the argument in the series. It is designed for
speed in generating a normal variate directly
as Phi^(-1)(U), but only provides accuracy to 6/7 digits,
(As I recall, the Taylor series was easy and fast only for the first few
terms.)

It is described in
``Normal (Gaussian) random variables in supercomputers'', (1991)
Journal of Supercomputing}, v 5, 49--55,
where the method is particularly suited for parallel computation.

George Marsaglia

 

[Bill Daly]

Very nice, and a definite improvement*! It speeds up the convergence
considerably which, quite possibly, accounts for more accurate results.
i.e. w/CVF on Wintel,

x phi(x)
0.123 0.5489464510164368
1.200 0.8849303297782917
2.400 0.9918024640754040
6.100 0.9999999994696567
-6.100 0.0000000005303433
-1.100 0.1356660609463827
7.200 0.9999999999996979

* Actually, it's an understatement considering the Syziphian effort at
LANL some odd thirty years ago. See, netlib\fn lib.

It is possible to compress the code further at the expense of clarity
(warning: this is not for the faint of heart):

double Phi(double x)
{long double s=x,t=0,b=x,x2=x*x,i=1;
while(s!=t) s=(t=s)+(b*=x2/(i+=2));
return .5+s*exp(-.5*x2-.91893853320467274178L);
}

This avoids the unnecessary computation of i+1 followed by the
unnecessary conversion of the result from int to long double,
eliminates the unnecessary variable pwr, and precomputes x*x.