Annotation of OpenXM_contrib/gnuplot/specfun.c, Revision 1.1.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