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

source code of include/boost/math/special_functions/detail/bessel_jy_zero.hpp