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

   1 /* specfunc/airy_zero.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_airy.h>
  26
  27 #include "error.h"
  28
  29 static const double zero_Ai[] = {
  30   0,
  31   -2.338107410459767039,
  32   -4.087949444130970617,
  33   -5.520559828095551059,
  34   -6.786708090071758999,
  35   -7.944133587120853123,
  36   -9.022650853340980380,
  37   -10.04017434155808593,
  38   -11.00852430373326289,
  39   -11.93601556323626252,
  40   -12.82877675286575720,
  41   -13.69148903521071793,
  42   -14.52782995177533498,
  43   -15.34075513597799686,
  44   -16.13268515694577144,
  45   -16.90563399742994263,
  46   -17.661300105697057509,
  47   -18.401132599207115416,
  48   -19.126380474246952144,
  49   -19.838129891721499701,
  50   -20.537332907677566360,
  51   -21.224829943642096955,
  52   -21.901367595585130707,
  53   -22.567612917496502831,
  54   -23.224165001121681061,
  55   -23.871564455535918567,
  56   -24.510301236589677490,
  57   -25.140821166148963748,
  58   -25.763531400982756459,
  59   -26.378805052137232374,
  60   -26.986985111606367686,
  61   -27.588387809882444812,
  62   -28.183305502632644923,
  63   -28.772009165237435382,
  64   -29.354750558766287963,
  65   -29.931764119086555913,
  66   -30.503268611418505287,
  67   -31.069468585183755604,
  68   -31.63055565801265934,
  69   -32.18670965295205069,
  70   -32.73809960900026913,
  71   -33.28488468190140188,
  72   -33.82721494950865194,
  73   -34.36523213386365906,
  74   -34.89907025034531210,
  75   -35.42885619274788846,
  76   -35.95471026189862926,
  77   -36.47674664437480896,
  78   -36.99507384699450161,
  79   -37.50979509200501613,
  80   -38.02100867725525443,
  81   -38.52880830509424882,
  82   -39.03328338327251391,
  83   -39.53451930072301805,
  84   -40.03259768075417603,
  85   -40.52759661388971821,
  86   -41.01959087233248966,
  87   -41.50865210780525018,
  88   -41.99484903432643004,
  89   -42.47824759730839188,
  90   -42.95891113021656009,
  91   -43.43690049989685412,
  92   -43.91227424156370168,
  93   -44.38508868433939023,
  94   -44.85539806814583243,
  95   -45.32325465267043011,
  96   -45.78870881905730086,
  97   -46.25180916491254629,
  98   -46.71260259315651633,
  99   -47.17113439520631705,
 100   -47.62744832892739292,
 101   -48.08158669175325711,
 102   -48.53359038933679845,
 103   -48.98349900006458366,
 104   -49.43135083573678341,
 105   -49.87718299868941729,
 106   -50.32103143561221860,
 107   -50.76293098829428522,
 108   -51.20291544151056412,
 109   -51.64101756824489758,
 110   -52.07726917242964943,
 111   -52.51170112936766183,
 112   -52.94434342398931824,
 113   -53.37522518708567514,
 114   -53.80437472964785717,
 115   -54.23181957543308298,
 116   -54.65758649186871111,
 117   -55.08170151939748312,
 118   -55.50418999935962251,
 119   -55.92507660050055598,
 120   -56.34438534418670066,
 121   -56.76213962840595327,
 122   -57.17836225062417808,
 123   -57.59307542956407782,
 124   -58.00630082596830627,
 125   -58.41805956240450934,
 126   -58.82837224216613231,
 127   -59.23725896731927534,
 128   -59.64473935594259360,
 129   -60.05083255860419805,
 130   -60.45555727411669871
 131 };
 132 static const size_t size_zero_Ai = sizeof(zero_Ai)/sizeof(double);
 133
 134
 135 static const double zero_Bi[] = {
 136   0,
 137   -1.173713222709127925,
 138   -3.271093302836352716,
 139   -4.830737841662015933,
 140   -6.169852128310251260,
 141   -7.376762079367763714,
 142   -8.491948846509388013,
 143   -9.538194379346238887,
 144   -10.52991350670535792,
 145   -11.47695355127877944,
 146   -12.38641713858273875,
 147   -13.26363952294180555,
 148   -14.11275680906865779,
 149   -14.93705741215416404,
 150   -15.739210351190482771,
 151   -16.521419550634379054,
 152   -17.285531624581242533,
 153   -18.033113287225001572,
 154   -18.765508284480081041,
 155   -19.483880132989234014,
 156   -20.189244785396202420,
 157   -20.882495994193175768,
 158   -21.564425284712977653,
 159   -22.235737881803385167,
 160   -22.897065554219793474,
 161   -23.548977079642448269,
 162   -24.191986850649000086,
 163   -24.826562012152892172,
 164   -25.453128427085131994,
 165   -26.072075698466804494,
 166   -26.683761425120990449,
 167   -27.288514830076298204,
 168   -27.886639871735962459,
 169   -28.478417925678661737,
 170   -29.064110107777650305,
 171   -29.643959295918396591,
 172   -30.218191897047274645,
 173   -30.787019397921766297,
 174   -31.350639731255585371,
 175   -31.90923848358456965,
 176   -32.46298996683685318,
 177   -33.01205817205683814,
 178   -33.55659762084006113,
 179   -34.09675412765602851,
 180   -34.63266548426775468,
 181   -35.16446207582101720,
 182   -35.69226743681080479,
 183   -36.21619875398748222,
 184   -36.73636732230120657,
 185   -37.25287895916828697,
 186   -37.76583438165180116,
 187   -38.27532955056003997,
 188   -38.78145598496327279,
 189   -39.28430105019802461,
 190   -39.78394822205711298,
 191   -40.28047732954369150,
 192   -40.77396477829068148,
 193   -41.26448375650675678,
 194   -41.75210442510106287,
 195   -42.23689409345656643,
 196   -42.71891738216253539,
 197   -43.19823637387693118,
 198   -43.67491075336673948,
 199   -44.14899793766617113,
 200   -44.62055319719727274,
 201   -45.08962976861312825,
 202   -45.55627896004907928,
 203   -46.02055024940102076,
 204   -46.48249137619078661,
 205   -46.94214842752602207,
 206   -47.39956591861496210,
 207   -47.85478686825452176,
 208   -48.30785286967246692,
 209   -48.75880415707066192,
 210   -49.20767966818603897,
 211   -49.65451710315861501,
 212   -50.09935297997125482,
 213   -50.54222268670364757,
 214   -50.98316053082286586,
 215   -51.42219978571468262,
 216   -51.85937273464332870,
 217   -52.29471071231240525,
 218   -52.72824414418606069,
 219   -53.16000258371716397,
 220   -53.59001474761792882,
 221   -54.01830854929815828,
 222   -54.44491113058688729,
 223   -54.86984889184461534,
 224   -55.29314752056546491,
 225   -55.71483201856140365,
 226   -56.13492672781406761,
 227   -56.55345535507366411,
 228   -56.97044099527886475,
 229   -57.38590615386647834,
 230   -57.79987276803497897,
 231   -58.21236222702161974,
 232   -58.62339539144885603,
 233   -59.03299261179210306,
 234   -59.44117374601743460,
 235   -59.84795817643466996,
 236   -60.25336482580837088
 237 };
 238 static const size_t size_zero_Bi = sizeof(zero_Bi)/sizeof(double);
 239
 240
 241 static const double zero_Aip[] = {
 242   0,
 243   -1.018792971647471089,
 244   -3.248197582179836738,
 245   -4.820099211178735639,
 246   -6.163307355639486822,
 247   -7.372177255047770177,
 248   -8.488486734019722133,
 249   -9.535449052433547471,
 250   -10.52766039695740728,
 251   -11.47505663348024529,
 252   -12.384788371845747325,
 253   -13.262218961665210382,
 254   -14.111501970462995282,
 255   -14.935937196720517467,
 256   -15.738201373692538303,
 257   -16.520503825433793542,
 258   -17.284695050216437357,
 259   -18.032344622504393395,
 260   -18.764798437665954740,
 261   -19.483221656567231178,
 262   -20.188631509463373154,
 263   -20.881922755516737701,
 264   -21.563887723198974958,
 265   -22.235232285348913331,
 266   -22.896588738874619001,
 267   -23.548526295928801574,
 268   -24.191559709526353841,
 269   -24.826156425921155001,
 270   -25.452742561777649948,
 271   -26.071707935173912515,
 272   -26.683410328322449767,
 273   -27.288179121523985029,
 274   -27.886318408768461192,
 275   -28.478109683102278108,
 276   -29.063814162638199090,
 277   -29.643674814632015921,
 278   -30.217918124468574603,
 279   -30.786755648012502519,
 280   -31.350385379083034671,
 281   -31.90899295843046320,
 282   -32.46275274623847982,
 283   -33.01182877663428709,
 284   -33.55637560978942190,
 285   -34.09653909480913771,
 286   -34.63245705463586589,
 287   -35.16425990255340758,
 288   -35.69207119851046870,
 289   -36.21600815233519918,
 290   -36.73618207994680321,
 291   -37.25269881785414827,
 292   -37.76565910053887108,
 293   -38.27515890473087933,
 294   -38.78128976408036876,
 295   -39.28413905729859644,
 296   -39.78379027246823278,
 297   -40.28032324990371935,
 298   -40.77381440566486637,
 299   -41.26433693758643383,
 300   -41.75196101547722703,
 301   -42.23675395695976012,
 302   -42.71878039026198233,
 303   -43.19810240513270670,
 304   -43.67477969292950869,
 305   -44.14886967681966886,
 306   -44.62042763293925724,
 307   -45.08950680327102630,
 308   -45.55615850092696446,
 309   -46.02043220845493728,
 310   -46.48237566972975615,
 311   -46.94203497593635633,
 312   -47.39945464610575493,
 313   -47.85467770262241617,
 314   -48.30774574208398774,
 315   -48.75869900186057804,
 316   -49.20757642267037247,
 317   -49.65441570746105074,
 318   -50.09925337686182515,
 319   -50.54212482144867502,
 320   -50.98306435104524282,
 321   -51.42210524126365311,
 322   -51.85927977747301469,
 323   -52.29461929636838876,
 324   -52.72815422529939506,
 325   -53.15991411950524351,
 326   -53.58992769739169611,
 327   -54.01822287397517367,
 328   -54.44482679260982599,
 329   -54.86976585510479430,
 330   -55.29306575033103518,
 331   -55.71475148140987392,
 332   -56.13484739156885235,
 333   -56.55337718874437424,
 334   -56.97036396900508167,
 335   -57.38583023886477265,
 336   -57.79979793654895377,
 337   -58.21228845227477578,
 338   -58.62332264760009139,
 339   -59.03292087389367419,
 340   -59.44110298997521892,
 341   -59.84788837897058171,
 342   -60.25329596442479317
 343 };
 344 static const size_t size_zero_Aip = sizeof(zero_Aip)/sizeof(double);
 345
 346
 347 static const double zero_Bip[] = {
 348   0,
 349   -2.294439682614123247,
 350   -4.073155089071828216,
 351   -5.512395729663599496,
 352   -6.781294445990305390,
 353   -7.940178689168578927,
 354   -9.019583358794239067,
 355   -10.037696334908545802,
 356   -11.006462667712289940,
 357   -11.934261645014844663,
 358   -12.827258309177217640,
 359   -13.690155826835049101,
 360   -14.526645763485711410,
 361   -15.339693082242404109,
 362   -16.131724782385900578,
 363   -16.904759411889649958,
 364   -17.660498743114976102,
 365   -18.400394367181703280,
 366   -19.125697156412638066,
 367   -19.837494718415910503,
 368   -20.536740241453273980,
 369   -21.224275044889266569,
 370   -21.900846445139208281,
 371   -22.567122080497200470,
 372   -23.223701521208962116,
 373   -23.871125771677973595,
 374   -24.509885117016242729,
 375   -25.140425655367878908,
 376   -25.763154776913454319,
 377   -26.378445791146615697,
 378   -26.986641859775034987,
 379   -27.588059359225600573,
 380   -28.182990771292975456,
 381   -28.771707180886056250,
 382   -29.354460444612957224,
 383   -29.931485082026055160,
 384   -30.502999931936645516,
 385   -31.069209608721234058,
 386   -31.63030578754333679,
 387   -32.18646834257807369,
 388   -32.73786635840274752,
 389   -33.28465903151424981,
 390   -33.82699647630635587,
 391   -34.36502044767239631,
 392   -34.89886499060196419,
 393   -35.42865702564380962,
 394   -35.95451687785511190,
 395   -36.47655875580547918,
 396   -36.99489118631672770,
 397   -37.50961740986809593,
 398   -38.02083574095788210
 399 };
 400 static const size_t size_zero_Bip = sizeof(zero_Bip)/sizeof(double);
 401
 402
 403
 404 /* [Abramowitz+Stegun, 10.4.105] */
 405 static double
 406 zero_f(double z)
 407 {
 408   const double pre = pow(z, 2.0/3.0);
 409   const double zi2 = 1.0/(z*z);
 410   const double zi4 = zi2 * zi2;
 411   const double t1  =  5.0/48.0 * zi2;
 412   const double t2  = -5.0/36.0 * zi4;
 413   const double t3  =  77125.0/82944.0 * zi4 * zi2;
 414   const double t4  = -108056875.0/6967296.0 * zi4 * zi4;
 415   return pre * (1.0 + t1 + t2 + t3 + t4);
 416
 417 }
 418 static double
 419 zero_g(double z)
 420 {
 421   const double pre = pow(z, 2.0/3.0);
 422   const double zi2 = 1.0/(z*z);
 423   const double zi4 = zi2 * zi2;
 424   const double t1  = -7.0/48.0 * zi2;
 425   const double t2  =  35.0/288.0 * zi4;
 426   const double t3  = -181223.0/207360.0 * zi4 * zi2;
 427   const double t4  =  18683371.0/1244160.0 * zi4 * zi4;
 428   return pre * (1.0 + t1 + t2 + t3 + t4);
 429 }
 430
 431
 432
 433 int
 434 gsl_sf_airy_zero_Ai_e(unsigned int s, gsl_sf_result * result)
 435 {
 436   /* CHECK_POINTER(result) */
 437
 438   if(s < 1) {
 439     DOMAIN_ERROR_MSG("s is less than 1", result);
 440   }
 441   else if(s < size_zero_Ai) {
 442     result->val = zero_Ai[s];
 443     result->err = GSL_DBL_EPSILON * fabs(result->val);
 444     return GSL_SUCCESS;
 445   }
 446   else {
 447     const double z = 3.0*M_PI/8.0 * (4.0*s - 1.0);
 448     const double f = zero_f(z);
 449     result->val = -f;
 450     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 451     return GSL_SUCCESS;
 452   }
 453 }
 454
 455
 456 int
 457 gsl_sf_airy_zero_Bi_e(unsigned int s, gsl_sf_result * result)
 458 {
 459   /* CHECK_POINTER(result) */
 460
 461   if(s < 1) {
 462     DOMAIN_ERROR_MSG("s is less than 1", result);
 463   }
 464   else if(s < size_zero_Bi) {
 465     result->val = zero_Bi[s];
 466     result->err = GSL_DBL_EPSILON * fabs(result->val);
 467     return GSL_SUCCESS;
 468   }
 469   else {
 470     const double z = 3.0*M_PI/8.0 * (4.0*s - 3.0);
 471     const double f = zero_f(z);
 472     result->val = -f;
 473     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 474     return GSL_SUCCESS;
 475   }
 476 }
 477
 478
 479 int
 480 gsl_sf_airy_zero_Ai_deriv_e(unsigned int s, gsl_sf_result * result)
 481 {
 482   /* CHECK_POINTER(result) */
 483
 484   if(s < 1) {
 485     DOMAIN_ERROR_MSG("s is less than 1", result);
 486   }
 487   else if(s < size_zero_Aip) {
 488     result->val = zero_Aip[s];
 489     result->err = GSL_DBL_EPSILON * fabs(result->val);
 490     return GSL_SUCCESS;
 491   }
 492   else {
 493     const double z = 3.0*M_PI/8.0 * (4.0*s - 3.0);
 494     const double g = zero_g(z);
 495     result->val = -g;
 496     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 497     return GSL_SUCCESS;
 498   }
 499 }
 500
 501
 502 int
 503 gsl_sf_airy_zero_Bi_deriv_e(unsigned int s, gsl_sf_result * result)
 504 {
 505   /* CHECK_POINTER(result) */
 506
 507   if(s < 1) {
 508     DOMAIN_ERROR_MSG("s is less than 1", result);
 509   }
 510   else if(s < size_zero_Bip) {
 511     result->val = zero_Bip[s];
 512     result->err = GSL_DBL_EPSILON * fabs(result->val);
 513     return GSL_SUCCESS;
 514   }
 515   else {
 516     const double z = 3.0*M_PI/8.0 * (4.0*s - 1.0);
 517     const double g = zero_g(z);
 518     result->val = -g;
 519     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 520     return GSL_SUCCESS;
 521   }
 522 }
 523
 524 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 525
 526 #include "eval.h"
 527
 528 double gsl_sf_airy_zero_Ai(unsigned int s)
 529 {
 530   EVAL_RESULT(gsl_sf_airy_zero_Ai_e(s, &result));
 531 }
 532
 533 double gsl_sf_airy_zero_Bi(unsigned int s)
 534 {
 535   EVAL_RESULT(gsl_sf_airy_zero_Bi_e(s, &result));
 536 }
 537
 538 double gsl_sf_airy_zero_Ai_deriv(unsigned int s)
 539 {
 540   EVAL_RESULT(gsl_sf_airy_zero_Ai_deriv_e(s, &result));
 541 }
 542
 543 double gsl_sf_airy_zero_Bi_deriv(unsigned int s)
 544 {
 545   EVAL_RESULT(gsl_sf_airy_zero_Bi_deriv_e(s, &result));
 546 }