123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627 |
- #ifndef BOOST_MATH_BESSEL_JY_ZERO_2013_01_18_HPP_
- #define BOOST_MATH_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)
- {
-
-
-
-
- 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:
- explicit equation_as_9_3_39_and_its_derivative(const T& zt) : zeta(zt) { }
- equation_as_9_3_39_and_its_derivative(const equation_as_9_3_39_and_its_derivative&) = default;
- boost::math::tuple<T, T> operator()(const T& z) const
- {
- BOOST_MATH_STD_USING
-
-
-
- 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&) = delete;
- const T zeta;
- };
- template<class T, class Policy>
- static T equation_as_9_5_26(const T& v, const T& ai_bi_root, const Policy& pol)
- {
- BOOST_MATH_STD_USING
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- const T v_pow_third(boost::math::cbrt(v, pol));
- const T v_pow_minus_two_thirds(T(1) / (v_pow_third * v_pow_third));
-
-
- const T zeta = v_pow_minus_two_thirds * (-ai_bi_root);
- const T zeta_sqrt = sqrt(zeta);
-
-
- const T b = -((((zeta * zeta_sqrt) * 2U) / 3U) + boost::math::constants::half_pi<T>());
-
- const T z_estimate = (-b + sqrt((b * b) - T(2))) / 2U;
-
-
- const T range_zmin = (std::max<T>)(z_estimate - T(1), T(1));
- const T range_zmax = z_estimate + T(1);
- const auto my_digits10 = static_cast<int>(static_cast<float>(boost::math::tools::digits<T>() * 0.301F));
-
-
- const auto iterations_allowed = static_cast<std::uintmax_t>((std::max)(12, my_digits10 * 2));
- std::uintmax_t iterations_used = iterations_allowed;
-
- 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);
-
- const T zsq_minus_one = (z * z) - T(1);
- const T zsq_minus_one_sqrt = sqrt(zsq_minus_one);
-
- 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);
-
- const T f1 = ((z * zeta_sqrt) * b0) / zsq_minus_one_sqrt;
-
- return (v * z) + (f1 / v);
- }
- namespace cyl_bessel_j_zero_detail
- {
- template<class T, class Policy>
- T equation_nist_10_21_40_a(const T& v, const Policy& pol)
- {
- const T v_pow_third(boost::math::cbrt(v, pol));
- 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) { }
- function_object_jv(const function_object_jv&) = default;
- 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&) = delete;
- };
- 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) { }
- function_object_jv_and_jv_prime(const function_object_jv_and_jv_prime&) = default;
- boost::math::tuple<T, T> operator()(const T& x) const
- {
-
-
-
-
- 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 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&) = delete;
- };
- 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
-
- T guess;
-
- if(v < 0)
- {
- if((m == 1) && (v > -0.5F))
- {
-
-
-
-
-
- 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;
- }
-
- const T vv(-v);
- const T vv_floor(floor(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)
- {
-
-
- 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;
-
- 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);
- }
-
- std::uintmax_t number_of_iterations(12U);
-
- 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)
- {
-
- if(v < 2.2F)
- {
-
-
-
-
-
- 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
- {
-
- guess = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::equation_nist_10_21_40_a(v, pol);
- }
- }
- else
- {
- if(v < 2.2F)
- {
-
- 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
- {
-
- const T airy_ai_root(boost::math::detail::airy_zero::airy_ai_zero_detail::initial_guess<T>(m, pol));
-
- guess = boost::math::detail::bessel_zero::equation_as_9_5_26(v, airy_ai_root, pol);
- }
- }
- return guess;
- }
- }
- namespace cyl_neumann_zero_detail
- {
- template<class T, class Policy>
- T equation_nist_10_21_40_b(const T& v, const Policy& pol)
- {
- const T v_pow_third(boost::math::cbrt(v, pol));
- 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) { }
- function_object_yv(const function_object_yv&) = default;
- 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&) = delete;
- };
- 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) { }
- function_object_yv_and_yv_prime(const function_object_yv_and_yv_prime&) = default;
- 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));
-
-
-
-
- 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 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&) = delete;
- };
- 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
-
- T guess;
-
- if(v < 0)
- {
-
- const T vv(-v);
- const T vv_floor(floor(vv));
-
-
-
-
-
-
-
-
- T root_hi;
- T root_lo;
- if(m == 1)
- {
-
-
-
-
- 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;
-
- 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;
- }
- }
-
- std::uintmax_t number_of_iterations(12U);
-
- 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)
- {
-
- if(v < 2.2F)
- {
-
-
-
-
-
- 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
- {
-
- guess = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::equation_nist_10_21_40_b(v, pol);
- }
- }
- else
- {
- if(v < 2.2F)
- {
-
- 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
- {
-
- const T airy_bi_root(boost::math::detail::airy_zero::airy_bi_zero_detail::initial_guess<T>(m, pol));
-
- guess = boost::math::detail::bessel_zero::equation_as_9_5_26(v, airy_bi_root, pol);
- }
- }
- return guess;
- }
- }
- }
- } } }
- #endif
|