1
// Copyright (c) 2013 Christopher Kormanyos
2
// Use, modification and distribution are subject to the
3
// Boost Software License, Version 1.0. (See accompanying file
4
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
6
// This work is based on an earlier work:
7
// "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
8
// in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
10
// This header contains implementation details for estimating the zeros
11
// of cylindrical Bessel and Neumann functions on the positive real axis.
12
// Support is included for both positive as well as negative order.
13
// Various methods are used to estimate the roots. These include
14
// empirical curve fitting and McMahon's asymptotic approximation
15
// for small order, uniform asymptotic expansion for large order,
16
// and iteration and root interlacing for negative order.
18
#ifndef _BESSEL_JY_ZERO_2013_01_18_HPP_
19
#define _BESSEL_JY_ZERO_2013_01_18_HPP_
22
#include <boost/math/constants/constants.hpp>
23
#include <boost/math/special_functions/math_fwd.hpp>
24
#include <boost/math/special_functions/cbrt.hpp>
25
#include <boost/math/special_functions/detail/airy_ai_bi_zero.hpp>
27
namespace boost { namespace math {
33
T equation_nist_10_21_19(const T& v, const T& a)
35
// Get the initial estimate of the m'th root of Jv or Yv.
36
// This subroutine is used for the order m with m > 1.
37
// The order m has been used to create the input parameter a.
39
// This is Eq. 10.21.19 in the NIST Handbook.
40
const T mu = (v * v) * 4U;
41
const T mu_minus_one = mu - T(1);
42
const T eight_a_inv = T(1) / (a * 8U);
43
const T eight_a_inv_squared = eight_a_inv * eight_a_inv;
45
const T term3 = ((mu_minus_one * 4U) * ((mu * 7U) - T(31U) )) / 3U;
46
const T term5 = ((mu_minus_one * 32U) * ((((mu * 83U) - T(982U) ) * mu) + T(3779U) )) / 15U;
47
const T term7 = ((mu_minus_one * 64U) * ((((((mu * 6949U) - T(153855UL)) * mu) + T(1585743UL)) * mu) - T(6277237UL))) / 105U;
49
return a + (((( - term7
50
* eight_a_inv_squared - term5)
51
* eight_a_inv_squared - term3)
52
* eight_a_inv_squared - mu_minus_one)
57
class equation_as_9_3_39_and_its_derivative
60
equation_as_9_3_39_and_its_derivative(const T& zt) : zeta(zt) { }
62
boost::math::tuple<T, T> operator()(const T& z) const
64
BOOST_MATH_STD_USING // ADL of std names, needed for acos, sqrt.
66
// Return the function of zeta that is implicitly defined
67
// in A&S Eq. 9.3.39 as a function of z. The function is
68
// returned along with its derivative with respect to z.
70
const T zsq_minus_one_sqrt = sqrt((z * z) - T(1));
74
- ( acos(T(1) / z) + ((T(2) / 3U) * (zeta * sqrt(zeta)))));
76
const T its_derivative(zsq_minus_one_sqrt / z);
78
return boost::math::tuple<T, T>(the_function, its_derivative);
82
const equation_as_9_3_39_and_its_derivative& operator=(const equation_as_9_3_39_and_its_derivative&);
87
static T equation_as_9_5_26(const T& v, const T& ai_bi_root)
89
BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt.
91
// Obtain the estimate of the m'th zero of Jv or Yv.
92
// The order m has been used to create the input parameter ai_bi_root.
93
// Here, v is larger than about 2.2. The estimate is computed
94
// from Abramowitz and Stegun Eqs. 9.5.22 and 9.5.26, page 371.
96
// The inversion of z as a function of zeta is mentioned in the text
97
// following A&S Eq. 9.5.26. Here, we accomplish the inversion by
98
// performing a Taylor expansion of Eq. 9.3.39 for large z to order 2
99
// and solving the resulting quadratic equation, thereby taking
100
// the positive root of the quadratic.
101
// In other words: (2/3)(-zeta)^(3/2) approx = z + 1/(2z) - pi/2.
102
// This leads to: z^2 - [(2/3)(-zeta)^(3/2) + pi/2]z + 1/2 = 0.
104
// With this initial estimate, Newton-Raphson iteration is used
105
// to refine the value of the estimate of the root of z
106
// as a function of zeta.
108
const T v_pow_third(boost::math::cbrt(v));
109
const T v_pow_minus_two_thirds(T(1) / (v_pow_third * v_pow_third));
111
// Obtain zeta using the order v combined with the m'th root of
112
// an airy function, as shown in A&S Eq. 9.5.22.
113
const T zeta = v_pow_minus_two_thirds * (-ai_bi_root);
115
const T zeta_sqrt = sqrt(zeta);
117
// Set up a quadratic equation based on the Taylor series
118
// expansion mentioned above.
119
const T b = -((((zeta * zeta_sqrt) * 2U) / 3U) + boost::math::constants::half_pi<T>());
121
// Solve the quadratic equation, taking the positive root.
122
const T z_estimate = (-b + sqrt((b * b) - T(2))) / 2U;
124
// Establish the range, the digits, and the iteration limit
125
// for the upcoming root-finding.
126
const T range_zmin = (std::max<T>)(z_estimate - T(1), T(1));
127
const T range_zmax = z_estimate + T(1);
129
const int digits2_of_t = int(float(std::numeric_limits<T>::digits)
130
* ( log(float(std::numeric_limits<T>::radix))
133
const int digits2_for_root = (std::min)(digits2_of_t, std::numeric_limits<double>::digits);
135
boost::uintmax_t iteration_count = boost::uintmax_t(std::numeric_limits<T>::digits10 * 2);
137
// Calculate the root of z as a function of zeta.
138
const T z = boost::math::tools::newton_raphson_iterate(
139
boost::math::detail::bessel_zero::equation_as_9_3_39_and_its_derivative<T>(zeta),
146
static_cast<void>(iteration_count);
148
// Continue with the implementation of A&S Eq. 9.3.39.
149
const T zsq_minus_one = (z * z) - T(1);
150
const T zsq_minus_one_sqrt = sqrt(zsq_minus_one);
152
// This is A&S Eq. 9.3.42.
153
const T b0_term_5_24 = T(5) / ((zsq_minus_one * zsq_minus_one_sqrt) * 24U);
154
const T b0_term_1_8 = T(1) / ( zsq_minus_one_sqrt * 8U);
155
const T b0_term_5_48 = T(5) / ((zeta * zeta) * 48U);
157
const T b0 = -b0_term_5_48 + ((b0_term_5_24 + b0_term_1_8) / zeta_sqrt);
159
// This is the second line of A&S Eq. 9.5.26 for f_k with k = 1.
160
const T f1 = ((z * zeta_sqrt) * b0) / zsq_minus_one_sqrt;
162
// This is A&S Eq. 9.5.22 expanded to k = 1 (i.e., one term in the series).
163
return (v * z) + (f1 / v);
166
namespace cyl_bessel_j_zero_detail
169
T equation_nist_10_21_40_a(const T& v)
171
const T v_pow_third(boost::math::cbrt(v));
172
const T v_pow_minus_two_thirds(T(1) / (v_pow_third * v_pow_third));
174
return v * ((((( + T(0.043)
175
* v_pow_minus_two_thirds - T(0.0908))
176
* v_pow_minus_two_thirds - T(0.00397))
177
* v_pow_minus_two_thirds + T(1.033150))
178
* v_pow_minus_two_thirds + T(1.8557571))
179
* v_pow_minus_two_thirds + T(1));
182
template<class T, class Policy>
183
class function_object_jv
186
function_object_jv(const T& v,
187
const Policy& pol) : my_v(v),
190
T operator()(const T& x) const
192
return boost::math::cyl_bessel_j(my_v, x, my_pol);
197
const Policy& my_pol;
198
const function_object_jv& operator=(const function_object_jv&);
201
template<class T, class Policy>
202
class function_object_jv_and_jv_prime
205
function_object_jv_and_jv_prime(const T& v,
206
const bool order_is_zero,
207
const Policy& pol) : my_v(v),
208
my_order_is_zero(order_is_zero),
211
boost::math::tuple<T, T> operator()(const T& x) const
213
// Obtain Jv(x) and Jv'(x).
214
// Chris's original code called the Bessel function implementation layer direct,
215
// but that circumvented optimizations for integer-orders. Call the documented
216
// top level functions instead, and let them sort out which implementation to use.
222
j_v = boost::math::cyl_bessel_j(0, x, my_pol);
223
j_v_prime = -boost::math::cyl_bessel_j(1, x, my_pol);
227
j_v = boost::math::cyl_bessel_j( my_v, x, my_pol);
228
const T j_v_m1 (boost::math::cyl_bessel_j(T(my_v - 1), x, my_pol));
229
j_v_prime = j_v_m1 - ((my_v * j_v) / x);
232
// Return a tuple containing both Jv(x) and Jv'(x).
233
return boost::math::make_tuple(j_v, j_v_prime);
238
const bool my_order_is_zero;
239
const Policy& my_pol;
240
const function_object_jv_and_jv_prime& operator=(const function_object_jv_and_jv_prime&);
243
template<class T> bool my_bisection_unreachable_tolerance(const T&, const T&) { return false; }
245
template<class T, class Policy>
246
T initial_guess(const T& v, const int m, const Policy& pol)
248
BOOST_MATH_STD_USING // ADL of std names, needed for floor.
250
// Compute an estimate of the m'th root of cyl_bessel_j.
254
// There is special handling for negative order.
257
if((m == 1) && (v > -0.5F))
259
// For small, negative v, use the results of empirical curve fitting.
260
// Mathematica(R) session for the coefficients:
261
// Table[{n, BesselJZero[n, 1]}, {n, -(1/2), 0, 1/10}]
263
// Fit[%, {n^0, n^1, n^2, n^3, n^4, n^5, n^6}, n]
264
guess = ((((( - T(0.2321156900729)
265
* v - T(0.1493247777488))
266
* v - T(0.15205419167239))
267
* v + T(0.07814930561249))
268
* v - T(0.17757573537688))
269
* v + T(1.542805677045663))
270
* v + T(2.40482555769577277);
275
// Create the positive order and extract its positive floor integer part.
277
const T vv_floor(floor(vv));
279
// The to-be-found root is bracketed by the roots of the
280
// Bessel function whose reflected, positive integer order
281
// is less than, but nearest to vv.
283
T root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(vv_floor, m, pol);
288
// The estimate of the first root for negative order is found using
289
// an adaptive range-searching algorithm.
290
root_lo = T(root_hi - 0.1F);
292
const bool hi_end_of_bracket_is_negative = (boost::math::cyl_bessel_j(v, root_hi, pol) < 0);
294
while((root_lo > boost::math::tools::epsilon<T>()))
296
const bool lo_end_of_bracket_is_negative = (boost::math::cyl_bessel_j(v, root_lo, pol) < 0);
298
if(hi_end_of_bracket_is_negative != lo_end_of_bracket_is_negative)
305
// Decrease the lower end of the bracket using an adaptive algorithm.
318
root_lo = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(vv_floor, m - 1, pol);
321
// Perform several steps of bisection iteration to refine the guess.
322
boost::uintmax_t number_of_iterations(12U);
324
// Do the bisection iteration.
325
const boost::math::tuple<T, T> guess_pair =
326
boost::math::tools::bisect(
327
boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::function_object_jv<T, Policy>(v, pol),
330
boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::my_bisection_unreachable_tolerance<T>,
331
number_of_iterations);
333
return (boost::math::get<0>(guess_pair) + boost::math::get<1>(guess_pair)) / 2U;
338
// Get the initial estimate of the first root.
342
// For small v, use the results of empirical curve fitting.
343
// Mathematica(R) session for the coefficients:
344
// Table[{n, BesselJZero[n, 1]}, {n, 0, 22/10, 1/10}]
346
// Fit[%, {n^0, n^1, n^2, n^3, n^4, n^5, n^6}, n]
347
guess = ((((( - T(0.0008342379046010)
348
* v + T(0.007590035637410))
349
* v - T(0.030640914772013))
350
* v + T(0.078232088020106))
351
* v - T(0.169668712590620))
352
* v + T(1.542187960073750))
353
* v + T(2.4048359915254634);
357
// For larger v, use the first line of Eqs. 10.21.40 in the NIST Handbook.
358
guess = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::equation_nist_10_21_40_a(v);
365
// Use Eq. 10.21.19 in the NIST Handbook.
366
const T a(((v + T(m * 2U)) - T(0.5)) * boost::math::constants::half_pi<T>());
368
guess = boost::math::detail::bessel_zero::equation_nist_10_21_19(v, a);
372
// Get an estimate of the m'th root of airy_ai.
373
const T airy_ai_root(boost::math::detail::airy_zero::airy_ai_zero_detail::initial_guess<T>(m));
375
// Use Eq. 9.5.26 in the A&S Handbook.
376
guess = boost::math::detail::bessel_zero::equation_as_9_5_26(v, airy_ai_root);
382
} // namespace cyl_bessel_j_zero_detail
384
namespace cyl_neumann_zero_detail
387
T equation_nist_10_21_40_b(const T& v)
389
const T v_pow_third(boost::math::cbrt(v));
390
const T v_pow_minus_two_thirds(T(1) / (v_pow_third * v_pow_third));
392
return v * ((((( - T(0.001)
393
* v_pow_minus_two_thirds - T(0.0060))
394
* v_pow_minus_two_thirds + T(0.01198))
395
* v_pow_minus_two_thirds + T(0.260351))
396
* v_pow_minus_two_thirds + T(0.9315768))
397
* v_pow_minus_two_thirds + T(1));
400
template<class T, class Policy>
401
class function_object_yv
404
function_object_yv(const T& v,
405
const Policy& pol) : my_v(v),
408
T operator()(const T& x) const
410
return boost::math::cyl_neumann(my_v, x, my_pol);
415
const Policy& my_pol;
416
const function_object_yv& operator=(const function_object_yv&);
419
template<class T, class Policy>
420
class function_object_yv_and_yv_prime
423
function_object_yv_and_yv_prime(const T& v,
424
const Policy& pol) : my_v(v),
427
boost::math::tuple<T, T> operator()(const T& x) const
429
const T half_epsilon(boost::math::tools::epsilon<T>() / 2U);
431
const bool order_is_zero = ((my_v > -half_epsilon) && (my_v < +half_epsilon));
433
// Obtain Yv(x) and Yv'(x).
434
// Chris's original code called the Bessel function implementation layer direct,
435
// but that circumvented optimizations for integer-orders. Call the documented
436
// top level functions instead, and let them sort out which implementation to use.
442
y_v = boost::math::cyl_neumann(0, x, my_pol);
443
y_v_prime = -boost::math::cyl_neumann(1, x, my_pol);
447
y_v = boost::math::cyl_neumann( my_v, x, my_pol);
448
const T y_v_m1 (boost::math::cyl_neumann(T(my_v - 1), x, my_pol));
449
y_v_prime = y_v_m1 - ((my_v * y_v) / x);
452
// Return a tuple containing both Yv(x) and Yv'(x).
453
return boost::math::make_tuple(y_v, y_v_prime);
458
const Policy& my_pol;
459
const function_object_yv_and_yv_prime& operator=(const function_object_yv_and_yv_prime&);
462
template<class T> bool my_bisection_unreachable_tolerance(const T&, const T&) { return false; }
464
template<class T, class Policy>
465
T initial_guess(const T& v, const int m, const Policy& pol)
467
BOOST_MATH_STD_USING // ADL of std names, needed for floor.
469
// Compute an estimate of the m'th root of cyl_neumann.
473
// There is special handling for negative order.
476
// Create the positive order and extract its positive floor and ceiling integer parts.
478
const T vv_floor(floor(vv));
480
// The to-be-found root is bracketed by the roots of the
481
// Bessel function whose reflected, positive integer order
482
// is less than, but nearest to vv.
484
// The special case of negative, half-integer order uses
485
// the relation between Yv and spherical Bessel functions
486
// in order to obtain the bracket for the root.
487
// In these special cases, cyl_neumann(-n/2, x) = sph_bessel_j(+n/2, x)
495
// The estimate of the first root for negative order is found using
496
// an adaptive range-searching algorithm.
497
// Take special precautions for the discontinuity at negative,
498
// half-integer orders and use different brackets above and below these.
499
if(T(vv - vv_floor) < 0.5F)
501
root_hi = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m, pol);
505
root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m, pol);
508
root_lo = T(root_hi - 0.1F);
510
const bool hi_end_of_bracket_is_negative = (boost::math::cyl_neumann(v, root_hi, pol) < 0);
512
while((root_lo > boost::math::tools::epsilon<T>()))
514
const bool lo_end_of_bracket_is_negative = (boost::math::cyl_neumann(v, root_lo, pol) < 0);
516
if(hi_end_of_bracket_is_negative != lo_end_of_bracket_is_negative)
523
// Decrease the lower end of the bracket using an adaptive algorithm.
536
if(T(vv - vv_floor) < 0.5F)
538
root_lo = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m - 1, pol);
539
root_hi = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m, pol);
545
root_lo = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m - 1, pol);
546
root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m, pol);
552
// Perform several steps of bisection iteration to refine the guess.
553
boost::uintmax_t number_of_iterations(12U);
555
// Do the bisection iteration.
556
const boost::math::tuple<T, T> guess_pair =
557
boost::math::tools::bisect(
558
boost::math::detail::bessel_zero::cyl_neumann_zero_detail::function_object_yv<T, Policy>(v, pol),
561
boost::math::detail::bessel_zero::cyl_neumann_zero_detail::my_bisection_unreachable_tolerance<T>,
562
number_of_iterations);
564
return (boost::math::get<0>(guess_pair) + boost::math::get<1>(guess_pair)) / 2U;
569
// Get the initial estimate of the first root.
573
// For small v, use the results of empirical curve fitting.
574
// Mathematica(R) session for the coefficients:
575
// Table[{n, BesselYZero[n, 1]}, {n, 0, 22/10, 1/10}]
577
// Fit[%, {n^0, n^1, n^2, n^3, n^4, n^5, n^6}, n]
578
guess = ((((( - T(0.0025095909235652)
579
* v + T(0.021291887049053))
580
* v - T(0.076487785486526))
581
* v + T(0.159110268115362))
582
* v - T(0.241681668765196))
583
* v + T(1.4437846310885244))
584
* v + T(0.89362115190200490);
588
// For larger v, use the second line of Eqs. 10.21.40 in the NIST Handbook.
589
guess = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::equation_nist_10_21_40_b(v);
596
// Use Eq. 10.21.19 in the NIST Handbook.
597
const T a(((v + T(m * 2U)) - T(1.5)) * boost::math::constants::half_pi<T>());
599
guess = boost::math::detail::bessel_zero::equation_nist_10_21_19(v, a);
603
// Get an estimate of the m'th root of airy_bi.
604
const T airy_bi_root(boost::math::detail::airy_zero::airy_bi_zero_detail::initial_guess<T>(m));
606
// Use Eq. 9.5.26 in the A&S Handbook.
607
guess = boost::math::detail::bessel_zero::equation_as_9_5_26(v, airy_bi_root);
613
} // namespace cyl_neumann_zero_detail
614
} // namespace bessel_zero
615
} } } // namespace boost::math::detail
617
#endif // _BESSEL_JY_ZERO_2013_01_18_HPP_
added
removed
boost/includes/boost/math/special_functions/detail/bessel_jy_zero.hpp
