1 /* GNUPLOT - specfun.c */
2
3 /*[
4 * Copyright 1986 - 1993, 1998, 2004 Thomas Williams, Colin Kelley
5 *
6 * Permission to use, copy, and distribute this software and its
7 * documentation for any purpose with or without fee is hereby granted,
8 * provided that the above copyright notice appear in all copies and
9 * that both that copyright notice and this permission notice appear
10 * in supporting documentation.
11 *
12 * Permission to modify the software is granted, but not the right to
13 * distribute the complete modified source code. Modifications are to
14 * be distributed as patches to the released version. Permission to
15 * distribute binaries produced by compiling modified sources is granted,
16 * provided you
17 * 1. distribute the corresponding source modifications from the
18 * released version in the form of a patch file along with the binaries,
19 * 2. add special version identification to distinguish your version
20 * in addition to the base release version number,
21 * 3. provide your name and address as the primary contact for the
22 * support of your modified version, and
23 * 4. retain our contact information in regard to use of the base
24 * software.
25 * Permission to distribute the released version of the source code along
26 * with corresponding source modifications in the form of a patch file is
27 * granted with same provisions 2 through 4 for binary distributions.
28 *
29 * This software is provided "as is" without express or implied warranty
30 * to the extent permitted by applicable law.
31 ]*/
32
33 /*
34 * AUTHORS
35 *
36 * Original Software:
37 * Jos van der Woude, jvdwoude@hut.nl
38 *
39 */
40
41 /* FIXME:
42 * plain comparisons of floating point numbers!
43 */
44
45 #include "specfun.h"
46 #include "stdfn.h"
47 #include "util.h"
48
49 #define ITMAX 200
50
51 #ifdef FLT_EPSILON
52 # define MACHEPS FLT_EPSILON /* 1.0E-08 */
53 #else
54 # define MACHEPS 1.0E-08
55 #endif
56
57 #ifndef E_MINEXP
58 /* AS239 value, e^-88 = 2^-127 */
59 #define E_MINEXP (-88.0)
60 #endif
61 #ifndef E_MAXEXP
62 #define E_MAXEXP (-E_MINEXP)
63 #endif
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 #if defined(GAMMA) && !HAVE_DECL_SIGNGAM
98 extern int signgam; /* this is not always declared in math.h */
99 #endif
100
101 /* Local function declarations, not visible outside this file */
102 static int mtherr __PROTO((char *, int));
103 static double polevl __PROTO((double x, const double coef[], int N));
104 static double p1evl __PROTO((double x, const double coef[], int N));
105 static double confrac __PROTO((double a, double b, double x));
106 static double ibeta __PROTO((double a, double b, double x));
107 static double igamma __PROTO((double a, double x));
108 static double ranf __PROTO((struct value * init));
109 static double inverse_error_func __PROTO((double p));
110 static double inverse_normal_func __PROTO((double p));
111 static double lambertw __PROTO((double x));
112 #if (0) /* Only used by low-precision Airy version */
113 static double airy_neg __PROTO(( double x ));
114 static double airy_pos __PROTO((double x));
115 #endif
116 #ifndef HAVE_LIBCERF
117 static double humlik __PROTO((double x, double y));
118 #endif
119 static double expint __PROTO((double n, double x));
120 #ifndef GAMMA
121 static int ISNAN __PROTO((double x));
122 static int ISFINITE __PROTO((double x));
123 static double lngamma __PROTO((double z));
124 #endif
125 #ifndef HAVE_ERF
126 static double erf __PROTO((double a));
127 #endif
128 #ifndef HAVE_ERFC
129 static double erfc __PROTO((double a));
130 #endif
131
132 /* Macros to configure routines taken from CEPHES: */
133
134 /* Unknown arithmetic, invokes coefficients given in
135 * normal decimal format. Beware of range boundary
136 * problems (MACHEP, MAXLOG, etc. in const.c) and
137 * roundoff problems in pow.c:
138 * (Sun SPARCstation)
139 */
140 #define UNK 1
141 #define MACHEP DBL_EPSILON
142 #define MAXNUM DBL_MAX
143
144 /* Define to support tiny denormal numbers, else undefine. */
145 #define DENORMAL 1
146
147 /* Define to ask for infinity support, else undefine. */
148 #define INFINITIES 1
149
150 /* Define to ask for support of numbers that are Not-a-Number,
151 else undefine. This may automatically define INFINITIES in some files. */
152 #define NANS 1
153
154 /* Define to distinguish between -0.0 and +0.0. */
155 #define MINUSZERO 1
156
157 /*
158 Cephes Math Library Release 2.0: April, 1987
159 Copyright 1984, 1987 by Stephen L. Moshier
160 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
161 */
162
163 static int merror = 0;
164
165 /* Notice: the order of appearance of the following messages cannot be bound
166 * to error codes defined in mconf.h or math.h or similar, as these files are
167 * not available on every platform. Thus, enumerate them explicitly.
168 */
169 #define MTHERR_DOMAIN 1
170 #define MTHERR_SING 2
171 #define MTHERR_OVERFLOW 3
172 #define MTHERR_UNDERFLOW 4
173 #define MTHERR_TLPREC 5
174 #define MTHERR_PLPREC 6
175
176 static int
mtherr(char * name,int code)177 mtherr(char *name, int code)
178 {
179 static const char *ermsg[7] = {
180 "unknown", /* error code 0 */
181 "domain", /* error code 1 */
182 "singularity", /* et seq. */
183 "overflow",
184 "underflow",
185 "total loss of precision",
186 "partial loss of precision"
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 by the code argument. */
199 if ((code <= 0) || (code >= 7))
200 code = 0;
201 printf("%s error\n", ermsg[code]);
202
203 /* Return to calling program */
204 return (0);
205 }
206
207 /* polevl.c
208 * p1evl.c
209 *
210 * Evaluate polynomial
211 *
212 *
213 *
214 * SYNOPSIS:
215 *
216 * int N;
217 * double x, y, coef[N+1], polevl[];
218 *
219 * y = polevl( x, coef, N );
220 *
221 *
222 *
223 * DESCRIPTION:
224 *
225 * Evaluates polynomial of degree N:
226 *
227 * 2 N
228 * y = C + C x + C x +...+ C x
229 * 0 1 2 N
230 *
231 * Coefficients are stored in reverse order:
232 *
233 * coef[0] = C , ..., coef[N] = C .
234 * N 0
235 *
236 * The function p1evl() assumes that coef[N] = 1.0 and is
237 * omitted from the array. Its calling arguments are
238 * otherwise the same as polevl().
239 *
240 *
241 * SPEED:
242 *
243 * In the interest of speed, there are no checks for out
244 * of bounds arithmetic. This routine is used by most of
245 * the functions in the library. Depending on available
246 * equipment features, the user may wish to rewrite the
247 * program in microcode or assembly language.
248 *
249 */
250
251 /*
252 Cephes Math Library Release 2.1: December, 1988
253 Copyright 1984, 1987, 1988 by Stephen L. Moshier
254 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
255 */
256 static double
polevl(double x,const double coef[],int N)257 polevl(double x, const double coef[], int N)
258 {
259 double ans;
260 int i;
261 const double *p;
262
263 p = coef;
264 ans = *p++;
265 i = N;
266
267 do
268 ans = ans * x + *p++;
269 while (--i);
270
271 return (ans);
272 }
273
274 /* N
275 * Evaluate polynomial when coefficient of x is 1.0.
276 * Otherwise same as polevl.
277 */
278 static double
p1evl(double x,const double coef[],int N)279 p1evl(double x, const double coef[], int N)
280 {
281 double ans;
282 const double *p;
283 int i;
284
285 p = coef;
286 ans = x + *p++;
287 i = N - 1;
288
289 do
290 ans = ans * x + *p++;
291 while (--i);
292
293 return (ans);
294 }
295
296 #ifndef GAMMA
297
298 /* Provide GAMMA function for those who do not already have one */
299
300 int sgngam;
301
302 static int
ISNAN(double x)303 ISNAN(double x)
304 {
305 volatile double a = x;
306
307 if (a != a)
308 return 1;
309 return 0;
310 }
311
312 static int
ISFINITE(double x)313 ISFINITE(double x)
314 {
315 volatile double a = x;
316
317 if (a < DBL_MAX)
318 return 1;
319 return 0;
320 }
321
322 double
lngamma(double x)323 lngamma(double x)
324 {
325 /* A[]: Stirling's formula expansion of log gamma
326 * B[], C[]: log gamma function between 2 and 3
327 */
328 #ifdef UNK
329 static const double A[] = {
330 8.11614167470508450300E-4,
331 -5.95061904284301438324E-4,
332 7.93650340457716943945E-4,
333 -2.77777777730099687205E-3,
334 8.33333333333331927722E-2
335 };
336 static const double B[] = {
337 -1.37825152569120859100E3,
338 -3.88016315134637840924E4,
339 -3.31612992738871184744E5,
340 -1.16237097492762307383E6,
341 -1.72173700820839662146E6,
342 -8.53555664245765465627E5
343 };
344 static const double C[] = {
345 /* 1.00000000000000000000E0, */
346 -3.51815701436523470549E2,
347 -1.70642106651881159223E4,
348 -2.20528590553854454839E5,
349 -1.13933444367982507207E6,
350 -2.53252307177582951285E6,
351 -2.01889141433532773231E6
352 };
353 /* log( sqrt( 2*pi ) ) */
354 static const double LS2PI = 0.91893853320467274178;
355 #define MAXLGM 2.556348e305
356 #endif /* UNK */
357
358 #ifdef DEC
359 static const unsigned short A[] = {
360 0035524, 0141201, 0034633, 0031405,
361 0135433, 0176755, 0126007, 0045030,
362 0035520, 0006371, 0003342, 0172730,
363 0136066, 0005540, 0132605, 0026407,
364 0037252, 0125252, 0125252, 0125132
365 };
366 static const unsigned short B[] = {
367 0142654, 0044014, 0077633, 0035410,
368 0144027, 0110641, 0125335, 0144760,
369 0144641, 0165637, 0142204, 0047447,
370 0145215, 0162027, 0146246, 0155211,
371 0145322, 0026110, 0010317, 0110130,
372 0145120, 0061472, 0120300, 0025363
373 };
374 static const unsigned short C[] = {
375 /*0040200,0000000,0000000,0000000*/
376 0142257, 0164150, 0163630, 0112622,
377 0143605, 0050153, 0156116, 0135272,
378 0144527, 0056045, 0145642, 0062332,
379 0145213, 0012063, 0106250, 0001025,
380 0145432, 0111254, 0044577, 0115142,
381 0145366, 0071133, 0050217, 0005122
382 };
383 /* log( sqrt( 2*pi ) ) */
384 static const unsigned short LS2P[] = {040153, 037616, 041445, 0172645,};
385 #define LS2PI *(double *)LS2P
386 #define MAXLGM 2.035093e36
387 #endif /* DEC */
388
389 #ifdef IBMPC
390 static const unsigned short A[] = {
391 0x6661, 0x2733, 0x9850, 0x3f4a,
392 0xe943, 0xb580, 0x7fbd, 0xbf43,
393 0x5ebb, 0x20dc, 0x019f, 0x3f4a,
394 0xa5a1, 0x16b0, 0xc16c, 0xbf66,
395 0x554b, 0x5555, 0x5555, 0x3fb5
396 };
397 static const unsigned short B[] = {
398 0x6761, 0x8ff3, 0x8901, 0xc095,
399 0xb93e, 0x355b, 0xf234, 0xc0e2,
400 0x89e5, 0xf890, 0x3d73, 0xc114,
401 0xdb51, 0xf994, 0xbc82, 0xc131,
402 0xf20b, 0x0219, 0x4589, 0xc13a,
403 0x055e, 0x5418, 0x0c67, 0xc12a
404 };
405 static const unsigned short C[] = {
406 /*0x0000,0x0000,0x0000,0x3ff0,*/
407 0x12b2, 0x1cf3, 0xfd0d, 0xc075,
408 0xd757, 0x7b89, 0xaa0d, 0xc0d0,
409 0x4c9b, 0xb974, 0xeb84, 0xc10a,
410 0x0043, 0x7195, 0x6286, 0xc131,
411 0xf34c, 0x892f, 0x5255, 0xc143,
412 0xe14a, 0x6a11, 0xce4b, 0xc13e
413 };
414 /* log( sqrt( 2*pi ) ) */
415 static const unsigned short LS2P[] = {
416 0xbeb5, 0xc864, 0x67f1, 0x3fed
417 };
418 #define LS2PI *(double *)LS2P
419 #define MAXLGM 2.556348e305
420 #endif /* IBMPC */
421
422 #ifdef MIEEE
423 static const unsigned short A[] = {
424 0x3f4a, 0x9850, 0x2733, 0x6661,
425 0xbf43, 0x7fbd, 0xb580, 0xe943,
426 0x3f4a, 0x019f, 0x20dc, 0x5ebb,
427 0xbf66, 0xc16c, 0x16b0, 0xa5a1,
428 0x3fb5, 0x5555, 0x5555, 0x554b
429 };
430 static const unsigned short B[] = {
431 0xc095, 0x8901, 0x8ff3, 0x6761,
432 0xc0e2, 0xf234, 0x355b, 0xb93e,
433 0xc114, 0x3d73, 0xf890, 0x89e5,
434 0xc131, 0xbc82, 0xf994, 0xdb51,
435 0xc13a, 0x4589, 0x0219, 0xf20b,
436 0xc12a, 0x0c67, 0x5418, 0x055e
437 };
438 static const unsigned short C[] = {
439 0xc075, 0xfd0d, 0x1cf3, 0x12b2,
440 0xc0d0, 0xaa0d, 0x7b89, 0xd757,
441 0xc10a, 0xeb84, 0xb974, 0x4c9b,
442 0xc131, 0x6286, 0x7195, 0x0043,
443 0xc143, 0x5255, 0x892f, 0xf34c,
444 0xc13e, 0xce4b, 0x6a11, 0xe14a
445 };
446 /* log( sqrt( 2*pi ) ) */
447 static const unsigned short LS2P[] = {
448 0x3fed, 0x67f1, 0xc864, 0xbeb5
449 };
450 #define LS2PI *(double *)LS2P
451 #define MAXLGM 2.556348e305
452 #endif /* MIEEE */
453
454 static const double LOGPI = 1.1447298858494001741434273513530587116472948129153;
455
456 double p, q, u, w, z;
457 int i;
458
459 sgngam = 1;
460 #ifdef NANS
461 if (ISNAN(x))
462 return (x);
463 #endif
464
465 #ifdef INFINITIES
466 if (!ISFINITE((x)))
467 return (DBL_MAX * DBL_MAX);
468 #endif
469
470 if (x < -34.0) {
471 q = -x;
472 w = lngamma(q); /* note this modifies sgngam! */
473 p = floor(q);
474 if (p == q) {
475 lgsing:
476 #ifdef INFINITIES
477 mtherr("lngamma", MTHERR_SING);
478 return (DBL_MAX * DBL_MAX);
479 #else
480 goto loverf;
481 #endif
482 }
483 i = p;
484 if ((i & 1) == 0)
485 sgngam = -1;
486 else
487 sgngam = 1;
488 z = q - p;
489 if (z > 0.5) {
490 p += 1.0;
491 z = p - q;
492 }
493 z = q * sin(PI * z);
494 if (z == 0.0)
495 goto lgsing;
496 /* z = log(PI) - log( z ) - w;*/
497 z = LOGPI - log(z) - w;
498 return (z);
499 }
500 if (x < 13.0) {
501 z = 1.0;
502 p = 0.0;
503 u = x;
504 while (u >= 3.0) {
505 p -= 1.0;
506 u = x + p;
507 z *= u;
508 }
509 while (u < 2.0) {
510 if (u == 0.0)
511 goto lgsing;
512 z /= u;
513 p += 1.0;
514 u = x + p;
515 }
516 if (z < 0.0) {
517 sgngam = -1;
518 z = -z;
519 } else
520 sgngam = 1;
521 if (u == 2.0)
522 return (log(z));
523 p -= 2.0;
524 x = x + p;
525 p = x * polevl(x, B, 5) / p1evl(x, C, 6);
526 return (log(z) + p);
527 }
528 if (x > MAXLGM) {
529 #ifdef INFINITIES
530 return (sgngam * (DBL_MAX * DBL_MAX));
531 #else
532 loverf:
533 mtherr("lngamma", MTHERR_OVERFLOW);
534 return (sgngam * MAXNUM);
535 #endif
536 }
537 q = (x - 0.5) * log(x) - x + LS2PI;
538 if (x > 1.0e8)
539 return (q);
540
541 p = 1.0 / (x * x);
542 if (x >= 1000.0)
543 q += ((7.9365079365079365079365e-4 * p
544 - 2.7777777777777777777778e-3) * p
545 + 0.0833333333333333333333) / x;
546 else
547 q += polevl(p, A, 4) / x;
548 return (q);
549 }
550
551 #define GAMMA(x) lngamma ((x))
552 /* HBB 20030816: must override name of sgngam so f_gamma() uses it */
553 #define signgam sgngam
554
555 #endif /* !GAMMA */
556
557 /*
558 * Make all the following internal routines f_whatever() perform
559 * autoconversion from string to numeric value.
560 */
561 #define pop(x) pop_or_convert_from_string(x)
562
563 void
f_erf(union argument * arg)564 f_erf(union argument *arg)
565 {
566 struct value a;
567 double x;
568
569 (void) arg; /* avoid -Wunused warning */
570 x = real(pop(&a));
571 x = erf(x);
572 push(Gcomplex(&a, x, 0.0));
573 }
574
575 void
f_erfc(union argument * arg)576 f_erfc(union argument *arg)
577 {
578 struct value a;
579 double x;
580
581 (void) arg; /* avoid -Wunused warning */
582 x = real(pop(&a));
583 x = erfc(x);
584 push(Gcomplex(&a, x, 0.0));
585 }
586
587 void
f_ibeta(union argument * arg)588 f_ibeta(union argument *arg)
589 {
590 struct value a;
591 double x;
592 double arg1;
593 double arg2;
594
595 (void) arg; /* avoid -Wunused warning */
596 x = real(pop(&a));
597 arg2 = real(pop(&a));
598 arg1 = real(pop(&a));
599
600 x = ibeta(arg1, arg2, x);
601 if (x == -1.0) {
602 undefined = TRUE;
603 push(Ginteger(&a, 0));
604 } else
605 push(Gcomplex(&a, x, 0.0));
606 }
607
f_igamma(union argument * arg)608 void f_igamma(union argument *arg)
609 {
610 struct value a;
611 double x;
612 double arg1;
613
614 (void) arg; /* avoid -Wunused warning */
615 x = real(pop(&a));
616 arg1 = real(pop(&a));
617
618 x = igamma(arg1, x);
619 if (x == -1.0) {
620 undefined = TRUE;
621 push(Ginteger(&a, 0));
622 } else
623 push(Gcomplex(&a, x, 0.0));
624 }
625
f_gamma(union argument * arg)626 void f_gamma(union argument *arg)
627 {
628 double y;
629 struct value a;
630
631 (void) arg; /* avoid -Wunused warning */
632 y = GAMMA(real(pop(&a)));
633 if (y > E_MAXEXP) {
634 undefined = TRUE;
635 push(Ginteger(&a, 0));
636 } else
637 push(Gcomplex(&a, signgam * gp_exp(y), 0.0));
638 }
639
f_lgamma(union argument * arg)640 void f_lgamma(union argument *arg)
641 {
642 struct value a;
643
644 (void) arg; /* avoid -Wunused warning */
645 push(Gcomplex(&a, GAMMA(real(pop(&a))), 0.0));
646 }
647
648 #ifndef BADRAND
649
f_rand(union argument * arg)650 void f_rand(union argument *arg)
651 {
652 struct value a;
653
654 (void) arg; /* avoid -Wunused warning */
655 push(Gcomplex(&a, ranf(pop(&a)), 0.0));
656 }
657
658 #else /* BADRAND */
659
660 /* Use only to observe the effect of a "bad" random number generator. */
f_rand(union argument * arg)661 void f_rand(union argument *arg)
662 {
663 struct value a;
664
665 (void) arg; /* avoid -Wunused warning */
666 static unsigned int y = 0;
667 unsigned int maxran = 1000;
668
669 (void) real(pop(&a));
670 y = (781 * y + 387) % maxran;
671
672 push(Gcomplex(&a, (double) y / maxran, 0.0));
673 }
674
675 #endif /* BADRAND */
676
677 /*
678 * Fallback implementation of the Faddeeva/Voigt function
679 * w(z) = exp(*-z^2) * erfc(-i*z)
680 * if not available from libcerf or some other library
681 */
682 #ifndef HAVE_LIBCERF
683 void
f_voigt(union argument * arg)684 f_voigt(union argument *arg)
685 {
686 struct value a;
687 double x,y;
688 (void) arg; /* avoid -Wunused warning */
689 y = real(pop(&a));
690 x = real(pop(&a));
691 push(Gcomplex(&a, humlik(x, y), 0.0));
692 }
693
694 /*
695 * Calculate the Voigt/Faddeeva function with relative error less than 10^(-4).
696 * (see http://www.atm.ox.ac.uk/user/wells/voigt.html)
697 *
698 * K(x,y) = \frac{y}{\pi} \int{\frac{e^{-t^2}}{(x-t)^2+y^2}}dt
699 *
700 * arguments:
701 * x, y - real and imaginary components of complex argument
702 * return value
703 * real value K(x,y)
704 *
705 * Algorithm: Josef Humlíček JQSRT 27 (1982) pp 437
706 * Fortran program by J.R. Wells JQSRT 62 (1999) pp 29-48.
707 * Translated to C++ with f2c program and modified by Marcin Wojdyr
708 * Minor adaptations from C++ to C by E. Stambulchik
709 * Adapted for gnuplot by Tommaso Vinci
710 */
humlik(double x,double y)711 static double humlik(double x, double y)
712 {
713
714 const double c[6] = { 1.0117281, -0.75197147, 0.012557727,
715 0.010022008, -2.4206814e-4, 5.0084806e-7 };
716 const double s[6] = { 1.393237, 0.23115241, -0.15535147,
717 0.0062183662, 9.1908299e-5, -6.2752596e-7 };
718 const double t[6] = { 0.31424038, 0.94778839, 1.5976826,
719 2.2795071, 3.020637, 3.8897249 };
720
721 const double rrtpi = 0.56418958; /* 1/SQRT(pi) */
722
723 double a0, d0, d2, e0, e2, e4, h0, h2, h4, h6,
724 p0, p2, p4, p6, p8, z0, z2, z4, z6, z8;
725 double mf[6], pf[6], mq[6], pq[6], xm[6], ym[6], xp[6], yp[6];
726 bool rg1, rg2, rg3;
727 double xlim0, xlim1, xlim2, xlim3, xlim4;
728 double yq, yrrtpi;
729 double abx, xq;
730 double k;
731
732 yq = y * y;
733 yrrtpi = y * rrtpi;
734 rg1 = true, rg2 = true, rg3 = true;
735 abx = fabs(x);
736 xq = abx * abx;
737
738 if (y >= 70.55)
739 return yrrtpi / (xq + yq);
740
741 xlim0 = sqrt(y * (40. - y * 3.6) + 15100.);
742 xlim1 = (y >= 8.425 ? 0. : sqrt(164. - y * (y * 1.8 + 4.3)));
743 xlim2 = 6.8 - y;
744 xlim3 = y * 2.4;
745 xlim4 = y * 18.1 + 1.65;
746 if (y <= 1e-6)
747 xlim2 = xlim1 = xlim0;
748
749 if (abx >= xlim0) /* Region 0 algorithm */
750 return yrrtpi / (xq + yq);
751
752 else if (abx >= xlim1) { /* Humlicek W4 Region 1 */
753 if (rg1) { /* First point in Region 1 */
754 rg1 = false;
755 a0 = yq + 0.5; /* Region 1 y-dependents */
756 d0 = a0 * a0;
757 d2 = yq + yq - 1.;
758 }
759 return rrtpi / (d0 + xq * (d2 + xq)) * y * (a0 + xq);
760 }
761
762 else if (abx > xlim2) { /* Humlicek W4 Region 2 */
763 if (rg2) { /* First point in Region 2 */
764 rg2 = false;
765 /* Region 2 y-dependents */
766 h0 = yq * (yq * (yq * (yq + 6.) + 10.5) + 4.5) + 0.5625;
767 h2 = yq * (yq * (yq * 4. + 6.) + 9.) - 4.5;
768 h4 = 10.5 - yq * (6. - yq * 6.);
769 h6 = yq * 4. - 6.;
770 e0 = yq * (yq * (yq + 5.5) + 8.25) + 1.875;
771 e2 = yq * (yq * 3. + 1.) + 5.25;
772 e4 = h6 * 0.75;
773 }
774 return rrtpi / (h0 + xq * (h2 + xq * (h4 + xq * (h6 + xq))))
775 * y * (e0 + xq * (e2 + xq * (e4 + xq)));
776 }
777
778 else if (abx < xlim3) { /* Humlicek W4 Region 3 */
779 if (rg3) { /* First point in Region 3 */
780 rg3 = false;
781 /* Region 3 y-dependents */
782 z0 = y * (y * (y * (y * (y * (y * (y * (y * (y * (y
783 + 13.3988) + 88.26741) + 369.1989) + 1074.409)
784 + 2256.981) + 3447.629) + 3764.966) + 2802.87)
785 + 1280.829) + 272.1014;
786 z2 = y * (y * (y * (y * (y * (y * (y * (y * 5. + 53.59518)
787 + 266.2987) + 793.4273) + 1549.675) + 2037.31)
788 + 1758.336) + 902.3066) + 211.678;
789 z4 = y * (y * (y * (y * (y * (y * 10. + 80.39278) + 269.2916)
790 + 479.2576) + 497.3014) + 308.1852) + 78.86585;
791 z6 = y * (y * (y * (y * 10. + 53.59518) + 92.75679)
792 + 55.02933) + 22.03523;
793 z8 = y * (y * 5. + 13.3988) + 1.49646;
794 p0 = y * (y * (y * (y * (y * (y * (y * (y * (y * 0.3183291
795 + 4.264678) + 27.93941) + 115.3772) + 328.2151) +
796 662.8097) + 946.897) + 919.4955) + 549.3954)
797 + 153.5168;
798 p2 = y * (y * (y * (y * (y * (y * (y * 1.2733163 + 12.79458)
799 + 56.81652) + 139.4665) + 189.773) + 124.5975)
800 - 1.322256) - 34.16955;
801 p4 = y * (y * (y * (y * (y * 1.9099744 + 12.79568)
802 + 29.81482) + 24.01655) + 10.46332) + 2.584042;
803 p6 = y * (y * (y * 1.273316 + 4.266322) + 0.9377051)
804 - 0.07272979;
805 p8 = y * .3183291 + 5.480304e-4;
806 }
807 return 1.7724538 / (z0 + xq * (z2 + xq * (z4 + xq * (z6 +
808 xq * (z8 + xq)))))
809 * (p0 + xq * (p2 + xq * (p4 + xq * (p6 + xq * p8))));
810 }
811
812 else { /* Humlicek CPF12 algorithm */
813 double ypy0 = y + 1.5;
814 double ypy0q = ypy0 * ypy0;
815 int j;
816 for (j = 0; j <= 5; ++j) {
817 double d = x - t[j];
818 mq[j] = d * d;
819 mf[j] = 1. / (mq[j] + ypy0q);
820 xm[j] = mf[j] * d;
821 ym[j] = mf[j] * ypy0;
822 d = x + t[j];
823 pq[j] = d * d;
824 pf[j] = 1. / (pq[j] + ypy0q);
825 xp[j] = pf[j] * d;
826 yp[j] = pf[j] * ypy0;
827 }
828 k = 0.;
829 if (abx <= xlim4) /* Humlicek CPF12 Region I */
830 for (j = 0; j <= 5; ++j)
831 k += c[j] * (ym[j]+yp[j]) - s[j] * (xm[j]-xp[j]);
832 else { /* Humlicek CPF12 Region II */
833 double yf = y + 3.;
834 for (j = 0; j <= 5; ++j)
835 k += (c[j] * (mq[j] * mf[j] - ym[j] * 1.5)
836 + s[j] * yf * xm[j]) / (mq[j] + 2.25)
837 + (c[j] * (pq[j] * pf[j] - yp[j] * 1.5)
838 - s[j] * yf * xp[j]) / (pq[j] + 2.25);
839 k = y * k + exp(-xq);
840 }
841 return k;
842 }
843 }
844 #endif /* libcerf not available */
845
846 /* ** ibeta.c
847 *
848 * DESCRIBE Approximate the incomplete beta function Ix(a, b).
849 *
850 * _
851 * |(a + b) /x (a-1) (b-1)
852 * Ix(a, b) = -_-------_--- * | t * (1 - t) dt (a,b > 0)
853 * |(a) * |(b) /0
854 *
855 *
856 *
857 * CALL p = ibeta(a, b, x)
858 *
859 * double a > 0
860 * double b > 0
861 * double x [0, 1]
862 *
863 * WARNING none
864 *
865 * RETURN double p [0, 1]
866 * -1.0 on error condition
867 *
868 * XREF lngamma()
869 *
870 * BUGS This approximation is only accurate on the domain
871 * x < (a-1)/(a+b-2)
872 *
873 * REFERENCE The continued fraction expansion as given by
874 * Abramowitz and Stegun (1964) is used.
875 *
876 * Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl
877 *
878 * Note: this function was translated from the Public Domain Fortran
879 * version available from http://lib.stat.cmu.edu/apstat/xxx
880 *
881 */
882
883 static double
ibeta(double a,double b,double x)884 ibeta(double a, double b, double x)
885 {
886 /* Test for admissibility of arguments */
887 if (a <= 0.0 || b <= 0.0)
888 return -1.0;
889 if (x < 0.0 || x > 1.0)
890 return -1.0;;
891
892 /* If x equals 0 or 1, return x as prob */
893 if (x == 0.0 || x == 1.0)
894 return x;
895
896 /* Swap a, b if necessary for more efficient evaluation */
897 if (a < x * (a + b)) {
898 double temp = confrac(b, a, 1.0 - x);
899 return (temp < 0.0) ? temp : 1.0 - temp;
900 } else {
901 return confrac(a, b, x);
902 }
903 }
904
905 static double
confrac(double a,double b,double x)906 confrac(double a, double b, double x)
907 {
908 double Alo = 0.0;
909 double Ahi;
910 double Aev;
911 double Aod;
912 double Blo = 1.0;
913 double Bhi = 1.0;
914 double Bod = 1.0;
915 double Bev = 1.0;
916 double f;
917 double fold;
918 double Apb = a + b;
919 double d;
920 int i;
921 int j;
922
923 /* Set up continued fraction expansion evaluation. */
924 Ahi = gp_exp(GAMMA(Apb) + a * log(x) + b * log(1.0 - x) -
925 GAMMA(a + 1.0) - GAMMA(b));
926
927 /*
928 * Continued fraction loop begins here. Evaluation continues until
929 * maximum iterations are exceeded, or convergence achieved.
930 */
931 for (i = 0, j = 1, f = Ahi; i <= ITMAX; i++, j++) {
932 d = a + j + i;
933 Aev = -(a + i) * (Apb + i) * x / d / (d - 1.0);
934 Aod = j * (b - j) * x / d / (d + 1.0);
935 Alo = Bev * Ahi + Aev * Alo;
936 Blo = Bev * Bhi + Aev * Blo;
937 Ahi = Bod * Alo + Aod * Ahi;
938 Bhi = Bod * Blo + Aod * Bhi;
939
940 if (fabs(Bhi) < MACHEPS)
941 Bhi = 0.0;
942
943 if (Bhi != 0.0) {
944 fold = f;
945 f = Ahi / Bhi;
946 if (fabs(f - fold) < fabs(f) * MACHEPS)
947 return f;
948 }
949 }
950
951 return -1.0;
952 }
953
954 /* ** igamma.c
955 *
956 * DESCRIBE Approximate the incomplete gamma function P(a, x).
957 *
958 * 1 /x -t (a-1)
959 * P(a, x) = -_--- * | e * t dt (a > 0)
960 * |(a) /0
961 *
962 * CALL p = igamma(a, x)
963 *
964 * double a > 0
965 * double x >= 0
966 *
967 * WARNING none
968 *
969 * RETURN double p [0, 1]
970 * -1.0 on error condition
971 *
972 * XREF lngamma()
973 *
974 * BUGS Values 0 <= x <= 1 may lead to inaccurate results.
975 *
976 * REFERENCE ALGORITHM AS239 APPL. STATIST. (1988) VOL. 37, NO. 3
977 *
978 * Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl
979 *
980 * Note: this function was translated from the Public Domain Fortran
981 * version available from http://lib.stat.cmu.edu/apstat/239
982 *
983 */
984
985 double
igamma(double a,double x)986 igamma(double a, double x)
987 {
988 double arg;
989 double aa;
990 double an;
991 double b;
992 int i;
993
994 /* Check that we have valid values for a and x */
995 if (x < 0.0 || a <= 0.0)
996 return -1.0;
997
998 /* Deal with special cases */
999 if (x == 0.0)
1000 return 0.0;
1001 if (x > XBIG)
1002 return 1.0;
1003
1004 /* Check value of factor arg */
1005 arg = a * log(x) - x - GAMMA(a + 1.0);
1006 /* HBB 20031006: removed a spurious check here */
1007 arg = gp_exp(arg);
1008
1009 /* Choose infinite series or continued fraction. */
1010
1011 if ((x > 1.0) && (x >= a + 2.0)) {
1012 /* Use a continued fraction expansion */
1013 double pn1, pn2, pn3, pn4, pn5, pn6;
1014 double rn;
1015 double rnold;
1016
1017 aa = 1.0 - a;
1018 b = aa + x + 1.0;
1019 pn1 = 1.0;
1020 pn2 = x;
1021 pn3 = x + 1.0;
1022 pn4 = x * b;
1023 rnold = pn3 / pn4;
1024
1025 for (i = 1; i <= ITMAX; i++) {
1026
1027 aa++;
1028 b += 2.0;
1029 an = aa * (double) i;
1030
1031 pn5 = b * pn3 - an * pn1;
1032 pn6 = b * pn4 - an * pn2;
1033
1034 if (pn6 != 0.0) {
1035
1036 rn = pn5 / pn6;
1037 if (fabs(rnold - rn) <= GPMIN(MACHEPS, MACHEPS * rn))
1038 return 1.0 - arg * rn * a;
1039
1040 rnold = rn;
1041 }
1042 pn1 = pn3;
1043 pn2 = pn4;
1044 pn3 = pn5;
1045 pn4 = pn6;
1046
1047 /* Re-scale terms in continued fraction if terms are large */
1048 if (fabs(pn5) >= OFLOW) {
1049
1050 pn1 /= OFLOW;
1051 pn2 /= OFLOW;
1052 pn3 /= OFLOW;
1053 pn4 /= OFLOW;
1054 }
1055 }
1056 } else {
1057 /* Use Pearson's series expansion. */
1058
1059 for (i = 0, aa = a, an = b = 1.0; i <= ITMAX; i++) {
1060
1061 aa++;
1062 an *= x / aa;
1063 b += an;
1064 if (an < b * MACHEPS)
1065 return arg * b;
1066 }
1067 }
1068 return -1.0;
1069 }
1070
1071
1072 /* ----------------------------------------------------------------
1073 Cummulative distribution function of the ChiSquare distribution
1074 ---------------------------------------------------------------- */
1075 double
chisq_cdf(int dof,double chisqr)1076 chisq_cdf(int dof, double chisqr)
1077 {
1078 if (dof <= 0)
1079 return not_a_number();
1080 if (chisqr <= 0.)
1081 return 0;
1082 return igamma(0.5 * dof, 0.5 * chisqr);
1083 }
1084
1085
1086 /***********************************************************************
1087 double ranf(double init)
1088 Random number generator as a Function
1089 Returns a random floating point number from a uniform distribution
1090 over 0 - 1 (endpoints of this interval are not returned) using a
1091 large integer generator.
1092 This is a transcription from Pascal to Fortran of routine
1093 Uniform_01 from the paper
1094 L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
1095 with Splitting Facilities." ACM Transactions on Mathematical
1096 Software, 17:98-111 (1991)
1097
1098 Generate Large Integer
1099 Returns a random integer following a uniform distribution over
1100 (1, 2147483562) using the generator.
1101 This is a transcription from Pascal to Fortran of routine
1102 Random from the paper
1103 L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
1104 with Splitting Facilities." ACM Transactions on Mathematical
1105 Software, 17:98-111 (1991)
1106 ***********************************************************************/
1107 static double
ranf(struct value * init)1108 ranf(struct value *init)
1109 {
1110 long k, z;
1111 static int firsttime = 1;
1112 static long seed1, seed2;
1113 static const long Xm1 = 2147483563L;
1114 static const long Xm2 = 2147483399L;
1115 static const long Xa1 = 40014L;
1116 static const long Xa2 = 40692L;
1117
1118 /* Seed values must be integer, but check for both values equal zero
1119 before casting for speed */
1120 if (real(init) != 0.0 || imag(init) != 0.0) {
1121
1122 /* Construct new seed values from input parameter */
1123 long seed1cvrt = real(init);
1124 long seed2cvrt = imag(init);
1125 if ( real(init) != (double)seed1cvrt ||
1126 imag(init) != (double)seed2cvrt ||
1127 seed1cvrt > 017777777777L ||
1128 seed2cvrt > 017777777777L ||
1129 (seed1cvrt <= 0 && seed2cvrt != 0) ||
1130 seed2cvrt < 0 )
1131 int_error(NO_CARET,"Illegal seed value");
1132 else if (seed1cvrt < 0)
1133 firsttime = 1;
1134 else {
1135 seed1 = seed1cvrt;
1136 seed2 = (seed2cvrt) ? seed2cvrt : seed1cvrt;
1137 firsttime = 0;
1138 }
1139 }
1140
1141 /* (Re)-Initialize seeds if necessary */
1142 if (firsttime) {
1143 firsttime = 0;
1144 seed1 = 1234567890L;
1145 seed2 = 1234567890L;
1146 }
1147
1148 FPRINTF((stderr,"ranf: seed = %lo %lo %ld %ld\n", seed1, seed2));
1149
1150 /* Generate pseudo random integers, which always end up positive */
1151 k = seed1 / 53668L;
1152 seed1 = Xa1 * (seed1 - k * 53668L) - k * 12211;
1153 if (seed1 < 0)
1154 seed1 += Xm1;
1155 k = seed2 / 52774L;
1156 seed2 = Xa2 * (seed2 - k * 52774L) - k * 3791;
1157 if (seed2 < 0)
1158 seed2 += Xm2;
1159 z = seed1 - seed2;
1160 if (z < 1)
1161 z += (Xm1 - 1);
1162
1163 /*
1164 * 4.656613057E-10 is 1/Xm1. Xm1 is set at the top of this file and is
1165 * currently 2147483563. If Xm1 changes, change this also.
1166 */
1167 return (double) 4.656613057E-10 *z;
1168 }
1169
1170 /* ----------------------------------------------------------------
1171 Following to specfun.c made by John Grosh (jgrosh@arl.mil)
1172 on 28 OCT 1992.
1173 ---------------------------------------------------------------- */
1174
1175 void
f_normal(union argument * arg)1176 f_normal(union argument *arg)
1177 { /* Normal or Gaussian Probability Function */
1178 struct value a;
1179 double x;
1180
1181 /* ref. Abramowitz and Stegun 1964, "Handbook of Mathematical
1182 Functions", Applied Mathematics Series, vol 55,
1183 Chapter 26, page 934, Eqn. 26.2.29 and Jos van der Woude
1184 code found above */
1185
1186 (void) arg; /* avoid -Wunused warning */
1187 x = real(pop(&a));
1188
1189 /* using erfc instead of erf produces accurate values for -38 < arg < -8 */
1190 if (x > -38) {
1191 x = 0.5 * SQRT_TWO * x;
1192 x = 0.5 * erfc(-x);
1193 } else {
1194 x = 0.0;
1195 }
1196 push(Gcomplex(&a, x, 0.0));
1197 }
1198
1199 void
f_inverse_normal(union argument * arg)1200 f_inverse_normal(union argument *arg)
1201 { /* Inverse normal distribution function */
1202 struct value a;
1203 double x;
1204
1205 (void) arg; /* avoid -Wunused warning */
1206 x = real(pop(&a));
1207
1208 if (x <= 0.0 || x >= 1.0) {
1209 undefined = TRUE;
1210 push(Gcomplex(&a, 0.0, 0.0));
1211 } else {
1212 push(Gcomplex(&a, inverse_normal_func(x), 0.0));
1213 }
1214 }
1215
1216
1217 void
f_inverse_erf(union argument * arg)1218 f_inverse_erf(union argument *arg)
1219 { /* Inverse error function */
1220 struct value a;
1221 double x;
1222
1223 (void) arg; /* avoid -Wunused warning */
1224 x = real(pop(&a));
1225
1226 if (fabs(x) >= 1.0) {
1227 undefined = TRUE;
1228 push(Gcomplex(&a, 0.0, 0.0));
1229 } else {
1230 push(Gcomplex(&a, inverse_error_func(x), 0.0));
1231 }
1232 }
1233
1234 /* ndtri.c
1235 *
1236 * Inverse of Normal distribution function
1237 *
1238 *
1239 *
1240 * SYNOPSIS:
1241 *
1242 * double x, y, ndtri();
1243 *
1244 * x = ndtri( y );
1245 *
1246 *
1247 *
1248 * DESCRIPTION:
1249 *
1250 * Returns the argument, x, for which the area under the
1251 * Gaussian probability density function (integrated from
1252 * minus infinity to x) is equal to y.
1253 *
1254 *
1255 * For small arguments 0 < y < exp(-2), the program computes
1256 * z = sqrt( -2.0 * log(y) ); then the approximation is
1257 * x = z - log(z)/z - (1/z) P(1/z) / Q(1/z).
1258 * There are two rational functions P/Q, one for 0 < y < exp(-32)
1259 * and the other for y up to exp(-2). For larger arguments,
1260 * w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
1261 *
1262 *
1263 * ACCURACY:
1264 *
1265 * Relative error:
1266 * arithmetic domain # trials peak rms
1267 * DEC 0.125, 1 5500 9.5e-17 2.1e-17
1268 * DEC 6e-39, 0.135 3500 5.7e-17 1.3e-17
1269 * IEEE 0.125, 1 20000 7.2e-16 1.3e-16
1270 * IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17
1271 *
1272 *
1273 * ERROR MESSAGES:
1274 *
1275 * message condition value returned
1276 * ndtri domain x <= 0 -DBL_MAX
1277 * ndtri domain x >= 1 DBL_MAX
1278 *
1279 */
1280
1281 /*
1282 Cephes Math Library Release 2.8: June, 2000
1283 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
1284 */
1285
1286 #ifdef UNK
1287 /* sqrt(2pi) */
1288 static double s2pi = 2.50662827463100050242E0;
1289 #endif
1290
1291 #ifdef DEC
1292 static unsigned short s2p[] = {0040440, 0066230, 0177661, 0034055};
1293 #define s2pi *(double *)s2p
1294 #endif
1295
1296 #ifdef IBMPC
1297 static unsigned short s2p[] = {0x2706, 0x1ff6, 0x0d93, 0x4004};
1298 #define s2pi *(double *)s2p
1299 #endif
1300
1301 #ifdef MIEEE
1302 static unsigned short s2p[] = {
1303 0x4004, 0x0d93, 0x1ff6, 0x2706
1304 };
1305 #define s2pi *(double *)s2p
1306 #endif
1307
1308 static double
inverse_normal_func(double y0)1309 inverse_normal_func(double y0)
1310 {
1311 /* approximation for 0 <= |y - 0.5| <= 3/8 */
1312 #ifdef UNK
1313 static const double P0[5] = {
1314 -5.99633501014107895267E1,
1315 9.80010754185999661536E1,
1316 -5.66762857469070293439E1,
1317 1.39312609387279679503E1,
1318 -1.23916583867381258016E0,
1319 };
1320 static const double Q0[8] = {
1321 /* 1.00000000000000000000E0,*/
1322 1.95448858338141759834E0,
1323 4.67627912898881538453E0,
1324 8.63602421390890590575E1,
1325 -2.25462687854119370527E2,
1326 2.00260212380060660359E2,
1327 -8.20372256168333339912E1,
1328 1.59056225126211695515E1,
1329 -1.18331621121330003142E0,
1330 };
1331 #endif
1332 #ifdef DEC
1333 static const unsigned short P0[20] = {
1334 0141557, 0155170, 0071360, 0120550,
1335 0041704, 0000214, 0172417, 0067307,
1336 0141542, 0132204, 0040066, 0156723,
1337 0041136, 0163161, 0157276, 0007747,
1338 0140236, 0116374, 0073666, 0051764,
1339 };
1340 static const unsigned short Q0[32] = {
1341 /*0040200,0000000,0000000,0000000,*/
1342 0040372, 0026256, 0110403, 0123707,
1343 0040625, 0122024, 0020277, 0026661,
1344 0041654, 0134161, 0124134, 0007244,
1345 0142141, 0073162, 0133021, 0131371,
1346 0042110, 0041235, 0043516, 0057767,
1347 0141644, 0011417, 0036155, 0137305,
1348 0041176, 0076556, 0004043, 0125430,
1349 0140227, 0073347, 0152776, 0067251,
1350 };
1351 #endif
1352 #ifdef IBMPC
1353 static const unsigned short P0[20] = {
1354 0x142d, 0x0e5e, 0xfb4f, 0xc04d,
1355 0xedd9, 0x9ea1, 0x8011, 0x4058,
1356 0xdbba, 0x8806, 0x5690, 0xc04c,
1357 0xc1fd, 0x3bd7, 0xdcce, 0x402b,
1358 0xca7e, 0x8ef6, 0xd39f, 0xbff3,
1359 };
1360 static const unsigned short Q0[36] = {
1361 /*0x0000,0x0000,0x0000,0x3ff0,*/
1362 0x74f9, 0xd220, 0x4595, 0x3fff,
1363 0xe5b6, 0x8417, 0xb482, 0x4012,
1364 0x81d4, 0x350b, 0x970e, 0x4055,
1365 0x365f, 0x56c2, 0x2ece, 0xc06c,
1366 0xcbff, 0xa8e9, 0x0853, 0x4069,
1367 0xb7d9, 0xe78d, 0x8261, 0xc054,
1368 0x7563, 0xc104, 0xcfad, 0x402f,
1369 0xcdd5, 0xfabf, 0xeedc, 0xbff2,
1370 };
1371 #endif
1372 #ifdef MIEEE
1373 static const unsigned short P0[20] = {
1374 0xc04d, 0xfb4f, 0x0e5e, 0x142d,
1375 0x4058, 0x8011, 0x9ea1, 0xedd9,
1376 0xc04c, 0x5690, 0x8806, 0xdbba,
1377 0x402b, 0xdcce, 0x3bd7, 0xc1fd,
1378 0xbff3, 0xd39f, 0x8ef6, 0xca7e,
1379 };
1380 static const unsigned short Q0[32] = {
1381 /*0x3ff0,0x0000,0x0000,0x0000,*/
1382 0x3fff, 0x4595, 0xd220, 0x74f9,
1383 0x4012, 0xb482, 0x8417, 0xe5b6,
1384 0x4055, 0x970e, 0x350b, 0x81d4,
1385 0xc06c, 0x2ece, 0x56c2, 0x365f,
1386 0x4069, 0x0853, 0xa8e9, 0xcbff,
1387 0xc054, 0x8261, 0xe78d, 0xb7d9,
1388 0x402f, 0xcfad, 0xc104, 0x7563,
1389 0xbff2, 0xeedc, 0xfabf, 0xcdd5,
1390 };
1391 #endif
1392
1393 /* Approximation for interval z = sqrt(-2 log y ) between 2 and 8
1394 * i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
1395 */
1396 #ifdef UNK
1397 static const double P1[9] = {
1398 4.05544892305962419923E0,
1399 3.15251094599893866154E1,
1400 5.71628192246421288162E1,
1401 4.40805073893200834700E1,
1402 1.46849561928858024014E1,
1403 2.18663306850790267539E0,
1404 -1.40256079171354495875E-1,
1405 -3.50424626827848203418E-2,
1406 -8.57456785154685413611E-4,
1407 };
1408 static const double Q1[8] = {
1409 /* 1.00000000000000000000E0,*/
1410 1.57799883256466749731E1,
1411 4.53907635128879210584E1,
1412 4.13172038254672030440E1,
1413 1.50425385692907503408E1,
1414 2.50464946208309415979E0,
1415 -1.42182922854787788574E-1,
1416 -3.80806407691578277194E-2,
1417 -9.33259480895457427372E-4,
1418 };
1419 #endif
1420 #ifdef DEC
1421 static const unsigned short P1[36] = {
1422 0040601, 0143074, 0150744, 0073326,
1423 0041374, 0031554, 0113253, 0146016,
1424 0041544, 0123272, 0012463, 0176771,
1425 0041460, 0051160, 0103560, 0156511,
1426 0041152, 0172624, 0117772, 0030755,
1427 0040413, 0170713, 0151545, 0176413,
1428 0137417, 0117512, 0022154, 0131671,
1429 0137017, 0104257, 0071432, 0007072,
1430 0135540, 0143363, 0063137, 0036166,
1431 };
1432 static const unsigned short Q1[32] = {
1433 /*0040200,0000000,0000000,0000000,*/
1434 0041174, 0075325, 0004736, 0120326,
1435 0041465, 0110044, 0047561, 0045567,
1436 0041445, 0042321, 0012142, 0030340,
1437 0041160, 0127074, 0166076, 0141051,
1438 0040440, 0046055, 0040745, 0150400,
1439 0137421, 0114146, 0067330, 0010621,
1440 0137033, 0175162, 0025555, 0114351,
1441 0135564, 0122773, 0145750, 0030357,
1442 };
1443 #endif
1444 #ifdef IBMPC
1445 static const unsigned short P1[36] = {
1446 0x8edb, 0x9a3c, 0x38c7, 0x4010,
1447 0x7982, 0x92d5, 0x866d, 0x403f,
1448 0x7fbf, 0x42a6, 0x94d7, 0x404c,
1449 0x1ba9, 0x10ee, 0x0a4e, 0x4046,
1450 0x463e, 0x93ff, 0x5eb2, 0x402d,
1451 0xbfa1, 0x7a6c, 0x7e39, 0x4001,
1452 0x9677, 0x448d, 0xf3e9, 0xbfc1,
1453 0x41c7, 0xee63, 0xf115, 0xbfa1,
1454 0xe78f, 0x6ccb, 0x18de, 0xbf4c,
1455 };
1456 static const unsigned short Q1[32] = {
1457 /*0x0000,0x0000,0x0000,0x3ff0,*/
1458 0xd41b, 0xa13b, 0x8f5a, 0x402f,
1459 0x296f, 0x89ee, 0xb204, 0x4046,
1460 0x461c, 0x228c, 0xa89a, 0x4044,
1461 0xd845, 0x9d87, 0x15c7, 0x402e,
1462 0xba20, 0xa83c, 0x0985, 0x4004,
1463 0x0232, 0xcddb, 0x330c, 0xbfc2,
1464 0xb31d, 0x456d, 0x7f4e, 0xbfa3,
1465 0x061e, 0x797d, 0x94bf, 0xbf4e,
1466 };
1467 #endif
1468 #ifdef MIEEE
1469 static const unsigned short P1[36] = {
1470 0x4010, 0x38c7, 0x9a3c, 0x8edb,
1471 0x403f, 0x866d, 0x92d5, 0x7982,
1472 0x404c, 0x94d7, 0x42a6, 0x7fbf,
1473 0x4046, 0x0a4e, 0x10ee, 0x1ba9,
1474 0x402d, 0x5eb2, 0x93ff, 0x463e,
1475 0x4001, 0x7e39, 0x7a6c, 0xbfa1,
1476 0xbfc1, 0xf3e9, 0x448d, 0x9677,
1477 0xbfa1, 0xf115, 0xee63, 0x41c7,
1478 0xbf4c, 0x18de, 0x6ccb, 0xe78f,
1479 };
1480 static const unsigned short Q1[32] = {
1481 /*0x3ff0,0x0000,0x0000,0x0000,*/
1482 0x402f, 0x8f5a, 0xa13b, 0xd41b,
1483 0x4046, 0xb204, 0x89ee, 0x296f,
1484 0x4044, 0xa89a, 0x228c, 0x461c,
1485 0x402e, 0x15c7, 0x9d87, 0xd845,
1486 0x4004, 0x0985, 0xa83c, 0xba20,
1487 0xbfc2, 0x330c, 0xcddb, 0x0232,
1488 0xbfa3, 0x7f4e, 0x456d, 0xb31d,
1489 0xbf4e, 0x94bf, 0x797d, 0x061e,
1490 };
1491 #endif
1492
1493 /* Approximation for interval z = sqrt(-2 log y ) between 8 and 64
1494 * i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
1495 */
1496
1497 #ifdef UNK
1498 static const double P2[9] = {
1499 3.23774891776946035970E0,
1500 6.91522889068984211695E0,
1501 3.93881025292474443415E0,
1502 1.33303460815807542389E0,
1503 2.01485389549179081538E-1,
1504 1.23716634817820021358E-2,
1505 3.01581553508235416007E-4,
1506 2.65806974686737550832E-6,
1507 6.23974539184983293730E-9,
1508 };
1509 static const double Q2[8] = {
1510 /* 1.00000000000000000000E0,*/
1511 6.02427039364742014255E0,
1512 3.67983563856160859403E0,
1513 1.37702099489081330271E0,
1514 2.16236993594496635890E-1,
1515 1.34204006088543189037E-2,
1516 3.28014464682127739104E-4,
1517 2.89247864745380683936E-6,
1518 6.79019408009981274425E-9,
1519 };
1520 #endif
1521 #ifdef DEC
1522 static const unsigned short P2[36] = {
1523 0040517, 0033507, 0036236, 0125641,
1524 0040735, 0044616, 0014473, 0140133,
1525 0040574, 0012567, 0114535, 0102541,
1526 0040252, 0120340, 0143474, 0150135,
1527 0037516, 0051057, 0115361, 0031211,
1528 0036512, 0131204, 0101511, 0125144,
1529 0035236, 0016627, 0043160, 0140216,
1530 0033462, 0060512, 0060141, 0010641,
1531 0031326, 0062541, 0101304, 0077706,
1532 };
1533 static const unsigned short Q2[32] = {
1534 /*0040200,0000000,0000000,0000000,*/
1535 0040700, 0143322, 0132137, 0040501,
1536 0040553, 0101155, 0053221, 0140257,
1537 0040260, 0041071, 0052573, 0010004,
1538 0037535, 0066472, 0177261, 0162330,
1539 0036533, 0160475, 0066666, 0036132,
1540 0035253, 0174533, 0027771, 0044027,
1541 0033502, 0016147, 0117666, 0063671,
1542 0031351, 0047455, 0141663, 0054751,
1543 };
1544 #endif
1545 #ifdef IBMPC
1546 static const unsigned short P2[36] = {
1547 0xd574, 0xe793, 0xe6e8, 0x4009,
1548 0x780b, 0xc327, 0xa931, 0x401b,
1549 0xb0ac, 0xf32b, 0x82ae, 0x400f,
1550 0x9a0c, 0x18e7, 0x541c, 0x3ff5,
1551 0x2651, 0xf35e, 0xca45, 0x3fc9,
1552 0x354d, 0x9069, 0x5650, 0x3f89,
1553 0x1812, 0xe8ce, 0xc3b2, 0x3f33,
1554 0x2234, 0x4c0c, 0x4c29, 0x3ec6,
1555 0x8ff9, 0x3058, 0xccac, 0x3e3a,
1556 };
1557 static const unsigned short Q2[32] = {
1558 /*0x0000,0x0000,0x0000,0x3ff0,*/
1559 0xe828, 0x568b, 0x18da, 0x4018,
1560 0x3816, 0xaad2, 0x704d, 0x400d,
1561 0x6200, 0x2aaf, 0x0847, 0x3ff6,
1562 0x3c9b, 0x5fd6, 0xada7, 0x3fcb,
1563 0xc78b, 0xadb6, 0x7c27, 0x3f8b,
1564 0x2903, 0x65ff, 0x7f2b, 0x3f35,
1565 0xccf7, 0xf3f6, 0x438c, 0x3ec8,
1566 0x6b3d, 0xb876, 0x29e5, 0x3e3d,
1567 };
1568 #endif
1569 #ifdef MIEEE
1570 static const unsigned short P2[36] = {
1571 0x4009, 0xe6e8, 0xe793, 0xd574,
1572 0x401b, 0xa931, 0xc327, 0x780b,
1573 0x400f, 0x82ae, 0xf32b, 0xb0ac,
1574 0x3ff5, 0x541c, 0x18e7, 0x9a0c,
1575 0x3fc9, 0xca45, 0xf35e, 0x2651,
1576 0x3f89, 0x5650, 0x9069, 0x354d,
1577 0x3f33, 0xc3b2, 0xe8ce, 0x1812,
1578 0x3ec6, 0x4c29, 0x4c0c, 0x2234,
1579 0x3e3a, 0xccac, 0x3058, 0x8ff9,
1580 };
1581 static const unsigned short Q2[32] = {
1582 /*0x3ff0,0x0000,0x0000,0x0000,*/
1583 0x4018, 0x18da, 0x568b, 0xe828,
1584 0x400d, 0x704d, 0xaad2, 0x3816,
1585 0x3ff6, 0x0847, 0x2aaf, 0x6200,
1586 0x3fcb, 0xada7, 0x5fd6, 0x3c9b,
1587 0x3f8b, 0x7c27, 0xadb6, 0xc78b,
1588 0x3f35, 0x7f2b, 0x65ff, 0x2903,
1589 0x3ec8, 0x438c, 0xf3f6, 0xccf7,
1590 0x3e3d, 0x29e5, 0xb876, 0x6b3d,
1591 };
1592 #endif
1593
1594 double x, y, z, y2, x0, x1;
1595 int code;
1596
1597 if (y0 <= 0.0) {
1598 mtherr("inverse_normal_func", MTHERR_DOMAIN);
1599 return (-DBL_MAX);
1600 }
1601 if (y0 >= 1.0) {
1602 mtherr("inverse_normal_func", MTHERR_DOMAIN);
1603 return (DBL_MAX);
1604 }
1605 code = 1;
1606 y = y0;
1607 if (y > (1.0 - 0.13533528323661269189)) { /* 0.135... = exp(-2) */
1608 y = 1.0 - y;
1609 code = 0;
1610 }
1611 if (y > 0.13533528323661269189) {
1612 y = y - 0.5;
1613 y2 = y * y;
1614 x = y + y * (y2 * polevl(y2, P0, 4) / p1evl(y2, Q0, 8));
1615 x = x * s2pi;
1616 return (x);
1617 }
1618 x = sqrt(-2.0 * log(y));
1619 x0 = x - log(x) / x;
1620
1621 z = 1.0 / x;
1622 if (x < 8.0) /* y > exp(-32) = 1.2664165549e-14 */
1623 x1 = z * polevl(z, P1, 8) / p1evl(z, Q1, 8);
1624 else
1625 x1 = z * polevl(z, P2, 8) / p1evl(z, Q2, 8);
1626 x = x0 - x1;
1627 if (code != 0)
1628 x = -x;
1629 return (x);
1630 }
1631
1632 /*
1633 Cephes Math Library Release 2.8: June, 2000
1634 Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
1635 */
1636
1637 #ifndef HAVE_ERFC
1638 /* erfc.c
1639 *
1640 * Complementary error function
1641 *
1642 *
1643 *
1644 * SYNOPSIS:
1645 *
1646 * double x, y, erfc();
1647 *
1648 * y = erfc( x );
1649 *
1650 *
1651 *
1652 * DESCRIPTION:
1653 *
1654 *
1655 * 1 - erf(x) =
1656 *
1657 * inf.
1658 * -
1659 * 2 | | 2
1660 * erfc(x) = -------- | exp( - t ) dt
1661 * sqrt(pi) | |
1662 * -
1663 * x
1664 *
1665 *
1666 * For small x, erfc(x) = 1 - erf(x); otherwise rational
1667 * approximations are computed.
1668 *
1669 *
1670 *
1671 * ACCURACY:
1672 *
1673 * Relative error:
1674 * arithmetic domain # trials peak rms
1675 * DEC 0, 9.2319 12000 5.1e-16 1.2e-16
1676 * IEEE 0,26.6417 30000 5.7e-14 1.5e-14
1677 *
1678 *
1679 * ERROR MESSAGES:
1680 *
1681 * message condition value returned
1682 * erfc underflow x > 9.231948545 (DEC) 0.0
1683 *
1684 *
1685 */
1686
1687 static double
erfc(double a)1688 erfc(double a)
1689 {
1690 #ifdef UNK
1691 static const double P[] = {
1692 2.46196981473530512524E-10,
1693 5.64189564831068821977E-1,
1694 7.46321056442269912687E0,
1695 4.86371970985681366614E1,
1696 1.96520832956077098242E2,
1697 5.26445194995477358631E2,
1698 9.34528527171957607540E2,
1699 1.02755188689515710272E3,
1700 5.57535335369399327526E2
1701 };
1702 static const double Q[] = {
1703 /* 1.00000000000000000000E0,*/
1704 1.32281951154744992508E1,
1705 8.67072140885989742329E1,
1706 3.54937778887819891062E2,
1707 9.75708501743205489753E2,
1708 1.82390916687909736289E3,
1709 2.24633760818710981792E3,
1710 1.65666309194161350182E3,
1711 5.57535340817727675546E2
1712 };
1713 static const double R[] = {
1714 5.64189583547755073984E-1,
1715 1.27536670759978104416E0,
1716 5.01905042251180477414E0,
1717 6.16021097993053585195E0,
1718 7.40974269950448939160E0,
1719 2.97886665372100240670E0
1720 };
1721 static const double S[] = {
1722 /* 1.00000000000000000000E0,*/
1723 2.26052863220117276590E0,
1724 9.39603524938001434673E0,
1725 1.20489539808096656605E1,
1726 1.70814450747565897222E1,
1727 9.60896809063285878198E0,
1728 3.36907645100081516050E0
1729 };
1730 #endif /* UNK */
1731
1732 #ifdef DEC
1733 static const unsigned short P[] = {
1734 0030207, 0054445, 0011173, 0021706,
1735 0040020, 0067272, 0030661, 0122075,
1736 0040756, 0151236, 0173053, 0067042,
1737 0041502, 0106175, 0062555, 0151457,
1738 0042104, 0102525, 0047401, 0003667,
1739 0042403, 0116176, 0011446, 0075303,
1740 0042551, 0120723, 0061641, 0123275,
1741 0042600, 0070651, 0007264, 0134516,
1742 0042413, 0061102, 0167507, 0176625
1743 };
1744 static const unsigned short Q[] = {
1745 /*0040200,0000000,0000000,0000000,*/
1746 0041123, 0123257, 0165741, 0017142,
1747 0041655, 0065027, 0173413, 0115450,
1748 0042261, 0074011, 0021573, 0004150,
1749 0042563, 0166530, 0013662, 0007200,
1750 0042743, 0176427, 0162443, 0105214,
1751 0043014, 0062546, 0153727, 0123772,
1752 0042717, 0012470, 0006227, 0067424,
1753 0042413, 0061103, 0003042, 0013254
1754 };
1755 static const unsigned short R[] = {
1756 0040020, 0067272, 0101024, 0155421,
1757 0040243, 0037467, 0056706, 0026462,
1758 0040640, 0116017, 0120665, 0034315,
1759 0040705, 0020162, 0143350, 0060137,
1760 0040755, 0016234, 0134304, 0130157,
1761 0040476, 0122700, 0051070, 0015473
1762 };
1763 static const unsigned short S[] = {
1764 /*0040200,0000000,0000000,0000000,*/
1765 0040420, 0126200, 0044276, 0070413,
1766 0041026, 0053051, 0007302, 0063746,
1767 0041100, 0144203, 0174051, 0061151,
1768 0041210, 0123314, 0126343, 0177646,
1769 0041031, 0137125, 0051431, 0033011,
1770 0040527, 0117362, 0152661, 0066201
1771 };
1772 #endif /* DEC */
1773
1774 #ifdef IBMPC
1775 static const unsigned short P[] = {
1776 0x6479, 0xa24f, 0xeb24, 0x3df0,
1777 0x3488, 0x4636, 0x0dd7, 0x3fe2,
1778 0x6dc4, 0xdec5, 0xda53, 0x401d,
1779 0xba66, 0xacad, 0x518f, 0x4048,
1780 0x20f7, 0xa9e0, 0x90aa, 0x4068,
1781 0xcf58, 0xc264, 0x738f, 0x4080,
1782 0x34d8, 0x6c74, 0x343a, 0x408d,
1783 0x972a, 0x21d6, 0x0e35, 0x4090,
1784 0xffb3, 0x5de8, 0x6c48, 0x4081
1785 };
1786 static const unsigned short Q[] = {
1787 /*0x0000,0x0000,0x0000,0x3ff0,*/
1788 0x23cc, 0xfd7c, 0x74d5, 0x402a,
1789 0x7365, 0xfee1, 0xad42, 0x4055,
1790 0x610d, 0x246f, 0x2f01, 0x4076,
1791 0x41d0, 0x02f6, 0x7dab, 0x408e,
1792 0x7151, 0xfca4, 0x7fa2, 0x409c,
1793 0xf4ff, 0xdafa, 0x8cac, 0x40a1,
1794 0xede2, 0x0192, 0xe2a7, 0x4099,
1795 0x42d6, 0x60c4, 0x6c48, 0x4081
1796 };
1797 static const unsigned short R[] = {
1798 0x9b62, 0x5042, 0x0dd7, 0x3fe2,
1799 0xc5a6, 0xebb8, 0x67e6, 0x3ff4,
1800 0xa71a, 0xf436, 0x1381, 0x4014,
1801 0x0c0c, 0x58dd, 0xa40e, 0x4018,
1802 0x960e, 0x9718, 0xa393, 0x401d,
1803 0x0367, 0x0a47, 0xd4b8, 0x4007
1804 };
1805 static const unsigned short S[] = {
1806 /*0x0000,0x0000,0x0000,0x3ff0,*/
1807 0xce21, 0x0917, 0x1590, 0x4002,
1808 0x4cfd, 0x21d8, 0xcac5, 0x4022,
1809 0x2c4d, 0x7f05, 0x1910, 0x4028,
1810 0x7ff5, 0x959c, 0x14d9, 0x4031,
1811 0x26c1, 0xaa63, 0x37ca, 0x4023,
1812 0x2d90, 0x5ab6, 0xf3de, 0x400a
1813 };
1814 #endif /* IBMPC */
1815
1816 #ifdef MIEEE
1817 static const unsigned short P[] = {
1818 0x3df0, 0xeb24, 0xa24f, 0x6479,
1819 0x3fe2, 0x0dd7, 0x4636, 0x3488,
1820 0x401d, 0xda53, 0xdec5, 0x6dc4,
1821 0x4048, 0x518f, 0xacad, 0xba66,
1822 0x4068, 0x90aa, 0xa9e0, 0x20f7,
1823 0x4080, 0x738f, 0xc264, 0xcf58,
1824 0x408d, 0x343a, 0x6c74, 0x34d8,
1825 0x4090, 0x0e35, 0x21d6, 0x972a,
1826 0x4081, 0x6c48, 0x5de8, 0xffb3
1827 };
1828 static const unsigned short Q[] = {
1829 0x402a, 0x74d5, 0xfd7c, 0x23cc,
1830 0x4055, 0xad42, 0xfee1, 0x7365,
1831 0x4076, 0x2f01, 0x246f, 0x610d,
1832 0x408e, 0x7dab, 0x02f6, 0x41d0,
1833 0x409c, 0x7fa2, 0xfca4, 0x7151,
1834 0x40a1, 0x8cac, 0xdafa, 0xf4ff,
1835 0x4099, 0xe2a7, 0x0192, 0xede2,
1836 0x4081, 0x6c48, 0x60c4, 0x42d6
1837 };
1838 static const unsigned short R[] = {
1839 0x3fe2, 0x0dd7, 0x5042, 0x9b62,
1840 0x3ff4, 0x67e6, 0xebb8, 0xc5a6,
1841 0x4014, 0x1381, 0xf436, 0xa71a,
1842 0x4018, 0xa40e, 0x58dd, 0x0c0c,
1843 0x401d, 0xa393, 0x9718, 0x960e,
1844 0x4007, 0xd4b8, 0x0a47, 0x0367
1845 };
1846 static const unsigned short S[] = {
1847 0x4002, 0x1590, 0x0917, 0xce21,
1848 0x4022, 0xcac5, 0x21d8, 0x4cfd,
1849 0x4028, 0x1910, 0x7f05, 0x2c4d,
1850 0x4031, 0x14d9, 0x959c, 0x7ff5,
1851 0x4023, 0x37ca, 0xaa63, 0x26c1,
1852 0x400a, 0xf3de, 0x5ab6, 0x2d90
1853 };
1854 #endif /* MIEEE */
1855
1856 double p, q, x, y, z;
1857
1858 if (a < 0.0)
1859 x = -a;
1860 else
1861 x = a;
1862
1863 if (x < 1.0)
1864 return (1.0 - erf(a));
1865
1866 z = -a * a;
1867
1868 if (z < DBL_MIN_10_EXP) {
1869 under:
1870 mtherr("erfc", MTHERR_UNDERFLOW);
1871 if (a < 0)
1872 return (2.0);
1873 else
1874 return (0.0);
1875 }
1876 z = exp(z);
1877
1878 if (x < 8.0) {
1879 p = polevl(x, P, 8);
1880 q = p1evl(x, Q, 8);
1881 } else {
1882 p = polevl(x, R, 5);
1883 q = p1evl(x, S, 6);
1884 }
1885 y = (z * p) / q;
1886
1887 if (a < 0)
1888 y = 2.0 - y;
1889
1890 if (y == 0.0)
1891 goto under;
1892
1893 return (y);
1894 }
1895 #endif /* !HAVE_ERFC */
1896
1897 #ifndef HAVE_ERF
1898 /* erf.c
1899 *
1900 * Error function
1901 *
1902 *
1903 *
1904 * SYNOPSIS:
1905 *
1906 * double x, y, erf();
1907 *
1908 * y = erf( x );
1909 *
1910 *
1911 *
1912 * DESCRIPTION:
1913 *
1914 * The integral is
1915 *
1916 * x
1917 * -
1918 * 2 | | 2
1919 * erf(x) = -------- | exp( - t ) dt.
1920 * sqrt(pi) | |
1921 * -
1922 * 0
1923 *
1924 * The magnitude of x is limited to 9.231948545 for DEC
1925 * arithmetic; 1 or -1 is returned outside this range.
1926 *
1927 * For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
1928 * erf(x) = 1 - erfc(x).
1929 *
1930 *
1931 *
1932 * ACCURACY:
1933 *
1934 * Relative error:
1935 * arithmetic domain # trials peak rms
1936 * DEC 0,1 14000 4.7e-17 1.5e-17
1937 * IEEE 0,1 30000 3.7e-16 1.0e-16
1938 *
1939 */
1940
1941 static double
erf(double x)1942 erf(double x)
1943 {
1944
1945 # ifdef UNK
1946 static const double T[] = {
1947 9.60497373987051638749E0,
1948 9.00260197203842689217E1,
1949 2.23200534594684319226E3,
1950 7.00332514112805075473E3,
1951 5.55923013010394962768E4
1952 };
1953 static const double U[] = {
1954 /* 1.00000000000000000000E0,*/
1955 3.35617141647503099647E1,
1956 5.21357949780152679795E2,
1957 4.59432382970980127987E3,
1958 2.26290000613890934246E4,
1959 4.92673942608635921086E4
1960 };
1961 # endif
1962
1963 # ifdef DEC
1964 static const unsigned short T[] = {
1965 0041031, 0126770, 0170672, 0166101,
1966 0041664, 0006522, 0072360, 0031770,
1967 0043013, 0100025, 0162641, 0126671,
1968 0043332, 0155231, 0161627, 0076200,
1969 0044131, 0024115, 0021020, 0117343
1970 };
1971 static const unsigned short U[] = {
1972 /*0040200,0000000,0000000,0000000,*/
1973 0041406, 0037461, 0177575, 0032714,
1974 0042402, 0053350, 0123061, 0153557,
1975 0043217, 0111227, 0032007, 0164217,
1976 0043660, 0145000, 0004013, 0160114,
1977 0044100, 0071544, 0167107, 0125471
1978 };
1979 # endif
1980
1981 # ifdef IBMPC
1982 static const unsigned short T[] = {
1983 0x5d88, 0x1e37, 0x35bf, 0x4023,
1984 0x067f, 0x4e9e, 0x81aa, 0x4056,
1985 0x35b7, 0xbcb4, 0x7002, 0x40a1,
1986 0xef90, 0x3c72, 0x5b53, 0x40bb,
1987 0x13dc, 0xa442, 0x2509, 0x40eb
1988 };
1989 static const unsigned short U[] = {
1990 /*0x0000,0x0000,0x0000,0x3ff0,*/
1991 0xa6ba, 0x3fef, 0xc7e6, 0x4040,
1992 0x3aee, 0x14c6, 0x4add, 0x4080,
1993 0xfd12, 0xe680, 0xf252, 0x40b1,
1994 0x7c0a, 0x0101, 0x1940, 0x40d6,
1995 0xf567, 0x9dc8, 0x0e6c, 0x40e8
1996 };
1997 # endif
1998
1999 # ifdef MIEEE
2000 static const unsigned short T[] = {
2001 0x4023, 0x35bf, 0x1e37, 0x5d88,
2002 0x4056, 0x81aa, 0x4e9e, 0x067f,
2003 0x40a1, 0x7002, 0xbcb4, 0x35b7,
2004 0x40bb, 0x5b53, 0x3c72, 0xef90,
2005 0x40eb, 0x2509, 0xa442, 0x13dc
2006 };
2007 static const unsigned short U[] = {
2008 0x4040, 0xc7e6, 0x3fef, 0xa6ba,
2009 0x4080, 0x4add, 0x14c6, 0x3aee,
2010 0x40b1, 0xf252, 0xe680, 0xfd12,
2011 0x40d6, 0x1940, 0x0101, 0x7c0a,
2012 0x40e8, 0x0e6c, 0x9dc8, 0xf567
2013 };
2014 # endif
2015
2016 double y, z;
2017
2018 if (fabs(x) > 1.0)
2019 return (1.0 - erfc(x));
2020 z = x * x;
2021 y = x * polevl(z, T, 4) / p1evl(z, U, 5);
2022 return (y);
2023 }
2024 #endif /* !HAVE_ERF */
2025
2026 /* ----------------------------------------------------------------
2027 Following function for the inverse error function is taken from
2028 NIST on 16. May 2002.
2029 Use Newton-Raphson correction also for range -1 to -y0 and
2030 add 3rd cycle to improve convergence - E A Merritt 21.10.2003
2031 ----------------------------------------------------------------
2032 */
2033
2034 static double
inverse_error_func(double y)2035 inverse_error_func(double y)
2036 {
2037 double x = 0.0; /* The output */
2038 double z = 0.0; /* Intermadiate variable */
2039 double y0 = 0.7; /* Central range variable */
2040
2041 /* Coefficients in rational approximations. */
2042 static const double a[4] = {
2043 0.886226899, -1.645349621, 0.914624893, -0.140543331
2044 };
2045 static const double b[4] = {
2046 -2.118377725, 1.442710462, -0.329097515, 0.012229801
2047 };
2048 static const double c[4] = {
2049 -1.970840454, -1.624906493, 3.429567803, 1.641345311
2050 };
2051 static const double d[2] = {
2052 3.543889200, 1.637067800
2053 };
2054
2055 if ((y < -1.0) || (1.0 < y)) {
2056 printf("inverse_error_func: The value out of the range of the function");
2057 x = log(-1.0);
2058 return (x);
2059 } else if ((y == -1.0) || (1.0 == y)) {
2060 x = -y * log(0.0);
2061 return (x);
2062 } else if ((-1.0 < y) && (y < -y0)) {
2063 z = sqrt(-log((1.0 + y) / 2.0));
2064 x = -(((c[3] * z + c[2]) * z + c[1]) * z + c[0]) / ((d[1] * z + d[0]) * z + 1.0);
2065 } else {
2066 if ((-y0 <= y) && (y <= y0)) {
2067 z = y * y;
2068 x = y * (((a[3] * z + a[2]) * z + a[1]) * z + a[0]) /
2069 ((((b[3] * z + b[3]) * z + b[1]) * z + b[0]) * z + 1.0);
2070 } else if ((y0 < y) && (y < 1.0)) {
2071 z = sqrt(-log((1.0 - y) / 2.0));
2072 x = (((c[3] * z + c[2]) * z + c[1]) * z + c[0]) / ((d[1] * z + d[0]) * z + 1.0);
2073 }
2074 }
2075 /* Three steps of Newton-Raphson correction to full accuracy. OK - four */
2076 x = x - (erf(x) - y) / (2.0 / sqrt(PI) * gp_exp(-x * x));
2077 x = x - (erf(x) - y) / (2.0 / sqrt(PI) * gp_exp(-x * x));
2078 x = x - (erf(x) - y) / (2.0 / sqrt(PI) * gp_exp(-x * x));
2079 x = x - (erf(x) - y) / (2.0 / sqrt(PI) * gp_exp(-x * x));
2080
2081 return (x);
2082 }
2083
2084
2085 /* Implementation of Lamberts W-function which is defined as
2086 * w(x)*e^(w(x))=x
2087 * Implementation by Gunter Kuhnle, gk@uni-leipzig.de
2088 * Algorithm originally developed by
2089 * KEITH BRIGGS, DEPARTMENT OF PLANT SCIENCES,
2090 * e-mail:kmb28@cam.ac.uk
2091 * http://epidem13.plantsci.cam.ac.uk/~kbriggs/W-ology.html */
2092
2093 static double
lambertw(double x)2094 lambertw(double x)
2095 {
2096 double p, e, t, w, eps;
2097 int i;
2098
2099 eps = MACHEPS;
2100
2101 if (x < -exp(-1))
2102 return -1; /* error, value undefined */
2103
2104 if (fabs(x) <= eps)
2105 return x;
2106
2107 if (x < 1) {
2108 p = sqrt(2.0 * (exp(1.0) * x + 1.0));
2109 w = -1.0 + p - p * p / 3.0 + 11.0 / 72.0 * p * p * p;
2110 } else {
2111 w = log(x);
2112 }
2113
2114 if (x > 3) {
2115 w = w - log(w);
2116 }
2117 for (i = 0; i < 20; i++) {
2118 e = gp_exp(w);
2119 t = w * e - x;
2120 t = t / (e * (w + 1.0) - 0.5 * (w + 2.0) * t / (w + 1.0));
2121 w = w - t;
2122 if (fabs(t) < eps * (1.0 + fabs(w)))
2123 return w;
2124 }
2125 return -1; /* error: iteration didn't converge */
2126 }
2127
2128 void
f_lambertw(union argument * arg)2129 f_lambertw(union argument *arg)
2130 {
2131 struct value a;
2132 double x;
2133
2134 (void) arg; /* avoid -Wunused warning */
2135 x = real(pop(&a));
2136
2137 x = lambertw(x);
2138 if (x <= -1)
2139 /* Error return from lambertw --> flag 'undefined' */
2140 undefined = TRUE;
2141
2142 push(Gcomplex(&a, x, 0.0));
2143 }
2144
2145 #if (0) /* This approximation of the airy function saves 2200 bytes relative
2146 * to the one from Cephes, but has low precision (2% relative error)
2147 */
2148 /* ------------------------------------------------------------
2149 Airy Function Ai(x)
2150
2151 After:
2152 "Two-Point Quasi-Fractional Approximations to the Airy Function Ai(x)"
2153 by Pablo Martin, Ricardo Perez, Antonio L. Guerrero
2154 Journal of Computational Physics 99, 337-340 (1992)
2155
2156 Beware of a misprint in equation (5) in this paper: The second term in
2157 parentheses must be multiplied by "x", as is clear from equation (3)
2158 and by comparison with equation (6). The implementation in this file
2159 uses the CORRECT formula (with the "x").
2160
2161 This is not a very high accuracy approximation, but sufficient for
2162 plotting and similar applications. Higher accuracy formulas are
2163 available, but are much more complicated (typically requiring iteration).
2164
2165 Added: janert (PKJ) 2009-09-05
2166 ------------------------------------------------------------ */
2167
2168 static double
airy_neg(double x)2169 airy_neg( double x ) {
2170 double z = sqrt( 0.37 + pow( fabs(x), 3.0 ) );
2171 double t = (2.0/3.0)*pow( fabs(x), 1.5 );
2172 double y = 0;
2173
2174 y += ( -0.043883564 + 0.3989422*z )*cos(t)/pow( z, 7.0/6.0 );
2175 y += x*( -0.013883003 - 0.3989422*z )*sin(t)/( pow( z, 5.0/6.0 ) * 1.5 * t );
2176
2177 return y;
2178 }
2179
2180 static double
airy_pos(double x)2181 airy_pos( double x ) {
2182 double z = sqrt( 0.0425 + pow( fabs(x), 3.0 ) );
2183 double y = 0;
2184
2185 y += (-0.002800908 + 0.326662423*z )/pow( z, 7.0/6.0 );
2186 y += x * ( -0.007232251 - 0.044567423*z )/pow( z, 11.0/6.0 );
2187 y *= exp(-(2.0/3.0)*z );
2188
2189 return y;
2190 }
2191
2192 void
f_airy(union argument * arg)2193 f_airy(union argument *arg)
2194 {
2195 struct value a;
2196 double x;
2197
2198 (void) arg; /* avoid -Wunused warning */
2199 x = real(pop(&a));
2200
2201 if( x < 0 ) {
2202 x = airy_neg(x);
2203 } else {
2204 x = airy_pos(x);
2205 }
2206
2207 push(Gcomplex(&a, x, 0.0));
2208 }
2209 #else /* Airy function from the Cephes library */
2210 /* airy.c
2211 *
2212 * Airy function
2213 *
2214 * SYNOPSIS:
2215 *
2216 * double x, ai, aip, bi, bip;
2217 * int airy();
2218 *
2219 * airy( x, _&ai, _&aip, _&bi, _&bip );
2220 *
2221 *
2222 * DESCRIPTION:
2223 *
2224 * Solution of the differential equation
2225 *
2226 * y"(x) = xy.
2227 *
2228 * The function returns the two independent solutions Ai, Bi
2229 * and their first derivatives Ai'(x), Bi'(x).
2230 *
2231 * Evaluation is by power series summation for small x,
2232 * by rational minimax approximations for large x.
2233 *
2234 *
2235 * ACCURACY:
2236 * Error criterion is absolute when function <= 1, relative
2237 * when function > 1, except * denotes relative error criterion.
2238 * For large negative x, the absolute error increases as x^1.5.
2239 * For large positive x, the relative error increases as x^1.5.
2240 *
2241 * Arithmetic domain function # trials peak rms
2242 * IEEE -10, 0 Ai 10000 1.6e-15 2.7e-16
2243 * IEEE 0, 10 Ai 10000 2.3e-14* 1.8e-15*
2244 * IEEE -10, 0 Ai' 10000 4.6e-15 7.6e-16
2245 * IEEE 0, 10 Ai' 10000 1.8e-14* 1.5e-15*
2246 * IEEE -10, 10 Bi 30000 4.2e-15 5.3e-16
2247 * IEEE -10, 10 Bi' 30000 4.9e-15 7.3e-16
2248 * DEC -10, 0 Ai 5000 1.7e-16 2.8e-17
2249 * DEC 0, 10 Ai 5000 2.1e-15* 1.7e-16*
2250 * DEC -10, 0 Ai' 5000 4.7e-16 7.8e-17
2251 * DEC 0, 10 Ai' 12000 1.8e-15* 1.5e-16*
2252 * DEC -10, 10 Bi 10000 5.5e-16 6.8e-17
2253 * DEC -10, 10 Bi' 7000 5.3e-16 8.7e-17
2254 *
2255 */
2256
2257 /*
2258 Cephes Math Library Release 2.8: June, 2000
2259 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
2260 */
2261
2262 static double c1 = 0.35502805388781723926;
2263 static double c2 = 0.258819403792806798405;
2264 static double sqrt3 = 1.732050807568877293527;
2265 static double sqpii = 5.64189583547756286948E-1;
2266
2267 #ifdef UNK
2268 #define MAXAIRY 25.77
2269 #endif
2270 #ifdef DEC
2271 #define MAXAIRY 25.77
2272 #endif
2273 #ifdef IBMPC
2274 #define MAXAIRY 103.892
2275 #endif
2276 #ifdef MIEEE
2277 #define MAXAIRY 103.892
2278 #endif
2279
2280
2281 #ifdef UNK
2282 static double AN[8] = {
2283 3.46538101525629032477E-1,
2284 1.20075952739645805542E1,
2285 7.62796053615234516538E1,
2286 1.68089224934630576269E2,
2287 1.59756391350164413639E2,
2288 7.05360906840444183113E1,
2289 1.40264691163389668864E1,
2290 9.99999999999999995305E-1,
2291 };
2292 static double AD[8] = {
2293 5.67594532638770212846E-1,
2294 1.47562562584847203173E1,
2295 8.45138970141474626562E1,
2296 1.77318088145400459522E2,
2297 1.64234692871529701831E2,
2298 7.14778400825575695274E1,
2299 1.40959135607834029598E1,
2300 1.00000000000000000470E0,
2301 };
2302 #endif
2303 #ifdef DEC
2304 static unsigned short AN[32] = {
2305 0037661,0066561,0024675,0131301,
2306 0041100,0017434,0034324,0101466,
2307 0041630,0107450,0067427,0007430,
2308 0042050,0013327,0071000,0034737,
2309 0042037,0140642,0156417,0167366,
2310 0041615,0011172,0075147,0051165,
2311 0041140,0066152,0160520,0075146,
2312 0040200,0000000,0000000,0000000,
2313 };
2314 static unsigned short AD[32] = {
2315 0040021,0046740,0011422,0064606,
2316 0041154,0014640,0024631,0062450,
2317 0041651,0003435,0101152,0106401,
2318 0042061,0050556,0034605,0136602,
2319 0042044,0036024,0152377,0151414,
2320 0041616,0172247,0072216,0115374,
2321 0041141,0104334,0124154,0166007,
2322 0040200,0000000,0000000,0000000,
2323 };
2324 #endif
2325 #ifdef IBMPC
2326 static unsigned short AN[32] = {
2327 0xb658,0x2537,0x2dae,0x3fd6,
2328 0x9067,0x871a,0x03e3,0x4028,
2329 0xe1e3,0x0de2,0x11e5,0x4053,
2330 0x073c,0xee40,0x02da,0x4065,
2331 0xfddf,0x5ba1,0xf834,0x4063,
2332 0xea4f,0x4f4c,0xa24f,0x4051,
2333 0x0f4d,0x5c2a,0x0d8d,0x402c,
2334 0x0000,0x0000,0x0000,0x3ff0,
2335 };
2336 static unsigned short AD[32] = {
2337 0x4d31,0x0262,0x29bc,0x3fe2,
2338 0x2ca5,0x0533,0x8334,0x402d,
2339 0x51a0,0xb04d,0x20e3,0x4055,
2340 0xb7b0,0xc730,0x2a2d,0x4066,
2341 0xfa61,0x9a9f,0x8782,0x4064,
2342 0xd35f,0xee91,0xde94,0x4051,
2343 0x9d81,0x950d,0x311b,0x402c,
2344 0x0000,0x0000,0x0000,0x3ff0,
2345 };
2346 #endif
2347 #ifdef MIEEE
2348 static unsigned short AN[32] = {
2349 0x3fd6,0x2dae,0x2537,0xb658,
2350 0x4028,0x03e3,0x871a,0x9067,
2351 0x4053,0x11e5,0x0de2,0xe1e3,
2352 0x4065,0x02da,0xee40,0x073c,
2353 0x4063,0xf834,0x5ba1,0xfddf,
2354 0x4051,0xa24f,0x4f4c,0xea4f,
2355 0x402c,0x0d8d,0x5c2a,0x0f4d,
2356 0x3ff0,0x0000,0x0000,0x0000,
2357 };
2358 static unsigned short AD[32] = {
2359 0x3fe2,0x29bc,0x0262,0x4d31,
2360 0x402d,0x8334,0x0533,0x2ca5,
2361 0x4055,0x20e3,0xb04d,0x51a0,
2362 0x4066,0x2a2d,0xc730,0xb7b0,
2363 0x4064,0x8782,0x9a9f,0xfa61,
2364 0x4051,0xde94,0xee91,0xd35f,
2365 0x402c,0x311b,0x950d,0x9d81,
2366 0x3ff0,0x0000,0x0000,0x0000,
2367 };
2368 #endif
2369
2370 #ifdef UNK
2371 static double APN[8] = {
2372 6.13759184814035759225E-1,
2373 1.47454670787755323881E1,
2374 8.20584123476060982430E1,
2375 1.71184781360976385540E2,
2376 1.59317847137141783523E2,
2377 6.99778599330103016170E1,
2378 1.39470856980481566958E1,
2379 1.00000000000000000550E0,
2380 };
2381 static double APD[8] = {
2382 3.34203677749736953049E-1,
2383 1.11810297306158156705E1,
2384 7.11727352147859965283E1,
2385 1.58778084372838313640E2,
2386 1.53206427475809220834E2,
2387 6.86752304592780337944E1,
2388 1.38498634758259442477E1,
2389 9.99999999999999994502E-1,
2390 };
2391 #endif
2392 #ifdef DEC
2393 static unsigned short APN[32] = {
2394 0040035,0017522,0065145,0054755,
2395 0041153,0166556,0161471,0057174,
2396 0041644,0016750,0034445,0046462,
2397 0042053,0027515,0152316,0046717,
2398 0042037,0050536,0067023,0023264,
2399 0041613,0172252,0007240,0131055,
2400 0041137,0023503,0052472,0002305,
2401 0040200,0000000,0000000,0000000,
2402 };
2403 static unsigned short APD[32] = {
2404 0037653,0016276,0112106,0126625,
2405 0041062,0162577,0067111,0111761,
2406 0041616,0054160,0140004,0137455,
2407 0042036,0143460,0104626,0157206,
2408 0042031,0032330,0067131,0114260,
2409 0041611,0054667,0147207,0134564,
2410 0041135,0114412,0070653,0146015,
2411 0040200,0000000,0000000,0000000,
2412 };
2413 #endif
2414 #ifdef IBMPC
2415 static unsigned short APN[32] = {
2416 0xab3e,0x4d4c,0xa3ea,0x3fe3,
2417 0x2bcf,0xdc67,0x7dad,0x402d,
2418 0xa9a6,0x0724,0x83bd,0x4054,
2419 0xc9ba,0xba99,0x65e9,0x4065,
2420 0x64d7,0xcdc2,0xea2b,0x4063,
2421 0x1646,0x41d4,0x7e95,0x4051,
2422 0x4099,0x6aa7,0xe4e8,0x402b,
2423 0x0000,0x0000,0x0000,0x3ff0,
2424 };
2425 static unsigned short APD[32] = {
2426 0xd5b3,0xd288,0x6397,0x3fd5,
2427 0x327e,0xedc9,0x5caf,0x4026,
2428 0x97e6,0x1800,0xcb0e,0x4051,
2429 0xdbd1,0x1132,0xd8e6,0x4063,
2430 0x3316,0x0dcb,0x269b,0x4063,
2431 0xf72f,0xf9d0,0x2b36,0x4051,
2432 0x7982,0x4e35,0xb321,0x402b,
2433 0x0000,0x0000,0x0000,0x3ff0,
2434 };
2435 #endif
2436 #ifdef MIEEE
2437 static unsigned short APN[32] = {
2438 0x3fe3,0xa3ea,0x4d4c,0xab3e,
2439 0x402d,0x7dad,0xdc67,0x2bcf,
2440 0x4054,0x83bd,0x0724,0xa9a6,
2441 0x4065,0x65e9,0xba99,0xc9ba,
2442 0x4063,0xea2b,0xcdc2,0x64d7,
2443 0x4051,0x7e95,0x41d4,0x1646,
2444 0x402b,0xe4e8,0x6aa7,0x4099,
2445 0x3ff0,0x0000,0x0000,0x0000,
2446 };
2447 static unsigned short APD[32] = {
2448 0x3fd5,0x6397,0xd288,0xd5b3,
2449 0x4026,0x5caf,0xedc9,0x327e,
2450 0x4051,0xcb0e,0x1800,0x97e6,
2451 0x4063,0xd8e6,0x1132,0xdbd1,
2452 0x4063,0x269b,0x0dcb,0x3316,
2453 0x4051,0x2b36,0xf9d0,0xf72f,
2454 0x402b,0xb321,0x4e35,0x7982,
2455 0x3ff0,0x0000,0x0000,0x0000,
2456 };
2457 #endif
2458
2459 #ifdef UNK
2460 static double BN16[5] = {
2461 -2.53240795869364152689E-1,
2462 5.75285167332467384228E-1,
2463 -3.29907036873225371650E-1,
2464 6.44404068948199951727E-2,
2465 -3.82519546641336734394E-3,
2466 };
2467 static double BD16[5] = {
2468 /* 1.00000000000000000000E0,*/
2469 -7.15685095054035237902E0,
2470 1.06039580715664694291E1,
2471 -5.23246636471251500874E0,
2472 9.57395864378383833152E-1,
2473 -5.50828147163549611107E-2,
2474 };
2475 #endif
2476 #ifdef DEC
2477 static unsigned short BN16[20] = {
2478 0137601,0124307,0010213,0035210,
2479 0040023,0042743,0101621,0016031,
2480 0137650,0164623,0036056,0074511,
2481 0037203,0174525,0000473,0142474,
2482 0136172,0130041,0066726,0064324,
2483 };
2484 static unsigned short BD16[20] = {
2485 /*0040200,0000000,0000000,0000000,*/
2486 0140745,0002354,0044335,0055276,
2487 0041051,0124717,0170130,0104013,
2488 0140647,0070135,0046473,0103501,
2489 0040165,0013745,0033324,0127766,
2490 0137141,0117204,0076164,0033107,
2491 };
2492 #endif
2493 #ifdef IBMPC
2494 static unsigned short BN16[20] = {
2495 0x6751,0xe211,0x3518,0xbfd0,
2496 0x2383,0x7072,0x68bc,0x3fe2,
2497 0xcf29,0x6785,0x1d32,0xbfd5,
2498 0x78a8,0xa027,0x7f2a,0x3fb0,
2499 0xcd1b,0x2dba,0x5604,0xbf6f,
2500 };
2501 static unsigned short BD16[20] = {
2502 /*0x0000,0x0000,0x0000,0x3ff0,*/
2503 0xab58,0x891b,0xa09d,0xc01c,
2504 0x1101,0xfe0b,0x3539,0x4025,
2505 0x70e8,0xa9a7,0xee0b,0xc014,
2506 0x95ff,0xa6da,0xa2fc,0x3fee,
2507 0x86c9,0x8f8e,0x33d0,0xbfac,
2508 };
2509 #endif
2510 #ifdef MIEEE
2511 static unsigned short BN16[20] = {
2512 0xbfd0,0x3518,0xe211,0x6751,
2513 0x3fe2,0x68bc,0x7072,0x2383,
2514 0xbfd5,0x1d32,0x6785,0xcf29,
2515 0x3fb0,0x7f2a,0xa027,0x78a8,
2516 0xbf6f,0x5604,0x2dba,0xcd1b,
2517 };
2518 static unsigned short BD16[20] = {
2519 /*0x3ff0,0x0000,0x0000,0x0000,*/
2520 0xc01c,0xa09d,0x891b,0xab58,
2521 0x4025,0x3539,0xfe0b,0x1101,
2522 0xc014,0xee0b,0xa9a7,0x70e8,
2523 0x3fee,0xa2fc,0xa6da,0x95ff,
2524 0xbfac,0x33d0,0x8f8e,0x86c9,
2525 };
2526 #endif
2527
2528 #ifdef UNK
2529 static double BPPN[5] = {
2530 4.65461162774651610328E-1,
2531 -1.08992173800493920734E0,
2532 6.38800117371827987759E-1,
2533 -1.26844349553102907034E-1,
2534 7.62487844342109852105E-3,
2535 };
2536 static double BPPD[5] = {
2537 /* 1.00000000000000000000E0,*/
2538 -8.70622787633159124240E0,
2539 1.38993162704553213172E1,
2540 -7.14116144616431159572E0,
2541 1.34008595960680518666E0,
2542 -7.84273211323341930448E-2,
2543 };
2544 #endif
2545 #ifdef DEC
2546 static unsigned short BPPN[20] = {
2547 0037756,0050354,0167531,0135731,
2548 0140213,0101216,0032767,0020375,
2549 0040043,0104147,0106312,0177632,
2550 0137401,0161574,0032015,0043714,
2551 0036371,0155035,0143165,0142262,
2552 };
2553 static unsigned short BPPD[20] = {
2554 /*0040200,0000000,0000000,0000000,*/
2555 0141013,0046265,0115005,0161053,
2556 0041136,0061631,0072445,0156131,
2557 0140744,0102145,0001127,0065304,
2558 0040253,0103757,0146453,0102513,
2559 0137240,0117200,0155402,0113500,
2560 };
2561 #endif
2562 #ifdef IBMPC
2563 static unsigned short BPPN[20] = {
2564 0x377b,0x9deb,0xca1d,0x3fdd,
2565 0xe420,0xc6be,0x7051,0xbff1,
2566 0x5ff3,0xf199,0x710c,0x3fe4,
2567 0xa8fa,0x8681,0x3c6f,0xbfc0,
2568 0xb896,0xb8ce,0x3b43,0x3f7f,
2569 };
2570 static unsigned short BPPD[20] = {
2571 /*0x0000,0x0000,0x0000,0x3ff0,*/
2572 0xbc45,0xb340,0x6996,0xc021,
2573 0xbb8b,0x2ea4,0xcc73,0x402b,
2574 0xed59,0xa04a,0x908c,0xc01c,
2575 0x70a9,0xf9a5,0x70fd,0x3ff5,
2576 0x52e8,0x1b60,0x13d0,0xbfb4,
2577 };
2578 #endif
2579 #ifdef MIEEE
2580 static unsigned short BPPN[20] = {
2581 0x3fdd,0xca1d,0x9deb,0x377b,
2582 0xbff1,0x7051,0xc6be,0xe420,
2583 0x3fe4,0x710c,0xf199,0x5ff3,
2584 0xbfc0,0x3c6f,0x8681,0xa8fa,
2585 0x3f7f,0x3b43,0xb8ce,0xb896,
2586 };
2587 static unsigned short BPPD[20] = {
2588 /*0x3ff0,0x0000,0x0000,0x0000,*/
2589 0xc021,0x6996,0xb340,0xbc45,
2590 0x402b,0xcc73,0x2ea4,0xbb8b,
2591 0xc01c,0x908c,0xa04a,0xed59,
2592 0x3ff5,0x70fd,0xf9a5,0x70a9,
2593 0xbfb4,0x13d0,0x1b60,0x52e8,
2594 };
2595 #endif
2596
2597 #ifdef UNK
2598 static double AFN[9] = {
2599 -1.31696323418331795333E-1,
2600 -6.26456544431912369773E-1,
2601 -6.93158036036933542233E-1,
2602 -2.79779981545119124951E-1,
2603 -4.91900132609500318020E-2,
2604 -4.06265923594885404393E-3,
2605 -1.59276496239262096340E-4,
2606 -2.77649108155232920844E-6,
2607 -1.67787698489114633780E-8,
2608 };
2609 static double AFD[9] = {
2610 /* 1.00000000000000000000E0,*/
2611 1.33560420706553243746E1,
2612 3.26825032795224613948E1,
2613 2.67367040941499554804E1,
2614 9.18707402907259625840E0,
2615 1.47529146771666414581E0,
2616 1.15687173795188044134E-1,
2617 4.40291641615211203805E-3,
2618 7.54720348287414296618E-5,
2619 4.51850092970580378464E-7,
2620 };
2621 #endif
2622 #ifdef DEC
2623 static unsigned short AFN[36] = {
2624 0137406,0155546,0124127,0033732,
2625 0140040,0057564,0141263,0041222,
2626 0140061,0071316,0013674,0175754,
2627 0137617,0037522,0056637,0120130,
2628 0137111,0075567,0121755,0166122,
2629 0136205,0020016,0043317,0002201,
2630 0135047,0001565,0075130,0002334,
2631 0133472,0051700,0165021,0131551,
2632 0131620,0020347,0132165,0013215,
2633 };
2634 static unsigned short AFD[36] = {
2635 /*0040200,0000000,0000000,0000000,*/
2636 0041125,0131131,0025627,0067623,
2637 0041402,0135342,0021703,0154315,
2638 0041325,0162305,0016671,0120175,
2639 0041022,0177101,0053114,0141632,
2640 0040274,0153131,0147364,0114306,
2641 0037354,0166545,0120042,0150530,
2642 0036220,0043127,0000727,0130273,
2643 0034636,0043275,0075667,0034733,
2644 0032762,0112715,0146250,0142474,
2645 };
2646 #endif
2647 #ifdef IBMPC
2648 static unsigned short AFN[36] = {
2649 0xe6fb,0xd50a,0xdb6c,0xbfc0,
2650 0x6852,0x9856,0x0bee,0xbfe4,
2651 0x9f7d,0xc2f7,0x2e59,0xbfe6,
2652 0xf40b,0x4bb3,0xe7ea,0xbfd1,
2653 0xbd8a,0xf47d,0x2f6e,0xbfa9,
2654 0xe090,0xc8d9,0xa401,0xbf70,
2655 0x009c,0xaf4b,0xe06e,0xbf24,
2656 0x366d,0x1d42,0x4a78,0xbec7,
2657 0xa2d2,0xf68e,0x041c,0xbe52,
2658 };
2659 static unsigned short AFD[36] = {
2660 /*0x0000,0x0000,0x0000,0x3ff0,*/
2661 0xedf2,0x2572,0xb64b,0x402a,
2662 0x7b1a,0x4478,0x575c,0x4040,
2663 0x3410,0xa3b7,0xbc98,0x403a,
2664 0x9873,0x2ac9,0x5fc8,0x4022,
2665 0x9319,0x39de,0x9acb,0x3ff7,
2666 0x5a2b,0xb404,0x9dac,0x3fbd,
2667 0xf617,0xe03a,0x08ca,0x3f72,
2668 0xe73b,0xaf76,0xc8d7,0x3f13,
2669 0x18a7,0xb995,0x52b9,0x3e9e,
2670 };
2671 #endif
2672 #ifdef MIEEE
2673 static unsigned short AFN[36] = {
2674 0xbfc0,0xdb6c,0xd50a,0xe6fb,
2675 0xbfe4,0x0bee,0x9856,0x6852,
2676 0xbfe6,0x2e59,0xc2f7,0x9f7d,
2677 0xbfd1,0xe7ea,0x4bb3,0xf40b,
2678 0xbfa9,0x2f6e,0xf47d,0xbd8a,
2679 0xbf70,0xa401,0xc8d9,0xe090,
2680 0xbf24,0xe06e,0xaf4b,0x009c,
2681 0xbec7,0x4a78,0x1d42,0x366d,
2682 0xbe52,0x041c,0xf68e,0xa2d2,
2683 };
2684 static unsigned short AFD[36] = {
2685 /*0x3ff0,0x0000,0x0000,0x0000,*/
2686 0x402a,0xb64b,0x2572,0xedf2,
2687 0x4040,0x575c,0x4478,0x7b1a,
2688 0x403a,0xbc98,0xa3b7,0x3410,
2689 0x4022,0x5fc8,0x2ac9,0x9873,
2690 0x3ff7,0x9acb,0x39de,0x9319,
2691 0x3fbd,0x9dac,0xb404,0x5a2b,
2692 0x3f72,0x08ca,0xe03a,0xf617,
2693 0x3f13,0xc8d7,0xaf76,0xe73b,
2694 0x3e9e,0x52b9,0xb995,0x18a7,
2695 };
2696 #endif
2697
2698 #ifdef UNK
2699 static double AGN[11] = {
2700 1.97339932091685679179E-2,
2701 3.91103029615688277255E-1,
2702 1.06579897599595591108E0,
2703 9.39169229816650230044E-1,
2704 3.51465656105547619242E-1,
2705 6.33888919628925490927E-2,
2706 5.85804113048388458567E-3,
2707 2.82851600836737019778E-4,
2708 6.98793669997260967291E-6,
2709 8.11789239554389293311E-8,
2710 3.41551784765923618484E-10,
2711 };
2712 static double AGD[10] = {
2713 /* 1.00000000000000000000E0,*/
2714 9.30892908077441974853E0,
2715 1.98352928718312140417E1,
2716 1.55646628932864612953E1,
2717 5.47686069422975497931E0,
2718 9.54293611618961883998E-1,
2719 8.64580826352392193095E-2,
2720 4.12656523824222607191E-3,
2721 1.01259085116509135510E-4,
2722 1.17166733214413521882E-6,
2723 4.91834570062930015649E-9,
2724 };
2725 #endif
2726 #ifdef DEC
2727 static unsigned short AGN[44] = {
2728 0036641,0124456,0167175,0157354,
2729 0037710,0037250,0001441,0136671,
2730 0040210,0066031,0150401,0123532,
2731 0040160,0066545,0003570,0153133,
2732 0037663,0171516,0072507,0170345,
2733 0037201,0151011,0007510,0045702,
2734 0036277,0172317,0104572,0101030,
2735 0035224,0045663,0000160,0136422,
2736 0033752,0074753,0047702,0135160,
2737 0032256,0052225,0156550,0107103,
2738 0030273,0142443,0166277,0071720,
2739 };
2740 static unsigned short AGD[40] = {
2741 /*0040200,0000000,0000000,0000000,*/
2742 0041024,0170537,0117253,0055003,
2743 0041236,0127256,0003570,0143240,
2744 0041171,0004333,0172476,0160645,
2745 0040657,0041161,0055716,0157161,
2746 0040164,0046226,0006257,0063431,
2747 0037261,0010357,0065445,0047563,
2748 0036207,0034043,0057434,0116732,
2749 0034724,0055416,0130035,0026377,
2750 0033235,0041056,0154071,0023502,
2751 0031250,0177071,0167254,0047242,
2752 };
2753 #endif
2754 #ifdef IBMPC
2755 static unsigned short AGN[44] = {
2756 0xbbde,0xddcf,0x3525,0x3f94,
2757 0x37b7,0x0064,0x07d5,0x3fd9,
2758 0x34eb,0x3a20,0x0d83,0x3ff1,
2759 0x1acb,0xa0ef,0x0dac,0x3fee,
2760 0xfe1d,0xcea8,0x7e69,0x3fd6,
2761 0x0978,0x21e9,0x3a41,0x3fb0,
2762 0x5043,0xf12f,0xfe99,0x3f77,
2763 0x17a2,0x600e,0x8976,0x3f32,
2764 0x574e,0x69f8,0x4f3d,0x3edd,
2765 0x11c8,0xbbad,0xca92,0x3e75,
2766 0xee7a,0x7d97,0x78a4,0x3df7,
2767 };
2768 static unsigned short AGD[40] = {
2769 /*0x0000,0x0000,0x0000,0x3ff0,*/
2770 0x6b40,0xf3d5,0x9e2b,0x4022,
2771 0x18d4,0xc0ef,0xd5d5,0x4033,
2772 0xdc35,0x7ea7,0x211b,0x402f,
2773 0xdbce,0x2b79,0xe84e,0x4015,
2774 0xece3,0xc195,0x8992,0x3fee,
2775 0xa9ee,0xed64,0x221d,0x3fb6,
2776 0x93bb,0x6be3,0xe704,0x3f70,
2777 0xa5a0,0xd603,0x8b61,0x3f1a,
2778 0x24e8,0xdb07,0xa845,0x3eb3,
2779 0x89d4,0x3dd5,0x1fc7,0x3e35,
2780 };
2781 #endif
2782 #ifdef MIEEE
2783 static unsigned short AGN[44] = {
2784 0x3f94,0x3525,0xddcf,0xbbde,
2785 0x3fd9,0x07d5,0x0064,0x37b7,
2786 0x3ff1,0x0d83,0x3a20,0x34eb,
2787 0x3fee,0x0dac,0xa0ef,0x1acb,
2788 0x3fd6,0x7e69,0xcea8,0xfe1d,
2789 0x3fb0,0x3a41,0x21e9,0x0978,
2790 0x3f77,0xfe99,0xf12f,0x5043,
2791 0x3f32,0x8976,0x600e,0x17a2,
2792 0x3edd,0x4f3d,0x69f8,0x574e,
2793 0x3e75,0xca92,0xbbad,0x11c8,
2794 0x3df7,0x78a4,0x7d97,0xee7a,
2795 };
2796 static unsigned short AGD[40] = {
2797 /*0x3ff0,0x0000,0x0000,0x0000,*/
2798 0x4022,0x9e2b,0xf3d5,0x6b40,
2799 0x4033,0xd5d5,0xc0ef,0x18d4,
2800 0x402f,0x211b,0x7ea7,0xdc35,
2801 0x4015,0xe84e,0x2b79,0xdbce,
2802 0x3fee,0x8992,0xc195,0xece3,
2803 0x3fb6,0x221d,0xed64,0xa9ee,
2804 0x3f70,0xe704,0x6be3,0x93bb,
2805 0x3f1a,0x8b61,0xd603,0xa5a0,
2806 0x3eb3,0xa845,0xdb07,0x24e8,
2807 0x3e35,0x1fc7,0x3dd5,0x89d4,
2808 };
2809 #endif
2810
2811 #ifdef UNK
2812 static double APFN[9] = {
2813 1.85365624022535566142E-1,
2814 8.86712188052584095637E-1,
2815 9.87391981747398547272E-1,
2816 4.01241082318003734092E-1,
2817 7.10304926289631174579E-2,
2818 5.90618657995661810071E-3,
2819 2.33051409401776799569E-4,
2820 4.08718778289035454598E-6,
2821 2.48379932900442457853E-8,
2822 };
2823 static double APFD[9] = {
2824 /* 1.00000000000000000000E0,*/
2825 1.47345854687502542552E1,
2826 3.75423933435489594466E1,
2827 3.14657751203046424330E1,
2828 1.09969125207298778536E1,
2829 1.78885054766999417817E0,
2830 1.41733275753662636873E-1,
2831 5.44066067017226003627E-3,
2832 9.39421290654511171663E-5,
2833 5.65978713036027009243E-7,
2834 };
2835 #endif
2836 #ifdef DEC
2837 static unsigned short APFN[36] = {
2838 0037475,0150174,0071752,0166651,
2839 0040142,0177621,0164246,0101757,
2840 0040174,0142670,0106760,0006573,
2841 0037715,0067570,0116274,0022404,
2842 0037221,0074157,0053341,0117207,
2843 0036301,0104257,0015075,0004777,
2844 0035164,0057502,0164034,0001313,
2845 0033611,0022254,0176000,0112565,
2846 0031725,0055523,0025153,0166057,
2847 };
2848 static unsigned short APFD[36] = {
2849 /*0040200,0000000,0000000,0000000,*/
2850 0041153,0140334,0130506,0061402,
2851 0041426,0025551,0024440,0070611,
2852 0041373,0134750,0047147,0176702,
2853 0041057,0171532,0105430,0017674,
2854 0040344,0174416,0001726,0047754,
2855 0037421,0021207,0020167,0136264,
2856 0036262,0043621,0151321,0124324,
2857 0034705,0001313,0163733,0016407,
2858 0033027,0166702,0150440,0170561,
2859 };
2860 #endif
2861 #ifdef IBMPC
2862 static unsigned short APFN[36] = {
2863 0x5db5,0x8e7d,0xba0f,0x3fc7,
2864 0xd07e,0x3d14,0x5ff2,0x3fec,
2865 0x01af,0x11be,0x98b7,0x3fef,
2866 0x84a1,0x1397,0xadef,0x3fd9,
2867 0x33d1,0xeadc,0x2f0d,0x3fb2,
2868 0xa140,0xe347,0x3115,0x3f78,
2869 0x8059,0x5d03,0x8be8,0x3f2e,
2870 0x12af,0x9f80,0x2495,0x3ed1,
2871 0x7d86,0x654d,0xab6a,0x3e5a,
2872 };
2873 static unsigned short APFD[36] = {
2874 /*0x0000,0x0000,0x0000,0x3ff0,*/
2875 0xcc60,0x9628,0x781b,0x402d,
2876 0x0e31,0x2524,0xc56d,0x4042,
2877 0xffb8,0x09cc,0x773d,0x403f,
2878 0x03f7,0x5163,0xfe6b,0x4025,
2879 0xc9fd,0xc07a,0x9f21,0x3ffc,
2880 0xf796,0xe40e,0x2450,0x3fc2,
2881 0x351a,0x3a5a,0x48f2,0x3f76,
2882 0x63a1,0x7cfb,0xa059,0x3f18,
2883 0x1e2e,0x5a24,0xfdb8,0x3ea2,
2884 };
2885 #endif
2886 #ifdef MIEEE
2887 static unsigned short APFN[36] = {
2888 0x3fc7,0xba0f,0x8e7d,0x5db5,
2889 0x3fec,0x5ff2,0x3d14,0xd07e,
2890 0x3fef,0x98b7,0x11be,0x01af,
2891 0x3fd9,0xadef,0x1397,0x84a1,
2892 0x3fb2,0x2f0d,0xeadc,0x33d1,
2893 0x3f78,0x3115,0xe347,0xa140,
2894 0x3f2e,0x8be8,0x5d03,0x8059,
2895 0x3ed1,0x2495,0x9f80,0x12af,
2896 0x3e5a,0xab6a,0x654d,0x7d86,
2897 };
2898 static unsigned short APFD[36] = {
2899 /*0x3ff0,0x0000,0x0000,0x0000,*/
2900 0x402d,0x781b,0x9628,0xcc60,
2901 0x4042,0xc56d,0x2524,0x0e31,
2902 0x403f,0x773d,0x09cc,0xffb8,
2903 0x4025,0xfe6b,0x5163,0x03f7,
2904 0x3ffc,0x9f21,0xc07a,0xc9fd,
2905 0x3fc2,0x2450,0xe40e,0xf796,
2906 0x3f76,0x48f2,0x3a5a,0x351a,
2907 0x3f18,0xa059,0x7cfb,0x63a1,
2908 0x3ea2,0xfdb8,0x5a24,0x1e2e,
2909 };
2910 #endif
2911
2912 #ifdef UNK
2913 static double APGN[11] = {
2914 -3.55615429033082288335E-2,
2915 -6.37311518129435504426E-1,
2916 -1.70856738884312371053E0,
2917 -1.50221872117316635393E0,
2918 -5.63606665822102676611E-1,
2919 -1.02101031120216891789E-1,
2920 -9.48396695961445269093E-3,
2921 -4.60325307486780994357E-4,
2922 -1.14300836484517375919E-5,
2923 -1.33415518685547420648E-7,
2924 -5.63803833958893494476E-10,
2925 };
2926 static double APGD[11] = {
2927 /* 1.00000000000000000000E0,*/
2928 9.85865801696130355144E0,
2929 2.16401867356585941885E1,
2930 1.73130776389749389525E1,
2931 6.17872175280828766327E0,
2932 1.08848694396321495475E0,
2933 9.95005543440888479402E-2,
2934 4.78468199683886610842E-3,
2935 1.18159633322838625562E-4,
2936 1.37480673554219441465E-6,
2937 5.79912514929147598821E-9,
2938 };
2939 #endif
2940 #ifdef DEC
2941 static unsigned short APGN[44] = {
2942 0137021,0124372,0176075,0075331,
2943 0140043,0023330,0177672,0161655,
2944 0140332,0131126,0010413,0171112,
2945 0140300,0044263,0175560,0054070,
2946 0140020,0044206,0142603,0073324,
2947 0137321,0015130,0066144,0144033,
2948 0136433,0061243,0175542,0103373,
2949 0135361,0053721,0020441,0053203,
2950 0134077,0141725,0160277,0130612,
2951 0132417,0040372,0100363,0060200,
2952 0130432,0175052,0171064,0034147,
2953 };
2954 static unsigned short APGD[40] = {
2955 /*0040200,0000000,0000000,0000000,*/
2956 0041035,0136420,0030124,0140220,
2957 0041255,0017432,0034447,0162256,
2958 0041212,0100456,0154544,0006321,
2959 0040705,0134026,0127154,0123414,
2960 0040213,0051612,0044470,0172607,
2961 0037313,0143362,0053273,0157051,
2962 0036234,0144322,0054536,0007264,
2963 0034767,0146170,0054265,0170342,
2964 0033270,0102777,0167362,0073631,
2965 0031307,0040644,0167103,0021763,
2966 };
2967 #endif
2968 #ifdef IBMPC
2969 static unsigned short APGN[44] = {
2970 0xaf5b,0x5f87,0x351f,0xbfa2,
2971 0x5c76,0x1ff7,0x64db,0xbfe4,
2972 0x7e49,0xc221,0x564a,0xbffb,
2973 0x0b07,0x7f6e,0x0916,0xbff8,
2974 0x6edb,0xd8b0,0x0910,0xbfe2,
2975 0x9903,0x0d8c,0x234b,0xbfba,
2976 0x50df,0x7f6c,0x6c54,0xbf83,
2977 0x2ad0,0x2424,0x2afa,0xbf3e,
2978 0xf631,0xbc17,0xf87a,0xbee7,
2979 0x6c10,0x501e,0xe81f,0xbe81,
2980 0x870d,0x5e46,0x5f45,0xbe03,
2981 };
2982 static unsigned short APGD[40] = {
2983 /*0x0000,0x0000,0x0000,0x3ff0,*/
2984 0x9812,0x060a,0xb7a2,0x4023,
2985 0xfc96,0x4724,0xa3e3,0x4035,
2986 0x819a,0xdb2c,0x5025,0x4031,
2987 0x94e2,0xd5cd,0xb702,0x4018,
2988 0x1eb1,0x4927,0x6a71,0x3ff1,
2989 0x7bc5,0x4ad7,0x78de,0x3fb9,
2990 0xc1d7,0x4b2b,0x991a,0x3f73,
2991 0xbe1c,0x0b16,0xf98f,0x3f1e,
2992 0x4ef3,0xfdde,0x10bf,0x3eb7,
2993 0x647e,0x9dc8,0xe834,0x3e38,
2994 };
2995 #endif
2996 #ifdef MIEEE
2997 static unsigned short APGN[44] = {
2998 0xbfa2,0x351f,0x5f87,0xaf5b,
2999 0xbfe4,0x64db,0x1ff7,0x5c76,
3000 0xbffb,0x564a,0xc221,0x7e49,
3001 0xbff8,0x0916,0x7f6e,0x0b07,
3002 0xbfe2,0x0910,0xd8b0,0x6edb,
3003 0xbfba,0x234b,0x0d8c,0x9903,
3004 0xbf83,0x6c54,0x7f6c,0x50df,
3005 0xbf3e,0x2afa,0x2424,0x2ad0,
3006 0xbee7,0xf87a,0xbc17,0xf631,
3007 0xbe81,0xe81f,0x501e,0x6c10,
3008 0xbe03,0x5f45,0x5e46,0x870d,
3009 };
3010 static unsigned short APGD[40] = {
3011 /*0x3ff0,0x0000,0x0000,0x0000,*/
3012 0x4023,0xb7a2,0x060a,0x9812,
3013 0x4035,0xa3e3,0x4724,0xfc96,
3014 0x4031,0x5025,0xdb2c,0x819a,
3015 0x4018,0xb702,0xd5cd,0x94e2,
3016 0x3ff1,0x6a71,0x4927,0x1eb1,
3017 0x3fb9,0x78de,0x4ad7,0x7bc5,
3018 0x3f73,0x991a,0x4b2b,0xc1d7,
3019 0x3f1e,0xf98f,0x0b16,0xbe1c,
3020 0x3eb7,0x10bf,0xfdde,0x4ef3,
3021 0x3e38,0xe834,0x9dc8,0x647e,
3022 };
3023 #endif
3024
3025 #ifdef ANSIPROT
3026 extern double fabs ( double );
3027 extern double exp ( double );
3028 extern double sqrt ( double );
3029 extern double polevl ( double, void *, int );
3030 extern double p1evl ( double, void *, int );
3031 extern double sin ( double );
3032 extern double cos ( double );
3033 #else
3034 double fabs(), exp(), sqrt();
3035 double polevl(), p1evl(), sin(), cos();
3036 #endif
3037
3038 int
airy(double x,double * ai,double * aip,double * bi,double * bip)3039 airy( double x, double *ai, double *aip, double *bi, double *bip )
3040 {
3041 double z, zz, t, f, g, uf, ug, k, zeta, theta;
3042 int domflg;
3043
3044 domflg = 0;
3045 if( x > MAXAIRY )
3046 {
3047 *ai = 0;
3048 *aip = 0;
3049 *bi = MAXNUM;
3050 *bip = MAXNUM;
3051 return(-1);
3052 }
3053
3054 if( x < -2.09 )
3055 {
3056 domflg = 15;
3057 t = sqrt(-x);
3058 zeta = -2.0 * x * t / 3.0;
3059 t = sqrt(t);
3060 k = sqpii / t;
3061 z = 1.0/zeta;
3062 zz = z * z;
3063 uf = 1.0 + zz * polevl( zz, AFN, 8 ) / p1evl( zz, AFD, 9 );
3064 ug = z * polevl( zz, AGN, 10 ) / p1evl( zz, AGD, 10 );
3065 theta = zeta + 0.25 * PI;
3066 f = sin( theta );
3067 g = cos( theta );
3068 *ai = k * (f * uf - g * ug);
3069 *bi = k * (g * uf + f * ug);
3070 uf = 1.0 + zz * polevl( zz, APFN, 8 ) / p1evl( zz, APFD, 9 );
3071 ug = z * polevl( zz, APGN, 10 ) / p1evl( zz, APGD, 10 );
3072 k = sqpii * t;
3073 *aip = -k * (g * uf + f * ug);
3074 *bip = k * (f * uf - g * ug);
3075 return(0);
3076 }
3077
3078 if( x >= 2.09 ) /* cbrt(9) */
3079 {
3080 domflg = 5;
3081 t = sqrt(x);
3082 zeta = 2.0 * x * t / 3.0;
3083 g = exp( zeta );
3084 t = sqrt(t);
3085 k = 2.0 * t * g;
3086 z = 1.0/zeta;
3087 f = polevl( z, AN, 7 ) / polevl( z, AD, 7 );
3088 *ai = sqpii * f / k;
3089 k = -0.5 * sqpii * t / g;
3090 f = polevl( z, APN, 7 ) / polevl( z, APD, 7 );
3091 *aip = f * k;
3092
3093 if( x > 8.3203353 ) /* zeta > 16 */
3094 {
3095 f = z * polevl( z, BN16, 4 ) / p1evl( z, BD16, 5 );
3096 k = sqpii * g;
3097 *bi = k * (1.0 + f) / t;
3098 f = z * polevl( z, BPPN, 4 ) / p1evl( z, BPPD, 5 );
3099 *bip = k * t * (1.0 + f);
3100 return(0);
3101 }
3102 }
3103
3104 f = 1.0;
3105 g = x;
3106 t = 1.0;
3107 uf = 1.0;
3108 ug = x;
3109 k = 1.0;
3110 z = x * x * x;
3111 while( t > MACHEP )
3112 {
3113 uf *= z;
3114 k += 1.0;
3115 uf /=k;
3116 ug *= z;
3117 k += 1.0;
3118 ug /=k;
3119 uf /=k;
3120 f += uf;
3121 k += 1.0;
3122 ug /=k;
3123 g += ug;
3124 t = fabs(uf/f);
3125 }
3126 uf = c1 * f;
3127 ug = c2 * g;
3128 if( (domflg & 1) == 0 )
3129 *ai = uf - ug;
3130 if( (domflg & 2) == 0 )
3131 *bi = sqrt3 * (uf + ug);
3132
3133 /* the deriviative of ai */
3134 k = 4.0;
3135 uf = x * x/2.0;
3136 ug = z/3.0;
3137 f = uf;
3138 g = 1.0 + ug;
3139 uf /= 3.0;
3140 t = 1.0;
3141
3142 while( t > MACHEP )
3143 {
3144 uf *= z;
3145 ug /=k;
3146 k += 1.0;
3147 ug *= z;
3148 uf /=k;
3149 f += uf;
3150 k += 1.0;
3151 ug /=k;
3152 uf /=k;
3153 g += ug;
3154 k += 1.0;
3155 t = fabs(ug/g);
3156 }
3157
3158 uf = c1 * f;
3159 ug = c2 * g;
3160 if( (domflg & 4) == 0 )
3161 *aip = uf - ug;
3162 if( (domflg & 8) == 0 )
3163 *bip = sqrt3 * (uf + ug);
3164 return(0);
3165 }
3166
3167 void
f_airy(union argument * arg)3168 f_airy(union argument *arg)
3169 {
3170 struct value a;
3171 double x;
3172 double ai, aip, bi, bip;
3173
3174 (void) arg; /* avoid -Wunused warning */
3175 x = real(pop(&a));
3176
3177 airy(x, &ai, &aip, &bi, &bip);
3178
3179 push(Gcomplex(&a, ai, 0.0));
3180 }
3181
3182 #endif /* End choice of low- or high-precision airy function */
3183
3184 /* ** expint.c
3185 *
3186 * DESCRIBE Approximate the exponential integral function
3187 *
3188 *
3189 * /inf -n -zt
3190 * E_n(z) = | t * e dt (n = 0, 1, 2, ...)
3191 * /1
3192 *
3193 *
3194 * CALL p = expint(n, z)
3195 *
3196 * double n >= 0
3197 * double z >= 0
3198 * also: n must be an integer
3199 * either z > 0 or n > 1
3200 *
3201 * WARNING none
3202 *
3203 * RETURN double p > 0
3204 * -1.0 on error condition
3205 *
3206 * REFERENCE Abramowitz and Stegun (1964)
3207 *
3208 * Copyright (c) 2010 James R. Van Zandt, jrvz@comcast.net
3209 */
3210
3211 static double
expint(double n,double z)3212 expint(double n, double z)
3213 {
3214 double y; /* the answer */
3215
3216 {
3217 /* Test for admissibility of arguments */
3218 double junk;
3219 if (n < 0 || z < 0 || modf(n,&junk))
3220 return -1.0;
3221 if (z == 0 && n < 2)
3222 return -1.0;
3223 }
3224
3225 /* special cases */
3226 if (n == 0) return exp(-z)/z;
3227 if (z == 0) return 1/(n-1);
3228
3229 /* for z=3, CF requires 36 terms and series requires 29 */
3230
3231 if (z > 3)
3232 { /* For large z, use continued fraction (Abramowitz & Stegun 5.1.22):
3233 E_n(z) = exp(-z)(1/(z+n/(1+1/(z+(n+1)/1+2/(z+...)))))
3234 The CF is valid and stable for z>0, and efficient for z>1 or so. */
3235 double n0, n1, n2, n3, d0, d1, d2, d3, y_prev=1;
3236 int i;
3237
3238 n0 = 0; n1 = 1;
3239 d0 = 1; d1 = z;
3240
3241 for (i=0; i<333; i++)
3242 { /* evaluate the CF "top down" using the recurrence
3243 relations for the numerators and denominators of
3244 successive convergents */
3245 n2 = n0*(n+i) + n1;
3246 d2 = d0*(n+i) + d1;
3247 n3 = n1*(1+i) + n2*z;
3248 d3 = d1*(1+i) + d2*z;
3249 y = n3/d3;
3250 if (y == y_prev) break;
3251 y_prev = y;
3252 n0 = n2; n1 = n3;
3253 d0 = d2; d1 = d3;
3254
3255 /* Re-scale terms in continued fraction if terms are large */
3256 if (d3 >= OFLOW) {
3257
3258 n0 /= OFLOW;
3259 n1 /= OFLOW;
3260 d0 /= OFLOW;
3261 d1 /= OFLOW;
3262 }
3263 }
3264 y = exp(-z)*y;
3265 }
3266
3267 else
3268 { /* For small z, use series (Abramowitz & Stegun 5.1.12):
3269 E_1(z) = -\gamma + \ln z +
3270 \sum_{m=1}^\infty { (-z)^m \over (m) m! }
3271 The series is valid for z>0, and efficient for z<4 or so. */
3272
3273 /* from Abramowitz & Stegun, Table 1.1 */
3274 double euler_constant = .577215664901532860606512;
3275
3276 double y_prev = 0;
3277 double t, m;
3278
3279 y = -euler_constant - log(z);
3280 t = 1;
3281 for (m = 1; m<333; m++)
3282 {
3283 t = -t*z/m;
3284 y = y - t/m;
3285 if (y == y_prev) break;
3286 y_prev = y;
3287 }
3288
3289 /* For n > 1, use recurrence relation (Abramowitz & Stegun 5.1.14):
3290 n E_{n+1}(z) + z E_n(z) = e^{-z}, n >= 1
3291 The recurrence is unstable for increasing n and z>4 or so,
3292 but okay for z<3. */
3293
3294 for (m=1; m<n; m++)
3295 y=(exp(-z) - z*y)/m;
3296 }
3297 return y;
3298 }
3299
3300 void
f_expint(union argument * arg)3301 f_expint(union argument *arg)
3302 {
3303 struct value a;
3304 double n, x;
3305
3306 (void) arg; /* avoid -Wunused warning */
3307 x = real(pop(&a));
3308 n = real(pop(&a));
3309
3310 x = expint(n, x);
3311 if (x <= -1)
3312 /* Error return from expint --> flag 'undefined' */
3313 undefined = TRUE;
3314
3315 push(Gcomplex(&a, x, 0.0));
3316 }
3317
3318