[BACK]Return to specfun.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / gnuplot

Diff for /OpenXM_contrib/gnuplot/Attic/specfun.c between version 1.1.1.2 and 1.1.1.3

version 1.1.1.2, 2000/01/22 14:16:01 version 1.1.1.3, 2003/09/15 07:09:26
Line 2 
Line 2 
 static char *RCSid = "$Id$";  static char *RCSid = "$Id$";
 #endif  #endif
   
   
 /* GNUPLOT - specfun.c */  /* GNUPLOT - specfun.c */
   
 /*[  /*[
Line 43  static char *RCSid = "$Id$";
Line 42  static char *RCSid = "$Id$";
  *   *
  */   */
   
   
 #include "plot.h"  #include "plot.h"
 #include "fnproto.h"  #include "fnproto.h"
   
   
 extern struct value stack[STACK_DEPTH];  extern struct value stack[STACK_DEPTH];
 extern int s_p;  extern int s_p;
 extern double zero;  extern double zero;
Line 107  static long Xa1 = 40014L;
Line 106  static long Xa1 = 40014L;
 static long Xa2 = 40692L;  static long Xa2 = 40692L;
   
 /* Local function declarations, not visible outside this file */  /* Local function declarations, not visible outside this file */
   static int mtherr(char *name, int code);
   static double polevl(double x, double coef[], int N);
   static double p1evl(double x, double coef[], int N);
 static double confrac __PROTO((double a, double b, double x));  static double confrac __PROTO((double a, double b, double x));
 static double ibeta __PROTO((double a, double b, double x));  static double ibeta __PROTO((double a, double b, double x));
 static double igamma __PROTO((double a, double x));  static double igamma __PROTO((double a, double x));
Line 114  static double ranf __PROTO((double init));
Line 116  static double ranf __PROTO((double init));
 static double inverse_normal_func __PROTO((double p));  static double inverse_normal_func __PROTO((double p));
 static double inverse_error_func __PROTO((double p));  static double inverse_error_func __PROTO((double p));
   
 #ifndef GAMMA  /* UNKnown arithmetic, invokes coefficients given in
 /* Provide GAMMA function for those who do not already have one */   * normal decimal format.  Beware of range boundary
 static double lngamma __PROTO((double z));   * problems (MACHEP, MAXLOG, etc. in const.c) and
 static double lgamneg __PROTO((double z));   * roundoff problems in pow.c:
 static double lgampos __PROTO((double z));   * (Sun SPARCstation)
    */
   #define UNK 1
   
 /**  /* If you define UNK, then be sure to set BIGENDIAN properly. */
  * from statlib, Thu Jan 23 15:02:27 EST 1992 ***  #ifdef FLOAT_WORDS_BIGENDIAN
  *  #define BIGENDIAN 1
  * This file contains two algorithms for the logarithm of the gamma function.  #else
  * Algorithm AS 245 is the faster (but longer) and gives an accuracy of about  #define BIGENDIAN 0
  * 10-12 significant decimal digits except for small regions around X = 1 and  #endif
  * X = 2, where the function goes to zero.  
  * The second algorithm is not part of the AS algorithms.   It is slower but  /* Define to support tiny denormal numbers, else undefine. */
  * gives 14 or more significant decimal digits accuracy, except around X = 1  #define DENORMAL 1
  * and X = 2.   The Lanczos series from which this algorithm is derived is  
  * interesting in that it is a convergent series approximation for the gamma  /* Define to ask for infinity support, else undefine. */
  * function, whereas the familiar series due to De Moivre (and usually wrongly  #define INFINITIES 1
  * called Stirling's approximation) is only an asymptotic approximation, as  
  * is the true and preferable approximation due to Stirling.  /* Define to ask for support of numbers that are Not-a-Number,
  *     else undefine.  This may automatically define INFINITIES in some files. */
  * Uses Lanczos-type approximation to ln(gamma) for z > 0. Reference: Lanczos,  #define NANS 1
  * C. 'A precision approximation of the gamma function', J. SIAM Numer.  
  * Anal., B, 1, 86-96, 1964. Accuracy: About 14 significant digits except for  /* Define to distinguish between -0.0 and +0.0.  */
  * small regions in the vicinity of 1 and 2.  #define MINUSZERO 1
  *  
  * Programmer: Alan Miller CSIRO Division of Mathematics & Statistics  /* Variable for error reporting.  See mtherr.c.  */
  *  extern int      merror;
  * Latest revision - 17 April 1988  
  *  /*
  * Additions: Translated from fortran to C, code added to handle values z < 0.  Cephes Math Library Release 2.0:  April, 1987
  * The global variable signgam contains the sign of the gamma function.  Copyright 1984, 1987 by Stephen L. Moshier
  *  Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  * IMPORTANT: The signgam variable contains garbage until AFTER the call to  */
  * lngamma().  
  *  #include <stdio.h>
  * Permission granted to distribute freely for non-commercial purposes only  
  * Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl  int             merror = 0;
   
   /* Notice: the order of appearance of the following
    * messages is bound to the error codes defined
    * in mconf.h.
  */   */
   static char    *ermsg[7] = {
       "unknown",                  /* error code 0 */
       "domain",                   /* error code 1 */
       "singularity",              /* et seq.      */
       "overflow",
       "underflow",
       "total loss of precision",
       "partial loss of precision"
   };
   
 /* Local data, not visible outside this file  #ifndef DOMAIN
    static double   a[] =  /* HBB 20021103: copied this from Cephes mconf.h */
    {  # define DOMAIN         1       /* argument domain error */
    0.9999999999995183E+00,  # define SING           2       /* argument singularity */
    0.6765203681218835E+03,  # define OVERFLOW       3       /* overflow range error */
    -.1259139216722289E+04,  # define UNDERFLOW      4       /* underflow range error */
    0.7713234287757674E+03,  # define TLOSS          5       /* total loss of precision */
    -.1766150291498386E+03,  # define PLOSS          6       /* partial loss of precision */
    0.1250734324009056E+02,  #endif
    -.1385710331296526E+00,  
    0.9934937113930748E-05,  
    0.1659470187408462E-06,  
    };   */  
   
 /* from Ray Toy */  
 static double GPFAR a[] =  
   static int mtherr(char *name, int code)
 {  {
     .99999999999980993227684700473478296744476168282198,  
     676.52036812188509856700919044401903816411251975244084,  
     -1259.13921672240287047156078755282840836424300664868028,  
     771.32342877765307884865282588943070775227268469602500,  
     -176.61502916214059906584551353999392943274507608117860,  
     12.50734327868690481445893685327104972970563021816420,  
     -.13857109526572011689554706984971501358032683492780,  
     .00000998436957801957085956266828104544089848531228,  
     .00000015056327351493115583383579667028994545044040,  
 };  
   
 static double lgamneg(z)  /* Display string passed by calling program,
 double z;   * which is supposed to be the name of the
    * function in which the error occurred:
    */
       printf("\n%s ", name);
   
   /* Set global error message word */
       merror = code;
   
   /* Display error message defined
    * by the code argument.
    */
       if ((code <= 0) || (code >= 7))
           code = 0;
       printf("%s error\n", ermsg[code]);
   
   /* Return to calling
    * program
    */
       return (0);
   }
   
   /*                                                      polevl.c
    *                                                      p1evl.c
    *
    *      Evaluate polynomial
    *
    *
    *
    * SYNOPSIS:
    *
    * int N;
    * double x, y, coef[N+1], polevl[];
    *
    * y = polevl( x, coef, N );
    *
    *
    *
    * DESCRIPTION:
    *
    * Evaluates polynomial of degree N:
    *
    *                     2          N
    * y  =  C  + C x + C x  +...+ C x
    *        0    1     2          N
    *
    * Coefficients are stored in reverse order:
    *
    * coef[0] = C  , ..., coef[N] = C  .
    *            N                   0
    *
    *  The function p1evl() assumes that coef[N] = 1.0 and is
    * omitted from the array.  Its calling arguments are
    * otherwise the same as polevl().
    *
    *
    * SPEED:
    *
    * In the interest of speed, there are no checks for out
    * of bounds arithmetic.  This routine is used by most of
    * the functions in the library.  Depending on available
    * equipment features, the user may wish to rewrite the
    * program in microcode or assembly language.
    *
    */
   
   /*
   Cephes Math Library Release 2.1:  December, 1988
   Copyright 1984, 1987, 1988 by Stephen L. Moshier
   Direct inquiries to 30 Frost Street, Cambridge, MA 02140
   */
   
   static double polevl(double x, double coef[], int N)
 {  {
     double tmp;      double          ans;
       int             i;
       double         *p;
   
     /* Use reflection formula, then call lgampos() */      p = coef;
     tmp = sin(z * PI);      ans = *p++;
       i = N;
   
     if (fabs(tmp) < MACHEPS) {      do
         tmp = 0.0;          ans = ans * x + *p++;
     } else if (tmp < 0.0) {      while (--i);
         tmp = -tmp;  
         signgam = -1;  
     }  
     return LNPI - lgampos(1.0 - z) - log(tmp);  
   
       return (ans);
 }  }
   
 static double lgampos(z)  /*                                                      p1evl() */
 double z;  /*                                          N
    * Evaluate polynomial when coefficient of x  is 1.0.
    * Otherwise same as polevl.
    */
   
   static double p1evl(double x, double coef[], int N)
 {  {
     double sum;      double          ans;
     double tmp;      double         *p;
     int i;      int             i;
   
     sum = a[0];      p = coef;
     for (i = 1, tmp = z; i < 9; i++) {      ans = x + *p++;
         sum += a[i] / tmp;      i = N - 1;
         tmp++;  
     }  
   
     return log(sum) + LNSQRT2PI - z - 6.5 + (z - 0.5) * log(z + 6.5);      do
           ans = ans * x + *p++;
       while (--i);
   
       return (ans);
 }  }
   
 static double lngamma(z)  #ifndef GAMMA
 double z;  /* Provide GAMMA function for those who do not already have one */
   static double lngamma __PROTO((double z));
   
   int             sgngam;
   /* A[]: Stirling's formula expansion of log gamma
    * B[], C[]: log gamma function between 2 and 3
    */
   #ifdef UNK
   static double   A[] = {
       8.11614167470508450300E-4,
       -5.95061904284301438324E-4,
       7.93650340457716943945E-4,
       -2.77777777730099687205E-3,
       8.33333333333331927722E-2
   };
   static double   B[] = {
       -1.37825152569120859100E3,
       -3.88016315134637840924E4,
       -3.31612992738871184744E5,
       -1.16237097492762307383E6,
       -1.72173700820839662146E6,
       -8.53555664245765465627E5
   };
   static double   C[] = {
   /* 1.00000000000000000000E0, */
       -3.51815701436523470549E2,
       -1.70642106651881159223E4,
       -2.20528590553854454839E5,
       -1.13933444367982507207E6,
       -2.53252307177582951285E6,
       -2.01889141433532773231E6
   };
   /* log( sqrt( 2*pi ) ) */
   static double   LS2PI = 0.91893853320467274178;
   #define MAXLGM 2.556348e305
   #endif
   #ifdef DEC
   static unsigned short A[] = {
       0035524, 0141201, 0034633, 0031405,
       0135433, 0176755, 0126007, 0045030,
       0035520, 0006371, 0003342, 0172730,
       0136066, 0005540, 0132605, 0026407,
       0037252, 0125252, 0125252, 0125132
   };
   static unsigned short B[] = {
       0142654, 0044014, 0077633, 0035410,
       0144027, 0110641, 0125335, 0144760,
       0144641, 0165637, 0142204, 0047447,
       0145215, 0162027, 0146246, 0155211,
       0145322, 0026110, 0010317, 0110130,
       0145120, 0061472, 0120300, 0025363
   };
   static unsigned short C[] = {
   /*0040200,0000000,0000000,0000000*/
       0142257, 0164150, 0163630, 0112622,
       0143605, 0050153, 0156116, 0135272,
       0144527, 0056045, 0145642, 0062332,
       0145213, 0012063, 0106250, 0001025,
       0145432, 0111254, 0044577, 0115142,
       0145366, 0071133, 0050217, 0005122
   };
   /* log( sqrt( 2*pi ) ) */
   static unsigned short LS2P[] = {040153, 037616, 041445, 0172645,};
   #define LS2PI *(double *)LS2P
   #define MAXLGM 2.035093e36
   #endif
   
   #ifdef IBMPC
   static unsigned short A[] = {
       0x6661, 0x2733, 0x9850, 0x3f4a,
       0xe943, 0xb580, 0x7fbd, 0xbf43,
       0x5ebb, 0x20dc, 0x019f, 0x3f4a,
       0xa5a1, 0x16b0, 0xc16c, 0xbf66,
       0x554b, 0x5555, 0x5555, 0x3fb5
   };
   static unsigned short B[] = {
       0x6761, 0x8ff3, 0x8901, 0xc095,
       0xb93e, 0x355b, 0xf234, 0xc0e2,
       0x89e5, 0xf890, 0x3d73, 0xc114,
       0xdb51, 0xf994, 0xbc82, 0xc131,
       0xf20b, 0x0219, 0x4589, 0xc13a,
       0x055e, 0x5418, 0x0c67, 0xc12a
   };
   static unsigned short C[] = {
   /*0x0000,0x0000,0x0000,0x3ff0,*/
       0x12b2, 0x1cf3, 0xfd0d, 0xc075,
       0xd757, 0x7b89, 0xaa0d, 0xc0d0,
       0x4c9b, 0xb974, 0xeb84, 0xc10a,
       0x0043, 0x7195, 0x6286, 0xc131,
       0xf34c, 0x892f, 0x5255, 0xc143,
       0xe14a, 0x6a11, 0xce4b, 0xc13e
   };
   /* log( sqrt( 2*pi ) ) */
   static unsigned short LS2P[] = {
       0xbeb5, 0xc864, 0x67f1, 0x3fed
   };
   #define LS2PI *(double *)LS2P
   #define MAXLGM 2.556348e305
   #endif
   
   #ifdef MIEEE
   static unsigned short A[] = {
       0x3f4a, 0x9850, 0x2733, 0x6661,
       0xbf43, 0x7fbd, 0xb580, 0xe943,
       0x3f4a, 0x019f, 0x20dc, 0x5ebb,
       0xbf66, 0xc16c, 0x16b0, 0xa5a1,
       0x3fb5, 0x5555, 0x5555, 0x554b
   };
   static unsigned short B[] = {
       0xc095, 0x8901, 0x8ff3, 0x6761,
       0xc0e2, 0xf234, 0x355b, 0xb93e,
       0xc114, 0x3d73, 0xf890, 0x89e5,
       0xc131, 0xbc82, 0xf994, 0xdb51,
       0xc13a, 0x4589, 0x0219, 0xf20b,
       0xc12a, 0x0c67, 0x5418, 0x055e
   };
   static unsigned short C[] = {
       0xc075, 0xfd0d, 0x1cf3, 0x12b2,
       0xc0d0, 0xaa0d, 0x7b89, 0xd757,
       0xc10a, 0xeb84, 0xb974, 0x4c9b,
       0xc131, 0x6286, 0x7195, 0x0043,
       0xc143, 0x5255, 0x892f, 0xf34c,
       0xc13e, 0xce4b, 0x6a11, 0xe14a
   };
   /* log( sqrt( 2*pi ) ) */
   static unsigned short LS2P[] = {
       0x3fed, 0x67f1, 0xc864, 0xbeb5
   };
   #define LS2PI *(double *)LS2P
   #define MAXLGM 2.556348e305
   #endif
   
   /* static const double PI = 3.1415926535897932384626433832795028841971693993751;*/
   static const double LOGPI = 1.1447298858494001741434273513530587116472948129153;
   
   int             ISNAN(double x)
 {  {
     signgam = 1;      volatile double a = x;
       if (a != a)
           return 1;
       return 0;
   }
   int             ISFINITE(double x)
   {
       volatile double a = x;
       if (a < DBL_MAX)
           return 1;
       return 0;
   }
   
     if (z <= 0.0)  double          lngamma(double x)
         return lgamneg(z);  {
       double          p,
                       q,
                       u,
                       w,
                       z;
       int             i;
   
       sgngam = 1;
   #ifdef NANS
       if (ISNAN(x))
           return (x);
   #endif
   
   #ifdef INFINITIES
       if (!ISFINITE((x)))
           return (DBL_MAX * DBL_MAX);
   #endif
   
       if (x < -34.0) {
           q = -x;
           w = lngamma(q);            /* note this modifies sgngam! */
           p = floor(q);
           if (p == q) {
       lgsing:
   #ifdef INFINITIES
               mtherr("lngamma", SING);
               return (DBL_MAX * DBL_MAX);
   #else
               goto loverf;
   #endif
           }
           i = p;
           if ((i & 1) == 0)
               sgngam = -1;
           else
               sgngam = 1;
           z = q - p;
           if (z > 0.5) {
               p += 1.0;
               z = p - q;
           }
           z = q * sin(PI * z);
           if (z == 0.0)
               goto lgsing;
   /*      z = log(PI) - log( z ) - w;*/
           z = LOGPI - log(z) - w;
           return (z);
       }
       if (x < 13.0) {
           z = 1.0;
           p = 0.0;
           u = x;
           while (u >= 3.0) {
               p -= 1.0;
               u = x + p;
               z *= u;
           }
           while (u < 2.0) {
               if (u == 0.0)
                   goto lgsing;
               z /= u;
               p += 1.0;
               u = x + p;
           }
           if (z < 0.0) {
               sgngam = -1;
               z = -z;
           } else
               sgngam = 1;
           if (u == 2.0)
               return (log(z));
           p -= 2.0;
           x = x + p;
           p = x * polevl(x, B, 5) / p1evl(x, C, 6);
           return (log(z) + p);
       }
       if (x > MAXLGM) {
   #ifdef INFINITIES
           return (sgngam * (DBL_MAX * DBL_MAX));
   #else
   loverf:
           mtherr("lngamma", OVERFLOW);
           return (sgngam * MAXNUM);
   #endif
       }
       q = (x - 0.5) * log(x) - x + LS2PI;
       if (x > 1.0e8)
           return (q);
   
       p = 1.0 / (x * x);
       if (x >= 1000.0)
           q += ((7.9365079365079365079365e-4 * p
                  - 2.7777777777777777777778e-3) * p
                 + 0.0833333333333333333333) / x;
     else      else
         return lgampos(z);          q += polevl(p, A, 4) / x;
       return (q);
 }  }
   
 # define GAMMA(x) lngamma ((x))  #define GAMMA(x) lngamma ((x))
 #endif /* !GAMMA */  #endif /* !GAMMA */
   
 void f_erf()  void f_erf()
Line 363  void f_rand()
Line 686  void f_rand()
   
 #endif /* BADRAND */  #endif /* BADRAND */
   
 /** ibeta.c  /* ** ibeta.c
  *   *
  *   DESCRIB   Approximate the incomplete beta function Ix(a, b).   *   DESCRIB   Approximate the incomplete beta function Ix(a, b).
  *   *
Line 392  void f_rand()
Line 715  void f_rand()
  *   REFERENCE The continued fraction expansion as given by   *   REFERENCE The continued fraction expansion as given by
  *             Abramowitz and Stegun (1964) is used.   *             Abramowitz and Stegun (1964) is used.
  *   *
  * Permission granted to distribute freely for non-commercial purposes only  
  * Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl   * Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl
    *
    * Note: this function was translated from the Public Domain Fortran
    *       version available from http://lib.stat.cmu.edu/apstat/xxx
    *
  */   */
   
 static double ibeta(a, b, x)  static double ibeta(a, b, x)
Line 462  double a, b, x;
Line 788  double a, b, x;
     return -1.0;      return -1.0;
 }  }
   
 /** igamma.c  /* ** igamma.c
  *   *
  *   DESCRIB   Approximate the incomplete gamma function P(a, x).   *   DESCRIB   Approximate the incomplete gamma function P(a, x).
  *   *
Line 486  double a, b, x;
Line 812  double a, b, x;
  *   *
  *   REFERENCE ALGORITHM AS239  APPL. STATIST. (1988) VOL. 37, NO. 3   *   REFERENCE ALGORITHM AS239  APPL. STATIST. (1988) VOL. 37, NO. 3
  *   *
  * Permission granted to distribute freely for non-commercial purposes only  
  * Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl   * Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl
    *
    * Note: this function was translated from the Public Domain Fortran
    *       version available from http://lib.stat.cmu.edu/apstat/239
    *
  */   */
   
 /* Global variables, not visible outside this file */  /* Global variables, not visible outside this file */
Line 704  void f_inverse_erf()
Line 1033  void f_inverse_erf()
     }      }
 }  }
   
 static double inverse_normal_func(p)  /*                                                      ndtri.c
 double p;   *
 {   *      Inverse of Normal distribution function
     /*   *
        Source: This routine was derived (using f2c) from the   *
        FORTRAN subroutine MDNRIS found in   *
        ACM Algorithm 602 obtained from netlib.   * SYNOPSIS:
    *
    * double x, y, ndtri();
    *
    * x = ndtri( y );
    *
    *
    *
    * DESCRIPTION:
    *
    * Returns the argument, x, for which the area under the
    * Gaussian probability density function (integrated from
    * minus infinity to x) is equal to y.
    *
    *
    * For small arguments 0 < y < exp(-2), the program computes
    * z = sqrt( -2.0 * log(y) );  then the approximation is
    * x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z).
    * There are two rational functions P/Q, one for 0 < y < exp(-32)
    * and the other for y up to exp(-2).  For larger arguments,
    * w = y - 0.5, and  x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
    *
    *
    * ACCURACY:
    *
    *                      Relative error:
    * arithmetic   domain        # trials      peak         rms
    *    DEC      0.125, 1         5500       9.5e-17     2.1e-17
    *    DEC      6e-39, 0.135     3500       5.7e-17     1.3e-17
    *    IEEE     0.125, 1        20000       7.2e-16     1.3e-16
    *    IEEE     3e-308, 0.135   50000       4.6e-16     9.8e-17
    *
    *
    * ERROR MESSAGES:
    *
    *   message         condition    value returned
    * ndtri domain       x <= 0        -DBL_MAX
    * ndtri domain       x >= 1         DBL_MAX
    *
    */
   
        MDNRIS code contains the 1978 Copyright  /*
        by IMSL, INC. .  Since MDNRIS has been  Cephes Math Library Release 2.8:  June, 2000
        submitted to netlib it may be used with  Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
        the restriction that it may only be  */
        used for noncommercial purposes and that  
        IMSL be acknowledged as the copyright-holder  
        of the code.  
      */  
   
     /* Initialized data */  #ifdef UNK
     static double eps = 1e-10;  /* sqrt(2pi) */
     static double g0 = 1.851159e-4;  static double   s2pi = 2.50662827463100050242E0;
     static double g1 = -.002028152;  #endif
     static double g2 = -.1498384;  
     static double g3 = .01078639;  
     static double h0 = .09952975;  
     static double h1 = .5211733;  
     static double h2 = -.06888301;  
     static double sqrt2 = 1.414213562373095;  
   
     /* Local variables */  #ifdef DEC
     static double a, w, x;  static unsigned short s2p[] = {0040440, 0066230, 0177661, 0034055};
     static double sd, wi, sn, y;  #define s2pi *(double *)s2p
   #endif
   
     /* Note: 0.0 < p < 1.0 */  #ifdef IBMPC
   static unsigned short s2p[] = {0x2706, 0x1ff6, 0x0d93, 0x4004};
   #define s2pi *(double *)s2p
   #endif
   
     /* p too small, compute y directly */  #ifdef MIEEE
     if (p <= eps) {  static unsigned short s2p[] = {
         a = p + p;      0x4004, 0x0d93, 0x1ff6, 0x2706
         w = sqrt(-(double) log(a + (a - a * a)));  };
   #define s2pi *(double *)s2p
   #endif
   
         /* use a rational function in 1.0 / w */  /* approximation for 0 <= |y - 0.5| <= 3/8 */
         wi = 1.0 / w;  #ifdef UNK
         sn = ((g3 * wi + g2) * wi + g1) * wi;  static double   P0[5] = {
         sd = ((wi + h2) * wi + h1) * wi + h0;      -5.99633501014107895267E1,
         y = w + w * (g0 + sn / sd);      9.80010754185999661536E1,
         y = -y * sqrt2;      -5.66762857469070293439E1,
     } else {      1.39312609387279679503E1,
         x = 1.0 - (p + p);      -1.23916583867381258016E0,
         y = inverse_error_func(x);  };
         y = -sqrt2 * y;  static double   Q0[8] = {
   /* 1.00000000000000000000E0,*/
       1.95448858338141759834E0,
       4.67627912898881538453E0,
       8.63602421390890590575E1,
       -2.25462687854119370527E2,
       2.00260212380060660359E2,
       -8.20372256168333339912E1,
       1.59056225126211695515E1,
       -1.18331621121330003142E0,
   };
   #endif
   #ifdef DEC
   static unsigned short P0[20] = {
       0141557, 0155170, 0071360, 0120550,
       0041704, 0000214, 0172417, 0067307,
       0141542, 0132204, 0040066, 0156723,
       0041136, 0163161, 0157276, 0007747,
       0140236, 0116374, 0073666, 0051764,
   };
   static unsigned short Q0[32] = {
   /*0040200,0000000,0000000,0000000,*/
       0040372, 0026256, 0110403, 0123707,
       0040625, 0122024, 0020277, 0026661,
       0041654, 0134161, 0124134, 0007244,
       0142141, 0073162, 0133021, 0131371,
       0042110, 0041235, 0043516, 0057767,
       0141644, 0011417, 0036155, 0137305,
       0041176, 0076556, 0004043, 0125430,
       0140227, 0073347, 0152776, 0067251,
   };
   #endif
   #ifdef IBMPC
   static unsigned short P0[20] = {
       0x142d, 0x0e5e, 0xfb4f, 0xc04d,
       0xedd9, 0x9ea1, 0x8011, 0x4058,
       0xdbba, 0x8806, 0x5690, 0xc04c,
       0xc1fd, 0x3bd7, 0xdcce, 0x402b,
       0xca7e, 0x8ef6, 0xd39f, 0xbff3,
   };
   static unsigned short Q0[36] = {
   /*0x0000,0x0000,0x0000,0x3ff0,*/
       0x74f9, 0xd220, 0x4595, 0x3fff,
       0xe5b6, 0x8417, 0xb482, 0x4012,
       0x81d4, 0x350b, 0x970e, 0x4055,
       0x365f, 0x56c2, 0x2ece, 0xc06c,
       0xcbff, 0xa8e9, 0x0853, 0x4069,
       0xb7d9, 0xe78d, 0x8261, 0xc054,
       0x7563, 0xc104, 0xcfad, 0x402f,
       0xcdd5, 0xfabf, 0xeedc, 0xbff2,
   };
   #endif
   #ifdef MIEEE
   static unsigned short P0[20] = {
       0xc04d, 0xfb4f, 0x0e5e, 0x142d,
       0x4058, 0x8011, 0x9ea1, 0xedd9,
       0xc04c, 0x5690, 0x8806, 0xdbba,
       0x402b, 0xdcce, 0x3bd7, 0xc1fd,
       0xbff3, 0xd39f, 0x8ef6, 0xca7e,
   };
   static unsigned short Q0[32] = {
   /*0x3ff0,0x0000,0x0000,0x0000,*/
       0x3fff, 0x4595, 0xd220, 0x74f9,
       0x4012, 0xb482, 0x8417, 0xe5b6,
       0x4055, 0x970e, 0x350b, 0x81d4,
       0xc06c, 0x2ece, 0x56c2, 0x365f,
       0x4069, 0x0853, 0xa8e9, 0xcbff,
       0xc054, 0x8261, 0xe78d, 0xb7d9,
       0x402f, 0xcfad, 0xc104, 0x7563,
       0xbff2, 0xeedc, 0xfabf, 0xcdd5,
   };
   #endif
   
   /* Approximation for interval z = sqrt(-2 log y ) between 2 and 8
    * i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
    */
   #ifdef UNK
   static double   P1[9] = {
       4.05544892305962419923E0,
       3.15251094599893866154E1,
       5.71628192246421288162E1,
       4.40805073893200834700E1,
       1.46849561928858024014E1,
       2.18663306850790267539E0,
       -1.40256079171354495875E-1,
       -3.50424626827848203418E-2,
       -8.57456785154685413611E-4,
   };
   static double   Q1[8] = {
   /*  1.00000000000000000000E0,*/
       1.57799883256466749731E1,
       4.53907635128879210584E1,
       4.13172038254672030440E1,
       1.50425385692907503408E1,
       2.50464946208309415979E0,
       -1.42182922854787788574E-1,
       -3.80806407691578277194E-2,
       -9.33259480895457427372E-4,
   };
   #endif
   #ifdef DEC
   static unsigned short P1[36] = {
       0040601, 0143074, 0150744, 0073326,
       0041374, 0031554, 0113253, 0146016,
       0041544, 0123272, 0012463, 0176771,
       0041460, 0051160, 0103560, 0156511,
       0041152, 0172624, 0117772, 0030755,
       0040413, 0170713, 0151545, 0176413,
       0137417, 0117512, 0022154, 0131671,
       0137017, 0104257, 0071432, 0007072,
       0135540, 0143363, 0063137, 0036166,
   };
   static unsigned short Q1[32] = {
   /*0040200,0000000,0000000,0000000,*/
       0041174, 0075325, 0004736, 0120326,
       0041465, 0110044, 0047561, 0045567,
       0041445, 0042321, 0012142, 0030340,
       0041160, 0127074, 0166076, 0141051,
       0040440, 0046055, 0040745, 0150400,
       0137421, 0114146, 0067330, 0010621,
       0137033, 0175162, 0025555, 0114351,
       0135564, 0122773, 0145750, 0030357,
   };
   #endif
   #ifdef IBMPC
   static unsigned short P1[36] = {
       0x8edb, 0x9a3c, 0x38c7, 0x4010,
       0x7982, 0x92d5, 0x866d, 0x403f,
       0x7fbf, 0x42a6, 0x94d7, 0x404c,
       0x1ba9, 0x10ee, 0x0a4e, 0x4046,
       0x463e, 0x93ff, 0x5eb2, 0x402d,
       0xbfa1, 0x7a6c, 0x7e39, 0x4001,
       0x9677, 0x448d, 0xf3e9, 0xbfc1,
       0x41c7, 0xee63, 0xf115, 0xbfa1,
       0xe78f, 0x6ccb, 0x18de, 0xbf4c,
   };
   static unsigned short Q1[32] = {
   /*0x0000,0x0000,0x0000,0x3ff0,*/
       0xd41b, 0xa13b, 0x8f5a, 0x402f,
       0x296f, 0x89ee, 0xb204, 0x4046,
       0x461c, 0x228c, 0xa89a, 0x4044,
       0xd845, 0x9d87, 0x15c7, 0x402e,
       0xba20, 0xa83c, 0x0985, 0x4004,
       0x0232, 0xcddb, 0x330c, 0xbfc2,
       0xb31d, 0x456d, 0x7f4e, 0xbfa3,
       0x061e, 0x797d, 0x94bf, 0xbf4e,
   };
   #endif
   #ifdef MIEEE
   static unsigned short P1[36] = {
       0x4010, 0x38c7, 0x9a3c, 0x8edb,
       0x403f, 0x866d, 0x92d5, 0x7982,
       0x404c, 0x94d7, 0x42a6, 0x7fbf,
       0x4046, 0x0a4e, 0x10ee, 0x1ba9,
       0x402d, 0x5eb2, 0x93ff, 0x463e,
       0x4001, 0x7e39, 0x7a6c, 0xbfa1,
       0xbfc1, 0xf3e9, 0x448d, 0x9677,
       0xbfa1, 0xf115, 0xee63, 0x41c7,
       0xbf4c, 0x18de, 0x6ccb, 0xe78f,
   };
   static unsigned short Q1[32] = {
   /*0x3ff0,0x0000,0x0000,0x0000,*/
       0x402f, 0x8f5a, 0xa13b, 0xd41b,
       0x4046, 0xb204, 0x89ee, 0x296f,
       0x4044, 0xa89a, 0x228c, 0x461c,
       0x402e, 0x15c7, 0x9d87, 0xd845,
       0x4004, 0x0985, 0xa83c, 0xba20,
       0xbfc2, 0x330c, 0xcddb, 0x0232,
       0xbfa3, 0x7f4e, 0x456d, 0xb31d,
       0xbf4e, 0x94bf, 0x797d, 0x061e,
   };
   #endif
   
   /* Approximation for interval z = sqrt(-2 log y ) between 8 and 64
    * i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
    */
   
   #ifdef UNK
   static double   P2[9] = {
       3.23774891776946035970E0,
       6.91522889068984211695E0,
       3.93881025292474443415E0,
       1.33303460815807542389E0,
       2.01485389549179081538E-1,
       1.23716634817820021358E-2,
       3.01581553508235416007E-4,
       2.65806974686737550832E-6,
       6.23974539184983293730E-9,
   };
   static double   Q2[8] = {
   /*  1.00000000000000000000E0,*/
       6.02427039364742014255E0,
       3.67983563856160859403E0,
       1.37702099489081330271E0,
       2.16236993594496635890E-1,
       1.34204006088543189037E-2,
       3.28014464682127739104E-4,
       2.89247864745380683936E-6,
       6.79019408009981274425E-9,
   };
   #endif
   #ifdef DEC
   static unsigned short P2[36] = {
       0040517, 0033507, 0036236, 0125641,
       0040735, 0044616, 0014473, 0140133,
       0040574, 0012567, 0114535, 0102541,
       0040252, 0120340, 0143474, 0150135,
       0037516, 0051057, 0115361, 0031211,
       0036512, 0131204, 0101511, 0125144,
       0035236, 0016627, 0043160, 0140216,
       0033462, 0060512, 0060141, 0010641,
       0031326, 0062541, 0101304, 0077706,
   };
   static unsigned short Q2[32] = {
   /*0040200,0000000,0000000,0000000,*/
       0040700, 0143322, 0132137, 0040501,
       0040553, 0101155, 0053221, 0140257,
       0040260, 0041071, 0052573, 0010004,
       0037535, 0066472, 0177261, 0162330,
       0036533, 0160475, 0066666, 0036132,
       0035253, 0174533, 0027771, 0044027,
       0033502, 0016147, 0117666, 0063671,
       0031351, 0047455, 0141663, 0054751,
   };
   #endif
   #ifdef IBMPC
   static unsigned short P2[36] = {
       0xd574, 0xe793, 0xe6e8, 0x4009,
       0x780b, 0xc327, 0xa931, 0x401b,
       0xb0ac, 0xf32b, 0x82ae, 0x400f,
       0x9a0c, 0x18e7, 0x541c, 0x3ff5,
       0x2651, 0xf35e, 0xca45, 0x3fc9,
       0x354d, 0x9069, 0x5650, 0x3f89,
       0x1812, 0xe8ce, 0xc3b2, 0x3f33,
       0x2234, 0x4c0c, 0x4c29, 0x3ec6,
       0x8ff9, 0x3058, 0xccac, 0x3e3a,
   };
   static unsigned short Q2[32] = {
   /*0x0000,0x0000,0x0000,0x3ff0,*/
       0xe828, 0x568b, 0x18da, 0x4018,
       0x3816, 0xaad2, 0x704d, 0x400d,
       0x6200, 0x2aaf, 0x0847, 0x3ff6,
       0x3c9b, 0x5fd6, 0xada7, 0x3fcb,
       0xc78b, 0xadb6, 0x7c27, 0x3f8b,
       0x2903, 0x65ff, 0x7f2b, 0x3f35,
       0xccf7, 0xf3f6, 0x438c, 0x3ec8,
       0x6b3d, 0xb876, 0x29e5, 0x3e3d,
   };
   #endif
   #ifdef MIEEE
   static unsigned short P2[36] = {
       0x4009, 0xe6e8, 0xe793, 0xd574,
       0x401b, 0xa931, 0xc327, 0x780b,
       0x400f, 0x82ae, 0xf32b, 0xb0ac,
       0x3ff5, 0x541c, 0x18e7, 0x9a0c,
       0x3fc9, 0xca45, 0xf35e, 0x2651,
       0x3f89, 0x5650, 0x9069, 0x354d,
       0x3f33, 0xc3b2, 0xe8ce, 0x1812,
       0x3ec6, 0x4c29, 0x4c0c, 0x2234,
       0x3e3a, 0xccac, 0x3058, 0x8ff9,
   };
   static unsigned short Q2[32] = {
   /*0x3ff0,0x0000,0x0000,0x0000,*/
       0x4018, 0x18da, 0x568b, 0xe828,
       0x400d, 0x704d, 0xaad2, 0x3816,
       0x3ff6, 0x0847, 0x2aaf, 0x6200,
       0x3fcb, 0xada7, 0x5fd6, 0x3c9b,
       0x3f8b, 0x7c27, 0xadb6, 0xc78b,
       0x3f35, 0x7f2b, 0x65ff, 0x2903,
       0x3ec8, 0x438c, 0xf3f6, 0xccf7,
       0x3e3d, 0x29e5, 0xb876, 0x6b3d,
   };
   #endif
   
   static double inverse_normal_func(double y0)
   {
       double          x,
                       y,
                       z,
                       y2,
                       x0,
                       x1;
       int             code;
   
       if (y0 <= 0.0) {
           mtherr("inverse_normal_func", DOMAIN);
           return (-DBL_MAX);
     }      }
     return (y);      if (y0 >= 1.0) {
           mtherr("inverse_normal_func", DOMAIN);
           return (DBL_MAX);
       }
       code = 1;
       y = y0;
       if (y > (1.0 - 0.13533528323661269189)) {   /* 0.135... = exp(-2) */
           y = 1.0 - y;
           code = 0;
       }
       if (y > 0.13533528323661269189) {
           y = y - 0.5;
           y2 = y * y;
           x = y + y * (y2 * polevl(y2, P0, 4) / p1evl(y2, Q0, 8));
           x = x * s2pi;
           return (x);
       }
       x = sqrt(-2.0 * log(y));
       x0 = x - log(x) / x;
   
       z = 1.0 / x;
       if (x < 8.0)                /* y > exp(-32) = 1.2664165549e-14 */
           x1 = z * polevl(z, P1, 8) / p1evl(z, Q1, 8);
       else
           x1 = z * polevl(z, P2, 8) / p1evl(z, Q2, 8);
       x = x0 - x1;
       if (code != 0)
           x = -x;
       return (x);
 }  }
   /*                                                      ndtr.c
    *
    *      Normal distribution function
    *
    *
    *
    * SYNOPSIS:
    *
    * double x, y, ndtr();
    *
    * y = ndtr( x );
    *
    *
    *
    * DESCRIPTION:
    *
    * Returns the area under the Gaussian probability density
    * function, integrated from minus infinity to x:
    *
    *                            x
    *                             -
    *                   1        | |          2
    *    ndtr(x)  = ---------    |    exp( - t /2 ) dt
    *               sqrt(2pi)  | |
    *                           -
    *                          -inf.
    *
    *             =  ( 1 + erf(z) ) / 2
    *             =  erfc(z) / 2
    *
    * where z = x/sqrt(2). Computation is via the functions
    * erf and erfc.
    *
    *
    * ACCURACY:
    *
    *                      Relative error:
    * arithmetic   domain     # trials      peak         rms
    *    DEC      -13,0         8000       2.1e-15     4.8e-16
    *    IEEE     -13,0        30000       3.4e-14     6.7e-15
    *
    *
    * ERROR MESSAGES:
    *
    *   message         condition         value returned
    * erfc underflow    x > 37.519379347       0.0
    *
    */
   
   /*                                                     erf.c
    *
    *      Error function
    *
    *
    *
    * SYNOPSIS:
    *
    * double x, y, erf();
    *
    * y = erf( x );
    *
    *
    *
    * DESCRIPTION:
    *
    * The integral is
    *
    *                           x
    *                            -
    *                 2         | |          2
    *   erf(x)  =  --------     |    exp( - t  ) dt.
    *              sqrt(pi)   | |
    *                          -
    *                           0
    *
    * The magnitude of x is limited to 9.231948545 for DEC
    * arithmetic; 1 or -1 is returned outside this range.
    *
    * For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
    * erf(x) = 1 - erfc(x).
    *
    *
    *
    * ACCURACY:
    *
    *                      Relative error:
    * arithmetic   domain     # trials      peak         rms
    *    DEC       0,1         14000       4.7e-17     1.5e-17
    *    IEEE      0,1         30000       3.7e-16     1.0e-16
    *
    */
   
 static double inverse_error_func(p)  /*                                                     erfc.c
 double p;   *
    *      Complementary error function
    *
    *
    *
    * SYNOPSIS:
    *
    * double x, y, erfc();
    *
    * y = erfc( x );
    *
    *
    *
    * DESCRIPTION:
    *
    *
    *  1 - erf(x) =
    *
    *                           inf.
    *                             -
    *                  2         | |          2
    *   erfc(x)  =  --------     |    exp( - t  ) dt
    *               sqrt(pi)   | |
    *                           -
    *                            x
    *
    *
    * For small x, erfc(x) = 1 - erf(x); otherwise rational
    * approximations are computed.
    *
    *
    *
    * ACCURACY:
    *
    *                      Relative error:
    * arithmetic   domain     # trials      peak         rms
    *    DEC       0, 9.2319   12000       5.1e-16     1.2e-16
    *    IEEE      0,26.6417   30000       5.7e-14     1.5e-14
    *
    *
    * ERROR MESSAGES:
    *
    *   message         condition              value returned
    * erfc underflow    x > 9.231948545 (DEC)       0.0
    *
    *
    */
   
   /*
   Cephes Math Library Release 2.8:  June, 2000
   Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
   */
   
   static const double SQRTH = 0.70710678118654752440084436210484903928483593768845;
   
   #ifdef UNK
   static double   P[] = {
       2.46196981473530512524E-10,
       5.64189564831068821977E-1,
       7.46321056442269912687E0,
       4.86371970985681366614E1,
       1.96520832956077098242E2,
       5.26445194995477358631E2,
       9.34528527171957607540E2,
       1.02755188689515710272E3,
       5.57535335369399327526E2
   };
   static double   Q[] = {
   /* 1.00000000000000000000E0,*/
       1.32281951154744992508E1,
       8.67072140885989742329E1,
       3.54937778887819891062E2,
       9.75708501743205489753E2,
       1.82390916687909736289E3,
       2.24633760818710981792E3,
       1.65666309194161350182E3,
       5.57535340817727675546E2
   };
   static double   R[] = {
       5.64189583547755073984E-1,
       1.27536670759978104416E0,
       5.01905042251180477414E0,
       6.16021097993053585195E0,
       7.40974269950448939160E0,
       2.97886665372100240670E0
   };
   static double   S[] = {
   /* 1.00000000000000000000E0,*/
       2.26052863220117276590E0,
       9.39603524938001434673E0,
       1.20489539808096656605E1,
       1.70814450747565897222E1,
       9.60896809063285878198E0,
       3.36907645100081516050E0
   };
   static double   T[] = {
       9.60497373987051638749E0,
       9.00260197203842689217E1,
       2.23200534594684319226E3,
       7.00332514112805075473E3,
       5.55923013010394962768E4
   };
   static double   U[] = {
   /* 1.00000000000000000000E0,*/
       3.35617141647503099647E1,
       5.21357949780152679795E2,
       4.59432382970980127987E3,
       2.26290000613890934246E4,
       4.92673942608635921086E4
   };
   
   #define UTHRESH 37.519379347
   #endif
   
   #ifdef DEC
   static unsigned short P[] = {
       0030207, 0054445, 0011173, 0021706,
       0040020, 0067272, 0030661, 0122075,
       0040756, 0151236, 0173053, 0067042,
       0041502, 0106175, 0062555, 0151457,
       0042104, 0102525, 0047401, 0003667,
       0042403, 0116176, 0011446, 0075303,
       0042551, 0120723, 0061641, 0123275,
       0042600, 0070651, 0007264, 0134516,
       0042413, 0061102, 0167507, 0176625
   };
   static unsigned short Q[] = {
   /*0040200,0000000,0000000,0000000,*/
       0041123, 0123257, 0165741, 0017142,
       0041655, 0065027, 0173413, 0115450,
       0042261, 0074011, 0021573, 0004150,
       0042563, 0166530, 0013662, 0007200,
       0042743, 0176427, 0162443, 0105214,
       0043014, 0062546, 0153727, 0123772,
       0042717, 0012470, 0006227, 0067424,
       0042413, 0061103, 0003042, 0013254
   };
   static unsigned short R[] = {
       0040020, 0067272, 0101024, 0155421,
       0040243, 0037467, 0056706, 0026462,
       0040640, 0116017, 0120665, 0034315,
       0040705, 0020162, 0143350, 0060137,
       0040755, 0016234, 0134304, 0130157,
       0040476, 0122700, 0051070, 0015473
   };
   static unsigned short S[] = {
   /*0040200,0000000,0000000,0000000,*/
       0040420, 0126200, 0044276, 0070413,
       0041026, 0053051, 0007302, 0063746,
       0041100, 0144203, 0174051, 0061151,
       0041210, 0123314, 0126343, 0177646,
       0041031, 0137125, 0051431, 0033011,
       0040527, 0117362, 0152661, 0066201
   };
   static unsigned short T[] = {
       0041031, 0126770, 0170672, 0166101,
       0041664, 0006522, 0072360, 0031770,
       0043013, 0100025, 0162641, 0126671,
       0043332, 0155231, 0161627, 0076200,
       0044131, 0024115, 0021020, 0117343
   };
   static unsigned short U[] = {
   /*0040200,0000000,0000000,0000000,*/
       0041406, 0037461, 0177575, 0032714,
       0042402, 0053350, 0123061, 0153557,
       0043217, 0111227, 0032007, 0164217,
       0043660, 0145000, 0004013, 0160114,
       0044100, 0071544, 0167107, 0125471
   };
   #define UTHRESH 14.0
   #endif
   
   #ifdef IBMPC
   static unsigned short P[] = {
       0x6479, 0xa24f, 0xeb24, 0x3df0,
       0x3488, 0x4636, 0x0dd7, 0x3fe2,
       0x6dc4, 0xdec5, 0xda53, 0x401d,
       0xba66, 0xacad, 0x518f, 0x4048,
       0x20f7, 0xa9e0, 0x90aa, 0x4068,
       0xcf58, 0xc264, 0x738f, 0x4080,
       0x34d8, 0x6c74, 0x343a, 0x408d,
       0x972a, 0x21d6, 0x0e35, 0x4090,
       0xffb3, 0x5de8, 0x6c48, 0x4081
   };
   static unsigned short Q[] = {
   /*0x0000,0x0000,0x0000,0x3ff0,*/
       0x23cc, 0xfd7c, 0x74d5, 0x402a,
       0x7365, 0xfee1, 0xad42, 0x4055,
       0x610d, 0x246f, 0x2f01, 0x4076,
       0x41d0, 0x02f6, 0x7dab, 0x408e,
       0x7151, 0xfca4, 0x7fa2, 0x409c,
       0xf4ff, 0xdafa, 0x8cac, 0x40a1,
       0xede2, 0x0192, 0xe2a7, 0x4099,
       0x42d6, 0x60c4, 0x6c48, 0x4081
   };
   static unsigned short R[] = {
       0x9b62, 0x5042, 0x0dd7, 0x3fe2,
       0xc5a6, 0xebb8, 0x67e6, 0x3ff4,
       0xa71a, 0xf436, 0x1381, 0x4014,
       0x0c0c, 0x58dd, 0xa40e, 0x4018,
       0x960e, 0x9718, 0xa393, 0x401d,
       0x0367, 0x0a47, 0xd4b8, 0x4007
   };
   static unsigned short S[] = {
   /*0x0000,0x0000,0x0000,0x3ff0,*/
       0xce21, 0x0917, 0x1590, 0x4002,
       0x4cfd, 0x21d8, 0xcac5, 0x4022,
       0x2c4d, 0x7f05, 0x1910, 0x4028,
       0x7ff5, 0x959c, 0x14d9, 0x4031,
       0x26c1, 0xaa63, 0x37ca, 0x4023,
       0x2d90, 0x5ab6, 0xf3de, 0x400a
   };
   static unsigned short T[] = {
       0x5d88, 0x1e37, 0x35bf, 0x4023,
       0x067f, 0x4e9e, 0x81aa, 0x4056,
       0x35b7, 0xbcb4, 0x7002, 0x40a1,
       0xef90, 0x3c72, 0x5b53, 0x40bb,
       0x13dc, 0xa442, 0x2509, 0x40eb
   };
   static unsigned short U[] = {
   /*0x0000,0x0000,0x0000,0x3ff0,*/
       0xa6ba, 0x3fef, 0xc7e6, 0x4040,
       0x3aee, 0x14c6, 0x4add, 0x4080,
       0xfd12, 0xe680, 0xf252, 0x40b1,
       0x7c0a, 0x0101, 0x1940, 0x40d6,
       0xf567, 0x9dc8, 0x0e6c, 0x40e8
   };
   #define UTHRESH 37.519379347
   #endif
   
   #ifdef MIEEE
   static unsigned short P[] = {
       0x3df0, 0xeb24, 0xa24f, 0x6479,
       0x3fe2, 0x0dd7, 0x4636, 0x3488,
       0x401d, 0xda53, 0xdec5, 0x6dc4,
       0x4048, 0x518f, 0xacad, 0xba66,
       0x4068, 0x90aa, 0xa9e0, 0x20f7,
       0x4080, 0x738f, 0xc264, 0xcf58,
       0x408d, 0x343a, 0x6c74, 0x34d8,
       0x4090, 0x0e35, 0x21d6, 0x972a,
       0x4081, 0x6c48, 0x5de8, 0xffb3
   };
   static unsigned short Q[] = {
       0x402a, 0x74d5, 0xfd7c, 0x23cc,
       0x4055, 0xad42, 0xfee1, 0x7365,
       0x4076, 0x2f01, 0x246f, 0x610d,
       0x408e, 0x7dab, 0x02f6, 0x41d0,
       0x409c, 0x7fa2, 0xfca4, 0x7151,
       0x40a1, 0x8cac, 0xdafa, 0xf4ff,
       0x4099, 0xe2a7, 0x0192, 0xede2,
       0x4081, 0x6c48, 0x60c4, 0x42d6
   };
   static unsigned short R[] = {
       0x3fe2, 0x0dd7, 0x5042, 0x9b62,
       0x3ff4, 0x67e6, 0xebb8, 0xc5a6,
       0x4014, 0x1381, 0xf436, 0xa71a,
       0x4018, 0xa40e, 0x58dd, 0x0c0c,
       0x401d, 0xa393, 0x9718, 0x960e,
       0x4007, 0xd4b8, 0x0a47, 0x0367
   };
   static unsigned short S[] = {
       0x4002, 0x1590, 0x0917, 0xce21,
       0x4022, 0xcac5, 0x21d8, 0x4cfd,
       0x4028, 0x1910, 0x7f05, 0x2c4d,
       0x4031, 0x14d9, 0x959c, 0x7ff5,
       0x4023, 0x37ca, 0xaa63, 0x26c1,
       0x400a, 0xf3de, 0x5ab6, 0x2d90
   };
   static unsigned short T[] = {
       0x4023, 0x35bf, 0x1e37, 0x5d88,
       0x4056, 0x81aa, 0x4e9e, 0x067f,
       0x40a1, 0x7002, 0xbcb4, 0x35b7,
       0x40bb, 0x5b53, 0x3c72, 0xef90,
       0x40eb, 0x2509, 0xa442, 0x13dc
   };
   static unsigned short U[] = {
       0x4040, 0xc7e6, 0x3fef, 0xa6ba,
       0x4080, 0x4add, 0x14c6, 0x3aee,
       0x40b1, 0xf252, 0xe680, 0xfd12,
       0x40d6, 0x1940, 0x0101, 0x7c0a,
       0x40e8, 0x0e6c, 0x9dc8, 0xf567
   };
   #define UTHRESH 37.519379347
   #endif
   /*
   double          ndtr(double a)
 {  {
     /*      double          x,
        Source: This routine was derived (using f2c) from the                      y,
        FORTRAN subroutine MERFI found in                      z;
        ACM Algorithm 602 obtained from netlib.  
   
        MDNRIS code contains the 1978 Copyright      x = a * SQRTH;
        by IMSL, INC. .  Since MERFI has been      z = fabs(x);
        submitted to netlib, it may be used with  
        the restriction that it may only be  
        used for noncommercial purposes and that  
        IMSL be acknowledged as the copyright-holder  
        of the code.  
      */  
   
       if (z < SQRTH)
           y = 0.5 + 0.5 * erf(x);
   
       else {
           y = 0.5 * erfc(z);
   
     /* Initialized data */          if (x > 0)
     static double a1 = -.5751703;              y = 1.0 - y;
     static double a2 = -1.896513;      }
     static double a3 = -.05496261;  
     static double b0 = -.113773;  
     static double b1 = -3.293474;  
     static double b2 = -2.374996;  
     static double b3 = -1.187515;  
     static double c0 = -.1146666;  
     static double c1 = -.1314774;  
     static double c2 = -.2368201;  
     static double c3 = .05073975;  
     static double d0 = -44.27977;  
     static double d1 = 21.98546;  
     static double d2 = -7.586103;  
     static double e0 = -.05668422;  
     static double e1 = .3937021;  
     static double e2 = -.3166501;  
     static double e3 = .06208963;  
     static double f0 = -6.266786;  
     static double f1 = 4.666263;  
     static double f2 = -2.962883;  
     static double g0 = 1.851159e-4;  
     static double g1 = -.002028152;  
     static double g2 = -.1498384;  
     static double g3 = .01078639;  
     static double h0 = .09952975;  
     static double h1 = .5211733;  
     static double h2 = -.06888301;  
   
     /* Local variables */      return (y);
     static double a, b, f, w, x, y, z, sigma, z2, sd, wi, sn;  }
   */
   double erf(double);
   
     x = p;  double          erfc(double a)
   {
       double          p,
                       q,
                       x,
                       y,
                       z;
   
     /* determine sign of x */      if (a < 0.0)
     if (x > 0)          x = -a;
         sigma = 1.0;  
     else      else
         sigma = -1.0;          x = a;
   
     /* Note: -1.0 < x < 1.0 */      if (x < 1.0)
           return (1.0 - erf(a));
   
     z = fabs(x);      z = -a * a;
   
     /* z between 0.0 and 0.85, approx. f by a      if (z < DBL_MIN_10_EXP) {
        rational function in z  */  under:
           mtherr("erfc", UNDERFLOW);
           if (a < 0)
               return (2.0);
           else
               return (0.0);
       }
       z = exp(z);
   
     if (z <= 0.85) {      if (x < 8.0) {
         z2 = z * z;          p = polevl(x, P, 8);
         f = z + z * (b0 + a1 * z2 / (b1 + z2 + a2          q = p1evl(x, Q, 8);
                                      / (b2 + z2 + a3 / (b3 + z2))));  
   
         /* z greater than 0.85 */  
     } else {      } else {
         a = 1.0 - z;          p = polevl(x, R, 5);
         b = z;          q = p1evl(x, S, 6);
       }
       y = (z * p) / q;
   
         /* reduced argument is in (0.85,1.0),      if (a < 0)
            obtain the transformed variable */          y = 2.0 - y;
   
         w = sqrt(-(double) log(a + a * b));      if (y == 0.0)
           goto under;
   
         /* w greater than 4.0, approx. f by a      return (y);
            rational function in 1.0 / w */  }
   
         if (w >= 4.0) {  double          erf(double x)
             wi = 1.0 / w;  {
             sn = ((g3 * wi + g2) * wi + g1) * wi;      double          y,
             sd = ((wi + h2) * wi + h1) * wi + h0;                      z;
             f = w + w * (g0 + sn / sd);  
   
             /* w between 2.5 and 4.0, approx.      if (fabs(x) > 1.0)
                f by a rational function in w */          return (1.0 - erfc(x));
       z = x * x;
       y = x * polevl(z, T, 4) / p1evl(z, U, 5);
       return (y);
   
         } else if (w < 4.0 && w > 2.5) {  }
             sn = ((e3 * w + e2) * w + e1) * w;  
             sd = ((w + f2) * w + f1) * w + f0;  
             f = w + w * (e0 + sn / sd);  
   
             /* w between 1.13222 and 2.5, approx. f by  static double inverse_error_func(double y)
                a rational function in w */  {
         } else if (w <= 2.5 && w > 1.13222) {      double          x = 0.0;    /* The output */
             sn = ((c3 * w + c2) * w + c1) * w;      double          z = 0.0;    /* Intermadiate variable */
             sd = ((w + d2) * w + d1) * w + d0;      double          y0 = 0.7;   /* Central range variable */
             f = w + w * (c0 + sn / sd);  
         }      /* Coefficients in rational approximations. */
       double          a[4] = {0.886226899, -1.645349621, 0.914624893, -0.140543331};
       double          b[4] = {-2.118377725, 1.442710462, -0.329097515, 0.012229801};
       double          c[4] = {-1.970840454, -1.624906493, 3.429567803, 1.641345311};
       double          d[2] = {3.543889200, 1.637067800};
   
       if ((y < -1.0) || (1.0 < y)) {
           printf("inverse_error_func: The value out of the range of the function");
           x = log(-1.0);
       } else if ((y == -1.0) || (1.0 == y)) {
           x = -y * log(0.0);
       } else if ((-1.0 < y) && (y < -y0)) {
           z = sqrt(-log((1.0 + y) / 2.0));
           x = -(((c[3] * z + c[2]) * z + c[1]) * z + c[0]) / ((d[1] * z + d[0]) * z + 1.0);
       } else {
           if ((-y0 < y) && (y < y0)) {
               z = y * y;
               x = y * (((a[3] * z + a[2]) * z + a[1]) * z + a[0]) /
                   ((((b[3] * z + b[3]) * z + b[1]) * z + b[0]) * z + 1.0);
           } else if ((y0 < y) && (y < 1.0)) {
               z = sqrt(-log((1.0 - y) / 2.0));
               x = (((c[3] * z + c[2]) * z + c[1]) * z + c[0]) / ((d[1] * z + d[0]) * z + 1.0);
           }
           /* Two steps of Newton-Raphson correction to full accuracy. */
           x = x - (erf(x) - y) / (2.0 / sqrt(PI) * exp(-x * x));
           x = x - (erf(x) - y) / (2.0 / sqrt(PI) * exp(-x * x));
     }      }
     y = sigma * f;      return (x);
     return (y);  
 }  }
Legend:
Removed from v.1.1.1.2  
changed lines
  Added in v.1.1.1.3
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>