dbessel.c

Go to the documentation of this file.

/**
 * @file dbessel.c
 * @brief Implements Bessel functions (I0, I1, J0, J1, K0, K1, Y0, Y1) in double precision.
 *
 * This file provides double-precision implementations of the modified and regular
 * Bessel functions. The code includes precise algorithms optimized for different
 * ranges of input values.
 *
 * @copyright 
 *   - (c) 1996 Takuya OOURA (email: ooura@mmm.t.u-tokyo.ac.jp)
 *   - (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
 *   - (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.
 *
 * @license 
 * This file is distributed under the terms of the Software License Agreement
 * found in the file LICENSE included with this distribution.
 *
 * @author Takuya OOURA, M. Borland, H. Shang, R. Soliday
 */
 
/*
    Copyright(C) 1996 Takuya OOURA (email: ooura@mmm.t.u-tokyo.ac.jp).
    You may use, copy, modify this code for any purpose and 
    without fee. You may distribute this ORIGINAL package.
*/
 
 
#include <mdb.h>
 
/**
 * @brief Computes the modified Bessel function of the first kind I0(x) in double precision.
 * 
 * @param x The input value.
 * @return The computed value of I0(x).
 */
double dbesi0(double x) {
  int k;
  double w, t, y;
  static double a[65] = {
    8.5246820682016865877e-11, 2.5966600546497407288e-9,
    7.9689994568640180274e-8, 1.9906710409667748239e-6,
    4.0312469446528002532e-5, 6.4499871606224265421e-4,
    0.0079012345761930579108, 0.071111111109207045212,
    0.444444444444724909, 1.7777777777777532045,
    4.0000000000000011182, 3.99999999999999998,
    1.0000000000000000001,
    1.1520919130377195927e-10, 2.2287613013610985225e-9,
    8.1903951930694585113e-8, 1.9821560631611544984e-6,
    4.0335461940910133184e-5, 6.4495330974432203401e-4,
    0.0079013012611467520626, 0.071111038160875566622,
    0.44444450319062699316, 1.7777777439146450067,
    4.0000000132337935071, 3.9999999968569015366,
    1.0000000003426703174,
    1.5476870780515238488e-10, 1.2685004214732975355e-9,
    9.2776861851114223267e-8, 1.9063070109379044378e-6,
    4.0698004389917945832e-5, 6.4370447244298070713e-4,
    0.0079044749458444976958, 0.071105052411749363882,
    0.44445280640924755082, 1.7777694934432109713,
    4.0000055808824003386, 3.9999977081165740932,
    1.0000004333949319118,
    2.0675200625006793075e-10, -6.1689554705125681442e-10,
    1.2436765915401571654e-7, 1.5830429403520613423e-6,
    4.2947227560776583326e-5, 6.3249861665073441312e-4,
    0.0079454472840953930811, 0.070994327785661860575,
    0.44467219586283000332, 1.7774588182255374745,
    4.0003038986252717972, 3.9998233869142057195,
    1.0000472932961288324,
    2.7475684794982708655e-10, -3.8991472076521332023e-9,
    1.9730170483976049388e-7, 5.9651531561967674521e-7,
    5.1992971474748995357e-5, 5.7327338675433770752e-4,
    0.0082293143836530412024, 0.069990934858728039037,
    0.44726764292723985087, 1.7726685170014087784,
    4.0062907863712704432, 3.9952750700487845355,
    1.0016354346654179322};
  static double b[70] = {
    6.7852367144945531383e-8, 4.6266061382821826854e-7,
    6.9703135812354071774e-6, 7.6637663462953234134e-5,
    7.9113515222612691636e-4, 0.0073401204731103808981,
    0.060677114958668837046, 0.43994941411651569622,
    2.7420017097661750609, 14.289661921740860534,
    59.820609640320710779, 188.78998681199150629,
    399.8731367825601118, 427.56411572180478514,
    1.8042097874891098754e-7, 1.2277164312044637357e-6,
    1.8484393221474274861e-5, 2.0293995900091309208e-4,
    0.0020918539850246207459, 0.019375315654033949297,
    0.15985869016767185908, 1.1565260527420641724,
    7.1896341224206072113, 37.354773811947484532,
    155.80993164266268457, 489.5211371158540918,
    1030.9147225169564806, 1093.5883545113746958,
    4.8017305613187493564e-7, 3.261317843912380074e-6,
    4.9073137508166159639e-5, 5.3806506676487583755e-4,
    0.0055387918291051866561, 0.051223717488786549025,
    0.42190298621367914765, 3.0463625987357355872,
    18.895299447327733204, 97.915189029455461554,
    407.13940115493494659, 1274.3088990480582632,
    2670.9883037012547506, 2815.7166284662544712,
    1.2789926338424623394e-6, 8.6718263067604918916e-6,
    1.3041508821299929489e-4, 0.001428224737372747892,
    0.014684070635768789378, 0.13561403190404185755,
    1.1152592585977393953, 8.0387088559465389038,
    49.761318895895479206, 257.2684232313529138,
    1066.8543146269566231, 3328.3874581009636362,
    6948.8586598121634874, 7288.4893398212481055,
    3.409350368197032893e-6, 2.3079025203103376076e-5,
    3.4691373283901830239e-4, 0.003794994977222908545,
    0.038974209677945602145, 0.3594948380414878371,
    2.9522878893539528226, 21.246564609514287056,
    131.28727387146173141, 677.38107093296675421,
    2802.3724744545046518, 8718.5731420798254081,
    18141.348781638832286, 18948.925349296308859};
  static double c[45] = {
    2.5568678676452702768e-15, 3.0393953792305924324e-14,
    6.3343751991094840009e-13, 1.5041298011833009649e-11,
    4.4569436918556541414e-10, 1.746393051427167951e-8,
    1.0059224011079852317e-6, 1.0729838945088577089e-4,
    0.05150322693642527738,
    5.2527963991711562216e-15, 7.202118481421005641e-15,
    7.2561421229904797156e-13, 1.482312146673104251e-11,
    4.4602670450376245434e-10, 1.7463600061788679671e-8,
    1.005922609132234756e-6, 1.0729838937545111487e-4,
    0.051503226936437300716,
    1.3365917359358069908e-14, -1.2932643065888544835e-13,
    1.7450199447905602915e-12, 1.0419051209056979788e-11,
    4.58047881980598326e-10, 1.7442405450073548966e-8,
    1.0059461453281292278e-6, 1.0729837434500161228e-4,
    0.051503226940658446941,
    5.3771611477352308649e-14, -1.1396193006413731702e-12,
    1.2858641335221653409e-11, -5.9802086004570057703e-11,
    7.3666894305929510222e-10, 1.6731837150730356448e-8,
    1.0070831435812128922e-6, 1.0729733111203704813e-4,
    0.051503227360726294675,
    3.7819492084858931093e-14, -4.8600496888588034879e-13,
    1.6898350504817224909e-12, 4.5884624327524255865e-11,
    1.2521615963377513729e-10, 1.8959658437754727957e-8,
    1.0020716710561353622e-6, 1.073037119856927559e-4,
    0.05150322383300230775};
 
  w = fabs(x);
  if (w < 8.5) {
    t = w * w * 0.0625;
    k = 13 * ((int)t);
    y = (((((((((((a[k] * t + a[k + 1]) * t +
                  a[k + 2]) *
                   t +
                 a[k + 3]) *
                  t +
                a[k + 4]) *
                 t +
               a[k + 5]) *
                t +
              a[k + 6]) *
               t +
             a[k + 7]) *
              t +
            a[k + 8]) *
             t +
           a[k + 9]) *
            t +
          a[k + 10]) *
           t +
         a[k + 11]) *
          t +
        a[k + 12];
  } else if (w < 12.5) {
    k = (int)w;
    t = w - k;
    k = 14 * (k - 8);
    y = ((((((((((((b[k] * t + b[k + 1]) * t +
                   b[k + 2]) *
                    t +
                  b[k + 3]) *
                   t +
                 b[k + 4]) *
                  t +
                b[k + 5]) *
                 t +
               b[k + 6]) *
                t +
              b[k + 7]) *
               t +
             b[k + 8]) *
              t +
            b[k + 9]) *
             t +
           b[k + 10]) *
            t +
          b[k + 11]) *
           t +
         b[k + 12]) *
          t +
        b[k + 13];
  } else {
    t = 60 / w;
    k = 9 * ((int)t);
    y = ((((((((c[k] * t + c[k + 1]) * t +
               c[k + 2]) *
                t +
              c[k + 3]) *
               t +
             c[k + 4]) *
              t +
            c[k + 5]) *
             t +
           c[k + 6]) *
            t +
          c[k + 7]) *
           t +
         c[k + 8]) *
        sqrt(t) * exp(w);
  }
  return y;
}
 
/**
 * @brief Computes the modified Bessel function of the first kind I1(x) in double precision.
 * 
 * @param x The input value.
 * @return The computed value of I1(x).
 */
double dbesi1(double x) {
  int k;
  double w, t, y;
  static double a[60] = {
    1.2787464404046789181e-10, 3.5705860060088241077e-9,
    9.961153761934733504e-8, 2.2395070088633043177e-6,
    4.0312466928887462346e-5, 5.6437387840203722356e-4,
    0.0059259259312934746096, 0.04444444444349900887,
    0.22222222222232042719, 0.66666666666666139867,
    1.0000000000000001106, 0.49999999999999999962,
    1.7281952384448634449e-10, 3.064720455997639013e-9,
    1.0237662138842827028e-7, 2.2299494417341498163e-6,
    4.0335364374929326943e-5, 5.6433440269141349899e-4,
    0.0059259754885893798654, 0.04444439941088039787,
    0.2222222511283502673, 0.66666665422146063244,
    1.0000000032274936821, 0.49999999961866867205,
    2.3216048939948030996e-10, 1.7443372702334489579e-9,
    1.1596478963485415499e-7, 2.1446755518623035147e-6,
    4.0697440347437076195e-5, 5.6324394900433192204e-4,
    0.0059283484996093060678, 0.044440673899150997921,
    0.2222263801685265786, 0.66666358151576732094,
    1.0000013834029985337, 0.49999971643129650249,
    3.1013758938255172562e-10, -8.4813676145611694984e-10,
    1.5544980187411802596e-7, 1.7811109378708045726e-6,
    4.2945322199060856985e-5, 5.5344850176852353639e-4,
    0.0059590327716950614802, 0.044371611097707060659,
    0.22233578241986401111, 0.6665474730046331531,
    1.0000756505206705927, 0.49997803664415994554,
    4.1214758313965020365e-10, -5.361331773534742944e-9,
    2.4661360807517345161e-7, 6.7144593918926723203e-7,
    5.1988027944945587571e-5, 5.0165568586065803067e-4,
    0.0061717530047005289953, 0.043745229577317251404,
    0.22363147971477747996, 0.6647546913111766024,
    1.0015686689447547657, 0.49941120439785391891};
  static double b[70] = {
    6.6324787943143095845e-8, 4.5125928898466638619e-7,
    6.7937793134877246623e-6, 7.4580507871505926302e-5,
    7.6866382927334005919e-4, 0.0071185174803491859307,
    0.058721838073486424416, 0.42473949281714196041,
    2.6396965606282079123, 13.710008536637016903,
    57.158647688180932003, 179.46182892089389037,
    377.57997362398478619, 399.87313678256009819,
    1.7652713206027939711e-7, 1.1988179244834708057e-6,
    1.8037851545747139231e-5, 1.9775785516370314656e-4,
    0.0020354870702829387283, 0.0188221641910322536,
    0.15500485219010424263, 1.119010001056057321,
    6.9391565185406617552, 35.948170579648649345,
    149.41909525103032616, 467.42979492780642582,
    979.04227423171290408, 1030.9147225169564443,
    4.7022299276154507603e-7, 3.1878571710170115972e-6,
    4.7940153875711448496e-5, 5.2496623508411440227e-4,
    0.0053968661134780824779, 0.049837081920693776234,
    0.40979593830387765545, 2.9533186922862948404,
    18.278176130722516369, 94.47649715018912107,
    391.66075612645333624, 1221.4182034643210345,
    2548.6177980961291004, 2670.9883037012546541,
    1.2535083724002034147e-6, 8.484587142065570825e-6,
    1.2753227372734042108e-4, 0.0013950105363562648921,
    0.014325473993765291906, 0.13212452778932829125,
    1.0849287786885151432, 7.8068089156260172673,
    48.232254570679165833, 248.80659424902394371,
    1029.0736929484210803, 3200.5629438795801652,
    6656.7749162019607914, 6948.8586598121632302,
    3.3439394490599745013e-6, 2.2600596902211837757e-5,
    3.3955927589987356838e-4, 0.0037105306061050972474,
    0.038065263634919156421, 0.35068223415665236079,
    2.8760027832105027316, 20.665999500843274339,
    127.47939148516390205, 656.43636874254000885,
    2709.524283793247992, 8407.1174233600734871,
    17437.146284159740233, 18141.3487816388316};
  static double c[45] = {
    -2.8849790431465382128e-15, -3.5125350943844774657e-14,
    -7.485086701370741975e-13, -1.8383904048277485153e-11,
    -5.7303556446977223342e-10, -2.4449502737311496525e-8,
    -1.6765373351766929724e-6, -3.2189516835265773471e-4,
    0.051503226936425277377,
    -5.8674306822281631119e-15, -9.4884898451194085565e-15,
    -8.503386513660036434e-13, -1.8142997866945285736e-11,
    -5.7340238386338193949e-10, -2.4449138101742183665e-8,
    -1.6765375646678855842e-6, -3.2189516826945356325e-4,
    0.051503226936412017608,
    -1.4723362506764340882e-14, 1.3945147385179042899e-13,
    -1.9618041857586930923e-12, -1.3343606394065121821e-11,
    -5.8649674606973244159e-10, -2.4426060539669553778e-8,
    -1.6765631828366988006e-6, -3.2189515191449587253e-4,
    0.051503226931820146445,
    -5.8203519372580372987e-14, 1.2266326995309845825e-12,
    -1.3921625844526453237e-11, 6.2228025878281625469e-11,
    -8.8636681342142794023e-10, -2.3661241616744818608e-8,
    -1.6777870960740520557e-6, -3.2189402882677074318e-4,
    0.051503226479551959376,
    -4.5801527369223291722e-14, 6.7998819697143727209e-13,
    -4.1624857909290468421e-12, -3.2849009406112440998e-11,
    -3.247827569043111827e-10, -2.5739209934053714983e-8,
    -1.6730566573215739195e-6, -3.2190010909008684076e-4,
    0.05150322986693207715};
 
  w = fabs(x);
  if (w < 8.5) {
    t = w * w * 0.0625;
    k = 12 * ((int)t);
    y = (((((((((((a[k] * t + a[k + 1]) * t +
                  a[k + 2]) *
                   t +
                 a[k + 3]) *
                  t +
                a[k + 4]) *
                 t +
               a[k + 5]) *
                t +
              a[k + 6]) *
               t +
             a[k + 7]) *
              t +
            a[k + 8]) *
             t +
           a[k + 9]) *
            t +
          a[k + 10]) *
           t +
         a[k + 11]) *
        w;
  } else if (w < 12.5) {
    k = (int)w;
    t = w - k;
    k = 14 * (k - 8);
    y = ((((((((((((b[k] * t + b[k + 1]) * t +
                   b[k + 2]) *
                    t +
                  b[k + 3]) *
                   t +
                 b[k + 4]) *
                  t +
                b[k + 5]) *
                 t +
               b[k + 6]) *
                t +
              b[k + 7]) *
               t +
             b[k + 8]) *
              t +
            b[k + 9]) *
             t +
           b[k + 10]) *
            t +
          b[k + 11]) *
           t +
         b[k + 12]) *
          t +
        b[k + 13];
  } else {
    t = 60 / w;
    k = 9 * ((int)t);
    y = ((((((((c[k] * t + c[k + 1]) * t +
               c[k + 2]) *
                t +
              c[k + 3]) *
               t +
             c[k + 4]) *
              t +
            c[k + 5]) *
             t +
           c[k + 6]) *
            t +
          c[k + 7]) *
           t +
         c[k + 8]) *
        sqrt(t) * exp(w);
  }
  return x < 0 ? -y : y;
}
 
/**
 * @brief Computes the regular Bessel function of the first kind J0(x) in double precision.
 * 
 * @param x The input value.
 * @return The computed value of J0(x).
 */
double dbesj0(double x) {
  int k;
  double w, t, y, v, theta;
  static double a[8] = {
    -2.3655394e-12, 4.70889868e-10,
    -6.78167892231e-8, 6.7816840038636e-6,
    -4.340277777716935e-4, 0.0156249999999992397,
    -0.2499999999999999638, 0.9999999999999999997};
  static double b[65] = {
    6.26681117e-11, -2.2270614428e-9,
    6.62981656302e-8, -1.6268486502196e-6,
    3.21978384111685e-5, -5.00523773331583e-4,
    0.0059060313537449816, -0.0505265323740109701,
    0.2936432097610503985, -1.0482565081091638637,
    1.9181123286040428113, -1.13191994752217001,
    -0.1965480952704682,
    4.57457332e-11, -1.5814772025e-9,
    4.55487446311e-8, -1.0735201286233e-6,
    2.02015179970014e-5, -2.942392368203808e-4,
    0.0031801987726150648, -0.0239875209742846362,
    0.1141447698973777641, -0.2766726722823530233,
    0.1088620480970941648, 0.5136514645381999197,
    -0.2100594022073706033,
    3.31366618e-11, -1.1119090229e-9,
    3.08823040363e-8, -6.956602653104e-7,
    1.23499947481762e-5, -1.66295194539618e-4,
    0.0016048663165678412, -0.0100785479932760966,
    0.0328996815223415274, -0.0056168761733860688,
    -0.2341096400274429386, 0.2551729256776404262,
    0.2288438186148935667,
    2.38007203e-11, -7.731046439e-10,
    2.06237001152e-8, -4.412291442285e-7,
    7.3107766249655e-6, -8.91749801028666e-5,
    7.34165451384135e-4, -0.0033303085445352071,
    0.0015425853045205717, 0.0521100583113136379,
    -0.1334447768979217815, -0.1401330292364750968,
    0.2685616168804818919,
    1.6935595e-11, -5.308092192e-10,
    1.35323005576e-8, -2.726650587978e-7,
    4.151324014176e-6, -4.43353052220157e-5,
    2.815740758993879e-4, -4.393235121629007e-4,
    -0.0067573531105799347, 0.0369141914660130814,
    0.0081673361942996237, -0.257338128589888186,
    0.0459580257102978932};
  static double c[70] = {
    -3.009451757e-11, -1.4958003844e-10,
    5.06854544776e-9, 1.863564222012e-8,
    -6.0304249068078e-7, -1.47686259937403e-6,
    4.714331342682714e-5, 6.286305481740818e-5,
    -0.00214137170594124344, -8.9157336676889788e-4,
    0.04508258728666024989, -0.00490362805828762224,
    -0.27312196367405374426, 0.04193925184293450356,
    -7.1245356e-12, -4.1170814825e-10,
    1.38012624364e-9, 5.704447670683e-8,
    -1.9026363528842e-7, -5.33925032409729e-6,
    1.736064885538091e-5, 3.0692619152608375e-4,
    -9.2598938200644367e-4, -0.00917934265960017663,
    0.02287952522866389076, 0.10545197546252853195,
    -0.16126443075752985095, -0.19392874768742235538,
    2.128344556e-11, -3.1053910272e-10,
    -3.34979293158e-9, 4.50723289505e-8,
    3.6437959146427e-7, -4.46421436266678e-6,
    -2.523429344576552e-5, 2.7519882931758163e-4,
    9.7185076358599358e-4, -0.00898326746345390692,
    -0.01665959196063987584, 0.11456933464891967814,
    0.07885001422733148815, -0.23664819446234712621,
    3.035295055e-11, 5.486066835e-11,
    -5.01026824811e-9, -5.0124684786e-9,
    5.8012340163034e-7, 1.6788922416169e-7,
    -4.373270270147275e-5, 1.183898532719802e-5,
    0.00189863342862291449, -0.0011375924956163613,
    -0.03846797195329871681, 0.02389746880951420335,
    0.22837862066532347461, -0.06765394811166522844,
    1.279875977e-11, 3.5925958103e-10,
    -2.28037105967e-9, -4.852770517176e-8,
    2.8696428000189e-7, 4.40131125178642e-6,
    -2.366617753349105e-5, -2.4412456252884129e-4,
    0.00113028178539430542, 0.0070847051391978908,
    -0.02526914792327618386, -0.08006137953480093426,
    0.16548380461475971846, 0.14688405470042110229};
  static double d[52] = {
    1.059601355592185731e-14, -2.71150591218550377e-13,
    8.6514809056201638e-12, -4.6264028554286627e-10,
    5.0815403835647104e-8, -1.76722552048141208e-5,
    0.16286750396763997378, 2.949651820598278873e-13,
    -8.818215611676125741e-12, 3.571119876162253451e-10,
    -2.63192412099371706e-8, 4.709502795656698909e-6,
    -0.005208333333333283282,
    7.18344107717531977e-15, -2.51623725588410308e-13,
    8.6017784918920604e-12, -4.6256876614290359e-10,
    5.0815343220437937e-8, -1.7672255176494197e-5,
    0.16286750396763433767, 2.2327570859680094777e-13,
    -8.464594853517051292e-12, 3.563766464349055183e-10,
    -2.631843986737892965e-8, 4.70950234228865941e-6,
    -0.0052083333332278466225,
    5.15413392842889366e-15, -2.27740238380640162e-13,
    8.4827767197609014e-12, -4.6224753682737618e-10,
    5.0814848128929134e-8, -1.7672254763876748e-5,
    0.16286750396748926663, 1.7316195320192170887e-13,
    -7.971122772293919646e-12, 3.544039469911895749e-10,
    -2.631443902081701081e-8, 4.709498228695400603e-6,
    -0.005208333331514365361,
    3.84653681453798517e-15, -2.04464520778789011e-13,
    8.3089298605177838e-12, -4.6155016158412096e-10,
    5.081326369646665e-8, -1.76722528311426167e-5,
    0.1628675039665006593, 1.3797879972460878797e-13,
    -7.448089381011684812e-12, 3.51273379710695978e-10,
    -2.630500895563592722e-8, 4.709483934775839193e-6,
    -0.0052083333227940760113};
 
  w = fabs(x);
  if (w < 1) {
    t = w * w;
    y = ((((((a[0] * t + a[1]) * t +
             a[2]) *
              t +
            a[3]) *
             t +
           a[4]) *
            t +
          a[5]) *
           t +
         a[6]) *
          t +
        a[7];
  } else if (w < 8.5) {
    t = w * w * 0.0625;
    k = (int)t;
    t -= k + 0.5;
    k *= 13;
    y = (((((((((((b[k] * t + b[k + 1]) * t +
                  b[k + 2]) *
                   t +
                 b[k + 3]) *
                  t +
                b[k + 4]) *
                 t +
               b[k + 5]) *
                t +
              b[k + 6]) *
               t +
             b[k + 7]) *
              t +
            b[k + 8]) *
             t +
           b[k + 9]) *
            t +
          b[k + 10]) *
           t +
         b[k + 11]) *
          t +
        b[k + 12];
  } else if (w < 12.5) {
    k = (int)w;
    t = w - (k + 0.5);
    k = 14 * (k - 8);
    y = ((((((((((((c[k] * t + c[k + 1]) * t +
                   c[k + 2]) *
                    t +
                  c[k + 3]) *
                   t +
                 c[k + 4]) *
                  t +
                c[k + 5]) *
                 t +
               c[k + 6]) *
                t +
              c[k + 7]) *
               t +
             c[k + 8]) *
              t +
            c[k + 9]) *
             t +
           c[k + 10]) *
            t +
          c[k + 11]) *
           t +
         c[k + 12]) *
          t +
        c[k + 13];
  } else {
    v = 24 / w;
    t = v * v;
    k = 13 * ((int)t);
    y = ((((((d[k] * t + d[k + 1]) * t +
             d[k + 2]) *
              t +
            d[k + 3]) *
             t +
           d[k + 4]) *
            t +
          d[k + 5]) *
           t +
         d[k + 6]) *
        sqrt(v);
    theta = (((((d[k + 7] * t + d[k + 8]) * t +
                d[k + 9]) *
                 t +
               d[k + 10]) *
                t +
              d[k + 11]) *
               t +
             d[k + 12]) *
              v -
            0.78539816339744830962;
    y *= cos(w + theta);
  }
  return y;
}
 
/**
 * @brief Computes the regular Bessel function of the first kind J1(x) in double precision.
 * 
 * @param x The input value.
 * @return The computed value of J1(x).
 */
double dbesj1(double x) {
  int k;
  double w, t, y, v, theta;
  static double a[8] = {
    -1.4810349e-13, 3.363594618e-11,
    -5.65140051697e-9, 6.7816840144764e-7,
    -5.425347222188379e-5, 0.00260416666666662438,
    -0.06249999999999999799, 0.49999999999999999998};
  static double b[65] = {
    2.43721316e-12, -9.400554763e-11,
    3.0605338998e-9, -8.287270492518e-8,
    1.83020515991344e-6, -3.219783841164382e-5,
    4.3795830161515318e-4, -0.00442952351530868999,
    0.03157908273375945955, -0.14682160488052520107,
    0.39309619054093640008, -0.4795280821510107028,
    0.1414899934402712514,
    1.82119257e-12, -6.862117678e-11,
    2.1732790836e-9, -5.69359291782e-8,
    1.20771046483277e-6, -2.020151799736374e-5,
    2.5745933218048448e-4, -0.00238514907946126334,
    0.01499220060892984289, -0.05707238494868888345,
    0.10375225210588234727, -0.02721551202427354117,
    -0.06420643306727498985,
    1.352611196e-12, -4.9706947875e-11,
    1.527944986332e-9, -3.8602878823401e-8,
    7.82618036237845e-7, -1.23499947484511e-5,
    1.45508295194426686e-4, -0.001203649737425854162,
    0.006299092495799005109, -0.016449840761170764763,
    0.002106328565019748701, 0.05852741000686073465,
    -0.031896615709705053191,
    9.97982124e-13, -3.5702556073e-11,
    1.062332772617e-9, -2.5779624221725e-8,
    4.96382962683556e-7, -7.310776625173004e-6,
    7.8028107569541842e-5, -5.50624088538081113e-4,
    0.002081442840335570371, -7.71292652260286633e-4,
    -0.019541271866742634199, 0.033361194224480445382,
    0.017516628654559387164,
    7.31050661e-13, -2.5404499912e-11,
    7.29360079088e-10, -1.6915375004937e-8,
    3.06748319652546e-7, -4.151324014331739e-6,
    3.8793392054271497e-5, -2.11180556924525773e-4,
    2.74577195102593786e-4, 0.003378676555289966782,
    -0.013842821799754920148, -0.002041834048574905921,
    0.032167266073736023299};
  static double c[70] = {
    -1.185964494e-11, 3.9110295657e-10,
    1.80385519493e-9, -5.575391345723e-8,
    -1.8635897017174e-7, 5.42738239401869e-6,
    1.181490114244279e-5, -3.300031939852107e-4,
    -3.7717832892725053e-4, 0.01070685852970608288,
    0.00356629346707622489, -0.13524776185998074716,
    0.00980725611657523952, 0.27312196367405374425,
    -3.029591097e-11, 9.259293559e-11,
    4.96321971223e-9, -1.518137078639e-8,
    -5.7045127595547e-7, 1.71237271302072e-6,
    4.271400348035384e-5, -1.2152454198713258e-4,
    -0.00184155714921474963, 0.00462994691003219055,
    0.03671737063840232452, -0.06863857568599167175,
    -0.21090395092505707655, 0.16126443075752985095,
    -2.19760208e-11, -2.7659100729e-10,
    3.74295124827e-9, 3.684765777023e-8,
    -4.5072801091574e-7, -3.27941630669276e-6,
    3.5713715545163e-5, 1.7664005411843533e-4,
    -0.00165119297594774104, -0.00485925381792986774,
    0.03593306985381680131, 0.04997877588191962563,
    -0.22913866929783936544, -0.07885001422733148814,
    5.16292316e-12, -3.9445956763e-10,
    -6.6220021263e-10, 5.511286218639e-8,
    5.01257940078e-8, -5.22111059203425e-6,
    -1.34311394455105e-6, 3.0612891890766805e-4,
    -7.103391195326182e-5, -0.00949316714311443491,
    0.00455036998246516948, 0.11540391585989614784,
    -0.04779493761902840455, -0.2283786206653234746,
    2.697817493e-11, -1.6633326949e-10,
    -4.3313486035e-9, 2.508404686362e-8,
    4.8528284780984e-7, -2.58267851112118e-6,
    -3.521049080466759e-5, 1.6566324273339952e-4,
    0.00146474737522491617, -0.00565140892697147306,
    -0.028338820556793004, 0.07580744376982855057,
    0.16012275906960187978, -0.16548380461475971845};
  static double d[52] = {
    -1.272346002224188092e-14, 3.370464692346669075e-13,
    -1.144940314335484869e-11, 6.863141561083429745e-10,
    -9.491933932960924159e-8, 5.301676561445687562e-5,
    0.162867503967639974, -3.652982212914147794e-13,
    1.151126750560028914e-11, -5.165585095674343486e-10,
    4.657991250060549892e-8, -1.186794704692706504e-5,
    0.01562499999999994026,
    -8.713069680903981555e-15, 3.140780373478474935e-13,
    -1.139089186076256597e-11, 6.862299023338785566e-10,
    -9.491926788274594674e-8, 5.301676558106268323e-5,
    0.162867503967646622, -2.792555727162752006e-13,
    1.108650207651756807e-11, -5.156745588549830981e-10,
    4.657894859077370979e-8, -1.186794650130550256e-5,
    0.01562499999987299901,
    -6.304859171204770696e-15, 2.857249044208791652e-13,
    -1.124956921556753188e-11, 6.858482894906716661e-10,
    -9.49186795351689846e-8, 5.301676509057781574e-5,
    0.1628675039678191167, -2.185193490132496053e-13,
    1.048820673697426074e-11, -5.132819367467680132e-10,
    4.65740943737299422e-8, -1.186794150862988921e-5,
    0.01562499999779270706,
    -4.74041720979200985e-15, 2.578715253644144182e-13,
    -1.104148898414138857e-11, 6.850134201626289183e-10,
    -9.49167823417491964e-8, 5.301676277588728159e-5,
    0.1628675039690033136, -1.75512205749384229e-13,
    9.848723331445182397e-12, -5.094535425482245697e-10,
    4.656255982268609304e-8, -1.186792402114394891e-5,
    0.01562499998712198636};
 
  w = fabs(x);
  if (w < 1) {
    t = w * w;
    y = (((((((a[0] * t + a[1]) * t +
              a[2]) *
               t +
             a[3]) *
              t +
            a[4]) *
             t +
           a[5]) *
            t +
          a[6]) *
           t +
         a[7]) *
        w;
  } else if (w < 8.5) {
    t = w * w * 0.0625;
    k = (int)t;
    t -= k + 0.5;
    k *= 13;
    y = ((((((((((((b[k] * t + b[k + 1]) * t +
                   b[k + 2]) *
                    t +
                  b[k + 3]) *
                   t +
                 b[k + 4]) *
                  t +
                b[k + 5]) *
                 t +
               b[k + 6]) *
                t +
              b[k + 7]) *
               t +
             b[k + 8]) *
              t +
            b[k + 9]) *
             t +
           b[k + 10]) *
            t +
          b[k + 11]) *
           t +
         b[k + 12]) *
        w;
  } else if (w < 12.5) {
    k = (int)w;
    t = w - (k + 0.5);
    k = 14 * (k - 8);
    y = ((((((((((((c[k] * t + c[k + 1]) * t +
                   c[k + 2]) *
                    t +
                  c[k + 3]) *
                   t +
                 c[k + 4]) *
                  t +
                c[k + 5]) *
                 t +
               c[k + 6]) *
                t +
              c[k + 7]) *
               t +
             c[k + 8]) *
              t +
            c[k + 9]) *
             t +
           c[k + 10]) *
            t +
          c[k + 11]) *
           t +
         c[k + 12]) *
          t +
        c[k + 13];
  } else {
    v = 24 / w;
    t = v * v;
    k = 13 * ((int)t);
    y = ((((((d[k] * t + d[k + 1]) * t +
             d[k + 2]) *
              t +
            d[k + 3]) *
             t +
           d[k + 4]) *
            t +
          d[k + 5]) *
           t +
         d[k + 6]) *
        sqrt(v);
    theta = (((((d[k + 7] * t + d[k + 8]) * t +
                d[k + 9]) *
                 t +
               d[k + 10]) *
                t +
              d[k + 11]) *
               t +
             d[k + 12]) *
              v -
            0.78539816339744830962;
    y *= sin(w + theta);
  }
  return x < 0 ? -y : y;
}
 
/**
 * @brief Computes the modified Bessel function of the second kind K0(x) in double precision.
 * 
 * @param x The input value.
 * @return The computed value of K0(x).
 */
double dbesk0(double x) {
  int k;
  double t, y;
  static double a[16] = {
    2.4307270476772195953e-12, 4.709166636378530437e-10,
    6.7816861334344265568e-8, 6.7816840204737508252e-6,
    4.3402777777915334676e-4, 0.015624999999999872796,
    0.25000000000000000448, 0.99999999999999999997,
    6.5878327432224993071e-12, 1.2083308769932888218e-9,
    1.6271062073716412046e-7, 1.4914719278555277887e-5,
    8.4603509071212245667e-4, 0.02524892993216233391,
    0.27898287891460312491, 0.11593151565841244874};
  static double b[112] = {
    -4.6430702971053162197e-13, 1.037793605956372823e-11,
    -1.0298475936392057807e-10, 5.3632747492333959219e-10,
    -2.1674628861036068105e-10, -2.3316071545820437669e-8,
    2.2557819578691704059e-7, -9.2325694638587080009e-7,
    -3.3569097781613661759e-6, 8.7355061305812582974e-5,
    -6.8021202111645760475e-4, 2.7434654781323362319e-4,
    0.10031787169953909561, 0.42102443824070833334,
    4.1447451117883103686e-12, -3.4026589638576604315e-11,
    9.3398790624638977468e-12, 1.518418175079985263e-9,
    -1.1364911665083029464e-8, 2.0619457602095915719e-8,
    3.043101803757224363e-7, -2.974973626447455551e-6,
    8.0143661611467038568e-6, 8.0937525149549218398e-5,
    -0.0010356346549612699886, 0.0028534806627578638795,
    0.097369634474060441807, 0.32175066577856452683,
    1.117088257074072752e-13, -8.2865909408297066068e-11,
    9.4656678749191182763e-10, -3.583201984184788338e-9,
    -9.5017955656904252761e-9, 1.5200595674883329093e-7,
    -3.866326257135605998e-7, -3.3350340828235103499e-6,
    2.9359886663960844231e-5, -1.1266401822556801563e-5,
    -0.0012113572742435576205, 0.0063158973673701376253,
    0.088291790250128171341, 0.22833982383240512262,
    -3.2880638807053948433e-11, 4.3194884830465283512e-10,
    -1.7455089683104033093e-9, -3.2437330799994764516e-9,
    4.7393655539139519778e-8, -1.1929265603456272466e-8,
    -1.3177845881013419388e-6, 3.3873375636197969526e-6,
    3.2729835880668256625e-5, -1.8367283883002494561e-4,
    -8.2830996454188084408e-4, 0.0095512732229514251931,
    0.072233832113719266702, 0.14753187103603405298,
    7.9998492614150860098e-11, -7.025734670268613949e-10,
    7.889882162708458627e-10, 1.1294796399671507085e-8,
    -1.1360539648638059137e-8, -3.0346309115270564487e-7,
    3.2235585426189451721e-7, 8.3575612102298214948e-6,
    -8.5169628089198208211e-6, -2.5740175232173357342e-4,
    1.246273401468915277e-4, 0.01068323286919220345,
    0.051515690033637395779, 0.085465862953544883657,
    -8.6111506537356531608e-11, 5.1862926131024597823e-10,
    7.5884324949371110022e-10, -6.4011975813006767417e-9,
    -4.1966181325111763156e-8, 9.1306285446881485314e-8,
    1.3573638315827954034e-6, 4.8683213252735694701e-7,
    -3.8805424608710197066e-5, -1.183898646868898061e-4,
    9.2796213947750964945e-4, 0.0089611057737319027776,
    0.031464453915862785606, 0.04426764808753663078,
    4.4400123834164610288e-11, -1.1411233140911074336e-10,
    -8.820067070246705983e-10, -1.9686735373323381456e-9,
    1.9921120728941773855e-8, 1.454397441858483474e-7,
    1.8238418041265854754e-8, -4.5363700392899066037e-6,
    -2.1688068222527688542e-5, 4.54960621666870347e-5,
    0.0010435238076080528284, 0.0058374528996419979931,
    0.01661121071042545585, 0.020756008367065750538,
    -6.5166519951106397214e-12, -5.857218285878853958e-11,
    1.5550375065815375404e-10, 1.9526509484993563229e-9,
    9.2637123346818426594e-9, -1.4136471501812055943e-8,
    -4.3024895710889717172e-7, -2.3235612243330592076e-6,
    4.0380616133862188804e-7, 9.2783767992909743602e-5,
    7.2964887597817095035e-4, 0.0031316245282223273413,
    0.0078028233022066428316, 0.0090014807263791058095};
  static double c[135] = {
    4.5161032649342790231e-11, -4.2774336988557091369e-11,
    6.0998467173896677777e-10, 1.9845167242599996944e-9,
    1.3097678767280215271e-8, 7.4505822268382641286e-8,
    4.2893920879106814989e-7, 2.3900851955655303104e-6,
    1.2533473009382380357e-5, 5.9693359063879871983e-5,
    2.477507066108730458e-4, 8.5106703131389516508e-4,
    0.0022500105115665788755, 0.00404461344545216346,
    0.0036910983340425942762,
    3.5732826433251464989e-12, -3.2906649482312266258e-12,
    7.0873811190464760555e-11, 2.955132058048417712e-10,
    2.2776940472505079894e-9, 1.5175463612815010036e-8,
    9.9462487812170164133e-8, 6.14487577978539011e-7,
    3.486953188290736075e-6, 1.7615836644757657443e-5,
    7.6373536037879531886e-5, 2.7098571871205999668e-4,
    7.3399047381788927036e-4, 0.0013439197177355085297,
    0.0012439943280131230863,
    3.6343547836242523646e-13, 9.7997961751276137602e-14,
    1.0184692699811569047e-11, 6.1495184828957652064e-11,
    5.0238328349302602543e-10, 3.7498626376004337661e-9,
    2.6689445483857236307e-8, 1.7591899737346368084e-7,
    1.0486448307010701679e-6, 5.4986458466257148573e-6,
    2.4521456351751345323e-5, 8.8900942259143832228e-5,
    2.448394771406830019e-4, 4.5418248688489693045e-4,
    4.2479574186923180694e-4,
    5.2460389348163395857e-14, 7.480206302650350354e-14,
    2.0012201610651998417e-12, 1.4887306044735163359e-11,
    1.294670541423294035e-10, 1.0391628915892803144e-9,
    7.8091180499677328456e-9, 5.3694223626907660084e-8,
    3.3063914804658509029e-7, 1.7776972424421486506e-6,
    8.0833148098458320202e-6, 2.975555630444881778e-5,
    8.2945928349220642178e-5, 1.5536921180500112883e-4,
    1.4647070522281538711e-4,
    9.7531436733955514559e-15, 2.4084291220447154982e-14,
    4.7654956400897494468e-13, 4.0200949504810597783e-12,
    3.6726577109162191533e-11, 3.0939005665422637601e-10,
    2.4122848979784500179e-9, 1.7071884462645525505e-8,
    1.0752238955654933405e-7, 5.8844190041189462347e-7,
    2.7136083303224014597e-6, 1.0102477728604441135e-5,
    2.8420490721532571809e-5, 5.3637016379451944413e-5,
    5.0881312956459247572e-5,
    2.173204986818937726e-15, 7.2720052142815590531e-15,
    1.28030837955368201e-13, 1.1696825543787717167e-12,
    1.1083298191597132094e-11, 9.6536661252658773139e-11,
    7.7242553835198536397e-10, 5.579836626711057562e-9,
    3.572134529654341437e-8, 1.9806931547193682466e-7,
    9.2312964655319555313e-7, 3.4666258590861079959e-6,
    9.8224698307751177077e-6, 1.8648773453825584428e-5,
    1.7780062316167651812e-5,
    5.5012463763851934112e-16, 2.2254763392767319419e-15,
    3.7187669817701214965e-14, 3.5819585377733489628e-13,
    3.4866061263191556694e-12, 3.110163345062965291e-11,
    2.5358235662235617663e-10, 1.8597629779492599046e-9,
    1.2052654739462999992e-8, 6.7501417351172136833e-8,
    3.1720052198654584574e-7, 1.1993651363602981832e-6,
    3.4179130317623363474e-6, 6.5208606745808860158e-6,
    6.2430205476536771454e-6,
    1.5225407517829491689e-16, 6.9834820025664405161e-16,
    1.1380182837138781431e-14, 1.1369488761077196511e-13,
    1.1291168681618466716e-12, 1.0250757630526871007e-11,
    8.4765287317253141514e-11, 6.2886627779402596211e-10,
    4.1142865598366029316e-9, 2.3223773435632014408e-8,
    1.0985095234166396934e-7, 4.1766260951820336228e-7,
    1.1958609263543792991e-6, 2.2907574647671878055e-6,
    2.2008253973114914005e-6,
    4.4863058691420695911e-17, 2.2437356594371819978e-16,
    3.6107964803015652759e-15, 3.7031193629853392081e-14,
    3.7341552790439784371e-13, 3.4355950129497564468e-12,
    2.8719942600171304499e-11, 2.1499646844509516453e-10,
    1.4171810843455227171e-9, 8.0501118772875784153e-9,
    3.8281889106330295876e-8, 1.4621673458431979989e-7,
    4.2029868696411098586e-7, 8.0785884122023473025e-7,
    7.7845438614204963209e-7};
  static double d[40] = {
    -7.9737703860537066166e-14, 1.9543834380466766627e-12,
    -4.7230794431646733538e-11, 1.4001773785771252004e-9,
    -5.4864553020583098585e-8, 3.1601984250143742772e-6,
    -3.3708783204090252161e-4, 0.16180215937964160437,
    -5.2593898374798632343e-14, 1.7725913926973236457e-12,
    -4.6672234858122387294e-11, 1.3991653503828889207e-9,
    -5.4863400156413929639e-8, 3.1601976099900075541e-6,
    -3.3708783171335864627e-4, 0.1618021593795843376,
    -3.6135496189875398132e-14, 1.5466239429618130284e-12,
    -4.5320259146602122624e-11, 1.3945974109459385552e-9,
    -5.4853994841172088787e-8, 3.1601858228022739196e-6,
    -3.370878233999830232e-4, 0.16180215937704286491,
    -2.5640663123518180635e-14, 1.3288079339404032671e-12,
    -4.3368537955908371563e-11, 1.3848103653102203186e-9,
    -5.4824335664256344123e-8, 3.1601315173126153586e-6,
    -3.370877677903569564e-4, 0.16180215935248373474,
    -1.8678321325292127767e-14, 1.1354310934105733311e-12,
    -4.1057197297998608931e-11, 1.369399096129635097e-9,
    -5.4762428935047089835e-8, 3.1599817092775027963e-6,
    -3.3708756559715893599e-4, 0.1618021592350814424};
 
  if (x < 0.86) {
    t = x * x;
    y = ((((((a[0] * t + a[1]) * t +
             a[2]) *
              t +
            a[3]) *
             t +
           a[4]) *
            t +
          a[5]) *
           t +
         a[6]) *
          t +
        a[7];
    y = ((((((a[8] * t + a[9]) * t +
             a[10]) *
              t +
            a[11]) *
             t +
           a[12]) *
            t +
          a[13]) *
           t +
         a[14]) *
          t +
        a[15] - y * log(x);
  } else if (x < 4.15) {
    t = x - 5 / x;
    k = (int)(t + 5);
    t = (k - 4) - t;
    k *= 14;
    y = ((((((((((((b[k] * t + b[k + 1]) * t +
                   b[k + 2]) *
                    t +
                  b[k + 3]) *
                   t +
                 b[k + 4]) *
                  t +
                b[k + 5]) *
                 t +
               b[k + 6]) *
                t +
              b[k + 7]) *
               t +
             b[k + 8]) *
              t +
            b[k + 9]) *
             t +
           b[k + 10]) *
            t +
          b[k + 11]) *
           t +
         b[k + 12]) *
          t +
        b[k + 13];
  } else if (x < 12.5) {
    k = (int)x;
    t = (k + 1) - x;
    k = 15 * (k - 4);
    y = (((((((((((((c[k] * t + c[k + 1]) * t +
                    c[k + 2]) *
                     t +
                   c[k + 3]) *
                    t +
                  c[k + 4]) *
                   t +
                 c[k + 5]) *
                  t +
                c[k + 6]) *
                 t +
               c[k + 7]) *
                t +
              c[k + 8]) *
               t +
             c[k + 9]) *
              t +
            c[k + 10]) *
             t +
           c[k + 11]) *
            t +
          c[k + 12]) *
           t +
         c[k + 13]) *
          t +
        c[k + 14];
  } else {
    t = 60 / x;
    k = 8 * ((int)t);
    y = (((((((d[k] * t + d[k + 1]) * t +
              d[k + 2]) *
               t +
             d[k + 3]) *
              t +
            d[k + 4]) *
             t +
           d[k + 5]) *
            t +
          d[k + 6]) *
           t +
         d[k + 7]) *
        sqrt(t) * exp(-x);
  }
  return y;
}
 
/**
 * @brief Computes the modified Bessel function of the second kind K1(x) in double precision.
 * 
 * @param x The input value.
 * @return The computed value of K1(x).
 */
double dbesk1(double x) {
  int k;
  double t, y, v;
  static double a[16] = {
    1.5151605362537935201e-13, 3.363790951353651035e-11,
    5.6514041131016827202e-9, 6.7816840255069534052e-7,
    5.4253472222259226487e-5, 0.0026041666666666637057,
    0.06250000000000000009, 0.5,
    -8.9790303384748696588e-11, -1.4029047449249185771e-8,
    -1.5592893752540998113e-6, -1.1253607018469017569e-4,
    -0.0046421827665011579173, -0.085370719728648409609,
    -0.3079657578292062966, 1.0000000000000000004};
  static double b[120] = {
    -9.4055461896630579928e-12, 3.1307934665844664773e-11,
    4.2005295001519243251e-10, -4.1636196779679820012e-9,
    1.4483026181700966164e-8, 1.1661000205428816914e-8,
    -3.5023996724943046209e-7, 1.4404279316339005012e-6,
    5.358156415715824208e-7, -3.5249754038612334639e-5,
    1.7150324075631641453e-4, -4.1276362611239191024e-5,
    -0.0046943110979636602591, 0.035085369853392357659,
    0.20063574339907819159,
    3.3998989888944034586e-11, 7.1558979072937373055e-11,
    -2.9226856932927698732e-9, 1.4591620256525610213e-8,
    -6.6141635609854161666e-9, -1.9991101838984472332e-7,
    5.9185836628873530572e-7, 1.9562880347358085687e-6,
    -1.5814366450418102764e-5, 7.6791682910944612028e-6,
    2.8354678948323983936e-4, -0.0010217932669418690641,
    -0.0032205661865210048433, 0.043497494842354644077,
    0.16110284302315089935,
    -2.0933987679665541827e-10, 7.9503322090520447316e-10,
    3.8000150948242575774e-9, -2.3076136195585571309e-8,
    -2.3902880302550799653e-8, 3.1914500937804377478e-7,
    3.2639909831082417694e-7, -5.3166994792995439449e-6,
    -3.1109524694269240094e-6, 9.2575906966353273247e-5,
    7.5063709094147644746e-7, -0.0017416491592625765379,
    0.0012138560335171676007, 0.045879687144659643175,
    0.11566544716132846709,
    3.1582384905164908749e-10, -1.9959561818098999516e-9,
    8.6959328920030927557e-10, 1.1642778282445577109e-8,
    4.3552264337818440471e-8, -1.5057982160481803238e-7,
    -1.0101729117980989857e-6, 7.7002002510177612013e-7,
    1.9580574235590194233e-5, 1.9358461980242834361e-5,
    -3.3932339942485532728e-4, -9.3416673584325090073e-4,
    0.0055800080455912847227, 0.038668683062477179235,
    0.072651643500517000658,
    -1.1554749629758510059e-10, 8.2270678758893273006e-10,
    -5.0211156951551538591e-10, -1.4929179050554858361e-9,
    -2.7107940791526366702e-8, -4.2204764086705349384e-8,
    3.7253098167927628867e-7, 2.4374697215363361156e-6,
    1.414194200690976837e-6, -4.8766389019473918231e-5,
    -2.1681387247526720978e-4, 2.9325729929653405236e-4,
    0.0064087534504827239815, 0.026054628289709454356,
    0.040156431128194184336,
    2.5506555170746221691e-11, -1.3521164018407978152e-10,
    -8.3281235274106699399e-11, -9.7764849575562351891e-10,
    3.4661828409940354542e-9, 3.9760633711791357544e-8,
    1.590290664550452993e-7, -1.4919441249454941275e-7,
    -5.3779684992094351263e-6, -2.7513862296246223142e-5,
    -9.7880089725297162007e-6, 7.0787668964515789714e-4,
    0.0046968199862345387583, 0.014745740181663320127,
    0.020048622219583455723,
    -3.4824483072529265585e-12, 1.5157161810563380451e-12,
    8.5303859696700686144e-12, 3.3455414203743741076e-10,
    2.0226016353844285376e-9, 5.312815400326633499e-9,
    -3.0799322316418042137e-8, -4.4455408272954712128e-7,
    -2.4293274626893384034e-6, -3.2129079340119038354e-6,
    5.922540368307538885e-5, 5.6822962576781683532e-4,
    0.0027152446516406682732, 0.0074075873691848838485,
    0.0093044450815739269849,
    -2.7683216166276377232e-13, 3.1986676777610155465e-12,
    9.4142986954031445666e-12, 6.7934609179456399334e-11,
    3.485152941147002933e-11, -2.5785248508896551557e-9,
    -2.8310220027112571258e-8, -1.6384131113072271115e-7,
    -3.2521663350596379097e-7, 4.038138875762230716e-6,
    5.1917606978077281001e-5, 3.3420027947470126154e-4,
    0.0013699550623118247094, 0.0034405619148342271096,
    0.0041042919106665762794};
  static double c[120] = {
    4.5281968025889407937e-12, 1.0806749918195271176e-11,
    9.6200972728717669027e-11, 5.721422706362526365e-10,
    3.6077804282954825099e-9, 2.2465236858536681852e-8,
    1.367696126430873523e-7, 7.9561767489531997361e-7,
    4.3014380065615550573e-6, 2.092171390555028559e-5,
    8.8079183950590176926e-5, 3.0549414408830252064e-4,
    8.1295715613927890473e-4, 0.0014679809476357079195,
    0.0013439197177355090057,
    7.6019964430402432637e-13, -2.261619859915827119e-13,
    1.7904450823779000744e-11, 9.1467054855312232717e-11,
    7.1378582044879519122e-10, 4.9925255415445769102e-9,
    3.3767315471315546644e-8, 2.1350774539167751457e-7,
    1.2314353082655232903e-6, 6.2918685053670619181e-6,
    2.7493229298777000013e-5, 9.8085825401369821771e-5,
    2.6670282677770444935e-4, 4.8967895428135985381e-4,
    4.5418248688489697144e-4,
    9.4180115230375147213e-14, 7.5943117003734061145e-14,
    3.0335730243874287654e-12, 2.0202796115462268051e-11,
    1.6839020189186971198e-10, 1.2907875663127201526e-9,
    9.354767612586579892e-9, 6.2471974110281880722e-8,
    3.7585985422997380441e-7, 1.9838348288114906484e-6,
    8.8884862203671982034e-6, 3.2333259238682810218e-5,
    8.9266668913380400243e-5, 1.6589185669844051903e-4,
    1.5536921180500113394e-4,
    1.5425475332301107271e-14, 2.8674534590132490434e-14,
    6.5078462279160216936e-13, 5.0939757793961391211e-12,
    4.497983746074897552e-11, 3.6662925847520171711e-10,
    2.7848878755089582413e-9, 1.929812005933947782e-8,
    1.1950323861976892013e-7, 6.4513432758147478287e-7,
    2.9422095033982461936e-6, 1.0854433321174584937e-5,
    3.0307433185818899481e-5, 5.684098144306501785e-5,
    5.3637016379451945253e-5,
    3.1077953698439839352e-15, 8.6899496170729520378e-15,
    1.6258562067326054104e-13, 1.4104842571366761537e-12,
    1.3019455544084110747e-11, 1.1070466372863950239e-10,
    8.6890603844230597917e-10, 6.1793722175049967488e-9,
    3.9058865943755615801e-8, 2.1432806981070368523e-7,
    9.9034657762983230155e-7, 3.6925185861895664251e-6,
    1.0399877577259449786e-5, 1.9644939661550210015e-5,
    1.8648773453825584597e-5,
    7.2831555285336286457e-16, 2.6077534095895783532e-15,
    4.4881202059263153495e-14, 4.1674329383944385626e-13,
    3.9760100480223728037e-12, 3.483597635535118301e-11,
    2.79932542127702497e-10, 2.0286513276830758107e-9,
    1.3018343087118439152e-8, 7.2315927974997999365e-8,
    3.3750708681924201599e-7, 1.2688020879407355571e-6,
    3.5980954090811587848e-6, 6.8358260635246667316e-6,
    6.5208606745808860557e-6,
    1.9026412343503745875e-16, 8.0073765508732553766e-16,
    1.3245754278055523992e-14, 1.2885201653055058502e-13,
    1.2600129301230402587e-12, 1.1283306843147549277e-11,
    9.2261481309646814329e-11, 6.7812033168299846818e-10,
    4.4020645304595102132e-9, 2.4685719238301517679e-8,
    1.1611886719473112509e-7, 4.3940380936523135466e-7,
    1.2529878285546791905e-6, 2.3917218527087570384e-6,
    2.290757464767187816e-6,
    5.3709522135744366512e-17, 2.5239221050372845433e-16,
    4.093314714589908336e-15, 4.1152784247617592367e-14,
    4.0998840572769381012e-13, 3.7319354625807158852e-12,
    3.0921671702920868014e-11, 2.2975898538634445343e-10,
    1.5049754445782364328e-9, 8.5030864719789148982e-9,
    4.025055939111842381e-8, 1.5312755642491878591e-7,
    4.3865020375297892208e-7, 8.4059737392822153101e-7,
    8.0785884122023473319e-7};
  static double d[40] = {
    9.2371554649979581914e-14, -2.3111336195788410887e-12,
    5.7728710326649832559e-11, -1.8002298130091372598e-9,
    7.6810375010517145638e-8, -5.2669973752193823306e-6,
    0.0010112634961227401357, 0.16180215937964160466,
    6.1381146507252683381e-14, -2.1034499679806301862e-12,
    5.7090233460448415278e-11, -1.7990724350642330817e-9,
    7.6809056078388019946e-8, -5.2669964425290062357e-6,
    0.001011263495747828339, 0.16180215937970716383,
    4.2458150578401296419e-14, -1.8435733128339016981e-12,
    5.5534955081564656595e-11, -1.7938162188526358466e-9,
    7.6798230945934117807e-8, -5.2669828728791921259e-6,
    0.0010112634861753356559, 0.16180215938263409582,
    3.0314798962267007518e-14, -1.5915009905364214455e-12,
    5.3275907427402047438e-11, -1.7824862013841369751e-9,
    7.676389035644707581e-8, -5.2669199860465945909e-6,
    0.0010112634217687349189, 0.16180215941108227283,
    2.2211515002229271212e-14, -1.3664088221521734796e-12,
    5.0585177270502341602e-11, -1.7645432205894533462e-9,
    7.6691805594577373698e-8, -5.2667455286976269634e-6,
    0.001011263186281097458, 0.16180215954783127877};
 
  if (x < 0.8) {
    t = x * x;
    y = (((((((a[0] * t + a[1]) * t +
              a[2]) *
               t +
             a[3]) *
              t +
            a[4]) *
             t +
           a[5]) *
            t +
          a[6]) *
           t +
         a[7]) *
        x;
    y = (((((((a[8] * t + a[9]) * t +
              a[10]) *
               t +
             a[11]) *
              t +
            a[12]) *
             t +
           a[13]) *
            t +
          a[14]) *
           t +
         a[15]) /
          x +
        y * log(x);
  } else if (x < 5.5) {
    v = 3 / x;
    t = x - v;
    k = (int)(t + 3);
    t = (k - 2) - t;
    k *= 15;
    y = ((((((((((((((b[k] * t + b[k + 1]) * t +
                     b[k + 2]) *
                      t +
                    b[k + 3]) *
                     t +
                   b[k + 4]) *
                    t +
                  b[k + 5]) *
                   t +
                 b[k + 6]) *
                  t +
                b[k + 7]) *
                 t +
               b[k + 8]) *
                t +
              b[k + 9]) *
               t +
             b[k + 10]) *
              t +
            b[k + 11]) *
             t +
           b[k + 12]) *
            t +
          b[k + 13]) *
           t +
         b[k + 14]) *
        v;
  } else if (x < 12.5) {
    k = (int)x;
    t = (k + 1) - x;
    k = 15 * (k - 5);
    y = (((((((((((((c[k] * t + c[k + 1]) * t +
                    c[k + 2]) *
                     t +
                   c[k + 3]) *
                    t +
                  c[k + 4]) *
                   t +
                 c[k + 5]) *
                  t +
                c[k + 6]) *
                 t +
               c[k + 7]) *
                t +
              c[k + 8]) *
               t +
             c[k + 9]) *
              t +
            c[k + 10]) *
             t +
           c[k + 11]) *
            t +
          c[k + 12]) *
           t +
         c[k + 13]) *
          t +
        c[k + 14];
  } else {
    t = 60 / x;
    k = 8 * ((int)t);
    y = (((((((d[k] * t + d[k + 1]) * t +
              d[k + 2]) *
               t +
             d[k + 3]) *
              t +
            d[k + 4]) *
             t +
           d[k + 5]) *
            t +
          d[k + 6]) *
           t +
         d[k + 7]) *
        sqrt(t) * exp(-x);
  }
  return y;
}
 
/**
 * @brief Computes the regular Bessel function of the second kind Y0(x) in double precision.
 * 
 * @param x The input value.
 * @return The computed value of Y0(x).
 */
double dbesy0(double x) {
  int k;
  double t, y, v, theta;
  static double a[16] = {
    -1.51249795e-12, 2.9979612902e-10,
    -4.317352912436e-8, 4.31735413787068e-6,
    -2.763106650893309e-4, 0.0099471839432433894,
    -0.15915494309189533339, 0.63661977236758134306,
    4.09490035e-12, -7.6925095943e-10,
    1.0358472550303e-7, -9.49500519343105e-6,
    5.3860266685948738e-4, -0.01607396802593822992,
    0.17760601686906713536, -0.07380429510868722524};
  static double b[112] = {
    -2.958527319e-11, -4.799424787e-11,
    7.0761365903e-10, 7.86916310954e-9,
    4.40632174633e-8, 1.047420705431e-7,
    -7.0359581447993e-7, -1.068225180236166e-5,
    -7.636033201484949e-5, -3.4318294473105478e-4,
    -5.0019381740765567e-4, 0.00898034680285547482,
    0.14783488851798565262, 0.01218096458136948754,
    -3.790394288e-11, -8.3012324083e-10,
    -4.16344143065e-9, -1.59750317942e-9,
    1.1897479804282e-7, 8.9462034665663e-7,
    2.7927785508031e-6, -7.10267818373816e-6,
    -1.4158631248747538e-4, -8.9512296835766028e-4,
    -0.00286419212921734238, 0.00448696141199596815,
    0.16247276853497028694, 0.16806523426268420517,
    4.9508820408e-10, 2.40256475451e-9,
    -3.4128702476e-9, -7.926555593644e-8,
    -2.7510687602104e-7, 8.2850014912691e-7,
    1.202838962100801e-5, 4.446440998573947e-5,
    -6.345799532264893e-5, -0.00153636349657759957,
    -0.00784582749467635199, -0.01091113270659651752,
    0.15855078058220179592, 0.33112075987570150227,
    -1.21381163034e-9, -1.85252445637e-9,
    2.782181715752e-8, 9.983716245339e-8,
    -5.1382328105461e-7, -4.13045112612025e-6,
    6.9219807212556e-7, 1.1346881914173059e-4,
    4.485714262231131e-4, -7.5862241142011713e-4,
    -0.01330170710405210993, -0.04337207915110066583,
    0.10708263075506099964, 0.46937172516119821029,
    1.23954313286e-9, -1.3150647347e-9,
    -3.208158247443e-8, 6.936208476862e-8,
    9.6326321605564e-7, -1.86813677726279e-6,
    -3.296596487362508e-5, -1.027532398838717e-5,
    8.8296805857821394e-4, 0.00288502916411875494,
    -0.00981551522729451135, -0.08175499367255420548,
    -0.0197090243242566078, 0.51958016003260732459,
    -5.0293922567e-10, 2.5905819115e-9,
    3.80069015768e-9, -1.257928760583e-7,
    2.9666593806065e-7, 5.32199410651341e-6,
    -1.466013674557952e-5, -2.1196179800800308e-4,
    1.488907515797525e-4, 0.00598670133725335804,
    0.00917434640154761191, -0.08592436550224288681,
    -0.1970089491453584694, 0.41202451505149651586,
    2.689416111e-11, -6.0225798842e-10,
    7.16632287642e-9, -3.038834206093e-8,
    -4.9829955736697e-7, 3.6467866107985e-6,
    2.599276311771812e-5, -1.6374249229498868e-4,
    -0.00112385259950125883, 0.00342431039748960465,
    0.03016663018029814938, -0.02432512430353342927,
    -0.31797371916576712163, 0.14418001571733803493,
    1.9463221e-12, -1.3451327719e-10,
    3.36581427068e-9, 2.163920016772e-8,
    -5.0868502638387e-7, -1.26562372408989e-6,
    3.554054033687039e-5, 7.490048474688881e-5,
    -0.0014238004200886133, -0.00355314274445060456,
    0.03042001885994619102, 0.07365497405664882878,
    -0.26882214651298246524, -0.16578636480856581212};
  static double c[126] = {
    1.4133441141e-10, -1.16239182574e-9,
    5.60813675548e-9, 4.361111919906e-8,
    -1.9980743680592e-7, -6.23828267379045e-6,
    2.659779111113435e-5, 3.4429768949740887e-4,
    -0.00128058451746064875, -0.01165906901727504569,
    0.03800023545011044967, 0.13079665132249175616,
    -0.30099732306965462304, -0.19470500862950453349,
    2.985070388e-11, -4.3935370028e-10,
    -2.70094988368e-9, 5.27013936718e-8,
    3.5596645009533e-7, -5.5939305510166e-6,
    -2.415386269153308e-5, 3.4989509210648774e-4,
    9.8204687587121096e-4, -0.01241793885175385487,
    -0.0139851984097438996, 0.16758045653582919005,
    0.02375823895638961834, -0.33948059288191103829,
    3.440543098e-11, -1.492311224e-11,
    -5.50307983816e-9, 2.88861169685e-9,
    6.5975576340131e-7, -6.4253433004282e-7,
    -5.095354686815779e-5, 6.35545499385876e-5,
    0.00231750171436134067, -0.00344140664586595266,
    -0.04796153689875432703, 0.06553727330877767204,
    0.27409127395927545297, -0.17324243491898233567,
    1.869953108e-11, 3.5678205287e-10,
    -3.23369443017e-9, -4.932837601739e-8,
    4.0636932259569e-7, 4.55294558480472e-6,
    -3.376596464782138e-5, -2.5758322585307069e-4,
    0.00167416927479443727, 0.00735300543429220468,
    -0.03904554680817849396, -0.07593187710651205994,
    0.25912851048611625179, 0.11731328614820863082,
    -1.168499781e-11, 4.049157062e-10,
    1.76042185153e-9, -5.796810963324e-8,
    -1.7664264701199e-7, 5.65232278687159e-6,
    1.060939366068237e-5, -3.4382320897779425e-4,
    -2.8786871040915476e-4, 0.0110369660250104046,
    0.00105742483935856071, -0.13664188676516064192,
    0.02616867939853747003, 0.27020510536578747599,
    -3.09714065e-11, 1.0275535029e-10,
    5.07546062007e-9, -1.697584561546e-8,
    -5.8222618994516e-7, 1.92483899143098e-6,
    4.3388818054202e-5, -1.3714074807064326e-4,
    -0.00184722933402281351, 0.0051736071998208299,
    0.03611657815299529568, -0.07491163418624572802,
    -0.20317989938720766824, 0.17121062620272384486,
    -2.289344046e-11, -2.7762757821e-10,
    3.9133021423e-9, 3.673304032817e-8,
    -4.7258696909223e-7, -3.23109113206792e-6,
    3.749284946160478e-5, 1.7039067180362504e-4,
    -0.00172637004639815403, -0.00454152531086626043,
    0.03717220530377418124, 0.04489395902785574534,
    -0.23370422835726857839, -0.06753037249787639679,
    4.73299038e-12, -4.0489249115e-10,
    -5.6747261733e-10, 5.656135976496e-8,
    3.564150059329e-8, -5.34821220288168e-6,
    7.581385461348e-8, 3.1190256463318867e-4,
    -1.4634564376714538e-4, -0.0095821355155047985,
    0.00624681472076521329, 0.11513529702572439703,
    -0.05794254714300082167, -0.22523211169118786537,
    2.732444065e-11, -1.7726475495e-10,
    -4.37562701999e-9, 2.68153688556e-8,
    4.8799100503315e-7, -2.76490187540644e-6,
    -3.513357925824063e-5, 1.7693155138571443e-4,
    0.00144533233498513085, -0.00599106092630544975,
    -0.02759437856689911694, 0.07945362316083459396,
    0.15383825653750118008, -0.17121430684466928734};
  static double d[52] = {
    1.059601355592185731e-14, -2.71150591218550377e-13,
    8.6514809056201638e-12, -4.6264028554286627e-10,
    5.0815403835647104e-8, -1.76722552048141208e-5,
    0.16286750396763997378, 2.949651820598278873e-13,
    -8.818215611676125741e-12, 3.571119876162253451e-10,
    -2.63192412099371706e-8, 4.709502795656698909e-6,
    -0.005208333333333283282,
    7.18344107717531977e-15, -2.51623725588410308e-13,
    8.6017784918920604e-12, -4.6256876614290359e-10,
    5.0815343220437937e-8, -1.7672255176494197e-5,
    0.16286750396763433767, 2.2327570859680094777e-13,
    -8.464594853517051292e-12, 3.563766464349055183e-10,
    -2.631843986737892965e-8, 4.70950234228865941e-6,
    -0.0052083333332278466225,
    5.15413392842889366e-15, -2.27740238380640162e-13,
    8.4827767197609014e-12, -4.6224753682737618e-10,
    5.0814848128929134e-8, -1.7672254763876748e-5,
    0.16286750396748926663, 1.7316195320192170887e-13,
    -7.971122772293919646e-12, 3.544039469911895749e-10,
    -2.631443902081701081e-8, 4.709498228695400603e-6,
    -0.005208333331514365361,
    3.84653681453798517e-15, -2.04464520778789011e-13,
    8.3089298605177838e-12, -4.6155016158412096e-10,
    5.081326369646665e-8, -1.76722528311426167e-5,
    0.1628675039665006593, 1.3797879972460878797e-13,
    -7.448089381011684812e-12, 3.51273379710695978e-10,
    -2.630500895563592722e-8, 4.709483934775839193e-6,
    -0.0052083333227940760113};
 
  if (x < 0.85) {
    t = x * x;
    y = ((((((a[0] * t + a[1]) * t +
             a[2]) *
              t +
            a[3]) *
             t +
           a[4]) *
            t +
          a[5]) *
           t +
         a[6]) *
          t +
        a[7];
    y = ((((((a[8] * t + a[9]) * t +
             a[10]) *
              t +
            a[11]) *
             t +
           a[12]) *
            t +
          a[13]) *
           t +
         a[14]) *
          t +
        a[15] + y * log(x);
  } else if (x < 4.5) {
    t = x - 4 / x;
    k = (int)(t + 4);
    t -= k - 3.5;
    k *= 14;
    y = ((((((((((((b[k] * t + b[k + 1]) * t +
                   b[k + 2]) *
                    t +
                  b[k + 3]) *
                   t +
                 b[k + 4]) *
                  t +
                b[k + 5]) *
                 t +
               b[k + 6]) *
                t +
              b[k + 7]) *
               t +
             b[k + 8]) *
              t +
            b[k + 9]) *
             t +
           b[k + 10]) *
            t +
          b[k + 11]) *
           t +
         b[k + 12]) *
          t +
        b[k + 13];
  } else if (x < 12.5) {
    k = (int)x;
    t = x - (k + 0.5);
    k = 14 * (k - 4);
    y = ((((((((((((c[k] * t + c[k + 1]) * t +
                   c[k + 2]) *
                    t +
                  c[k + 3]) *
                   t +
                 c[k + 4]) *
                  t +
                c[k + 5]) *
                 t +
               c[k + 6]) *
                t +
              c[k + 7]) *
               t +
             c[k + 8]) *
              t +
            c[k + 9]) *
             t +
           c[k + 10]) *
            t +
          c[k + 11]) *
           t +
         c[k + 12]) *
          t +
        c[k + 13];
  } else {
    v = 24 / x;
    t = v * v;
    k = 13 * ((int)t);
    y = ((((((d[k] * t + d[k + 1]) * t +
             d[k + 2]) *
              t +
            d[k + 3]) *
             t +
           d[k + 4]) *
            t +
          d[k + 5]) *
           t +
         d[k + 6]) *
        sqrt(v);
    theta = (((((d[k + 7] * t + d[k + 8]) * t +
                d[k + 9]) *
                 t +
               d[k + 10]) *
                t +
              d[k + 11]) *
               t +
             d[k + 12]) *
              v -
            0.78539816339744830962;
    y *= sin(x + theta);
  }
  return y;
}
 
/**
 * @brief Computes the regular Bessel function of the second kind Y1(x) in double precision.
 * 
 * @param x The input value.
 * @return The computed value of Y1(x).
 */
double dbesy1(double x) {
  int k;
  double t, y, v, theta;
  static double a[16] = {
    -9.46499e-14, 2.141432796e-11,
    -3.5977944347e-9, 4.3173541397185e-7,
    -3.453883313621946e-5, 0.00165786399054057257,
    -0.03978873577297383381, 0.31830988618379067154,
    -5.571931463e-11, 8.93098810365e-9,
    -9.9267351748541e-7, 7.164268731546281e-5,
    -0.00295530533604595764, 0.05434868816050718549,
    -0.1960570906462388432, -0.63661977236758134367};
  static double b[112] = {
    2.553685555e-11, 7.384080449e-11,
    -2.9605446927e-10, -5.2081739106e-9,
    -3.676554001376e-8, -1.532652533627e-7,
    -1.3967947147386e-7, 3.88730181348638e-6,
    4.052240780599446e-5, 2.5133005806772436e-4,
    0.00109380629039004777, 0.00298157276348996763,
    4.7188983238582062e-4, -0.19643830264444486066,
    1.056916126e-11, 6.2160122556e-10,
    3.85741197273e-9, 7.491469406e-9,
    -6.306874425799e-8, -6.869833657799e-7,
    -3.29785491054157e-6, -5.57893716372629e-6,
    4.634569030169549e-5, 4.9086790937752985e-4,
    0.00256488729380268318, 0.00822571374211646414,
    0.01094512145046992425, -0.19159562910398869317,
    -4.0211480479e-10, -2.29721816531e-9,
    6.2931855684e-10, 6.509879184236e-8,
    3.1589903600713e-7, -3.093448289001e-8,
    -8.45532004628483e-6, -4.930719004112743e-5,
    -9.865914153735679e-5, 4.7059803263198826e-4,
    0.00472271367216891688, 0.01915518502418162917,
    0.0372239136138415874, -0.1693323078896687811,
    1.10728466147e-9, 2.05608441891e-9,
    -2.536199750468e-8, -1.1751665325851e-7,
    3.404534385077e-7, 4.36083390507752e-6,
    8.34725302600388e-6, -7.154971448642437e-5,
    -5.1956165332380432e-4, -0.00100811886858456098,
    0.00436854257164588225, 0.03425644258611852849,
    0.0907397117201493499, -0.10791597118914788657,
    -1.17749768422e-9, 1.68403859797e-9,
    3.471030171126e-8, -6.227816248757e-8,
    -1.18405395385178e-6, 1.4478514453927e-7,
    3.630744090017806e-5, 1.0481065165516454e-4,
    -5.2483993430622472e-4, -0.00407014291708523893,
    -0.00574311514714191566, 0.03534768318503458091,
    0.16539360869993020163, 0.01986240425870529979,
    4.008050552e-10, -3.36044652853e-9,
    -9.2387130705e-10, 1.7840823064933e-7,
    -2.2643988836534e-7, -8.08896250262948e-6,
    -1.69403516892975e-6, 2.6617475073041589e-4,
    7.2562705080786576e-4, -0.00399208199997748377,
    -0.02399816920965700841, -0.0092557474455754244,
    0.2008281402282066699, 0.21040564991987672576,
    6.692975871e-11, 1.672202705e-10,
    -1.433097981251e-8, 5.455077276306e-8,
    1.04954046999028e-6, -3.39572219821681e-6,
    -5.541577161092338e-5, 4.344113487321887e-5,
    0.00184798594297387512, 0.0030104150846815045,
    -0.02789739141907399514, -0.09421628519129373841,
    0.09950212445529771576, 0.3749697583874640408,
    4.528583204e-11, 6.6861092459e-10,
    -9.2464816857e-9, -8.002058622489e-8,
    8.7578676979688e-7, 6.28456150216647e-6,
    -4.276239654834506e-5, -3.4630842915911896e-4,
    8.9431202169021474e-4, 0.01086380107342193819,
    0.00145563493724746238, -0.14193383817937981478,
    -0.15148562296986494725, 0.35720232689066512346};
  static double c[126] = {
    3.201620498e-10, -1.77431944127e-9,
    1.367520951871e-8, -6.171755212447e-8,
    -4.3602779932091e-7, 1.79827173442894e-6,
    4.990624965075363e-5, -1.8618453816525107e-4,
    -0.00206578613620093339, 0.00640292258731756063,
    0.04663627606907821596, -0.11400070635033154332,
    -0.26159330264498333979, 0.30099732306965462346,
    -5.9774708e-13, -3.8827648808e-10,
    5.27269271362e-9, 2.971059814012e-8,
    -5.2701406514027e-7, -3.20369808906187e-6,
    4.475144442564492e-5, 1.6907703884518866e-4,
    -0.00209937055264007564, -0.00491023437935628693,
    0.0496717554070154514, 0.04195559522923170315,
    -0.33516091307165838036, -0.02375823895638961835,
    2.83380103e-12, -4.4715450235e-10,
    1.769519961e-10, 6.05337984007e-8,
    -2.888550814408e-8, -5.93780185020349e-6,
    5.14027455732102e-6, 3.5667482807472357e-4,
    -3.813272996260773e-4, -0.01158750857180657934,
    0.01376562658346365929, 0.14388461069626297876,
    -0.1310745466175553429, -0.27409127395927545296,
    2.683697058e-11, -2.4302161327e-10,
    -4.30151236236e-9, 3.557058903199e-8,
    4.9328952592926e-7, -3.65732389065382e-6,
    -3.642356546467712e-5, 2.3636175253326716e-4,
    0.00154549935517002108, -0.00837084637397210913,
    -0.02941202173717025197, 0.11713664042453548043,
    0.15186375421302413107, -0.25912851048611625179,
    2.914434604e-11, 1.5185121893e-10,
    -4.88084673389e-9, -1.936460341189e-8,
    5.7968735781303e-7, 1.58978381365925e-6,
    -4.521858314881102e-5, -7.426575562367423e-5,
    0.00206293925392279862, 0.00143934355204571641,
    -0.04414786410004317489, -0.00317227451807568106,
    0.273283773530321296, -0.02616867939853747003,
    6.44550655e-12, 4.0249587279e-10,
    -1.23789833341e-9, -5.58299757877e-8,
    1.6975984093144e-7, 5.24003568623094e-6,
    -1.539871212028104e-5, -3.0372172637669852e-4,
    8.2284448843625174e-4, 0.00923614667011392611,
    -0.02069442879928366384, -0.10834973445898588438,
    0.14982326837249145873, 0.20317989938720766823,
    -2.115691215e-11, 2.9752035314e-10,
    3.34739862269e-9, -4.304625868398e-8,
    -3.6733494871207e-7, 4.25328270524107e-6,
    2.584872967637478e-5, -2.6244994622929809e-4,
    -0.00102234403086242668, 0.00863185023199066935,
    0.01816610124346617164, -0.11151661591132254182,
    -0.0897879180557114995, 0.23370422835726857838,
    -2.957977086e-11, -6.150420279e-11,
    4.88089472193e-9, 6.24218182859e-9,
    -5.6561995267846e-7, -3.2077350100271e-7,
    4.27856984896483e-5, -5.3069698280036e-7,
    -0.0018714153878560023, 7.3172821883575326e-4,
    0.03832854206202077374, -0.01874044416229564036,
    -0.2302705940514488064, 0.05794254714300082167,
    -1.203179184e-11, -3.5509887475e-10,
    2.13620090322e-9, 4.813181550803e-8,
    -2.6815627351126e-7, -4.39191902440609e-6,
    2.211921535574542e-5, 2.45935054805247e-4,
    -0.00106158930833741898, -0.00722666167492552731,
    0.02396424370522244155, 0.08278313570069734843,
    -0.15890724632166919294, -0.15383825653750118007};
  static double d[52] = {
    -1.272346002224188092e-14, 3.370464692346669075e-13,
    -1.144940314335484869e-11, 6.863141561083429745e-10,
    -9.491933932960924159e-8, 5.301676561445687562e-5,
    0.162867503967639974, -3.652982212914147794e-13,
    1.151126750560028914e-11, -5.165585095674343486e-10,
    4.657991250060549892e-8, -1.186794704692706504e-5,
    0.01562499999999994026,
    -8.713069680903981555e-15, 3.140780373478474935e-13,
    -1.139089186076256597e-11, 6.862299023338785566e-10,
    -9.491926788274594674e-8, 5.301676558106268323e-5,
    0.162867503967646622, -2.792555727162752006e-13,
    1.108650207651756807e-11, -5.156745588549830981e-10,
    4.657894859077370979e-8, -1.186794650130550256e-5,
    0.01562499999987299901,
    -6.304859171204770696e-15, 2.857249044208791652e-13,
    -1.124956921556753188e-11, 6.858482894906716661e-10,
    -9.49186795351689846e-8, 5.301676509057781574e-5,
    0.1628675039678191167, -2.185193490132496053e-13,
    1.048820673697426074e-11, -5.132819367467680132e-10,
    4.65740943737299422e-8, -1.186794150862988921e-5,
    0.01562499999779270706,
    -4.74041720979200985e-15, 2.578715253644144182e-13,
    -1.104148898414138857e-11, 6.850134201626289183e-10,
    -9.49167823417491964e-8, 5.301676277588728159e-5,
    0.1628675039690033136, -1.75512205749384229e-13,
    9.848723331445182397e-12, -5.094535425482245697e-10,
    4.656255982268609304e-8, -1.186792402114394891e-5,
    0.01562499998712198636};
 
  if (x < 0.85) {
    t = x * x;
    y = (((((((a[0] * t + a[1]) * t +
              a[2]) *
               t +
             a[3]) *
              t +
            a[4]) *
             t +
           a[5]) *
            t +
          a[6]) *
           t +
         a[7]) *
        x;
    y = (((((((a[8] * t + a[9]) * t +
              a[10]) *
               t +
             a[11]) *
              t +
            a[12]) *
             t +
           a[13]) *
            t +
          a[14]) *
           t +
         a[15]) /
          x +
        y * log(x);
  } else if (x < 4.15) {
    v = 4 / x;
    t = x - v;
    k = (int)(t + 4);
    t -= k - 3.5;
    k *= 14;
    y = (((((((((((((b[k] * t + b[k + 1]) * t +
                    b[k + 2]) *
                     t +
                   b[k + 3]) *
                    t +
                  b[k + 4]) *
                   t +
                 b[k + 5]) *
                  t +
                b[k + 6]) *
                 t +
               b[k + 7]) *
                t +
              b[k + 8]) *
               t +
             b[k + 9]) *
              t +
            b[k + 10]) *
             t +
           b[k + 11]) *
            t +
          b[k + 12]) *
           t +
         b[k + 13]) *
        v;
  } else if (x < 12.5) {
    k = (int)x;
    t = x - (k + 0.5);
    k = 14 * (k - 4);
    y = ((((((((((((c[k] * t + c[k + 1]) * t +
                   c[k + 2]) *
                    t +
                  c[k + 3]) *
                   t +
                 c[k + 4]) *
                  t +
                c[k + 5]) *
                 t +
               c[k + 6]) *
                t +
              c[k + 7]) *
               t +
             c[k + 8]) *
              t +
            c[k + 9]) *
             t +
           c[k + 10]) *
            t +
          c[k + 11]) *
           t +
         c[k + 12]) *
          t +
        c[k + 13];
  } else {
    v = 24 / x;
    t = v * v;
    k = 13 * ((int)t);
    y = ((((((d[k] * t + d[k + 1]) * t +
             d[k + 2]) *
              t +
            d[k + 3]) *
             t +
           d[k + 4]) *
            t +
          d[k + 5]) *
           t +
         d[k + 6]) *
        sqrt(v);
    theta = (((((d[k + 7] * t + d[k + 8]) * t +
                d[k + 9]) *
                 t +
               d[k + 10]) *
                t +
              d[k + 11]) *
               t +
             d[k + 12]) *
              v -
            0.78539816339744830962;
    y *= -cos(x + theta);
  }
  return y;
}
 
/**
 * @brief Evaluates a Chebyshev series at a given point.
 *
 * This function evaluates the Chebyshev series using the Clenshaw algorithm.
 * 
 * @param a Lower bound of the interval.
 * @param b Upper bound of the interval.
 * @param c[] Array of Chebyshev coefficients.
 * @param m Number of coefficients in the array.
 * @param x The point at which the series is to be evaluated. Must be within [a, b].
 * @return The value of the Chebyshev series at x.
 * @note If x is outside the interval [a, b], the function will terminate with an error message.
 */
double chebev(double a, double b, double c[], int m, double x) {
  double d = 0.0, dd = 0.0, sv, y, y2;
  int j;
 
  if ((x - a) * (x - b) > 0.0) {
    fprintf(stderr, "x not in range in routine chebev\n");
    exit(1);
  }
  y2 = 2.0 * (y = (2.0 * x - a - b) / (b - a));
  for (j = m - 1; j >= 1; j--) {
    sv = d;
    d = y2 * d - dd + c[j];
    dd = sv;
  }
  return y * d - dd + 0.5 * c[0];
}
 
#define NUSE1 5
#define NUSE2 5
 
/**
 * @brief Computes auxiliary Bessel functions \f$\gamma_1\f$, \f$\gamma_2\f$, \f$\gamma_+\f$, and \f$\gamma_-\f$.
 *
 * This function uses precomputed Chebyshev coefficients to calculate the auxiliary Bessel functions.
 * 
 * @param x The input value for which the Bessel functions are computed.
 * @param gam1 Pointer to store the computed value of \f$\gamma_1\f$.
 * @param gam2 Pointer to store the computed value of \f$\gamma_2\f$.
 * @param gampl Pointer to store the computed value of \f$\gamma_+\f$.
 * @param gammi Pointer to store the computed value of \f$\gamma_-\f$.
 */
void beschb(double x, double *gam1, double *gam2, double *gampl, double *gammi) {
  double xx;
  static double c1[] = {
    -1.142022680371172e0, 6.516511267076e-3,
    3.08709017308e-4, -3.470626964e-6, 6.943764e-9,
    3.6780e-11, -1.36e-13};
  static double c2[] = {
    1.843740587300906e0, -0.076852840844786e0,
    1.271927136655e-3, -4.971736704e-6, -3.3126120e-8,
    2.42310e-10, -1.70e-13, -1.0e-15};
 
  xx = 8.0 * x * x - 1.0;
  *gam1 = chebev(-1.0, 1.0, c1, NUSE1, xx);
  *gam2 = chebev(-1.0, 1.0, c2, NUSE2, xx);
  *gampl = *gam2 - x * (*gam1);
  *gammi = *gam2 + x * (*gam1);
}