diff options
Diffstat (limited to 'gsl-1.9/specfunc/bessel_zero.c')
-rw-r--r-- | gsl-1.9/specfunc/bessel_zero.c | 1219 |
1 files changed, 1219 insertions, 0 deletions
diff --git a/gsl-1.9/specfunc/bessel_zero.c b/gsl-1.9/specfunc/bessel_zero.c new file mode 100644 index 0000000..7ff8f7e --- /dev/null +++ b/gsl-1.9/specfunc/bessel_zero.c @@ -0,0 +1,1219 @@ +/* specfunc/bessel_zero.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* Author: G. Jungman */ + +#include <config.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_sf_airy.h> +#include <gsl/gsl_sf_pow_int.h> +#include <gsl/gsl_sf_bessel.h> + +#include "error.h" + +#include "bessel_olver.h" + +/* For Chebyshev expansions of the roots as functions of nu, + * see [G. Nemeth, Mathematical Approximation of Special Functions]. + * This gives the fits for all nu and s <= 10. + * I made the fits for other values of s myself [GJ]. + */ + +/* Chebyshev expansion: j_{nu,1} = c_k T_k*(nu/2), nu <= 2 */ +static const double coef_jnu1_a[] = { + 3.801775243633476, + 1.360704737511120, + -0.030707710261106, + 0.004526823746202, + -0.000808682832134, + 0.000159218792489, + -0.000033225189761, + 0.000007205599763, + -0.000001606110397, + 0.000000365439424, + -0.000000084498039, + 0.000000019793815, + -0.000000004687054, + 0.000000001120052, + -0.000000000269767, + 0.000000000065420, + -0.000000000015961, + 0.000000000003914, + -0.000000000000965, + 0.000000000000239, + -0.000000000000059, + 0.000000000000015, + -0.000000000000004, + 0.000000000000001 +}; + + +/* Chebyshev expansion: j_{nu,1} = nu c_k T_k*((2/nu)^(2/3)), nu >= 2 */ +static const double coef_jnu1_b[] = { + 1.735063412537096, + 0.784478100951978, + 0.048881473180370, + -0.000578279783021, + -0.000038984957864, + 0.000005758297879, + -0.000000327583229, + -0.000000003853878, + 0.000000002284653, + -0.000000000153079, + -0.000000000000895, + 0.000000000000283, + 0.000000000000043, + 0.000000000000010, + -0.000000000000003 +}; + + +/* Chebyshev expansion: j_{nu,2} = c_k T_k*(nu/2), nu <= 2 */ +static const double coef_jnu2_a[] = { + 6.992370244046161, + 1.446379282056534, + -0.023458616207293, + 0.002172149448700, + -0.000246262775620, + 0.000030990180959, + -0.000004154183047, + 0.000000580766328, + -0.000000083648175, + 0.000000012317355, + -0.000000001844887, + 0.000000000280076, + -0.000000000042986, + 0.000000000006658, + -0.000000000001039, + 0.000000000000163, + -0.000000000000026, + 0.000000000000004, + -0.000000000000001 +}; + + +/* Chebyshev expansion: j_{nu,2} = nu c_k T_k*((2/nu)^(2/3)), nu >= 2 */ +static const double coef_jnu2_b[] = { + 2.465611864263400, + 1.607952988471069, + 0.138758034431497, + -0.003687791182054, + -0.000051276007868, + 0.000045113570749, + -0.000007579172152, + 0.000000736469208, + -0.000000011118527, + -0.000000011919884, + 0.000000002696788, + -0.000000000314488, + 0.000000000008124, + 0.000000000005211, + -0.000000000001292, + 0.000000000000158, + -0.000000000000004, + -0.000000000000003, + 0.000000000000001 +}; + + +/* Chebyshev expansion: j_{nu,3} = c_k T_k*(nu/3), nu <= 3 */ +static const double coef_jnu3_a[] = { + 10.869647065239236, + 2.177524286141710, + -0.034822817125293, + 0.003167249102413, + -0.000353960349344, + 0.000044039086085, + -0.000005851380981, + 0.000000812575483, + -0.000000116463617, + 0.000000017091246, + -0.000000002554376, + 0.000000000387335, + -0.000000000059428, + 0.000000000009207, + -0.000000000001438, + 0.000000000000226, + -0.000000000000036, + 0.000000000000006, + -0.000000000000001 +}; + + +/* Chebyshev expansion: j_{nu,3} = nu c_k T_k*((3/nu)^(2/3)), nu >= 3 */ +static const double coef_jnu3_b[] = { + 2.522816775173244, + 1.673199424973720, + 0.146431617506314, + -0.004049001763912, + -0.000039517767244, + 0.000048781729288, + -0.000008729705695, + 0.000000928737310, + -0.000000028388244, + -0.000000012927432, + 0.000000003441008, + -0.000000000471695, + 0.000000000025590, + 0.000000000005502, + -0.000000000001881, + 0.000000000000295, + -0.000000000000020, + -0.000000000000003, + 0.000000000000001 +}; + + +/* Chebyshev expansion: j_{nu,4} = c_k T_k*(nu/4), nu <= 4 */ +static const double coef_jnu4_a[] = { + 14.750310252773009, + 2.908010932941708, + -0.046093293420315, + 0.004147172321412, + -0.000459092310473, + 0.000056646951906, + -0.000007472351546, + 0.000001031210065, + -0.000000147008137, + 0.000000021475218, + -0.000000003197208, + 0.000000000483249, + -0.000000000073946, + 0.000000000011431, + -0.000000000001782, + 0.000000000000280, + -0.000000000000044, + 0.000000000000007, + -0.000000000000001 +}; + + +/* Chebyshev expansion: j_{nu,4} = nu c_k T_k*((4/nu)^(2/3)), nu >= 4 */ +static const double coef_jnu4_b[] = { + 2.551681323117914, + 1.706177978336572, + 0.150357658406131, + -0.004234001378590, + -0.000033854229898, + 0.000050763551485, + -0.000009337464057, + 0.000001029717834, + -0.000000037474196, + -0.000000013450153, + 0.000000003836180, + -0.000000000557404, + 0.000000000035748, + 0.000000000005487, + -0.000000000002187, + 0.000000000000374, + -0.000000000000031, + -0.000000000000003, + 0.000000000000001 +}; + + + +/* Chebyshev expansion: j_{nu,5} = c_k T_k*(nu/5), nu <= 5 */ +static const double coef_jnu5_a[] = { + 18.632261081028211, + 3.638249012596966, + -0.057329705998828, + 0.005121709126820, + -0.000563325259487, + 0.000069100826174, + -0.000009066603030, + 0.000001245181383, + -0.000000176737282, + 0.000000025716695, + -0.000000003815184, + 0.000000000574839, + -0.000000000087715, + 0.000000000013526, + -0.000000000002104, + 0.000000000000330, + -0.000000000000052, + 0.000000000000008, + -0.000000000000001 +}; + + +/* Chebyshev expansion: j_{nu,5} = nu c_k T_k*((5/nu)^(2/3)), nu >= 5 */ +/* FIXME: There is something wrong with this fit, in about the + * 9th or 10th decimal place. + */ +static const double coef_jnu5_b[] = { + 2.569079487591442, + 1.726073360882134, + 0.152740776809531, + -0.004346449660148, + -0.000030512461856, + 0.000052000821080, + -0.000009713343981, + 0.000001091997863, + -0.000000043061707, + -0.000000013779413, + 0.000000004082870, + -0.000000000611259, + 0.000000000042242, + 0.000000000005448, + -0.000000000002377, + 0.000000000000424, + -0.000000000000038, + -0.000000000000002, + 0.000000000000002 +}; + + +/* Chebyshev expansion: j_{nu,6} = c_k T_k*(nu/6), nu <= 6 */ +static const double coef_jnu6_a[] = { + 22.514836143374042, + 4.368367257557198, + -0.068550155285562, + 0.006093776505822, + -0.000667152784957, + 0.000081486022398, + -0.000010649011647, + 0.000001457089679, + -0.000000206105082, + 0.000000029894724, + -0.000000004422012, + 0.000000000664471, + -0.000000000101140, + 0.000000000015561, + -0.000000000002416, + 0.000000000000378, + -0.000000000000060, + 0.000000000000009, + -0.000000000000002 +}; + + +/* Chebyshev expansion: j_{nu,6} = nu c_k T_k*((6/nu)^(2/3)), nu >= 6 */ +static const double coef_jnu6_b[] = { + 2.580710285494837, + 1.739380728566154, + 0.154340696401691, + -0.004422028860168, + -0.000028305272624, + 0.000052845975269, + -0.000009968794373, + 0.000001134252926, + -0.000000046841241, + -0.000000014007555, + 0.000000004251816, + -0.000000000648213, + 0.000000000046728, + 0.000000000005414, + -0.000000000002508, + 0.000000000000459, + -0.000000000000043, + -0.000000000000002, + 0.000000000000002 +}; + + +/* Chebyshev expansion: j_{nu,7} = c_k T_k*(nu/7), nu <= 7 */ +static const double coef_jnu7_a[] = { + 26.397760539730869, + 5.098418721711790, + -0.079761896398948, + 0.007064521280487, + -0.000770766522482, + 0.000093835449636, + -0.000012225308542, + 0.000001667939800, + -0.000000235288157, + 0.000000034040347, + -0.000000005023142, + 0.000000000753101, + -0.000000000114389, + 0.000000000017564, + -0.000000000002722, + 0.000000000000425, + -0.000000000000067, + 0.000000000000011, + -0.000000000000002 +}; + + +/* Chebyshev expansion: j_{nu,7} = nu c_k T_k*((7/nu)^(2/3)), nu >= 7 */ +static const double coef_jnu7_b[] = { + 2.589033335856773, + 1.748907007612678, + 0.155488900387653, + -0.004476317805688, + -0.000026737952924, + 0.000053459680946, + -0.000010153699240, + 0.000001164804272, + -0.000000049566917, + -0.000000014175403, + 0.000000004374840, + -0.000000000675135, + 0.000000000050004, + 0.000000000005387, + -0.000000000002603, + 0.000000000000485, + -0.000000000000047, + -0.000000000000002, + 0.000000000000002 +}; + + +/* Chebyshev expansion: j_{nu,8} = c_k T_k*(nu/8), nu <= 8 */ +static const double coef_jnu8_a[] = { + 30.280900001606662, + 5.828429205461221, + -0.090968381181069, + 0.008034479731033, + -0.000874254899080, + 0.000106164151611, + -0.000013798098749, + 0.000001878187386, + -0.000000264366627, + 0.000000038167685, + -0.000000005621060, + 0.000000000841165, + -0.000000000127538, + 0.000000000019550, + -0.000000000003025, + 0.000000000000472, + -0.000000000000074, + 0.000000000000012, + -0.000000000000002 +}; + + +/* Chebyshev expansion: j_{nu,8} = nu c_k T_k*((8/nu)^(2/3)), nu >= 8 */ +static const double coef_jnu8_b[] = { + 2.595283877150078, + 1.756063044986928, + 0.156352972371030, + -0.004517201896761, + -0.000025567187878, + 0.000053925472558, + -0.000010293734486, + 0.000001187923085, + -0.000000051625122, + -0.000000014304212, + 0.000000004468450, + -0.000000000695620, + 0.000000000052500, + 0.000000000005367, + -0.000000000002676, + 0.000000000000505, + -0.000000000000050, + -0.000000000000002, + 0.000000000000002 +}; + + +/* Chebyshev expansion: j_{nu,9} = c_k T_k*(nu/9), nu <= 9 */ +static const double coef_jnu9_a[] = { + 34.164181213238386, + 6.558412747925228, + -0.102171455365016, + 0.009003934361201, + -0.000977663914535, + 0.000118479876579, + -0.000015368714220, + 0.000002088064285, + -0.000000293381154, + 0.000000042283900, + -0.000000006217033, + 0.000000000928887, + -0.000000000140627, + 0.000000000021526, + -0.000000000003326, + 0.000000000000518, + -0.000000000000081, + 0.000000000000013, + -0.000000000000002 +}; + + +/* Chebyshev expansion: j_{nu,9} = nu c_k T_k*((9/nu)^(2/3)), nu >= 9 */ +static const double coef_jnu9_b[] = { + 2.600150240905079, + 1.761635491694032, + 0.157026743724010, + -0.004549100368716, + -0.000024659248617, + 0.000054291035068, + -0.000010403464334, + 0.000001206027524, + -0.000000053234089, + -0.000000014406241, + 0.000000004542078, + -0.000000000711728, + 0.000000000054464, + 0.000000000005350, + -0.000000000002733, + 0.000000000000521, + -0.000000000000052, + -0.000000000000002, + 0.000000000000002 +}; + + +/* Chebyshev expansion: j_{nu,10} = c_k T_k*(nu/10), nu <= 10 */ +static const double coef_jnu10_a[] = { + 38.047560766184647, + 7.288377637926008, + -0.113372193277897, + 0.009973047509098, + -0.001081019701335, + 0.000130786983847, + -0.000016937898538, + 0.000002297699179, + -0.000000322354218, + 0.000000046392941, + -0.000000006811759, + 0.000000001016395, + -0.000000000153677, + 0.000000000023486, + -0.000000000003616, + 0.000000000000561, + -0.000000000000095, + 0.000000000000027, + -0.000000000000013, + 0.000000000000005 +}; + + +/* Chebyshev expansion: j_{nu,10} = nu c_k T_k*((10/nu)^(2/3)), nu >= 10 */ +static const double coef_jnu10_b[] = { + 2.604046346867949, + 1.766097596481182, + 0.157566834446511, + -0.004574682244089, + -0.000023934500688, + 0.000054585558231, + -0.000010491765415, + 0.000001220589364, + -0.000000054526331, + -0.000000014489078, + 0.000000004601510, + -0.000000000724727, + 0.000000000056049, + 0.000000000005337, + -0.000000000002779, + 0.000000000000533, + -0.000000000000054, + -0.000000000000002, + 0.000000000000002 +}; + + +/* Chebyshev expansion: j_{nu,11} = c_k T_k*(nu/22), nu <= 22 */ +static const double coef_jnu11_a[] = { + 49.5054081076848637, + 15.33692279367165101, + -0.33677234163517130, + 0.04623235772920729, + -0.00781084960665093, + 0.00147217395434708, + -0.00029695043846867, + 0.00006273356860235, + -0.00001370575125628, + 3.07171282012e-6, + -7.0235041249e-7, + 1.6320559339e-7, + -3.843117306e-8, + 9.15083800e-9, + -2.19957642e-9, + 5.3301703e-10, + -1.3007541e-10, + 3.193827e-11, + -7.88605e-12, + 1.95918e-12, + -4.9020e-13, + 1.2207e-13, + -2.820e-14, + 5.25e-15, + -1.88e-15, + 2.80e-15, + -2.45e-15 +}; + + +/* Chebyshev expansion: j_{nu,12} = c_k T_k*(nu/24), nu <= 24 */ +static const double coef_jnu12_a[] = { + 54.0787833216641519, + 16.7336367772863598, + -0.36718411124537953, + 0.05035523375053820, + -0.00849884978867533, + 0.00160027692813434, + -0.00032248114889921, + 0.00006806354127199, + -0.00001485665901339, + 3.32668783672e-6, + -7.5998952729e-7, + 1.7644939709e-7, + -4.151538210e-8, + 9.87722772e-9, + -2.37230133e-9, + 5.7442875e-10, + -1.4007767e-10, + 3.437166e-11, + -8.48215e-12, + 2.10554e-12, + -5.2623e-13, + 1.3189e-13, + -3.175e-14, + 5.73e-15, + 5.6e-16, + -8.7e-16, + -6.5e-16 +}; + + +/* Chebyshev expansion: j_{nu,13} = c_k T_k*(nu/26), nu <= 26 */ +static const double coef_jnu13_a[] = { + 58.6521941921708890, + 18.1303398137970284, + -0.39759381380126650, + 0.05447765240465494, + -0.00918674227679980, + 0.00172835361420579, + -0.00034800528297612, + 0.00007339183835188, + -0.00001600713368099, + 3.58154960392e-6, + -8.1759873497e-7, + 1.8968523220e-7, + -4.459745253e-8, + 1.060304419e-8, + -2.54487624e-9, + 6.1580214e-10, + -1.5006751e-10, + 3.679707e-11, + -9.07159e-12, + 2.24713e-12, + -5.5943e-13, + 1.4069e-13, + -3.679e-14, + 1.119e-14, + -4.99e-15, + 3.43e-15, + -2.85e-15, + 2.3e-15, + -1.7e-15, + 8.7e-16 +}; + + +/* Chebyshev expansion: j_{nu,14} = c_k T_k*(nu/28), nu <= 28 */ +static const double coef_jnu14_a[] = { + 63.2256329577315566, + 19.5270342832914901, + -0.42800190567884337, + 0.05859971627729398, + -0.00987455163523582, + 0.00185641011402081, + -0.00037352439419968, + 0.00007871886257265, + -0.00001715728110045, + 3.83632624437e-6, + -8.7518558668e-7, + 2.0291515353e-7, + -4.767795233e-8, + 1.132844415e-8, + -2.71734219e-9, + 6.5714886e-10, + -1.6005342e-10, + 3.922557e-11, + -9.66637e-12, + 2.39379e-12, + -5.9541e-13, + 1.4868e-13, + -3.726e-14, + 9.37e-15, + -2.36e-15, + 6.0e-16 +}; + + +/* Chebyshev expansion: j_{nu,15} = c_k T_k*(nu/30), nu <= 30 */ +static const double coef_jnu15_a[] = { + 67.7990939565631635, + 20.9237219226859859, + -0.45840871823085836, + 0.06272149946755639, + -0.01056229551143042, + 0.00198445078693100, + -0.00039903958650729, + 0.00008404489865469, + -0.00001830717574922, + 4.09103745566e-6, + -9.3275533309e-7, + 2.1614056403e-7, + -5.075725222e-8, + 1.205352081e-8, + -2.88971837e-9, + 6.9846848e-10, + -1.7002946e-10, + 4.164941e-11, + -1.025859e-11, + 2.53921e-12, + -6.3128e-13, + 1.5757e-13, + -3.947e-14, + 9.92e-15, + -2.50e-15, + 6.3e-16 +}; + + +/* Chebyshev expansion: j_{nu,16} = c_k T_k*(nu/32), nu <= 32 */ +static const double coef_jnu16_a[] = { + 72.3725729616724770, + 22.32040402918608585, + -0.48881449782358690, + 0.06684305681828766, + -0.01124998690363398, + 0.00211247882775445, + -0.00042455166484632, + 0.00008937015316346, + -0.00001945687139551, + 4.34569739281e-6, + -9.9031173548e-7, + 2.2936247195e-7, + -5.383562595e-8, + 1.277835103e-8, + -3.06202860e-9, + 7.3977037e-10, + -1.8000071e-10, + 4.407196e-11, + -1.085046e-11, + 2.68453e-12, + -6.6712e-13, + 1.6644e-13, + -4.168e-14, + 1.047e-14, + -2.64e-15, + 6.7e-16 +}; + + +/* Chebyshev expansion: j_{nu,17} = c_k T_k*(nu/34), nu <= 34 */ +static const double coef_jnu17_a[] = { + 76.9460667535209549, + 23.71708159112252670, + -0.51921943142405352, + 0.07096442978067622, + -0.01193763559341369, + 0.00224049662974902, + -0.00045006122941781, + 0.00009469477941684, + -0.00002060640777107, + 4.60031647195e-6, + -1.04785755046e-6, + 2.4258161247e-7, + -5.691327087e-8, + 1.350298805e-8, + -3.23428733e-9, + 7.8105847e-10, + -1.8996825e-10, + 4.649350e-11, + -1.144205e-11, + 2.82979e-12, + -7.0294e-13, + 1.7531e-13, + -4.388e-14, + 1.102e-14, + -2.78e-15, + 7.0e-16 +}; + + +/* Chebyshev expansion: j_{nu,18} = c_k T_k*(nu/36), nu <= 36 */ +static const double coef_jnu18_a[] = { + 81.5195728368096659, + 25.11375537470259305, + -0.54962366347317668, + 0.07508565026117689, + -0.01262524908033818, + 0.00236850602019778, + -0.00047556873651929, + 0.00010001889347161, + -0.00002175581482429, + 4.85490251239e-6, + -1.10539483940e-6, + 2.5579853343e-7, + -5.999033352e-8, + 1.422747129e-8, + -3.40650521e-9, + 8.2233565e-10, + -1.9993286e-10, + 4.891426e-11, + -1.203343e-11, + 2.97498e-12, + -7.3875e-13, + 1.8418e-13, + -4.608e-14, + 1.157e-14, + -2.91e-15, + 7.4e-16 +}; + + +/* Chebyshev expansion: j_{nu,19} = c_k T_k*(nu/38), nu <= 38 */ +static const double coef_jnu19_a[] = { + 86.0930892477047512, + 26.51042598308271729, + -0.58002730731948358, + 0.07920674321589394, + -0.01331283320930301, + 0.00249650841778073, + -0.00050107453900793, + 0.00010534258471335, + -0.00002290511552874, + 5.10946148897e-6, + -1.16292517157e-6, + 2.6901365037e-7, + -6.306692473e-8, + 1.495183048e-8, + -3.57869025e-9, + 8.6360410e-10, + -2.0989514e-10, + 5.133439e-11, + -1.262465e-11, + 3.12013e-12, + -7.7455e-13, + 1.9304e-13, + -4.829e-14, + 1.212e-14, + -3.05e-15, + 7.7e-16 +}; + + +/* Chebyshev expansion: j_{nu,20} = c_k T_k*(nu/40), nu <= 40 */ +static const double coef_jnu20_a[] = { + 90.6666144195163770, + 27.9070938975436823, + -0.61043045315390591, + 0.08332772844325554, + -0.01400039260208282, + 0.00262450494035660, + -0.00052657891389470, + 0.00011066592304919, + -0.00002405432778364, + 5.36399803946e-6, + -1.22044976064e-6, + 2.8222728362e-7, + -6.614312964e-8, + 1.567608839e-8, + -3.75084856e-9, + 9.0486546e-10, + -2.1985553e-10, + 5.375401e-11, + -1.321572e-11, + 3.26524e-12, + -8.1033e-13, + 2.0190e-13, + -5.049e-14, + 1.267e-14, + -3.19e-15, + 8.0e-16, + -2.0e-16 +}; + + +static const double * coef_jnu_a[] = { + 0, + coef_jnu1_a, + coef_jnu2_a, + coef_jnu3_a, + coef_jnu4_a, + coef_jnu5_a, + coef_jnu6_a, + coef_jnu7_a, + coef_jnu8_a, + coef_jnu9_a, + coef_jnu10_a, + coef_jnu11_a, + coef_jnu12_a, + coef_jnu13_a, + coef_jnu14_a, + coef_jnu15_a, + coef_jnu16_a, + coef_jnu17_a, + coef_jnu18_a, + coef_jnu19_a, + coef_jnu20_a +}; + +static const size_t size_jnu_a[] = { + 0, + sizeof(coef_jnu1_a)/sizeof(double), + sizeof(coef_jnu2_a)/sizeof(double), + sizeof(coef_jnu3_a)/sizeof(double), + sizeof(coef_jnu4_a)/sizeof(double), + sizeof(coef_jnu5_a)/sizeof(double), + sizeof(coef_jnu6_a)/sizeof(double), + sizeof(coef_jnu7_a)/sizeof(double), + sizeof(coef_jnu8_a)/sizeof(double), + sizeof(coef_jnu9_a)/sizeof(double), + sizeof(coef_jnu10_a)/sizeof(double), + sizeof(coef_jnu11_a)/sizeof(double), + sizeof(coef_jnu12_a)/sizeof(double), + sizeof(coef_jnu13_a)/sizeof(double), + sizeof(coef_jnu14_a)/sizeof(double), + sizeof(coef_jnu15_a)/sizeof(double), + sizeof(coef_jnu16_a)/sizeof(double), + sizeof(coef_jnu17_a)/sizeof(double), + sizeof(coef_jnu18_a)/sizeof(double), + sizeof(coef_jnu19_a)/sizeof(double), + sizeof(coef_jnu20_a)/sizeof(double) +}; + + +static const double * coef_jnu_b[] = { + 0, + coef_jnu1_b, + coef_jnu2_b, + coef_jnu3_b, + coef_jnu4_b, + coef_jnu5_b, + coef_jnu6_b, + coef_jnu7_b, + coef_jnu8_b, + coef_jnu9_b, + coef_jnu10_b +}; + +static const size_t size_jnu_b[] = { + 0, + sizeof(coef_jnu1_b)/sizeof(double), + sizeof(coef_jnu2_b)/sizeof(double), + sizeof(coef_jnu3_b)/sizeof(double), + sizeof(coef_jnu4_b)/sizeof(double), + sizeof(coef_jnu5_b)/sizeof(double), + sizeof(coef_jnu6_b)/sizeof(double), + sizeof(coef_jnu7_b)/sizeof(double), + sizeof(coef_jnu8_b)/sizeof(double), + sizeof(coef_jnu9_b)/sizeof(double), + sizeof(coef_jnu10_b)/sizeof(double) +}; + + + +/* Evaluate Clenshaw recurrence for + * a T* Chebyshev series. + * sizeof(c) = N+1 + */ +static double +clenshaw(const double * c, int N, double u) +{ + double B_np1 = 0.0; + double B_n = c[N]; + double B_nm1; + int n; + for(n=N; n>0; n--) { + B_nm1 = 2.0*(2.0*u-1.0) * B_n - B_np1 + c[n-1]; + B_np1 = B_n; + B_n = B_nm1; + } + return B_n - (2.0*u-1.0)*B_np1; +} + + + +/* correction terms to leading McMahon expansion + * [Abramowitz+Stegun 9.5.12] + * [Olver, Royal Society Math. Tables, v. 7] + * We factor out a beta, so that this is a multiplicative + * correction: + * j_{nu,s} = beta(s,nu) * mcmahon_correction(nu, beta(s,nu)) + * macmahon_correction --> 1 as s --> Inf + */ +static double +mcmahon_correction(const double mu, const double beta) +{ + const double eb = 8.0*beta; + const double ebsq = eb*eb; + + if(mu < GSL_DBL_EPSILON) { + /* Prevent division by zero below. */ + const double term1 = 1.0/ebsq; + const double term2 = -4.0*31.0/(3*ebsq*ebsq); + const double term3 = 32.0*3779.0/(15.0*ebsq*ebsq*ebsq); + const double term4 = -64.0*6277237.0/(105.0*ebsq*ebsq*ebsq*ebsq); + const double term5 = 512.0*2092163573.0/(315.0*ebsq*ebsq*ebsq*ebsq*ebsq); + return 1.0 + 8.0*(term1 + term2 + term3 + term4 + term5); + } + else { + /* Here we do things in terms of 1/mu, which + * is purely to prevent overflow in the very + * unlikely case that mu is really big. + */ + const double mi = 1.0/mu; + const double r = mu/ebsq; + const double n2 = 4.0/3.0 * (7.0 - 31.0*mi); + const double n3 = 32.0/15.0 * (83.0 + (-982.0 + 3779.0*mi)*mi); + const double n4 = 64.0/105.0 * (6949.0 + (-153855.0 + (1585743.0 - 6277237.0*mi)*mi)*mi); + const double n5 = 512.0/315.0 * (70197.0 + (-2479316.0 + (48010494.0 + (-512062548.0 + 2092163573.0*mi)*mi)*mi)*mi); + const double n6 = 2048.0/3465.0 * (5592657.0 + (-287149133.0 + (8903961290.0 + (-179289628602.0 + (1982611456181.0 - 8249725736393.0*mi)*mi)*mi)*mi)*mi); + const double term1 = (1.0 - mi) * r; + const double term2 = term1 * n2 * r; + const double term3 = term1 * n3 * r*r; + const double term4 = term1 * n4 * r*r*r; + const double term5 = term1 * n5 * r*r*r*r; + const double term6 = term1 * n6 * r*r*r*r*r; + return 1.0 - 8.0*(term1 + term2 + term3 + term4 + term5 + term6); + } +} + + +/* Assumes z >= 1.0 */ +static double +olver_b0(double z, double minus_zeta) +{ + if(z < 1.02) { + const double a = 1.0-z; + const double c0 = 0.0179988721413553309252458658183; + const double c1 = 0.0111992982212877614645974276203; + const double c2 = 0.0059404069786014304317781160605; + const double c3 = 0.0028676724516390040844556450173; + const double c4 = 0.0012339189052567271708525111185; + const double c5 = 0.0004169250674535178764734660248; + const double c6 = 0.0000330173385085949806952777365; + const double c7 = -0.0001318076238578203009990106425; + const double c8 = -0.0001906870370050847239813945647; + return c0 + a*(c1 + a*(c2 + a*(c3 + a*(c4 + a*(c5 + a*(c6 + a*(c7 + a*c8))))))); + } + else { + const double abs_zeta = minus_zeta; + const double t = 1.0/(z*sqrt(1.0 - 1.0/(z*z))); + return -5.0/(48.0*abs_zeta*abs_zeta) + t*(3.0 + 5.0*t*t)/(24.0*sqrt(abs_zeta)); + } +} + + +inline +static double +olver_f1(double z, double minus_zeta) +{ + const double b0 = olver_b0(z, minus_zeta); + const double h2 = sqrt(4.0*minus_zeta/(z*z-1.0)); /* FIXME */ + return 0.5 * z * h2 * b0; +} + + +int +gsl_sf_bessel_zero_J0_e(unsigned int s, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if(s == 0){ + result->val = 0.0; + result->err = 0.0; + GSL_ERROR ("error", GSL_EINVAL); + } + else { + /* See [F. Lether, J. Comp. Appl .Math. 67, 167 (1996)]. */ + + static const double P[] = { 1567450796.0/12539606369.0, + 8903660.0/2365861.0, + 10747040.0/536751.0, + 17590991.0/1696654.0 + }; + static const double Q[] = { 1.0, + 29354255.0/954518.0, + 76900001.0/431847.0, + 67237052.0/442411.0 + }; + + const double beta = (s - 0.25) * M_PI; + const double bi2 = 1.0/(beta*beta); + const double R33num = P[0] + bi2 * (P[1] + bi2 * (P[2] + P[3] * bi2)); + const double R33den = Q[0] + bi2 * (Q[1] + bi2 * (Q[2] + Q[3] * bi2)); + const double R33 = R33num/R33den; + result->val = beta + R33/beta; + result->err = fabs(3.0e-15 * result->val); + return GSL_SUCCESS; + } +} + + +int +gsl_sf_bessel_zero_J1_e(unsigned int s, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if(s == 0) { + result->val = 0.0; + result->err = 0.0; + return GSL_SUCCESS; + } + else { + /* See [M. Branders et al., J. Comp. Phys. 42, 403 (1981)]. */ + + static const double a[] = { -0.362804405737084, + 0.120341279038597, + 0.439454547101171e-01, + 0.159340088474713e-02 + }; + static const double b[] = { 1.0, + -0.325641790801361, + -0.117453445968927, + -0.424906902601794e-02 + }; + + const double beta = (s + 0.25) * M_PI; + const double bi2 = 1.0/(beta*beta); + const double Rnum = a[3] + bi2 * (a[2] + bi2 * (a[1] + bi2 * a[0])); + const double Rden = b[3] + bi2 * (b[2] + bi2 * (b[1] + bi2 * b[0])); + const double R = Rnum/Rden; + result->val = beta * (1.0 + R*bi2); + result->err = fabs(2.0e-14 * result->val); + return GSL_SUCCESS; + } +} + + +int +gsl_sf_bessel_zero_Jnu_e(double nu, unsigned int s, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if(nu <= -1.0) { + DOMAIN_ERROR(result); + } + else if(s == 0) { + result->val = 0.0; + result->err = 0.0; + if (nu == 0.0) { + GSL_ERROR ("no zero-th root for nu = 0.0", GSL_EINVAL); + } + return GSL_SUCCESS; + } + else if(nu < 0.0) { + /* This can be done, I'm just lazy now. */ + result->val = 0.0; + result->err = 0.0; + GSL_ERROR("unimplemented", GSL_EUNIMPL); + } + else if(s == 1) { + /* Chebyshev fits for the first positive zero. + * For some reason Nemeth made this different from the others. + */ + if(nu < 2.0) { + const double * c = coef_jnu_a[s]; + const size_t L = size_jnu_a[s]; + const double arg = nu/2.0; + const double chb = clenshaw(c, L-1, arg); + result->val = chb; + result->err = 2.0e-15 * result->val; + } + else { + const double * c = coef_jnu_b[s]; + const size_t L = size_jnu_b[s]; + const double arg = pow(2.0/nu, 2.0/3.0); + const double chb = clenshaw(c, L-1, arg); + result->val = nu * chb; + result->err = 2.0e-15 * result->val; + } + return GSL_SUCCESS; + } + else if(s <= 10) { + /* Chebyshev fits for the first 10 positive zeros. */ + if(nu < s) { + const double * c = coef_jnu_a[s]; + const size_t L = size_jnu_a[s]; + const double arg = nu/s; + const double chb = clenshaw(c, L-1, arg); + result->val = chb; + result->err = 2.0e-15 * result->val; + } + else { + const double * c = coef_jnu_b[s]; + const size_t L = size_jnu_b[s]; + const double arg = pow(s/nu, 2.0/3.0); + const double chb = clenshaw(c, L-1, arg); + result->val = nu * chb; + result->err = 2.0e-15 * result->val; + + /* FIXME: truth in advertising for the screwed up + * s = 5 fit. Need to fix that. + */ + if(s == 5) { + result->err *= 5.0e+06; + } + } + return GSL_SUCCESS; + } + else if(s > 0.5*nu && s <= 20) { + /* Chebyshev fits for 10 < s <= 20. */ + const double * c = coef_jnu_a[s]; + const size_t L = size_jnu_a[s]; + const double arg = nu/(2.0*s); + const double chb = clenshaw(c, L-1, arg); + result->val = chb; + result->err = 4.0e-15 * chb; + return GSL_SUCCESS; + } + else if(s > 2.0 * nu) { + /* McMahon expansion if s is large compared to nu. */ + const double beta = (s + 0.5*nu - 0.25) * M_PI; + const double mc = mcmahon_correction(4.0*nu*nu, beta); + gsl_sf_result rat12; + gsl_sf_pow_int_e(nu/beta, 14, &rat12); + result->val = beta * mc; + result->err = 4.0 * fabs(beta) * rat12.val; + result->err += 4.0 * fabs(GSL_DBL_EPSILON * result->val); + return GSL_SUCCESS; + } + else { + /* Olver uniform asymptotic. */ + gsl_sf_result as; + const int stat_as = gsl_sf_airy_zero_Ai_e(s, &as); + const double minus_zeta = -pow(nu,-2.0/3.0) * as.val; + const double z = gsl_sf_bessel_Olver_zofmzeta(minus_zeta); + const double f1 = olver_f1(z, minus_zeta); + result->val = nu * (z + f1/(nu*nu)); + result->err = 0.001/(nu*nu*nu); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return stat_as; + } +} + + +/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ + +#include "eval.h" + +double gsl_sf_bessel_zero_J0(unsigned int s) +{ + EVAL_RESULT(gsl_sf_bessel_zero_J0_e(s, &result)); +} + +double gsl_sf_bessel_zero_J1(unsigned int s) +{ + EVAL_RESULT(gsl_sf_bessel_zero_J1_e(s, &result)); +} + +double gsl_sf_bessel_zero_Jnu(double nu, unsigned int s) +{ + EVAL_RESULT(gsl_sf_bessel_zero_Jnu_e(nu, s, &result)); +} |