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

   1 /* specfunc/bessel_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 #include <gsl/gsl_sf_pow_int.h>
  27 #include <gsl/gsl_sf_bessel.h>
  28
  29 #include "error.h"
  30
  31 #include "bessel_olver.h"
  32
  33 /* For Chebyshev expansions of the roots as functions of nu,
  34  * see [G. Nemeth, Mathematical Approximation of Special Functions].
  35  * This gives the fits for all nu and s <= 10.
  36  * I made the fits for other values of s myself [GJ].
  37  */
  38
  39 /* Chebyshev expansion: j_{nu,1} = c_k T_k*(nu/2), nu <= 2 */
  40 static const double coef_jnu1_a[] = {
  41   3.801775243633476,
  42   1.360704737511120,
  43  -0.030707710261106,
  44   0.004526823746202,
  45  -0.000808682832134,
  46   0.000159218792489,
  47  -0.000033225189761,
  48   0.000007205599763,
  49  -0.000001606110397,
  50   0.000000365439424,
  51  -0.000000084498039,
  52   0.000000019793815,
  53  -0.000000004687054,
  54   0.000000001120052,
  55  -0.000000000269767,
  56   0.000000000065420,
  57  -0.000000000015961,
  58   0.000000000003914,
  59  -0.000000000000965,
  60   0.000000000000239,
  61  -0.000000000000059,
  62   0.000000000000015,
  63  -0.000000000000004,
  64   0.000000000000001
  65 };
  66
  67
  68 /* Chebyshev expansion: j_{nu,1} = nu c_k T_k*((2/nu)^(2/3)), nu >= 2 */
  69 static const double coef_jnu1_b[] = {
  70   1.735063412537096,
  71   0.784478100951978,
  72   0.048881473180370,
  73  -0.000578279783021,
  74  -0.000038984957864,
  75   0.000005758297879,
  76  -0.000000327583229,
  77  -0.000000003853878,
  78   0.000000002284653,
  79  -0.000000000153079,
  80  -0.000000000000895,
  81   0.000000000000283,
  82   0.000000000000043,
  83   0.000000000000010,
  84  -0.000000000000003
  85 };
  86
  87
  88 /* Chebyshev expansion: j_{nu,2} = c_k T_k*(nu/2), nu <= 2 */
  89 static const double coef_jnu2_a[] = {
  90   6.992370244046161,
  91   1.446379282056534,
  92  -0.023458616207293,
  93   0.002172149448700,
  94  -0.000246262775620,
  95   0.000030990180959,
  96  -0.000004154183047,
  97   0.000000580766328,
  98  -0.000000083648175,
  99   0.000000012317355,
 100  -0.000000001844887,
 101   0.000000000280076,
 102  -0.000000000042986,
 103   0.000000000006658,
 104  -0.000000000001039,
 105   0.000000000000163,
 106  -0.000000000000026,
 107   0.000000000000004,
 108  -0.000000000000001
 109 };
 110
 111
 112 /* Chebyshev expansion: j_{nu,2} = nu c_k T_k*((2/nu)^(2/3)), nu >= 2 */
 113 static const double coef_jnu2_b[] = {
 114   2.465611864263400,
 115   1.607952988471069,
 116   0.138758034431497,
 117  -0.003687791182054,
 118  -0.000051276007868,
 119   0.000045113570749,
 120  -0.000007579172152,
 121   0.000000736469208,
 122  -0.000000011118527,
 123  -0.000000011919884,
 124   0.000000002696788,
 125  -0.000000000314488,
 126   0.000000000008124,
 127   0.000000000005211,
 128  -0.000000000001292,
 129   0.000000000000158,
 130  -0.000000000000004,
 131  -0.000000000000003,
 132   0.000000000000001
 133 };
 134
 135
 136 /* Chebyshev expansion: j_{nu,3} = c_k T_k*(nu/3), nu <= 3 */
 137 static const double coef_jnu3_a[] = {
 138   10.869647065239236,
 139    2.177524286141710,
 140   -0.034822817125293,
 141    0.003167249102413,
 142   -0.000353960349344,
 143    0.000044039086085,
 144   -0.000005851380981,
 145    0.000000812575483,
 146   -0.000000116463617,
 147    0.000000017091246,
 148   -0.000000002554376,
 149    0.000000000387335,
 150   -0.000000000059428,
 151    0.000000000009207,
 152   -0.000000000001438,
 153    0.000000000000226,
 154   -0.000000000000036,
 155    0.000000000000006,
 156   -0.000000000000001
 157 };
 158
 159
 160 /* Chebyshev expansion: j_{nu,3} = nu c_k T_k*((3/nu)^(2/3)), nu >= 3 */
 161 static const double coef_jnu3_b[] = {
 162   2.522816775173244,
 163   1.673199424973720,
 164   0.146431617506314,
 165  -0.004049001763912,
 166  -0.000039517767244,
 167   0.000048781729288,
 168  -0.000008729705695,
 169   0.000000928737310,
 170  -0.000000028388244,
 171  -0.000000012927432,
 172   0.000000003441008,
 173  -0.000000000471695,
 174   0.000000000025590,
 175   0.000000000005502,
 176  -0.000000000001881,
 177   0.000000000000295,
 178  -0.000000000000020,
 179  -0.000000000000003,
 180   0.000000000000001
 181 };
 182
 183
 184 /* Chebyshev expansion: j_{nu,4} = c_k T_k*(nu/4), nu <= 4 */
 185 static const double coef_jnu4_a[] = {
 186   14.750310252773009,
 187    2.908010932941708,
 188   -0.046093293420315,
 189    0.004147172321412,
 190   -0.000459092310473,
 191    0.000056646951906,
 192   -0.000007472351546,
 193    0.000001031210065,
 194   -0.000000147008137,
 195    0.000000021475218,
 196   -0.000000003197208,
 197    0.000000000483249,
 198   -0.000000000073946,
 199    0.000000000011431,
 200   -0.000000000001782,
 201    0.000000000000280,
 202   -0.000000000000044,
 203    0.000000000000007,
 204   -0.000000000000001
 205 };
 206
 207
 208 /* Chebyshev expansion: j_{nu,4} = nu c_k T_k*((4/nu)^(2/3)), nu >= 4 */
 209 static const double coef_jnu4_b[] = {
 210   2.551681323117914,
 211   1.706177978336572,
 212   0.150357658406131,
 213  -0.004234001378590,
 214  -0.000033854229898,
 215   0.000050763551485,
 216  -0.000009337464057,
 217   0.000001029717834,
 218  -0.000000037474196,
 219  -0.000000013450153,
 220   0.000000003836180,
 221  -0.000000000557404,
 222   0.000000000035748,
 223   0.000000000005487,
 224  -0.000000000002187,
 225   0.000000000000374,
 226  -0.000000000000031,
 227  -0.000000000000003,
 228   0.000000000000001
 229 };
 230
 231
 232
 233 /* Chebyshev expansion: j_{nu,5} = c_k T_k*(nu/5), nu <= 5 */
 234 static const double coef_jnu5_a[] = {
 235   18.632261081028211,
 236    3.638249012596966,
 237   -0.057329705998828,
 238    0.005121709126820,
 239   -0.000563325259487,
 240    0.000069100826174,
 241   -0.000009066603030,
 242    0.000001245181383,
 243   -0.000000176737282,
 244    0.000000025716695,
 245   -0.000000003815184,
 246    0.000000000574839,
 247   -0.000000000087715,
 248    0.000000000013526,
 249   -0.000000000002104,
 250    0.000000000000330,
 251   -0.000000000000052,
 252    0.000000000000008,
 253   -0.000000000000001
 254 };
 255
 256
 257 /* Chebyshev expansion: j_{nu,5} = nu c_k T_k*((5/nu)^(2/3)), nu >= 5 */
 258 /* FIXME: There is something wrong with this fit, in about the
 259  * 9th or 10th decimal place.
 260  */
 261 static const double coef_jnu5_b[] = {
 262   2.569079487591442,
 263   1.726073360882134,
 264   0.152740776809531,
 265  -0.004346449660148,
 266  -0.000030512461856,
 267   0.000052000821080,
 268  -0.000009713343981,
 269   0.000001091997863,
 270  -0.000000043061707,
 271  -0.000000013779413,
 272   0.000000004082870,
 273  -0.000000000611259,
 274   0.000000000042242,
 275   0.000000000005448,
 276  -0.000000000002377,
 277   0.000000000000424,
 278  -0.000000000000038,
 279  -0.000000000000002,
 280   0.000000000000002
 281 };
 282
 283
 284 /* Chebyshev expansion: j_{nu,6} = c_k T_k*(nu/6), nu <= 6 */
 285 static const double coef_jnu6_a[] = {
 286   22.514836143374042,
 287    4.368367257557198,
 288   -0.068550155285562,
 289    0.006093776505822,
 290   -0.000667152784957,
 291    0.000081486022398,
 292   -0.000010649011647,
 293    0.000001457089679,
 294   -0.000000206105082,
 295    0.000000029894724,
 296   -0.000000004422012,
 297    0.000000000664471,
 298   -0.000000000101140,
 299    0.000000000015561,
 300   -0.000000000002416,
 301    0.000000000000378,
 302   -0.000000000000060,
 303    0.000000000000009,
 304   -0.000000000000002
 305 };
 306
 307
 308 /* Chebyshev expansion: j_{nu,6} = nu c_k T_k*((6/nu)^(2/3)), nu >= 6 */
 309 static const double coef_jnu6_b[] = {
 310   2.580710285494837,
 311   1.739380728566154,
 312   0.154340696401691,
 313  -0.004422028860168,
 314  -0.000028305272624,
 315   0.000052845975269,
 316  -0.000009968794373,
 317   0.000001134252926,
 318  -0.000000046841241,
 319  -0.000000014007555,
 320   0.000000004251816,
 321  -0.000000000648213,
 322   0.000000000046728,
 323   0.000000000005414,
 324  -0.000000000002508,
 325   0.000000000000459,
 326  -0.000000000000043,
 327  -0.000000000000002,
 328   0.000000000000002
 329 };
 330
 331
 332 /* Chebyshev expansion: j_{nu,7} = c_k T_k*(nu/7), nu <= 7 */
 333 static const double coef_jnu7_a[] = {
 334   26.397760539730869,
 335    5.098418721711790,
 336   -0.079761896398948,
 337    0.007064521280487,
 338   -0.000770766522482,
 339    0.000093835449636,
 340   -0.000012225308542,
 341    0.000001667939800,
 342   -0.000000235288157,
 343    0.000000034040347,
 344   -0.000000005023142,
 345    0.000000000753101,
 346   -0.000000000114389,
 347    0.000000000017564,
 348   -0.000000000002722,
 349    0.000000000000425,
 350   -0.000000000000067,
 351    0.000000000000011,
 352   -0.000000000000002
 353 };
 354
 355
 356 /* Chebyshev expansion: j_{nu,7} = nu c_k T_k*((7/nu)^(2/3)), nu >= 7 */
 357 static const double coef_jnu7_b[] = {
 358   2.589033335856773,
 359   1.748907007612678,
 360   0.155488900387653,
 361  -0.004476317805688,
 362  -0.000026737952924,
 363   0.000053459680946,
 364  -0.000010153699240,
 365   0.000001164804272,
 366  -0.000000049566917,
 367  -0.000000014175403,
 368   0.000000004374840,
 369  -0.000000000675135,
 370   0.000000000050004,
 371   0.000000000005387,
 372  -0.000000000002603,
 373   0.000000000000485,
 374  -0.000000000000047,
 375  -0.000000000000002,
 376   0.000000000000002
 377 };
 378
 379
 380 /* Chebyshev expansion: j_{nu,8} = c_k T_k*(nu/8), nu <= 8 */
 381 static const double coef_jnu8_a[] = {
 382   30.280900001606662,
 383    5.828429205461221,
 384   -0.090968381181069,
 385    0.008034479731033,
 386   -0.000874254899080,
 387    0.000106164151611,
 388   -0.000013798098749,
 389    0.000001878187386,
 390   -0.000000264366627,
 391    0.000000038167685,
 392   -0.000000005621060,
 393    0.000000000841165,
 394   -0.000000000127538,
 395    0.000000000019550,
 396   -0.000000000003025,
 397    0.000000000000472,
 398   -0.000000000000074,
 399    0.000000000000012,
 400   -0.000000000000002
 401 };
 402
 403
 404 /* Chebyshev expansion: j_{nu,8} = nu c_k T_k*((8/nu)^(2/3)), nu >= 8 */
 405 static const double coef_jnu8_b[] = {
 406   2.595283877150078,
 407   1.756063044986928,
 408   0.156352972371030,
 409  -0.004517201896761,
 410  -0.000025567187878,
 411   0.000053925472558,
 412  -0.000010293734486,
 413   0.000001187923085,
 414  -0.000000051625122,
 415  -0.000000014304212,
 416   0.000000004468450,
 417  -0.000000000695620,
 418   0.000000000052500,
 419   0.000000000005367,
 420  -0.000000000002676,
 421   0.000000000000505,
 422  -0.000000000000050,
 423  -0.000000000000002,
 424   0.000000000000002
 425 };
 426
 427
 428 /* Chebyshev expansion: j_{nu,9} = c_k T_k*(nu/9), nu <= 9 */
 429 static const double coef_jnu9_a[] = {
 430   34.164181213238386,
 431    6.558412747925228,
 432   -0.102171455365016,
 433    0.009003934361201,
 434   -0.000977663914535,
 435    0.000118479876579,
 436   -0.000015368714220,
 437    0.000002088064285,
 438   -0.000000293381154,
 439    0.000000042283900,
 440   -0.000000006217033,
 441    0.000000000928887,
 442   -0.000000000140627,
 443    0.000000000021526,
 444   -0.000000000003326,
 445    0.000000000000518,
 446   -0.000000000000081,
 447    0.000000000000013,
 448   -0.000000000000002
 449 };
 450
 451
 452 /* Chebyshev expansion: j_{nu,9} = nu c_k T_k*((9/nu)^(2/3)), nu >= 9 */
 453 static const double coef_jnu9_b[] = {
 454   2.600150240905079,
 455   1.761635491694032,
 456   0.157026743724010,
 457  -0.004549100368716,
 458  -0.000024659248617,
 459   0.000054291035068,
 460  -0.000010403464334,
 461   0.000001206027524,
 462  -0.000000053234089,
 463  -0.000000014406241,
 464   0.000000004542078,
 465  -0.000000000711728,
 466   0.000000000054464,
 467   0.000000000005350,
 468  -0.000000000002733,
 469   0.000000000000521,
 470  -0.000000000000052,
 471  -0.000000000000002,
 472   0.000000000000002
 473 };
 474
 475
 476 /* Chebyshev expansion: j_{nu,10} = c_k T_k*(nu/10), nu <= 10 */
 477 static const double coef_jnu10_a[] = {
 478   38.047560766184647,
 479    7.288377637926008,
 480   -0.113372193277897,
 481    0.009973047509098,
 482   -0.001081019701335,
 483    0.000130786983847,
 484   -0.000016937898538,
 485    0.000002297699179,
 486   -0.000000322354218,
 487    0.000000046392941,
 488   -0.000000006811759,
 489    0.000000001016395,
 490   -0.000000000153677,
 491    0.000000000023486,
 492   -0.000000000003616,
 493    0.000000000000561,
 494   -0.000000000000095,
 495    0.000000000000027,
 496   -0.000000000000013,
 497    0.000000000000005
 498 };
 499
 500
 501 /* Chebyshev expansion: j_{nu,10} = nu c_k T_k*((10/nu)^(2/3)), nu >= 10 */
 502 static const double coef_jnu10_b[] = {
 503   2.604046346867949,
 504   1.766097596481182,
 505   0.157566834446511,
 506  -0.004574682244089,
 507  -0.000023934500688,
 508   0.000054585558231,
 509  -0.000010491765415,
 510   0.000001220589364,
 511  -0.000000054526331,
 512  -0.000000014489078,
 513   0.000000004601510,
 514  -0.000000000724727,
 515   0.000000000056049,
 516   0.000000000005337,
 517  -0.000000000002779,
 518   0.000000000000533,
 519  -0.000000000000054,
 520  -0.000000000000002,
 521   0.000000000000002
 522 };
 523
 524
 525 /* Chebyshev expansion: j_{nu,11} = c_k T_k*(nu/22), nu <= 22 */
 526 static const double coef_jnu11_a[] = {
 527   49.5054081076848637,
 528   15.33692279367165101,
 529  -0.33677234163517130,
 530   0.04623235772920729,
 531  -0.00781084960665093,
 532   0.00147217395434708,
 533  -0.00029695043846867,
 534   0.00006273356860235,
 535  -0.00001370575125628,
 536   3.07171282012e-6,
 537  -7.0235041249e-7,
 538   1.6320559339e-7,
 539  -3.843117306e-8,
 540   9.15083800e-9,
 541  -2.19957642e-9,
 542   5.3301703e-10,
 543  -1.3007541e-10,
 544   3.193827e-11,
 545  -7.88605e-12,
 546   1.95918e-12,
 547  -4.9020e-13,
 548   1.2207e-13,
 549  -2.820e-14,
 550   5.25e-15,
 551  -1.88e-15,
 552   2.80e-15,
 553  -2.45e-15
 554 };
 555
 556
 557 /* Chebyshev expansion: j_{nu,12} = c_k T_k*(nu/24), nu <= 24 */
 558 static const double coef_jnu12_a[] = {
 559   54.0787833216641519,
 560   16.7336367772863598,
 561  -0.36718411124537953,
 562   0.05035523375053820,
 563  -0.00849884978867533,
 564   0.00160027692813434,
 565  -0.00032248114889921,
 566   0.00006806354127199,
 567  -0.00001485665901339,
 568   3.32668783672e-6,
 569  -7.5998952729e-7,
 570   1.7644939709e-7,
 571  -4.151538210e-8,
 572   9.87722772e-9,
 573  -2.37230133e-9,
 574   5.7442875e-10,
 575  -1.4007767e-10,
 576   3.437166e-11,
 577  -8.48215e-12,
 578   2.10554e-12,
 579  -5.2623e-13,
 580   1.3189e-13,
 581  -3.175e-14,
 582   5.73e-15,
 583   5.6e-16,
 584  -8.7e-16,
 585  -6.5e-16
 586 };
 587
 588
 589 /* Chebyshev expansion: j_{nu,13} = c_k T_k*(nu/26), nu <= 26 */
 590 static const double coef_jnu13_a[] = {
 591   58.6521941921708890,
 592   18.1303398137970284,
 593  -0.39759381380126650,
 594   0.05447765240465494,
 595  -0.00918674227679980,
 596   0.00172835361420579,
 597  -0.00034800528297612,
 598   0.00007339183835188,
 599  -0.00001600713368099,
 600   3.58154960392e-6,
 601  -8.1759873497e-7,
 602   1.8968523220e-7,
 603  -4.459745253e-8,
 604   1.060304419e-8,
 605  -2.54487624e-9,
 606   6.1580214e-10,
 607  -1.5006751e-10,
 608   3.679707e-11,
 609  -9.07159e-12,
 610   2.24713e-12,
 611  -5.5943e-13,
 612   1.4069e-13,
 613  -3.679e-14,
 614   1.119e-14,
 615  -4.99e-15,
 616   3.43e-15,
 617  -2.85e-15,
 618   2.3e-15,
 619  -1.7e-15,
 620   8.7e-16
 621 };
 622
 623
 624 /* Chebyshev expansion: j_{nu,14} = c_k T_k*(nu/28), nu <= 28 */
 625 static const double coef_jnu14_a[] = {
 626   63.2256329577315566,
 627   19.5270342832914901,
 628  -0.42800190567884337,
 629   0.05859971627729398,
 630  -0.00987455163523582,
 631   0.00185641011402081,
 632  -0.00037352439419968,
 633   0.00007871886257265,
 634  -0.00001715728110045,
 635   3.83632624437e-6,
 636  -8.7518558668e-7,
 637   2.0291515353e-7,
 638  -4.767795233e-8,
 639   1.132844415e-8,
 640  -2.71734219e-9,
 641   6.5714886e-10,
 642  -1.6005342e-10,
 643   3.922557e-11,
 644  -9.66637e-12,
 645   2.39379e-12,
 646  -5.9541e-13,
 647   1.4868e-13,
 648  -3.726e-14,
 649   9.37e-15,
 650  -2.36e-15,
 651   6.0e-16
 652 };
 653
 654
 655 /* Chebyshev expansion: j_{nu,15} = c_k T_k*(nu/30), nu <= 30 */
 656 static const double coef_jnu15_a[] = {
 657   67.7990939565631635,
 658   20.9237219226859859,
 659  -0.45840871823085836,
 660   0.06272149946755639,
 661  -0.01056229551143042,
 662   0.00198445078693100,
 663  -0.00039903958650729,
 664   0.00008404489865469,
 665  -0.00001830717574922,
 666   4.09103745566e-6,
 667  -9.3275533309e-7,
 668   2.1614056403e-7,
 669  -5.075725222e-8,
 670   1.205352081e-8,
 671  -2.88971837e-9,
 672   6.9846848e-10,
 673  -1.7002946e-10,
 674   4.164941e-11,
 675  -1.025859e-11,
 676   2.53921e-12,
 677  -6.3128e-13,
 678   1.5757e-13,
 679  -3.947e-14,
 680   9.92e-15,
 681  -2.50e-15,
 682   6.3e-16
 683 };
 684
 685
 686 /* Chebyshev expansion: j_{nu,16} = c_k T_k*(nu/32), nu <= 32 */
 687 static const double coef_jnu16_a[] = {
 688   72.3725729616724770,
 689   22.32040402918608585,
 690  -0.48881449782358690,
 691   0.06684305681828766,
 692  -0.01124998690363398,
 693   0.00211247882775445,
 694  -0.00042455166484632,
 695   0.00008937015316346,
 696  -0.00001945687139551,
 697   4.34569739281e-6,
 698  -9.9031173548e-7,
 699   2.2936247195e-7,
 700  -5.383562595e-8,
 701   1.277835103e-8,
 702  -3.06202860e-9,
 703   7.3977037e-10,
 704  -1.8000071e-10,
 705   4.407196e-11,
 706  -1.085046e-11,
 707   2.68453e-12,
 708  -6.6712e-13,
 709   1.6644e-13,
 710  -4.168e-14,
 711   1.047e-14,
 712  -2.64e-15,
 713   6.7e-16
 714 };
 715
 716
 717 /* Chebyshev expansion: j_{nu,17} = c_k T_k*(nu/34), nu <= 34 */
 718 static const double coef_jnu17_a[] = {
 719   76.9460667535209549,
 720   23.71708159112252670,
 721  -0.51921943142405352,
 722   0.07096442978067622,
 723  -0.01193763559341369,
 724   0.00224049662974902,
 725  -0.00045006122941781,
 726   0.00009469477941684,
 727  -0.00002060640777107,
 728   4.60031647195e-6,
 729  -1.04785755046e-6,
 730   2.4258161247e-7,
 731  -5.691327087e-8,
 732   1.350298805e-8,
 733  -3.23428733e-9,
 734   7.8105847e-10,
 735  -1.8996825e-10,
 736   4.649350e-11,
 737  -1.144205e-11,
 738   2.82979e-12,
 739  -7.0294e-13,
 740   1.7531e-13,
 741  -4.388e-14,
 742   1.102e-14,
 743  -2.78e-15,
 744   7.0e-16
 745 };
 746
 747
 748 /* Chebyshev expansion: j_{nu,18} = c_k T_k*(nu/36), nu <= 36 */
 749 static const double coef_jnu18_a[] = {
 750   81.5195728368096659,
 751   25.11375537470259305,
 752  -0.54962366347317668,
 753   0.07508565026117689,
 754  -0.01262524908033818,
 755   0.00236850602019778,
 756  -0.00047556873651929,
 757   0.00010001889347161,
 758  -0.00002175581482429,
 759   4.85490251239e-6,
 760  -1.10539483940e-6,
 761   2.5579853343e-7,
 762  -5.999033352e-8,
 763   1.422747129e-8,
 764  -3.40650521e-9,
 765   8.2233565e-10,
 766  -1.9993286e-10,
 767   4.891426e-11,
 768  -1.203343e-11,
 769   2.97498e-12,
 770  -7.3875e-13,
 771   1.8418e-13,
 772  -4.608e-14,
 773   1.157e-14,
 774  -2.91e-15,
 775   7.4e-16
 776 };
 777
 778
 779 /* Chebyshev expansion: j_{nu,19} = c_k T_k*(nu/38), nu <= 38 */
 780 static const double coef_jnu19_a[] = {
 781   86.0930892477047512,
 782   26.51042598308271729,
 783  -0.58002730731948358,
 784   0.07920674321589394,
 785  -0.01331283320930301,
 786   0.00249650841778073,
 787  -0.00050107453900793,
 788   0.00010534258471335,
 789  -0.00002290511552874,
 790   5.10946148897e-6,
 791  -1.16292517157e-6,
 792   2.6901365037e-7,
 793  -6.306692473e-8,
 794   1.495183048e-8,
 795  -3.57869025e-9,
 796   8.6360410e-10,
 797  -2.0989514e-10,
 798   5.133439e-11,
 799  -1.262465e-11,
 800   3.12013e-12,
 801  -7.7455e-13,
 802   1.9304e-13,
 803  -4.829e-14,
 804   1.212e-14,
 805  -3.05e-15,
 806   7.7e-16
 807 };
 808
 809
 810 /* Chebyshev expansion: j_{nu,20} = c_k T_k*(nu/40), nu <= 40 */
 811 static const double coef_jnu20_a[] = {
 812   90.6666144195163770,
 813   27.9070938975436823,
 814  -0.61043045315390591,
 815   0.08332772844325554,
 816  -0.01400039260208282,
 817   0.00262450494035660,
 818  -0.00052657891389470,
 819   0.00011066592304919,
 820  -0.00002405432778364,
 821   5.36399803946e-6,
 822  -1.22044976064e-6,
 823   2.8222728362e-7,
 824  -6.614312964e-8,
 825   1.567608839e-8,
 826  -3.75084856e-9,
 827   9.0486546e-10,
 828  -2.1985553e-10,
 829   5.375401e-11,
 830  -1.321572e-11,
 831   3.26524e-12,
 832  -8.1033e-13,
 833   2.0190e-13,
 834  -5.049e-14,
 835   1.267e-14,
 836  -3.19e-15,
 837   8.0e-16,
 838  -2.0e-16
 839 };
 840
 841
 842 static const double * coef_jnu_a[] = {
 843   0,
 844   coef_jnu1_a,
 845   coef_jnu2_a,
 846   coef_jnu3_a,
 847   coef_jnu4_a,
 848   coef_jnu5_a,
 849   coef_jnu6_a,
 850   coef_jnu7_a,
 851   coef_jnu8_a,
 852   coef_jnu9_a,
 853   coef_jnu10_a,
 854   coef_jnu11_a,
 855   coef_jnu12_a,
 856   coef_jnu13_a,
 857   coef_jnu14_a,
 858   coef_jnu15_a,
 859   coef_jnu16_a,
 860   coef_jnu17_a,
 861   coef_jnu18_a,
 862   coef_jnu19_a,
 863   coef_jnu20_a
 864 };
 865
 866 static const size_t size_jnu_a[] = {
 867   0,
 868   sizeof(coef_jnu1_a)/sizeof(double),
 869   sizeof(coef_jnu2_a)/sizeof(double),
 870   sizeof(coef_jnu3_a)/sizeof(double),
 871   sizeof(coef_jnu4_a)/sizeof(double),
 872   sizeof(coef_jnu5_a)/sizeof(double),
 873   sizeof(coef_jnu6_a)/sizeof(double),
 874   sizeof(coef_jnu7_a)/sizeof(double),
 875   sizeof(coef_jnu8_a)/sizeof(double),
 876   sizeof(coef_jnu9_a)/sizeof(double),
 877   sizeof(coef_jnu10_a)/sizeof(double),
 878   sizeof(coef_jnu11_a)/sizeof(double),
 879   sizeof(coef_jnu12_a)/sizeof(double),
 880   sizeof(coef_jnu13_a)/sizeof(double),
 881   sizeof(coef_jnu14_a)/sizeof(double),
 882   sizeof(coef_jnu15_a)/sizeof(double),
 883   sizeof(coef_jnu16_a)/sizeof(double),
 884   sizeof(coef_jnu17_a)/sizeof(double),
 885   sizeof(coef_jnu18_a)/sizeof(double),
 886   sizeof(coef_jnu19_a)/sizeof(double),
 887   sizeof(coef_jnu20_a)/sizeof(double)
 888 };
 889
 890
 891 static const double * coef_jnu_b[] = {
 892   0,
 893   coef_jnu1_b,
 894   coef_jnu2_b,
 895   coef_jnu3_b,
 896   coef_jnu4_b,
 897   coef_jnu5_b,
 898   coef_jnu6_b,
 899   coef_jnu7_b,
 900   coef_jnu8_b,
 901   coef_jnu9_b,
 902   coef_jnu10_b
 903 };
 904
 905 static const size_t size_jnu_b[] = {
 906   0,
 907   sizeof(coef_jnu1_b)/sizeof(double),
 908   sizeof(coef_jnu2_b)/sizeof(double),
 909   sizeof(coef_jnu3_b)/sizeof(double),
 910   sizeof(coef_jnu4_b)/sizeof(double),
 911   sizeof(coef_jnu5_b)/sizeof(double),
 912   sizeof(coef_jnu6_b)/sizeof(double),
 913   sizeof(coef_jnu7_b)/sizeof(double),
 914   sizeof(coef_jnu8_b)/sizeof(double),
 915   sizeof(coef_jnu9_b)/sizeof(double),
 916   sizeof(coef_jnu10_b)/sizeof(double)
 917 };
 918
 919
 920
 921 /* Evaluate Clenshaw recurrence for
 922  * a T* Chebyshev series.
 923  * sizeof(c) = N+1
 924  */
 925 static double
 926 clenshaw(const double * c, int N, double u)
 927 {
 928   double B_np1 = 0.0;
 929   double B_n   = c[N];
 930   double B_nm1;
 931   int n;
 932   for(n=N; n>0; n--) {
 933     B_nm1 = 2.0*(2.0*u-1.0) * B_n - B_np1 + c[n-1];
 934     B_np1 = B_n;
 935     B_n   = B_nm1;
 936   }
 937   return B_n - (2.0*u-1.0)*B_np1;
 938 }
 939
 940
 941
 942 /* correction terms to leading McMahon expansion
 943  * [Abramowitz+Stegun 9.5.12]
 944  * [Olver, Royal Society Math. Tables, v. 7]
 945  * We factor out a beta, so that this is a multiplicative
 946  * correction:
 947  *   j_{nu,s} = beta(s,nu) * mcmahon_correction(nu, beta(s,nu))
 948  *   macmahon_correction --> 1 as s --> Inf
 949  */
 950 static double
 951 mcmahon_correction(const double mu, const double beta)
 952 {
 953   const double eb   = 8.0*beta;
 954   const double ebsq = eb*eb;
 955
 956   if(mu < GSL_DBL_EPSILON) {
 957     /* Prevent division by zero below. */
 958     const double term1 =  1.0/ebsq;
 959     const double term2 = -4.0*31.0/(3*ebsq*ebsq);
 960     const double term3 =  32.0*3779.0/(15.0*ebsq*ebsq*ebsq);
 961     const double term4 = -64.0*6277237.0/(105.0*ebsq*ebsq*ebsq*ebsq);
 962     const double term5 =  512.0*2092163573.0/(315.0*ebsq*ebsq*ebsq*ebsq*ebsq);
 963     return 1.0 + 8.0*(term1 + term2 + term3 + term4 + term5);
 964   }
 965   else {
 966     /* Here we do things in terms of 1/mu, which
 967      * is purely to prevent overflow in the very
 968      * unlikely case that mu is really big.
 969      */
 970     const double mi   = 1.0/mu;
 971     const double r  = mu/ebsq;
 972     const double n2 = 4.0/3.0    * (7.0 - 31.0*mi);
 973     const double n3 = 32.0/15.0  * (83.0 + (-982.0 + 3779.0*mi)*mi);
 974     const double n4 = 64.0/105.0 * (6949.0 + (-153855.0 + (1585743.0 - 6277237.0*mi)*mi)*mi);
 975     const double n5 = 512.0/315.0 * (70197.0 + (-2479316.0 + (48010494.0 + (-512062548.0 + 2092163573.0*mi)*mi)*mi)*mi);
 976     const double n6 = 2048.0/3465.0 * (5592657.0 + (-287149133.0 + (8903961290.0 + (-179289628602.0 + (1982611456181.0 - 8249725736393.0*mi)*mi)*mi)*mi)*mi);
 977     const double term1 = (1.0 - mi) * r;
 978     const double term2 = term1 * n2 * r;
 979     const double term3 = term1 * n3 * r*r;
 980     const double term4 = term1 * n4 * r*r*r;
 981     const double term5 = term1 * n5 * r*r*r*r;
 982     const double term6 = term1 * n6 * r*r*r*r*r;
 983     return 1.0 - 8.0*(term1 + term2 + term3 + term4 + term5 + term6);
 984   }
 985 }
 986
 987
 988 /* Assumes z >= 1.0 */
 989 static double
 990 olver_b0(double z, double minus_zeta)
 991 {
 992   if(z < 1.02) {
 993     const double a = 1.0-z;
 994     const double c0 =  0.0179988721413553309252458658183;
 995     const double c1 =  0.0111992982212877614645974276203;
 996     const double c2 =  0.0059404069786014304317781160605;
 997     const double c3 =  0.0028676724516390040844556450173;
 998     const double c4 =  0.0012339189052567271708525111185;
 999     const double c5 =  0.0004169250674535178764734660248;
1000     const double c6 =  0.0000330173385085949806952777365;
1001     const double c7 = -0.0001318076238578203009990106425;
1002     const double c8 = -0.0001906870370050847239813945647;
1003     return c0 + a*(c1 + a*(c2 + a*(c3 + a*(c4 + a*(c5 + a*(c6 + a*(c7 + a*c8)))))));
1004   }
1005   else {
1006     const double abs_zeta = minus_zeta;
1007     const double t = 1.0/(z*sqrt(1.0 - 1.0/(z*z)));
1008     return -5.0/(48.0*abs_zeta*abs_zeta) + t*(3.0 + 5.0*t*t)/(24.0*sqrt(abs_zeta));
1009   }
1010 }
1011
1012
1013 inline
1014 static double
1015 olver_f1(double z, double minus_zeta)
1016 {
1017   const double b0 = olver_b0(z, minus_zeta);
1018   const double h2 = sqrt(4.0*minus_zeta/(z*z-1.0)); /* FIXME */
1019   return 0.5 * z * h2 * b0;
1020 }
1021
1022
1023 int
1024 gsl_sf_bessel_zero_J0_e(unsigned int s, gsl_sf_result * result)
1025 {
1026   /* CHECK_POINTER(result) */
1027
1028   if(s == 0){
1029     result->val = 0.0;
1030     result->err = 0.0;
1031     GSL_ERROR ("error", GSL_EINVAL);
1032   }
1033   else {
1034     /* See [F. Lether, J. Comp. Appl .Math. 67, 167 (1996)]. */
1035
1036     static const double P[] = { 1567450796.0/12539606369.0,
1037                                 8903660.0/2365861.0,
1038                                 10747040.0/536751.0,
1039                                 17590991.0/1696654.0
1040                               };
1041     static const double Q[] = { 1.0,
1042                                 29354255.0/954518.0,
1043                                 76900001.0/431847.0,
1044                                 67237052.0/442411.0
1045                               };
1046
1047     const double beta = (s - 0.25) * M_PI;
1048     const double bi2  = 1.0/(beta*beta);
1049     const double R33num = P[0] + bi2 * (P[1] + bi2 * (P[2] + P[3] * bi2));
1050     const double R33den = Q[0] + bi2 * (Q[1] + bi2 * (Q[2] + Q[3] * bi2));
1051     const double R33 = R33num/R33den;
1052     result->val = beta + R33/beta;
1053     result->err = fabs(3.0e-15 * result->val);
1054     return GSL_SUCCESS;
1055   }
1056 }
1057
1058
1059 int
1060 gsl_sf_bessel_zero_J1_e(unsigned int s, gsl_sf_result * result)
1061 {
1062   /* CHECK_POINTER(result) */
1063
1064   if(s == 0) {
1065     result->val = 0.0;
1066     result->err = 0.0;
1067     return GSL_SUCCESS;
1068   }
1069   else {
1070     /* See [M. Branders et al., J. Comp. Phys. 42, 403 (1981)]. */
1071
1072     static const double a[] = { -0.362804405737084,
1073                                  0.120341279038597,
1074                                  0.439454547101171e-01,
1075                                  0.159340088474713e-02
1076                               };
1077     static const double b[] = {  1.0,
1078                                 -0.325641790801361,
1079                                 -0.117453445968927,
1080                                 -0.424906902601794e-02
1081                               };
1082
1083     const double beta = (s + 0.25) * M_PI;
1084     const double bi2  = 1.0/(beta*beta);
1085     const double Rnum = a[3] + bi2 * (a[2] + bi2 * (a[1] + bi2 * a[0]));
1086     const double Rden = b[3] + bi2 * (b[2] + bi2 * (b[1] + bi2 * b[0]));
1087     const double R = Rnum/Rden;
1088     result->val = beta * (1.0 + R*bi2);
1089     result->err = fabs(2.0e-14 * result->val);
1090     return GSL_SUCCESS;
1091   }
1092 }
1093
1094
1095 int
1096 gsl_sf_bessel_zero_Jnu_e(double nu, unsigned int s, gsl_sf_result * result)
1097 {
1098   /* CHECK_POINTER(result) */
1099
1100   if(nu <= -1.0) {
1101     DOMAIN_ERROR(result);
1102   }
1103   else if(s == 0) {
1104     result->val = 0.0;
1105     result->err = 0.0;
1106     if (nu == 0.0) {
1107       GSL_ERROR ("no zero-th root for nu = 0.0", GSL_EINVAL);
1108     }
1109     return GSL_SUCCESS;
1110   }
1111   else if(nu < 0.0) {
1112     /* This can be done, I'm just lazy now. */
1113     result->val = 0.0;
1114     result->err = 0.0;
1115     GSL_ERROR("unimplemented", GSL_EUNIMPL);
1116   }
1117   else if(s == 1) {
1118     /* Chebyshev fits for the first positive zero.
1119      * For some reason Nemeth made this different from the others.
1120      */
1121     if(nu < 2.0) {
1122       const double * c = coef_jnu_a[s];
1123       const size_t   L = size_jnu_a[s];
1124       const double arg = nu/2.0;
1125       const double chb = clenshaw(c, L-1, arg);
1126       result->val = chb;
1127       result->err = 2.0e-15 * result->val;
1128     }
1129     else {
1130       const double * c = coef_jnu_b[s];
1131       const size_t   L = size_jnu_b[s];
1132       const double arg = pow(2.0/nu, 2.0/3.0);
1133       const double chb = clenshaw(c, L-1, arg);
1134       result->val = nu * chb;
1135       result->err = 2.0e-15 * result->val;
1136     }
1137     return GSL_SUCCESS;
1138   }
1139   else if(s <= 10) {
1140     /* Chebyshev fits for the first 10 positive zeros. */
1141     if(nu < s) {
1142       const double * c = coef_jnu_a[s];
1143       const size_t   L = size_jnu_a[s];
1144       const double arg = nu/s;
1145       const double chb = clenshaw(c, L-1, arg);
1146       result->val = chb;
1147       result->err = 2.0e-15 * result->val;
1148     }
1149     else {
1150       const double * c = coef_jnu_b[s];
1151       const size_t   L = size_jnu_b[s];
1152       const double arg = pow(s/nu, 2.0/3.0);
1153       const double chb = clenshaw(c, L-1, arg);
1154       result->val = nu * chb;
1155       result->err = 2.0e-15 * result->val;
1156
1157       /* FIXME: truth in advertising for the screwed up
1158        * s = 5 fit. Need to fix that.
1159        */
1160       if(s == 5) {
1161         result->err *= 5.0e+06;
1162       }
1163     }
1164     return GSL_SUCCESS;
1165   }
1166   else if(s > 0.5*nu && s <= 20) {
1167     /* Chebyshev fits for 10 < s <= 20. */
1168     const double * c = coef_jnu_a[s];
1169     const size_t   L = size_jnu_a[s];
1170     const double arg = nu/(2.0*s);
1171     const double chb = clenshaw(c, L-1, arg);
1172     result->val = chb;
1173     result->err = 4.0e-15 * chb;
1174     return GSL_SUCCESS;
1175   }
1176   else if(s > 2.0 * nu) {
1177     /* McMahon expansion if s is large compared to nu. */
1178     const double beta = (s + 0.5*nu - 0.25) * M_PI;
1179     const double mc   = mcmahon_correction(4.0*nu*nu, beta);
1180     gsl_sf_result rat12;
1181     gsl_sf_pow_int_e(nu/beta, 14, &rat12);
1182     result->val  = beta * mc;
1183     result->err  = 4.0 * fabs(beta) * rat12.val;
1184     result->err += 4.0 * fabs(GSL_DBL_EPSILON * result->val);
1185     return GSL_SUCCESS;
1186   }
1187   else {
1188     /* Olver uniform asymptotic. */
1189     gsl_sf_result as;
1190     const int stat_as = gsl_sf_airy_zero_Ai_e(s, &as);
1191     const double minus_zeta = -pow(nu,-2.0/3.0) * as.val;
1192     const double z  = gsl_sf_bessel_Olver_zofmzeta(minus_zeta);
1193     const double f1 = olver_f1(z, minus_zeta);
1194     result->val  = nu * (z + f1/(nu*nu));
1195     result->err  = 0.001/(nu*nu*nu);
1196     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
1197     return stat_as;
1198   }
1199 }
1200
1201
1202 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
1203
1204 #include "eval.h"
1205
1206 double gsl_sf_bessel_zero_J0(unsigned int s)
1207 {
1208   EVAL_RESULT(gsl_sf_bessel_zero_J0_e(s, &result));
1209 }
1210
1211 double gsl_sf_bessel_zero_J1(unsigned int s)
1212 {
1213   EVAL_RESULT(gsl_sf_bessel_zero_J1_e(s, &result));
1214 }
1215
1216 double gsl_sf_bessel_zero_Jnu(double nu, unsigned int s)
1217 {
1218   EVAL_RESULT(gsl_sf_bessel_zero_Jnu_e(nu, s, &result));
1219 }