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

   1 /* specfunc/airy_der.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_exp.h>
  26 #include <gsl/gsl_sf_airy.h>
  27
  28 #include "error.h"
  29
  30 #include "chebyshev.h"
  31 #include "cheb_eval_mode.c"
  32
  33 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  34
  35
  36 /* based on SLATEC aide.f, bide.f, aid.f, bid.f, r9admp.f */
  37
  38 /*
  39  series for aif on the interval -1.00000e+00 to  1.00000e+00
  40                                         with weighted error   5.22e-18
  41                                          log weighted error  17.28
  42                                significant figures required  16.01
  43                                     decimal places required  17.73
  44 */
  45 static double aif_data[8] = {
  46    0.10527461226531408809,
  47    0.01183613628152997844,
  48    0.00012328104173225664,
  49    0.00000062261225638140,
  50    0.00000000185298887844,
  51    0.00000000000363328873,
  52    0.00000000000000504622,
  53    0.00000000000000000522
  54 };
  55 static cheb_series aif_cs = {
  56   aif_data,
  57   7,
  58   -1, 1,
  59   7
  60 };
  61
  62 /*
  63  series for aig on the interval -1.00000e+00 to  1.00000e+00
  64                                         with weighted error   3.14e-19
  65                                          log weighted error  18.50
  66                                significant figures required  17.44
  67                                     decimal places required  18.98
  68 */
  69 static double aig_data[9] = {
  70    0.021233878150918666852,
  71    0.086315930335214406752,
  72    0.001797594720383231358,
  73    0.000014265499875550693,
  74    0.000000059437995283683,
  75    0.000000000152403366479,
  76    0.000000000000264587660,
  77    0.000000000000000331562,
  78    0.000000000000000000314
  79 };
  80 static cheb_series aig_cs = {
  81   aig_data,
  82   8,
  83   -1, 1,
  84   8
  85 };
  86
  87 /*
  88  series for aip2 on the interval  0.00000e+00 to  1.25000e-01
  89                                         with weighted error   2.15e-17
  90                                          log weighted error  16.67
  91                                significant figures required  14.27
  92                                     decimal places required  17.26
  93 */
  94 static double aip2_data[15] = {
  95     0.0065457691989713757,
  96     0.0023833724120774592,
  97    -0.0000430700770220586,
  98     0.0000015629125858629,
  99    -0.0000000815417186163,
 100     0.0000000054103738057,
 101    -0.0000000004284130883,
 102     0.0000000000389497963,
 103    -0.0000000000039623161,
 104     0.0000000000004428184,
 105    -0.0000000000000536297,
 106     0.0000000000000069650,
 107    -0.0000000000000009620,
 108     0.0000000000000001403,
 109    -0.0000000000000000215
 110 };
 111 static cheb_series aip2_cs = {
 112   aip2_data,
 113   14,
 114   -1, 1,
 115   9
 116 };
 117
 118 /*
 119  series for aip1 on the interval  1.25000e-01 to  1.00000e+00
 120                                         with weighted error   2.60e-17
 121                                          log weighted error  16.58
 122                                significant figures required  14.91
 123                                     decimal places required  17.28
 124 */
 125 static double aip1_data[25] = {
 126     0.0358865097808301538,
 127     0.0114668575627764899,
 128    -0.0007592073583861400,
 129     0.0000869517610893841,
 130    -0.0000128237294298592,
 131     0.0000022062695681038,
 132    -0.0000004222295185921,
 133     0.0000000874686415726,
 134    -0.0000000192773588418,
 135     0.0000000044668460054,
 136    -0.0000000010790108052,
 137     0.0000000002700029447,
 138    -0.0000000000696480108,
 139     0.0000000000184489907,
 140    -0.0000000000050027817,
 141     0.0000000000013852243,
 142    -0.0000000000003908218,
 143     0.0000000000001121536,
 144    -0.0000000000000326862,
 145     0.0000000000000096619,
 146    -0.0000000000000028935,
 147     0.0000000000000008770,
 148    -0.0000000000000002688,
 149     0.0000000000000000832,
 150    -0.0000000000000000260
 151 };
 152 static cheb_series aip1_cs = {
 153   aip1_data,
 154   24,
 155   -1, 1,
 156   14
 157 };
 158
 159
 160 /*
 161  series for bif on the interval -1.00000e+00 to  1.00000e+00
 162                                         with weighted error   9.05e-18
 163                                          log weighted error  17.04
 164                                significant figures required  15.83
 165                                     decimal places required  17.49
 166 */
 167 static double bif_data[8] = {
 168    0.1153536790828570243,
 169    0.0205007894049192875,
 170    0.0002135290278902876,
 171    0.0000010783960614677,
 172    0.0000000032094708833,
 173    0.0000000000062930407,
 174    0.0000000000000087403,
 175    0.0000000000000000090
 176 };
 177 static cheb_series bif_cs = {
 178   bif_data,
 179   7,
 180   -1, 1,
 181   7
 182 };
 183
 184 /*
 185  series for big on the interval -1.00000e+00 to  1.00000e+00
 186                                         with weighted error   5.44e-19
 187                                          log weighted error  18.26
 188                                significant figures required  17.46
 189                                     decimal places required  18.74
 190 */
 191 static double big_data[9] = {
 192    -0.097196440416443537390,
 193     0.149503576843167066571,
 194     0.003113525387121326042,
 195     0.000024708570579821297,
 196     0.000000102949627731379,
 197     0.000000000263970373987,
 198     0.000000000000458279271,
 199     0.000000000000000574283,
 200     0.000000000000000000544
 201 };
 202 static cheb_series big_cs = {
 203   big_data,
 204   8,
 205   -1, 1,
 206   8
 207 };
 208
 209 /*
 210  series for bif2 on the interval  1.00000e+00 to  8.00000e+00
 211                                         with weighted error   3.82e-19
 212                                          log weighted error  18.42
 213                                significant figures required  17.68
 214                                     decimal places required  18.92
 215 */
 216 static double bif2_data[10] = {
 217    0.323493987603522033521,
 218    0.086297871535563559139,
 219    0.002994025552655397426,
 220    0.000051430528364661637,
 221    0.000000525840250036811,
 222    0.000000003561751373958,
 223    0.000000000017146864007,
 224    0.000000000000061663520,
 225    0.000000000000000171911,
 226    0.000000000000000000382
 227 };
 228 static cheb_series bif2_cs = {
 229   bif2_data,
 230   9,
 231   -1, 1,
 232   9
 233 };
 234
 235 /*
 236  series for big2 on the interval  1.00000e+00 to  8.00000e+00
 237                                         with weighted error   3.35e-17
 238                                          log weighted error  16.48
 239                                significant figures required  16.52
 240                                     decimal places required  16.98
 241 */
 242 static double big2_data[10] = {
 243    1.6062999463621294578,
 244    0.7449088819876088652,
 245    0.0470138738610277380,
 246    0.0012284422062548239,
 247    0.0000173222412256624,
 248    0.0000001521901652368,
 249    0.0000000009113560249,
 250    0.0000000000039547918,
 251    0.0000000000000130017,
 252    0.0000000000000000335
 253 };
 254 static cheb_series big2_cs = {
 255   big2_data,
 256   9,
 257   -1, 1,
 258   9
 259 };
 260
 261 /*
 262  series for bip2 on the interval  0.00000e+00 to  1.25000e-01
 263                                         with weighted error   2.07e-18
 264                                          log weighted error  17.69
 265                                significant figures required  16.51
 266                                     decimal places required  18.42
 267 */
 268 static double bip2_data[29] = {
 269     -0.13269705443526630495,
 270     -0.00568443626045977481,
 271     -0.00015643601119611610,
 272     -0.00001136737203679562,
 273     -0.00000143464350991284,
 274     -0.00000018098531185164,
 275      0.00000000926177343611,
 276      0.00000001710005490721,
 277      0.00000000476698163504,
 278     -0.00000000035195022023,
 279     -0.00000000058890614316,
 280     -0.00000000006678499608,
 281      0.00000000006395565102,
 282      0.00000000001554529427,
 283     -0.00000000000792397000,
 284     -0.00000000000258326243,
 285      0.00000000000121655048,
 286      0.00000000000038707207,
 287     -0.00000000000022487045,
 288     -0.00000000000004953477,
 289      0.00000000000004563782,
 290      0.00000000000000332998,
 291     -0.00000000000000921750,
 292      0.00000000000000094157,
 293      0.00000000000000167154,
 294     -0.00000000000000055134,
 295     -0.00000000000000022369,
 296      0.00000000000000017487,
 297      0.00000000000000000207
 298 };
 299 static cheb_series bip2_cs = {
 300   bip2_data,
 301   28,
 302   -1, 1,
 303   14
 304 };
 305
 306 /*
 307  series for bip1 on the interval  1.25000e-01 to  3.53553e-01
 308                                         with weighted error   1.86e-17
 309                                          log weighted error  16.73
 310                                significant figures required  15.67
 311                                     decimal places required  17.42
 312 */
 313 static double bip1_data[24] = {
 314    -0.1729187351079553719,
 315    -0.0149358492984694364,
 316    -0.0005471104951678566,
 317     0.0001537966292958408,
 318     0.0000154353476192179,
 319    -0.0000065434113851906,
 320     0.0000003728082407879,
 321     0.0000002072078388189,
 322    -0.0000000658173336470,
 323     0.0000000074926746354,
 324     0.0000000011101336884,
 325    -0.0000000007265140553,
 326     0.0000000001782723560,
 327    -0.0000000000217346352,
 328    -0.0000000000020302035,
 329     0.0000000000019311827,
 330    -0.0000000000006044953,
 331     0.0000000000001209450,
 332    -0.0000000000000125109,
 333    -0.0000000000000019917,
 334     0.0000000000000015154,
 335    -0.0000000000000004977,
 336     0.0000000000000001155,
 337    -0.0000000000000000186
 338 };
 339 static cheb_series bip1_cs = {
 340   bip1_data,
 341   23,
 342   -1, 1,
 343   13
 344 };
 345
 346 /*
 347  series for an22 on the interval -1.00000e+00 to -1.25000e-01
 348                                         with weighted error   3.30e-17
 349                                          log weighted error  16.48
 350                                significant figures required  14.95
 351                                     decimal places required  17.24
 352 */
 353 static double an22_data[33] = {
 354     0.0537418629629794329,
 355    -0.0126661435859883193,
 356    -0.0011924334106593007,
 357    -0.0002032327627275655,
 358    -0.0000446468963075164,
 359    -0.0000113359036053123,
 360    -0.0000031641352378546,
 361    -0.0000009446708886149,
 362    -0.0000002966562236472,
 363    -0.0000000969118892024,
 364    -0.0000000326822538653,
 365    -0.0000000113144618964,
 366    -0.0000000040042691002,
 367    -0.0000000014440333684,
 368    -0.0000000005292853746,
 369    -0.0000000001967763374,
 370    -0.0000000000740800096,
 371    -0.0000000000282016314,
 372    -0.0000000000108440066,
 373    -0.0000000000042074801,
 374    -0.0000000000016459150,
 375    -0.0000000000006486827,
 376    -0.0000000000002574095,
 377    -0.0000000000001027889,
 378    -0.0000000000000412846,
 379    -0.0000000000000166711,
 380    -0.0000000000000067657,
 381    -0.0000000000000027585,
 382    -0.0000000000000011296,
 383    -0.0000000000000004645,
 384    -0.0000000000000001917,
 385    -0.0000000000000000794,
 386    -0.0000000000000000330
 387 };
 388 static cheb_series an22_cs = {
 389   an22_data,
 390   32,
 391   -1, 1,
 392   18
 393 };
 394
 395 /*
 396  series for an21 on the interval -1.25000e-01 to -1.56250e-02
 397                                         with weighted error   3.43e-17
 398                                          log weighted error  16.47
 399                                significant figures required  14.48
 400                                     decimal places required  17.16
 401 */
 402 static double an21_data[24] = {
 403     0.0198313155263169394,
 404    -0.0029376249067087533,
 405    -0.0001136260695958196,
 406    -0.0000100554451087156,
 407    -0.0000013048787116563,
 408    -0.0000002123881993151,
 409    -0.0000000402270833384,
 410    -0.0000000084996745953,
 411    -0.0000000019514839426,
 412    -0.0000000004783865344,
 413    -0.0000000001236733992,
 414    -0.0000000000334137486,
 415    -0.0000000000093702824,
 416    -0.0000000000027130128,
 417    -0.0000000000008075954,
 418    -0.0000000000002463214,
 419    -0.0000000000000767656,
 420    -0.0000000000000243883,
 421    -0.0000000000000078831,
 422    -0.0000000000000025882,
 423    -0.0000000000000008619,
 424    -0.0000000000000002908,
 425    -0.0000000000000000993,
 426    -0.0000000000000000343
 427 };
 428 static cheb_series an21_cs = {
 429   an21_data,
 430   23,
 431   -1, 1,
 432   12
 433 };
 434
 435 /*
 436  series for an20 on the interval -1.56250e-02 to  0.00000e+00
 437                                         with weighted error   4.41e-17
 438                                          log weighted error  16.36
 439                                significant figures required  14.16
 440                                     decimal places required  16.96
 441 */
 442 static double an20_data[16] = {
 443     0.0126732217145738027,
 444    -0.0005212847072615621,
 445    -0.0000052672111140370,
 446    -0.0000001628202185026,
 447    -0.0000000090991442687,
 448    -0.0000000007438647126,
 449    -0.0000000000795494752,
 450    -0.0000000000104050944,
 451    -0.0000000000015932426,
 452    -0.0000000000002770648,
 453    -0.0000000000000535343,
 454    -0.0000000000000113062,
 455    -0.0000000000000025772,
 456    -0.0000000000000006278,
 457    -0.0000000000000001621,
 458    -0.0000000000000000441
 459 };
 460 static cheb_series an20_cs = {
 461   an20_data,
 462   15,
 463   -1, 1,
 464   8
 465 };
 466
 467 /*
 468  series for aph2 on the interval -1.00000e+00 to -1.25000e-01
 469                                         with weighted error   2.94e-17
 470                                          log weighted error  16.53
 471                                significant figures required  15.58
 472                                     decimal places required  17.28
 473 */
 474 static double aph2_data[32] = {
 475    -0.2057088719781465107,
 476     0.0422196961357771922,
 477     0.0020482560511207275,
 478     0.0002607800735165006,
 479     0.0000474824268004729,
 480     0.0000105102756431612,
 481     0.0000026353534014668,
 482     0.0000007208824863499,
 483     0.0000002103236664473,
 484     0.0000000644975634555,
 485     0.0000000205802377264,
 486     0.0000000067836273921,
 487     0.0000000022974015284,
 488     0.0000000007961306765,
 489     0.0000000002813860610,
 490     0.0000000001011749057,
 491     0.0000000000369306738,
 492     0.0000000000136615066,
 493     0.0000000000051142751,
 494     0.0000000000019351689,
 495     0.0000000000007393607,
 496     0.0000000000002849792,
 497     0.0000000000001107281,
 498     0.0000000000000433412,
 499     0.0000000000000170801,
 500     0.0000000000000067733,
 501     0.0000000000000027017,
 502     0.0000000000000010835,
 503     0.0000000000000004367,
 504     0.0000000000000001769,
 505     0.0000000000000000719,
 506     0.0000000000000000294
 507 };
 508 static cheb_series aph2_cs = {
 509   aph2_data,
 510   31,
 511   -1, 1,
 512   16
 513 };
 514
 515 /*
 516  series for aph1 on the interval -1.25000e-01 to -1.56250e-02
 517                                         with weighted error   6.38e-17
 518                                          log weighted error  16.20
 519                                significant figures required  14.91
 520                                     decimal places required  16.87
 521 */
 522 static double aph1_data[22] = {
 523   -0.1024172908077571694,
 524    0.0071697275146591248,
 525    0.0001209959363122329,
 526    0.0000073361512841220,
 527    0.0000007535382954272,
 528    0.0000001041478171741,
 529    0.0000000174358728519,
 530    0.0000000033399795033,
 531    0.0000000007073075174,
 532    0.0000000001619187515,
 533    0.0000000000394539982,
 534    0.0000000000101192282,
 535    0.0000000000027092778,
 536    0.0000000000007523806,
 537    0.0000000000002156369,
 538    0.0000000000000635283,
 539    0.0000000000000191757,
 540    0.0000000000000059143,
 541    0.0000000000000018597,
 542    0.0000000000000005950,
 543    0.0000000000000001934,
 544    0.0000000000000000638
 545 };
 546 static cheb_series aph1_cs = {
 547   aph1_data,
 548   21,
 549   -1, 1,
 550   10
 551 };
 552
 553 /*
 554  series for aph0 on the interval -1.56250e-02 to  0.00000e+00
 555                                         with weighted error   2.29e-17
 556                                          log weighted error  16.64
 557                                significant figures required  15.27
 558                                     decimal places required  17.23
 559 */
 560 static double aph0_data[15] = {
 561  -0.0855849241130933257,
 562   0.0011214378867065261,
 563   0.0000042721029353664,
 564   0.0000000817607381483,
 565   0.0000000033907645000,
 566   0.0000000002253264423,
 567   0.0000000000206284209,
 568   0.0000000000023858763,
 569   0.0000000000003301618,
 570   0.0000000000000527010,
 571   0.0000000000000094555,
 572   0.0000000000000018709,
 573   0.0000000000000004024,
 574   0.0000000000000000930,
 575   0.0000000000000000229
 576 };
 577 static cheb_series aph0_cs = {
 578   aph0_data,
 579   14,
 580   -1, 1,
 581   7
 582 };
 583
 584
 585 static
 586 int
 587 airy_deriv_mod_phase(const double x, gsl_mode_t mode,
 588                      gsl_sf_result * ampl, gsl_sf_result * phi)
 589 {
 590   const double pi34 = 2.356194490192344928847;
 591   gsl_sf_result result_a;
 592   gsl_sf_result result_p;
 593   double a, p;
 594   double sqx;
 595   double x32;
 596
 597   if(x <= -4.0) {
 598     double z = 128.0/(x*x*x) + 1.0;
 599     cheb_eval_mode_e(&an20_cs, z, mode, &result_a);
 600     cheb_eval_mode_e(&aph0_cs, z, mode, &result_p);
 601   }
 602   else if(x <= -2.0) {
 603     double z = (128.0/(x*x*x) + 9.0) / 7.0;
 604     cheb_eval_mode_e(&an21_cs, z, mode, &result_a);
 605     cheb_eval_mode_e(&aph1_cs, z, mode, &result_p);
 606   }
 607   else if(x <= -1.0) {
 608     double z = (16.0/(x*x*x) + 9.0) / 7.0;
 609     cheb_eval_mode_e(&an22_cs, z, mode, &result_a);
 610     cheb_eval_mode_e(&aph2_cs, z, mode, &result_p);
 611   }
 612   else {
 613     ampl->val = 0.0;
 614     ampl->err = 0.0;
 615     phi->val  = 0.0;
 616     phi->err  = 0.0;
 617     GSL_ERROR ("x is greater than 1.0", GSL_EDOM);
 618   }
 619
 620   a =  0.3125 + result_a.val;
 621   p = -0.625  + result_p.val;
 622
 623   sqx = sqrt(-x);
 624   x32   = x*sqx;
 625
 626   ampl->val = sqrt(a * sqx);
 627   ampl->err = fabs(ampl->val) * (GSL_DBL_EPSILON + fabs(result_a.err/result_a.val));
 628   phi->val  = pi34 - x * sqx * p;
 629   phi->err = fabs(phi->val) * (GSL_DBL_EPSILON + fabs(result_p.err/result_p.val));
 630
 631   return GSL_SUCCESS;
 632 }
 633
 634
 635 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 636
 637 int
 638 gsl_sf_airy_Ai_deriv_scaled_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
 639 {
 640   /* CHECK_POINTER(result) */
 641
 642   if(x < -1.0) {
 643     gsl_sf_result a;
 644     gsl_sf_result p;
 645     int status_ap = airy_deriv_mod_phase(x, mode, &a, &p);
 646     double c    = cos(p.val);
 647     result->val  = a.val * c;
 648     result->err  = fabs(result->val * p.err) + fabs(c * a.err);
 649     result->err += GSL_DBL_EPSILON * fabs(result->val);
 650     return status_ap;
 651   }
 652   else if(x <= 1.0) {
 653     const double x3 = x*x*x;
 654     const double x2 = x*x;
 655     gsl_sf_result result_c0;
 656     gsl_sf_result result_c1;
 657     cheb_eval_mode_e(&aif_cs, x3, mode, &result_c0);
 658     cheb_eval_mode_e(&aig_cs, x3, mode, &result_c1);
 659
 660     result->val  = (x2*(0.125 + result_c0.val) - result_c1.val) - 0.25;
 661     result->err  = fabs(x2*result_c0.val) + result_c1.err;
 662     result->err += GSL_DBL_EPSILON * fabs(result->val);
 663
 664     if(x > GSL_ROOT3_DBL_EPSILON*GSL_ROOT3_DBL_EPSILON) {
 665       /* scale only if x is positive */
 666       double s = exp(2.0*x*sqrt(x)/3.0);
 667       result->val *= s;
 668       result->err *= s;
 669     }
 670
 671     return GSL_SUCCESS;
 672   }
 673   else if(x <= 4.0) {
 674     const double sqrtx = sqrt(x);
 675     const double z = (16.0/(x*sqrtx) - 9.0)/7.0;
 676     const double s = sqrt(sqrtx);
 677     gsl_sf_result result_c0;
 678     cheb_eval_mode_e(&aip1_cs, z, mode, &result_c0);
 679     result->val  = -(0.28125 + result_c0.val) * s;
 680     result->err  = result_c0.err * s;
 681     result->err += GSL_DBL_EPSILON * fabs(result->val);
 682     return GSL_SUCCESS;
 683   }
 684   else {
 685     const double sqrtx = sqrt(x);
 686     const double z = 16.0/(x*sqrtx) - 1.0;
 687     const double s = sqrt(sqrtx);
 688     gsl_sf_result result_c0;
 689     cheb_eval_mode_e(&aip2_cs, z, mode, &result_c0);
 690     result->val  = -(0.28125 + result_c0.val) * s;
 691     result->err  = result_c0.err * s;
 692     result->err += GSL_DBL_EPSILON * fabs(result->val);
 693     return GSL_SUCCESS;
 694   }
 695 }
 696
 697
 698 int
 699 gsl_sf_airy_Ai_deriv_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
 700 {
 701   /* CHECK_POINTER(result) */
 702
 703   if(x < -1.0) {
 704     gsl_sf_result a;
 705     gsl_sf_result p;
 706     int status_ap = airy_deriv_mod_phase(x, mode, &a, &p);
 707     double c    = cos(p.val);
 708     result->val  = a.val * c;
 709     result->err  = fabs(result->val * p.err) + fabs(c * a.err);
 710     result->err += GSL_DBL_EPSILON * fabs(result->val);
 711     return status_ap;
 712   }
 713   else if(x < 1.0) {
 714     const double x3 = x*x*x;
 715     gsl_sf_result result_c1;
 716     gsl_sf_result result_c2;
 717     cheb_eval_mode_e(&aif_cs, x3, mode, &result_c1);
 718     cheb_eval_mode_e(&aig_cs, x3, mode, &result_c2);
 719     result->val  = (x*x*(0.125 + result_c1.val) - result_c2.val) - 0.25;
 720     result->err  = fabs(x*x*result_c1.err) + result_c2.err;
 721     result->err += GSL_DBL_EPSILON * fabs(result->val);
 722     return GSL_SUCCESS;
 723   }
 724   else if(x*x*x < 9.0/4.0 * GSL_LOG_DBL_MIN*GSL_LOG_DBL_MIN) {
 725     gsl_sf_result result_aps;
 726     const double arg = -2.0*x*sqrt(x)/3.0;
 727     const int stat_a = gsl_sf_airy_Ai_deriv_scaled_e(x, mode, &result_aps);
 728     const int stat_e = gsl_sf_exp_mult_err_e(arg, 1.5*fabs(arg*GSL_DBL_EPSILON),
 729                                                 result_aps.val, result_aps.err,
 730                                                 result);
 731     return GSL_ERROR_SELECT_2(stat_e, stat_a);
 732   }
 733   else {
 734     UNDERFLOW_ERROR(result);
 735   }
 736 }
 737
 738
 739 int
 740 gsl_sf_airy_Bi_deriv_scaled_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
 741 {
 742   const double atr =  8.7506905708484345;   /* 16./(sqrt(8)-1) */
 743   const double btr = -2.0938363213560543;   /* -(sqrt(8)+1)/(sqrt(8)-1) */
 744
 745   /* CHECK_POINTER(result) */
 746
 747   if(x < -1.0) {
 748     gsl_sf_result a;
 749     gsl_sf_result p;
 750     int status_ap = airy_deriv_mod_phase(x, mode, &a, &p);
 751     double s     = sin(p.val);
 752     result->val  = a.val * s;
 753     result->err  = fabs(result->val * p.err) + fabs(s * a.err);
 754     result->err += GSL_DBL_EPSILON * fabs(result->val);
 755     return status_ap;
 756   }
 757   else if(x < 1.0) {
 758     const double x3 = x*x*x;
 759     const double x2 = x*x;
 760     gsl_sf_result result_c1;
 761     gsl_sf_result result_c2;
 762     cheb_eval_mode_e(&bif_cs, x3, mode, &result_c1);
 763     cheb_eval_mode_e(&big_cs, x3, mode, &result_c2);
 764     result->val  = x2 * (result_c1.val + 0.25) + result_c2.val + 0.5;
 765     result->err  = x2 * result_c1.err + result_c2.err;
 766     result->err += GSL_DBL_EPSILON * fabs(result->val);
 767
 768     if(x > GSL_ROOT3_DBL_EPSILON*GSL_ROOT3_DBL_EPSILON) {
 769       /* scale only if x is positive */
 770       const double s = exp(-2.0*x*sqrt(x)/3.0);
 771       result->val *= s;
 772       result->err *= s;
 773     }
 774
 775     return GSL_SUCCESS;
 776   }
 777   else if(x < 2.0) {
 778     const double z = (2.0*x*x*x - 9.0) / 7.0;
 779     const double s = exp(-2.0*x*sqrt(x)/3.0);
 780     gsl_sf_result result_c0;
 781     gsl_sf_result result_c1;
 782     cheb_eval_mode_e(&bif2_cs, z, mode, &result_c0);
 783     cheb_eval_mode_e(&big2_cs, z, mode, &result_c1);
 784     result->val  = s * (x*x * (0.25 + result_c0.val) + 0.5 + result_c1.val);
 785     result->err  = s * (x*x * result_c0.err + result_c1.err);
 786     result->err += GSL_DBL_EPSILON * fabs(result->val);
 787     return GSL_SUCCESS;
 788   }
 789   else if(x < 4.0) {
 790     const double sqrtx = sqrt(x);
 791     const double z = atr/(x*sqrtx) + btr;
 792     const double s = sqrt(sqrtx);
 793     gsl_sf_result result_c0;
 794     cheb_eval_mode_e(&bip1_cs, z, mode, &result_c0);
 795     result->val  = s * (0.625 + result_c0.val);
 796     result->err  = s * result_c0.err;
 797     result->err += GSL_DBL_EPSILON * fabs(result->val);
 798     return GSL_SUCCESS;
 799   }
 800   else {
 801     const double sqrtx = sqrt(x);
 802     const double z = 16.0/(x*sqrtx) - 1.0;
 803     const double s = sqrt(sqrtx);
 804     gsl_sf_result result_c0;
 805     cheb_eval_mode_e(&bip2_cs, z, mode, &result_c0);
 806     result->val  = s * (0.625 + result_c0.val);
 807     result->err  = s * result_c0.err;
 808     result->err += GSL_DBL_EPSILON * fabs(result->val);
 809     return GSL_SUCCESS;
 810   }
 811 }
 812
 813
 814 int
 815 gsl_sf_airy_Bi_deriv_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
 816 {
 817   /* CHECK_POINTER(result) */
 818
 819   if(x < -1.0) {
 820     gsl_sf_result a;
 821     gsl_sf_result p;
 822     int status_ap = airy_deriv_mod_phase(x, mode, &a, &p);
 823     double s    = sin(p.val);
 824     result->val  = a.val * s;
 825     result->err  = fabs(result->val * p.err) + fabs(s * a.err);
 826     result->err += GSL_DBL_EPSILON * fabs(result->val);
 827     return status_ap;
 828   }
 829   else if(x < 1.0) {
 830     const double x3 = x*x*x;
 831     const double x2 = x*x;
 832     gsl_sf_result result_c1;
 833     gsl_sf_result result_c2;
 834     cheb_eval_mode_e(&bif_cs, x3, mode, &result_c1);
 835     cheb_eval_mode_e(&big_cs, x3, mode, &result_c2);
 836     result->val  = x2 * (result_c1.val + 0.25) + result_c2.val + 0.5;
 837     result->err  = x2 * result_c1.err + result_c2.err;
 838     result->err += GSL_DBL_EPSILON * fabs(result->val);
 839     return GSL_SUCCESS;
 840   }
 841   else if(x < 2.0) {
 842     const double z = (2.0*x*x*x - 9.0) / 7.0;
 843     gsl_sf_result result_c1;
 844     gsl_sf_result result_c2;
 845     cheb_eval_mode_e(&bif2_cs, z, mode, &result_c1);
 846     cheb_eval_mode_e(&big2_cs, z, mode, &result_c2);
 847     result->val  = x*x * (result_c1.val + 0.25) + result_c2.val + 0.5;
 848     result->err  = x*x * result_c1.err + result_c2.err;
 849     result->err += GSL_DBL_EPSILON * fabs(result->val);
 850     return GSL_SUCCESS;
 851   }
 852   else if(x < GSL_ROOT3_DBL_MAX*GSL_ROOT3_DBL_MAX) {
 853     gsl_sf_result result_bps;
 854     const double arg = 2.0*(x*sqrt(x)/3.0);
 855     int stat_b = gsl_sf_airy_Bi_deriv_scaled_e(x, mode, &result_bps);
 856     int stat_e = gsl_sf_exp_mult_err_e(arg, 1.5*fabs(arg*GSL_DBL_EPSILON),
 857                                           result_bps.val, result_bps.err,
 858                                           result);
 859     return GSL_ERROR_SELECT_2(stat_e, stat_b);
 860   }
 861   else {
 862     OVERFLOW_ERROR(result);
 863   }
 864 }
 865
 866 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 867
 868 #include "eval.h"
 869
 870 double gsl_sf_airy_Ai_deriv_scaled(const double x, gsl_mode_t mode)
 871 {
 872   EVAL_RESULT(gsl_sf_airy_Ai_deriv_scaled_e(x, mode, &result));
 873 }
 874
 875 double gsl_sf_airy_Ai_deriv(const double x, gsl_mode_t mode)
 876 {
 877   EVAL_RESULT(gsl_sf_airy_Ai_deriv_e(x, mode, &result));
 878 }
 879
 880 double gsl_sf_airy_Bi_deriv_scaled(const double x, gsl_mode_t mode)
 881 {
 882   EVAL_RESULT(gsl_sf_airy_Bi_deriv_scaled_e(x, mode, &result));
 883 }
 884
 885 double gsl_sf_airy_Bi_deriv(const double x, gsl_mode_t mode)
 886 {
 887   EVAL_RESULT(gsl_sf_airy_Bi_deriv_e(x, mode, &result));
 888 }