git.maemo.org Git - gnuplot/blob - src/specfun.c

   1 #ifndef lint
   2 static char *RCSid() { return RCSid("$Id: specfun.c,v 1.34.2.1 2008/05/21 18:14:23 sfeam Exp $"); }
   3 #endif
   4
   5 /* GNUPLOT - specfun.c */
   6
   7 /*[
   8  * Copyright 1986 - 1993, 1998, 2004   Thomas Williams, Colin Kelley
   9  *
  10  * Permission to use, copy, and distribute this software and its
  11  * documentation for any purpose with or without fee is hereby granted,
  12  * provided that the above copyright notice appear in all copies and
  13  * that both that copyright notice and this permission notice appear
  14  * in supporting documentation.
  15  *
  16  * Permission to modify the software is granted, but not the right to
  17  * distribute the complete modified source code.  Modifications are to
  18  * be distributed as patches to the released version.  Permission to
  19  * distribute binaries produced by compiling modified sources is granted,
  20  * provided you
  21  *   1. distribute the corresponding source modifications from the
  22  *    released version in the form of a patch file along with the binaries,
  23  *   2. add special version identification to distinguish your version
  24  *    in addition to the base release version number,
  25  *   3. provide your name and address as the primary contact for the
  26  *    support of your modified version, and
  27  *   4. retain our contact information in regard to use of the base
  28  *    software.
  29  * Permission to distribute the released version of the source code along
  30  * with corresponding source modifications in the form of a patch file is
  31  * granted with same provisions 2 through 4 for binary distributions.
  32  *
  33  * This software is provided "as is" without express or implied warranty
  34  * to the extent permitted by applicable law.
  35 ]*/
  36
  37 /*
  38  * AUTHORS
  39  *
  40  *   Original Software:
  41  *   Jos van der Woude, jvdwoude@hut.nl
  42  *
  43  */
  44
  45 /* FIXME:
  46  * plain comparisons of floating point numbers!
  47  */
  48
  49 #include "specfun.h"
  50 #include "stdfn.h"
  51 #include "util.h"
  52
  53 #define ITMAX   200
  54
  55 #ifdef FLT_EPSILON
  56 # define MACHEPS FLT_EPSILON    /* 1.0E-08 */
  57 #else
  58 # define MACHEPS 1.0E-08
  59 #endif
  60
  61 #ifndef E_MINEXP
  62 /* AS239 value, e^-88 = 2^-127 */
  63 #define E_MINEXP  (-88.0)
  64 #endif
  65 #ifndef E_MAXEXP
  66 #define E_MAXEXP (-E_MINEXP)
  67 #endif
  68
  69 #ifdef FLT_MAX
  70 # define OFLOW   FLT_MAX                /* 1.0E+37 */
  71 #else
  72 # define OFLOW   1.0E+37
  73 #endif
  74
  75 /* AS239 value for igamma(a,x>=XBIG) = 1.0 */
  76 #define XBIG    1.0E+08
  77
  78 /*
  79  * Mathematical constants
  80  */
  81 #define LNPI 1.14472988584940016
  82 #define LNSQRT2PI 0.9189385332046727
  83 #ifdef PI
  84 # undef PI
  85 #endif
  86 #define PI 3.14159265358979323846
  87 #define PNT68 0.6796875
  88 #define SQRT_TWO 1.41421356237309504880168872420969809  /* JG */
  89
  90 /* Prefer lgamma */
  91 #ifndef GAMMA
  92 # ifdef HAVE_LGAMMA
  93 #  define GAMMA(x) lgamma (x)
  94 # elif defined(HAVE_GAMMA)
  95 #  define GAMMA(x) gamma (x)
  96 # else
  97 #  undef GAMMA
  98 # endif
  99 #endif
 100
 101 #if defined(GAMMA) && !HAVE_DECL_SIGNGAM
 102 extern int signgam;             /* this is not always declared in math.h */
 103 #endif
 104
 105 /* Local function declarations, not visible outside this file */
 106 static int mtherr __PROTO((char *, int));
 107 static double polevl __PROTO((double x, const double coef[], int N));
 108 static double p1evl __PROTO((double x, const double coef[], int N));
 109 static double confrac __PROTO((double a, double b, double x));
 110 static double ibeta __PROTO((double a, double b, double x));
 111 static double igamma __PROTO((double a, double x));
 112 static double ranf __PROTO((struct value * init));
 113 static double inverse_normal_func __PROTO((double p));
 114 static double inverse_error_func __PROTO((double p));
 115 static double lambertw __PROTO((double x));
 116 #ifndef GAMMA
 117 static int ISNAN __PROTO((double x));
 118 static int ISFINITE __PROTO((double x));
 119 static double lngamma __PROTO((double z));
 120 #endif
 121 #ifndef HAVE_ERF
 122 static double erf __PROTO((double a));
 123 #endif
 124 #ifndef HAVE_ERFC
 125 static double erfc __PROTO((double a));
 126 #endif
 127
 128 /* Macros to configure routines taken from CEPHES: */
 129
 130 /* Unknown arithmetic, invokes coefficients given in
 131  * normal decimal format.  Beware of range boundary
 132  * problems (MACHEP, MAXLOG, etc. in const.c) and
 133  * roundoff problems in pow.c:
 134  * (Sun SPARCstation)
 135  */
 136 #define UNK 1
 137
 138 /* Define to support tiny denormal numbers, else undefine. */
 139 #define DENORMAL 1
 140
 141 /* Define to ask for infinity support, else undefine. */
 142 #define INFINITIES 1
 143
 144 /* Define to ask for support of numbers that are Not-a-Number,
 145    else undefine.  This may automatically define INFINITIES in some files. */
 146 #define NANS 1
 147
 148 /* Define to distinguish between -0.0 and +0.0.  */
 149 #define MINUSZERO 1
 150
 151 /*
 152 Cephes Math Library Release 2.0:  April, 1987
 153 Copyright 1984, 1987 by Stephen L. Moshier
 154 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 155 */
 156
 157 static int merror = 0;
 158
 159 /* Notice: the order of appearance of the following messages cannot be bound
 160  * to error codes defined in mconf.h or math.h or similar, as these files are
 161  * not available on every platform. Thus, enumerate them explicitly.
 162  */
 163 #define MTHERR_DOMAIN    1
 164 #define MTHERR_SING      2
 165 #define MTHERR_OVERFLOW  3
 166 #define MTHERR_UNDERFLOW 4
 167 #define MTHERR_TLPREC    5
 168 #define MTHERR_PLPREC    6
 169
 170 static int
 171 mtherr(char *name, int code)
 172 {
 173     static const char *ermsg[7] = {
 174         "unknown",                  /* error code 0 */
 175         "domain",                   /* error code 1 */
 176         "singularity",              /* et seq.      */
 177         "overflow",
 178         "underflow",
 179         "total loss of precision",
 180         "partial loss of precision"
 181     };
 182
 183     /* Display string passed by calling program,
 184      * which is supposed to be the name of the
 185      * function in which the error occurred:
 186      */
 187     printf("\n%s ", name);
 188
 189     /* Set global error message word */
 190     merror = code;
 191
 192     /* Display error message defined by the code argument.  */
 193     if ((code <= 0) || (code >= 7))
 194         code = 0;
 195     printf("%s error\n", ermsg[code]);
 196
 197     /* Return to calling program */
 198     return (0);
 199 }
 200
 201 /*                                                      polevl.c
 202  *                                                      p1evl.c
 203  *
 204  *      Evaluate polynomial
 205  *
 206  *
 207  *
 208  * SYNOPSIS:
 209  *
 210  * int N;
 211  * double x, y, coef[N+1], polevl[];
 212  *
 213  * y = polevl( x, coef, N );
 214  *
 215  *
 216  *
 217  * DESCRIPTION:
 218  *
 219  * Evaluates polynomial of degree N:
 220  *
 221  *                     2          N
 222  * y  =  C  + C x + C x  +...+ C x
 223  *        0    1     2          N
 224  *
 225  * Coefficients are stored in reverse order:
 226  *
 227  * coef[0] = C  , ..., coef[N] = C  .
 228  *            N                   0
 229  *
 230  *  The function p1evl() assumes that coef[N] = 1.0 and is
 231  * omitted from the array.  Its calling arguments are
 232  * otherwise the same as polevl().
 233  *
 234  *
 235  * SPEED:
 236  *
 237  * In the interest of speed, there are no checks for out
 238  * of bounds arithmetic.  This routine is used by most of
 239  * the functions in the library.  Depending on available
 240  * equipment features, the user may wish to rewrite the
 241  * program in microcode or assembly language.
 242  *
 243  */
 244
 245 /*
 246 Cephes Math Library Release 2.1:  December, 1988
 247 Copyright 1984, 1987, 1988 by Stephen L. Moshier
 248 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 249 */
 250 static double
 251 polevl(double x, const double coef[], int N)
 252 {
 253     double          ans;
 254     int             i;
 255     const double    *p;
 256
 257     p = coef;
 258     ans = *p++;
 259     i = N;
 260
 261     do
 262         ans = ans * x + *p++;
 263     while (--i);
 264
 265     return (ans);
 266 }
 267
 268 /*                                          N
 269  * Evaluate polynomial when coefficient of x  is 1.0.
 270  * Otherwise same as polevl.
 271  */
 272 static double
 273 p1evl(double x, const double coef[], int N)
 274 {
 275     double              ans;
 276     const double        *p;
 277     int                 i;
 278
 279     p = coef;
 280     ans = x + *p++;
 281     i = N - 1;
 282
 283     do
 284         ans = ans * x + *p++;
 285     while (--i);
 286
 287     return (ans);
 288 }
 289
 290 #ifndef GAMMA
 291
 292 /* Provide GAMMA function for those who do not already have one */
 293
 294 int             sgngam;
 295
 296 static int
 297 ISNAN(double x)
 298 {
 299     volatile double a = x;
 300
 301     if (a != a)
 302         return 1;
 303     return 0;
 304 }
 305
 306 static int
 307 ISFINITE(double x)
 308 {
 309     volatile double a = x;
 310
 311     if (a < DBL_MAX)
 312         return 1;
 313     return 0;
 314 }
 315
 316 double
 317 lngamma(double x)
 318 {
 319     /* A[]: Stirling's formula expansion of log gamma
 320      * B[], C[]: log gamma function between 2 and 3
 321      */
 322 #ifdef UNK
 323     static const double   A[] = {
 324         8.11614167470508450300E-4,
 325         -5.95061904284301438324E-4,
 326         7.93650340457716943945E-4,
 327         -2.77777777730099687205E-3,
 328         8.33333333333331927722E-2
 329     };
 330     static const double   B[] = {
 331         -1.37825152569120859100E3,
 332         -3.88016315134637840924E4,
 333         -3.31612992738871184744E5,
 334         -1.16237097492762307383E6,
 335         -1.72173700820839662146E6,
 336         -8.53555664245765465627E5
 337     };
 338     static const double   C[] = {
 339         /* 1.00000000000000000000E0, */
 340         -3.51815701436523470549E2,
 341         -1.70642106651881159223E4,
 342         -2.20528590553854454839E5,
 343         -1.13933444367982507207E6,
 344         -2.53252307177582951285E6,
 345         -2.01889141433532773231E6
 346     };
 347     /* log( sqrt( 2*pi ) ) */
 348     static const double   LS2PI = 0.91893853320467274178;
 349 #define MAXLGM 2.556348e305
 350 #endif /* UNK */
 351
 352 #ifdef DEC
 353     static const unsigned short A[] = {
 354         0035524, 0141201, 0034633, 0031405,
 355         0135433, 0176755, 0126007, 0045030,
 356         0035520, 0006371, 0003342, 0172730,
 357         0136066, 0005540, 0132605, 0026407,
 358         0037252, 0125252, 0125252, 0125132
 359     };
 360     static const unsigned short B[] = {
 361         0142654, 0044014, 0077633, 0035410,
 362         0144027, 0110641, 0125335, 0144760,
 363         0144641, 0165637, 0142204, 0047447,
 364         0145215, 0162027, 0146246, 0155211,
 365         0145322, 0026110, 0010317, 0110130,
 366         0145120, 0061472, 0120300, 0025363
 367     };
 368     static const unsigned short C[] = {
 369         /*0040200,0000000,0000000,0000000*/
 370         0142257, 0164150, 0163630, 0112622,
 371         0143605, 0050153, 0156116, 0135272,
 372         0144527, 0056045, 0145642, 0062332,
 373         0145213, 0012063, 0106250, 0001025,
 374         0145432, 0111254, 0044577, 0115142,
 375         0145366, 0071133, 0050217, 0005122
 376     };
 377     /* log( sqrt( 2*pi ) ) */
 378     static const unsigned short LS2P[] = {040153, 037616, 041445, 0172645,};
 379 #define LS2PI *(double *)LS2P
 380 #define MAXLGM 2.035093e36
 381 #endif /* DEC */
 382
 383 #ifdef IBMPC
 384     static const unsigned short A[] = {
 385         0x6661, 0x2733, 0x9850, 0x3f4a,
 386         0xe943, 0xb580, 0x7fbd, 0xbf43,
 387         0x5ebb, 0x20dc, 0x019f, 0x3f4a,
 388         0xa5a1, 0x16b0, 0xc16c, 0xbf66,
 389         0x554b, 0x5555, 0x5555, 0x3fb5
 390     };
 391     static const unsigned short B[] = {
 392         0x6761, 0x8ff3, 0x8901, 0xc095,
 393         0xb93e, 0x355b, 0xf234, 0xc0e2,
 394         0x89e5, 0xf890, 0x3d73, 0xc114,
 395         0xdb51, 0xf994, 0xbc82, 0xc131,
 396         0xf20b, 0x0219, 0x4589, 0xc13a,
 397         0x055e, 0x5418, 0x0c67, 0xc12a
 398     };
 399     static const unsigned short C[] = {
 400         /*0x0000,0x0000,0x0000,0x3ff0,*/
 401         0x12b2, 0x1cf3, 0xfd0d, 0xc075,
 402         0xd757, 0x7b89, 0xaa0d, 0xc0d0,
 403         0x4c9b, 0xb974, 0xeb84, 0xc10a,
 404         0x0043, 0x7195, 0x6286, 0xc131,
 405         0xf34c, 0x892f, 0x5255, 0xc143,
 406         0xe14a, 0x6a11, 0xce4b, 0xc13e
 407     };
 408     /* log( sqrt( 2*pi ) ) */
 409     static const unsigned short LS2P[] = {
 410         0xbeb5, 0xc864, 0x67f1, 0x3fed
 411     };
 412 #define LS2PI *(double *)LS2P
 413 #define MAXLGM 2.556348e305
 414 #endif /* IBMPC */
 415
 416 #ifdef MIEEE
 417     static const unsigned short A[] = {
 418         0x3f4a, 0x9850, 0x2733, 0x6661,
 419         0xbf43, 0x7fbd, 0xb580, 0xe943,
 420         0x3f4a, 0x019f, 0x20dc, 0x5ebb,
 421         0xbf66, 0xc16c, 0x16b0, 0xa5a1,
 422         0x3fb5, 0x5555, 0x5555, 0x554b
 423     };
 424     static const unsigned short B[] = {
 425         0xc095, 0x8901, 0x8ff3, 0x6761,
 426         0xc0e2, 0xf234, 0x355b, 0xb93e,
 427         0xc114, 0x3d73, 0xf890, 0x89e5,
 428         0xc131, 0xbc82, 0xf994, 0xdb51,
 429         0xc13a, 0x4589, 0x0219, 0xf20b,
 430         0xc12a, 0x0c67, 0x5418, 0x055e
 431     };
 432     static const unsigned short C[] = {
 433         0xc075, 0xfd0d, 0x1cf3, 0x12b2,
 434         0xc0d0, 0xaa0d, 0x7b89, 0xd757,
 435         0xc10a, 0xeb84, 0xb974, 0x4c9b,
 436         0xc131, 0x6286, 0x7195, 0x0043,
 437         0xc143, 0x5255, 0x892f, 0xf34c,
 438         0xc13e, 0xce4b, 0x6a11, 0xe14a
 439     };
 440     /* log( sqrt( 2*pi ) ) */
 441     static const unsigned short LS2P[] = {
 442         0x3fed, 0x67f1, 0xc864, 0xbeb5
 443     };
 444 #define LS2PI *(double *)LS2P
 445 #define MAXLGM 2.556348e305
 446 #endif /* MIEEE */
 447
 448     static const double LOGPI = 1.1447298858494001741434273513530587116472948129153;
 449
 450     double          p, q, u, w, z;
 451     int             i;
 452
 453     sgngam = 1;
 454 #ifdef NANS
 455     if (ISNAN(x))
 456         return (x);
 457 #endif
 458
 459 #ifdef INFINITIES
 460     if (!ISFINITE((x)))
 461         return (DBL_MAX * DBL_MAX);
 462 #endif
 463
 464     if (x < -34.0) {
 465         q = -x;
 466         w = lngamma(q);            /* note this modifies sgngam! */
 467         p = floor(q);
 468         if (p == q) {
 469         lgsing:
 470 #ifdef INFINITIES
 471             mtherr("lngamma", MTHERR_SING);
 472             return (DBL_MAX * DBL_MAX);
 473 #else
 474             goto loverf;
 475 #endif
 476         }
 477         i = p;
 478         if ((i & 1) == 0)
 479             sgngam = -1;
 480         else
 481             sgngam = 1;
 482         z = q - p;
 483         if (z > 0.5) {
 484             p += 1.0;
 485             z = p - q;
 486         }
 487         z = q * sin(PI * z);
 488         if (z == 0.0)
 489             goto lgsing;
 490         /*      z = log(PI) - log( z ) - w;*/
 491         z = LOGPI - log(z) - w;
 492         return (z);
 493     }
 494     if (x < 13.0) {
 495         z = 1.0;
 496         p = 0.0;
 497         u = x;
 498         while (u >= 3.0) {
 499             p -= 1.0;
 500             u = x + p;
 501             z *= u;
 502         }
 503         while (u < 2.0) {
 504             if (u == 0.0)
 505                 goto lgsing;
 506             z /= u;
 507             p += 1.0;
 508             u = x + p;
 509         }
 510         if (z < 0.0) {
 511             sgngam = -1;
 512             z = -z;
 513         } else
 514             sgngam = 1;
 515         if (u == 2.0)
 516             return (log(z));
 517         p -= 2.0;
 518         x = x + p;
 519         p = x * polevl(x, B, 5) / p1evl(x, C, 6);
 520         return (log(z) + p);
 521     }
 522     if (x > MAXLGM) {
 523 #ifdef INFINITIES
 524         return (sgngam * (DBL_MAX * DBL_MAX));
 525 #else
 526     loverf:
 527         mtherr("lngamma", MTHERR_OVERFLOW);
 528         return (sgngam * MAXNUM);
 529 #endif
 530     }
 531     q = (x - 0.5) * log(x) - x + LS2PI;
 532     if (x > 1.0e8)
 533         return (q);
 534
 535     p = 1.0 / (x * x);
 536     if (x >= 1000.0)
 537         q += ((7.9365079365079365079365e-4 * p
 538                - 2.7777777777777777777778e-3) * p
 539               + 0.0833333333333333333333) / x;
 540     else
 541         q += polevl(p, A, 4) / x;
 542     return (q);
 543 }
 544
 545 #define GAMMA(x) lngamma ((x))
 546 /* HBB 20030816: must override name of sgngam so f_gamma() uses it */
 547 #define signgam sgngam
 548
 549 #endif /* !GAMMA */
 550
 551 void
 552 f_erf(union argument *arg)
 553 {
 554     struct value a;
 555     double x;
 556
 557     (void) arg;                         /* avoid -Wunused warning */
 558     x = real(pop(&a));
 559     x = erf(x);
 560     push(Gcomplex(&a, x, 0.0));
 561 }
 562
 563 void
 564 f_erfc(union argument *arg)
 565 {
 566     struct value a;
 567     double x;
 568
 569     (void) arg;                         /* avoid -Wunused warning */
 570     x = real(pop(&a));
 571     x = erfc(x);
 572     push(Gcomplex(&a, x, 0.0));
 573 }
 574
 575 void
 576 f_ibeta(union argument *arg)
 577 {
 578     struct value a;
 579     double x;
 580     double arg1;
 581     double arg2;
 582
 583     (void) arg;                         /* avoid -Wunused warning */
 584     x = real(pop(&a));
 585     arg2 = real(pop(&a));
 586     arg1 = real(pop(&a));
 587
 588     x = ibeta(arg1, arg2, x);
 589     if (x == -1.0) {
 590         undefined = TRUE;
 591         push(Ginteger(&a, 0));
 592     } else
 593         push(Gcomplex(&a, x, 0.0));
 594 }
 595
 596 void f_igamma(union argument *arg)
 597 {
 598     struct value a;
 599     double x;
 600     double arg1;
 601
 602     (void) arg;                         /* avoid -Wunused warning */
 603     x = real(pop(&a));
 604     arg1 = real(pop(&a));
 605
 606     x = igamma(arg1, x);
 607     if (x == -1.0) {
 608         undefined = TRUE;
 609         push(Ginteger(&a, 0));
 610     } else
 611         push(Gcomplex(&a, x, 0.0));
 612 }
 613
 614 void f_gamma(union argument *arg)
 615 {
 616     double y;
 617     struct value a;
 618
 619     (void) arg;                         /* avoid -Wunused warning */
 620     y = GAMMA(real(pop(&a)));
 621     if (y > E_MAXEXP) {
 622         undefined = TRUE;
 623         push(Ginteger(&a, 0));
 624     } else
 625         push(Gcomplex(&a, signgam * gp_exp(y), 0.0));
 626 }
 627
 628 void f_lgamma(union argument *arg)
 629 {
 630     struct value a;
 631
 632     (void) arg;                         /* avoid -Wunused warning */
 633     push(Gcomplex(&a, GAMMA(real(pop(&a))), 0.0));
 634 }
 635
 636 #ifndef BADRAND
 637
 638 void f_rand(union argument *arg)
 639 {
 640     struct value a;
 641
 642     (void) arg;                         /* avoid -Wunused warning */
 643     push(Gcomplex(&a, ranf(pop(&a)), 0.0));
 644 }
 645
 646 #else /* BADRAND */
 647
 648 /* Use only to observe the effect of a "bad" random number generator. */
 649 void f_rand(union argument *arg)
 650 {
 651     struct value a;
 652
 653     (void) arg;                         /* avoid -Wunused warning */
 654     static unsigned int y = 0;
 655     unsigned int maxran = 1000;
 656
 657     (void) real(pop(&a));
 658     y = (781 * y + 387) % maxran;
 659
 660     push(Gcomplex(&a, (double) y / maxran, 0.0));
 661 }
 662
 663 #endif /* BADRAND */
 664
 665 /* ** ibeta.c
 666  *
 667  *   DESCRIBE  Approximate the incomplete beta function Ix(a, b).
 668  *
 669  *                           _
 670  *                          |(a + b)     /x  (a-1)         (b-1)
 671  *             Ix(a, b) = -_-------_--- * |  t     * (1 - t)     dt (a,b > 0)
 672  *                        |(a) * |(b)   /0
 673  *
 674  *
 675  *
 676  *   CALL      p = ibeta(a, b, x)
 677  *
 678  *             double    a    > 0
 679  *             double    b    > 0
 680  *             double    x    [0, 1]
 681  *
 682  *   WARNING   none
 683  *
 684  *   RETURN    double    p    [0, 1]
 685  *                            -1.0 on error condition
 686  *
 687  *   XREF      lngamma()
 688  *
 689  *   BUGS      none
 690  *
 691  *   REFERENCE The continued fraction expansion as given by
 692  *             Abramowitz and Stegun (1964) is used.
 693  *
 694  * Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl
 695  *
 696  * Note: this function was translated from the Public Domain Fortran
 697  *       version available from http://lib.stat.cmu.edu/apstat/xxx
 698  *
 699  */
 700
 701 static double
 702 ibeta(double a, double b, double x)
 703 {
 704     /* Test for admissibility of arguments */
 705     if (a <= 0.0 || b <= 0.0)
 706         return -1.0;
 707     if (x < 0.0 || x > 1.0)
 708         return -1.0;;
 709
 710     /* If x equals 0 or 1, return x as prob */
 711     if (x == 0.0 || x == 1.0)
 712         return x;
 713
 714     /* Swap a, b if necessary for more efficient evaluation */
 715     return a < x * (a + b) ? 1.0 - confrac(b, a, 1.0 - x) : confrac(a, b, x);
 716 }
 717
 718 static double
 719 confrac(double a, double b, double x)
 720 {
 721     double Alo = 0.0;
 722     double Ahi;
 723     double Aev;
 724     double Aod;
 725     double Blo = 1.0;
 726     double Bhi = 1.0;
 727     double Bod = 1.0;
 728     double Bev = 1.0;
 729     double f;
 730     double fold;
 731     double Apb = a + b;
 732     double d;
 733     int i;
 734     int j;
 735
 736     /* Set up continued fraction expansion evaluation. */
 737     Ahi = gp_exp(GAMMA(Apb) + a * log(x) + b * log(1.0 - x) -
 738                  GAMMA(a + 1.0) - GAMMA(b));
 739
 740     /*
 741      * Continued fraction loop begins here. Evaluation continues until
 742      * maximum iterations are exceeded, or convergence achieved.
 743      */
 744     for (i = 0, j = 1, f = Ahi; i <= ITMAX; i++, j++) {
 745         d = a + j + i;
 746         Aev = -(a + i) * (Apb + i) * x / d / (d - 1.0);
 747         Aod = j * (b - j) * x / d / (d + 1.0);
 748         Alo = Bev * Ahi + Aev * Alo;
 749         Blo = Bev * Bhi + Aev * Blo;
 750         Ahi = Bod * Alo + Aod * Ahi;
 751         Bhi = Bod * Blo + Aod * Bhi;
 752
 753         if (fabs(Bhi) < MACHEPS)
 754             Bhi = 0.0;
 755
 756         if (Bhi != 0.0) {
 757             fold = f;
 758             f = Ahi / Bhi;
 759             if (fabs(f - fold) < fabs(f) * MACHEPS)
 760                 return f;
 761         }
 762     }
 763
 764     return -1.0;
 765 }
 766
 767 /* ** igamma.c
 768  *
 769  *   DESCRIBE  Approximate the incomplete gamma function P(a, x).
 770  *
 771  *                         1     /x  -t   (a-1)
 772  *             P(a, x) = -_--- * |  e  * t     dt      (a > 0)
 773  *                       |(a)   /0
 774  *
 775  *   CALL      p = igamma(a, x)
 776  *
 777  *             double    a    >  0
 778  *             double    x    >= 0
 779  *
 780  *   WARNING   none
 781  *
 782  *   RETURN    double    p    [0, 1]
 783  *                            -1.0 on error condition
 784  *
 785  *   XREF      lngamma()
 786  *
 787  *   BUGS      Values 0 <= x <= 1 may lead to inaccurate results.
 788  *
 789  *   REFERENCE ALGORITHM AS239  APPL. STATIST. (1988) VOL. 37, NO. 3
 790  *
 791  * Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl
 792  *
 793  * Note: this function was translated from the Public Domain Fortran
 794  *       version available from http://lib.stat.cmu.edu/apstat/239
 795  *
 796  */
 797
 798 static double
 799 igamma(double a, double x)
 800 {
 801     double arg;
 802     double aa;
 803     double an;
 804     double b;
 805     int i;
 806
 807     /* Check that we have valid values for a and x */
 808     if (x < 0.0 || a <= 0.0)
 809         return -1.0;
 810
 811     /* Deal with special cases */
 812     if (x == 0.0)
 813         return 0.0;
 814     if (x > XBIG)
 815         return 1.0;
 816
 817     /* Check value of factor arg */
 818     arg = a * log(x) - x - GAMMA(a + 1.0);
 819     /* HBB 20031006: removed a spurious check here */
 820     arg = gp_exp(arg);
 821
 822     /* Choose infinite series or continued fraction. */
 823
 824     if ((x > 1.0) && (x >= a + 2.0)) {
 825         /* Use a continued fraction expansion */
 826         double pn1, pn2, pn3, pn4, pn5, pn6;
 827         double rn;
 828         double rnold;
 829
 830         aa = 1.0 - a;
 831         b = aa + x + 1.0;
 832         pn1 = 1.0;
 833         pn2 = x;
 834         pn3 = x + 1.0;
 835         pn4 = x * b;
 836         rnold = pn3 / pn4;
 837
 838         for (i = 1; i <= ITMAX; i++) {
 839
 840             aa++;
 841             b += 2.0;
 842             an = aa * (double) i;
 843
 844             pn5 = b * pn3 - an * pn1;
 845             pn6 = b * pn4 - an * pn2;
 846
 847             if (pn6 != 0.0) {
 848
 849                 rn = pn5 / pn6;
 850                 if (fabs(rnold - rn) <= GPMIN(MACHEPS, MACHEPS * rn))
 851                     return 1.0 - arg * rn * a;
 852
 853                 rnold = rn;
 854             }
 855             pn1 = pn3;
 856             pn2 = pn4;
 857             pn3 = pn5;
 858             pn4 = pn6;
 859
 860             /* Re-scale terms in continued fraction if terms are large */
 861             if (fabs(pn5) >= OFLOW) {
 862
 863                 pn1 /= OFLOW;
 864                 pn2 /= OFLOW;
 865                 pn3 /= OFLOW;
 866                 pn4 /= OFLOW;
 867             }
 868         }
 869     } else {
 870         /* Use Pearson's series expansion. */
 871
 872         for (i = 0, aa = a, an = b = 1.0; i <= ITMAX; i++) {
 873
 874             aa++;
 875             an *= x / aa;
 876             b += an;
 877             if (an < b * MACHEPS)
 878                 return arg * b;
 879         }
 880     }
 881     return -1.0;
 882 }
 883
 884 /***********************************************************************
 885      double ranf(double init)
 886                 Random number generator as a Function
 887      Returns a random floating point number from a uniform distribution
 888      over 0 - 1 (endpoints of this interval are not returned) using a
 889      large integer generator.
 890      This is a transcription from Pascal to Fortran of routine
 891      Uniform_01 from the paper
 892      L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
 893      with Splitting Facilities." ACM Transactions on Mathematical
 894      Software, 17:98-111 (1991)
 895
 896                Generate Large Integer
 897      Returns a random integer following a uniform distribution over
 898      (1, 2147483562) using the generator.
 899      This is a transcription from Pascal to Fortran of routine
 900      Random from the paper
 901      L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
 902      with Splitting Facilities." ACM Transactions on Mathematical
 903      Software, 17:98-111 (1991)
 904 ***********************************************************************/
 905 static double
 906 ranf(struct value *init)
 907 {
 908     long k, z;
 909     static int firsttime = 1;
 910     static long seed1, seed2;
 911     static const long Xm1 = 2147483563L;
 912     static const long Xm2 = 2147483399L;
 913     static const long Xa1 = 40014L;
 914     static const long Xa2 = 40692L;
 915
 916
 917     /* (Re)-Initialize seeds if necessary */
 918     if ( real(init) < 0.0 || firsttime == 1) {
 919         firsttime = 0;
 920         seed1 = 1234567890L;
 921         seed2 = 1234567890L;
 922     }
 923
 924     /* Construct new seed values from input parameter */
 925     /* FIXME: Ideally we should allow all 64 bits of seed to be set */
 926     if (real(init) > 0.0) {
 927         if (real(init) >= (double)(017777777777UL))
 928             int_error(NO_CARET,"Illegal seed value");
 929         if (imag(init) >= (double)(017777777777UL))
 930             int_error(NO_CARET,"Illegal seed value");
 931         seed1 = (int)real(init);
 932         seed2 = (int)imag(init);
 933         if (seed2 == 0)
 934             seed2 = seed1;
 935     }
 936     FPRINTF((stderr,"ranf: seed = %lo %lo        %ld %ld\n", seed1,seed2));
 937
 938     /* Generate pseudo random integers */
 939     k = seed1 / 53668L;
 940     seed1 = Xa1 * (seed1 - k * 53668L) - k * 12211;
 941     if (seed1 < 0)
 942         seed1 += Xm1;
 943     k = seed2 / 52774L;
 944     seed2 = Xa2 * (seed2 - k * 52774L) - k * 3791;
 945     if (seed2 < 0)
 946         seed2 += Xm2;
 947     z = seed1 - seed2;
 948     if (z < 1)
 949         z += (Xm1 - 1);
 950
 951     /*
 952      * 4.656613057E-10 is 1/Xm1.  Xm1 is set at the top of this file and is
 953      * currently 2147483563. If Xm1 changes, change this also.
 954      */
 955     return (double) 4.656613057E-10 *z;
 956 }
 957
 958 /* ----------------------------------------------------------------
 959    Following to specfun.c made by John Grosh (jgrosh@arl.mil)
 960    on 28 OCT 1992.
 961    ---------------------------------------------------------------- */
 962
 963 void
 964 f_normal(union argument *arg)
 965 {                               /* Normal or Gaussian Probability Function */
 966     struct value a;
 967     double x;
 968
 969     /* ref. Abramowitz and Stegun 1964, "Handbook of Mathematical
 970        Functions", Applied Mathematics Series, vol 55,
 971        Chapter 26, page 934, Eqn. 26.2.29 and Jos van der Woude
 972        code found above */
 973
 974     (void) arg;                         /* avoid -Wunused warning */
 975     x = real(pop(&a));
 976
 977     x = 0.5 * SQRT_TWO * x;
 978     x = 0.5 * (1.0 + erf(x));
 979     push(Gcomplex(&a, x, 0.0));
 980 }
 981
 982 void
 983 f_inverse_normal(union argument *arg)
 984 {                               /* Inverse normal distribution function */
 985     struct value a;
 986     double x;
 987
 988     (void) arg;                 /* avoid -Wunused warning */
 989     x = real(pop(&a));
 990
 991     if (x <= 0.0 || x >= 1.0) {
 992         undefined = TRUE;
 993         push(Gcomplex(&a, 0.0, 0.0));
 994     } else {
 995         push(Gcomplex(&a, inverse_normal_func(x), 0.0));
 996     }
 997 }
 998
 999
1000 void
1001 f_inverse_erf(union argument *arg)
1002 {                               /* Inverse error function */
1003     struct value a;
1004     double x;
1005
1006     (void) arg;                         /* avoid -Wunused warning */
1007     x = real(pop(&a));
1008
1009     if (fabs(x) >= 1.0) {
1010         undefined = TRUE;
1011         push(Gcomplex(&a, 0.0, 0.0));
1012     } else {
1013         push(Gcomplex(&a, inverse_error_func(x), 0.0));
1014     }
1015 }
1016
1017 /*                                                      ndtri.c
1018  *
1019  *      Inverse of Normal distribution function
1020  *
1021  *
1022  *
1023  * SYNOPSIS:
1024  *
1025  * double x, y, ndtri();
1026  *
1027  * x = ndtri( y );
1028  *
1029  *
1030  *
1031  * DESCRIPTION:
1032  *
1033  * Returns the argument, x, for which the area under the
1034  * Gaussian probability density function (integrated from
1035  * minus infinity to x) is equal to y.
1036  *
1037  *
1038  * For small arguments 0 < y < exp(-2), the program computes
1039  * z = sqrt( -2.0 * log(y) );  then the approximation is
1040  * x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z).
1041  * There are two rational functions P/Q, one for 0 < y < exp(-32)
1042  * and the other for y up to exp(-2).  For larger arguments,
1043  * w = y - 0.5, and  x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
1044  *
1045  *
1046  * ACCURACY:
1047  *
1048  *                      Relative error:
1049  * arithmetic   domain        # trials      peak         rms
1050  *    DEC      0.125, 1         5500       9.5e-17     2.1e-17
1051  *    DEC      6e-39, 0.135     3500       5.7e-17     1.3e-17
1052  *    IEEE     0.125, 1        20000       7.2e-16     1.3e-16
1053  *    IEEE     3e-308, 0.135   50000       4.6e-16     9.8e-17
1054  *
1055  *
1056  * ERROR MESSAGES:
1057  *
1058  *   message         condition    value returned
1059  * ndtri domain       x <= 0        -DBL_MAX
1060  * ndtri domain       x >= 1         DBL_MAX
1061  *
1062  */
1063
1064 /*
1065 Cephes Math Library Release 2.8:  June, 2000
1066 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
1067 */
1068
1069 #ifdef UNK
1070 /* sqrt(2pi) */
1071 static double   s2pi = 2.50662827463100050242E0;
1072 #endif
1073
1074 #ifdef DEC
1075 static unsigned short s2p[] = {0040440, 0066230, 0177661, 0034055};
1076 #define s2pi *(double *)s2p
1077 #endif
1078
1079 #ifdef IBMPC
1080 static unsigned short s2p[] = {0x2706, 0x1ff6, 0x0d93, 0x4004};
1081 #define s2pi *(double *)s2p
1082 #endif
1083
1084 #ifdef MIEEE
1085 static unsigned short s2p[] = {
1086     0x4004, 0x0d93, 0x1ff6, 0x2706
1087 };
1088 #define s2pi *(double *)s2p
1089 #endif
1090
1091 static double
1092 inverse_normal_func(double y0)
1093 {
1094     /* approximation for 0 <= |y - 0.5| <= 3/8 */
1095 #ifdef UNK
1096     static const double   P0[5] = {
1097         -5.99633501014107895267E1,
1098         9.80010754185999661536E1,
1099         -5.66762857469070293439E1,
1100         1.39312609387279679503E1,
1101         -1.23916583867381258016E0,
1102     };
1103     static const double   Q0[8] = {
1104         /* 1.00000000000000000000E0,*/
1105         1.95448858338141759834E0,
1106         4.67627912898881538453E0,
1107         8.63602421390890590575E1,
1108         -2.25462687854119370527E2,
1109         2.00260212380060660359E2,
1110         -8.20372256168333339912E1,
1111         1.59056225126211695515E1,
1112         -1.18331621121330003142E0,
1113     };
1114 #endif
1115 #ifdef DEC
1116     static const unsigned short P0[20] = {
1117         0141557, 0155170, 0071360, 0120550,
1118         0041704, 0000214, 0172417, 0067307,
1119         0141542, 0132204, 0040066, 0156723,
1120         0041136, 0163161, 0157276, 0007747,
1121         0140236, 0116374, 0073666, 0051764,
1122     };
1123     static const unsigned short Q0[32] = {
1124         /*0040200,0000000,0000000,0000000,*/
1125         0040372, 0026256, 0110403, 0123707,
1126         0040625, 0122024, 0020277, 0026661,
1127         0041654, 0134161, 0124134, 0007244,
1128         0142141, 0073162, 0133021, 0131371,
1129         0042110, 0041235, 0043516, 0057767,
1130         0141644, 0011417, 0036155, 0137305,
1131         0041176, 0076556, 0004043, 0125430,
1132         0140227, 0073347, 0152776, 0067251,
1133     };
1134 #endif
1135 #ifdef IBMPC
1136     static const unsigned short P0[20] = {
1137         0x142d, 0x0e5e, 0xfb4f, 0xc04d,
1138         0xedd9, 0x9ea1, 0x8011, 0x4058,
1139         0xdbba, 0x8806, 0x5690, 0xc04c,
1140         0xc1fd, 0x3bd7, 0xdcce, 0x402b,
1141         0xca7e, 0x8ef6, 0xd39f, 0xbff3,
1142     };
1143     static const unsigned short Q0[36] = {
1144         /*0x0000,0x0000,0x0000,0x3ff0,*/
1145         0x74f9, 0xd220, 0x4595, 0x3fff,
1146         0xe5b6, 0x8417, 0xb482, 0x4012,
1147         0x81d4, 0x350b, 0x970e, 0x4055,
1148         0x365f, 0x56c2, 0x2ece, 0xc06c,
1149         0xcbff, 0xa8e9, 0x0853, 0x4069,
1150         0xb7d9, 0xe78d, 0x8261, 0xc054,
1151         0x7563, 0xc104, 0xcfad, 0x402f,
1152         0xcdd5, 0xfabf, 0xeedc, 0xbff2,
1153     };
1154 #endif
1155 #ifdef MIEEE
1156     static const unsigned short P0[20] = {
1157         0xc04d, 0xfb4f, 0x0e5e, 0x142d,
1158         0x4058, 0x8011, 0x9ea1, 0xedd9,
1159         0xc04c, 0x5690, 0x8806, 0xdbba,
1160         0x402b, 0xdcce, 0x3bd7, 0xc1fd,
1161         0xbff3, 0xd39f, 0x8ef6, 0xca7e,
1162     };
1163     static const unsigned short Q0[32] = {
1164         /*0x3ff0,0x0000,0x0000,0x0000,*/
1165         0x3fff, 0x4595, 0xd220, 0x74f9,
1166         0x4012, 0xb482, 0x8417, 0xe5b6,
1167         0x4055, 0x970e, 0x350b, 0x81d4,
1168         0xc06c, 0x2ece, 0x56c2, 0x365f,
1169         0x4069, 0x0853, 0xa8e9, 0xcbff,
1170         0xc054, 0x8261, 0xe78d, 0xb7d9,
1171         0x402f, 0xcfad, 0xc104, 0x7563,
1172         0xbff2, 0xeedc, 0xfabf, 0xcdd5,
1173     };
1174 #endif
1175
1176     /* Approximation for interval z = sqrt(-2 log y ) between 2 and 8
1177      * i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
1178      */
1179 #ifdef UNK
1180     static const double   P1[9] = {
1181         4.05544892305962419923E0,
1182         3.15251094599893866154E1,
1183         5.71628192246421288162E1,
1184         4.40805073893200834700E1,
1185         1.46849561928858024014E1,
1186         2.18663306850790267539E0,
1187         -1.40256079171354495875E-1,
1188         -3.50424626827848203418E-2,
1189         -8.57456785154685413611E-4,
1190     };
1191     static const double   Q1[8] = {
1192         /*  1.00000000000000000000E0,*/
1193         1.57799883256466749731E1,
1194         4.53907635128879210584E1,
1195         4.13172038254672030440E1,
1196         1.50425385692907503408E1,
1197         2.50464946208309415979E0,
1198         -1.42182922854787788574E-1,
1199         -3.80806407691578277194E-2,
1200         -9.33259480895457427372E-4,
1201     };
1202 #endif
1203 #ifdef DEC
1204     static const unsigned short P1[36] = {
1205         0040601, 0143074, 0150744, 0073326,
1206         0041374, 0031554, 0113253, 0146016,
1207         0041544, 0123272, 0012463, 0176771,
1208         0041460, 0051160, 0103560, 0156511,
1209         0041152, 0172624, 0117772, 0030755,
1210         0040413, 0170713, 0151545, 0176413,
1211         0137417, 0117512, 0022154, 0131671,
1212         0137017, 0104257, 0071432, 0007072,
1213         0135540, 0143363, 0063137, 0036166,
1214     };
1215     static const unsigned short Q1[32] = {
1216         /*0040200,0000000,0000000,0000000,*/
1217         0041174, 0075325, 0004736, 0120326,
1218         0041465, 0110044, 0047561, 0045567,
1219         0041445, 0042321, 0012142, 0030340,
1220         0041160, 0127074, 0166076, 0141051,
1221         0040440, 0046055, 0040745, 0150400,
1222         0137421, 0114146, 0067330, 0010621,
1223         0137033, 0175162, 0025555, 0114351,
1224         0135564, 0122773, 0145750, 0030357,
1225     };
1226 #endif
1227 #ifdef IBMPC
1228     static const unsigned short P1[36] = {
1229         0x8edb, 0x9a3c, 0x38c7, 0x4010,
1230         0x7982, 0x92d5, 0x866d, 0x403f,
1231         0x7fbf, 0x42a6, 0x94d7, 0x404c,
1232         0x1ba9, 0x10ee, 0x0a4e, 0x4046,
1233         0x463e, 0x93ff, 0x5eb2, 0x402d,
1234         0xbfa1, 0x7a6c, 0x7e39, 0x4001,
1235         0x9677, 0x448d, 0xf3e9, 0xbfc1,
1236         0x41c7, 0xee63, 0xf115, 0xbfa1,
1237         0xe78f, 0x6ccb, 0x18de, 0xbf4c,
1238     };
1239     static const unsigned short Q1[32] = {
1240         /*0x0000,0x0000,0x0000,0x3ff0,*/
1241         0xd41b, 0xa13b, 0x8f5a, 0x402f,
1242         0x296f, 0x89ee, 0xb204, 0x4046,
1243         0x461c, 0x228c, 0xa89a, 0x4044,
1244         0xd845, 0x9d87, 0x15c7, 0x402e,
1245         0xba20, 0xa83c, 0x0985, 0x4004,
1246         0x0232, 0xcddb, 0x330c, 0xbfc2,
1247         0xb31d, 0x456d, 0x7f4e, 0xbfa3,
1248         0x061e, 0x797d, 0x94bf, 0xbf4e,
1249     };
1250 #endif
1251 #ifdef MIEEE
1252     static const unsigned short P1[36] = {
1253         0x4010, 0x38c7, 0x9a3c, 0x8edb,
1254         0x403f, 0x866d, 0x92d5, 0x7982,
1255         0x404c, 0x94d7, 0x42a6, 0x7fbf,
1256         0x4046, 0x0a4e, 0x10ee, 0x1ba9,
1257         0x402d, 0x5eb2, 0x93ff, 0x463e,
1258         0x4001, 0x7e39, 0x7a6c, 0xbfa1,
1259         0xbfc1, 0xf3e9, 0x448d, 0x9677,
1260         0xbfa1, 0xf115, 0xee63, 0x41c7,
1261         0xbf4c, 0x18de, 0x6ccb, 0xe78f,
1262     };
1263     static const unsigned short Q1[32] = {
1264         /*0x3ff0,0x0000,0x0000,0x0000,*/
1265         0x402f, 0x8f5a, 0xa13b, 0xd41b,
1266         0x4046, 0xb204, 0x89ee, 0x296f,
1267         0x4044, 0xa89a, 0x228c, 0x461c,
1268         0x402e, 0x15c7, 0x9d87, 0xd845,
1269         0x4004, 0x0985, 0xa83c, 0xba20,
1270         0xbfc2, 0x330c, 0xcddb, 0x0232,
1271         0xbfa3, 0x7f4e, 0x456d, 0xb31d,
1272         0xbf4e, 0x94bf, 0x797d, 0x061e,
1273     };
1274 #endif
1275
1276     /* Approximation for interval z = sqrt(-2 log y ) between 8 and 64
1277      * i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
1278      */
1279
1280 #ifdef UNK
1281     static const double   P2[9] = {
1282         3.23774891776946035970E0,
1283         6.91522889068984211695E0,
1284         3.93881025292474443415E0,
1285         1.33303460815807542389E0,
1286         2.01485389549179081538E-1,
1287         1.23716634817820021358E-2,
1288         3.01581553508235416007E-4,
1289         2.65806974686737550832E-6,
1290         6.23974539184983293730E-9,
1291     };
1292     static const double   Q2[8] = {
1293         /*  1.00000000000000000000E0,*/
1294         6.02427039364742014255E0,
1295         3.67983563856160859403E0,
1296         1.37702099489081330271E0,
1297         2.16236993594496635890E-1,
1298         1.34204006088543189037E-2,
1299         3.28014464682127739104E-4,
1300         2.89247864745380683936E-6,
1301         6.79019408009981274425E-9,
1302     };
1303 #endif
1304 #ifdef DEC
1305     static const unsigned short P2[36] = {
1306         0040517, 0033507, 0036236, 0125641,
1307         0040735, 0044616, 0014473, 0140133,
1308         0040574, 0012567, 0114535, 0102541,
1309         0040252, 0120340, 0143474, 0150135,
1310         0037516, 0051057, 0115361, 0031211,
1311         0036512, 0131204, 0101511, 0125144,
1312         0035236, 0016627, 0043160, 0140216,
1313         0033462, 0060512, 0060141, 0010641,
1314         0031326, 0062541, 0101304, 0077706,
1315     };
1316     static const unsigned short Q2[32] = {
1317         /*0040200,0000000,0000000,0000000,*/
1318         0040700, 0143322, 0132137, 0040501,
1319         0040553, 0101155, 0053221, 0140257,
1320         0040260, 0041071, 0052573, 0010004,
1321         0037535, 0066472, 0177261, 0162330,
1322         0036533, 0160475, 0066666, 0036132,
1323         0035253, 0174533, 0027771, 0044027,
1324         0033502, 0016147, 0117666, 0063671,
1325         0031351, 0047455, 0141663, 0054751,
1326     };
1327 #endif
1328 #ifdef IBMPC
1329     static const unsigned short P2[36] = {
1330         0xd574, 0xe793, 0xe6e8, 0x4009,
1331         0x780b, 0xc327, 0xa931, 0x401b,
1332         0xb0ac, 0xf32b, 0x82ae, 0x400f,
1333         0x9a0c, 0x18e7, 0x541c, 0x3ff5,
1334         0x2651, 0xf35e, 0xca45, 0x3fc9,
1335         0x354d, 0x9069, 0x5650, 0x3f89,
1336         0x1812, 0xe8ce, 0xc3b2, 0x3f33,
1337         0x2234, 0x4c0c, 0x4c29, 0x3ec6,
1338         0x8ff9, 0x3058, 0xccac, 0x3e3a,
1339     };
1340     static const unsigned short Q2[32] = {
1341         /*0x0000,0x0000,0x0000,0x3ff0,*/
1342         0xe828, 0x568b, 0x18da, 0x4018,
1343         0x3816, 0xaad2, 0x704d, 0x400d,
1344         0x6200, 0x2aaf, 0x0847, 0x3ff6,
1345         0x3c9b, 0x5fd6, 0xada7, 0x3fcb,
1346         0xc78b, 0xadb6, 0x7c27, 0x3f8b,
1347         0x2903, 0x65ff, 0x7f2b, 0x3f35,
1348         0xccf7, 0xf3f6, 0x438c, 0x3ec8,
1349         0x6b3d, 0xb876, 0x29e5, 0x3e3d,
1350     };
1351 #endif
1352 #ifdef MIEEE
1353     static const unsigned short P2[36] = {
1354         0x4009, 0xe6e8, 0xe793, 0xd574,
1355         0x401b, 0xa931, 0xc327, 0x780b,
1356         0x400f, 0x82ae, 0xf32b, 0xb0ac,
1357         0x3ff5, 0x541c, 0x18e7, 0x9a0c,
1358         0x3fc9, 0xca45, 0xf35e, 0x2651,
1359         0x3f89, 0x5650, 0x9069, 0x354d,
1360         0x3f33, 0xc3b2, 0xe8ce, 0x1812,
1361         0x3ec6, 0x4c29, 0x4c0c, 0x2234,
1362         0x3e3a, 0xccac, 0x3058, 0x8ff9,
1363     };
1364     static const unsigned short Q2[32] = {
1365         /*0x3ff0,0x0000,0x0000,0x0000,*/
1366         0x4018, 0x18da, 0x568b, 0xe828,
1367         0x400d, 0x704d, 0xaad2, 0x3816,
1368         0x3ff6, 0x0847, 0x2aaf, 0x6200,
1369         0x3fcb, 0xada7, 0x5fd6, 0x3c9b,
1370         0x3f8b, 0x7c27, 0xadb6, 0xc78b,
1371         0x3f35, 0x7f2b, 0x65ff, 0x2903,
1372         0x3ec8, 0x438c, 0xf3f6, 0xccf7,
1373         0x3e3d, 0x29e5, 0xb876, 0x6b3d,
1374     };
1375 #endif
1376
1377     double          x, y, z, y2, x0, x1;
1378     int             code;
1379
1380     if (y0 <= 0.0) {
1381         mtherr("inverse_normal_func", MTHERR_DOMAIN);
1382         return (-DBL_MAX);
1383     }
1384     if (y0 >= 1.0) {
1385         mtherr("inverse_normal_func", MTHERR_DOMAIN);
1386         return (DBL_MAX);
1387     }
1388     code = 1;
1389     y = y0;
1390     if (y > (1.0 - 0.13533528323661269189)) {   /* 0.135... = exp(-2) */
1391         y = 1.0 - y;
1392         code = 0;
1393     }
1394     if (y > 0.13533528323661269189) {
1395         y = y - 0.5;
1396         y2 = y * y;
1397         x = y + y * (y2 * polevl(y2, P0, 4) / p1evl(y2, Q0, 8));
1398         x = x * s2pi;
1399         return (x);
1400     }
1401     x = sqrt(-2.0 * log(y));
1402     x0 = x - log(x) / x;
1403
1404     z = 1.0 / x;
1405     if (x < 8.0)                /* y > exp(-32) = 1.2664165549e-14 */
1406         x1 = z * polevl(z, P1, 8) / p1evl(z, Q1, 8);
1407     else
1408         x1 = z * polevl(z, P2, 8) / p1evl(z, Q2, 8);
1409     x = x0 - x1;
1410     if (code != 0)
1411         x = -x;
1412     return (x);
1413 }
1414
1415 /*
1416 Cephes Math Library Release 2.8:  June, 2000
1417 Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
1418 */
1419
1420 #ifndef HAVE_ERFC
1421 /*                                                     erfc.c
1422  *
1423  *      Complementary error function
1424  *
1425  *
1426  *
1427  * SYNOPSIS:
1428  *
1429  * double x, y, erfc();
1430  *
1431  * y = erfc( x );
1432  *
1433  *
1434  *
1435  * DESCRIPTION:
1436  *
1437  *
1438  *  1 - erf(x) =
1439  *
1440  *                           inf.
1441  *                             -
1442  *                  2         | |          2
1443  *   erfc(x)  =  --------     |    exp( - t  ) dt
1444  *               sqrt(pi)   | |
1445  *                           -
1446  *                            x
1447  *
1448  *
1449  * For small x, erfc(x) = 1 - erf(x); otherwise rational
1450  * approximations are computed.
1451  *
1452  *
1453  *
1454  * ACCURACY:
1455  *
1456  *                      Relative error:
1457  * arithmetic   domain     # trials      peak         rms
1458  *    DEC       0, 9.2319   12000       5.1e-16     1.2e-16
1459  *    IEEE      0,26.6417   30000       5.7e-14     1.5e-14
1460  *
1461  *
1462  * ERROR MESSAGES:
1463  *
1464  *   message         condition              value returned
1465  * erfc underflow    x > 9.231948545 (DEC)       0.0
1466  *
1467  *
1468  */
1469
1470 static double
1471 erfc(double a)
1472 {
1473 #ifdef UNK
1474     static const double   P[] = {
1475         2.46196981473530512524E-10,
1476         5.64189564831068821977E-1,
1477         7.46321056442269912687E0,
1478         4.86371970985681366614E1,
1479         1.96520832956077098242E2,
1480         5.26445194995477358631E2,
1481         9.34528527171957607540E2,
1482         1.02755188689515710272E3,
1483         5.57535335369399327526E2
1484     };
1485     static const double   Q[] = {
1486         /* 1.00000000000000000000E0,*/
1487         1.32281951154744992508E1,
1488         8.67072140885989742329E1,
1489         3.54937778887819891062E2,
1490         9.75708501743205489753E2,
1491         1.82390916687909736289E3,
1492         2.24633760818710981792E3,
1493         1.65666309194161350182E3,
1494         5.57535340817727675546E2
1495     };
1496     static const double   R[] = {
1497         5.64189583547755073984E-1,
1498         1.27536670759978104416E0,
1499         5.01905042251180477414E0,
1500         6.16021097993053585195E0,
1501         7.40974269950448939160E0,
1502         2.97886665372100240670E0
1503     };
1504     static const double   S[] = {
1505         /* 1.00000000000000000000E0,*/
1506         2.26052863220117276590E0,
1507         9.39603524938001434673E0,
1508         1.20489539808096656605E1,
1509         1.70814450747565897222E1,
1510         9.60896809063285878198E0,
1511         3.36907645100081516050E0
1512     };
1513 #endif /* UNK */
1514
1515 #ifdef DEC
1516     static const unsigned short P[] = {
1517         0030207, 0054445, 0011173, 0021706,
1518         0040020, 0067272, 0030661, 0122075,
1519         0040756, 0151236, 0173053, 0067042,
1520         0041502, 0106175, 0062555, 0151457,
1521         0042104, 0102525, 0047401, 0003667,
1522         0042403, 0116176, 0011446, 0075303,
1523         0042551, 0120723, 0061641, 0123275,
1524         0042600, 0070651, 0007264, 0134516,
1525         0042413, 0061102, 0167507, 0176625
1526     };
1527     static const unsigned short Q[] = {
1528         /*0040200,0000000,0000000,0000000,*/
1529         0041123, 0123257, 0165741, 0017142,
1530         0041655, 0065027, 0173413, 0115450,
1531         0042261, 0074011, 0021573, 0004150,
1532         0042563, 0166530, 0013662, 0007200,
1533         0042743, 0176427, 0162443, 0105214,
1534         0043014, 0062546, 0153727, 0123772,
1535         0042717, 0012470, 0006227, 0067424,
1536         0042413, 0061103, 0003042, 0013254
1537     };
1538     static const unsigned short R[] = {
1539         0040020, 0067272, 0101024, 0155421,
1540         0040243, 0037467, 0056706, 0026462,
1541         0040640, 0116017, 0120665, 0034315,
1542         0040705, 0020162, 0143350, 0060137,
1543         0040755, 0016234, 0134304, 0130157,
1544         0040476, 0122700, 0051070, 0015473
1545     };
1546     static const unsigned short S[] = {
1547         /*0040200,0000000,0000000,0000000,*/
1548         0040420, 0126200, 0044276, 0070413,
1549         0041026, 0053051, 0007302, 0063746,
1550         0041100, 0144203, 0174051, 0061151,
1551         0041210, 0123314, 0126343, 0177646,
1552         0041031, 0137125, 0051431, 0033011,
1553         0040527, 0117362, 0152661, 0066201
1554     };
1555 #endif /* DEC */
1556
1557 #ifdef IBMPC
1558     static const unsigned short P[] = {
1559         0x6479, 0xa24f, 0xeb24, 0x3df0,
1560         0x3488, 0x4636, 0x0dd7, 0x3fe2,
1561         0x6dc4, 0xdec5, 0xda53, 0x401d,
1562         0xba66, 0xacad, 0x518f, 0x4048,
1563         0x20f7, 0xa9e0, 0x90aa, 0x4068,
1564         0xcf58, 0xc264, 0x738f, 0x4080,
1565         0x34d8, 0x6c74, 0x343a, 0x408d,
1566         0x972a, 0x21d6, 0x0e35, 0x4090,
1567         0xffb3, 0x5de8, 0x6c48, 0x4081
1568     };
1569     static const unsigned short Q[] = {
1570         /*0x0000,0x0000,0x0000,0x3ff0,*/
1571         0x23cc, 0xfd7c, 0x74d5, 0x402a,
1572         0x7365, 0xfee1, 0xad42, 0x4055,
1573         0x610d, 0x246f, 0x2f01, 0x4076,
1574         0x41d0, 0x02f6, 0x7dab, 0x408e,
1575         0x7151, 0xfca4, 0x7fa2, 0x409c,
1576         0xf4ff, 0xdafa, 0x8cac, 0x40a1,
1577         0xede2, 0x0192, 0xe2a7, 0x4099,
1578         0x42d6, 0x60c4, 0x6c48, 0x4081
1579     };
1580     static const unsigned short R[] = {
1581         0x9b62, 0x5042, 0x0dd7, 0x3fe2,
1582         0xc5a6, 0xebb8, 0x67e6, 0x3ff4,
1583         0xa71a, 0xf436, 0x1381, 0x4014,
1584         0x0c0c, 0x58dd, 0xa40e, 0x4018,
1585         0x960e, 0x9718, 0xa393, 0x401d,
1586         0x0367, 0x0a47, 0xd4b8, 0x4007
1587     };
1588     static const unsigned short S[] = {
1589         /*0x0000,0x0000,0x0000,0x3ff0,*/
1590         0xce21, 0x0917, 0x1590, 0x4002,
1591         0x4cfd, 0x21d8, 0xcac5, 0x4022,
1592         0x2c4d, 0x7f05, 0x1910, 0x4028,
1593         0x7ff5, 0x959c, 0x14d9, 0x4031,
1594         0x26c1, 0xaa63, 0x37ca, 0x4023,
1595         0x2d90, 0x5ab6, 0xf3de, 0x400a
1596     };
1597 #endif /* IBMPC */
1598
1599 #ifdef MIEEE
1600     static const unsigned short P[] = {
1601         0x3df0, 0xeb24, 0xa24f, 0x6479,
1602         0x3fe2, 0x0dd7, 0x4636, 0x3488,
1603         0x401d, 0xda53, 0xdec5, 0x6dc4,
1604         0x4048, 0x518f, 0xacad, 0xba66,
1605         0x4068, 0x90aa, 0xa9e0, 0x20f7,
1606         0x4080, 0x738f, 0xc264, 0xcf58,
1607         0x408d, 0x343a, 0x6c74, 0x34d8,
1608         0x4090, 0x0e35, 0x21d6, 0x972a,
1609         0x4081, 0x6c48, 0x5de8, 0xffb3
1610     };
1611     static const unsigned short Q[] = {
1612         0x402a, 0x74d5, 0xfd7c, 0x23cc,
1613         0x4055, 0xad42, 0xfee1, 0x7365,
1614         0x4076, 0x2f01, 0x246f, 0x610d,
1615         0x408e, 0x7dab, 0x02f6, 0x41d0,
1616         0x409c, 0x7fa2, 0xfca4, 0x7151,
1617         0x40a1, 0x8cac, 0xdafa, 0xf4ff,
1618         0x4099, 0xe2a7, 0x0192, 0xede2,
1619         0x4081, 0x6c48, 0x60c4, 0x42d6
1620     };
1621     static const unsigned short R[] = {
1622         0x3fe2, 0x0dd7, 0x5042, 0x9b62,
1623         0x3ff4, 0x67e6, 0xebb8, 0xc5a6,
1624         0x4014, 0x1381, 0xf436, 0xa71a,
1625         0x4018, 0xa40e, 0x58dd, 0x0c0c,
1626         0x401d, 0xa393, 0x9718, 0x960e,
1627         0x4007, 0xd4b8, 0x0a47, 0x0367
1628     };
1629     static const unsigned short S[] = {
1630         0x4002, 0x1590, 0x0917, 0xce21,
1631         0x4022, 0xcac5, 0x21d8, 0x4cfd,
1632         0x4028, 0x1910, 0x7f05, 0x2c4d,
1633         0x4031, 0x14d9, 0x959c, 0x7ff5,
1634         0x4023, 0x37ca, 0xaa63, 0x26c1,
1635         0x400a, 0xf3de, 0x5ab6, 0x2d90
1636     };
1637 #endif /* MIEEE */
1638
1639     double p, q, x, y, z;
1640
1641     if (a < 0.0)
1642         x = -a;
1643     else
1644         x = a;
1645
1646     if (x < 1.0)
1647         return (1.0 - erf(a));
1648
1649     z = -a * a;
1650
1651     if (z < DBL_MIN_10_EXP) {
1652     under:
1653         mtherr("erfc", MTHERR_UNDERFLOW);
1654         if (a < 0)
1655             return (2.0);
1656         else
1657             return (0.0);
1658     }
1659     z = exp(z);
1660
1661     if (x < 8.0) {
1662         p = polevl(x, P, 8);
1663         q = p1evl(x, Q, 8);
1664     } else {
1665         p = polevl(x, R, 5);
1666         q = p1evl(x, S, 6);
1667     }
1668     y = (z * p) / q;
1669
1670     if (a < 0)
1671         y = 2.0 - y;
1672
1673     if (y == 0.0)
1674         goto under;
1675
1676     return (y);
1677 }
1678 #endif /* !HAVE_ERFC */
1679
1680 #ifndef HAVE_ERF
1681 /*                                                     erf.c
1682  *
1683  *      Error function
1684  *
1685  *
1686  *
1687  * SYNOPSIS:
1688  *
1689  * double x, y, erf();
1690  *
1691  * y = erf( x );
1692  *
1693  *
1694  *
1695  * DESCRIPTION:
1696  *
1697  * The integral is
1698  *
1699  *                           x
1700  *                            -
1701  *                 2         | |          2
1702  *   erf(x)  =  --------     |    exp( - t  ) dt.
1703  *              sqrt(pi)   | |
1704  *                          -
1705  *                           0
1706  *
1707  * The magnitude of x is limited to 9.231948545 for DEC
1708  * arithmetic; 1 or -1 is returned outside this range.
1709  *
1710  * For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
1711  * erf(x) = 1 - erfc(x).
1712  *
1713  *
1714  *
1715  * ACCURACY:
1716  *
1717  *                      Relative error:
1718  * arithmetic   domain     # trials      peak         rms
1719  *    DEC       0,1         14000       4.7e-17     1.5e-17
1720  *    IEEE      0,1         30000       3.7e-16     1.0e-16
1721  *
1722  */
1723
1724 static double
1725 erf(double x)
1726 {
1727
1728 # ifdef UNK
1729     static const double   T[] = {
1730         9.60497373987051638749E0,
1731         9.00260197203842689217E1,
1732         2.23200534594684319226E3,
1733         7.00332514112805075473E3,
1734         5.55923013010394962768E4
1735     };
1736     static const double   U[] = {
1737         /* 1.00000000000000000000E0,*/
1738         3.35617141647503099647E1,
1739         5.21357949780152679795E2,
1740         4.59432382970980127987E3,
1741         2.26290000613890934246E4,
1742         4.92673942608635921086E4
1743     };
1744 # endif
1745
1746 # ifdef DEC
1747     static const unsigned short T[] = {
1748         0041031, 0126770, 0170672, 0166101,
1749         0041664, 0006522, 0072360, 0031770,
1750         0043013, 0100025, 0162641, 0126671,
1751         0043332, 0155231, 0161627, 0076200,
1752         0044131, 0024115, 0021020, 0117343
1753     };
1754     static const unsigned short U[] = {
1755         /*0040200,0000000,0000000,0000000,*/
1756         0041406, 0037461, 0177575, 0032714,
1757         0042402, 0053350, 0123061, 0153557,
1758         0043217, 0111227, 0032007, 0164217,
1759         0043660, 0145000, 0004013, 0160114,
1760         0044100, 0071544, 0167107, 0125471
1761     };
1762 # endif
1763
1764 # ifdef IBMPC
1765     static const unsigned short T[] = {
1766         0x5d88, 0x1e37, 0x35bf, 0x4023,
1767         0x067f, 0x4e9e, 0x81aa, 0x4056,
1768         0x35b7, 0xbcb4, 0x7002, 0x40a1,
1769         0xef90, 0x3c72, 0x5b53, 0x40bb,
1770         0x13dc, 0xa442, 0x2509, 0x40eb
1771     };
1772     static const unsigned short U[] = {
1773         /*0x0000,0x0000,0x0000,0x3ff0,*/
1774         0xa6ba, 0x3fef, 0xc7e6, 0x4040,
1775         0x3aee, 0x14c6, 0x4add, 0x4080,
1776         0xfd12, 0xe680, 0xf252, 0x40b1,
1777         0x7c0a, 0x0101, 0x1940, 0x40d6,
1778         0xf567, 0x9dc8, 0x0e6c, 0x40e8
1779     };
1780 # endif
1781
1782 # ifdef MIEEE
1783     static const unsigned short T[] = {
1784         0x4023, 0x35bf, 0x1e37, 0x5d88,
1785         0x4056, 0x81aa, 0x4e9e, 0x067f,
1786         0x40a1, 0x7002, 0xbcb4, 0x35b7,
1787         0x40bb, 0x5b53, 0x3c72, 0xef90,
1788         0x40eb, 0x2509, 0xa442, 0x13dc
1789     };
1790     static const unsigned short U[] = {
1791         0x4040, 0xc7e6, 0x3fef, 0xa6ba,
1792         0x4080, 0x4add, 0x14c6, 0x3aee,
1793         0x40b1, 0xf252, 0xe680, 0xfd12,
1794         0x40d6, 0x1940, 0x0101, 0x7c0a,
1795         0x40e8, 0x0e6c, 0x9dc8, 0xf567
1796     };
1797 # endif
1798
1799     double y, z;
1800
1801     if (fabs(x) > 1.0)
1802         return (1.0 - erfc(x));
1803     z = x * x;
1804     y = x * polevl(z, T, 4) / p1evl(z, U, 5);
1805     return (y);
1806 }
1807 #endif /* !HAVE_ERF */
1808
1809 /* ----------------------------------------------------------------
1810    Following function for the inverse error function is taken from
1811    NIST on 16. May 2002.
1812    Use Newton-Raphson correction also for range -1 to -y0 and
1813    add 3rd cycle to improve convergence -  E A Merritt 21.10.2003
1814    ----------------------------------------------------------------
1815  */
1816
1817 static double
1818 inverse_error_func(double y)
1819 {
1820     double x = 0.0;    /* The output */
1821     double z = 0.0;    /* Intermadiate variable */
1822     double y0 = 0.7;   /* Central range variable */
1823
1824     /* Coefficients in rational approximations. */
1825     static const double a[4] = {
1826         0.886226899, -1.645349621, 0.914624893, -0.140543331
1827     };
1828     static const double b[4] = {
1829         -2.118377725, 1.442710462, -0.329097515, 0.012229801
1830     };
1831     static const double c[4] = {
1832         -1.970840454, -1.624906493, 3.429567803, 1.641345311
1833     };
1834     static const double d[2] = {
1835         3.543889200, 1.637067800
1836     };
1837
1838     if ((y < -1.0) || (1.0 < y)) {
1839         printf("inverse_error_func: The value out of the range of the function");
1840         x = log(-1.0);
1841         return (x);
1842     } else if ((y == -1.0) || (1.0 == y)) {
1843         x = -y * log(0.0);
1844         return (x);
1845     } else if ((-1.0 < y) && (y < -y0)) {
1846         z = sqrt(-log((1.0 + y) / 2.0));
1847         x = -(((c[3] * z + c[2]) * z + c[1]) * z + c[0]) / ((d[1] * z + d[0]) * z + 1.0);
1848     } else {
1849         if ((-y0 <= y) && (y <= y0)) {
1850             z = y * y;
1851             x = y * (((a[3] * z + a[2]) * z + a[1]) * z + a[0]) /
1852                 ((((b[3] * z + b[3]) * z + b[1]) * z + b[0]) * z + 1.0);
1853         } else if ((y0 < y) && (y < 1.0)) {
1854             z = sqrt(-log((1.0 - y) / 2.0));
1855             x = (((c[3] * z + c[2]) * z + c[1]) * z + c[0]) / ((d[1] * z + d[0]) * z + 1.0);
1856         }
1857     }
1858     /* Three steps of Newton-Raphson correction to full accuracy. OK - four */
1859     x = x - (erf(x) - y) / (2.0 / sqrt(PI) * gp_exp(-x * x));
1860     x = x - (erf(x) - y) / (2.0 / sqrt(PI) * gp_exp(-x * x));
1861     x = x - (erf(x) - y) / (2.0 / sqrt(PI) * gp_exp(-x * x));
1862     x = x - (erf(x) - y) / (2.0 / sqrt(PI) * gp_exp(-x * x));
1863
1864     return (x);
1865 }
1866
1867
1868 /* Implementation of Lamberts W-function which is defined as
1869  * w(x)*e^(w(x))=x
1870  * Implementation by Gunter Kuhnle, gk@uni-leipzig.de
1871  * Algorithm originally developed by
1872  * KEITH BRIGGS, DEPARTMENT OF PLANT SCIENCES,
1873  * e-mail:kmb28@cam.ac.uk
1874  * http://epidem13.plantsci.cam.ac.uk/~kbriggs/W-ology.html */
1875
1876 static double
1877 lambertw(double x)
1878 {
1879     double p, e, t, w, eps;
1880     int i;
1881
1882     eps = MACHEPS;
1883
1884     if (x < -exp(-1))
1885         return -1;              /* error, value undefined */
1886
1887     if (fabs(x) <= eps)
1888         return x;
1889
1890     if (x < 1) {
1891         p = sqrt(2.0 * (exp(1.0) * x + 1.0));
1892         w = -1.0 + p - p * p / 3.0 + 11.0 / 72.0 * p * p * p;
1893     } else {
1894         w = log(x);
1895     }
1896
1897     if (x > 3) {
1898         w = w - log(w);
1899     }
1900     for (i = 0; i < 20; i++) {
1901         e = gp_exp(w);
1902         t = w * e - x;
1903         t = t / (e * (w + 1.0) - 0.5 * (w + 2.0) * t / (w + 1.0));
1904         w = w - t;
1905         if (fabs(t) < eps * (1.0 + fabs(w)))
1906             return w;
1907     }
1908     return -1;                 /* error: iteration didn't converge */
1909 }
1910
1911 void
1912 f_lambertw(union argument *arg)
1913 {
1914     struct value a;
1915     double x;
1916
1917     (void) arg;                        /* avoid -Wunused warning */
1918     x = real(pop(&a));
1919
1920     x = lambertw(x);
1921     if (x <= -1)
1922         /* Error return from lambertw --> flag 'undefined' */
1923         undefined = TRUE;
1924
1925     push(Gcomplex(&a, x, 0.0));
1926 }
1927