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) |
5 | // |
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 |
9 | // |
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. |
17 | // |
18 | #ifndef BOOST_MATH_BESSEL_JY_ZERO_2013_01_18_HPP_ |
19 | #define BOOST_MATH_BESSEL_JY_ZERO_2013_01_18_HPP_ |
20 | |
21 | #include <algorithm> |
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> |
26 | |
27 | namespace boost { namespace math { |
28 | namespace detail |
29 | { |
30 | namespace bessel_zero |
31 | { |
32 | template<class T> |
33 | T equation_nist_10_21_19(const T& v, const T& a) |
34 | { |
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. |
38 | |
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; |
44 | |
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; |
48 | |
49 | return a + (((( - term7 |
50 | * eight_a_inv_squared - term5) |
51 | * eight_a_inv_squared - term3) |
52 | * eight_a_inv_squared - mu_minus_one) |
53 | * eight_a_inv); |
54 | } |
55 | |
56 | template<typename T> |
57 | class equation_as_9_3_39_and_its_derivative |
58 | { |
59 | public: |
60 | equation_as_9_3_39_and_its_derivative(const T& zt) : zeta(zt) { } |
61 | |
62 | boost::math::tuple<T, T> operator()(const T& z) const |
63 | { |
64 | BOOST_MATH_STD_USING // ADL of std names, needed for acos, sqrt. |
65 | |
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. |
69 | |
70 | const T zsq_minus_one_sqrt = sqrt((z * z) - T(1)); |
71 | |
72 | const T the_function( |
73 | zsq_minus_one_sqrt |
74 | - ( acos(T(1) / z) + ((T(2) / 3U) * (zeta * sqrt(zeta))))); |
75 | |
76 | const T its_derivative(zsq_minus_one_sqrt / z); |
77 | |
78 | return boost::math::tuple<T, T>(the_function, its_derivative); |
79 | } |
80 | |
81 | private: |
82 | const equation_as_9_3_39_and_its_derivative& operator=(const equation_as_9_3_39_and_its_derivative&); |
83 | const T zeta; |
84 | }; |
85 | |
86 | template<class T> |
87 | static T equation_as_9_5_26(const T& v, const T& ai_bi_root) |
88 | { |
89 | BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt. |
90 | |
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. |
95 | // |
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. |
103 | // |
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. |
107 | |
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)); |
110 | |
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); |
114 | |
115 | const T zeta_sqrt = sqrt(zeta); |
116 | |
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>()); |
120 | |
121 | // Solve the quadratic equation, taking the positive root. |
122 | const T z_estimate = (-b + sqrt((b * b) - T(2))) / 2U; |
123 | |
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); |
128 | |
129 | const int my_digits10 = static_cast<int>(static_cast<float>(boost::math::tools::digits<T>() * 0.301F)); |
130 | |
131 | // Select the maximum allowed iterations based on the number |
132 | // of decimal digits in the numeric type T, being at least 12. |
133 | const boost::uintmax_t iterations_allowed = static_cast<boost::uintmax_t>((std::max)(a: 12, b: my_digits10 * 2)); |
134 | |
135 | boost::uintmax_t iterations_used = iterations_allowed; |
136 | |
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), |
140 | z_estimate, |
141 | range_zmin, |
142 | range_zmax, |
143 | (std::min)(boost::math::tools::digits<T>(), boost::math::tools::digits<float>()), |
144 | iterations_used); |
145 | |
146 | static_cast<void>(iterations_used); |
147 | |
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); |
151 | |
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); |
156 | |
157 | const T b0 = -b0_term_5_48 + ((b0_term_5_24 + b0_term_1_8) / zeta_sqrt); |
158 | |
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; |
161 | |
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); |
164 | } |
165 | |
166 | namespace cyl_bessel_j_zero_detail |
167 | { |
168 | template<class T> |
169 | T equation_nist_10_21_40_a(const T& v) |
170 | { |
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)); |
173 | |
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)); |
180 | } |
181 | |
182 | template<class T, class Policy> |
183 | class function_object_jv |
184 | { |
185 | public: |
186 | function_object_jv(const T& v, |
187 | const Policy& pol) : my_v(v), |
188 | my_pol(pol) { } |
189 | |
190 | T operator()(const T& x) const |
191 | { |
192 | return boost::math::cyl_bessel_j(my_v, x, my_pol); |
193 | } |
194 | |
195 | private: |
196 | const T my_v; |
197 | const Policy& my_pol; |
198 | const function_object_jv& operator=(const function_object_jv&); |
199 | }; |
200 | |
201 | template<class T, class Policy> |
202 | class function_object_jv_and_jv_prime |
203 | { |
204 | public: |
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), |
209 | my_pol(pol) { } |
210 | |
211 | boost::math::tuple<T, T> operator()(const T& x) const |
212 | { |
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. |
217 | T j_v; |
218 | T j_v_prime; |
219 | |
220 | if(my_order_is_zero) |
221 | { |
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); |
224 | } |
225 | else |
226 | { |
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); |
230 | } |
231 | |
232 | // Return a tuple containing both Jv(x) and Jv'(x). |
233 | return boost::math::make_tuple(j_v, j_v_prime); |
234 | } |
235 | |
236 | private: |
237 | const T my_v; |
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&); |
241 | }; |
242 | |
243 | template<class T> bool my_bisection_unreachable_tolerance(const T&, const T&) { return false; } |
244 | |
245 | template<class T, class Policy> |
246 | T initial_guess(const T& v, const int m, const Policy& pol) |
247 | { |
248 | BOOST_MATH_STD_USING // ADL of std names, needed for floor. |
249 | |
250 | // Compute an estimate of the m'th root of cyl_bessel_j. |
251 | |
252 | T guess; |
253 | |
254 | // There is special handling for negative order. |
255 | if(v < 0) |
256 | { |
257 | if((m == 1) && (v > -0.5F)) |
258 | { |
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}] |
262 | // N[%, 20] |
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); |
271 | |
272 | return guess; |
273 | } |
274 | |
275 | // Create the positive order and extract its positive floor integer part. |
276 | const T vv(-v); |
277 | const T vv_floor(floor(vv)); |
278 | |
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. |
282 | |
283 | T root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(vv_floor, m, pol); |
284 | T root_lo; |
285 | |
286 | if(m == 1) |
287 | { |
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); |
291 | |
292 | const bool hi_end_of_bracket_is_negative = (boost::math::cyl_bessel_j(v, root_hi, pol) < 0); |
293 | |
294 | while((root_lo > boost::math::tools::epsilon<T>())) |
295 | { |
296 | const bool lo_end_of_bracket_is_negative = (boost::math::cyl_bessel_j(v, root_lo, pol) < 0); |
297 | |
298 | if(hi_end_of_bracket_is_negative != lo_end_of_bracket_is_negative) |
299 | { |
300 | break; |
301 | } |
302 | |
303 | root_hi = root_lo; |
304 | |
305 | // Decrease the lower end of the bracket using an adaptive algorithm. |
306 | if(root_lo > 0.5F) |
307 | { |
308 | root_lo -= 0.5F; |
309 | } |
310 | else |
311 | { |
312 | root_lo *= 0.75F; |
313 | } |
314 | } |
315 | } |
316 | else |
317 | { |
318 | root_lo = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(vv_floor, m - 1, pol); |
319 | } |
320 | |
321 | // Perform several steps of bisection iteration to refine the guess. |
322 | boost::uintmax_t number_of_iterations(12U); |
323 | |
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), |
328 | root_lo, |
329 | root_hi, |
330 | boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::my_bisection_unreachable_tolerance<T>, |
331 | number_of_iterations); |
332 | |
333 | return (boost::math::get<0>(guess_pair) + boost::math::get<1>(guess_pair)) / 2U; |
334 | } |
335 | |
336 | if(m == 1U) |
337 | { |
338 | // Get the initial estimate of the first root. |
339 | |
340 | if(v < 2.2F) |
341 | { |
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}] |
345 | // N[%, 20] |
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); |
354 | } |
355 | else |
356 | { |
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); |
359 | } |
360 | } |
361 | else |
362 | { |
363 | if(v < 2.2F) |
364 | { |
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>()); |
367 | |
368 | guess = boost::math::detail::bessel_zero::equation_nist_10_21_19(v, a); |
369 | } |
370 | else |
371 | { |
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)); |
374 | |
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); |
377 | } |
378 | } |
379 | |
380 | return guess; |
381 | } |
382 | } // namespace cyl_bessel_j_zero_detail |
383 | |
384 | namespace cyl_neumann_zero_detail |
385 | { |
386 | template<class T> |
387 | T equation_nist_10_21_40_b(const T& v) |
388 | { |
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)); |
391 | |
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)); |
398 | } |
399 | |
400 | template<class T, class Policy> |
401 | class function_object_yv |
402 | { |
403 | public: |
404 | function_object_yv(const T& v, |
405 | const Policy& pol) : my_v(v), |
406 | my_pol(pol) { } |
407 | |
408 | T operator()(const T& x) const |
409 | { |
410 | return boost::math::cyl_neumann(my_v, x, my_pol); |
411 | } |
412 | |
413 | private: |
414 | const T my_v; |
415 | const Policy& my_pol; |
416 | const function_object_yv& operator=(const function_object_yv&); |
417 | }; |
418 | |
419 | template<class T, class Policy> |
420 | class function_object_yv_and_yv_prime |
421 | { |
422 | public: |
423 | function_object_yv_and_yv_prime(const T& v, |
424 | const Policy& pol) : my_v(v), |
425 | my_pol(pol) { } |
426 | |
427 | boost::math::tuple<T, T> operator()(const T& x) const |
428 | { |
429 | const T half_epsilon(boost::math::tools::epsilon<T>() / 2U); |
430 | |
431 | const bool order_is_zero = ((my_v > -half_epsilon) && (my_v < +half_epsilon)); |
432 | |
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. |
437 | T y_v; |
438 | T y_v_prime; |
439 | |
440 | if(order_is_zero) |
441 | { |
442 | y_v = boost::math::cyl_neumann(0, x, my_pol); |
443 | y_v_prime = -boost::math::cyl_neumann(1, x, my_pol); |
444 | } |
445 | else |
446 | { |
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); |
450 | } |
451 | |
452 | // Return a tuple containing both Yv(x) and Yv'(x). |
453 | return boost::math::make_tuple(y_v, y_v_prime); |
454 | } |
455 | |
456 | private: |
457 | const T my_v; |
458 | const Policy& my_pol; |
459 | const function_object_yv_and_yv_prime& operator=(const function_object_yv_and_yv_prime&); |
460 | }; |
461 | |
462 | template<class T> bool my_bisection_unreachable_tolerance(const T&, const T&) { return false; } |
463 | |
464 | template<class T, class Policy> |
465 | T initial_guess(const T& v, const int m, const Policy& pol) |
466 | { |
467 | BOOST_MATH_STD_USING // ADL of std names, needed for floor. |
468 | |
469 | // Compute an estimate of the m'th root of cyl_neumann. |
470 | |
471 | T guess; |
472 | |
473 | // There is special handling for negative order. |
474 | if(v < 0) |
475 | { |
476 | // Create the positive order and extract its positive floor and ceiling integer parts. |
477 | const T vv(-v); |
478 | const T vv_floor(floor(vv)); |
479 | |
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. |
483 | |
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) |
488 | // for v = -n/2. |
489 | |
490 | T root_hi; |
491 | T root_lo; |
492 | |
493 | if(m == 1) |
494 | { |
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) |
500 | { |
501 | root_hi = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m, pol); |
502 | } |
503 | else |
504 | { |
505 | root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m, pol); |
506 | } |
507 | |
508 | root_lo = T(root_hi - 0.1F); |
509 | |
510 | const bool hi_end_of_bracket_is_negative = (boost::math::cyl_neumann(v, root_hi, pol) < 0); |
511 | |
512 | while((root_lo > boost::math::tools::epsilon<T>())) |
513 | { |
514 | const bool lo_end_of_bracket_is_negative = (boost::math::cyl_neumann(v, root_lo, pol) < 0); |
515 | |
516 | if(hi_end_of_bracket_is_negative != lo_end_of_bracket_is_negative) |
517 | { |
518 | break; |
519 | } |
520 | |
521 | root_hi = root_lo; |
522 | |
523 | // Decrease the lower end of the bracket using an adaptive algorithm. |
524 | if(root_lo > 0.5F) |
525 | { |
526 | root_lo -= 0.5F; |
527 | } |
528 | else |
529 | { |
530 | root_lo *= 0.75F; |
531 | } |
532 | } |
533 | } |
534 | else |
535 | { |
536 | if(T(vv - vv_floor) < 0.5F) |
537 | { |
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); |
540 | root_lo += 0.01F; |
541 | root_hi += 0.01F; |
542 | } |
543 | else |
544 | { |
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); |
547 | root_lo += 0.01F; |
548 | root_hi += 0.01F; |
549 | } |
550 | } |
551 | |
552 | // Perform several steps of bisection iteration to refine the guess. |
553 | boost::uintmax_t number_of_iterations(12U); |
554 | |
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), |
559 | root_lo, |
560 | root_hi, |
561 | boost::math::detail::bessel_zero::cyl_neumann_zero_detail::my_bisection_unreachable_tolerance<T>, |
562 | number_of_iterations); |
563 | |
564 | return (boost::math::get<0>(guess_pair) + boost::math::get<1>(guess_pair)) / 2U; |
565 | } |
566 | |
567 | if(m == 1U) |
568 | { |
569 | // Get the initial estimate of the first root. |
570 | |
571 | if(v < 2.2F) |
572 | { |
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}] |
576 | // N[%, 20] |
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); |
585 | } |
586 | else |
587 | { |
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); |
590 | } |
591 | } |
592 | else |
593 | { |
594 | if(v < 2.2F) |
595 | { |
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>()); |
598 | |
599 | guess = boost::math::detail::bessel_zero::equation_nist_10_21_19(v, a); |
600 | } |
601 | else |
602 | { |
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)); |
605 | |
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); |
608 | } |
609 | } |
610 | |
611 | return guess; |
612 | } |
613 | } // namespace cyl_neumann_zero_detail |
614 | } // namespace bessel_zero |
615 | } } } // namespace boost::math::detail |
616 | |
617 | #endif // BOOST_MATH_BESSEL_JY_ZERO_2013_01_18_HPP_ |
618 | |