]>
Commit | Line | Data |
---|---|---|
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 _BESSEL_JY_ZERO_2013_01_18_HPP_ | |
19 | #define _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)(12, 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 // _BESSEL_JY_ZERO_2013_01_18_HPP_ |