blob f8f218a359d59c36a1059137d664c77657fd11db
1 #include "xmath.h"
3 #define zero 0.0
4 #define half 0.5
5 #define one 1.0
6 #define split 0.42e0
7 #define a0 2.50662823884e0
8 #define a1 -18.61500062529e0
9 #define a2 41.39119773534e0
10 #define a3 -25.44106049637e0
11 #define b1 -8.47351093090e0
12 #define b2 23.08336743743e0
13 #define b3 -21.06224101826e0
14 #define b4 3.13082909833e0
16 #define c0 -2.78718931138e0
17 #define c1 -2.29796479134e0
18 #define c2 4.85014127135e0
19 #define c3 2.32121276858e0
20 #define d1 3.54388924762e0
21 #define d2 1.63706781897e0
25 c Algorithm as 111 Applied statistics (1977), vol 26 no 1 page 121
26 c Produces normal deviate corresponding to lower tail area of p
27 c the hash sums are the sums of the moduli of the coefficients
28 c they nave no inherent meanings but are incuded for use in
29 c checking transcriptions. Functions abs,alog and sqrt are used.
32 derived from AS111 fortran version
35 double ppnd(p, ifault)
36 double p;
37 int *ifault;
39 double q,r,ppn;
41 *ifault = 0;
42 q = p - half;
43 if( fabs(q) <= split) {
44 r = q*q;
45 ppn = q * (((a3 * r + a2) * r + a1) * r + a0)
46 / ((((b4 * r + b3) * r + b2) * r + b1) * r + one);
48 else {
49 r = p;
50 if(q > zero) r = one - p;
51 if(r <= zero) {
52 *ifault = 1;
53 return(0.0);
55 r = sqrt(-log(r));
56 ppn = (((c3*r+c2)*r + c1) * r + c0) / ((d2 * r + d1) * r + one);
57 if( q < zero) ppn = -ppn;
59 return(ppn);