Rapid evaluation of the inverse of the normal distribution function (Q1324541)

This is an interesting article with direct application in generating normal random variable by computer programs. The suggested applications are related to Monte Carlo simulation based on massively parallel systems or supercomputers. The idea is to replace larger programs with complicated computations and with difficulties in accuracy controlling by simpler arithmetic programs that use tabled constants. These seem to be the normal evolution since memory becomes cheaper and cheaper. The authors compute the inverse of the cPhi function \[ cPhi(x) = (2/\pi)^{1/2} \int^ \infty_ x \exp (-t^ 2/2) dt=u, \] using a uniform random variable as input and the truncated Taylor series development of it. In order to increase the speed the coefficients of the truncated Taylor series \[ x(u_ 0+h) = x(u_ 0) + x'(u_ 0) \cdot h + {1 \over 2} x''(u_ 0) \cdot h^ 2 + {1 \over 6} x'''(u_ 0) \cdot h^ 3, \] are predetermined for 1024 points. And here comes another bright idea: the 1024 points are chosen based on the representation of the uniform random variable in modern computers as floating point variable of the form: \(u=2^{-k} ((1/2) + (j/64)) + 2^{-k} \cdot (m/2^{24})\) with \(0 \leq k<32\), \(0 \leq j<32\) and \(0 \leq m<2^{18}\) and considering 32 bit representation. With this assumptions and the truncation to the third power of \(h\) of the Taylor series, the authors show that the error does not exceed the limit of single precision accuracy. Furthermore the calculations are speeded up based on reducing multiplications. A number of FORTRAN programs are also presented in order to evaluate the complementary normal distribution function cPhi (several versions) with great accuracy, create the constant tables, and generate the normal distribution variable. These simple programs give the user the possibility to completely control the accuracy.