Annotation of OpenXM_contrib/gnuplot/specfun.c, Revision 1.1
1.1 ! maekawa 1: #ifndef lint
! 2: static char *RCSid = "$Id: specfun.c,v 1.21 1998/04/14 00:16:20 drd Exp $";
! 3: #endif
! 4:
! 5:
! 6: /* GNUPLOT - specfun.c */
! 7:
! 8: /*[
! 9: * Copyright 1986 - 1993, 1998 Thomas Williams, Colin Kelley
! 10: *
! 11: * Permission to use, copy, and distribute this software and its
! 12: * documentation for any purpose with or without fee is hereby granted,
! 13: * provided that the above copyright notice appear in all copies and
! 14: * that both that copyright notice and this permission notice appear
! 15: * in supporting documentation.
! 16: *
! 17: * Permission to modify the software is granted, but not the right to
! 18: * distribute the complete modified source code. Modifications are to
! 19: * be distributed as patches to the released version. Permission to
! 20: * distribute binaries produced by compiling modified sources is granted,
! 21: * provided you
! 22: * 1. distribute the corresponding source modifications from the
! 23: * released version in the form of a patch file along with the binaries,
! 24: * 2. add special version identification to distinguish your version
! 25: * in addition to the base release version number,
! 26: * 3. provide your name and address as the primary contact for the
! 27: * support of your modified version, and
! 28: * 4. retain our contact information in regard to use of the base
! 29: * software.
! 30: * Permission to distribute the released version of the source code along
! 31: * with corresponding source modifications in the form of a patch file is
! 32: * granted with same provisions 2 through 4 for binary distributions.
! 33: *
! 34: * This software is provided "as is" without express or implied warranty
! 35: * to the extent permitted by applicable law.
! 36: ]*/
! 37:
! 38: /*
! 39: * AUTHORS
! 40: *
! 41: * Original Software:
! 42: * Jos van der Woude, jvdwoude@hut.nl
! 43: *
! 44: */
! 45:
! 46: #include "plot.h"
! 47: #include "fnproto.h"
! 48:
! 49:
! 50: extern struct value stack[STACK_DEPTH];
! 51: extern int s_p;
! 52: extern double zero;
! 53:
! 54: #define ITMAX 200
! 55:
! 56: #ifdef FLT_EPSILON
! 57: # define MACHEPS FLT_EPSILON /* 1.0E-08 */
! 58: #else
! 59: # define MACHEPS 1.0E-08
! 60: #endif
! 61:
! 62: /* AS239 value, e^-88 = 2^-127 */
! 63: #define MINEXP -88.0
! 64:
! 65: #ifdef FLT_MAX
! 66: # define OFLOW FLT_MAX /* 1.0E+37 */
! 67: #else
! 68: # define OFLOW 1.0E+37
! 69: #endif
! 70:
! 71: /* AS239 value for igamma(a,x>=XBIG) = 1.0 */
! 72: #define XBIG 1.0E+08
! 73:
! 74: /*
! 75: * Mathematical constants
! 76: */
! 77: #define LNPI 1.14472988584940016
! 78: #define LNSQRT2PI 0.9189385332046727
! 79: #ifdef PI
! 80: # undef PI
! 81: #endif
! 82: #define PI 3.14159265358979323846
! 83: #define PNT68 0.6796875
! 84: #define SQRT_TWO 1.41421356237309504880168872420969809 /* JG */
! 85:
! 86: /* Prefer lgamma */
! 87: #ifndef GAMMA
! 88: # ifdef HAVE_LGAMMA
! 89: # define GAMMA(x) lgamma (x)
! 90: # elif defined(HAVE_GAMMA)
! 91: # define GAMMA(x) gamma (x)
! 92: # else
! 93: # undef GAMMA
! 94: # endif
! 95: #endif
! 96:
! 97: #ifndef GAMMA
! 98: int signgam = 0;
! 99: #else
! 100: extern int signgam; /* this is not always declared in math.h */
! 101: #endif
! 102:
! 103: /* Global variables, not visible outside this file */
! 104: static long Xm1 = 2147483563L;
! 105: static long Xm2 = 2147483399L;
! 106: static long Xa1 = 40014L;
! 107: static long Xa2 = 40692L;
! 108:
! 109: /* Local function declarations, not visible outside this file */
! 110: static double confrac __PROTO((double a, double b, double x));
! 111: static double ibeta __PROTO((double a, double b, double x));
! 112: static double igamma __PROTO((double a, double x));
! 113: static double ranf __PROTO((double init));
! 114: static double inverse_normal_func __PROTO((double p));
! 115: static double inverse_error_func __PROTO((double p));
! 116:
! 117: #ifndef GAMMA
! 118: /* Provide GAMMA function for those who do not already have one */
! 119: static double lngamma __PROTO((double z));
! 120: static double lgamneg __PROTO((double z));
! 121: static double lgampos __PROTO((double z));
! 122:
! 123: /**
! 124: * from statlib, Thu Jan 23 15:02:27 EST 1992 ***
! 125: *
! 126: * This file contains two algorithms for the logarithm of the gamma function.
! 127: * Algorithm AS 245 is the faster (but longer) and gives an accuracy of about
! 128: * 10-12 significant decimal digits except for small regions around X = 1 and
! 129: * X = 2, where the function goes to zero.
! 130: * The second algorithm is not part of the AS algorithms. It is slower but
! 131: * gives 14 or more significant decimal digits accuracy, except around X = 1
! 132: * and X = 2. The Lanczos series from which this algorithm is derived is
! 133: * interesting in that it is a convergent series approximation for the gamma
! 134: * function, whereas the familiar series due to De Moivre (and usually wrongly
! 135: * called Stirling's approximation) is only an asymptotic approximation, as
! 136: * is the true and preferable approximation due to Stirling.
! 137: *
! 138: * Uses Lanczos-type approximation to ln(gamma) for z > 0. Reference: Lanczos,
! 139: * C. 'A precision approximation of the gamma function', J. SIAM Numer.
! 140: * Anal., B, 1, 86-96, 1964. Accuracy: About 14 significant digits except for
! 141: * small regions in the vicinity of 1 and 2.
! 142: *
! 143: * Programmer: Alan Miller CSIRO Division of Mathematics & Statistics
! 144: *
! 145: * Latest revision - 17 April 1988
! 146: *
! 147: * Additions: Translated from fortran to C, code added to handle values z < 0.
! 148: * The global variable signgam contains the sign of the gamma function.
! 149: *
! 150: * IMPORTANT: The signgam variable contains garbage until AFTER the call to
! 151: * lngamma().
! 152: *
! 153: * Permission granted to distribute freely for non-commercial purposes only
! 154: * Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl
! 155: */
! 156:
! 157: /* Local data, not visible outside this file
! 158: static double a[] =
! 159: {
! 160: 0.9999999999995183E+00,
! 161: 0.6765203681218835E+03,
! 162: -.1259139216722289E+04,
! 163: 0.7713234287757674E+03,
! 164: -.1766150291498386E+03,
! 165: 0.1250734324009056E+02,
! 166: -.1385710331296526E+00,
! 167: 0.9934937113930748E-05,
! 168: 0.1659470187408462E-06,
! 169: }; */
! 170:
! 171: /* from Ray Toy */
! 172: static double GPFAR a[] =
! 173: {
! 174: .99999999999980993227684700473478296744476168282198,
! 175: 676.52036812188509856700919044401903816411251975244084,
! 176: -1259.13921672240287047156078755282840836424300664868028,
! 177: 771.32342877765307884865282588943070775227268469602500,
! 178: -176.61502916214059906584551353999392943274507608117860,
! 179: 12.50734327868690481445893685327104972970563021816420,
! 180: -.13857109526572011689554706984971501358032683492780,
! 181: .00000998436957801957085956266828104544089848531228,
! 182: .00000015056327351493115583383579667028994545044040,
! 183: };
! 184:
! 185: static double lgamneg(z)
! 186: double z;
! 187: {
! 188: double tmp;
! 189:
! 190: /* Use reflection formula, then call lgampos() */
! 191: tmp = sin(z * PI);
! 192:
! 193: if (fabs(tmp) < MACHEPS) {
! 194: tmp = 0.0;
! 195: } else if (tmp < 0.0) {
! 196: tmp = -tmp;
! 197: signgam = -1;
! 198: }
! 199: return LNPI - lgampos(1.0 - z) - log(tmp);
! 200:
! 201: }
! 202:
! 203: static double lgampos(z)
! 204: double z;
! 205: {
! 206: double sum;
! 207: double tmp;
! 208: int i;
! 209:
! 210: sum = a[0];
! 211: for (i = 1, tmp = z; i < 9; i++) {
! 212: sum += a[i] / tmp;
! 213: tmp++;
! 214: }
! 215:
! 216: return log(sum) + LNSQRT2PI - z - 6.5 + (z - 0.5) * log(z + 6.5);
! 217: }
! 218:
! 219: static double lngamma(z)
! 220: double z;
! 221: {
! 222: signgam = 1;
! 223:
! 224: if (z <= 0.0)
! 225: return lgamneg(z);
! 226: else
! 227: return lgampos(z);
! 228: }
! 229:
! 230: # define GAMMA(x) lngamma ((x))
! 231: #endif /* !GAMMA */
! 232:
! 233: void f_erf()
! 234: {
! 235: struct value a;
! 236: double x;
! 237:
! 238: x = real(pop(&a));
! 239:
! 240: #ifdef HAVE_ERF
! 241: x = erf(x);
! 242: #else
! 243: {
! 244: int fsign;
! 245: fsign = x >= 0 ? 1 : 0;
! 246: x = igamma(0.5, (x)*(x));
! 247: if (x == -1.0) {
! 248: undefined = TRUE;
! 249: x = 0.0;
! 250: } else {
! 251: if (fsign == 0)
! 252: x = -x;
! 253: }
! 254: }
! 255: #endif
! 256: push(Gcomplex(&a, x, 0.0));
! 257: }
! 258:
! 259: void f_erfc()
! 260: {
! 261: struct value a;
! 262: double x;
! 263:
! 264: x = real(pop(&a));
! 265: #ifdef HAVE_ERFC
! 266: x = erfc(x);
! 267: #else
! 268: {
! 269: int fsign;
! 270: fsign = x >= 0 ? 1 : 0;
! 271: x = igamma(0.5, (x)*(x));
! 272: if (x == 1.0) {
! 273: undefined = TRUE;
! 274: x = 0.0;
! 275: } else {
! 276: x = fsign > 0 ? 1.0 - x : 1.0 + x ;
! 277: }
! 278: }
! 279: #endif
! 280: push(Gcomplex(&a, x, 0.0));
! 281: }
! 282:
! 283: void f_ibeta()
! 284: {
! 285: struct value a;
! 286: double x;
! 287: double arg1;
! 288: double arg2;
! 289:
! 290: x = real(pop(&a));
! 291: arg2 = real(pop(&a));
! 292: arg1 = real(pop(&a));
! 293:
! 294: x = ibeta(arg1, arg2, x);
! 295: if (x == -1.0) {
! 296: undefined = TRUE;
! 297: push(Ginteger(&a, 0));
! 298: } else
! 299: push(Gcomplex(&a, x, 0.0));
! 300: }
! 301:
! 302: void f_igamma()
! 303: {
! 304: struct value a;
! 305: double x;
! 306: double arg1;
! 307:
! 308: x = real(pop(&a));
! 309: arg1 = real(pop(&a));
! 310:
! 311: x = igamma(arg1, x);
! 312: if (x == -1.0) {
! 313: undefined = TRUE;
! 314: push(Ginteger(&a, 0));
! 315: } else
! 316: push(Gcomplex(&a, x, 0.0));
! 317: }
! 318:
! 319: void f_gamma()
! 320: {
! 321: register double y;
! 322: struct value a;
! 323:
! 324: y = GAMMA(real(pop(&a)));
! 325: if (y > 88.0) {
! 326: undefined = TRUE;
! 327: push(Ginteger(&a, 0));
! 328: } else
! 329: push(Gcomplex(&a, signgam * gp_exp(y), 0.0));
! 330: }
! 331:
! 332: void f_lgamma()
! 333: {
! 334: struct value a;
! 335:
! 336: push(Gcomplex(&a, GAMMA(real(pop(&a))), 0.0));
! 337: }
! 338:
! 339: #ifndef BADRAND
! 340:
! 341: void f_rand()
! 342: {
! 343: struct value a;
! 344:
! 345: push(Gcomplex(&a, ranf(real(pop(&a))), 0.0));
! 346: }
! 347:
! 348: #else /* BADRAND */
! 349:
! 350: /* Use only to observe the effect of a "bad" random number generator. */
! 351: void f_rand()
! 352: {
! 353: struct value a;
! 354:
! 355: static unsigned int y = 0;
! 356: unsigned int maxran = 1000;
! 357:
! 358: (void) real(pop(&a));
! 359: y = (781 * y + 387) % maxran;
! 360:
! 361: push(Gcomplex(&a, (double) y / maxran, 0.0));
! 362: }
! 363:
! 364: #endif /* BADRAND */
! 365:
! 366: /** ibeta.c
! 367: *
! 368: * DESCRIB Approximate the incomplete beta function Ix(a, b).
! 369: *
! 370: * _
! 371: * |(a + b) /x (a-1) (b-1)
! 372: * Ix(a, b) = -_-------_--- * | t * (1 - t) dt (a,b > 0)
! 373: * |(a) * |(b) /0
! 374: *
! 375: *
! 376: *
! 377: * CALL p = ibeta(a, b, x)
! 378: *
! 379: * double a > 0
! 380: * double b > 0
! 381: * double x [0, 1]
! 382: *
! 383: * WARNING none
! 384: *
! 385: * RETURN double p [0, 1]
! 386: * -1.0 on error condition
! 387: *
! 388: * XREF lngamma()
! 389: *
! 390: * BUGS none
! 391: *
! 392: * REFERENCE The continued fraction expansion as given by
! 393: * Abramowitz and Stegun (1964) is used.
! 394: *
! 395: * Permission granted to distribute freely for non-commercial purposes only
! 396: * Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl
! 397: */
! 398:
! 399: static double ibeta(a, b, x)
! 400: double a, b, x;
! 401: {
! 402: /* Test for admissibility of arguments */
! 403: if (a <= 0.0 || b <= 0.0)
! 404: return -1.0;
! 405: if (x < 0.0 || x > 1.0)
! 406: return -1.0;;
! 407:
! 408: /* If x equals 0 or 1, return x as prob */
! 409: if (x == 0.0 || x == 1.0)
! 410: return x;
! 411:
! 412: /* Swap a, b if necessarry for more efficient evaluation */
! 413: return a < x * (a + b) ? 1.0 - confrac(b, a, 1.0 - x) : confrac(a, b, x);
! 414: }
! 415:
! 416: static double confrac(a, b, x)
! 417: double a, b, x;
! 418: {
! 419: double Alo = 0.0;
! 420: double Ahi;
! 421: double Aev;
! 422: double Aod;
! 423: double Blo = 1.0;
! 424: double Bhi = 1.0;
! 425: double Bod = 1.0;
! 426: double Bev = 1.0;
! 427: double f;
! 428: double fold;
! 429: double Apb = a + b;
! 430: double d;
! 431: int i;
! 432: int j;
! 433:
! 434: /* Set up continued fraction expansion evaluation. */
! 435: Ahi = gp_exp(GAMMA(Apb) + a * log(x) + b * log(1.0 - x) -
! 436: GAMMA(a + 1.0) - GAMMA(b));
! 437:
! 438: /*
! 439: * Continued fraction loop begins here. Evaluation continues until
! 440: * maximum iterations are exceeded, or convergence achieved.
! 441: */
! 442: for (i = 0, j = 1, f = Ahi; i <= ITMAX; i++, j++) {
! 443: d = a + j + i;
! 444: Aev = -(a + i) * (Apb + i) * x / d / (d - 1.0);
! 445: Aod = j * (b - j) * x / d / (d + 1.0);
! 446: Alo = Bev * Ahi + Aev * Alo;
! 447: Blo = Bev * Bhi + Aev * Blo;
! 448: Ahi = Bod * Alo + Aod * Ahi;
! 449: Bhi = Bod * Blo + Aod * Bhi;
! 450:
! 451: if (fabs(Bhi) < MACHEPS)
! 452: Bhi = 0.0;
! 453:
! 454: if (Bhi != 0.0) {
! 455: fold = f;
! 456: f = Ahi / Bhi;
! 457: if (fabs(f - fold) < fabs(f) * MACHEPS)
! 458: return f;
! 459: }
! 460: }
! 461:
! 462: return -1.0;
! 463: }
! 464:
! 465: /** igamma.c
! 466: *
! 467: * DESCRIB Approximate the incomplete gamma function P(a, x).
! 468: *
! 469: * 1 /x -t (a-1)
! 470: * P(a, x) = -_--- * | e * t dt (a > 0)
! 471: * |(a) /0
! 472: *
! 473: * CALL p = igamma(a, x)
! 474: *
! 475: * double a > 0
! 476: * double x >= 0
! 477: *
! 478: * WARNING none
! 479: *
! 480: * RETURN double p [0, 1]
! 481: * -1.0 on error condition
! 482: *
! 483: * XREF lngamma()
! 484: *
! 485: * BUGS Values 0 <= x <= 1 may lead to inaccurate results.
! 486: *
! 487: * REFERENCE ALGORITHM AS239 APPL. STATIST. (1988) VOL. 37, NO. 3
! 488: *
! 489: * Permission granted to distribute freely for non-commercial purposes only
! 490: * Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl
! 491: */
! 492:
! 493: /* Global variables, not visible outside this file */
! 494: static double pn1, pn2, pn3, pn4, pn5, pn6;
! 495:
! 496: static double igamma(a, x)
! 497: double a, x;
! 498: {
! 499: double arg;
! 500: double aa;
! 501: double an;
! 502: double b;
! 503: int i;
! 504:
! 505: /* Check that we have valid values for a and x */
! 506: if (x < 0.0 || a <= 0.0)
! 507: return -1.0;
! 508:
! 509: /* Deal with special cases */
! 510: if (x == 0.0)
! 511: return 0.0;
! 512: if (x > XBIG)
! 513: return 1.0;
! 514:
! 515: /* Check value of factor arg */
! 516: arg = a * log(x) - x - GAMMA(a + 1.0);
! 517: if (arg < MINEXP)
! 518: return -1.0;
! 519: arg = gp_exp(arg);
! 520:
! 521: /* Choose infinite series or continued fraction. */
! 522:
! 523: if ((x > 1.0) && (x >= a + 2.0)) {
! 524: /* Use a continued fraction expansion */
! 525:
! 526: double rn;
! 527: double rnold;
! 528:
! 529: aa = 1.0 - a;
! 530: b = aa + x + 1.0;
! 531: pn1 = 1.0;
! 532: pn2 = x;
! 533: pn3 = x + 1.0;
! 534: pn4 = x * b;
! 535: rnold = pn3 / pn4;
! 536:
! 537: for (i = 1; i <= ITMAX; i++) {
! 538:
! 539: aa++;
! 540: b += 2.0;
! 541: an = aa * (double) i;
! 542:
! 543: pn5 = b * pn3 - an * pn1;
! 544: pn6 = b * pn4 - an * pn2;
! 545:
! 546: if (pn6 != 0.0) {
! 547:
! 548: rn = pn5 / pn6;
! 549: if (fabs(rnold - rn) <= GPMIN(MACHEPS, MACHEPS * rn))
! 550: return 1.0 - arg * rn * a;
! 551:
! 552: rnold = rn;
! 553: }
! 554: pn1 = pn3;
! 555: pn2 = pn4;
! 556: pn3 = pn5;
! 557: pn4 = pn6;
! 558:
! 559: /* Re-scale terms in continued fraction if terms are large */
! 560: if (fabs(pn5) >= OFLOW) {
! 561:
! 562: pn1 /= OFLOW;
! 563: pn2 /= OFLOW;
! 564: pn3 /= OFLOW;
! 565: pn4 /= OFLOW;
! 566: }
! 567: }
! 568: } else {
! 569: /* Use Pearson's series expansion. */
! 570:
! 571: for (i = 0, aa = a, an = b = 1.0; i <= ITMAX; i++) {
! 572:
! 573: aa++;
! 574: an *= x / aa;
! 575: b += an;
! 576: if (an < b * MACHEPS)
! 577: return arg * b;
! 578: }
! 579: }
! 580: return -1.0;
! 581: }
! 582:
! 583: /***********************************************************************
! 584: double ranf(double init)
! 585: RANDom number generator as a Function
! 586: Returns a random floating point number from a uniform distribution
! 587: over 0 - 1 (endpoints of this interval are not returned) using a
! 588: large integer generator.
! 589: This is a transcription from Pascal to Fortran of routine
! 590: Uniform_01 from the paper
! 591: L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
! 592: with Splitting Facilities." ACM Transactions on Mathematical
! 593: Software, 17:98-111 (1991)
! 594:
! 595: GeNerate LarGe Integer
! 596: Returns a random integer following a uniform distribution over
! 597: (1, 2147483562) using the generator.
! 598: This is a transcription from Pascal to Fortran of routine
! 599: Random from the paper
! 600: L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
! 601: with Splitting Facilities." ACM Transactions on Mathematical
! 602: Software, 17:98-111 (1991)
! 603: ***********************************************************************/
! 604: static double ranf(init)
! 605: double init;
! 606: {
! 607:
! 608: long k, z;
! 609: static int firsttime = 1;
! 610: static long s1, s2;
! 611:
! 612: /* (Re)-Initialize seeds if necessary */
! 613: if (init < 0.0 || firsttime == 1) {
! 614: firsttime = 0;
! 615: s1 = 1234567890L;
! 616: s2 = 1234567890L;
! 617: }
! 618: /* Generate pseudo random integers */
! 619: k = s1 / 53668L;
! 620: s1 = Xa1 * (s1 - k * 53668L) - k * 12211;
! 621: if (s1 < 0)
! 622: s1 += Xm1;
! 623: k = s2 / 52774L;
! 624: s2 = Xa2 * (s2 - k * 52774L) - k * 3791;
! 625: if (s2 < 0)
! 626: s2 += Xm2;
! 627: z = s1 - s2;
! 628: if (z < 1)
! 629: z += (Xm1 - 1);
! 630:
! 631: /*
! 632: * 4.656613057E-10 is 1/Xm1. Xm1 is set at the top of this file and is
! 633: * currently 2147483563. If Xm1 changes, change this also.
! 634: */
! 635: return (double) 4.656613057E-10 *z;
! 636: }
! 637:
! 638: /* ----------------------------------------------------------------
! 639: Following to specfun.c made by John Grosh (jgrosh@arl.mil)
! 640: on 28 OCT 1992.
! 641: ---------------------------------------------------------------- */
! 642:
! 643: void f_normal()
! 644: { /* Normal or Gaussian Probability Function */
! 645: struct value a;
! 646: double x;
! 647:
! 648: /* ref. Abramowitz and Stegun 1964, "Handbook of Mathematical
! 649: Functions", Applied Mathematics Series, vol 55,
! 650: Chapter 26, page 934, Eqn. 26.2.29 and Jos van der Woude
! 651: code found above */
! 652:
! 653: x = real(pop(&a));
! 654:
! 655: x = 0.5 * SQRT_TWO * x;
! 656: #ifdef HAVE_ERF
! 657: x = 0.5 * (1.0 + erf(x));
! 658: #else
! 659: {
! 660: int fsign;
! 661: fsign = x >= 0 ? 1 : 0;
! 662: x = igamma(0.5, (x)*(x));
! 663: if (x == 1.0) {
! 664: undefined = TRUE;
! 665: x = 0.0;
! 666: } else {
! 667: if (fsign == 0)
! 668: x = -(x);
! 669: x = 0.5 * (1.0 + x);
! 670: }
! 671: }
! 672: #endif
! 673: push(Gcomplex(&a, x, 0.0));
! 674: }
! 675:
! 676: void f_inverse_normal()
! 677: { /* Inverse normal distribution function */
! 678: struct value a;
! 679: double x;
! 680:
! 681: x = real(pop(&a));
! 682:
! 683: if (x <= 0.0 || x >= 1.0) {
! 684: undefined = TRUE;
! 685: push(Gcomplex(&a, 0.0, 0.0));
! 686: } else {
! 687: push(Gcomplex(&a, inverse_normal_func(x), 0.0));
! 688: }
! 689: }
! 690:
! 691:
! 692: void f_inverse_erf()
! 693: { /* Inverse error function */
! 694: struct value a;
! 695: double x;
! 696:
! 697: x = real(pop(&a));
! 698:
! 699: if (fabs(x) >= 1.0) {
! 700: undefined = TRUE;
! 701: push(Gcomplex(&a, 0.0, 0.0));
! 702: } else {
! 703: push(Gcomplex(&a, inverse_error_func(x), 0.0));
! 704: }
! 705: }
! 706:
! 707: static double inverse_normal_func(p)
! 708: double p;
! 709: {
! 710: /*
! 711: Source: This routine was derived (using f2c) from the
! 712: FORTRAN subroutine MDNRIS found in
! 713: ACM Algorithm 602 obtained from netlib.
! 714:
! 715: MDNRIS code contains the 1978 Copyright
! 716: by IMSL, INC. . Since MDNRIS has been
! 717: submitted to netlib it may be used with
! 718: the restriction that it may only be
! 719: used for noncommercial purposes and that
! 720: IMSL be acknowledged as the copyright-holder
! 721: of the code.
! 722: */
! 723:
! 724: /* Initialized data */
! 725: static double eps = 1e-10;
! 726: static double g0 = 1.851159e-4;
! 727: static double g1 = -.002028152;
! 728: static double g2 = -.1498384;
! 729: static double g3 = .01078639;
! 730: static double h0 = .09952975;
! 731: static double h1 = .5211733;
! 732: static double h2 = -.06888301;
! 733: static double sqrt2 = 1.414213562373095;
! 734:
! 735: /* Local variables */
! 736: static double a, w, x;
! 737: static double sd, wi, sn, y;
! 738:
! 739: /* Note: 0.0 < p < 1.0 */
! 740:
! 741: /* p too small, compute y directly */
! 742: if (p <= eps) {
! 743: a = p + p;
! 744: w = sqrt(-(double) log(a + (a - a * a)));
! 745:
! 746: /* use a rational function in 1.0 / w */
! 747: wi = 1.0 / w;
! 748: sn = ((g3 * wi + g2) * wi + g1) * wi;
! 749: sd = ((wi + h2) * wi + h1) * wi + h0;
! 750: y = w + w * (g0 + sn / sd);
! 751: y = -y * sqrt2;
! 752: } else {
! 753: x = 1.0 - (p + p);
! 754: y = inverse_error_func(x);
! 755: y = -sqrt2 * y;
! 756: }
! 757: return (y);
! 758: }
! 759:
! 760:
! 761: static double inverse_error_func(p)
! 762: double p;
! 763: {
! 764: /*
! 765: Source: This routine was derived (using f2c) from the
! 766: FORTRAN subroutine MERFI found in
! 767: ACM Algorithm 602 obtained from netlib.
! 768:
! 769: MDNRIS code contains the 1978 Copyright
! 770: by IMSL, INC. . Since MERFI has been
! 771: submitted to netlib, it may be used with
! 772: the restriction that it may only be
! 773: used for noncommercial purposes and that
! 774: IMSL be acknowledged as the copyright-holder
! 775: of the code.
! 776: */
! 777:
! 778:
! 779:
! 780: /* Initialized data */
! 781: static double a1 = -.5751703;
! 782: static double a2 = -1.896513;
! 783: static double a3 = -.05496261;
! 784: static double b0 = -.113773;
! 785: static double b1 = -3.293474;
! 786: static double b2 = -2.374996;
! 787: static double b3 = -1.187515;
! 788: static double c0 = -.1146666;
! 789: static double c1 = -.1314774;
! 790: static double c2 = -.2368201;
! 791: static double c3 = .05073975;
! 792: static double d0 = -44.27977;
! 793: static double d1 = 21.98546;
! 794: static double d2 = -7.586103;
! 795: static double e0 = -.05668422;
! 796: static double e1 = .3937021;
! 797: static double e2 = -.3166501;
! 798: static double e3 = .06208963;
! 799: static double f0 = -6.266786;
! 800: static double f1 = 4.666263;
! 801: static double f2 = -2.962883;
! 802: static double g0 = 1.851159e-4;
! 803: static double g1 = -.002028152;
! 804: static double g2 = -.1498384;
! 805: static double g3 = .01078639;
! 806: static double h0 = .09952975;
! 807: static double h1 = .5211733;
! 808: static double h2 = -.06888301;
! 809:
! 810: /* Local variables */
! 811: static double a, b, f, w, x, y, z, sigma, z2, sd, wi, sn;
! 812:
! 813: x = p;
! 814:
! 815: /* determine sign of x */
! 816: if (x > 0)
! 817: sigma = 1.0;
! 818: else
! 819: sigma = -1.0;
! 820:
! 821: /* Note: -1.0 < x < 1.0 */
! 822:
! 823: z = fabs(x);
! 824:
! 825: /* z between 0.0 and 0.85, approx. f by a
! 826: rational function in z */
! 827:
! 828: if (z <= 0.85) {
! 829: z2 = z * z;
! 830: f = z + z * (b0 + a1 * z2 / (b1 + z2 + a2
! 831: / (b2 + z2 + a3 / (b3 + z2))));
! 832:
! 833: /* z greater than 0.85 */
! 834: } else {
! 835: a = 1.0 - z;
! 836: b = z;
! 837:
! 838: /* reduced argument is in (0.85,1.0),
! 839: obtain the transformed variable */
! 840:
! 841: w = sqrt(-(double) log(a + a * b));
! 842:
! 843: /* w greater than 4.0, approx. f by a
! 844: rational function in 1.0 / w */
! 845:
! 846: if (w >= 4.0) {
! 847: wi = 1.0 / w;
! 848: sn = ((g3 * wi + g2) * wi + g1) * wi;
! 849: sd = ((wi + h2) * wi + h1) * wi + h0;
! 850: f = w + w * (g0 + sn / sd);
! 851:
! 852: /* w between 2.5 and 4.0, approx.
! 853: f by a rational function in w */
! 854:
! 855: } else if (w < 4.0 && w > 2.5) {
! 856: sn = ((e3 * w + e2) * w + e1) * w;
! 857: sd = ((w + f2) * w + f1) * w + f0;
! 858: f = w + w * (e0 + sn / sd);
! 859:
! 860: /* w between 1.13222 and 2.5, approx. f by
! 861: a rational function in w */
! 862: } else if (w <= 2.5 && w > 1.13222) {
! 863: sn = ((c3 * w + c2) * w + c1) * w;
! 864: sd = ((w + d2) * w + d1) * w + d0;
! 865: f = w + w * (c0 + sn / sd);
! 866: }
! 867: }
! 868: y = sigma * f;
! 869: return (y);
! 870: }
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