lib/ppnd.c

   1 #include "xmath.h"
   2
   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
  15
  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
  22
  23 /*
  24 c
  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.
  30 c
  31
  32 derived from AS111 fortran version
  33 */
  34
  35 double ppnd(p, ifault)
  36         double p;
  37         int *ifault;
  38 {
  39   double q,r,ppn;
  40
  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);
  47   }
  48   else {
  49     r = p;
  50     if(q > zero) r = one - p;
  51     if(r <= zero) {
  52       *ifault = 1;
  53       return(0.0);
  54     }
  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;
  58   }
  59   return(ppn);
  60 }