boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-28 07:22:12 +00:00
Code Issues Projects Releases Wiki Activity

Tidy up policy usage and error handling in Bessel functions.

Browse Source 
Change zero-finder functors to call top level Bessel functions.

[SVN r82792]
This commit is contained in:
John Maddock
2013-02-09 11:56:52 +00:00
parent 4d1e810132
commit 0e2463bfd7
2 changed files with 82 additions and 24 deletions
Download Patch File Download Diff File Expand all files Collapse all files

88
include/boost/math/special_functions/bessel.hpp
View File

@@ -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<T>(v, m);
// Select the maximum allowed iterations, being at least 24.
boost::uintmax_t number_of_iterations = (std::max)(24, int(std::numeric_limits<T>::digits10));
boost::uintmax_t number_of_iterations = policies::get_max_root_iterations<Policy>();
// 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<T>::digits)
* ( log(float(std::numeric_limits<T>::radix))
/ log(2.0F)));
const int my_digits2 = policies::digits<T, Policy>();
// 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<void>(number_of_iterations);
if(number_of_iterations >= policies::get_max_root_iterations<Policy>())
{
policies::raise_evaluation_error<T>(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<T>(v, m);
// Select the maximum allowed iterations, being at least 24.
boost::uintmax_t number_of_iterations = (std::max)(24, int(std::numeric_limits<T>::digits10));
boost::uintmax_t number_of_iterations = policies::get_max_root_iterations<Policy>();
// 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<T>::digits)
* ( log(float(std::numeric_limits<T>::radix))
/ log(2.0F)));
const int my_digits2 = policies::digits<T, Policy>();
// 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<void>(number_of_iterations);
if(number_of_iterations >= policies::get_max_root_iterations<Policy>())
{
policies::raise_evaluation_error<T>(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<T1, T2, Policy>::result_type cyl_bessel_j(
typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_j_imp<value_type>(v, static_cast<value_type>(x), tag_type(), pol), "boost::math::cyl_bessel_j<%1%>(%1%,%1%)");
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_j_imp<value_type>(v, static_cast<value_type>(x), tag_type(), forwarding_policy()), "boost::math::cyl_bessel_j<%1%>(%1%,%1%)");
}
template <class T1, class T2>
@@ -472,7 +482,13 @@ inline typename detail::bessel_traits<T, T, Policy>::result_type sph_bessel(unsi
BOOST_FPU_EXCEPTION_GUARD
typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
return policies::checked_narrowing_cast<result_type, Policy>(detail::sph_bessel_j_imp<value_type>(v, static_cast<value_type>(x), pol), "boost::math::sph_bessel<%1%>(%1%,%1%)");
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
return policies::checked_narrowing_cast<result_type, Policy>(detail::sph_bessel_j_imp<value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::sph_bessel<%1%>(%1%,%1%)");
}
template <class T>
@@ -487,7 +503,13 @@ inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_i(
BOOST_FPU_EXCEPTION_GUARD
typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_i_imp<value_type>(v, static_cast<value_type>(x), pol), "boost::math::cyl_bessel_i<%1%>(%1%,%1%)");
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_i_imp<value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_i<%1%>(%1%,%1%)");
}
template <class T1, class T2>
@@ -503,7 +525,13 @@ inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_k(
typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_k_imp<value_type>(v, static_cast<value_type>(x), tag_type(), pol), "boost::math::cyl_bessel_k<%1%>(%1%,%1%)");
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_k_imp<value_type>(v, static_cast<value_type>(x), tag_type(), forwarding_policy()), "boost::math::cyl_bessel_k<%1%>(%1%,%1%)");
}
template <class T1, class T2>
@@ -519,7 +547,13 @@ inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_neumann(T
typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_neumann_imp<value_type>(v, static_cast<value_type>(x), tag_type(), pol), "boost::math::cyl_neumann<%1%>(%1%,%1%)");
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_neumann_imp<value_type>(v, static_cast<value_type>(x), tag_type(), forwarding_policy()), "boost::math::cyl_neumann<%1%>(%1%,%1%)");
}
template <class T1, class T2>
@@ -534,7 +568,13 @@ inline typename detail::bessel_traits<T, T, Policy>::result_type sph_neumann(uns
BOOST_FPU_EXCEPTION_GUARD
typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
return policies::checked_narrowing_cast<result_type, Policy>(detail::sph_neumann_imp<value_type>(v, static_cast<value_type>(x), pol), "boost::math::sph_neumann<%1%>(%1%,%1%)");
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
return policies::checked_narrowing_cast<result_type, Policy>(detail::sph_neumann_imp<value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::sph_neumann<%1%>(%1%,%1%)");
}
template <class T>
@@ -549,8 +589,14 @@ inline typename detail::bessel_traits<T, T, Policy>::result_type cyl_bessel_j_ze
BOOST_FPU_EXCEPTION_GUARD
typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<value_type>::is_integer, "Order must be a floating-point type.");
return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_j_zero_imp<value_type>(v, m, pol), "boost::math::cyl_bessel_j_zero<%1%>(%1%,%1%)");
return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_j_zero_imp<value_type>(v, m, forwarding_policy()), "boost::math::cyl_bessel_j_zero<%1%>(%1%,%1%)");
}
template <class T>
@@ -591,8 +637,14 @@ inline typename detail::bessel_traits<T, T, Policy>::result_type cyl_neumann_zer
BOOST_FPU_EXCEPTION_GUARD
typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<value_type>::is_integer, "Order must be a floating-point type.");
return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_neumann_zero_imp<value_type>(v, m, pol), "boost::math::cyl_neumann_zero<%1%>(%1%,%1%)");
return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_neumann_zero_imp<value_type>(v, m, forwarding_policy()), "boost::math::cyl_neumann_zero<%1%>(%1%,%1%)");
}
template <class T>

18
include/boost/math/special_functions/detail/bessel_jy_zero.hpp
View File

@@ -257,18 +257,21 @@
&& (my_v < +boost::math::tools::epsilon<T>()));
// 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<T, T> 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).