1 /* specfunc/airy.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 /* Author: G. Jungman */
21
22 #include "gsl__config.h"
23 #include "gsl_math.h"
24 #include "gsl_errno.h"
25 #include "gsl_sf_trig.h"
26 #include "gsl_sf_airy.h"
27
28 #include "gsl_specfunc__error.h"
29 #include "gsl_specfunc__check.h"
30
31 #include "gsl_specfunc__chebyshev.h"
32 #include "gsl_specfunc__cheb_eval_mode.c"
33
34 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
35
36
37 /* chebyshev expansions for Airy modulus and phase
38 based on SLATEC r9aimp()
39
40 Series for AM21 on the interval -1.25000D-01 to 0.
41 with weighted error 2.89E-17
42 log weighted error 16.54
43 significant figures required 14.15
44 decimal places required 17.34
45
46 Series for ATH1 on the interval -1.25000D-01 to 0.
47 with weighted error 2.53E-17
48 log weighted error 16.60
49 significant figures required 15.15
50 decimal places required 17.38
51
52 Series for AM22 on the interval -1.00000D+00 to -1.25000D-01
53 with weighted error 2.99E-17
54 log weighted error 16.52
55 significant figures required 14.57
56 decimal places required 17.28
57
58 Series for ATH2 on the interval -1.00000D+00 to -1.25000D-01
59 with weighted error 2.57E-17
60 log weighted error 16.59
61 significant figures required 15.07
62 decimal places required 17.34
63 */
64
65 static double am21_data[37] = {
66 0.0065809191761485,
67 0.0023675984685722,
68 0.0001324741670371,
69 0.0000157600904043,
70 0.0000027529702663,
71 0.0000006102679017,
72 0.0000001595088468,
73 0.0000000471033947,
74 0.0000000152933871,
75 0.0000000053590722,
76 0.0000000020000910,
77 0.0000000007872292,
78 0.0000000003243103,
79 0.0000000001390106,
80 0.0000000000617011,
81 0.0000000000282491,
82 0.0000000000132979,
83 0.0000000000064188,
84 0.0000000000031697,
85 0.0000000000015981,
86 0.0000000000008213,
87 0.0000000000004296,
88 0.0000000000002284,
89 0.0000000000001232,
90 0.0000000000000675,
91 0.0000000000000374,
92 0.0000000000000210,
93 0.0000000000000119,
94 0.0000000000000068,
95 0.0000000000000039,
96 0.0000000000000023,
97 0.0000000000000013,
98 0.0000000000000008,
99 0.0000000000000005,
100 0.0000000000000003,
101 0.0000000000000001,
102 0.0000000000000001
103 };
104 static cheb_series am21_cs = {
105 am21_data,
106 36,
107 -1, 1,
108 20
109 };
110
111 static double ath1_data[36] = {
112 -0.07125837815669365,
113 -0.00590471979831451,
114 -0.00012114544069499,
115 -0.00000988608542270,
116 -0.00000138084097352,
117 -0.00000026142640172,
118 -0.00000006050432589,
119 -0.00000001618436223,
120 -0.00000000483464911,
121 -0.00000000157655272,
122 -0.00000000055231518,
123 -0.00000000020545441,
124 -0.00000000008043412,
125 -0.00000000003291252,
126 -0.00000000001399875,
127 -0.00000000000616151,
128 -0.00000000000279614,
129 -0.00000000000130428,
130 -0.00000000000062373,
131 -0.00000000000030512,
132 -0.00000000000015239,
133 -0.00000000000007758,
134 -0.00000000000004020,
135 -0.00000000000002117,
136 -0.00000000000001132,
137 -0.00000000000000614,
138 -0.00000000000000337,
139 -0.00000000000000188,
140 -0.00000000000000105,
141 -0.00000000000000060,
142 -0.00000000000000034,
143 -0.00000000000000020,
144 -0.00000000000000011,
145 -0.00000000000000007,
146 -0.00000000000000004,
147 -0.00000000000000002
148 };
149 static cheb_series ath1_cs = {
150 ath1_data,
151 35,
152 -1, 1,
153 15
154 };
155
156 static double am22_data[33] = {
157 -0.01562844480625341,
158 0.00778336445239681,
159 0.00086705777047718,
160 0.00015696627315611,
161 0.00003563962571432,
162 0.00000924598335425,
163 0.00000262110161850,
164 0.00000079188221651,
165 0.00000025104152792,
166 0.00000008265223206,
167 0.00000002805711662,
168 0.00000000976821090,
169 0.00000000347407923,
170 0.00000000125828132,
171 0.00000000046298826,
172 0.00000000017272825,
173 0.00000000006523192,
174 0.00000000002490471,
175 0.00000000000960156,
176 0.00000000000373448,
177 0.00000000000146417,
178 0.00000000000057826,
179 0.00000000000022991,
180 0.00000000000009197,
181 0.00000000000003700,
182 0.00000000000001496,
183 0.00000000000000608,
184 0.00000000000000248,
185 0.00000000000000101,
186 0.00000000000000041,
187 0.00000000000000017,
188 0.00000000000000007,
189 0.00000000000000002
190 };
191 static cheb_series am22_cs = {
192 am22_data,
193 32,
194 -1, 1,
195 15
196 };
197
198 static double ath2_data[32] = {
199 0.00440527345871877,
200 -0.03042919452318455,
201 -0.00138565328377179,
202 -0.00018044439089549,
203 -0.00003380847108327,
204 -0.00000767818353522,
205 -0.00000196783944371,
206 -0.00000054837271158,
207 -0.00000016254615505,
208 -0.00000005053049981,
209 -0.00000001631580701,
210 -0.00000000543420411,
211 -0.00000000185739855,
212 -0.00000000064895120,
213 -0.00000000023105948,
214 -0.00000000008363282,
215 -0.00000000003071196,
216 -0.00000000001142367,
217 -0.00000000000429811,
218 -0.00000000000163389,
219 -0.00000000000062693,
220 -0.00000000000024260,
221 -0.00000000000009461,
222 -0.00000000000003716,
223 -0.00000000000001469,
224 -0.00000000000000584,
225 -0.00000000000000233,
226 -0.00000000000000093,
227 -0.00000000000000037,
228 -0.00000000000000015,
229 -0.00000000000000006,
230 -0.00000000000000002
231 };
232 static cheb_series ath2_cs = {
233 ath2_data,
234 31,
235 -1, 1,
236 16
237 };
238
239
240 /* Airy modulus and phase for x < -1 */
241 static
242 int
airy_mod_phase(const double x,gsl_mode_t mode,gsl_sf_result * mod,gsl_sf_result * phase)243 airy_mod_phase(const double x, gsl_mode_t mode, gsl_sf_result * mod, gsl_sf_result * phase)
244 {
245 gsl_sf_result result_m;
246 gsl_sf_result result_p;
247 double m, p;
248 double sqx;
249
250 if(x < -2.0) {
251 double z = 16.0/(x*x*x) + 1.0;
252 cheb_eval_mode_e(&am21_cs, z, mode, &result_m);
253 cheb_eval_mode_e(&ath1_cs, z, mode, &result_p);
254 }
255 else if(x <= -1.0) {
256 double z = (16.0/(x*x*x) + 9.0)/7.0;
257 cheb_eval_mode_e(&am22_cs, z, mode, &result_m);
258 cheb_eval_mode_e(&ath2_cs, z, mode, &result_p);
259 }
260 else {
261 mod->val = 0.0;
262 mod->err = 0.0;
263 phase->val = 0.0;
264 phase->err = 0.0;
265 GSL_ERROR ("x is greater than 1.0", GSL_EDOM);
266 }
267
268 m = 0.3125 + result_m.val;
269 p = -0.625 + result_p.val;
270
271 sqx = sqrt(-x);
272
273 mod->val = sqrt(m/sqx);
274 mod->err = fabs(mod->val) * (GSL_DBL_EPSILON + fabs(result_m.err/result_m.val));
275 phase->val = M_PI_4 - x*sqx * p;
276 phase->err = fabs(phase->val) * (GSL_DBL_EPSILON + fabs(result_p.err/result_p.val));
277
278 return GSL_SUCCESS;
279 }
280
281
282
283 /* Chebyshev series for Ai(x) with x in [-1,1]
284 based on SLATEC ai(x)
285
286 series for aif on the interval -1.00000d+00 to 1.00000d+00
287 with weighted error 1.09e-19
288 log weighted error 18.96
289 significant figures required 17.76
290 decimal places required 19.44
291
292 series for aig on the interval -1.00000d+00 to 1.00000d+00
293 with weighted error 1.51e-17
294 log weighted error 16.82
295 significant figures required 15.19
296 decimal places required 17.27
297 */
298 static double ai_data_f[9] = {
299 -0.03797135849666999750,
300 0.05919188853726363857,
301 0.00098629280577279975,
302 0.00000684884381907656,
303 0.00000002594202596219,
304 0.00000000006176612774,
305 0.00000000000010092454,
306 0.00000000000000012014,
307 0.00000000000000000010
308 };
309 static cheb_series aif_cs = {
310 ai_data_f,
311 8,
312 -1, 1,
313 8
314 };
315
316 static double ai_data_g[8] = {
317 0.01815236558116127,
318 0.02157256316601076,
319 0.00025678356987483,
320 0.00000142652141197,
321 0.00000000457211492,
322 0.00000000000952517,
323 0.00000000000001392,
324 0.00000000000000001
325 };
326 static cheb_series aig_cs = {
327 ai_data_g,
328 7,
329 -1, 1,
330 7
331 };
332
333
334 /* Chebvyshev series for Bi(x) with x in [-1,1]
335 based on SLATEC bi(x)
336
337 series for bif on the interval -1.00000d+00 to 1.00000d+00
338 with weighted error 1.88e-19
339 log weighted error 18.72
340 significant figures required 17.74
341 decimal places required 19.20
342
343 series for big on the interval -1.00000d+00 to 1.00000d+00
344 with weighted error 2.61e-17
345 log weighted error 16.58
346 significant figures required 15.17
347 decimal places required 17.03
348 */
349 static double data_bif[9] = {
350 -0.01673021647198664948,
351 0.10252335834249445610,
352 0.00170830925073815165,
353 0.00001186254546774468,
354 0.00000004493290701779,
355 0.00000000010698207143,
356 0.00000000000017480643,
357 0.00000000000000020810,
358 0.00000000000000000018
359 };
360 static cheb_series bif_cs = {
361 data_bif,
362 8,
363 -1, 1,
364 8
365 };
366
367 static double data_big[8] = {
368 0.02246622324857452,
369 0.03736477545301955,
370 0.00044476218957212,
371 0.00000247080756363,
372 0.00000000791913533,
373 0.00000000001649807,
374 0.00000000000002411,
375 0.00000000000000002
376 };
377 static cheb_series big_cs = {
378 data_big,
379 7,
380 -1, 1,
381 7
382 };
383
384
385 /* Chebyshev series for Bi(x) with x in [1,8]
386 based on SLATEC bi(x)
387 */
388 static double data_bif2[10] = {
389 0.0998457269381604100,
390 0.4786249778630055380,
391 0.0251552119604330118,
392 0.0005820693885232645,
393 0.0000074997659644377,
394 0.0000000613460287034,
395 0.0000000003462753885,
396 0.0000000000014288910,
397 0.0000000000000044962,
398 0.0000000000000000111
399 };
400 static cheb_series bif2_cs = {
401 data_bif2,
402 9,
403 -1, 1,
404 9
405 };
406
407 static double data_big2[10] = {
408 0.033305662145514340,
409 0.161309215123197068,
410 0.0063190073096134286,
411 0.0001187904568162517,
412 0.0000013045345886200,
413 0.0000000093741259955,
414 0.0000000000474580188,
415 0.0000000000001783107,
416 0.0000000000000005167,
417 0.0000000000000000011
418 };
419 static cheb_series big2_cs = {
420 data_big2,
421 9,
422 -1, 1,
423 9
424 };
425
426
427 /* chebyshev for Ai(x) asymptotic factor
428 based on SLATEC aie()
429
430 Series for AIP on the interval 0. to 1.00000D+00
431 with weighted error 5.10E-17
432 log weighted error 16.29
433 significant figures required 14.41
434 decimal places required 17.06
435
436 [GJ] Sun Apr 19 18:14:31 EDT 1998
437 There was something wrong with these coefficients. I was getting
438 errors after 3 or 4 digits. So I recomputed this table. Now I get
439 double precision agreement with Mathematica. But it does not seem
440 possible that the small differences here would account for the
441 original discrepancy. There must have been something wrong with my
442 original usage...
443 */
444 static double data_aip[36] = {
445 -0.0187519297793867540198,
446 -0.0091443848250055004725,
447 0.0009010457337825074652,
448 -0.0001394184127221491507,
449 0.0000273815815785209370,
450 -0.0000062750421119959424,
451 0.0000016064844184831521,
452 -0.0000004476392158510354,
453 0.0000001334635874651668,
454 -0.0000000420735334263215,
455 0.0000000139021990246364,
456 -0.0000000047831848068048,
457 0.0000000017047897907465,
458 -0.0000000006268389576018,
459 0.0000000002369824276612,
460 -0.0000000000918641139267,
461 0.0000000000364278543037,
462 -0.0000000000147475551725,
463 0.0000000000060851006556,
464 -0.0000000000025552772234,
465 0.0000000000010906187250,
466 -0.0000000000004725870319,
467 0.0000000000002076969064,
468 -0.0000000000000924976214,
469 0.0000000000000417096723,
470 -0.0000000000000190299093,
471 0.0000000000000087790676,
472 -0.0000000000000040927557,
473 0.0000000000000019271068,
474 -0.0000000000000009160199,
475 0.0000000000000004393567,
476 -0.0000000000000002125503,
477 0.0000000000000001036735,
478 -0.0000000000000000509642,
479 0.0000000000000000252377,
480 -0.0000000000000000125793
481 /*
482 -.0187519297793868
483 -.0091443848250055,
484 .0009010457337825,
485 -.0001394184127221,
486 .0000273815815785,
487 -.0000062750421119,
488 .0000016064844184,
489 -.0000004476392158,
490 .0000001334635874,
491 -.0000000420735334,
492 .0000000139021990,
493 -.0000000047831848,
494 .0000000017047897,
495 -.0000000006268389,
496 .0000000002369824,
497 -.0000000000918641,
498 .0000000000364278,
499 -.0000000000147475,
500 .0000000000060851,
501 -.0000000000025552,
502 .0000000000010906,
503 -.0000000000004725,
504 .0000000000002076,
505 -.0000000000000924,
506 .0000000000000417,
507 -.0000000000000190,
508 .0000000000000087,
509 -.0000000000000040,
510 .0000000000000019,
511 -.0000000000000009,
512 .0000000000000004,
513 -.0000000000000002,
514 .0000000000000001,
515 -.0000000000000000
516 */
517 };
518 static cheb_series aip_cs = {
519 data_aip,
520 35,
521 -1, 1,
522 17
523 };
524
525
526 /* chebyshev for Bi(x) asymptotic factor
527 based on SLATEC bie()
528
529 Series for BIP on the interval 1.25000D-01 to 3.53553D-01
530 with weighted error 1.91E-17
531 log weighted error 16.72
532 significant figures required 15.35
533 decimal places required 17.41
534
535 Series for BIP2 on the interval 0. to 1.25000D-01
536 with weighted error 1.05E-18
537 log weighted error 17.98
538 significant figures required 16.74
539 decimal places required 18.71
540 */
541 static double data_bip[24] = {
542 -0.08322047477943447,
543 0.01146118927371174,
544 0.00042896440718911,
545 -0.00014906639379950,
546 -0.00001307659726787,
547 0.00000632759839610,
548 -0.00000042226696982,
549 -0.00000019147186298,
550 0.00000006453106284,
551 -0.00000000784485467,
552 -0.00000000096077216,
553 0.00000000070004713,
554 -0.00000000017731789,
555 0.00000000002272089,
556 0.00000000000165404,
557 -0.00000000000185171,
558 0.00000000000059576,
559 -0.00000000000012194,
560 0.00000000000001334,
561 0.00000000000000172,
562 -0.00000000000000145,
563 0.00000000000000049,
564 -0.00000000000000011,
565 0.00000000000000001
566 };
567 static cheb_series bip_cs = {
568 data_bip,
569 23,
570 -1, 1,
571 14
572 };
573
574 static double data_bip2[29] = {
575 -0.113596737585988679,
576 0.0041381473947881595,
577 0.0001353470622119332,
578 0.0000104273166530153,
579 0.0000013474954767849,
580 0.0000001696537405438,
581 -0.0000000100965008656,
582 -0.0000000167291194937,
583 -0.0000000045815364485,
584 0.0000000003736681366,
585 0.0000000005766930320,
586 0.0000000000621812650,
587 -0.0000000000632941202,
588 -0.0000000000149150479,
589 0.0000000000078896213,
590 0.0000000000024960513,
591 -0.0000000000012130075,
592 -0.0000000000003740493,
593 0.0000000000002237727,
594 0.0000000000000474902,
595 -0.0000000000000452616,
596 -0.0000000000000030172,
597 0.0000000000000091058,
598 -0.0000000000000009814,
599 -0.0000000000000016429,
600 0.0000000000000005533,
601 0.0000000000000002175,
602 -0.0000000000000001737,
603 -0.0000000000000000010
604 };
605 static cheb_series bip2_cs = {
606 data_bip2,
607 28,
608 -1, 1,
609 10
610 };
611
612
613 /* assumes x >= 1.0 */
614 inline static int
airy_aie(const double x,gsl_mode_t mode,gsl_sf_result * result)615 airy_aie(const double x, gsl_mode_t mode, gsl_sf_result * result)
616 {
617 double sqx = sqrt(x);
618 double z = 2.0/(x*sqx) - 1.0;
619 double y = sqrt(sqx);
620 gsl_sf_result result_c;
621 cheb_eval_mode_e(&aip_cs, z, mode, &result_c);
622 result->val = (0.28125 + result_c.val)/y;
623 result->err = result_c.err/y + GSL_DBL_EPSILON * fabs(result->val);
624 return GSL_SUCCESS;
625 }
626
627 /* assumes x >= 2.0 */
airy_bie(const double x,gsl_mode_t mode,gsl_sf_result * result)628 static int airy_bie(const double x, gsl_mode_t mode, gsl_sf_result * result)
629 {
630 const double ATR = 8.7506905708484345;
631 const double BTR = -2.0938363213560543;
632
633 if(x < 4.0) {
634 double sqx = sqrt(x);
635 double z = ATR/(x*sqx) + BTR;
636 double y = sqrt(sqx);
637 gsl_sf_result result_c;
638 cheb_eval_mode_e(&bip_cs, z, mode, &result_c);
639 result->val = (0.625 + result_c.val)/y;
640 result->err = result_c.err/y + GSL_DBL_EPSILON * fabs(result->val);
641 }
642 else {
643 double sqx = sqrt(x);
644 double z = 16.0/(x*sqx) - 1.0;
645 double y = sqrt(sqx);
646 gsl_sf_result result_c;
647 cheb_eval_mode_e(&bip2_cs, z, mode, &result_c);
648 result->val = (0.625 + result_c.val)/y;
649 result->err = result_c.err/y + GSL_DBL_EPSILON * fabs(result->val);
650 }
651
652 return GSL_SUCCESS;
653 }
654
655
656 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
657
658 int
gsl_sf_airy_Ai_e(const double x,const gsl_mode_t mode,gsl_sf_result * result)659 gsl_sf_airy_Ai_e(const double x, const gsl_mode_t mode, gsl_sf_result * result)
660 {
661 /* CHECK_POINTER(result) */
662
663 if(x < -1.0) {
664 gsl_sf_result mod;
665 gsl_sf_result theta;
666 gsl_sf_result cos_result;
667 int stat_mp = airy_mod_phase(x, mode, &mod, &theta);
668 int stat_cos = gsl_sf_cos_err_e(theta.val, theta.err, &cos_result);
669 result->val = mod.val * cos_result.val;
670 result->err = fabs(mod.val * cos_result.err) + fabs(cos_result.val * mod.err);
671 result->err += GSL_DBL_EPSILON * fabs(result->val);
672 return GSL_ERROR_SELECT_2(stat_mp, stat_cos);
673 }
674 else if(x <= 1.0) {
675 const double z = x*x*x;
676 gsl_sf_result result_c0;
677 gsl_sf_result result_c1;
678 cheb_eval_mode_e(&aif_cs, z, mode, &result_c0);
679 cheb_eval_mode_e(&aig_cs, z, mode, &result_c1);
680 result->val = 0.375 + (result_c0.val - x*(0.25 + result_c1.val));
681 result->err = result_c0.err + fabs(x*result_c1.err);
682 result->err += GSL_DBL_EPSILON * fabs(result->val);
683 return GSL_SUCCESS;
684 }
685 else {
686 double x32 = x * sqrt(x);
687 double s = exp(-2.0*x32/3.0);
688 gsl_sf_result result_aie;
689 int stat_aie = airy_aie(x, mode, &result_aie);
690 result->val = result_aie.val * s;
691 result->err = result_aie.err * s + result->val * x32 * GSL_DBL_EPSILON;
692 result->err += GSL_DBL_EPSILON * fabs(result->val);
693 CHECK_UNDERFLOW(result);
694 return stat_aie;
695 }
696 }
697
698
699 int
gsl_sf_airy_Ai_scaled_e(const double x,gsl_mode_t mode,gsl_sf_result * result)700 gsl_sf_airy_Ai_scaled_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
701 {
702 /* CHECK_POINTER(result) */
703
704 if(x < -1.0) {
705 gsl_sf_result mod;
706 gsl_sf_result theta;
707 gsl_sf_result cos_result;
708 int stat_mp = airy_mod_phase(x, mode, &mod, &theta);
709 int stat_cos = gsl_sf_cos_err_e(theta.val, theta.err, &cos_result);
710 result->val = mod.val * cos_result.val;
711 result->err = fabs(mod.val * cos_result.err) + fabs(cos_result.val * mod.err);
712 result->err += GSL_DBL_EPSILON * fabs(result->val);
713 return GSL_ERROR_SELECT_2(stat_mp, stat_cos);
714 }
715 else if(x <= 1.0) {
716 const double z = x*x*x;
717 gsl_sf_result result_c0;
718 gsl_sf_result result_c1;
719 cheb_eval_mode_e(&aif_cs, z, mode, &result_c0);
720 cheb_eval_mode_e(&aig_cs, z, mode, &result_c1);
721 result->val = 0.375 + (result_c0.val - x*(0.25 + result_c1.val));
722 result->err = result_c0.err + fabs(x*result_c1.err);
723 result->err += GSL_DBL_EPSILON * fabs(result->val);
724
725 if(x > 0.0) {
726 const double scale = exp(2.0/3.0 * sqrt(z));
727 result->val *= scale;
728 result->err *= scale;
729 }
730
731 return GSL_SUCCESS;
732 }
733 else {
734 return airy_aie(x, mode, result);
735 }
736 }
737
738
gsl_sf_airy_Bi_e(const double x,gsl_mode_t mode,gsl_sf_result * result)739 int gsl_sf_airy_Bi_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
740 {
741 /* CHECK_POINTER(result) */
742 if(x < -1.0) {
743 gsl_sf_result mod;
744 gsl_sf_result theta;
745 gsl_sf_result sin_result;
746 int stat_mp = airy_mod_phase(x, mode, &mod, &theta);
747 int stat_sin = gsl_sf_sin_err_e(theta.val, theta.err, &sin_result);
748 result->val = mod.val * sin_result.val;
749 result->err = fabs(mod.val * sin_result.err) + fabs(sin_result.val * mod.err);
750 result->err += GSL_DBL_EPSILON * fabs(result->val);
751 return GSL_ERROR_SELECT_2(stat_mp, stat_sin);
752 }
753 else if(x < 1.0) {
754 const double z = x*x*x;
755 gsl_sf_result result_c0;
756 gsl_sf_result result_c1;
757 cheb_eval_mode_e(&bif_cs, z, mode, &result_c0);
758 cheb_eval_mode_e(&big_cs, z, mode, &result_c1);
759 result->val = 0.625 + result_c0.val + x*(0.4375 + result_c1.val);
760 result->err = result_c0.err + fabs(x * result_c1.err);
761 result->err += GSL_DBL_EPSILON * fabs(result->val);
762 return GSL_SUCCESS;
763 }
764 else if(x <= 2.0) {
765 const double z = (2.0*x*x*x - 9.0)/7.0;
766 gsl_sf_result result_c0;
767 gsl_sf_result result_c1;
768 cheb_eval_mode_e(&bif2_cs, z, mode, &result_c0);
769 cheb_eval_mode_e(&big2_cs, z, mode, &result_c1);
770 result->val = 1.125 + result_c0.val + x*(0.625 + result_c1.val);
771 result->err = result_c0.err + fabs(x * result_c1.err);
772 result->err += GSL_DBL_EPSILON * fabs(result->val);
773 return GSL_SUCCESS;
774 }
775 else {
776 const double y = 2.0*x*sqrt(x)/3.0;
777 const double s = exp(y);
778
779 if(y > GSL_LOG_DBL_MAX - 1.0) {
780 OVERFLOW_ERROR(result);
781 }
782 else {
783 gsl_sf_result result_bie;
784 int stat_bie = airy_bie(x, mode, &result_bie);
785 result->val = result_bie.val * s;
786 result->err = result_bie.err * s + fabs(1.5*y * (GSL_DBL_EPSILON * result->val));
787 result->err += GSL_DBL_EPSILON * fabs(result->val);
788 return stat_bie;
789 }
790 }
791 }
792
793
794 int
gsl_sf_airy_Bi_scaled_e(const double x,gsl_mode_t mode,gsl_sf_result * result)795 gsl_sf_airy_Bi_scaled_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
796 {
797 /* CHECK_POINTER(result) */
798
799 if(x < -1.0) {
800 gsl_sf_result mod;
801 gsl_sf_result theta;
802 gsl_sf_result sin_result;
803 int stat_mp = airy_mod_phase(x, mode, &mod, &theta);
804 int stat_sin = gsl_sf_sin_err_e(theta.val, theta.err, &sin_result);
805 result->val = mod.val * sin_result.val;
806 result->err = fabs(mod.val * sin_result.err) + fabs(sin_result.val * mod.err);
807 result->err += GSL_DBL_EPSILON * fabs(result->val);
808 return GSL_ERROR_SELECT_2(stat_mp, stat_sin);
809 }
810 else if(x < 1.0) {
811 const double z = x*x*x;
812 gsl_sf_result result_c0;
813 gsl_sf_result result_c1;
814 cheb_eval_mode_e(&bif_cs, z, mode, &result_c0);
815 cheb_eval_mode_e(&big_cs, z, mode, &result_c1);
816 result->val = 0.625 + result_c0.val + x*(0.4375 + result_c1.val);
817 result->err = result_c0.err + fabs(x * result_c1.err);
818 result->err += GSL_DBL_EPSILON * fabs(result->val);
819 if(x > 0.0) {
820 const double scale = exp(-2.0/3.0 * sqrt(z));
821 result->val *= scale;
822 result->err *= scale;
823 }
824 return GSL_SUCCESS;
825 }
826 else if(x <= 2.0) {
827 const double x3 = x*x*x;
828 const double z = (2.0*x3 - 9.0)/7.0;
829 const double s = exp(-2.0/3.0 * sqrt(x3));
830 gsl_sf_result result_c0;
831 gsl_sf_result result_c1;
832 cheb_eval_mode_e(&bif2_cs, z, mode, &result_c0);
833 cheb_eval_mode_e(&big2_cs, z, mode, &result_c1);
834 result->val = s * (1.125 + result_c0.val + x*(0.625 + result_c1.val));
835 result->err = s * (result_c0.err + fabs(x * result_c1.err));
836 result->err += GSL_DBL_EPSILON * fabs(result->val);
837 return GSL_SUCCESS;
838 }
839 else {
840 return airy_bie(x, mode, result);
841 }
842 }
843
844
845 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
846
847 #include "gsl_specfunc__eval.h"
848
gsl_sf_airy_Ai(const double x,gsl_mode_t mode)849 double gsl_sf_airy_Ai(const double x, gsl_mode_t mode)
850 {
851 EVAL_RESULT(gsl_sf_airy_Ai_e(x, mode, &result));
852 }
853
gsl_sf_airy_Ai_scaled(const double x,gsl_mode_t mode)854 double gsl_sf_airy_Ai_scaled(const double x, gsl_mode_t mode)
855 {
856 EVAL_RESULT(gsl_sf_airy_Ai_scaled_e(x, mode, &result));
857 }
858
gsl_sf_airy_Bi(const double x,gsl_mode_t mode)859 double gsl_sf_airy_Bi(const double x, gsl_mode_t mode)
860 {
861 EVAL_RESULT(gsl_sf_airy_Bi_e(x, mode, &result));
862 }
863
gsl_sf_airy_Bi_scaled(const double x,gsl_mode_t mode)864 double gsl_sf_airy_Bi_scaled(const double x, gsl_mode_t mode)
865 {
866 EVAL_RESULT(gsl_sf_airy_Bi_scaled_e(x, mode, &result));
867 }
868
869
870