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