lib/nor.c

   1 # include <stdio.h>
   2 # include "xmath.h"
   3
   4 #define P10 242.66795523053175
   5 #define P11 21.979261618294152
   6 #define P12 6.9963834886191355
   7 #define P13 -.035609843701815385
   8 #define Q10 215.05887586986120
   9 #define Q11 91.164905404514901
  10 #define Q12 15.082797630407787
  11 #define Q13 1.0
  12
  13 #define P20 300.4592610201616005
  14 #define P21 451.9189537118729422
  15 #define P22 339.3208167343436870
  16 #define P23 152.9892850469404039
  17 #define P24 43.16222722205673530
  18 #define P25 7.211758250883093659
  19 #define P26 .5641955174789739711
  20 #define P27 -.0000001368648573827167067
  21 #define Q20 300.4592609569832933
  22 #define Q21 790.9509253278980272
  23 #define Q22 931.3540948506096211
  24 #define Q23 638.9802644656311665
  25 #define Q24 277.5854447439876434
  26 #define Q25 77.00015293522947295
  27 #define Q26 12.78272731962942351
  28 #define Q27 1.0
  29
  30 #define P30 -.00299610707703542174
  31 #define P31 -.0494730910623250734
  32 #define P32 -.226956593539686930
  33 #define P33 -.278661308609647788
  34 #define P34 -.0223192459734184686
  35 #define Q30 .0106209230528467918
  36 #define Q31 .191308926107829841
  37 #define Q32 1.05167510706793207
  38 #define Q33 1.98733201817135256
  39 #define Q34 1.0
  40
  41 #define SQRT2 1.414213562373095049
  42 #define SQRTPI 1.772453850905516027
  43
  44 void
  45 normbase(double *x, double *phi)
  46 {
  47   int sn;
  48   double R1, R2, y, y2, y3, y4, y5, y6, y7, erf, erfc, z, z2, z3, z4;
  49
  50   y = *x / SQRT2;
  51   if (y < 0) {
  52     y = -y;
  53     sn = -1;
  54   } else {
  55     sn = 1;
  56   }
  57   y2 = y * y;
  58   y4 = y2 * y2;
  59   y6 = y4 * y2;
  60
  61   if(y < 0.46875) {
  62     R1 = P10 + P11 * y2 + P12 * y4 + P13 * y6;
  63     R2 = Q10 + Q11 * y2 + Q12 * y4 + Q13 * y6;
  64     erf = y * R1 / R2;
  65     if(sn == 1) {
  66       *phi = 0.5 + 0.5*erf;
  67     } else {
  68       *phi = 0.5 - 0.5*erf;
  69     }
  70   } else {
  71     if(y < 4.0) {
  72       y3 = y2 * y;
  73       y5 = y4 * y;
  74       y7 = y6 * y;
  75       R1 = P20 + P21 * y + P22 * y2 + P23 * y3 + P24 * y4 + P25 * y5 +
  76         P26 * y6 + P27 * y7;
  77       R2 = Q20 + Q21 * y + Q22 * y2 + Q23 * y3 + Q24 * y4 + Q25 * y5 +
  78         Q26 * y6 + Q27 * y7;
  79       erfc = exp(-y2) * R1 / R2;
  80       if (sn == 1) {
  81         *phi = 1.0 - (0.5 * erfc);
  82       } else {
  83         *phi = 0.5 * erfc;
  84       }
  85     }  else {
  86       z = y4;
  87       z2 = z * z;
  88       z3 = z2 * z;
  89       z4 = z2 * z2;
  90       R1 = P30 + P31 * z + P32 * z2 + P33 * z3 + P34 * z4;
  91       R2 = Q30 + Q31 * z + Q32 * z2 + Q33 * z3 + Q34 * z4;
  92       erfc = (exp(-y2)/y) * (1.0 / SQRTPI + R1 / (R2 * y2));
  93       if(sn == 1) {
  94         *phi = 1.0 - ( 0.5 * erfc);
  95       } else {
  96         *phi = (0.5 * erfc);
  97       }
  98     }
  99   }
 100 }