From 0e2463bfd787ca91d5e293e581ea3aaa5e81932e Mon Sep 17 00:00:00 2001 From: John Maddock Date: Sat, 9 Feb 2013 11:56:52 +0000 Subject: [PATCH] Tidy up policy usage and error handling in Bessel functions. Change zero-finder functors to call top level Bessel functions. [SVN r82792] --- .../boost/math/special_functions/bessel.hpp | 88 +++++++++++++++---- .../detail/bessel_jy_zero.hpp | 18 ++-- 2 files changed, 82 insertions(+), 24 deletions(-) diff --git a/include/boost/math/special_functions/bessel.hpp b/include/boost/math/special_functions/bessel.hpp index 78b5e4ad2..284afa63d 100644 --- a/include/boost/math/special_functions/bessel.hpp +++ b/include/boost/math/special_functions/bessel.hpp @@ -386,13 +386,11 @@ inline T cyl_bessel_j_zero_imp(T v, unsigned m, const Policy& pol) const T guess_root = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(v, m); // Select the maximum allowed iterations, being at least 24. - boost::uintmax_t number_of_iterations = (std::max)(24, int(std::numeric_limits::digits10)); + boost::uintmax_t number_of_iterations = policies::get_max_root_iterations(); // Select the desired number of binary digits of precision. // Account for the radix of number representations having non-two radix! - const int my_digits2 = int(float(std::numeric_limits::digits) - * ( log(float(std::numeric_limits::radix)) - / log(2.0F))); + const int my_digits2 = policies::digits(); // Perform the root-finding using Newton-Raphson iteration from Boost.Math. const T jvm = @@ -404,7 +402,11 @@ inline T cyl_bessel_j_zero_imp(T v, unsigned m, const Policy& pol) my_digits2, number_of_iterations); - static_cast(number_of_iterations); + if(number_of_iterations >= policies::get_max_root_iterations()) + { + policies::raise_evaluation_error(function, "Unable to locate root in a reasonable time:" + " Current best guess is %1%", jvm, Policy()); + } return jvm; } @@ -425,13 +427,11 @@ inline T cyl_neumann_zero_imp(T v, unsigned m, const Policy& pol) const T guess_root = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(v, m); // Select the maximum allowed iterations, being at least 24. - boost::uintmax_t number_of_iterations = (std::max)(24, int(std::numeric_limits::digits10)); + boost::uintmax_t number_of_iterations = policies::get_max_root_iterations(); // Select the desired number of binary digits of precision. // Account for the radix of number representations having non-two radix! - const int my_digits2 = int(float(std::numeric_limits::digits) - * ( log(float(std::numeric_limits::radix)) - / log(2.0F))); + const int my_digits2 = policies::digits(); // Perform the root-finding using Newton-Raphson iteration from Boost.Math. const T yvm = @@ -443,7 +443,11 @@ inline T cyl_neumann_zero_imp(T v, unsigned m, const Policy& pol) my_digits2, number_of_iterations); - static_cast(number_of_iterations); + if(number_of_iterations >= policies::get_max_root_iterations()) + { + policies::raise_evaluation_error(function, "Unable to locate root in a reasonable time:" + " Current best guess is %1%", yvm, Policy()); + } return yvm; } @@ -457,7 +461,13 @@ inline typename detail::bessel_traits::result_type cyl_bessel_j( typedef typename detail::bessel_traits::result_type result_type; typedef typename detail::bessel_traits::optimisation_tag tag_type; typedef typename policies::evaluation::type value_type; - return policies::checked_narrowing_cast(detail::cyl_bessel_j_imp(v, static_cast(x), tag_type(), pol), "boost::math::cyl_bessel_j<%1%>(%1%,%1%)"); + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + return policies::checked_narrowing_cast(detail::cyl_bessel_j_imp(v, static_cast(x), tag_type(), forwarding_policy()), "boost::math::cyl_bessel_j<%1%>(%1%,%1%)"); } template @@ -472,7 +482,13 @@ inline typename detail::bessel_traits::result_type sph_bessel(unsi BOOST_FPU_EXCEPTION_GUARD typedef typename detail::bessel_traits::result_type result_type; typedef typename policies::evaluation::type value_type; - return policies::checked_narrowing_cast(detail::sph_bessel_j_imp(v, static_cast(x), pol), "boost::math::sph_bessel<%1%>(%1%,%1%)"); + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + return policies::checked_narrowing_cast(detail::sph_bessel_j_imp(v, static_cast(x), forwarding_policy()), "boost::math::sph_bessel<%1%>(%1%,%1%)"); } template @@ -487,7 +503,13 @@ inline typename detail::bessel_traits::result_type cyl_bessel_i( BOOST_FPU_EXCEPTION_GUARD typedef typename detail::bessel_traits::result_type result_type; typedef typename policies::evaluation::type value_type; - return policies::checked_narrowing_cast(detail::cyl_bessel_i_imp(v, static_cast(x), pol), "boost::math::cyl_bessel_i<%1%>(%1%,%1%)"); + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + return policies::checked_narrowing_cast(detail::cyl_bessel_i_imp(v, static_cast(x), forwarding_policy()), "boost::math::cyl_bessel_i<%1%>(%1%,%1%)"); } template @@ -503,7 +525,13 @@ inline typename detail::bessel_traits::result_type cyl_bessel_k( typedef typename detail::bessel_traits::result_type result_type; typedef typename detail::bessel_traits::optimisation_tag tag_type; typedef typename policies::evaluation::type value_type; - return policies::checked_narrowing_cast(detail::cyl_bessel_k_imp(v, static_cast(x), tag_type(), pol), "boost::math::cyl_bessel_k<%1%>(%1%,%1%)"); + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + return policies::checked_narrowing_cast(detail::cyl_bessel_k_imp(v, static_cast(x), tag_type(), forwarding_policy()), "boost::math::cyl_bessel_k<%1%>(%1%,%1%)"); } template @@ -519,7 +547,13 @@ inline typename detail::bessel_traits::result_type cyl_neumann(T typedef typename detail::bessel_traits::result_type result_type; typedef typename detail::bessel_traits::optimisation_tag tag_type; typedef typename policies::evaluation::type value_type; - return policies::checked_narrowing_cast(detail::cyl_neumann_imp(v, static_cast(x), tag_type(), pol), "boost::math::cyl_neumann<%1%>(%1%,%1%)"); + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + return policies::checked_narrowing_cast(detail::cyl_neumann_imp(v, static_cast(x), tag_type(), forwarding_policy()), "boost::math::cyl_neumann<%1%>(%1%,%1%)"); } template @@ -534,7 +568,13 @@ inline typename detail::bessel_traits::result_type sph_neumann(uns BOOST_FPU_EXCEPTION_GUARD typedef typename detail::bessel_traits::result_type result_type; typedef typename policies::evaluation::type value_type; - return policies::checked_narrowing_cast(detail::sph_neumann_imp(v, static_cast(x), pol), "boost::math::sph_neumann<%1%>(%1%,%1%)"); + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + return policies::checked_narrowing_cast(detail::sph_neumann_imp(v, static_cast(x), forwarding_policy()), "boost::math::sph_neumann<%1%>(%1%,%1%)"); } template @@ -549,8 +589,14 @@ inline typename detail::bessel_traits::result_type cyl_bessel_j_ze BOOST_FPU_EXCEPTION_GUARD typedef typename detail::bessel_traits::result_type result_type; typedef typename policies::evaluation::type value_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits::is_integer, "Order must be a floating-point type."); - return policies::checked_narrowing_cast(detail::cyl_bessel_j_zero_imp(v, m, pol), "boost::math::cyl_bessel_j_zero<%1%>(%1%,%1%)"); + return policies::checked_narrowing_cast(detail::cyl_bessel_j_zero_imp(v, m, forwarding_policy()), "boost::math::cyl_bessel_j_zero<%1%>(%1%,%1%)"); } template @@ -591,8 +637,14 @@ inline typename detail::bessel_traits::result_type cyl_neumann_zer BOOST_FPU_EXCEPTION_GUARD typedef typename detail::bessel_traits::result_type result_type; typedef typename policies::evaluation::type value_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits::is_integer, "Order must be a floating-point type."); - return policies::checked_narrowing_cast(detail::cyl_neumann_zero_imp(v, m, pol), "boost::math::cyl_neumann_zero<%1%>(%1%,%1%)"); + return policies::checked_narrowing_cast(detail::cyl_neumann_zero_imp(v, m, forwarding_policy()), "boost::math::cyl_neumann_zero<%1%>(%1%,%1%)"); } template diff --git a/include/boost/math/special_functions/detail/bessel_jy_zero.hpp b/include/boost/math/special_functions/detail/bessel_jy_zero.hpp index 05693a725..d0363e363 100644 --- a/include/boost/math/special_functions/detail/bessel_jy_zero.hpp +++ b/include/boost/math/special_functions/detail/bessel_jy_zero.hpp @@ -257,18 +257,21 @@ && (my_v < +boost::math::tools::epsilon())); // 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(order_is_zero) { - j_v = boost::math::detail::cyl_bessel_j_imp(T(0), x, boost::math::detail::bessel_no_int_tag(), my_pol); - j_v_prime = -boost::math::detail::cyl_bessel_j_imp(T(1), x, boost::math::detail::bessel_no_int_tag(), my_pol); + 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::detail::cyl_bessel_j_imp( my_v, x, boost::math::detail::bessel_no_int_tag(), my_pol); - const T j_v_m1 (boost::math::detail::cyl_bessel_j_imp(T(my_v - 1), x, boost::math::detail::bessel_no_int_tag(), my_pol)); + 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); } @@ -370,8 +373,11 @@ boost::math::tuple operator()(const T& x) const { // Obtain Yv(x) and Yv'(x). - const T y_v (boost::math::detail::cyl_neumann_imp( my_v, x, boost::math::detail::bessel_no_int_tag(), my_pol)); - const T y_v_m1 (boost::math::detail::cyl_neumann_imp(T(my_v - 1), x, boost::math::detail::bessel_no_int_tag(), my_pol)); + // 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. + const T 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)); const T y_v_prime(y_v_m1 - ((my_v * y_v) / x)); // Return a tuple containing both Yv(x) and Yv'(x).