diff options
Diffstat (limited to 'boost/math/special_functions/detail/bessel_jy_zero.hpp')
-rw-r--r-- | boost/math/special_functions/detail/bessel_jy_zero.hpp | 617 |
1 files changed, 617 insertions, 0 deletions
diff --git a/boost/math/special_functions/detail/bessel_jy_zero.hpp b/boost/math/special_functions/detail/bessel_jy_zero.hpp new file mode 100644 index 0000000000..ecd8696eee --- /dev/null +++ b/boost/math/special_functions/detail/bessel_jy_zero.hpp @@ -0,0 +1,617 @@ +// Copyright (c) 2013 Christopher Kormanyos +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +// +// This work is based on an earlier work: +// "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations", +// in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469 +// +// This header contains implementation details for estimating the zeros +// of cylindrical Bessel and Neumann functions on the positive real axis. +// Support is included for both positive as well as negative order. +// Various methods are used to estimate the roots. These include +// empirical curve fitting and McMahon's asymptotic approximation +// for small order, uniform asymptotic expansion for large order, +// and iteration and root interlacing for negative order. +// +#ifndef _BESSEL_JY_ZERO_2013_01_18_HPP_ + #define _BESSEL_JY_ZERO_2013_01_18_HPP_ + + #include <algorithm> + #include <boost/math/constants/constants.hpp> + #include <boost/math/special_functions/math_fwd.hpp> + #include <boost/math/special_functions/cbrt.hpp> + #include <boost/math/special_functions/detail/airy_ai_bi_zero.hpp> + + namespace boost { namespace math { + namespace detail + { + namespace bessel_zero + { + template<class T> + T equation_nist_10_21_19(const T& v, const T& a) + { + // Get the initial estimate of the m'th root of Jv or Yv. + // This subroutine is used for the order m with m > 1. + // The order m has been used to create the input parameter a. + + // This is Eq. 10.21.19 in the NIST Handbook. + const T mu = (v * v) * 4U; + const T mu_minus_one = mu - T(1); + const T eight_a_inv = T(1) / (a * 8U); + const T eight_a_inv_squared = eight_a_inv * eight_a_inv; + + const T term3 = ((mu_minus_one * 4U) * ((mu * 7U) - T(31U) )) / 3U; + const T term5 = ((mu_minus_one * 32U) * ((((mu * 83U) - T(982U) ) * mu) + T(3779U) )) / 15U; + const T term7 = ((mu_minus_one * 64U) * ((((((mu * 6949U) - T(153855UL)) * mu) + T(1585743UL)) * mu) - T(6277237UL))) / 105U; + + return a + (((( - term7 + * eight_a_inv_squared - term5) + * eight_a_inv_squared - term3) + * eight_a_inv_squared - mu_minus_one) + * eight_a_inv); + } + + template<typename T> + class equation_as_9_3_39_and_its_derivative + { + public: + equation_as_9_3_39_and_its_derivative(const T& zt) : zeta(zt) { } + + boost::math::tuple<T, T> operator()(const T& z) const + { + BOOST_MATH_STD_USING // ADL of std names, needed for acos, sqrt. + + // Return the function of zeta that is implicitly defined + // in A&S Eq. 9.3.39 as a function of z. The function is + // returned along with its derivative with respect to z. + + const T zsq_minus_one_sqrt = sqrt((z * z) - T(1)); + + const T the_function( + zsq_minus_one_sqrt + - ( acos(T(1) / z) + ((T(2) / 3U) * (zeta * sqrt(zeta))))); + + const T its_derivative(zsq_minus_one_sqrt / z); + + return boost::math::tuple<T, T>(the_function, its_derivative); + } + + private: + const equation_as_9_3_39_and_its_derivative& operator=(const equation_as_9_3_39_and_its_derivative&); + const T zeta; + }; + + template<class T> + static T equation_as_9_5_26(const T& v, const T& ai_bi_root) + { + BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt. + + // Obtain the estimate of the m'th zero of Jv or Yv. + // The order m has been used to create the input parameter ai_bi_root. + // Here, v is larger than about 2.2. The estimate is computed + // from Abramowitz and Stegun Eqs. 9.5.22 and 9.5.26, page 371. + // + // The inversion of z as a function of zeta is mentioned in the text + // following A&S Eq. 9.5.26. Here, we accomplish the inversion by + // performing a Taylor expansion of Eq. 9.3.39 for large z to order 2 + // and solving the resulting quadratic equation, thereby taking + // the positive root of the quadratic. + // In other words: (2/3)(-zeta)^(3/2) approx = z + 1/(2z) - pi/2. + // This leads to: z^2 - [(2/3)(-zeta)^(3/2) + pi/2]z + 1/2 = 0. + // + // With this initial estimate, Newton-Raphson iteration is used + // to refine the value of the estimate of the root of z + // as a function of zeta. + + const T v_pow_third(boost::math::cbrt(v)); + const T v_pow_minus_two_thirds(T(1) / (v_pow_third * v_pow_third)); + + // Obtain zeta using the order v combined with the m'th root of + // an airy function, as shown in A&S Eq. 9.5.22. + const T zeta = v_pow_minus_two_thirds * (-ai_bi_root); + + const T zeta_sqrt = sqrt(zeta); + + // Set up a quadratic equation based on the Taylor series + // expansion mentioned above. + const T b = -((((zeta * zeta_sqrt) * 2U) / 3U) + boost::math::constants::half_pi<T>()); + + // Solve the quadratic equation, taking the positive root. + const T z_estimate = (-b + sqrt((b * b) - T(2))) / 2U; + + // Establish the range, the digits, and the iteration limit + // for the upcoming root-finding. + const T range_zmin = (std::max<T>)(z_estimate - T(1), T(1)); + const T range_zmax = z_estimate + T(1); + + const int my_digits10 = static_cast<int>(static_cast<float>(boost::math::tools::digits<T>() * 0.301F)); + + // Select the maximum allowed iterations based on the number + // of decimal digits in the numeric type T, being at least 12. + const boost::uintmax_t iterations_allowed = static_cast<boost::uintmax_t>((std::max)(12, my_digits10 * 2)); + + boost::uintmax_t iterations_used = iterations_allowed; + + // Calculate the root of z as a function of zeta. + const T z = boost::math::tools::newton_raphson_iterate( + boost::math::detail::bessel_zero::equation_as_9_3_39_and_its_derivative<T>(zeta), + z_estimate, + range_zmin, + range_zmax, + (std::min)(boost::math::tools::digits<T>(), boost::math::tools::digits<float>()), + iterations_used); + + static_cast<void>(iterations_used); + + // Continue with the implementation of A&S Eq. 9.3.39. + const T zsq_minus_one = (z * z) - T(1); + const T zsq_minus_one_sqrt = sqrt(zsq_minus_one); + + // This is A&S Eq. 9.3.42. + const T b0_term_5_24 = T(5) / ((zsq_minus_one * zsq_minus_one_sqrt) * 24U); + const T b0_term_1_8 = T(1) / ( zsq_minus_one_sqrt * 8U); + const T b0_term_5_48 = T(5) / ((zeta * zeta) * 48U); + + const T b0 = -b0_term_5_48 + ((b0_term_5_24 + b0_term_1_8) / zeta_sqrt); + + // This is the second line of A&S Eq. 9.5.26 for f_k with k = 1. + const T f1 = ((z * zeta_sqrt) * b0) / zsq_minus_one_sqrt; + + // This is A&S Eq. 9.5.22 expanded to k = 1 (i.e., one term in the series). + return (v * z) + (f1 / v); + } + + namespace cyl_bessel_j_zero_detail + { + template<class T> + T equation_nist_10_21_40_a(const T& v) + { + const T v_pow_third(boost::math::cbrt(v)); + const T v_pow_minus_two_thirds(T(1) / (v_pow_third * v_pow_third)); + + return v * ((((( + T(0.043) + * v_pow_minus_two_thirds - T(0.0908)) + * v_pow_minus_two_thirds - T(0.00397)) + * v_pow_minus_two_thirds + T(1.033150)) + * v_pow_minus_two_thirds + T(1.8557571)) + * v_pow_minus_two_thirds + T(1)); + } + + template<class T, class Policy> + class function_object_jv + { + public: + function_object_jv(const T& v, + const Policy& pol) : my_v(v), + my_pol(pol) { } + + T operator()(const T& x) const + { + return boost::math::cyl_bessel_j(my_v, x, my_pol); + } + + private: + const T my_v; + const Policy& my_pol; + const function_object_jv& operator=(const function_object_jv&); + }; + + template<class T, class Policy> + class function_object_jv_and_jv_prime + { + public: + function_object_jv_and_jv_prime(const T& v, + const bool order_is_zero, + const Policy& pol) : my_v(v), + my_order_is_zero(order_is_zero), + my_pol(pol) { } + + boost::math::tuple<T, T> operator()(const T& x) const + { + // Obtain Jv(x) and Jv'(x). + // Chris's original code called the Bessel function implementation layer direct, + // but that circumvented optimizations for integer-orders. Call the documented + // top level functions instead, and let them sort out which implementation to use. + T j_v; + T j_v_prime; + + if(my_order_is_zero) + { + j_v = boost::math::cyl_bessel_j(0, x, my_pol); + j_v_prime = -boost::math::cyl_bessel_j(1, x, my_pol); + } + else + { + j_v = boost::math::cyl_bessel_j( my_v, x, my_pol); + const T j_v_m1 (boost::math::cyl_bessel_j(T(my_v - 1), x, my_pol)); + j_v_prime = j_v_m1 - ((my_v * j_v) / x); + } + + // Return a tuple containing both Jv(x) and Jv'(x). + return boost::math::make_tuple(j_v, j_v_prime); + } + + private: + const T my_v; + const bool my_order_is_zero; + const Policy& my_pol; + const function_object_jv_and_jv_prime& operator=(const function_object_jv_and_jv_prime&); + }; + + template<class T> bool my_bisection_unreachable_tolerance(const T&, const T&) { return false; } + + template<class T, class Policy> + T initial_guess(const T& v, const int m, const Policy& pol) + { + BOOST_MATH_STD_USING // ADL of std names, needed for floor. + + // Compute an estimate of the m'th root of cyl_bessel_j. + + T guess; + + // There is special handling for negative order. + if(v < 0) + { + if((m == 1) && (v > -0.5F)) + { + // For small, negative v, use the results of empirical curve fitting. + // Mathematica(R) session for the coefficients: + // Table[{n, BesselJZero[n, 1]}, {n, -(1/2), 0, 1/10}] + // N[%, 20] + // Fit[%, {n^0, n^1, n^2, n^3, n^4, n^5, n^6}, n] + guess = ((((( - T(0.2321156900729) + * v - T(0.1493247777488)) + * v - T(0.15205419167239)) + * v + T(0.07814930561249)) + * v - T(0.17757573537688)) + * v + T(1.542805677045663)) + * v + T(2.40482555769577277); + + return guess; + } + + // Create the positive order and extract its positive floor integer part. + const T vv(-v); + const T vv_floor(floor(vv)); + + // The to-be-found root is bracketed by the roots of the + // Bessel function whose reflected, positive integer order + // is less than, but nearest to vv. + + T root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(vv_floor, m, pol); + T root_lo; + + if(m == 1) + { + // The estimate of the first root for negative order is found using + // an adaptive range-searching algorithm. + root_lo = T(root_hi - 0.1F); + + const bool hi_end_of_bracket_is_negative = (boost::math::cyl_bessel_j(v, root_hi, pol) < 0); + + while((root_lo > boost::math::tools::epsilon<T>())) + { + const bool lo_end_of_bracket_is_negative = (boost::math::cyl_bessel_j(v, root_lo, pol) < 0); + + if(hi_end_of_bracket_is_negative != lo_end_of_bracket_is_negative) + { + break; + } + + root_hi = root_lo; + + // Decrease the lower end of the bracket using an adaptive algorithm. + if(root_lo > 0.5F) + { + root_lo -= 0.5F; + } + else + { + root_lo *= 0.75F; + } + } + } + else + { + root_lo = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(vv_floor, m - 1, pol); + } + + // Perform several steps of bisection iteration to refine the guess. + boost::uintmax_t number_of_iterations(12U); + + // Do the bisection iteration. + const boost::math::tuple<T, T> guess_pair = + boost::math::tools::bisect( + boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::function_object_jv<T, Policy>(v, pol), + root_lo, + root_hi, + boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::my_bisection_unreachable_tolerance<T>, + number_of_iterations); + + return (boost::math::get<0>(guess_pair) + boost::math::get<1>(guess_pair)) / 2U; + } + + if(m == 1U) + { + // Get the initial estimate of the first root. + + if(v < 2.2F) + { + // For small v, use the results of empirical curve fitting. + // Mathematica(R) session for the coefficients: + // Table[{n, BesselJZero[n, 1]}, {n, 0, 22/10, 1/10}] + // N[%, 20] + // Fit[%, {n^0, n^1, n^2, n^3, n^4, n^5, n^6}, n] + guess = ((((( - T(0.0008342379046010) + * v + T(0.007590035637410)) + * v - T(0.030640914772013)) + * v + T(0.078232088020106)) + * v - T(0.169668712590620)) + * v + T(1.542187960073750)) + * v + T(2.4048359915254634); + } + else + { + // For larger v, use the first line of Eqs. 10.21.40 in the NIST Handbook. + guess = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::equation_nist_10_21_40_a(v); + } + } + else + { + if(v < 2.2F) + { + // Use Eq. 10.21.19 in the NIST Handbook. + const T a(((v + T(m * 2U)) - T(0.5)) * boost::math::constants::half_pi<T>()); + + guess = boost::math::detail::bessel_zero::equation_nist_10_21_19(v, a); + } + else + { + // Get an estimate of the m'th root of airy_ai. + const T airy_ai_root(boost::math::detail::airy_zero::airy_ai_zero_detail::initial_guess<T>(m)); + + // Use Eq. 9.5.26 in the A&S Handbook. + guess = boost::math::detail::bessel_zero::equation_as_9_5_26(v, airy_ai_root); + } + } + + return guess; + } + } // namespace cyl_bessel_j_zero_detail + + namespace cyl_neumann_zero_detail + { + template<class T> + T equation_nist_10_21_40_b(const T& v) + { + const T v_pow_third(boost::math::cbrt(v)); + const T v_pow_minus_two_thirds(T(1) / (v_pow_third * v_pow_third)); + + return v * ((((( - T(0.001) + * v_pow_minus_two_thirds - T(0.0060)) + * v_pow_minus_two_thirds + T(0.01198)) + * v_pow_minus_two_thirds + T(0.260351)) + * v_pow_minus_two_thirds + T(0.9315768)) + * v_pow_minus_two_thirds + T(1)); + } + + template<class T, class Policy> + class function_object_yv + { + public: + function_object_yv(const T& v, + const Policy& pol) : my_v(v), + my_pol(pol) { } + + T operator()(const T& x) const + { + return boost::math::cyl_neumann(my_v, x, my_pol); + } + + private: + const T my_v; + const Policy& my_pol; + const function_object_yv& operator=(const function_object_yv&); + }; + + template<class T, class Policy> + class function_object_yv_and_yv_prime + { + public: + function_object_yv_and_yv_prime(const T& v, + const Policy& pol) : my_v(v), + my_pol(pol) { } + + boost::math::tuple<T, T> operator()(const T& x) const + { + const T half_epsilon(boost::math::tools::epsilon<T>() / 2U); + + const bool order_is_zero = ((my_v > -half_epsilon) && (my_v < +half_epsilon)); + + // Obtain Yv(x) and Yv'(x). + // Chris's original code called the Bessel function implementation layer direct, + // but that circumvented optimizations for integer-orders. Call the documented + // top level functions instead, and let them sort out which implementation to use. + T y_v; + T y_v_prime; + + if(order_is_zero) + { + y_v = boost::math::cyl_neumann(0, x, my_pol); + y_v_prime = -boost::math::cyl_neumann(1, x, my_pol); + } + else + { + y_v = boost::math::cyl_neumann( my_v, x, my_pol); + const T y_v_m1 (boost::math::cyl_neumann(T(my_v - 1), x, my_pol)); + y_v_prime = y_v_m1 - ((my_v * y_v) / x); + } + + // Return a tuple containing both Yv(x) and Yv'(x). + return boost::math::make_tuple(y_v, y_v_prime); + } + + private: + const T my_v; + const Policy& my_pol; + const function_object_yv_and_yv_prime& operator=(const function_object_yv_and_yv_prime&); + }; + + template<class T> bool my_bisection_unreachable_tolerance(const T&, const T&) { return false; } + + template<class T, class Policy> + T initial_guess(const T& v, const int m, const Policy& pol) + { + BOOST_MATH_STD_USING // ADL of std names, needed for floor. + + // Compute an estimate of the m'th root of cyl_neumann. + + T guess; + + // There is special handling for negative order. + if(v < 0) + { + // Create the positive order and extract its positive floor and ceiling integer parts. + const T vv(-v); + const T vv_floor(floor(vv)); + + // The to-be-found root is bracketed by the roots of the + // Bessel function whose reflected, positive integer order + // is less than, but nearest to vv. + + // The special case of negative, half-integer order uses + // the relation between Yv and spherical Bessel functions + // in order to obtain the bracket for the root. + // In these special cases, cyl_neumann(-n/2, x) = sph_bessel_j(+n/2, x) + // for v = -n/2. + + T root_hi; + T root_lo; + + if(m == 1) + { + // The estimate of the first root for negative order is found using + // an adaptive range-searching algorithm. + // Take special precautions for the discontinuity at negative, + // half-integer orders and use different brackets above and below these. + if(T(vv - vv_floor) < 0.5F) + { + root_hi = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m, pol); + } + else + { + root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m, pol); + } + + root_lo = T(root_hi - 0.1F); + + const bool hi_end_of_bracket_is_negative = (boost::math::cyl_neumann(v, root_hi, pol) < 0); + + while((root_lo > boost::math::tools::epsilon<T>())) + { + const bool lo_end_of_bracket_is_negative = (boost::math::cyl_neumann(v, root_lo, pol) < 0); + + if(hi_end_of_bracket_is_negative != lo_end_of_bracket_is_negative) + { + break; + } + + root_hi = root_lo; + + // Decrease the lower end of the bracket using an adaptive algorithm. + if(root_lo > 0.5F) + { + root_lo -= 0.5F; + } + else + { + root_lo *= 0.75F; + } + } + } + else + { + if(T(vv - vv_floor) < 0.5F) + { + root_lo = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m - 1, pol); + root_hi = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m, pol); + root_lo += 0.01F; + root_hi += 0.01F; + } + else + { + root_lo = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m - 1, pol); + root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m, pol); + root_lo += 0.01F; + root_hi += 0.01F; + } + } + + // Perform several steps of bisection iteration to refine the guess. + boost::uintmax_t number_of_iterations(12U); + + // Do the bisection iteration. + const boost::math::tuple<T, T> guess_pair = + boost::math::tools::bisect( + boost::math::detail::bessel_zero::cyl_neumann_zero_detail::function_object_yv<T, Policy>(v, pol), + root_lo, + root_hi, + boost::math::detail::bessel_zero::cyl_neumann_zero_detail::my_bisection_unreachable_tolerance<T>, + number_of_iterations); + + return (boost::math::get<0>(guess_pair) + boost::math::get<1>(guess_pair)) / 2U; + } + + if(m == 1U) + { + // Get the initial estimate of the first root. + + if(v < 2.2F) + { + // For small v, use the results of empirical curve fitting. + // Mathematica(R) session for the coefficients: + // Table[{n, BesselYZero[n, 1]}, {n, 0, 22/10, 1/10}] + // N[%, 20] + // Fit[%, {n^0, n^1, n^2, n^3, n^4, n^5, n^6}, n] + guess = ((((( - T(0.0025095909235652) + * v + T(0.021291887049053)) + * v - T(0.076487785486526)) + * v + T(0.159110268115362)) + * v - T(0.241681668765196)) + * v + T(1.4437846310885244)) + * v + T(0.89362115190200490); + } + else + { + // For larger v, use the second line of Eqs. 10.21.40 in the NIST Handbook. + guess = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::equation_nist_10_21_40_b(v); + } + } + else + { + if(v < 2.2F) + { + // Use Eq. 10.21.19 in the NIST Handbook. + const T a(((v + T(m * 2U)) - T(1.5)) * boost::math::constants::half_pi<T>()); + + guess = boost::math::detail::bessel_zero::equation_nist_10_21_19(v, a); + } + else + { + // Get an estimate of the m'th root of airy_bi. + const T airy_bi_root(boost::math::detail::airy_zero::airy_bi_zero_detail::initial_guess<T>(m)); + + // Use Eq. 9.5.26 in the A&S Handbook. + guess = boost::math::detail::bessel_zero::equation_as_9_5_26(v, airy_bi_root); + } + } + + return guess; + } + } // namespace cyl_neumann_zero_detail + } // namespace bessel_zero + } } } // namespace boost::math::detail + +#endif // _BESSEL_JY_ZERO_2013_01_18_HPP_ |