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 explicit equation_as_9_3_39_and_its_derivative(const T& zt) : zeta(zt) { }
61
62 equation_as_9_3_39_and_its_derivative(const equation_as_9_3_39_and_its_derivative&) = default;
63
64 boost::math::tuple<T, T> operator()(const T& z) const
65 {
66 BOOST_MATH_STD_USING // ADL of std names, needed for acos, sqrt.
67
68 // Return the function of zeta that is implicitly defined
69 // in A&S Eq. 9.3.39 as a function of z. The function is
70 // returned along with its derivative with respect to z.
71
72 const T zsq_minus_one_sqrt = sqrt((z * z) - T(1));
73
74 const T the_function(
75 zsq_minus_one_sqrt
76 - ( acos(T(1) / z) + ((T(2) / 3U) * (zeta * sqrt(zeta)))));
77
78 const T its_derivative(zsq_minus_one_sqrt / z);
79
80 return boost::math::tuple<T, T>(the_function, its_derivative);
81 }
82
83 private:
84 const equation_as_9_3_39_and_its_derivative& operator=(const equation_as_9_3_39_and_its_derivative&) = delete;
85 const T zeta;
86 };
87
88 template<class T, class Policy>
89 static T equation_as_9_5_26(const T& v, const T& ai_bi_root, const Policy& pol)
90 {
91 BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt.
92
93 // Obtain the estimate of the m'th zero of Jv or Yv.
94 // The order m has been used to create the input parameter ai_bi_root.
95 // Here, v is larger than about 2.2. The estimate is computed
96 // from Abramowitz and Stegun Eqs. 9.5.22 and 9.5.26, page 371.
97 //
98 // The inversion of z as a function of zeta is mentioned in the text
99 // following A&S Eq. 9.5.26. Here, we accomplish the inversion by
100 // performing a Taylor expansion of Eq. 9.3.39 for large z to order 2
101 // and solving the resulting quadratic equation, thereby taking
102 // the positive root of the quadratic.
103 // In other words: (2/3)(-zeta)^(3/2) approx = z + 1/(2z) - pi/2.
104 // This leads to: z^2 - [(2/3)(-zeta)^(3/2) + pi/2]z + 1/2 = 0.
105 //
106 // With this initial estimate, Newton-Raphson iteration is used
107 // to refine the value of the estimate of the root of z
108 // as a function of zeta.
109
110 const T v_pow_third(boost::math::cbrt(v, pol));
111 const T v_pow_minus_two_thirds(T(1) / (v_pow_third * v_pow_third));
112
113 // Obtain zeta using the order v combined with the m'th root of
114 // an airy function, as shown in A&S Eq. 9.5.22.
115 const T zeta = v_pow_minus_two_thirds * (-ai_bi_root);
116
117 const T zeta_sqrt = sqrt(zeta);
118
119 // Set up a quadratic equation based on the Taylor series
120 // expansion mentioned above.
121 const T b = -((((zeta * zeta_sqrt) * 2U) / 3U) + boost::math::constants::half_pi<T>());
122
123 // Solve the quadratic equation, taking the positive root.
124 const T z_estimate = (-b + sqrt((b * b) - T(2))) / 2U;
125
126 // Establish the range, the digits, and the iteration limit
127 // for the upcoming root-finding.
128 const T range_zmin = (std::max<T>)(z_estimate - T(1), T(1));
129 const T range_zmax = z_estimate + T(1);
130
131 const auto my_digits10 = static_cast<int>(static_cast<float>(boost::math::tools::digits<T>() * 0.301F));
132
133 // Select the maximum allowed iterations based on the number
134 // of decimal digits in the numeric type T, being at least 12.
135 const auto iterations_allowed = static_cast<std::uintmax_t>((std::max)(a: 12, b: my_digits10 * 2));
136
137 std::uintmax_t iterations_used = iterations_allowed;
138
139 // Calculate the root of z as a function of zeta.
140 const T z = boost::math::tools::newton_raphson_iterate(
141 boost::math::detail::bessel_zero::equation_as_9_3_39_and_its_derivative<T>(zeta),
142 z_estimate,
143 range_zmin,
144 range_zmax,
145 (std::min)(boost::math::tools::digits<T>(), boost::math::tools::digits<float>()),
146 iterations_used);
147
148 static_cast<void>(iterations_used);
149
150 // Continue with the implementation of A&S Eq. 9.3.39.
151 const T zsq_minus_one = (z * z) - T(1);
152 const T zsq_minus_one_sqrt = sqrt(zsq_minus_one);
153
154 // This is A&S Eq. 9.3.42.
155 const T b0_term_5_24 = T(5) / ((zsq_minus_one * zsq_minus_one_sqrt) * 24U);
156 const T b0_term_1_8 = T(1) / ( zsq_minus_one_sqrt * 8U);
157 const T b0_term_5_48 = T(5) / ((zeta * zeta) * 48U);
158
159 const T b0 = -b0_term_5_48 + ((b0_term_5_24 + b0_term_1_8) / zeta_sqrt);
160
161 // This is the second line of A&S Eq. 9.5.26 for f_k with k = 1.
162 const T f1 = ((z * zeta_sqrt) * b0) / zsq_minus_one_sqrt;
163
164 // This is A&S Eq. 9.5.22 expanded to k = 1 (i.e., one term in the series).
165 return (v * z) + (f1 / v);
166 }
167
168 namespace cyl_bessel_j_zero_detail
169 {
170 template<class T, class Policy>
171 T equation_nist_10_21_40_a(const T& v, const Policy& pol)
172 {
173 const T v_pow_third(boost::math::cbrt(v, pol));
174 const T v_pow_minus_two_thirds(T(1) / (v_pow_third * v_pow_third));
175
176 return v * ((((( + T(0.043)
177 * v_pow_minus_two_thirds - T(0.0908))
178 * v_pow_minus_two_thirds - T(0.00397))
179 * v_pow_minus_two_thirds + T(1.033150))
180 * v_pow_minus_two_thirds + T(1.8557571))
181 * v_pow_minus_two_thirds + T(1));
182 }
183
184 template<class T, class Policy>
185 class function_object_jv
186 {
187 public:
188 function_object_jv(const T& v,
189 const Policy& pol) : my_v(v),
190 my_pol(pol) { }
191
192 function_object_jv(const function_object_jv&) = default;
193
194 T operator()(const T& x) const
195 {
196 return boost::math::cyl_bessel_j(my_v, x, my_pol);
197 }
198
199 private:
200 const T my_v;
201 const Policy& my_pol;
202 const function_object_jv& operator=(const function_object_jv&) = delete;
203 };
204
205 template<class T, class Policy>
206 class function_object_jv_and_jv_prime
207 {
208 public:
209 function_object_jv_and_jv_prime(const T& v,
210 const bool order_is_zero,
211 const Policy& pol) : my_v(v),
212 my_order_is_zero(order_is_zero),
213 my_pol(pol) { }
214
215 function_object_jv_and_jv_prime(const function_object_jv_and_jv_prime&) = default;
216
217 boost::math::tuple<T, T> operator()(const T& x) const
218 {
219 // Obtain Jv(x) and Jv'(x).
220 // Chris's original code called the Bessel function implementation layer direct,
221 // but that circumvented optimizations for integer-orders. Call the documented
222 // top level functions instead, and let them sort out which implementation to use.
223 T j_v;
224 T j_v_prime;
225
226 if(my_order_is_zero)
227 {
228 j_v = boost::math::cyl_bessel_j(0, x, my_pol);
229 j_v_prime = -boost::math::cyl_bessel_j(1, x, my_pol);
230 }
231 else
232 {
233 j_v = boost::math::cyl_bessel_j( my_v, x, my_pol);
234 const T j_v_m1 (boost::math::cyl_bessel_j(T(my_v - 1), x, my_pol));
235 j_v_prime = j_v_m1 - ((my_v * j_v) / x);
236 }
237
238 // Return a tuple containing both Jv(x) and Jv'(x).
239 return boost::math::make_tuple(j_v, j_v_prime);
240 }
241
242 private:
243 const T my_v;
244 const bool my_order_is_zero;
245 const Policy& my_pol;
246 const function_object_jv_and_jv_prime& operator=(const function_object_jv_and_jv_prime&) = delete;
247 };
248
249 template<class T> bool my_bisection_unreachable_tolerance(const T&, const T&) { return false; }
250
251 template<class T, class Policy>
252 T initial_guess(const T& v, const int m, const Policy& pol)
253 {
254 BOOST_MATH_STD_USING // ADL of std names, needed for floor.
255
256 // Compute an estimate of the m'th root of cyl_bessel_j.
257
258 T guess;
259
260 // There is special handling for negative order.
261 if(v < 0)
262 {
263 if((m == 1) && (v > -0.5F))
264 {
265 // For small, negative v, use the results of empirical curve fitting.
266 // Mathematica(R) session for the coefficients:
267 // Table[{n, BesselJZero[n, 1]}, {n, -(1/2), 0, 1/10}]
268 // N[%, 20]
269 // Fit[%, {n^0, n^1, n^2, n^3, n^4, n^5, n^6}, n]
270 guess = ((((( - T(0.2321156900729)
271 * v - T(0.1493247777488))
272 * v - T(0.15205419167239))
273 * v + T(0.07814930561249))
274 * v - T(0.17757573537688))
275 * v + T(1.542805677045663))
276 * v + T(2.40482555769577277);
277
278 return guess;
279 }
280
281 // Create the positive order and extract its positive floor integer part.
282 const T vv(-v);
283 const T vv_floor(floor(vv));
284
285 // The to-be-found root is bracketed by the roots of the
286 // Bessel function whose reflected, positive integer order
287 // is less than, but nearest to vv.
288
289 T root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(vv_floor, m, pol);
290 T root_lo;
291
292 if(m == 1)
293 {
294 // The estimate of the first root for negative order is found using
295 // an adaptive range-searching algorithm.
296 root_lo = T(root_hi - 0.1F);
297
298 const bool hi_end_of_bracket_is_negative = (boost::math::cyl_bessel_j(v, root_hi, pol) < 0);
299
300 while((root_lo > boost::math::tools::epsilon<T>()))
301 {
302 const bool lo_end_of_bracket_is_negative = (boost::math::cyl_bessel_j(v, root_lo, pol) < 0);
303
304 if(hi_end_of_bracket_is_negative != lo_end_of_bracket_is_negative)
305 {
306 break;
307 }
308
309 root_hi = root_lo;
310
311 // Decrease the lower end of the bracket using an adaptive algorithm.
312 if(root_lo > 0.5F)
313 {
314 root_lo -= 0.5F;
315 }
316 else
317 {
318 root_lo *= 0.75F; // LCOV_EXCL_LINE probably unreachable, but hard to prove?
319 }
320 }
321 }
322 else
323 {
324 root_lo = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(vv_floor, m - 1, pol);
325 }
326
327 // Perform several steps of bisection iteration to refine the guess.
328 std::uintmax_t number_of_iterations(12U);
329
330 // Do the bisection iteration.
331 const boost::math::tuple<T, T> guess_pair =
332 boost::math::tools::bisect(
333 boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::function_object_jv<T, Policy>(v, pol),
334 root_lo,
335 root_hi,
336 boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::my_bisection_unreachable_tolerance<T>,
337 number_of_iterations);
338
339 return (boost::math::get<0>(guess_pair) + boost::math::get<1>(guess_pair)) / 2U;
340 }
341
342 if(m == 1U)
343 {
344 // Get the initial estimate of the first root.
345
346 if(v < 2.2F)
347 {
348 // For small v, use the results of empirical curve fitting.
349 // Mathematica(R) session for the coefficients:
350 // Table[{n, BesselJZero[n, 1]}, {n, 0, 22/10, 1/10}]
351 // N[%, 20]
352 // Fit[%, {n^0, n^1, n^2, n^3, n^4, n^5, n^6}, n]
353 guess = ((((( - T(0.0008342379046010)
354 * v + T(0.007590035637410))
355 * v - T(0.030640914772013))
356 * v + T(0.078232088020106))
357 * v - T(0.169668712590620))
358 * v + T(1.542187960073750))
359 * v + T(2.4048359915254634);
360 }
361 else
362 {
363 // For larger v, use the first line of Eqs. 10.21.40 in the NIST Handbook.
364 guess = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::equation_nist_10_21_40_a(v, pol);
365 }
366 }
367 else
368 {
369 if(v < 2.2F)
370 {
371 // Use Eq. 10.21.19 in the NIST Handbook.
372 const T a(((v + T(m * 2U)) - T(0.5)) * boost::math::constants::half_pi<T>());
373
374 guess = boost::math::detail::bessel_zero::equation_nist_10_21_19(v, a);
375 }
376 else
377 {
378 // Get an estimate of the m'th root of airy_ai.
379 const T airy_ai_root(boost::math::detail::airy_zero::airy_ai_zero_detail::initial_guess<T>(m, pol));
380
381 // Use Eq. 9.5.26 in the A&S Handbook.
382 guess = boost::math::detail::bessel_zero::equation_as_9_5_26(v, airy_ai_root, pol);
383 }
384 }
385
386 return guess;
387 }
388 } // namespace cyl_bessel_j_zero_detail
389
390 namespace cyl_neumann_zero_detail
391 {
392 template<class T, class Policy>
393 T equation_nist_10_21_40_b(const T& v, const Policy& pol)
394 {
395 const T v_pow_third(boost::math::cbrt(v, pol));
396 const T v_pow_minus_two_thirds(T(1) / (v_pow_third * v_pow_third));
397
398 return v * ((((( - T(0.001)
399 * v_pow_minus_two_thirds - T(0.0060))
400 * v_pow_minus_two_thirds + T(0.01198))
401 * v_pow_minus_two_thirds + T(0.260351))
402 * v_pow_minus_two_thirds + T(0.9315768))
403 * v_pow_minus_two_thirds + T(1));
404 }
405
406 template<class T, class Policy>
407 class function_object_yv
408 {
409 public:
410 function_object_yv(const T& v,
411 const Policy& pol) : my_v(v),
412 my_pol(pol) { }
413
414 function_object_yv(const function_object_yv&) = default;
415
416 T operator()(const T& x) const
417 {
418 return boost::math::cyl_neumann(my_v, x, my_pol);
419 }
420
421 private:
422 const T my_v;
423 const Policy& my_pol;
424 const function_object_yv& operator=(const function_object_yv&) = delete;
425 };
426
427 template<class T, class Policy>
428 class function_object_yv_and_yv_prime
429 {
430 public:
431 function_object_yv_and_yv_prime(const T& v,
432 const Policy& pol) : my_v(v),
433 my_pol(pol) { }
434
435 function_object_yv_and_yv_prime(const function_object_yv_and_yv_prime&) = default;
436
437 boost::math::tuple<T, T> operator()(const T& x) const
438 {
439 const T half_epsilon(boost::math::tools::epsilon<T>() / 2U);
440
441 const bool order_is_zero = ((my_v > -half_epsilon) && (my_v < +half_epsilon));
442
443 // Obtain Yv(x) and Yv'(x).
444 // Chris's original code called the Bessel function implementation layer direct,
445 // but that circumvented optimizations for integer-orders. Call the documented
446 // top level functions instead, and let them sort out which implementation to use.
447 T y_v;
448 T y_v_prime;
449
450 if(order_is_zero)
451 {
452 y_v = boost::math::cyl_neumann(0, x, my_pol);
453 y_v_prime = -boost::math::cyl_neumann(1, x, my_pol);
454 }
455 else
456 {
457 y_v = boost::math::cyl_neumann( my_v, x, my_pol);
458 const T y_v_m1 (boost::math::cyl_neumann(T(my_v - 1), x, my_pol));
459 y_v_prime = y_v_m1 - ((my_v * y_v) / x);
460 }
461
462 // Return a tuple containing both Yv(x) and Yv'(x).
463 return boost::math::make_tuple(y_v, y_v_prime);
464 }
465
466 private:
467 const T my_v;
468 const Policy& my_pol;
469 const function_object_yv_and_yv_prime& operator=(const function_object_yv_and_yv_prime&) = delete;
470 };
471
472 template<class T> bool my_bisection_unreachable_tolerance(const T&, const T&) { return false; }
473
474 template<class T, class Policy>
475 T initial_guess(const T& v, const int m, const Policy& pol)
476 {
477 BOOST_MATH_STD_USING // ADL of std names, needed for floor.
478
479 // Compute an estimate of the m'th root of cyl_neumann.
480
481 T guess;
482
483 // There is special handling for negative order.
484 if(v < 0)
485 {
486 // Create the positive order and extract its positive floor and ceiling integer parts.
487 const T vv(-v);
488 const T vv_floor(floor(vv));
489
490 // The to-be-found root is bracketed by the roots of the
491 // Bessel function whose reflected, positive integer order
492 // is less than, but nearest to vv.
493
494 // The special case of negative, half-integer order uses
495 // the relation between Yv and spherical Bessel functions
496 // in order to obtain the bracket for the root.
497 // In these special cases, cyl_neumann(-n/2, x) = sph_bessel_j(+n/2, x)
498 // for v = -n/2.
499
500 T root_hi;
501 T root_lo;
502
503 if(m == 1)
504 {
505 // The estimate of the first root for negative order is found using
506 // an adaptive range-searching algorithm.
507 // Take special precautions for the discontinuity at negative,
508 // half-integer orders and use different brackets above and below these.
509 if(T(vv - vv_floor) < 0.5F)
510 {
511 root_hi = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m, pol);
512 }
513 else
514 {
515 root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m, pol);
516 }
517
518 root_lo = T(root_hi - 0.1F);
519
520 const bool hi_end_of_bracket_is_negative = (boost::math::cyl_neumann(v, root_hi, pol) < 0);
521
522 while((root_lo > boost::math::tools::epsilon<T>()))
523 {
524 const bool lo_end_of_bracket_is_negative = (boost::math::cyl_neumann(v, root_lo, pol) < 0);
525
526 if(hi_end_of_bracket_is_negative != lo_end_of_bracket_is_negative)
527 {
528 break;
529 }
530
531 root_hi = root_lo;
532
533 // Decrease the lower end of the bracket using an adaptive algorithm.
534 if(root_lo > 0.5F)
535 {
536 root_lo -= 0.5F;
537 }
538 else
539 {
540 root_lo *= 0.75F; // LCOV_EXCL_LINE probably unreachable, but hard to prove?
541 }
542 }
543 }
544 else
545 {
546 if(T(vv - vv_floor) < 0.5F)
547 {
548 root_lo = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m - 1, pol);
549 root_hi = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m, pol);
550 root_lo += 0.01F;
551 root_hi += 0.01F;
552 }
553 else
554 {
555 root_lo = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m - 1, pol);
556 root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m, pol);
557 root_lo += 0.01F;
558 root_hi += 0.01F;
559 }
560 }
561
562 // Perform several steps of bisection iteration to refine the guess.
563 std::uintmax_t number_of_iterations(12U);
564
565 // Do the bisection iteration.
566 const boost::math::tuple<T, T> guess_pair =
567 boost::math::tools::bisect(
568 boost::math::detail::bessel_zero::cyl_neumann_zero_detail::function_object_yv<T, Policy>(v, pol),
569 root_lo,
570 root_hi,
571 boost::math::detail::bessel_zero::cyl_neumann_zero_detail::my_bisection_unreachable_tolerance<T>,
572 number_of_iterations);
573
574 return (boost::math::get<0>(guess_pair) + boost::math::get<1>(guess_pair)) / 2U;
575 }
576
577 if(m == 1U)
578 {
579 // Get the initial estimate of the first root.
580
581 if(v < 2.2F)
582 {
583 // For small v, use the results of empirical curve fitting.
584 // Mathematica(R) session for the coefficients:
585 // Table[{n, BesselYZero[n, 1]}, {n, 0, 22/10, 1/10}]
586 // N[%, 20]
587 // Fit[%, {n^0, n^1, n^2, n^3, n^4, n^5, n^6}, n]
588 guess = ((((( - T(0.0025095909235652)
589 * v + T(0.021291887049053))
590 * v - T(0.076487785486526))
591 * v + T(0.159110268115362))
592 * v - T(0.241681668765196))
593 * v + T(1.4437846310885244))
594 * v + T(0.89362115190200490);
595 }
596 else
597 {
598 // For larger v, use the second line of Eqs. 10.21.40 in the NIST Handbook.
599 guess = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::equation_nist_10_21_40_b(v, pol);
600 }
601 }
602 else
603 {
604 if(v < 2.2F)
605 {
606 // Use Eq. 10.21.19 in the NIST Handbook.
607 const T a(((v + T(m * 2U)) - T(1.5)) * boost::math::constants::half_pi<T>());
608
609 guess = boost::math::detail::bessel_zero::equation_nist_10_21_19(v, a);
610 }
611 else
612 {
613 // Get an estimate of the m'th root of airy_bi.
614 const T airy_bi_root(boost::math::detail::airy_zero::airy_bi_zero_detail::initial_guess<T>(m, pol));
615
616 // Use Eq. 9.5.26 in the A&S Handbook.
617 guess = boost::math::detail::bessel_zero::equation_as_9_5_26(v, airy_bi_root, pol);
618 }
619 }
620
621 return guess;
622 }
623 } // namespace cyl_neumann_zero_detail
624 } // namespace bessel_zero
625 } } } // namespace boost::math::detail
626
627#endif // BOOST_MATH_BESSEL_JY_ZERO_2013_01_18_HPP_
628

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