Re: 'erf' function in C

From: George Marsaglia (geo_at_stat.fsu.edu)
Date: 02/23/04 

Date: Mon, 23 Feb 2004 11:34:53 -0500

"Deane Yang" <deane_yang@yahoo.com> wrote in message
news:36ab6f9c.0402221143.13262bfa@posting.google.com...
> 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

Relevant Pages

  • Re: erf function in C
    ... Would you still recommend this? ... Is there another fast but portable algorithm with sufficient ... > should serve very well for solving the equation ... the Taylor series was easy and fast only for the first few ...
    (comp.lang.c)
  • Re: Algorithm for calculating matrix exponential
    ... I am a programmer and I will like most of all to see an algorithm, ... Taylor series approach is probably a good one. ... diagonalizing the matrix should be considered. ...
    (sci.math)
  • Re: opponents of taylor and lhospital ?
    ... like l'hospitals rule nor taylor series. ... There is in fact an _algorithm_ for evaluating limits of a wide class ... precisely it works by computing higher-rank valuations in certain ...
    (sci.math)