src/c/nor.c

   1 /* Double precision version of routines in ERF from the netlib SPECFUN */
   2 /* library. Translated by f2c and modified. */
   3
   4 #include "xlisp.h"
   5 #include "xlstat.h"
   6 #include "linalg.h"
   7
   8 static VOID calerf P3C(double, arg, double *, result, int, jint)
   9 {
  10   /* ------------------------------------------------------------------ */
  11   /* This packet evaluates  erf(x),  erfc(x),  and  exp(x*x)*erfc(x) */
  12   /*   for a real argument  x.  It contains three FUNCTION type */
  13   /*   subprograms: ERF, ERFC, and ERFCX (or DERF, DERFC, and DERFCX), */
  14   /*   and one SUBROUTINE type subprogram, CALERF.  The calling */
  15   /*   statements for the primary entries are: */
  16
  17   /*                   Y=ERF(X)     (or   Y=DERF(X)), */
  18
  19   /*                   Y=ERFC(X)    (or   Y=DERFC(X)), */
  20   /*   and */
  21   /*                   Y=ERFCX(X)   (or   Y=DERFCX(X)). */
  22
  23   /*   The routine  CALERF  is intended for internal packet use only, */
  24   /*   all computations within the packet being concentrated in this */
  25   /*   routine.  The function subprograms invoke  CALERF  with the */
  26   /*   statement */
  27
  28   /*          CALL CALERF(ARG,RESULT,JINT) */
  29
  30   /*   where the parameter usage is as follows */
  31
  32   /*      Function                     Parameters for CALERF */
  33   /*       call              ARG                  Result          JINT */
  34
  35   /*     ERF(ARG)      ANY REAL ARGUMENT         ERF(ARG)          0 */
  36   /*     ERFC(ARG)     ABS(ARG) .LT. XBIG        ERFC(ARG)         1 */
  37   /*     ERFCX(ARG)    XNEG .LT. ARG .LT. XMAX   ERFCX(ARG)        2 */
  38
  39   /*   The main computation evaluates near-minimax approximations */
  40   /*   from "Rational Chebyshev approximations for the error function" */
  41   /*   by W. J. Cody, Math. Comp., 1969, PP. 631-638.  This */
  42   /*   transportable program uses rational functions that theoretically */
  43   /*   approximate  erf(x)  and  erfc(x)  to at least 18 significant */
  44   /*   decimal digits.  The accuracy achieved depends on the arithmetic */
  45   /*   system, the compiler, the intrinsic functions, and proper */
  46   /*   selection of the machine-dependent constants. */
  47
  48   /* ******************************************************************* */
  49   /* ******************************************************************* */
  50
  51   /* Explanation of machine-dependent constants */
  52
  53   /*   XMIN   = the smallest positive floating-point number. */
  54   /*   XINF   = the largest positive finite floating-point number. */
  55   /*   XNEG   = the largest negative argument acceptable to ERFCX; */
  56   /*            the negative of the solution to the equation */
  57   /*            2*exp(x*x) = XINF. */
  58   /*   XSMALL = argument below which erf(x) may be represented by */
  59   /*            2*x/sqrt(pi)  and above which  x*x  will not underflow. */
  60   /*            A conservative value is the largest machine number X */
  61   /*            such that   1.0 + X = 1.0   to machine precision. */
  62   /*   XBIG   = largest argument acceptable to ERFC;  solution to */
  63   /*            the equation:  W(x) * (1-0.5/x**2) = XMIN,  where */
  64   /*            W(x) = exp(-x*x)/[x*sqrt(pi)]. */
  65   /*   XHUGE  = argument above which  1.0 - 1/(2*x*x) = 1.0  to */
  66   /*            machine precision.  A conservative value is */
  67   /*            1/[2*sqrt(XSMALL)] */
  68   /*   XMAX   = largest acceptable argument to ERFCX; the minimum */
  69   /*            of XINF and 1/[sqrt(pi)*XMIN]. */
  70
  71   /*   Approximate values for some important machines are: */
  72
  73   /*                          XMIN       XINF        XNEG     XSMALL */
  74
  75   /*  CDC 7600      (S.P.)  3.13E-294   1.26E+322   -27.220  7.11E-15 */
  76   /*  CRAY-1        (S.P.)  4.58E-2467  5.45E+2465  -75.345  7.11E-15 */
  77   /*  IEEE (IBM/XT, */
  78   /*    SUN, etc.)  (S.P.)  1.18E-38    3.40E+38     -9.382  5.96E-8 */
  79   /*  IEEE (IBM/XT, */
  80   /*    SUN, etc.)  (D.P.)  2.23D-308   1.79D+308   -26.628  1.11D-16 */
  81   /*  IBM 195       (D.P.)  5.40D-79    7.23E+75    -13.190  1.39D-17 */
  82   /*  UNIVAC 1108   (D.P.)  2.78D-309   8.98D+307   -26.615  1.73D-18 */
  83   /*  VAX D-Format  (D.P.)  2.94D-39    1.70D+38     -9.345  1.39D-17 */
  84   /*  VAX G-Format  (D.P.)  5.56D-309   8.98D+307   -26.615  1.11D-16 */
  85
  86
  87   /*                          XBIG       XHUGE       XMAX */
  88
  89   /*  CDC 7600      (S.P.)  25.922      8.39E+6     1.80X+293 */
  90   /*  CRAY-1        (S.P.)  75.326      8.39E+6     5.45E+2465 */
  91   /*  IEEE (IBM/XT, */
  92   /*    SUN, etc.)  (S.P.)   9.194      2.90E+3     4.79E+37 */
  93   /*  IEEE (IBM/XT, */
  94   /*    SUN, etc.)  (D.P.)  26.543      6.71D+7     2.53D+307 */
  95   /*  IBM 195       (D.P.)  13.306      1.90D+8     7.23E+75 */
  96   /*  UNIVAC 1108   (D.P.)  26.582      5.37D+8     8.98D+307 */
  97   /*  VAX D-Format  (D.P.)   9.269      1.90D+8     1.70D+38 */
  98   /*  VAX G-Format  (D.P.)  26.569      6.71D+7     8.98D+307 */
  99
 100   /* ******************************************************************* */
 101   /* ******************************************************************* */
 102
 103   /* Error returns */
 104
 105   /*  The program returns  ERFC = 0      for  ARG .GE. XBIG; */
 106
 107   /*                       ERFCX = XINF  for  ARG .LT. XNEG; */
 108   /*      and */
 109   /*                       ERFCX = 0     for  ARG .GE. XMAX. */
 110
 111
 112   /* Intrinsic functions required are: */
 113
 114   /*     ABS, AINT, EXP */
 115
 116
 117   /*  Author: W. J. Cody */
 118   /*          Mathematics and Computer Science Division */
 119   /*          Argonne National Laboratory */
 120   /*          Argonne, IL 60439 */
 121
 122   /*  Latest modification: March 19, 1990 */
 123
 124   /* ------------------------------------------------------------------ */
 125   double xden, xnum;
 126   int i;
 127   double x, y, del, ysq;
 128   /* ------------------------------------------------------------------ */
 129   /*  Mathematical constants */
 130   /* ------------------------------------------------------------------ */
 131   static double four = 4.;
 132   static double one = 1.;
 133   static double half = .5;
 134   static double two = 2.;
 135   static double zero = 0.;
 136   static double sqrpi = .56418958354775628695;
 137   static double thresh = .46875;
 138   static double sixten = 16.;
 139   /* ------------------------------------------------------------------ */
 140   /*  Machine-dependent constants */
 141   /* ------------------------------------------------------------------ */
 142 #ifdef IEEEFP
 143   static double xinf = 1.79e308;
 144   static double xneg = -26.628;
 145   static double xsmall = 1.11e-16;
 146   static double xbig = 26.543;
 147   static double xhuge = 6.71e7;
 148   static double xmax = 2.53e307;
 149 #else
 150 #ifdef CRAYCC
 151   static double xinf = 5.45e2465;
 152   static double xneg = -75.345;
 153   static double xsmall = 7.11e-15;
 154   static double xbig = 75.326;
 155   static double xhuge = 8.39e6;
 156   static double xmax = 5.45e2465;
 157 #else /* use IBM 196 values */
 158   static double xinf = 7.23e75;
 159   static double xneg = -13.190;
 160   static double xsmall = 1.39e-17;
 161   static double xbig = 13.306;
 162   static double xhuge = 1.90e8;
 163   static double xmax = 7.23e75;
 164 #endif /* CRAYCC */
 165 #endif /* IEEEFP */
 166   /* ------------------------------------------------------------------ */
 167   /*  Coefficients for approximation to  erf  in first interval */
 168   /* ------------------------------------------------------------------ */
 169   static double a[5] = { 3.1611237438705656,113.864154151050156,
 170                            377.485237685302021,3209.37758913846947,
 171                            .185777706184603153 };
 172   static double b[4] = { 23.6012909523441209,244.024637934444173,
 173                            1282.61652607737228,2844.23683343917062 };
 174   /* ------------------------------------------------------------------ */
 175   /*  Coefficients for approximation to  erfc  in second interval */
 176   /* ------------------------------------------------------------------ */
 177   static double c[9] = { .564188496988670089,8.88314979438837594,
 178                            66.1191906371416295,298.635138197400131,
 179                            881.95222124176909,1712.04761263407058,
 180                            2051.07837782607147,1230.33935479799725,
 181                            2.15311535474403846e-8 };
 182   static double d[8] = { 15.7449261107098347,117.693950891312499,
 183                            537.181101862009858,1621.38957456669019,
 184                            3290.79923573345963,4362.61909014324716,
 185                            3439.36767414372164,1230.33935480374942 };
 186   /* ------------------------------------------------------------------ */
 187   /*  Coefficients for approximation to  erfc  in third interval */
 188   /* ------------------------------------------------------------------ */
 189   static double p[6] = { .305326634961232344,.360344899949804439,
 190                            .125781726111229246,.0160837851487422766,
 191                            6.58749161529837803e-4,.0163153871373020978 };
 192   static double q[5] = { 2.56852019228982242,1.87295284992346047,
 193                            .527905102951428412,.0605183413124413191,
 194                            .00233520497626869185 };
 195   /* ------------------------------------------------------------------ */
 196   x = arg;
 197   y = abs(x);
 198   if (y <= thresh) {
 199     /* ------------------------------------------------------------------ */
 200     /*  Evaluate  erf  for  |X| <= 0.46875 */
 201     /* ------------------------------------------------------------------ */
 202     ysq = zero;
 203     if (y > xsmall) {
 204       ysq = y * y;
 205     }
 206     xnum = a[4] * ysq;
 207     xden = ysq;
 208     for (i = 1; i <= 3; ++i) {
 209       xnum = (xnum + a[i - 1]) * ysq;
 210       xden = (xden + b[i - 1]) * ysq;
 211     }
 212     *result = x * (xnum + a[3]) / (xden + b[3]);
 213     if (jint != 0) {
 214       *result = one - *result;
 215     }
 216     if (jint == 2) {
 217       *result = exp(ysq) * *result;
 218     }
 219     return;
 220     /* ------------------------------------------------------------------ */
 221     /*  Evaluate  erfc  for 0.46875 <= |X| <= 4.0 */
 222     /* ------------------------------------------------------------------ */
 223   }
 224   else if (y <= four) {
 225     xnum = c[8] * y;
 226     xden = y;
 227     for (i = 1; i <= 7; ++i) {
 228       xnum = (xnum + c[i - 1]) * y;
 229       xden = (xden + d[i - 1]) * y;
 230     }
 231     *result = (xnum + c[7]) / (xden + d[7]);
 232     if (jint != 2) {
 233       ysq = floor(y * sixten) / sixten;
 234       del = (y - ysq) * (y + ysq);
 235       *result = exp(-ysq * ysq) * exp(-del) * *result;
 236     }
 237     /* ------------------------------------------------------------------ */
 238     /*  Evaluate  erfc  for |X| > 4.0 */
 239     /* ------------------------------------------------------------------ */
 240   }
 241   else {
 242     *result = zero;
 243     if (y >= xbig) {
 244       if (jint != 2 || y >= xmax) {
 245         goto L300;
 246       }
 247       if (y >= xhuge) {
 248         *result = sqrpi / y;
 249         goto L300;
 250       }
 251     }
 252     ysq = one / (y * y);
 253     xnum = p[5] * ysq;
 254     xden = ysq;
 255     for (i = 1; i <= 4; ++i) {
 256       xnum = (xnum + p[i - 1]) * ysq;
 257       xden = (xden + q[i - 1]) * ysq;
 258     }
 259     *result = ysq * (xnum + p[4]) / (xden + q[4]);
 260     *result = (sqrpi - *result) / y;
 261     if (jint != 2) {
 262       ysq = floor(y * sixten) / sixten;
 263       del = (y - ysq) * (y + ysq);
 264       *result = exp(-ysq * ysq) * exp(-del) * *result;
 265     }
 266   }
 267   /* ------------------------------------------------------------------ */
 268   /*  Fix up for negative argument, erf, etc. */
 269   /* ------------------------------------------------------------------ */
 270  L300:
 271   if (jint == 0) {
 272     *result = half - *result + half;
 273     if (x < zero) {
 274       *result = -(*result);
 275     }
 276   }
 277   else if (jint == 1) {
 278     if (x < zero) {
 279       *result = two - *result;
 280     }
 281   }
 282   else {
 283     if (x < zero) {
 284       if (x < xneg) {
 285         *result = xinf;
 286       }
 287       else {
 288         ysq = -floor(-(x * sixten)) / sixten;
 289         del = (x - ysq) * (x + ysq);
 290         y = exp(ysq * ysq) * exp(del);
 291         *result = y + y - *result;
 292       }
 293     }
 294   }
 295   return;
 296 }
 297
 298 #ifdef DEBUG
 299 double derf_ P1C(double *, x)
 300 {
 301   int jint;
 302   double result;
 303
 304   /* This subprogram computes approximate values for erf(x). */
 305   /*   (see comments heading CALERF). */
 306   /*   Author/date: W. J. Cody, January 8, 1985 */
 307
 308   jint = 0;
 309   calerf(*x, &result, jint);
 310   return result;
 311 }
 312
 313 double derfc_ P1C(double *, x)
 314 {
 315   int jint;
 316   double result;
 317
 318   /* This subprogram computes approximate values for erfc(x). */
 319   /*   (see comments heading CALERF). */
 320   /*   Author/date: W. J. Cody, January 8, 1985 */
 321
 322   jint = 1;
 323   calerf(*x, &result, jint);
 324   return result;
 325 }
 326
 327 double derfcx_ (double *, x)
 328 {
 329   int jint;
 330   double result;
 331
 332   /* This subprogram computes approximate values for exp(x*x) * erfc(x). */
 333   /*   (see comments heading CALERF). */
 334   /*   Author/date: W. J. Cody, March 30, 1987 */
 335
 336   jint = 2;
 337   calerf(*x, &result, jint);
 338   return result;
 339 }
 340 #endif /* DEBUG */
 341
 342 #define SQRT2 1.414213562373095049
 343
 344 VOID normbase P2C(double *, x, double *, phi)
 345 {
 346   double y;
 347
 348   y = *x / SQRT2;
 349   if (y < 0) {
 350     calerf(-y, phi, 1);
 351     *phi = 0.5 * *phi;
 352   }
 353   else {
 354     calerf(y, phi, 1);
 355     *phi = 1.0 - 0.5 * *phi;
 356   }
 357 }