ceph/src/boost/libs/math/doc/sf/bessel_jy.qbk

   1
   2 [section:bessel_first Bessel Functions of the First and Second Kinds]
   3
   4 [h4 Synopsis]
   5
   6 `#include <boost/math/special_functions/bessel.hpp>`
   7
   8    template <class T1, class T2>
   9    ``__sf_result`` cyl_bessel_j(T1 v, T2 x);
  10
  11    template <class T1, class T2, class ``__Policy``>
  12    ``__sf_result`` cyl_bessel_j(T1 v, T2 x, const ``__Policy``&);
  13
  14    template <class T1, class T2>
  15    ``__sf_result`` cyl_neumann(T1 v, T2 x);
  16
  17    template <class T1, class T2, class ``__Policy``>
  18    ``__sf_result`` cyl_neumann(T1 v, T2 x, const ``__Policy``&);
  19
  20
  21 [h4 Description]
  22
  23 The functions __cyl_bessel_j and __cyl_neumann return the result of the
  24 Bessel functions of the first and second kinds respectively:
  25
  26 cyl_bessel_j(v, x) = J[sub v](x)
  27
  28 cyl_neumann(v, x) = Y[sub v](x) = N[sub v](x)
  29
  30 where:
  31
  32 [equation bessel2]
  33
  34 [equation bessel3]
  35
  36 The return type of these functions is computed using the __arg_promotion_rules
  37 when T1 and T2 are different types.  The functions are also optimised for the
  38 relatively common case that T1 is an integer.
  39
  40 [optional_policy]
  41
  42 The functions return the result of __domain_error whenever the result is
  43 undefined or complex.  For __cyl_bessel_j this occurs when `x < 0` and v is not
  44 an integer, or when `x == 0` and `v != 0`.  For __cyl_neumann this occurs
  45 when `x <= 0`.
  46
  47 The following graph illustrates the cyclic nature of J[sub v]:
  48
  49 [graph cyl_bessel_j]
  50
  51 The following graph shows the behaviour of Y[sub v]: this is also
  52 cyclic for large /x/, but tends to -[infin][space] for small /x/:
  53
  54 [graph cyl_neumann]
  55
  56 [h4 Testing]
  57
  58 There are two sets of test values: spot values calculated using
  59 [@http://functions.wolfram.com functions.wolfram.com],
  60 and a much larger set of tests computed using
  61 a simplified version of this implementation
  62 (with all the special case handling removed).
  63
  64 [h4 Accuracy]
  65
  66 The following tables show how the accuracy of these functions
  67 varies on various platforms, along with comparisons to other
  68 libraries.  Note that the cyclic nature of these
  69 functions means that they have an infinite number of irrational
  70 roots: in general these functions have arbitrarily large /relative/
  71 errors when the arguments are sufficiently close to a root.  Of
  72 course the absolute error in such cases is always small.
  73 Note that only results for the widest floating-point type on the
  74 system are given as narrower types have __zero_error.  All values
  75 are relative errors in units of epsilon.  Most of the gross errors
  76 exhibited by other libraries occur for very large arguments - you will
  77 need to drill down into the actual program output if you need more
  78 information on this.
  79
  80 [table_cyl_bessel_j_integer_orders_]
  81
  82 [table_cyl_bessel_j]
  83
  84 [table_cyl_neumann_integer_orders_]
  85
  86 [table_cyl_neumann]
  87
  88 Note that for large /x/ these functions are largely dependent on
  89 the accuracy of the `std::sin` and `std::cos` functions.
  90
  91 Comparison to GSL and __cephes is interesting: both __cephes and this library optimise
  92 the integer order case - leading to identical results - simply using the general
  93 case is for the most part slightly more accurate though, as noted by the
  94 better accuracy of GSL in the integer argument cases.  This implementation tends to
  95 perform much better when the arguments become large, __cephes in particular produces
  96 some remarkably inaccurate results with some of the test data (no significant figures
  97 correct), and even GSL performs badly with some inputs to J[sub v].  Note that
  98 by way of double-checking these results, the worst performing __cephes and GSL cases
  99 were recomputed using [@http://functions.wolfram.com functions.wolfram.com],
 100 and the result checked against our test data: no errors in the test data were found.
 101
 102 [h4 Implementation]
 103
 104 The implementation is mostly about filtering off various special cases:
 105
 106 When /x/ is negative, then the order /v/ must be an integer or the
 107 result is a domain error.  If the order is an integer then the function
 108 is odd for odd orders and even for even orders, so we reflect to /x > 0/.
 109
 110 When the order /v/ is negative then the reflection formulae can be used to
 111 move to /v > 0/:
 112
 113 [equation bessel9]
 114
 115 [equation bessel10]
 116
 117 Note that if the order is an integer, then these formulae reduce to:
 118
 119 J[sub -n] = (-1)[super n]J[sub n]
 120
 121 Y[sub -n] = (-1)[super n]Y[sub n]
 122
 123 However, in general, a negative order implies that we will need to compute
 124 both J and Y.
 125
 126 When /x/ is large compared to the order /v/ then the asymptotic expansions
 127 for large /x/ in M. Abramowitz and I.A. Stegun,
 128 ['Handbook of Mathematical Functions] 9.2.19 are used
 129 (these were found to be more reliable
 130 than those in A&S 9.2.5).
 131
 132 When the order /v/ is an integer the method first relates the result
 133 to J[sub 0], J[sub 1], Y[sub 0][space] and Y[sub 1][space] using either
 134 forwards or backwards recurrence (Miller's algorithm) depending upon which is stable.
 135 The values for J[sub 0], J[sub 1], Y[sub 0][space] and Y[sub 1][space] are
 136 calculated using the rational minimax approximations on
 137 root-bracketing intervals for small ['|x|] and Hankel asymptotic
 138 expansion for large ['|x|]. The coefficients are from:
 139
 140 W.J. Cody, ['ALGORITHM 715: SPECFUN - A Portable FORTRAN Package of
 141 Special Function Routines and Test Drivers], ACM Transactions on Mathematical
 142 Software, vol 19, 22 (1993).
 143
 144 and
 145
 146 J.F. Hart et al, ['Computer Approximations], John Wiley & Sons, New York, 1968.
 147
 148 These approximations are accurate to around 19 decimal digits: therefore
 149 these methods are not used when type T has more than 64 binary digits.
 150
 151 When /x/ is smaller than machine epsilon then the following approximations for
 152 Y[sub 0](x), Y[sub 1](x), Y[sub 2](x) and Y[sub n](x) can be used
 153 (see: [@http://functions.wolfram.com/03.03.06.0037.01 http://functions.wolfram.com/03.03.06.0037.01],
 154 [@http://functions.wolfram.com/03.03.06.0038.01 http://functions.wolfram.com/03.03.06.0038.01],
 155 [@http://functions.wolfram.com/03.03.06.0039.01 http://functions.wolfram.com/03.03.06.0039.01]
 156 and [@http://functions.wolfram.com/03.03.06.0040.01 http://functions.wolfram.com/03.03.06.0040.01]):
 157
 158 [equation bessel_y0_small_z]
 159
 160 [equation bessel_y1_small_z]
 161
 162 [equation bessel_y2_small_z]
 163
 164 [equation bessel_yn_small_z]
 165
 166 When /x/ is small compared to /v/ and /v/ is not an integer, then the following
 167 series approximation can be used for Y[sub v](x), this is also an area where other
 168 approximations are often too slow to converge to be used
 169 (see [@http://functions.wolfram.com/03.03.06.0034.01 http://functions.wolfram.com/03.03.06.0034.01]):
 170
 171 [equation bessel_yv_small_z]
 172
 173 When /x/ is small compared to /v/, J[sub v]x[space] is best computed directly from the series:
 174
 175 [equation bessel2]
 176
 177 In the general case we compute J[sub v][space] and
 178 Y[sub v][space] simultaneously.
 179
 180 To get the initial values, let
 181 [mu][space] = [nu] - floor([nu] + 1/2), then [mu][space] is the fractional part
 182 of [nu][space] such that
 183 |[mu]| <= 1/2 (we need this for convergence later). The idea is to
 184 calculate J[sub [mu]](x), J[sub [mu]+1](x), Y[sub [mu]](x), Y[sub [mu]+1](x)
 185 and use them to obtain J[sub [nu]](x), Y[sub [nu]](x).
 186
 187 The algorithm is called Steed's method, which needs two
 188 continued fractions as well as the Wronskian:
 189
 190 [equation bessel8]
 191
 192 [equation bessel11]
 193
 194 [equation bessel12]
 195
 196 See: F.S. Acton, ['Numerical Methods that Work],
 197     The Mathematical Association of America, Washington, 1997.
 198
 199 The continued fractions are computed using the modified Lentz's method
 200 (W.J. Lentz, ['Generating Bessel functions in Mie scattering calculations
 201 using continued fractions], Applied Optics, vol 15, 668 (1976)).
 202 Their convergence rates depend on ['x], therefore we need
 203 different strategies for large ['x] and small ['x].
 204
 205 ['x > v], CF1 needs O(['x]) iterations to converge, CF2 converges rapidly
 206
 207 ['x <= v], CF1 converges rapidly, CF2 fails to converge when ['x] [^->] 0
 208
 209 When ['x] is large (['x] > 2), both continued fractions converge (CF1
 210 may be slow for really large ['x]). J[sub [mu]], J[sub [mu]+1],
 211 Y[sub [mu]], Y[sub [mu]+1] can be calculated by
 212
 213 [equation bessel13]
 214
 215 where
 216
 217 [equation bessel14]
 218
 219 J[sub [nu]] and Y[sub [mu]] are then calculated using backward
 220 (Miller's algorithm) and forward recurrence respectively.
 221
 222 When ['x] is small (['x] <= 2), CF2 convergence may fail (but CF1
 223 works very well). The solution here is Temme's series:
 224
 225 [equation bessel15]
 226
 227 where
 228
 229 [equation bessel16]
 230
 231 g[sub k][space] and h[sub k][space]
 232 are also computed by recursions (involving gamma functions), but the
 233 formulas are a little complicated, readers are refered to
 234 N.M. Temme, ['On the numerical evaluation of the ordinary Bessel function
 235 of the second kind], Journal of Computational Physics, vol 21, 343 (1976).
 236 Note Temme's series converge only for |[mu]| <= 1/2.
 237
 238 As the previous case, Y[sub [nu]][space] is calculated from the forward
 239 recurrence, so is Y[sub [nu]+1]. With these two
 240 values and f[sub [nu]], the Wronskian yields J[sub [nu]](x) directly
 241 without backward recurrence.
 242
 243 [endsect]
 244
 245 [section:bessel_root Finding Zeros of Bessel Functions of the First and Second Kinds]
 246
 247 [h4 Synopsis]
 248
 249 `#include <boost/math/special_functions/bessel.hpp>`
 250
 251 Functions for obtaining both a single zero or root of the Bessel function,
 252 and placing multiple zeros into a container like `std::vector`
 253 by providing an output iterator.
 254
 255 The signature of the single value functions are:
 256
 257   template <class T>
 258   T cyl_bessel_j_zero(
 259            T v,            // Floating-point value for Jv.
 260            int m);         // 1-based index of zero.
 261
 262   template <class T>
 263   T cyl_neumann_zero(
 264            T v,            // Floating-point value for Jv.
 265            int m);         // 1-based index of zero.
 266
 267 and for multiple zeros:
 268
 269  template <class T, class OutputIterator>
 270  OutputIterator cyl_bessel_j_zero(
 271                       T v,                       // Floating-point value for Jv.
 272                       int start_index,           // 1-based index of first zero.
 273                       unsigned number_of_zeros,  // How many zeros to generate.
 274                       OutputIterator out_it);    // Destination for zeros.
 275
 276  template <class T, class OutputIterator>
 277  OutputIterator cyl_neumann_zero(
 278                       T v,                       // Floating-point value for Jv.
 279                       int start_index,           // 1-based index of zero.
 280                       unsigned number_of_zeros,  // How many zeros to generate
 281                       OutputIterator out_it);    // Destination for zeros.
 282
 283 There are also versions which allow control of the __policy_section for error handling and precision.
 284
 285   template <class T>
 286   T cyl_bessel_j_zero(
 287            T v,            // Floating-point value for Jv.
 288            int m,          // 1-based index of zero.
 289            const Policy&); // Policy to use.
 290
 291   template <class T>
 292   T cyl_neumann_zero(
 293            T v,            // Floating-point value for Jv.
 294            int m,          // 1-based index of zero.
 295            const Policy&); // Policy to use.
 296
 297
 298  template <class T, class OutputIterator>
 299  OutputIterator cyl_bessel_j_zero(
 300                       T v,                       // Floating-point value for Jv.
 301                       int start_index,           // 1-based index of first zero.
 302                       unsigned number_of_zeros,  // How many zeros to generate.
 303                       OutputIterator out_it,     // Destination for zeros.
 304                       const Policy& pol);        // Policy to use.
 305
 306  template <class T, class OutputIterator>
 307  OutputIterator cyl_neumann_zero(
 308                       T v,                       // Floating-point value for Jv.
 309                       int start_index,           // 1-based index of zero.
 310                       unsigned number_of_zeros,  // How many zeros to generate.
 311                       OutputIterator out_it,     // Destination for zeros.
 312                       const Policy& pol);        // Policy to use.
 313
 314 [h4 Description]
 315
 316 Every real order [nu] cylindrical Bessel and Neumann functions have an infinite
 317 number of zeros on the positive real axis. The real zeros on the positive real
 318 axis can be found by solving for the roots of
 319
 320 [emquad] ['J[sub [nu]](j[sub [nu], m]) = 0]
 321
 322 [emquad] ['Y[sub [nu]](y[sub [nu], m]) = 0]
 323
 324 Here, ['j[sub [nu], m]] represents the ['m[super th]]
 325 root of the cylindrical Bessel function of order ['[nu]],
 326 and ['y[sub [nu], m]] represents the ['m[super th]]
 327 root of the cylindrical Neumann function of order ['[nu]].
 328
 329 The zeros or roots (values of `x` where the function crosses the horizontal `y = 0` axis)
 330 of the Bessel and Neumann functions are computed by two functions,
 331 `cyl_bessel_j_zero` and `cyl_neumann_zero`.
 332
 333 In each case the index or rank of the zero
 334 returned is 1-based, which is to say:
 335
 336    cyl_bessel_j_zero(v, 1);
 337
 338 returns the first zero of Bessel J.
 339
 340 Passing an `start_index <= 0` results in a `std::domain_error` being raised.
 341
 342 For certain parameters, however, the zero'th root is defined and
 343 it has a value of zero. For example, the zero'th root
 344 of `J[sub v](x)` is defined and it has a value of zero for all
 345 values of `v > 0` and for negative integer values of `v = -n`.
 346 Similar cases are described in the implementation details below.
 347
 348 The order `v` of `J` can be positive, negative and zero for the `cyl_bessel_j`
 349 and `cyl_neumann` functions, but not infinite nor NaN.
 350
 351 [graph bessel_j_zeros]
 352
 353 [graph neumann_y_zeros]
 354
 355 [h4 Examples of finding Bessel and Neumann zeros]
 356
 357 [import ../../example/bessel_zeros_example_1.cpp]
 358
 359 [bessel_zeros_example_1]
 360 [bessel_zeros_example_2]
 361
 362 [import ../../example/bessel_zeros_interator_example.cpp]
 363
 364 [bessel_zeros_iterator_example_1]
 365 [bessel_zeros_iterator_example_2]
 366
 367 [import ../../example/neumann_zeros_example_1.cpp]
 368
 369 [neumann_zeros_example_1]
 370 [neumann_zeros_example_2]
 371
 372 [import ../../example/bessel_errors_example.cpp]
 373
 374 [bessel_errors_example_1]
 375 [bessel_errors_example_2]
 376
 377 The full code (and output) for these examples is at
 378 [@../../example/bessel_zeros_example_1.cpp Bessel zeros],
 379 [@../../example/bessel_zeros_interator_example.cpp Bessel zeros iterator],
 380 [@../../example/neumann_zeros_example_1.cpp Neumann zeros],
 381 [@../../example/bessel_errors_example.cpp  Bessel error messages].
 382
 383 [h3 Implementation]
 384
 385 Various methods are used to compute initial estimates
 386 for ['j[sub [nu], m]] and ['y[sub [nu], m]] ; these are described in detail below.
 387
 388 After finding the initial estimate of a given root,
 389 its precision is subsequently refined to the desired level
 390 using Newton-Raphson iteration from Boost.Math's __root_finding_with_derivatives
 391 utilities combined with the functions __cyl_bessel_j  and __cyl_neumann.
 392
 393 Newton iteration requires both ['J[sub [nu]](x)] or ['Y[sub [nu]](x)]
 394 as well as its derivative. The derivatives of ['J[sub [nu]](x)] and ['Y[sub [nu]](x)]
 395 with respect to  ['x] are given by __Abramowitz_Stegun.
 396 In particular,
 397
 398 [emquad] ['d/[sub dx] ['J[sub [nu]](x)] = ['J[sub [nu]-1](x)] - [nu] J[sub [nu]](x)] / x
 399
 400 [emquad] ['d/[sub dx] ['Y[sub [nu]](x)] = ['Y[sub [nu]-1](x)] - [nu] Y[sub [nu]](x)] / x
 401
 402 Enumeration of the rank of a root (in other words the index of a root)
 403 begins with one and counts up, in other words
 404 ['m,=1,2,3,[ellipsis]] The value of the first root is always greater than zero.
 405
 406 For certain special parameters, cylindrical Bessel functions
 407 and cylindrical Neumann functions have a root at the origin. For example,
 408 ['J[sub [nu]](x)] has a root at the origin for every positive order
 409 ['[nu] > 0], and for every negative integer order
 410 ['[nu] = -n] with ['n [isin] [negative] [super +]] and ['n [ne] 0].
 411
 412 In addition, ['Y[sub [nu]](x)] has a root at the origin
 413 for every negative half-integer order ['[nu] = -n/2], with ['n [isin] [negative] [super +]] and
 414 and ['n [ne] 0].
 415
 416 For these special parameter values, the origin with
 417 a value of ['x = 0] is provided as the ['0[super th]]
 418 root generated by `cyl_bessel_j_zero()`
 419 and `cyl_neumann_zero()`.
 420
 421 When calculating initial estimates for the roots
 422 of Bessel functions, a distinction is made between
 423 positive order and negative order, and different
 424 methods are used for these. In addition, different algorithms
 425 are used for the first root ['m = 1] and
 426 for subsequent roots with higher rank ['m [ge] 2].
 427 Furthermore, estimates of the roots for Bessel functions
 428 with order above and below a cutoff at ['[nu] = 2.2]
 429 are calculated with different methods.
 430
 431 Calculations of the estimates of ['j[sub [nu],1]] and ['y[sub [nu],1]]
 432 with ['0 [le] [nu] < 2.2] use empirically tabulated values.
 433 The coefficients for these have been generated by a
 434 computer algebra system.
 435
 436 Calculations of the estimates of ['j[sub [nu],1]] and ['y[sub [nu],1]]
 437 with ['[nu][ge] 2.2] use Eqs.9.5.14 and 9.5.15 in __Abramowitz_Stegun.
 438
 439 In particular,
 440
 441 [emquad] ['j[sub [nu],1] [cong] [nu] + 1.85575 [nu][super [frac13]] + 1.033150 [nu][super -[frac13]] - 0.00397 [nu][super -1] - 0.0908 [nu][super -5/3] + 0.043 [nu][super -7/3] + [ellipsis]]
 442
 443 and
 444
 445 [emquad] ['y[sub [nu],1] [cong] [nu] + 0.93157 [nu][super [frac13]] + 0.26035 [nu][super -[frac13]] + 0.01198 [nu][super -1] - 0.0060 [nu][super -5/3] - 0.001 [nu][super -7/3] + [ellipsis]]
 446
 447 Calculations of the estimates of ['j[sub [nu], m]]  and  ['y[sub [nu], m]]
 448 with rank ['m > 2] and ['0 [le] [nu] < 2.2]  use
 449 McMahon's approximation, as described in M. Abramowitz and I. A. Stegan, Section 9.5 and 9.5.12.
 450 In particular,
 451
 452 [emquad] ['j[sub [nu],m], y[sub [nu],m] [cong] [beta] - ([mu]-1) / 8[beta]]
 453
 454 [emquad] [emquad] [emquad] ['- 4([mu]-1)(7[mu] - 31) / 3(8[beta])[super 3]]
 455
 456 [emquad] [emquad] [emquad] ['-32([mu]-1)(83[mu][sup2] - 982[mu] + 3779) / 15(8[beta])[super 5]]
 457
 458 [emquad] [emquad] [emquad] ['-64([mu]-1)(6949[mu][super 3] - 153855[mu][sup2] + 1585743[mu]- 6277237) / 105(8a)[super 7]]
 459
 460 [emquad] [emquad] [emquad] ['- [ellipsis]]       [emquad] [emquad]                                              (5)
 461
 462 where ['[mu] = 4[nu][super 2]] and ['[beta] = (m + [frac12][nu] - [frac14])[pi]]
 463 for ['j[sub [nu],m]] and
 464 ['[beta] = (m + [frac12][nu] -[frac34])[pi] for ['y[sub [nu],m]]].
 465
 466 Calculations of the estimates of ['j[sub [nu], m]]  and  ['y[sub [nu], m]]
 467 with ['[nu] [ge] 2.2] use
 468 one term in the asymptotic expansion given in
 469 Eq.9.5.22 and top line of Eq.9.5.26 combined with Eq. 9.3.39,
 470 all in __Abramowitz_Stegun explicit and easy-to-understand treatment
 471 for asymptotic expansion of zeros.
 472 The latter two equations are expressed for argument ['(x)] greater than one.
 473 (Olver also gives the series form of the equations in
 474 [@http://dlmf.nist.gov/10.21#vi [sect]10.21(vi) McMahon's Asymptotic Expansions for Large Zeros] - using slightly different variable names).
 475
 476 In summary,
 477
 478 [emquad] ['j[sub [nu], m] [sim] [nu]x(-[zeta]) + f[sub 1](-[zeta]/[nu])]
 479
 480 where ['-[zeta] = [nu][super -2/3]a[sub m]] and ['a[sub m]] is
 481 the absolute value of the ['m[super th]] root of ['Ai(x)] on the negative real axis.
 482
 483 Here ['x = x(-[zeta])] is the inverse of the function
 484
 485 [emquad] ['[frac23](-[zeta])[super 3/2] = [radic](x[sup2] - 1) - cos[supminus][sup1](1/x)]  [emquad] [emquad]   (7)
 486
 487 Furthermore,
 488
 489 [emquad] ['f[sub 1](-[zeta]) = [frac12]x(-[zeta]) {h(-[zeta])}[sup2] [sdot] b[sub 0](-[zeta])]
 490
 491 where
 492
 493 [emquad] ['h(-[zeta]) = {4(-[zeta]) / (x[sup2] - 1)}[super 4]]
 494
 495 and
 496
 497 [emquad] ['b[sub 0](-[zeta]) = -5/(48[zeta][sup2]) + 1/(-[zeta])[super [frac12]] [sdot] { 5/(24(x[super 2]-1)[super 3/2]) + 1/(8(x[super 2]-1)[super [frac12])]}]
 498
 499 When solving for ['x(-[zeta])] in Eq. 7 above,
 500 the right-hand-side is expanded to order 2 in
 501 a Taylor series for large ['x]. This results in
 502
 503 [emquad] ['[frac23](-[zeta])[super 3/2] [approx] x + 1/2x - [pi]/2]
 504
 505 The positive root of the resulting quadratic equation
 506 is used to find an initial estimate ['x(-[zeta])].
 507 This initial estimate is subsequently refined with
 508 several steps of Newton-Raphson iteration
 509 in Eq. 7.
 510
 511 Estimates of the roots of cylindrical Bessel functions
 512 of negative order on the positive real axis are found
 513 using interlacing relations. For example, the ['m[super th]]
 514 root of the cylindrical Bessel function ['j[sub -[nu],m]]
 515 is bracketed by the ['m[super th]] root and the
 516 ['(m+1)[super th]] root of the Bessel function of
 517 corresponding positive integer order. In other words,
 518
 519 [emquad]  ['j[sub n[nu],m]] <  ['j[sub -[nu],m]] <  ['j[sub n[nu],m+1]]
 520
 521 where ['m > 1] and ['n[sub [nu]]] represents the integral
 522 floor of the absolute value of ['|-[nu]|].
 523
 524 Similar bracketing relations are used to find estimates
 525 of the roots of Neumann functions of negative order,
 526 whereby a discontinuity at every negative half-integer
 527 order needs to be handled.
 528
 529 Bracketing relations do not hold for the first root
 530 of cylindrical Bessel functions and cylindrical Neumann
 531 functions with negative order. Therefore, iterative algorithms
 532 combined with root-finding via bisection are used
 533 to localize ['j[sub -[nu],1]] and ['y[sub -[nu],1]].
 534
 535 [h3 Testing]
 536
 537 The precision of evaluation of zeros was tested at 50 decimal digits using `cpp_dec_float_50`
 538 and found identical with spot values computed by __WolframAlpha.
 539
 540 [endsect]  [/section:bessel Finding Zeros of Bessel Functions of the First and Second Kinds]
 541
 542 [/
 543   Copyright 2006, 2013 John Maddock, Paul A. Bristow, Xiaogang Zhang and Christopher Kormanyos.
 544
 545   Distributed under the Boost Software License, Version 1.0.
 546   (See accompanying file LICENSE_1_0.txt or copy at
 547   http://www.boost.org/LICENSE_1_0.txt).
 548 ]