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