htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/airy.c

   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 <config.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_sf_trig.h>
  26 #include <gsl/gsl_sf_airy.h>
  27
  28 #include "error.h"
  29 #include "check.h"
  30
  31 #include "chebyshev.h"
  32 #include "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
 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
 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 */
 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
 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
 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
 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
 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 "eval.h"
 848
 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
 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
 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
 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