Changeset 57

Show
Ignore:
Timestamp:
03/13/06 05:36:42 (3 years ago)
Author:
Don Clugston
Message:

Added more test cases to Bessel functions. Changed their names.

Total coverage:
1219 cov, 57 uncov -- 95.5% covered.

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • trunk/bessel.d

    r56 r57  
    11 
     2// Bessel functions are like the planet Mars: 
     3// "I have no use for it -- but It's Nice To Know It's There". 
    24module bessel; 
    35 
    …… 
    911 
    1012import std.math; 
     13 
     14 
     15/* 
     16sqrt(j0^2(1/x^2) + y0^2(1/x^2)) = z P(z**2)/Q(z**2) 
     17z(x) = 1/sqrt(x) 
     18Relative error 
     19n=7, d=7 
     20Peak error =  1.80e-20 
     21Relative error spread =  5.1e-2 
     22*/ 
     23 
     24const real modulusn[] = [ 
     25 
     26   0x1.154700ea96e79656p-7, // 0.0084618334268988678395 
     27   0x1.72244b6e998cd6fp-4,  // 0.090366644531602000525 
     28   0x1.6ebccf42e9c19fd2p-1, // 0.71628425304232057209 
     29   0x1.6bd844e89cbd639ap+1, // 2.8425375114252161456 
     30   0x1.e812b377c75ebc96p+2, // 7.6261414212908496305 
     31   0x1.46d69ca24ce76686p+3, // 10.213697735779743439 
     32   0x1.b756f7234cc67146p+2, // 6.8646829457021346239 
     33   0x1.943a7471eaa50ab2p-2  // 0.39475423760692244614 
     34]; 
     35 
     36const real modulusd[] = [ 
     37   0x1.5b84007c37011506p-7, // 0.010605335461541217703 
     38   0x1.cfe76758639bdab4p-4, // 0.11325779313322123049 
     39   0x1.cbfa09bf71bafc7ep-1, // 0.8983920141407590632 
     40   0x1.c8eafb3836f2eeb4p+1, // 3.5696710609899109018 
     41   0x1.339db78060eb706ep+3, // 9.6130025393862137883 
     42   0x1.a06530916be8bc7ap+3, // 13.012352260614782615 
     43   0x1.23bfe7f67a54893p+3,  // 9.117176038171821116 
     44   0x1p+0   // 1 
     45]; 
     46 
     47/* 
     48atan(y0(x)/j0(x)) = x - pi/4 + x P(x**2)/Q(x**2) 
     49Absolute error 
     50n=5, d=6 
     51Peak error =  2.80e-21 
     52Relative error spread =  5.5e-1 
     53*/ 
     54 
     55const real phasen[] = [ 
     56   -0x1.ccbaf3865bb0985ep-22,   // -4.2908850907731129636e-07 
     57   -0x1.3a6b175e64bdb82ep-14,   // -7.4963170368299641512e-05 
     58   -0x1.06124b5310cdca28p-8,    // -0.0039988931558269906429 
     59   -0x1.3cebb7ab09cf1b14p-4,    // -0.077373235181415168822 
     60   -0x1.00156ed209b43c6p-1, // -0.50016351999224936946 
     61   -0x1.78aa9ba4254ca20cp-1     // -0.73567663553935715194 
     62]; 
     63 
     64const real phased[] = [ 
     65   0x1.ccbaf3865bb09918p-19,    // 3.4327080726184903899e-06 
     66   0x1.3b5b0e12900d58b8p-11,    // 0.00060149323173421904041 
     67   0x1.0897373ff9906f7ep-5, // 0.032298667821850250481 
     68   0x1.450a5b8c552ade4ap-1, // 0.63484464729352451028 
     69   0x1.123e263e7f0e96d2p+2, // 4.2850432977977361182 
     70   0x1.d82ecca5654be7d2p+2, // 7.3778564086143760725 
     71   0x1p+0   // 1 
     72]; 
     73 
     74/*** 
     75 *  Bessel function of order zero 
     76 * 
     77 * Returns Bessel function of first kind, order zero of the argument. 
     78 * 
     79 * The domain is divided into the intervals [0, 9] and 
     80 * (9, infinity). In the first interval the rational approximation 
     81 * is (x^2 - r^2) (x^2 - s^2) (x^2 - t^2) P7(x^2) / Q8(x^2), 
     82 * where r, s, t are the first three zeros of the function. 
     83 * In the second interval the expansion is in terms of the 
     84 * modulus M0(x) = sqrt(J0(x)^2 + Y0(x)^2) and phase  P0(x) 
     85 * = atan(Y0(x)/J0(x)).  M0 is approximated by sqrt(1/x)P7(1/x)/Q7(1/x). 
     86 * The approximation to J0 is M0 * cos(x -  pi/4 + 1/x P5(1/x^2)/Q6(1/x^2)). 
     87 * 
     88 * 
     89 * ACCURACY: 
     90 * 
     91 *                      Absolute error: 
     92 * arithmetic   domain      # trials      peak         rms 
     93 *    IEEE      0, 30       100000      2.8e-19      7.4e-20 
     94 */ 
     95real cylBessel_j0(real x) 
     96{ 
    1197 
    1298/* 
    …… 
    42128]; 
    43129 
    44 /* 
    45 sqrt(j0^2(1/x^2) + y0^2(1/x^2)) = z P(z**2)/Q(z**2) 
    46 z(x) = 1/sqrt(x) 
    47 Relative error 
    48 n=7, d=7 
    49 Peak error =  1.80e-20 
    50 Relative error spread =  5.1e-2 
    51 */ 
    52  
    53  
    54 const real modulusn[] = [ 
    55  
    56    0x1.154700ea96e79656p-7, // 0.0084618334268988678395 
    57    0x1.72244b6e998cd6fp-4,  // 0.090366644531602000525 
    58    0x1.6ebccf42e9c19fd2p-1, // 0.71628425304232057209 
    59    0x1.6bd844e89cbd639ap+1, // 2.8425375114252161456 
    60    0x1.e812b377c75ebc96p+2, // 7.6261414212908496305 
    61    0x1.46d69ca24ce76686p+3, // 10.213697735779743439 
    62    0x1.b756f7234cc67146p+2, // 6.8646829457021346239 
    63    0x1.943a7471eaa50ab2p-2  // 0.39475423760692244614 
    64 ]; 
    65  
    66 const real modulusd[] = [ 
    67    0x1.5b84007c37011506p-7, // 0.010605335461541217703 
    68    0x1.cfe76758639bdab4p-4, // 0.11325779313322123049 
    69    0x1.cbfa09bf71bafc7ep-1, // 0.8983920141407590632 
    70    0x1.c8eafb3836f2eeb4p+1, // 3.5696710609899109018 
    71    0x1.339db78060eb706ep+3, // 9.6130025393862137883 
    72    0x1.a06530916be8bc7ap+3, // 13.012352260614782615 
    73    0x1.23bfe7f67a54893p+3,  // 9.117176038171821116 
    74    0x1p+0   // 1 
    75 ]; 
    76  
    77 /* 
    78 atan(y0(x)/j0(x)) = x - pi/4 + x P(x**2)/Q(x**2) 
    79 Absolute error 
    80 n=5, d=6 
    81 Peak error =  2.80e-21 
    82 Relative error spread =  5.5e-1 
    83 */ 
    84  
    85 const real phasen[] = [ 
    86    -0x1.ccbaf3865bb0985ep-22,   // -4.2908850907731129636e-07 
    87    -0x1.3a6b175e64bdb82ep-14,   // -7.4963170368299641512e-05 
    88    -0x1.06124b5310cdca28p-8,    // -0.0039988931558269906429 
    89    -0x1.3cebb7ab09cf1b14p-4,    // -0.077373235181415168822 
    90    -0x1.00156ed209b43c6p-1, // -0.50016351999224936946 
    91    -0x1.78aa9ba4254ca20cp-1     // -0.73567663553935715194 
    92 ]; 
    93  
    94 const real phased[] = [ 
    95    0x1.ccbaf3865bb09918p-19,    // 3.4327080726184903899e-06 
    96    0x1.3b5b0e12900d58b8p-11,    // 0.00060149323173421904041 
    97    0x1.0897373ff9906f7ep-5, // 0.032298667821850250481 
    98    0x1.450a5b8c552ade4ap-1, // 0.63484464729352451028 
    99    0x1.123e263e7f0e96d2p+2, // 4.2850432977977361182 
    100    0x1.d82ecca5654be7d2p+2, // 7.3778564086143760725 
    101    0x1p+0   // 1 
    102 ]; 
     130 
     131 
     132    real xx, y, z, modulus, phase; 
     133 
     134    xx = x * x; 
     135    if ( xx < 81.0L ) { 
     136        const real [] JZ = [5.783185962946784521176L, 30.47126234366208639908L, 7.488700679069518344489e1L]; 
     137        y = (xx - JZ[0]) * (xx - JZ[1]) * (xx - JZ[2]); 
     138        y *= poly( xx, j0n) / poly( xx, j0d); 
     139        return y; 
     140    } 
     141 
     142    y = fabs(x); 
     143    xx = 1.0/xx; 
     144    phase = poly( xx, phasen) / poly( xx, phased); 
     145 
     146    z = 1.0/y; 
     147    modulus = poly( z, modulusn) / poly( z, modulusd); 
     148 
     149    y = modulus * cos( y -  PI_4 + z*phase) / sqrt(y); 
     150    return y; 
     151} 
     152 
     153/** 
     154 *  Bessel function of the second kind, order zero 
     155 * Also known as the cylindrical Neumann function, order zero. 
     156 * 
     157 * Returns Bessel function of the second kind, of order 
     158 * zero, of the argument. 
     159 * 
     160 * The domain is divided into the intervals [0, 5>, [5,9> and 
     161 * [9, infinity). In the first interval a rational approximation 
     162 * R(x) is employed to compute y0(x)  = R(x) + 2/pi * log(x) * j0(x). 
     163 * 
     164 * In the second interval, the approximation is 
     165 *     (x - p)(x - q)(x - r)(x - s)P7(x)/Q7(x) 
     166 * where p, q, r, s are zeros of y0(x). 
     167 * 
     168 * The third interval uses the same approximations to modulus 
     169 * and phase as j0(x), whence y0(x) = modulus * sin(phase). 
     170 * 
     171 * ACCURACY: 
     172 * 
     173 *  Absolute error, when y0(x) < 1; else relative error: 
     174 * 
     175 * arithmetic   domain     # trials      peak         rms 
     176 *    IEEE      0, 30       100000      3.4e-19     7.6e-20 
     177 */ 
     178real cylBessel_y0(real x) 
     179{ 
    103180 
    104181 
    …… 
    172249 
    173250 
    174 /*** 
    175  *  Bessel function of order zero 
    176  * 
    177  * Returns Bessel function of first kind, order zero of the argument. 
    178  * 
    179  * The domain is divided into the intervals [0, 9] and 
    180  * (9, infinity). In the first interval the rational approximation 
    181  * is (x^2 - r^2) (x^2 - s^2) (x^2 - t^2) P7(x^2) / Q8(x^2), 
    182  * where r, s, t are the first three zeros of the function. 
    183  * In the second interval the expansion is in terms of the 
    184  * modulus M0(x) = sqrt(J0(x)^2 + Y0(x)^2) and phase  P0(x) 
    185  * = atan(Y0(x)/J0(x)).  M0 is approximated by sqrt(1/x)P7(1/x)/Q7(1/x). 
    186  * The approximation to J0 is M0 * cos(x -  pi/4 + 1/x P5(1/x^2)/Q6(1/x^2)). 
    187  * 
    188  * 
    189  * ACCURACY: 
    190  * 
    191  *                      Absolute error: 
    192  * arithmetic   domain      # trials      peak         rms 
    193  *    IEEE      0, 30       100000      2.8e-19      7.4e-20 
    194  */ 
    195 real j0(real x) 
    196 { 
    197  
    198251    real xx, y, z, modulus, phase; 
    199  
    200     xx = x * x; 
    201     if ( xx < 81.0L ) { 
    202         const real JZ1 = 5.783185962946784521176L; 
    203         const real JZ2 = 30.47126234366208639908L; 
    204         const real JZ3 = 7.488700679069518344489e1L; 
    205         y = (xx - JZ1) * (xx - JZ2) * (xx -JZ3); 
    206         y *= poly( xx, j0n) / poly( xx, j0d); 
    207         return y; 
    208     } 
    209  
    210     y = fabs(x); 
    211     xx = 1.0/xx; 
    212     phase = poly( xx, phasen) / poly( xx, phased); 
    213  
    214     z = 1.0/y; 
    215     modulus = poly( z, modulusn) / poly( z, modulusd); 
    216  
    217     y = modulus * cos( y -  PI_4 + z*phase) / sqrt(y); 
    218     return y; 
    219 } 
    220  
    221 /** 
    222  *  Bessel function of the second kind, order zero 
    223  * 
    224  * Returns Bessel function of the second kind, of order 
    225  * zero, of the argument. 
    226  * 
    227  * The domain is divided into the intervals [0, 5>, [5,9> and 
    228  * [9, infinity). In the first interval a rational approximation 
    229  * R(x) is employed to compute y0(x)  = R(x) + 2/pi * log(x) * j0(x). 
    230  * 
    231  * In the second interval, the approximation is 
    232  *     (x - p)(x - q)(x - r)(x - s)P7(x)/Q7(x) 
    233  * where p, q, r, s are zeros of y0(x). 
    234  * 
    235  * The third interval uses the same approximations to modulus 
    236  * and phase as j0(x), whence y0(x) = modulus * sin(phase). 
    237  * 
    238  * ACCURACY: 
    239  * 
    240  *  Absolute error, when y0(x) < 1; else relative error: 
    241  * 
    242  * arithmetic   domain     # trials      peak         rms 
    243  *    IEEE      0, 30       100000      3.4e-19     7.6e-20 
    244  */ 
    245 real y0(real x) 
    246 { 
    247     real xx, y, z, modulus, phase; 
    248  
     252     
    249253    if ( x < 0.0 ) return -real.max; 
    250254    xx = x * x; 
    251255    if ( xx < 81.0L ) { 
    252256        if ( xx < 20.25L ) { 
    253             y = 2*PI * log(x) * j0(x); 
     257            y = M_2_PI * log(x) * cylBessel_j0(x); 
    254258            y += poly( xx, y0n) / poly( xx, y0d); 
    255259        } else { 
    256 const real Y0Z1 = 3.957678419314857868376e0L; 
    257 const real Y0Z2 = 7.086051060301772697624e0L; 
    258 const real Y0Z3 = 1.022234504349641701900e1L; 
    259 const real Y0Z4 = 1.336109747387276347827e1L; 
    260             y = (x - Y0Z1)*(x - Y0Z2)*(x - Y0Z3)*(x - Y0Z4); 
     260            const real [] Y0Z = [3.957678419314857868376e0L, 7.086051060301772697624e0L, 
     261                1.022234504349641701900e1L, 1.336109747387276347827e1L]; 
     262            y = (x - Y0Z[0])*(x - Y0Z[1])*(x - Y0Z[2])*(x - Y0Z[3]); 
    261263            y *= poly( x, y059n) / poly( x, y059d); 
    262264        } 
    …… 
    276278 
    277279unittest { 
    278 // argument, exactval, extra precision, derivative. Use absolute error. 
    279 assert(feqrel(j0(8.0L), 1.71646118164062500000E-1L + 4.68897349140609086941E-6L)>=real.mant_dig-3); 
    280 assert(feqrel(j0(2.0L), 2.23876953125000000000E-1L+  1.38260162356680518275E-5L)>=real.mant_dig-3); 
    281 assert(feqrel(j0(-2.0L), 2.23876953125000000000E-1L+  1.38260162356680518275E-5L)>=real.mant_dig-3); 
    282 assert(feqrel(y0(8.0L), 2.23510742187500000000E-1L + 1.07472000662205273234E-5L)>=real.mant_dig-4); 
    283  
    284  
    285 /+ 
    286 y  {"j0", j0l, 8.0L, 1.71646118164062500000E-1L, 
    287    4.68897349140609086941E-6L, -2.34636346853914624381E-1, -4}, 
    288   {"j0", j0l, 4.54541015625L, -3.09783935546875000000E-1L, 
    289    7.07472668157686463367E-6L, 2.42993657373627558460E-1L, -4}, 
    290   {"j0", j0l, 2.85711669921875L, -2.07901000976562500000E-1L, 
    291    1.15237285263902751582E-5L, -3.90402225324501311651E-1L, -4}, 
    292   {"j0", j0l, 2.0L, 2.23876953125000000000E-1L, 
    293    1.38260162356680518275E-5L, -5.76724807756873387202E-1L, -4}, 
    294   {"j0", j0l, 1.16415321826934814453125e-10L, 9.99984741210937500000E-1L, 
    295    1.52587890624999966119E-5L, 9.99999999999999999997E-1L, -4}, 
    296   {"j0", j0l, -2.0L, 2.23876953125000000000E-1L, 
    297    1.38260162356680518275E-5L, 5.76724807756873387202E-1L, -4}, 
    298 y  {"y0", y0l, 8.0L, 2.23510742187500000000E-1L, 
    299    1.07472000662205273234E-5L, 1.58060461731247494256E-1L, -4}, 
    300   {"y0", y0l, 4.54541015625L, -2.08114624023437500000E-1L, 
    301    1.45018823856668874574E-5L, -2.88887645307401250876E-1L, -4}, 
    302   {"y0", y0l, 2.85711669921875L, 4.20303344726562500000E-1L, 
    303    1.32781607563122276008E-5L, -2.82488638474982469213E-1, -4}, 
    304   {"y0", y0l, 2.0L, 5.10360717773437500000E-1L, 
    305    1.49548763076195966066E-5L, 1.07032431540937546888E-1L, -4}, 
    306   {"y0", y0l, 1.16415321826934814453125e-10L, -1.46357574462890625000E1L, 
    307    3.54110537011061127637E-6L, 5.46852220461145271913E9L, -4}, 
    308   {"j1", j1l, 8.0L, 2.34634399414062500000E-1L, 
    309    1.94743985212438127665E-6L,1.42321263780814578043E-1, -4}, 
    310   {"j1", j1l, 4.54541015625L, -2.42996215820312500000E-1L, 
    311    2.55844668494153980076E-6L, -2.56317734136211337012E-1, -4}, 
    312   {"j1", j1l, 2.85711669921875L, 3.90396118164062500000E-1L, 
    313    6.10716043881165077013E-6L, -3.44531507106757980441E-1L, -4}, 
    314   {"j1", j1l, 2.0L, 5.76721191406250000000E-1L, 
    315    3.61635062338720244824E-6L,  -6.44716247372010255494E-2L, -4}, 
    316   {"j1", j1l, 1.16415321826934814453125e-10L, 
    317    5.820677273504770710133016109466552734375e-11L, 
    318    8.881784197001251337312921818461805735896e-16L, 
    319    4.99999999999999999997E-1L, -4}, 
    320   {"j1", j1l, -2.0L, -5.76721191406250000000E-1L, 
    321    -3.61635062338720244824E-6L,  -6.44716247372010255494E-2L, -4}, 
    322   {"y1", y1l, 8.0L, -1.58065795898437500000E-1L, 
    323    5.33416719000574444473E-6L, 2.43279047103972157309E-1L, -4}, 
    324   {"y1", y1l, 4.54541015625L, 2.88879394531250000000E-1L, 
    325    8.25077615125087585195E-6L, -2.71656024771791736625E-1L, -4}, 
    326   {"y1", y1l, 2.85711669921875L, 2.82485961914062500000E-1, 
    327    2.67656091996921314433E-6L, 3.21444694221532719737E-1, -4}, 
    328   {"y1", y1l, 2.0L, -1.07040405273437500000E-1L, 
    329    7.97373249995311162923E-6L, 5.63891888420213893041E-1, -4}, 
    330   {"y1", y1l, 1.16415321826934814453125e-10L, -5.46852220500000000000E9L,  
    331    3.88547280871200700671E-1L, 4.69742480525120196168E19L, -4}, 
    332 +/    
     280  // argument, result1, result2, derivative. Correct result is result1+result2. 
     281const real [4][] j0_test_points = [ 
     282    [8.0L, 1.71646118164062500000E-1L, 4.68897349140609086941E-6L, -2.34636346853914624381E-1L], 
     283    [4.54541015625L, -3.09783935546875000000E-1L, 7.07472668157686463367E-6L, 2.42993657373627558460E-1L], 
     284    [2.85711669921875L, -2.07901000976562500000E-1L, 1.15237285263902751582E-5L, -3.90402225324501311651E-1L], 
     285    [2.0L, 2.23876953125000000000E-1L, 1.38260162356680518275E-5L, -5.76724807756873387202E-1L], 
     286    [1.16415321826934814453125e-10L, 9.99984741210937500000E-1L, 1.52587890624999966119E-5L, 
     287        9.99999999999999999997E-1L], 
     288    [-2.0L, 2.23876953125000000000E-1L, 
     289        1.38260162356680518275E-5L, 5.76724807756873387202E-1L] 
     290]; 
     291 
     292const real [4][] y0_test_points = [ 
     293    [ 8.0L, 2.23510742187500000000E-1L, 1.07472000662205273234E-5L, 1.58060461731247494256E-1L], 
     294    [4.54541015625L, -2.08114624023437500000E-1L, 1.45018823856668874574E-5L, -2.88887645307401250876E-1L], 
     295    [2.85711669921875L, 4.20303344726562500000E-1L, 1.32781607563122276008E-5L, -2.82488638474982469213E-1], 
     296    [2.0L, 5.10360717773437500000E-1L, 1.49548763076195966066E-5L, 1.07032431540937546888E-1L], 
     297    [1.16415321826934814453125e-10L, -1.46357574462890625000E1L, 3.54110537011061127637E-6L, 
     298        5.46852220461145271913E9L]   
     299]; 
     300 
     301const real [4][] j1_test_points = [ 
     302  [ 8.0L, 2.34634399414062500000E-1L, 1.94743985212438127665E-6L,1.42321263780814578043E-1], 
     303  [4.54541015625L, -2.42996215820312500000E-1L, 2.55844668494153980076E-6L, -2.56317734136211337012E-1], 
     304  [2.85711669921875L, 3.90396118164062500000E-1L, 6.10716043881165077013E-6L, -3.44531507106757980441E-1L], 
     305  [2.0L, 5.76721191406250000000E-1L, 3.61635062338720244824E-6L,  -6.44716247372010255494E-2L], 
     306  [1.16415321826934814453125e-10L, 5.820677273504770710133016109466552734375e-11L, 
     307   8.881784197001251337312921818461805735896e-16L, 4.99999999999999999997E-1L], 
     308  [-2.0L, -5.76721191406250000000E-1L, -3.61635062338720244824E-6L, -6.44716247372010255494E-2L] 
     309]; 
     310 
     311const real [4][] y1_test_points = [ 
     312    [8.0L, -1.58065795898437500000E-1L, 
     313        5.33416719000574444473E-6L, 2.43279047103972157309E-1L], 
     314    [4.54541015625L, 2.88879394531250000000E-1L, 
     315        8.25077615125087585195E-6L, -2.71656024771791736625E-1L], 
     316    [2.85711669921875L, 2.82485961914062500000E-1, 
     317        2.67656091996921314433E-6L, 3.21444694221532719737E-1], 
     318    [2.0L, -1.07040405273437500000E-1L, 
     319        7.97373249995311162923E-6L, 5.63891888420213893041E-1], 
     320    [1.16415321826934814453125e-10L, -5.46852220500000000000E9L,  
     321        3.88547280871200700671E-1L, 4.69742480525120196168E19L] 
     322]; 
     323 
     324    foreach(real [4] t; j0_test_points) { 
     325        assert(feqrel(cylBessel_j0(t[0]), t[1]+t[2]) >=real.mant_dig-3); 
     326    } 
     327 
     328    foreach(real [4] t; y0_test_points) { 
     329        assert(feqrel(cylBessel_y0(t[0]), t[1]+t[2]) >=real.mant_dig-4); 
     330    } 
     331     
     332    // Values from MS Excel, of doubtful accuracy. 
     333    assert(fabs(-0.060_409_940_421_649 - cylBessel_j0(173.5)) < 0.000_000_000_1); 
     334    assert(fabs(-0.044_733_447_576_5866 - cylBessel_y0(313.25)) < 0.000_000_000_1); 
     335 
    333336} 

  • trunk/docs/bessel.html

    r56 r57  
    66    <!-- Generated by Ddoc from bessel.d --> 
    77<br><br> 
    8 <dl><dt><big>real <u>j0</u>(real <i>x</i>); 
     8<dl><dt><big>real <u>cylBessel_j0</u>(real <i>x</i>); 
    99</big></dt> 
    1010<dd>Bessel function of order zero 
    …… 
    3333 
    3434</dd> 
    35 <dt><big>real <u>y0</u>(real <i>x</i>); 
     35<dt><big>real <u>cylBessel_y0</u>(real <i>x</i>); 
    3636</big></dt> 
    3737<dd>Bessel function of the second kind, order zero 
     38 Also known as the cylindrical Neumann function, order zero. 
    3839<br><br> 
    3940Returns Bessel function of the second kind, of order 
    …… 
    4344 The domain is divided into the intervals [0, 5&gt;, [5,9&gt; and 
    4445 [9, infinity<br><br>. In the first interval a rational approximation 
    45  R(<i>x</i>) is employed to compute <u>y0</u>(<i>x</i>)  = R(<i>x</i>) + 2/pi * log(<i>x</i>) * j0(<i>x</i>). 
     46 R(<i>x</i>) is employed to compute y0(<i>x</i>)  = R(<i>x</i>) + 2/pi * log(<i>x</i>) * j0(<i>x</i>). 
    4647<br><br> 
    4748 
    4849 In the second interval, the approximation is 
    4950     (<i>x</i> - p)(<i>x</i> - q)(<i>x</i> - r)(<i>x</i> - s)P7(<i>x</i>)/Q7(<i>x</i>) 
    50  where p, q, r, s are zeros of <u>y0</u>(<i>x</i>). 
     51 where p, q, r, s are zeros of y0(<i>x</i>). 
    5152<br><br> 
    5253 
    5354 The third interval uses the same approximations to modulus 
    54  and phase as j0(<i>x</i>), whence <u>y0</u>(<i>x</i>) = modulus * sin(phase). 
     55 and phase as j0(<i>x</i>), whence y0(<i>x</i>) = modulus * sin(phase). 
    5556 
    5657 
    5758<b>ACCURACY:</b><br> 
    58 Absolute error, when <u>y0</u>(<i>x</i>) &lt; 1; else relative error: 
     59Absolute error, when y0(<i>x</i>) &lt; 1; else relative error: 
    5960<br><br> 
    6061