'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]

"James Harlacher" <jp****@hotmail .com> a écrit dans le message de
news:82******** *************** ***@posting.goo gle.com...

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.1283791670955 125738961589031 21545171688
* 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.84506 291151)
* 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.0000000000000 0000000e-01, /* 0x3FE00000, 0x00000000 */
one = 1.0000000000000 0000000e+00, /* 0x3FF00000, 0x00000000 */
two = 2.0000000000000 0000000e+00, /* 0x40000000, 0x00000000 */
/* c = (float)0.845062 91151 */
erx = 8.4506291151046 7529297e-01, /* 0x3FEB0AC1, 0x60000000 */
/*
* Coefficients for approximation to erf on [0,0.84375]
*/
efx = 1.2837916709551 2586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
efx8= 1.0270333367641 0069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
pp0 = 1.2837916709551 2558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
pp1 = -3.2504210724700 1499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
pp2 = -2.8481749575598 5104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
pp3 = -5.7702702964894 4159157e-03, /* 0xBF77A291, 0x236668E4 */
pp4 = -2.3763016656650 1626084e-05, /* 0xBEF8EAD6, 0x120016AC */
qq1 = 3.9791722395915 5352819e-01, /* 0x3FD97779, 0xCDDADC09 */
qq2 = 6.5022249988767 2944485e-02, /* 0x3FB0A54C, 0x5536CEBA */
qq3 = 5.0813062818757 6562776e-03, /* 0x3F74D022, 0xC4D36B0F */
qq4 = 1.3249473800432 1644526e-04, /* 0x3F215DC9, 0x221C1A10 */
qq5 = -3.9602282787753 6812320e-06, /* 0xBED09C43, 0x42A26120 */
/*
* Coefficients for approximation to erf in [0.84375,1.25]
*/
pa0 = -2.3621185607526 5944077e-03, /* 0xBF6359B8, 0xBEF77538 */
pa1 = 4.1485611868374 8331666e-01, /* 0x3FDA8D00, 0xAD92B34D */
pa2 = -3.7220787603570 1323847e-01, /* 0xBFD7D240, 0xFBB8C3F1 */
pa3 = 3.1834661990116 1753674e-01, /* 0x3FD45FCA, 0x805120E4 */
pa4 = -1.1089469428239 6677476e-01, /* 0xBFBC6398, 0x3D3E28EC */
pa5 = 3.5478304325618 2359371e-02, /* 0x3FA22A36, 0x599795EB */
pa6 = -2.1663755948687 9084300e-03, /* 0xBF61BF38, 0x0A96073F */
qa1 = 1.0642088040084 4228286e-01, /* 0x3FBB3E66, 0x18EEE323 */
qa2 = 5.4039791770217 1048937e-01, /* 0x3FE14AF0, 0x92EB6F33 */
qa3 = 7.1828654414196 2662868e-02, /* 0x3FB2635C, 0xD99FE9A7 */
qa4 = 1.2617121980876 1642112e-01, /* 0x3FC02660, 0xE763351F */
qa5 = 1.3637083912029 0507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */
qa6 = 1.1984499846799 1074170e-02, /* 0x3F888B54, 0x5735151D */
/*
* Coefficients for approximation to erfc in [1.25,1/0.35]
*/
ra0 = -9.8649440348471 4822705e-03, /* 0xBF843412, 0x600D6435 */
ra1 = -6.9385857270718 1764372e-01, /* 0xBFE63416, 0xE4BA7360 */
ra2 = -1.0558626225323 2909814e+01, /* 0xC0251E04, 0x41B0E726 */
ra3 = -6.2375332450326 0060396e+01, /* 0xC04F300A, 0xE4CBA38D */
ra4 = -1.6239666946257 3470355e+02, /* 0xC0644CB1, 0x84282266 */
ra5 = -1.8460509290671 1035994e+02, /* 0xC067135C, 0xEBCCABB2 */
ra6 = -8.1287435506306 5934246e+01, /* 0xC0545265, 0x57E4D2F2 */
ra7 = -9.8143293441691 4548592e+00, /* 0xC023A0EF, 0xC69AC25C */
sa1 = 1.9651271667439 2571292e+01, /* 0x4033A6B9, 0xBD707687 */
sa2 = 1.3765775414351 9042600e+02, /* 0x4061350C, 0x526AE721 */
sa3 = 4.3456587747522 9228821e+02, /* 0x407B290D, 0xD58A1A71 */
sa4 = 6.4538727173326 7880336e+02, /* 0x40842B19, 0x21EC2868 */
sa5 = 4.2900814002756 7833386e+02, /* 0x407AD021, 0x57700314 */
sa6 = 1.0863500554177 9435134e+02, /* 0x405B28A3, 0xEE48AE2C */
sa7 = 6.5702497703192 8170135e+00, /* 0x401A47EF, 0x8E484A93 */
sa8 = -6.0424415214858 0987438e-02, /* 0xBFAEEFF2, 0xEE749A62 */
/*
* Coefficients for approximation to erfc in [1/.35,28]
*/
rb0 = -9.8649429247000 9928597e-03, /* 0xBF843412, 0x39E86F4A */
rb1 = -7.9928323768052 3006574e-01, /* 0xBFE993BA, 0x70C285DE */
rb2 = -1.7757954917754 7519889e+01, /* 0xC031C209, 0x555F995A */
rb3 = -1.6063638485582 1916062e+02, /* 0xC064145D, 0x43C5ED98 */
rb4 = -6.3756644336838 9627722e+02, /* 0xC083EC88, 0x1375F228 */
rb5 = -1.0250951316110 7724954e+03, /* 0xC0900461, 0x6A2E5992 */
rb6 = -4.8351919160865 1397019e+02, /* 0xC07E384E, 0x9BDC383F */
sb1 = 3.0338060743482 4582924e+01, /* 0x403E568B, 0x261D5190 */
sb2 = 3.2579251299657 3918826e+02, /* 0x40745CAE, 0x221B9F0A */
sb3 = 1.5367295860844 3695994e+03, /* 0x409802EB, 0x189D5118 */
sb4 = 3.1998582195085 9553908e+03, /* 0x40A8FFB7, 0x688C246A */
sb5 = 2.5530504064331 6442583e+03, /* 0x40A3F219, 0xCEDF3BE6 */
sb6 = 4.7452854120695 5367215e+02, /* 0x407DA874, 0xE79FE763 */
sb7 = -2.2440952446585 8183362e+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>=0x7ff000 00) { /* 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+ef x8*x); /*avoid underflow */
return x + efx*x;
}
z = x*x;
r = pp0+z*(pp1+z*(p p2+z*(pp3+z*pp4 )));
s = one+z*(qq1+z*(q q2+z*(qq3+z*(qq 4+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*(p a2+s*(pa3+s*(pa 4+s*(pa5+s*pa6) ))));
Q = one+s*(qa1+s*(q a2+s*(qa3+s*(qa 4+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*ra 7))))));
S=one+s*(sa1+s* (sa2+s*(sa3+s*( sa4+s*(
sa5+s*(sa6+s*(s a7+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*sb 7))))));
}
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>=0x7ff000 00) { /* erfc(nan)=nan */
/* erfc(+-inf)=0,2 */
return (double)(((unsi gned)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*(p p2+z*(pp3+z*pp4 )));
s = one+z*(qq1+z*(q q2+z*(qq3+z*(qq 4+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*(p a2+s*(pa3+s*(pa 4+s*(pa5+s*pa6) ))));
Q = one+s*(qa1+s*(q a2+s*(qa3+s*(qa 4+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*ra 7))))));
S=one+s*(sa1+s* (sa2+s*(sa3+s*( sa4+s*(
sa5+s*(sa6+s*(s a7+s*sa8))))))) ;
} else { /* |x| >= 1/.35 ~ 2.857143 */
if(hx<0&&ix>=0x 40180000) 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*sb 7))))));
}
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]

"James Harlacher" <jp****@hotmail .com> wrote in message
news:82******** *************** ***@posting.goo gle.com...

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]

James Harlacher wrote:

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]

"James Harlacher" <jp****@hotmail .com> wrote in message
news:82******** *************** ***@posting.goo gle.com...

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.25331 41373155L,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-.91893853320467 274178L);
}

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.5489464510164 369
.54894645101643 675909
1.2 0.8849303297782 918
.88493032977829 173198
2.4 0.9918024640754 040
.99180246407540 387055
6.1 0.9999999994696 578
.99999999946965 767370
-6.1 0.0000000005303 427
.00000000053034 232638
-1.1 0.1356660609463 828
.13566606094638 267517
7.2 0.9999999999996 990
.99999999999969 893721
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]

On Mon, 16 Feb 2004 11:13:48 -0500, "George Marsaglia"
<ge*@stat.fsu.e du> wrote:

"James Harlacher" <jp****@hotmail .com> wrote in message
news:82******* *************** ****@posting.go ogle.com...

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.5489464510164 369
.54894645101643 675909
1.2 0.8849303297782 918
.88493032977829 173198
2.4 0.9918024640754 040
.99180246407540 387055
6.1 0.9999999994696 578
.99999999946965 767370
-6.1 0.0000000005303 427
.00000000053034 232638
-1.1 0.1356660609463 828
.13566606094638 267517
7.2 0.9999999999996 990
.99999999999969 893721
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

[...]

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

With kind regards
Matthias Kläy
--
www.kcc.ch

[P.J. Plauger]

"George Marsaglia" <ge*@stat.fsu.e du> wrote in message
news:7J******** ************@co mcast.com...

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]

George Marsaglia wrote:

"James Harlacher" <jp****@hotmail .com> wrote in message
news:82******** *************** ***@posting.goo gle.com...

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.

.... 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]

Matthias Klaey wrote:

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.wolfr am.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]

Sidney Cadot <si****@jigsaw. nl> writes:

Matthias Klaey wrote:

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.wolfr am.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=Standar dForm" part of the URL to
"format=OutputF orm". 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