|
|
Boost-Commit : |
From: johnmaddock_at_[hidden]
Date: 2007-07-13 08:49:00
Author: johnmaddock
Date: 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
New Revision: 7420
URL: http://svn.boost.org/trac/boost/changeset/7420
Log:
Added policy framework to the bessel functions, plus a few others.
Removed:
sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_jv.hpp
sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_kv.hpp
sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_yv.hpp
Text files modified:
sandbox/math_toolkit/policy/boost/math/policy/error_handling.hpp | 6
sandbox/math_toolkit/policy/boost/math/special_functions/acosh.hpp | 21 +
sandbox/math_toolkit/policy/boost/math/special_functions/asinh.hpp | 8
sandbox/math_toolkit/policy/boost/math/special_functions/atanh.hpp | 35 ++-
sandbox/math_toolkit/policy/boost/math/special_functions/bessel.hpp | 292 ++++++++++++++++++++++-----------------
sandbox/math_toolkit/policy/boost/math/special_functions/cos_pi.hpp | 6
sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_ik.hpp | 48 +++---
sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_jn.hpp | 6
sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_jy.hpp | 52 +++---
sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_jy_asym.hpp | 36 ++--
sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_k0.hpp | 14 +
sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_k1.hpp | 14 +
sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_kn.hpp | 22 +-
sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_y0.hpp | 14 +
sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_y1.hpp | 10
sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_yn.hpp | 22 +-
sandbox/math_toolkit/policy/boost/math/special_functions/factorials.hpp | 73 +++++++--
sandbox/math_toolkit/policy/boost/math/special_functions/hypot.hpp | 18 +
sandbox/math_toolkit/policy/boost/math/special_functions/sin_pi.hpp | 6
sandbox/math_toolkit/policy/boost/math/special_functions/sinc.hpp | 14 +
sandbox/math_toolkit/policy/boost/math/tools/test.hpp | 16 +
sandbox/math_toolkit/policy/libs/math/test/hypot_test.cpp | 1
sandbox/math_toolkit/policy/libs/math/test/test_bessel_i.cpp | 1
sandbox/math_toolkit/policy/libs/math/test/test_bessel_k.cpp | 1
sandbox/math_toolkit/policy/libs/math/test/test_bessel_y.cpp | 1
25 files changed, 449 insertions(+), 288 deletions(-)
Modified: sandbox/math_toolkit/policy/boost/math/policy/error_handling.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/policy/error_handling.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/policy/error_handling.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -181,7 +181,7 @@
const char* message,
const ::boost::math::policy::overflow_error< ::boost::math::policy::throw_on_error>&)
{
- raise_error<std::overflow_error, T>(function, message);
+ raise_error<std::overflow_error, T>(function, message ? message : "numeric overflow");
// we never get here:
return std::numeric_limits<T>::has_infinity ? std::numeric_limits<T>::infinity() : boost::math::tools::max_value<T>();
}
@@ -224,7 +224,7 @@
const char* message,
const ::boost::math::policy::underflow_error< ::boost::math::policy::throw_on_error>&)
{
- raise_error<std::underflow_error, T>(function, message);
+ raise_error<std::underflow_error, T>(function, message ? message : "numeric underflow");
// we never get here:
return 0;
}
@@ -268,7 +268,7 @@
const T& val,
const ::boost::math::policy::denorm_error< ::boost::math::policy::throw_on_error>&)
{
- raise_error<std::underflow_error, T>(function, message);
+ raise_error<std::underflow_error, T>(function, message ? message : "denormalised result");
// we never get here:
return 0;
}
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/acosh.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/acosh.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/acosh.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -14,7 +14,7 @@
#include <cmath>
#include <boost/config.hpp>
#include <boost/math/tools/precision.hpp>
-#include <boost/math/tools/error_handling.hpp>
+#include <boost/math/policy/error_handling.hpp>
#include <boost/math/special_functions/math_fwd.hpp>
// This is the inverse of the hyperbolic cosine function.
@@ -36,8 +36,8 @@
using ::std::numeric_limits;
#endif
- template<typename T>
- inline T acosh_imp(const T x)
+ template<typename T, typename Policy>
+ inline T acosh_imp(const T x, const Policy& pol)
{
using ::std::abs;
using ::std::sqrt;
@@ -52,9 +52,9 @@
if(x < one)
{
- return tools::domain_error<T>(
- BOOST_CURRENT_FUNCTION,
- "acosh requires x >= 1, but got x = %1%.", x);
+ return policy::raise_domain_error<T>(
+ "boost::math::acosh<%1%>(%1%)",
+ "acosh requires x >= 1, but got x = %1%.", x, pol);
}
else if (x >= taylor_n_bound)
{
@@ -88,12 +88,19 @@
}
}
+ template<typename T, typename Policy>
+ inline typename tools::promote_args<T>::type acosh(const T x, const Policy& pol)
+ {
+ typedef typename tools::promote_args<T>::type result_type;
+ return detail::acosh_imp(
+ static_cast<result_type>(x), pol);
+ }
template<typename T>
inline typename tools::promote_args<T>::type acosh(const T x)
{
typedef typename tools::promote_args<T>::type result_type;
return detail::acosh_imp(
- static_cast<result_type>(x));
+ static_cast<result_type>(x), policy::policy<>());
}
}
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/asinh.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/asinh.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/asinh.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -14,7 +14,6 @@
#include <cmath>
#include <boost/config.hpp>
#include <boost/math/tools/precision.hpp>
-#include <boost/math/tools/error_handling.hpp>
#include <boost/math/special_functions/math_fwd.hpp>
// This is the inverse of the hyperbolic sine function.
@@ -99,6 +98,13 @@
return detail::asinh_imp(
static_cast<result_type>(x));
}
+ template<typename T, typename Policy>
+ inline typename tools::promote_args<T>::type asinh(const T x, const Policy&)
+ {
+ typedef typename tools::promote_args<T>::type result_type;
+ return detail::asinh_imp(
+ static_cast<result_type>(x));
+ }
}
}
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/atanh.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/atanh.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/atanh.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -14,7 +14,7 @@
#include <cmath>
#include <boost/config.hpp>
#include <boost/math/tools/precision.hpp>
-#include <boost/math/tools/error_handling.hpp>
+#include <boost/math/policy/error_handling.hpp>
#include <boost/math/special_functions/math_fwd.hpp>
// This is the inverse of the hyperbolic tangent function.
@@ -38,8 +38,8 @@
// This is the main fare
- template<typename T>
- inline T atanh_imp(const T x)
+ template<typename T, typename Policy>
+ inline T atanh_imp(const T x, const Policy& pol)
{
using ::std::abs;
using ::std::sqrt;
@@ -52,30 +52,30 @@
static T const taylor_2_bound = sqrt(tools::epsilon<T>());
static T const taylor_n_bound = sqrt(taylor_2_bound);
+
+ static const char* function = "boost::math::atanh<%1%>(%1%)";
if (x < -one)
{
- return tools::domain_error<T>(
- BOOST_CURRENT_FUNCTION,
- "atanh requires x >= -1, but got x = %1%.", x);
+ return policy::raise_domain_error<T>(
+ function,
+ "atanh requires x >= -1, but got x = %1%.", x, pol);
}
else if (x < -one + tools::epsilon<T>())
{
// -Infinity:
- return -tools::overflow_error<T>(
- BOOST_CURRENT_FUNCTION);
+ return -policy::raise_overflow_error<T>(function, 0, pol);
}
else if (x > one - tools::epsilon<T>())
{
// Infinity:
- return -tools::overflow_error<T>(
- BOOST_CURRENT_FUNCTION);
+ return -policy::raise_overflow_error<T>(function, 0, pol);
}
else if (x > +one)
{
- return tools::domain_error<T>(
- BOOST_CURRENT_FUNCTION,
- "atanh requires x <= 1, but got x = %1%.", x);
+ return policy::raise_domain_error<T>(
+ function,
+ "atanh requires x <= 1, but got x = %1%.", x, pol);
}
else if (abs(x) >= taylor_n_bound)
{
@@ -99,12 +99,19 @@
}
}
+ template<typename T, typename Policy>
+ inline typename tools::promote_args<T>::type atanh(const T x, const Policy& pol)
+ {
+ typedef typename tools::promote_args<T>::type result_type;
+ return detail::atanh_imp(
+ static_cast<result_type>(x), pol);
+ }
template<typename T>
inline typename tools::promote_args<T>::type atanh(const T x)
{
typedef typename tools::promote_args<T>::type result_type;
return detail::atanh_imp(
- static_cast<result_type>(x));
+ static_cast<result_type>(x), policy::policy<>());
}
}
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/bessel.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/bessel.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/bessel.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -32,17 +32,19 @@
typedef mpl::int_<1> bessel_maybe_int_tag; // Maybe integer optimisation.
typedef mpl::int_<2> bessel_int_tag; // Definite integer optimistaion.
-template <class T1, class T2>
+template <class T1, class T2, class Policy>
struct bessel_traits
{
typedef typename tools::promote_args<
T1, T2
>::type result_type;
- typedef typename mpl::if_c<
- (std::numeric_limits<result_type>::is_specialized == 0)
- || (std::numeric_limits<result_type>::digits == 0)
- || (std::numeric_limits<result_type>::digits > 64),
+ typedef typename policy::precision<result_type, Policy>::type precision_type;
+
+ typedef typename mpl::if_<
+ mpl::or_<
+ mpl::less_equal<precision_type, mpl::int_<0> >,
+ mpl::greater<precision_type, mpl::int_<64> > >,
bessel_no_int_tag,
typename mpl::if_<
is_integral<T1>,
@@ -52,7 +54,7 @@
>::type optimisation_tag;
};
-template <class T>
+template <class T, class Policy>
struct bessel_j_small_z_series_term
{
typedef T result_type;
@@ -62,7 +64,7 @@
{
using namespace std;
mult = x / 2;
- term = pow(mult, v) / boost::math::tgamma(v+1);
+ term = pow(mult, v) / boost::math::tgamma(v+1, Policy());
mult *= -mult;
}
T operator()()
@@ -79,7 +81,7 @@
T term;
};
-template <class T>
+template <class T, class Policy>
struct sph_bessel_j_small_z_series_term
{
typedef T result_type;
@@ -89,7 +91,7 @@
{
using namespace std;
mult = x / 2;
- term = pow(mult, T(v)) / boost::math::tgamma(v+1+T(0.5f));
+ term = pow(mult, T(v)) / boost::math::tgamma(v+1+T(0.5f), Policy());
mult *= -mult;
}
T operator()()
@@ -106,83 +108,84 @@
T term;
};
-template <class T>
-inline T bessel_j_small_z_series(T v, T x)
+template <class T, class Policy>
+inline T bessel_j_small_z_series(T v, T x, const Policy& pol)
{
- bessel_j_small_z_series_term<T> s(v, x);
+ bessel_j_small_z_series_term<T, Policy> s(v, x);
boost::uintmax_t max_iter = BOOST_MATH_MAX_ITER;
- T result = boost::math::tools::sum_series(s, boost::math::tools::digits<T>(), max_iter);
- tools::check_series_iterations(BOOST_CURRENT_FUNCTION, max_iter);
+ T result = boost::math::tools::sum_series(s, boost::math::policy::digits<T, Policy>(), max_iter);
+ policy::check_series_iterations("boost::math::bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
return result;
}
-template <class T>
-inline T sph_bessel_j_small_z_series(unsigned v, T x)
+template <class T, class Policy>
+inline T sph_bessel_j_small_z_series(unsigned v, T x, const Policy& pol)
{
using namespace std; // ADL of std names
- sph_bessel_j_small_z_series_term<T> s(v, x);
+ sph_bessel_j_small_z_series_term<T, Policy> s(v, x);
boost::uintmax_t max_iter = BOOST_MATH_MAX_ITER;
- T result = boost::math::tools::sum_series(s, boost::math::tools::digits<T>(), max_iter);
- tools::check_series_iterations(BOOST_CURRENT_FUNCTION, max_iter);
+ T result = boost::math::tools::sum_series(s, boost::math::policy::digits<T, Policy>(), max_iter);
+ policy::check_series_iterations("boost::math::sph_bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
return result * sqrt(constants::pi<T>() / 4);
}
-template <class T>
-T cyl_bessel_j_imp(T v, T x, const bessel_no_int_tag& t)
+template <class T, class Policy>
+T cyl_bessel_j_imp(T v, T x, const bessel_no_int_tag& t, const Policy& pol)
{
using namespace std;
+ static const char* function = "boost::math::bessel_j<%1%>(%1%,%1%)";
if(x < 0)
{
// better have integer v:
if(floor(v) == v)
{
- T r = cyl_bessel_j_imp(v, -x, t);
+ T r = cyl_bessel_j_imp(v, -x, t, pol);
if(tools::real_cast<int>(v) & 1)
r = -r;
return r;
}
else
- return tools::domain_error<T>(
- BOOST_CURRENT_FUNCTION,
- "Got x = %1%, but we need x >= 0", x);
+ return policy::raise_domain_error<T>(
+ function,
+ "Got x = %1%, but we need x >= 0", x, pol);
}
if(x == 0)
return (v == 0) ? 1 : (v > 0) ? 0 :
- tools::domain_error<T>(
- BOOST_CURRENT_FUNCTION,
- "Got v = %1%, but require v >= 0 or a negative integer: the result would be complex.", v);
+ policy::raise_domain_error<T>(
+ function,
+ "Got v = %1%, but require v >= 0 or a negative integer: the result would be complex.", v, pol);
if((v >= 0) && ((x < 1) || (v > x * x / 4)))
{
- return bessel_j_small_z_series(v, x);
+ return bessel_j_small_z_series(v, x, pol);
}
T j, y;
- bessel_jy(v, x, &j, &y, need_j);
+ bessel_jy(v, x, &j, &y, need_j, pol);
return j;
}
-template <class T>
-inline T cyl_bessel_j_imp(T v, T x, const bessel_maybe_int_tag&)
+template <class T, class Policy>
+inline T cyl_bessel_j_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol)
{
using namespace std; // ADL of std names.
- typedef typename bessel_asymptotic_tag<T>::type tag_type;
+ typedef typename bessel_asymptotic_tag<T, Policy>::type tag_type;
if((fabs(v) < 200) && (floor(v) == v))
{
if(fabs(x) > asymptotic_bessel_j_limit<T>(v, tag_type()))
return asymptotic_bessel_j_large_x_2(v, x);
else
- return bessel_jn(tools::real_cast<int>(v), x);
+ return bessel_jn(tools::real_cast<int>(v), x, pol);
}
- return cyl_bessel_j_imp(v, x, bessel_no_int_tag());
+ return cyl_bessel_j_imp(v, x, bessel_no_int_tag(), pol);
}
-template <class T>
-inline T cyl_bessel_j_imp(int v, T x, const bessel_int_tag&)
+template <class T, class Policy>
+inline T cyl_bessel_j_imp(int v, T x, const bessel_int_tag&, const Policy& pol)
{
using namespace std;
- typedef typename bessel_asymptotic_tag<T>::type tag_type;
+ typedef typename bessel_asymptotic_tag<T, Policy>::type tag_type;
if(fabs(x) > asymptotic_bessel_j_limit<T>(abs(v), tag_type()))
{
T r = asymptotic_bessel_j_large_x_2(static_cast<T>(abs(v)), x);
@@ -191,37 +194,37 @@
return r;
}
else
- return bessel_jn(v, x);
+ return bessel_jn(v, x, pol);
}
-template <class T>
-inline T sph_bessel_j_imp(unsigned n, T x)
+template <class T, class Policy>
+inline T sph_bessel_j_imp(unsigned n, T x, const Policy& pol)
{
using namespace std; // ADL of std names
if(x < 0)
- return tools::domain_error<T>(
- BOOST_CURRENT_FUNCTION,
- "Got x = %1%, but function requires x > 0.", x);
+ return policy::raise_domain_error<T>(
+ "boost::math::sph_bessel_j<%1%>(%1%,%1%)",
+ "Got x = %1%, but function requires x > 0.", x, pol);
//
// Special case, n == 0 resolves down to the sinus cardinal of x:
//
if(n == 0)
- return boost::math::sinc_pi(x);
+ return boost::math::sinc_pi(x, pol);
//
// When x is small we may end up with 0/0, use series evaluation
// instead, especially as it converges rapidly:
//
if(x < 1)
- return sph_bessel_j_small_z_series(n, x);
+ return sph_bessel_j_small_z_series(n, x, pol);
//
// Default case is just a naive evaluation of the definition:
//
return sqrt(constants::pi<T>() / (2 * x))
- * cyl_bessel_j_imp(T(n)+T(0.5f), x, bessel_no_int_tag());
+ * cyl_bessel_j_imp(T(n)+T(0.5f), x, bessel_no_int_tag(), pol);
}
-template <class T>
-T cyl_bessel_i_imp(T v, T x)
+template <class T, class Policy>
+T cyl_bessel_i_imp(T v, T x, const Policy& pol)
{
//
// This handles all the bessel I functions, note that we don't optimise
@@ -235,15 +238,15 @@
// better have integer v:
if(floor(v) == v)
{
- T r = cyl_bessel_i_imp(v, -x);
+ T r = cyl_bessel_i_imp(v, -x, pol);
if(tools::real_cast<int>(v) & 1)
r = -r;
return r;
}
else
- return tools::domain_error<T>(
- BOOST_CURRENT_FUNCTION,
- "Got x = %1%, but we need x >= 0", x);
+ return policy::raise_domain_error<T>(
+ "boost::math::cyl_bessel_i<%1%>(%1%,%1%)",
+ "Got x = %1%, but we need x >= 0", x, pol);
}
if(x == 0)
{
@@ -255,7 +258,7 @@
T e = exp(x / 2);
return e * (e / sqrt(2 * x * constants::pi<T>()));
}
- if(tools::digits<T>() <= 64)
+ if(policy::digits<T, Policy>() <= 64)
{
if(v == 0)
{
@@ -267,77 +270,79 @@
}
}
T I, K;
- bessel_ik(v, x, &I, &K, need_i);
+ bessel_ik(v, x, &I, &K, need_i, pol);
return I;
}
-template <class T>
-inline T cyl_bessel_k_imp(T v, T x, const bessel_no_int_tag& t)
+template <class T, class Policy>
+inline T cyl_bessel_k_imp(T v, T x, const bessel_no_int_tag& t, const Policy& pol)
{
+ static const char* function = "boost::math::cyl_bessel_k<%1%>(%1%,%1%)";
using namespace std;
if(x < 0)
{
- return tools::domain_error<T>(
- BOOST_CURRENT_FUNCTION,
- "Got x = %1%, but we need x > 0", x);
+ return policy::raise_domain_error<T>(
+ function,
+ "Got x = %1%, but we need x > 0", x, pol);
}
if(x == 0)
{
- return (v == 0) ? tools::overflow_error<T>(BOOST_CURRENT_FUNCTION)
- : tools::domain_error<T>(
- BOOST_CURRENT_FUNCTION,
- "Got x = %1%, but we need x > 0", x);
+ return (v == 0) ? policy::raise_overflow_error<T>(function, 0, pol)
+ : policy::raise_domain_error<T>(
+ function,
+ "Got x = %1%, but we need x > 0", x, pol);
}
T I, K;
- bessel_ik(v, x, &I, &K, need_k);
+ bessel_ik(v, x, &I, &K, need_k, pol);
return K;
}
-template <class T>
-inline T cyl_bessel_k_imp(T v, T x, const bessel_maybe_int_tag&)
+template <class T, class Policy>
+inline T cyl_bessel_k_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol)
{
using namespace std;
if((floor(v) == v))
{
- return bessel_kn(tools::real_cast<int>(v), x);
+ return bessel_kn(tools::real_cast<int>(v), x, pol);
}
- return cyl_bessel_k_imp(v, x, bessel_no_int_tag());
+ return cyl_bessel_k_imp(v, x, bessel_no_int_tag(), pol);
}
-template <class T>
-inline T cyl_bessel_k_imp(int v, T x, const bessel_int_tag&)
+template <class T, class Policy>
+inline T cyl_bessel_k_imp(int v, T x, const bessel_int_tag&, const Policy& pol)
{
- return bessel_kn(v, x);
+ return bessel_kn(v, x, pol);
}
-template <class T>
-inline T cyl_neumann_imp(T v, T x, const bessel_no_int_tag&)
+template <class T, class Policy>
+inline T cyl_neumann_imp(T v, T x, const bessel_no_int_tag&, const Policy& pol)
{
+ static const char* function = "boost::math::cyl_neumann<%1%>(%1%,%1%)";
if(x <= 0)
{
return (v == 0) && (x == 0) ?
- tools::overflow_error<T>(BOOST_CURRENT_FUNCTION)
- : tools::domain_error<T>(
- BOOST_CURRENT_FUNCTION,
- "Got x = %1%, but result is complex for x <= 0", x);
+ policy::raise_overflow_error<T>(function, 0, pol)
+ : policy::raise_domain_error<T>(
+ function,
+ "Got x = %1%, but result is complex for x <= 0", x, pol);
}
T j, y;
- bessel_jy(v, x, &j, &y, need_y);
+ bessel_jy(v, x, &j, &y, need_y, pol);
//
// Post evaluation check for internal overflow during evaluation,
// can occur when x is small and v is large, in which case the result
// is -INF:
//
if(!(boost::math::isfinite)(y))
- return -tools::overflow_error<T>(BOOST_CURRENT_FUNCTION);
+ return -policy::raise_overflow_error<T>(function, 0, pol);
return y;
}
-template <class T>
-inline T cyl_neumann_imp(T v, T x, const bessel_maybe_int_tag&)
+template <class T, class Policy>
+inline T cyl_neumann_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol)
{
using namespace std;
- typedef typename bessel_asymptotic_tag<T>::type tag_type;
+ typedef typename bessel_asymptotic_tag<T, Policy>::type tag_type;
if(floor(v) == v)
{
if((fabs(x) > asymptotic_bessel_y_limit<T>(tag_type())) && (fabs(x) > 5 * abs(v)))
@@ -348,16 +353,16 @@
return r;
}
else
- return bessel_yn(tools::real_cast<int>(v), x);
+ return bessel_yn(tools::real_cast<int>(v), x, pol);
}
- return cyl_neumann_imp<T>(v, x, bessel_no_int_tag());
+ return cyl_neumann_imp<T>(v, x, bessel_no_int_tag(), pol);
}
-template <class T>
-inline T cyl_neumann_imp(int v, T x, const bessel_int_tag&)
+template <class T, class Policy>
+inline T cyl_neumann_imp(int v, T x, const bessel_int_tag&, const Policy& pol)
{
using namespace std;
- typedef typename bessel_asymptotic_tag<T>::type tag_type;
+ typedef typename bessel_asymptotic_tag<T, Policy>::type tag_type;
if((fabs(x) > asymptotic_bessel_y_limit<T>(tag_type())) && (fabs(x) > 5 * abs(v)))
{
T r = asymptotic_bessel_y_large_x_2(static_cast<T>(abs(v)), x);
@@ -366,91 +371,128 @@
return r;
}
else
- return bessel_yn(tools::real_cast<int>(v), x);
+ return bessel_yn(tools::real_cast<int>(v), x, pol);
}
-template <class T>
-inline T sph_neumann_imp(unsigned v, T x)
+template <class T, class Policy>
+inline T sph_neumann_imp(unsigned v, T x, const Policy& pol)
{
using namespace std; // ADL of std names
+ static const char* function = "boost::math::sph_neumann<%1%>(%1%,%1%)";
//
// Nothing much to do here but check for errors, and
// evaluate the function's definition directly:
//
if(x < 0)
- return tools::domain_error<T>(
- BOOST_CURRENT_FUNCTION,
- "Got x = %1%, but function requires x > 0.", x);
+ return policy::raise_domain_error<T>(
+ function,
+ "Got x = %1%, but function requires x > 0.", x, pol);
if(x < 2 * tools::min_value<T>())
- return -tools::overflow_error<T>(BOOST_CURRENT_FUNCTION);
+ return -policy::raise_overflow_error<T>(function, 0, pol);
- T result = cyl_neumann_imp(T(v)+0.5f, x, bessel_no_int_tag());
+ T result = cyl_neumann_imp(T(v)+0.5f, x, bessel_no_int_tag(), pol);
T tx = sqrt(constants::pi<T>() / (2 * x));
if((tx > 1) && (tools::max_value<T>() / tx < result))
- return -tools::overflow_error<T>(BOOST_CURRENT_FUNCTION);
+ return -policy::raise_overflow_error<T>(function, 0, pol);
return result * tx;
}
} // namespace detail
+template <class T1, class T2, class Policy>
+inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_j(T1 v, T2 x, const Policy& pol)
+{
+ BOOST_FPU_EXCEPTION_GUARD
+ 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 policy::evaluation<result_type, Policy>::type value_type;
+ return policy::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%)");
+}
+
template <class T1, class T2>
-inline typename detail::bessel_traits<T1, T2>::result_type cyl_bessel_j(T1 v, T2 x)
+inline typename detail::bessel_traits<T1, T2, policy::policy<> >::result_type cyl_bessel_j(T1 v, T2 x)
+{
+ return cyl_bessel_j(v, x, policy::policy<>());
+}
+
+template <class T, class Policy>
+inline typename detail::bessel_traits<T, T, Policy>::result_type sph_bessel(unsigned v, T x, const Policy& pol)
{
BOOST_FPU_EXCEPTION_GUARD
- typedef typename detail::bessel_traits<T1, T2>::result_type result_type;
- typedef typename detail::bessel_traits<T1, T2>::optimisation_tag tag_type;
- typedef typename tools::evaluation<result_type>::type value_type;
- return tools::checked_narrowing_cast<result_type>(detail::cyl_bessel_j_imp<value_type>(v, static_cast<value_type>(x), tag_type()), BOOST_CURRENT_FUNCTION);
+ typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
+ typedef typename policy::evaluation<result_type, Policy>::type value_type;
+ return policy::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%)");
}
template <class T>
-inline typename detail::bessel_traits<T, T>::result_type sph_bessel(unsigned v, T x)
+inline typename detail::bessel_traits<T, T, policy::policy<> >::result_type sph_bessel(unsigned v, T x)
+{
+ return sph_bessel(v, x, policy::policy<>());
+}
+
+template <class T1, class T2, class Policy>
+inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_i(T1 v, T2 x, const Policy& pol)
{
BOOST_FPU_EXCEPTION_GUARD
- typedef typename detail::bessel_traits<T, T>::result_type result_type;
- typedef typename tools::evaluation<result_type>::type value_type;
- return tools::checked_narrowing_cast<result_type>(detail::sph_bessel_j_imp<value_type>(v, static_cast<value_type>(x)), BOOST_CURRENT_FUNCTION);
+ typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
+ typedef typename policy::evaluation<result_type, Policy>::type value_type;
+ return policy::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%)");
}
template <class T1, class T2>
-inline typename detail::bessel_traits<T1, T2>::result_type cyl_bessel_i(T1 v, T2 x)
+inline typename detail::bessel_traits<T1, T2, policy::policy<> >::result_type cyl_bessel_i(T1 v, T2 x)
+{
+ return cyl_bessel_i(v, x, policy::policy<>());
+}
+
+template <class T1, class T2, class Policy>
+inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_k(T1 v, T2 x, const Policy& pol)
{
BOOST_FPU_EXCEPTION_GUARD
- typedef typename detail::bessel_traits<T1, T2>::result_type result_type;
- typedef typename tools::evaluation<result_type>::type value_type;
- return tools::checked_narrowing_cast<result_type>(detail::cyl_bessel_i_imp<value_type>(v, static_cast<value_type>(x)), BOOST_CURRENT_FUNCTION);
+ 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 policy::evaluation<result_type, Policy>::type value_type;
+ return policy::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%)");
}
template <class T1, class T2>
-inline typename detail::bessel_traits<T1, T2>::result_type cyl_bessel_k(T1 v, T2 x)
+inline typename detail::bessel_traits<T1, T2, policy::policy<> >::result_type cyl_bessel_k(T1 v, T2 x)
+{
+ return cyl_bessel_k(v, x, policy::policy<>());
+}
+
+template <class T1, class T2, class Policy>
+inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_neumann(T1 v, T2 x, const Policy& pol)
{
BOOST_FPU_EXCEPTION_GUARD
- typedef typename detail::bessel_traits<T1, T2>::result_type result_type;
- typedef typename detail::bessel_traits<T1, T2>::optimisation_tag tag_type;
- typedef typename tools::evaluation<result_type>::type value_type;
- return tools::checked_narrowing_cast<result_type>(detail::cyl_bessel_k_imp<value_type>(v, static_cast<value_type>(x), tag_type()), BOOST_CURRENT_FUNCTION);
+ 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 policy::evaluation<result_type, Policy>::type value_type;
+ return policy::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%)");
}
template <class T1, class T2>
-inline typename detail::bessel_traits<T1, T2>::result_type cyl_neumann(T1 v, T2 x)
+inline typename detail::bessel_traits<T1, T2, policy::policy<> >::result_type cyl_neumann(T1 v, T2 x)
+{
+ return cyl_neumann(v, x, policy::policy<>());
+}
+
+template <class T, class Policy>
+inline typename detail::bessel_traits<T, T, Policy>::result_type sph_neumann(unsigned v, T x, const Policy& pol)
{
BOOST_FPU_EXCEPTION_GUARD
- typedef typename detail::bessel_traits<T1, T2>::result_type result_type;
- typedef typename detail::bessel_traits<T1, T2>::optimisation_tag tag_type;
- typedef typename tools::evaluation<result_type>::type value_type;
- return tools::checked_narrowing_cast<result_type>(detail::cyl_neumann_imp<value_type>(v, static_cast<value_type>(x), tag_type()), BOOST_CURRENT_FUNCTION);
+ typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
+ typedef typename policy::evaluation<result_type, Policy>::type value_type;
+ return policy::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%)");
}
template <class T>
-inline typename detail::bessel_traits<T, T>::result_type sph_neumann(unsigned v, T x)
+inline typename detail::bessel_traits<T, T, policy::policy<> >::result_type sph_neumann(unsigned v, T x)
{
- BOOST_FPU_EXCEPTION_GUARD
- typedef typename detail::bessel_traits<T, T>::result_type result_type;
- typedef typename tools::evaluation<result_type>::type value_type;
- return tools::checked_narrowing_cast<result_type>(detail::sph_neumann_imp<value_type>(v, static_cast<value_type>(x)), BOOST_CURRENT_FUNCTION);
+ return sph_neumann(v, x, policy::policy<>());
}
} // namespace math
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/cos_pi.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/cos_pi.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/cos_pi.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -40,6 +40,12 @@
return invert ? -rem : rem;
}
+template <class T, class Policy>
+inline T cos_pi(T x, const Policy&)
+{
+ return cos_pi(x);
+}
+
} // namespace math
} // namespace boost
#endif
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_ik.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_ik.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_ik.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -9,7 +9,7 @@
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/special_functions/sin_pi.hpp>
#include <boost/math/constants/constants.hpp>
-#include <boost/math/tools/error_handling.hpp>
+#include <boost/math/policy/error_handling.hpp>
#include <boost/math/tools/config.hpp>
// Modified Bessel functions of the first and second kind of fractional order
@@ -20,8 +20,8 @@
// Calculate K(v, x) and K(v+1, x) by method analogous to
// Temme, Journal of Computational Physics, vol 21, 343 (1976)
-template <typename T>
-int temme_ik(T v, T x, T* K, T* K1)
+template <typename T, typename Policy>
+int temme_ik(T v, T x, T* K, T* K1, const Policy& pol)
{
T f, h, p, q, coef, sum, sum1, tolerance;
T a, b, c, d, sigma, gamma1, gamma2;
@@ -37,8 +37,8 @@
BOOST_ASSERT(abs(x) <= 2);
BOOST_ASSERT(abs(v) <= 0.5f);
- T gp = tgamma1pm1(v);
- T gm = tgamma1pm1(-v);
+ T gp = tgamma1pm1(v, pol);
+ T gm = tgamma1pm1(-v, pol);
a = log(x / 2);
b = exp(v * a);
@@ -76,7 +76,7 @@
break;
}
}
- check_series_iterations(BOOST_CURRENT_FUNCTION, k);
+ policy::check_series_iterations("boost::math::bessel_ik<%1%>(%1%,%1%)", k, pol);
*K = sum;
*K1 = 2 * sum1 / x;
@@ -86,8 +86,8 @@
// Evaluate continued fraction fv = I_(v+1) / I_v, derived from
// Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73
-template <typename T>
-int CF1_ik(T v, T x, T* fv)
+template <typename T, typename Policy>
+int CF1_ik(T v, T x, T* fv, const Policy& pol)
{
T C, D, f, a, b, delta, tiny, tolerance;
int k;
@@ -119,7 +119,7 @@
break;
}
}
- tools::check_series_iterations(BOOST_CURRENT_FUNCTION, k);
+ policy::check_series_iterations("boost::math::bessel_ik<%1%>(%1%,%1%)", k, pol);
*fv = f;
@@ -129,8 +129,8 @@
// Calculate K(v, x) and K(v+1, x) by evaluating continued fraction
// z1 / z0 = U(v+1.5, 2v+1, 2x) / U(v+0.5, 2v+1, 2x), see
// Thompson and Barnett, Computer Physics Communications, vol 47, 245 (1987)
-template <typename T>
-int CF2_ik(T v, T x, T* Kv, T* Kv1)
+template <typename T, typename Policy>
+int CF2_ik(T v, T x, T* Kv, T* Kv1, const Policy& pol)
{
using namespace std;
using namespace boost::math::constants;
@@ -177,7 +177,7 @@
break;
}
}
- tools::check_series_iterations(BOOST_CURRENT_FUNCTION, k);
+ policy::check_series_iterations("boost::math::bessel_ik<%1%>(%1%,%1%)", k, pol);
*Kv = sqrt(pi<T>() / (2 * x)) * exp(-x) / S;
*Kv1 = *Kv * (0.5f + v + x + (v * v - 0.25f) * f) / x;
@@ -192,8 +192,8 @@
// Compute I(v, x) and K(v, x) simultaneously by Temme's method, see
// Temme, Journal of Computational Physics, vol 19, 324 (1975)
-template <typename T>
-int bessel_ik(T v, T x, T* I, T* K, int kind = need_i | need_k)
+template <typename T, typename Policy>
+int bessel_ik(T v, T x, T* I, T* K, int kind, const Policy& pol)
{
// Kv1 = K_(v+1), fv = I_(v+1) / I_v
// Ku1 = K_(u+1), fu = I_(u+1) / I_u
@@ -206,6 +206,8 @@
using namespace boost::math::tools;
using namespace boost::math::constants;
+ static const char* function = "boost::math::bessel_ik<%1%>(%1%,%1%)";
+
if (v < 0)
{
reflect = true;
@@ -217,8 +219,8 @@
if (x < 0)
{
- *I = *K = domain_error<T>(BOOST_CURRENT_FUNCTION,
- "Got x = %1% but real argument x must be non-negative, complex number result not supported.", x);
+ *I = *K = policy::raise_domain_error<T>(function,
+ "Got x = %1% but real argument x must be non-negative, complex number result not supported.", x, pol);
return 1;
}
if (x == 0)
@@ -226,7 +228,7 @@
Iv = (v == 0) ? static_cast<T>(1) : static_cast<T>(0);
if(kind & need_k)
{
- Kv = tools::overflow_error<T>(BOOST_CURRENT_FUNCTION);
+ Kv = policy::raise_overflow_error<T>(function, 0, pol);
}
else
{
@@ -236,9 +238,9 @@
if(reflect && (kind & need_i))
{
T z = (u + n % 2);
- Iv = sin_pi(z) == 0 ?
+ Iv = sin_pi(z, pol) == 0 ?
Iv :
- tools::overflow_error<T>(BOOST_CURRENT_FUNCTION); // reflection formula
+ policy::raise_overflow_error<T>(function, 0, pol); // reflection formula
}
*I = Iv;
@@ -250,11 +252,11 @@
W = 1 / x; // Wronskian
if (x <= 2) // x in (0, 2]
{
- temme_ik(u, x, &Ku, &Ku1); // Temme series
+ temme_ik(u, x, &Ku, &Ku1, pol); // Temme series
}
else // x in (2, \infty)
{
- CF2_ik(u, x, &Ku, &Ku1); // continued fraction CF2_ik
+ CF2_ik(u, x, &Ku, &Ku1, pol); // continued fraction CF2_ik
}
prev = Ku;
current = Ku1;
@@ -279,11 +281,11 @@
// x case instead. Note that the asymptotic expansion
// isn't very accurate - so it's deliberately very hard
// to get here - probably we're going to overflow:
- Iv = asymptotic_bessel_i_large_x(v, x);
+ Iv = asymptotic_bessel_i_large_x(v, x, pol);
}
else
{
- CF1_ik(v, x, &fv); // continued fraction CF1_ik
+ CF1_ik(v, x, &fv, pol); // continued fraction CF1_ik
Iv = W / (Kv * fv + Kv1); // Wronskian relation
}
}
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_jn.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_jn.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_jn.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -17,8 +17,8 @@
namespace boost { namespace math { namespace detail{
-template <typename T>
-T bessel_jn(int n, T x)
+template <typename T, typename Policy>
+T bessel_jn(int n, T x, const Policy& pol)
{
T value(0), factor, current, prev, next;
@@ -62,7 +62,7 @@
{
T fn; int s; // fn = J_(n+1) / J_n
// |x| <= n, fast convergence for continued fraction CF1
- boost::math::detail::CF1_jy(static_cast<T>(n), x, &fn, &s);
+ boost::math::detail::CF1_jy(static_cast<T>(n), x, &fn, &s, pol);
// tiny initial value to prevent overflow
T init = sqrt(tools::min_value<T>());
prev = fn * init;
Deleted: sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_jv.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_jv.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
+++ (empty file)
@@ -1,122 +0,0 @@
-// Copyright (c) 2006 Xiaogang Zhang
-// 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)
-
-#ifndef BOOST_MATH_BESSEL_JV_HPP
-#define BOOST_MATH_BESSEL_JV_HPP
-
-#include <boost/math/special_functions/bessel_jn.hpp>
-#include <boost/math/special_functions/bessel_jy.hpp>
-#include <boost/math/constants/constants.hpp>
-
-// Bessel and Spherical Bessel functions of the first kind
-
-namespace boost { namespace math {
-
-// Bessel function of the first kind of any order
-template <typename T>
-inline T bessel_jv(T v, T x)
-{
- int n = static_cast<int>(v);
- if (v == n)
- {
- return bessel_jn(n, x); // v is integer
- }
- else
- {
- T J, Y;
- bessel_jy(v, x, &J, &Y);
- return J;
- }
-}
-
-// Spherical Bessel function of the first kind of non-negative order
-template <typename T>
-inline T bessel_sjv(unsigned n, T x)
-{
- using namespace std;
- using namespace boost::math::constants;
-
- T v = n + 0.5L;
- if (x == 0)
- {
- return (n == 0) ? static_cast<T>(1) : static_cast<T>(0);
- }
- else
- {
- return sqrt(0.5L * pi<T>() / x) * bessel_jv(v, x);
- }
-}
-
-// -------------------- TR1 functions --------------------
-
-inline float cyl_bessel_jf(float nu, float x)
-{
- if (nu >= 128)
- {
- std::cout << "Warning: cyl_bessel_jf(nu, x), nu >= 128, "
- << "result is implementation defined according to C++ TR1"
- << std::endl;
- }
- return bessel_jv<float>(nu, x);
-}
-
-inline double cyl_bessel_j(double nu, double x)
-{
- if (nu >= 128)
- {
- std::cout << "Warning: cyl_bessel_j(nu, x), nu >= 128, "
- << "result is implementation defined according to C++ TR1"
- << std::endl;
- }
- return bessel_jv<double>(nu, x);
-}
-
-inline long double cyl_bessel_jl(long double nu, long double x)
-{
- if (nu >= 128)
- {
- std::cout << "Warning: cyl_bessel_jl(nu, x), nu >= 128, "
- << "result is implementation defined according to C++ TR1"
- << std::endl;
- }
- return bessel_jv<long double>(nu, x);
-}
-
-inline float sph_besself(unsigned n, float x)
-{
- if (n >= 128)
- {
- std::cout << "Warning: sph_besself(n, x), n >= 128, "
- << "result is implementation defined according to C++ TR1"
- << std::endl;
- }
- return bessel_sjv<float>(n, x);
-}
-
-inline double sph_bessel(unsigned n, double x)
-{
- if (n >= 128)
- {
- std::cout << "Warning: sph_bessel(n, x), n >= 128, "
- << "result is implementation defined according to C++ TR1"
- << std::endl;
- }
- return bessel_sjv<double>(n, x);
-}
-
-inline long double sph_bessell(unsigned n, long double x)
-{
- if (n >= 128)
- {
- std::cout << "Warning: sph_bessell(n, x), n >= 128, "
- << "result is implementation defined according to C++ TR1"
- << std::endl;
- }
- return bessel_sjv<long double>(n, x);
-}
-
-}} // namespaces
-
-#endif // BOOST_MATH_BESSEL_JV_HPP
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_jy.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_jy.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_jy.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -15,7 +15,7 @@
#include <boost/math/special_functions/detail/simple_complex.hpp>
#include <boost/math/special_functions/detail/bessel_jy_asym.hpp>
#include <boost/math/constants/constants.hpp>
-#include <boost/math/tools/error_handling.hpp>
+#include <boost/math/policy/error_handling.hpp>
#include <boost/mpl/if.hpp>
#include <boost/type_traits/is_floating_point.hpp>
#include <complex>
@@ -28,8 +28,8 @@
// Calculate Y(v, x) and Y(v+1, x) by Temme's method, see
// Temme, Journal of Computational Physics, vol 21, 343 (1976)
-template <typename T>
-int temme_jy(T v, T x, T* Y, T* Y1)
+template <typename T, typename Policy>
+int temme_jy(T v, T x, T* Y, T* Y1, const Policy& pol)
{
T g, h, p, q, f, coef, sum, sum1, tolerance;
T a, d, e, sigma;
@@ -41,10 +41,10 @@
BOOST_ASSERT(fabs(v) <= 0.5f); // precondition for using this routine
- T gp = tgamma1pm1(v);
- T gm = tgamma1pm1(-v);
- T spv = sin_pi(v);
- T spv2 = sin_pi(v/2);
+ T gp = tgamma1pm1(v, pol);
+ T gm = tgamma1pm1(-v, pol);
+ T spv = sin_pi(v, pol);
+ T spv2 = sin_pi(v/2, pol);
T xp = pow(x/2, v);
a = log(x / 2);
@@ -88,7 +88,7 @@
break;
}
}
- check_series_iterations(BOOST_CURRENT_FUNCTION, k);
+ policy::check_series_iterations("boost::math::bessel_jy<%1%>(%1%,%1%)", k, pol);
*Y = -sum;
*Y1 = -2 * sum1 / x;
@@ -97,8 +97,8 @@
// Evaluate continued fraction fv = J_(v+1) / J_v, see
// Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73
-template <typename T>
-int CF1_jy(T v, T x, T* fv, int* sign)
+template <typename T, typename Policy>
+int CF1_jy(T v, T x, T* fv, int* sign, const Policy& pol)
{
T C, D, f, a, b, delta, tiny, tolerance;
int k, s = 1;
@@ -129,7 +129,7 @@
if (abs(delta - 1.0L) < tolerance)
{ break; }
}
- tools::check_series_iterations(BOOST_CURRENT_FUNCTION, k / 100);
+ policy::check_series_iterations("boost::math::bessel_jy<%1%>(%1%,%1%)", k / 100, pol);
*fv = -f;
*sign = s; // sign of denominator
@@ -145,8 +145,8 @@
// Evaluate continued fraction p + iq = (J' + iY') / (J + iY), see
// Press et al, Numerical Recipes in C, 2nd edition, 1992
-template <typename T>
-int CF2_jy(T v, T x, T* p, T* q)
+template <typename T, typename Policy>
+int CF2_jy(T v, T x, T* p, T* q, const Policy& pol)
{
using namespace std;
@@ -183,7 +183,7 @@
f *= delta;
if (abs(delta - one) < tolerance) { break; }
}
- tools::check_series_iterations(BOOST_CURRENT_FUNCTION, k);
+ policy::check_series_iterations("boost::math::bessel_jy<%1%>(%1%,%1%)", k, pol);
*p = real(f);
*q = imag(f);
@@ -197,8 +197,8 @@
// Compute J(v, x) and Y(v, x) simultaneously by Steed's method, see
// Barnett et al, Computer Physics Communications, vol 8, 377 (1974)
-template <typename T>
-int bessel_jy(T v, T x, T* J, T* Y, int kind = need_j|need_y)
+template <typename T, typename Policy>
+int bessel_jy(T v, T x, T* J, T* Y, int kind, const Policy& pol)
{
BOOST_ASSERT(x >= 0);
@@ -207,6 +207,8 @@
bool reflect = false;
int n, k, s;
+ static const char* function = "boost::math::bessel_jy<%1%>(%1%,%1%)";
+
using namespace std;
using namespace boost::math::tools;
using namespace boost::math::constants;
@@ -222,8 +224,8 @@
if (x == 0)
{
- *J = *Y = tools::overflow_error<T>(
- BOOST_CURRENT_FUNCTION);
+ *J = *Y = policy::raise_overflow_error<T>(
+ function, 0, pol);
return 1;
}
@@ -231,7 +233,7 @@
W = T(2) / (x * pi<T>()); // Wronskian
if (x <= 2) // x in (0, 2]
{
- if(temme_jy(u, x, &Yu, &Yu1)) // Temme series
+ if(temme_jy(u, x, &Yu, &Yu1, pol)) // Temme series
{
// domain error:
*J = *Y = Yu;
@@ -249,7 +251,7 @@
Yv1 = current;
if(kind&need_j)
{
- CF1_jy(v, x, &fv, &s); // continued fraction CF1_jy
+ CF1_jy(v, x, &fv, &s, pol); // continued fraction CF1_jy
Jv = W / (Yv * fv - Yv1); // Wronskian relation
}
else
@@ -259,7 +261,7 @@
{
// Get Y(u, x):
// define tag type that will dispatch to right limits:
- typedef typename bessel_asymptotic_tag<T>::type tag_type;
+ typedef typename bessel_asymptotic_tag<T, Policy>::type tag_type;
T lim;
switch(kind)
@@ -294,7 +296,7 @@
}
else
{
- CF1_jy(v, x, &fv, &s);
+ CF1_jy(v, x, &fv, &s, pol);
// tiny initial value to prevent overflow
T init = sqrt(tools::min_value<T>());
prev = fv * s * init;
@@ -308,7 +310,7 @@
T ratio = (s * init) / current; // scaling ratio
// can also call CF1_jy() to get fu, not much difference in precision
fu = prev / current;
- CF2_jy(u, x, &p, &q); // continued fraction CF2_jy
+ CF2_jy(u, x, &p, &q, pol); // continued fraction CF2_jy
T t = u / x - fu; // t = J'/J
gamma = (p - t) / q;
Ju = sign(current) * sqrt(W / (q + gamma * (p - t)));
@@ -338,8 +340,8 @@
if (reflect)
{
T z = (u + n % 2);
- *J = cos_pi(z) * Jv - sin_pi(z) * Yv; // reflection formula
- *Y = sin_pi(z) * Jv + cos_pi(z) * Yv;
+ *J = cos_pi(z, pol) * Jv - sin_pi(z, pol) * Yv; // reflection formula
+ *Y = sin_pi(z, pol) * Jv + cos_pi(z, pol) * Yv;
}
else
{
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_jy_asym.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_jy_asym.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_jy_asym.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -178,18 +178,20 @@
return 1245243 /*3128000*/;
}
-template <class T>
+template <class T, class Policy>
struct bessel_asymptotic_tag
{
- typedef typename mpl::if_c<
- (std::numeric_limits<T>::digits == 0
- || std::numeric_limits<T>::digits > 113),
+ typedef typename policy::precision<T, Policy>::type precision_type;
+ typedef typename mpl::if_<
+ mpl::or_<
+ mpl::equal_to<precision_type, mpl::int_<0> >,
+ mpl::greater<precision_type, mpl::int_<113> > >,
mpl::int_<0>,
- typename mpl::if_c<
- (std::numeric_limits<T>::digits > 64),
+ typename mpl::if_<
+ mpl::greater<precision_type, mpl::int_<64> >,
mpl::int_<113>,
- typename mpl::if_c<
- (std::numeric_limits<T>::digits > 53),
+ typename mpl::if_<
+ mpl::greater<precision_type, mpl::int_<53> >,
mpl::int_<64>,
mpl::int_<53>
>::type
@@ -227,14 +229,14 @@
return v2 * 39154 /*85700*/;
}
-template <class T>
-void temme_asyptotic_y_small_x(T v, T x, T* Y, T* Y1)
+template <class T, class Policy>
+void temme_asyptotic_y_small_x(T v, T x, T* Y, T* Y1, const Policy& pol)
{
T c = 1;
- T p = (v / sin_pi(v)) * pow(x / 2, -v) / tgamma(1 - v);
- T q = (v / sin_pi(v)) * pow(x / 2, v) / tgamma(1 + v);
+ T p = (v / sin_pi(v, pol)) * pow(x / 2, -v) / tgamma(1 - v, pol);
+ T q = (v / sin_pi(v, pol)) * pow(x / 2, v) / tgamma(1 + v, pol);
T f = (p - q) / v;
- T g_prefix = sin_pi(v / 2);
+ T g_prefix = sin_pi(v / 2, pol);
g_prefix *= g_prefix * 2 / v;
T g = f + g_prefix * q;
T h = p;
@@ -248,7 +250,7 @@
p /= k - v;
q /= k + v;
c *= c_mult / k;
- T c1 = pow(-x * x / 4, k) / factorial<T>(k);
+ T c1 = pow(-x * x / 4, k) / factorial<T>(k, pol);
g = f + g_prefix * q;
h = -k * g + p;
y += c * g;
@@ -261,8 +263,8 @@
*Y1 = (-2 / x) * y1;
}
-template <class T>
-T asymptotic_bessel_i_large_x(T v, T x)
+template <class T, class Policy>
+T asymptotic_bessel_i_large_x(T v, T x, const Policy& pol)
{
using namespace std; // ADL of std names
T s = 1;
@@ -287,7 +289,7 @@
s = e * (e * s / sqrt(2 * x * constants::pi<T>()));
return (boost::math::isfinite)(s) ?
- s : tools::overflow_error<T>(BOOST_CURRENT_FUNCTION);
+ s : policy::raise_overflow_error<T>("boost::math::asymptotic_bessel_i_large_x<%1%>(%1%,%1%)", 0, pol);
}
}}} // namespaces
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_k0.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_k0.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_k0.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -7,7 +7,7 @@
#define BOOST_MATH_BESSEL_K0_HPP
#include <boost/math/tools/rational.hpp>
-#include <boost/math/tools/error_handling.hpp>
+#include <boost/math/policy/error_handling.hpp>
#include <boost/assert.hpp>
// Modified Bessel function of the second kind of order zero
@@ -16,8 +16,8 @@
namespace boost { namespace math { namespace detail{
-template <typename T>
-T bessel_k0(T x)
+template <typename T, typename Policy>
+T bessel_k0(T x, const Policy& pol)
{
static const T P1[] = {
static_cast<T>(2.4708152720399552679e+03L),
@@ -75,14 +75,16 @@
using namespace std;
using namespace boost::math::tools;
+ static const char* function = "boost::math::bessel_k0<%1%>(%1%,%1%)";
+
if (x < 0)
{
- return domain_error<T>(BOOST_CURRENT_FUNCTION,
- "Got x = %1%, but argument x must be non-negative, complex number result not supported", x);
+ return policy::raise_domain_error<T>(function,
+ "Got x = %1%, but argument x must be non-negative, complex number result not supported", x, pol);
}
if (x == 0)
{
- return overflow_error<T>(BOOST_CURRENT_FUNCTION);
+ return policy::raise_overflow_error<T>(function, 0, pol);
}
if (x <= 1) // x in (0, 1]
{
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_k1.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_k1.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_k1.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -7,7 +7,7 @@
#define BOOST_MATH_BESSEL_K1_HPP
#include <boost/math/tools/rational.hpp>
-#include <boost/math/tools/error_handling.hpp>
+#include <boost/math/policy/error_handling.hpp>
#include <boost/assert.hpp>
// Modified Bessel function of the second kind of order one
@@ -16,8 +16,8 @@
namespace boost { namespace math { namespace detail{
-template <typename T>
-T bessel_k1(T x)
+template <typename T, typename Policy>
+T bessel_k1(T x, const Policy& pol)
{
static const T P1[] = {
static_cast<T>(-2.2149374878243304548e+06L),
@@ -77,14 +77,16 @@
using namespace std;
using namespace boost::math::tools;
+ static const char* function = "boost::math::bessel_k1<%1%>(%1%,%1%)";
+
if (x < 0)
{
- return domain_error<T>(BOOST_CURRENT_FUNCTION,
- "Got x = %1%, but argument x must be non-negative, complex number result not supported.", x);
+ return policy::raise_domain_error<T>(function,
+ "Got x = %1%, but argument x must be non-negative, complex number result not supported.", x, pol);
}
if (x == 0)
{
- return overflow_error<T>(BOOST_CURRENT_FUNCTION);
+ return policy::raise_overflow_error<T>(function, 0, pol);
}
if (x <= 1) // x in (0, 1]
{
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_kn.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_kn.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_kn.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -8,45 +8,47 @@
#include <boost/math/special_functions/detail/bessel_k0.hpp>
#include <boost/math/special_functions/detail/bessel_k1.hpp>
-#include <boost/math/tools/error_handling.hpp>
+#include <boost/math/policy/error_handling.hpp>
// Modified Bessel function of the second kind of integer order
// K_n(z) is the dominant solution, forward recurrence always OK (though unstable)
namespace boost { namespace math { namespace detail{
-template <typename T>
-T bessel_kn(int n, T x)
+template <typename T, typename Policy>
+T bessel_kn(int n, T x, const Policy& pol)
{
T value, current, prev;
using namespace boost::math::tools;
+ static const char* function = "boost::math::bessel_kn<%1%>(%1%,%1%)";
+
if (x < 0)
{
- return domain_error<T>(BOOST_CURRENT_FUNCTION,
- "Got x = %1%, but argument x must be non-negative, complex number result not supported.", x);
+ return policy::raise_domain_error<T>(function,
+ "Got x = %1%, but argument x must be non-negative, complex number result not supported.", x, pol);
}
if (x == 0)
{
- return overflow_error<T>(BOOST_CURRENT_FUNCTION);
+ return policy::raise_overflow_error<T>(function, 0, pol);
}
if (n == 0)
{
- return bessel_k0(x);
+ return bessel_k0(x, pol);
}
if (n == 1)
{
- return bessel_k1(x);
+ return bessel_k1(x, pol);
}
if (n < 0)
{
n = -n; // K_{-n}(z) = K_n(z)
}
- prev = bessel_k0(x);
- current = bessel_k1(x);
+ prev = bessel_k0(x, pol);
+ current = bessel_k1(x, pol);
for (int k = 1; k < n; k++) // n >= 2
{
value = 2 * k * current / x + prev;
Deleted: sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_kv.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_kv.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
+++ (empty file)
@@ -1,91 +0,0 @@
-// Copyright (c) 2006 Xiaogang Zhang
-// 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)
-
-#ifndef BOOST_MATH_BESSEL_KV_HPP
-#define BOOST_MATH_BESSEL_KV_HPP
-
-#include <boost/math/special_functions/bessel_kn.hpp>
-#include <boost/math/special_functions/bessel_ik.hpp>
-#include <boost/math/constants/constants.hpp>
-
-// Modified Bessel and Spherical Bessel functions of the second kind
-
-namespace boost { namespace math {
-
-// Modified Bessel function of the second kind of any order
-template <typename T>
-inline T bessel_kv(T v, T x)
-{
- int n = static_cast<int>(v);
- if (v == n)
- {
- return bessel_kn(n, x); // v is integer
- }
- else
- {
- T I, K;
- bessel_ik(v, x, &I, &K);
- return K;
- }
-}
-
-// Modified Spherical Bessel function of the second kind of non-negative order
-template <typename T>
-inline T bessel_skv(unsigned n, T x)
-{
- using namespace std;
- using namespace boost::math::tools;
- using namespace boost::math::constants;
-
- T v = n + 0.5L;
- if (x == 0)
- {
- return overflow_error<T>("boost::math::bessel_skv(n, x)",
- "infinity occurred but not supported");
- }
- else
- {
- return sqrt(0.5L * pi<T>() / x) * bessel_kv(v, x);
- }
-}
-
-// -------------------- TR1 functions --------------------
-
-inline float cyl_bessel_kf(float nu, float x)
-{
- if (nu >= 128)
- {
- std::cout << "Warning: cyl_bessel_kf(nu, x), nu >= 128, "
- << "result is implementation defined according to C++ TR1"
- << std::endl;
- }
- return bessel_kv<float>(nu, x);
-}
-
-inline double cyl_bessel_k(double nu, double x)
-{
- if (nu >= 128)
- {
- std::cout << "Warning: cyl_bessel_k(nu, x), nu >= 128, "
- << "result is implementation defined according to C++ TR1"
- << std::endl;
- }
- return bessel_kv<double>(nu, x);
-}
-
-inline long double cyl_bessel_kl(long double nu, long double x)
-{
- if (nu >= 128)
- {
- std::cout << "Warning: cyl_bessel_kl(nu, x), nu >= 128, "
- << "result is implementation defined according to C++ TR1"
- << std::endl;
- }
- return bessel_kv<long double>(nu, x);
-}
-
-}} // namespaces
-
-#endif // BOOST_MATH_BESSEL_KV_HPP
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_y0.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_y0.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_y0.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -9,7 +9,7 @@
#include <boost/math/special_functions/detail/bessel_j0.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/rational.hpp>
-#include <boost/math/tools/error_handling.hpp>
+#include <boost/math/policy/error_handling.hpp>
#include <boost/assert.hpp>
// Bessel function of the second kind of order zero
@@ -18,8 +18,8 @@
namespace boost { namespace math { namespace detail{
-template <typename T>
-T bessel_y0(T x)
+template <typename T, typename Policy>
+T bessel_y0(T x, const Policy& pol)
{
static const T P1[] = {
static_cast<T>(1.0723538782003176831e+11L),
@@ -123,14 +123,16 @@
using namespace boost::math::tools;
using namespace boost::math::constants;
+ static const char* function = "boost::math::bessel_y0<%1%>(%1%,%1%)";
+
if (x < 0)
{
- return domain_error<T>(BOOST_CURRENT_FUNCTION,
- "Got x = %1% but x must be non-negative, complex result not supported.", x);
+ return policy::raise_domain_error<T>(function,
+ "Got x = %1% but x must be non-negative, complex result not supported.", x, pol);
}
if (x == 0)
{
- return -overflow_error<T>(BOOST_CURRENT_FUNCTION);
+ return -policy::raise_overflow_error<T>(function, 0, pol);
}
if (x <= 3) // x in (0, 3]
{
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_y1.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_y1.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_y1.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -9,7 +9,7 @@
#include <boost/math/special_functions/detail/bessel_j1.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/rational.hpp>
-#include <boost/math/tools/error_handling.hpp>
+#include <boost/math/policy/error_handling.hpp>
#include <boost/assert.hpp>
// Bessel function of the second kind of order one
@@ -18,8 +18,8 @@
namespace boost { namespace math { namespace detail{
-template <typename T>
-T bessel_y1(T x)
+template <typename T, typename Policy>
+T bessel_y1(T x, const Policy& pol)
{
static const T P1[] = {
static_cast<T>(4.0535726612579544093e+13L),
@@ -112,8 +112,8 @@
if (x <= 0)
{
- return domain_error<T>(BOOST_CURRENT_FUNCTION,
- "Got x == %1%, but x must be > 0, complex result not supported.", x);
+ return policy::raise_domain_error<T>("bost::math::bessel_y1<%1%>(%1%,%1%)",
+ "Got x == %1%, but x must be > 0, complex result not supported.", x, pol);
}
if (x <= 4) // x in (0, 4]
{
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_yn.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_yn.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_yn.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -8,37 +8,39 @@
#include <boost/math/special_functions/detail/bessel_y0.hpp>
#include <boost/math/special_functions/detail/bessel_y1.hpp>
-#include <boost/math/tools/error_handling.hpp>
+#include <boost/math/policy/error_handling.hpp>
// Bessel function of the second kind of integer order
// Y_n(z) is the dominant solution, forward recurrence always OK (though unstable)
namespace boost { namespace math { namespace detail{
-template <typename T>
-T bessel_yn(int n, T x)
+template <typename T, typename Policy>
+T bessel_yn(int n, T x, const Policy& pol)
{
T value, factor, current, prev;
using namespace boost::math::tools;
+ static const char* function = "boost::math::bessel_yn<%1%>(%1%,%1%)";
+
if ((x == 0) && (n == 0))
{
- return -overflow_error<T>(BOOST_CURRENT_FUNCTION);
+ return -policy::raise_overflow_error<T>(function, 0, pol);
}
if (x <= 0)
{
- return domain_error<T>(BOOST_CURRENT_FUNCTION,
- "Got x = %1%, but x must be > 0, complex result not supported.", x);
+ return policy::raise_domain_error<T>(function,
+ "Got x = %1%, but x must be > 0, complex result not supported.", x, pol);
}
if (n == 0)
{
- return bessel_y0(x);
+ return bessel_y0(x, pol);
}
if (n == 1)
{
- return bessel_y1(x);
+ return bessel_y1(x, pol);
}
if (n < 0)
{
@@ -50,8 +52,8 @@
factor = 1;
}
- prev = bessel_y0(x);
- current = bessel_y1(x);
+ prev = bessel_y0(x, pol);
+ current = bessel_y1(x, pol);
for (int k = 1; k < n; k++) // n >= 2
{
value = 2 * k * current / x - prev;
Deleted: sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_yv.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/detail/bessel_yv.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
+++ (empty file)
@@ -1,125 +0,0 @@
-// Copyright (c) 2006 Xiaogang Zhang
-// 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)
-
-#ifndef BOOST_MATH_BESSEL_YV_HPP
-#define BOOST_MATH_BESSEL_YV_HPP
-
-#include <boost/math/special_functions/bessel_yn.hpp>
-#include <boost/math/special_functions/bessel_jy.hpp>
-#include <boost/math/constants/constants.hpp>
-#include <boost/math/tools/error_handling.hpp>
-
-// Bessel and Spherical Bessel functions of the second kind
-
-namespace boost { namespace math {
-
-// Bessel function of the second kind of any order
-template <typename T>
-inline T bessel_yv(T v, T x)
-{
- int n = static_cast<int>(v);
- if (v == n)
- {
- return bessel_yn(n, x); // v is integer
- }
- else
- {
- T J, Y;
- bessel_jy(v, x, &J, &Y);
- return Y;
- }
-}
-
-// Spherical Bessel function of the second kind of non-negative order
-template <typename T>
-inline T bessel_syv(unsigned n, T x)
-{
- using namespace std;
- using namespace boost::math::tools;
- using namespace boost::math::constants;
-
- T v = n + 0.5L;
- if (x == 0)
- {
- return -overflow_error<T>("boost::math::bessel_syv(n, x)",
- "infinity occurred but not supported");
- }
- else
- {
- return sqrt(0.5L * pi<T>() / x) * bessel_yv(v, x);
- }
-}
-
-// -------------------- TR1 functions --------------------
-
-inline float cyl_neumannf(float nu, float x)
-{
- if (nu >= 128)
- {
- std::cout << "Warning: cyl_neumannf(nu, x), nu >= 128, "
- << "result is implementation defined according to C++ TR1"
- << std::endl;
- }
- return bessel_yv<float>(nu, x);
-}
-
-inline double cyl_neumann(double nu, double x)
-{
- if (nu >= 128)
- {
- std::cout << "Warning: cyl_neumann(nu, x), nu >= 128, "
- << "result is implementation defined according to C++ TR1"
- << std::endl;
- }
- return bessel_yv<double>(nu, x);
-}
-
-inline long double cyl_neumannl(long double nu, long double x)
-{
- if (nu >= 128)
- {
- std::cout << "Warning: cyl_neumannl(nu, x), nu >= 128, "
- << "result is implementation defined according to C++ TR1"
- << std::endl;
- }
- return bessel_yv<long double>(nu, x);
-}
-
-inline float sph_neumannf(unsigned n, float x)
-{
- if (n >= 128)
- {
- std::cout << "Warning: sph_neumannf(n, x), n >= 128, "
- << "result is implementation defined according to C++ TR1"
- << std::endl;
- }
- return bessel_syv<float>(n, x);
-}
-
-inline double sph_neumann(unsigned n, double x)
-{
- if (n >= 128)
- {
- std::cout << "Warning: sph_neumann(n, x), n >= 128, "
- << "result is implementation defined according to C++ TR1"
- << std::endl;
- }
- return bessel_syv<double>(n, x);
-}
-
-inline long double sph_neumannl(unsigned n, long double x)
-{
- if (n >= 128)
- {
- std::cout << "Warning: sph_neumannl(n, x), n >= 128, "
- << "result is implementation defined according to C++ TR1"
- << std::endl;
- }
- return bessel_syv<long double>(n, x);
-}
-
-}} // namespaces
-
-#endif // BOOST_MATH_BESSEL_YV_HPP
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/factorials.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/factorials.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/factorials.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -23,19 +23,26 @@
namespace boost { namespace math
{
-template <class T>
-inline T factorial(unsigned i)
+template <class T, class Policy>
+inline T factorial(unsigned i, const Policy& pol)
{
using namespace std; // Aid ADL for floor.
if(i <= max_factorial<T>::value)
return unchecked_factorial<T>(i);
- T result = boost::math::tgamma(static_cast<T>(i+1));
+ T result = boost::math::tgamma(static_cast<T>(i+1), pol);
if(result > tools::max_value<T>())
return result; // Overflowed value! (But tgamma will have signalled the error already).
return floor(result + 0.5);
}
+template <class T>
+inline T factorial(unsigned i)
+{
+ return factorial<T>(i, policy::policy<>());
+}
+/*
+// Can't have these in a policy enabled world?
template<>
inline float factorial<float>(unsigned i)
{
@@ -51,9 +58,9 @@
return unchecked_factorial<double>(i);
return tools::overflow_error<double>(BOOST_CURRENT_FUNCTION);
}
-
-template <class T>
-T double_factorial(unsigned i)
+*/
+template <class T, class Policy>
+T double_factorial(unsigned i, const Policy& pol)
{
using namespace std; // ADL lookup of std names
if(i & 1)
@@ -68,7 +75,7 @@
// Fallthrough: i is too large to use table lookup, try the
// gamma function instead.
//
- T result = boost::math::tgamma(static_cast<T>(i) / 2 + 1) / sqrt(constants::pi<T>());
+ T result = boost::math::tgamma(static_cast<T>(i) / 2 + 1, pol) / sqrt(constants::pi<T>());
if(ldexp(tools::max_value<T>(), -static_cast<int>(i+1) / 2) > result)
return ceil(result * ldexp(T(1), (i+1) / 2) - 0.5f);
}
@@ -76,20 +83,26 @@
{
// even i:
unsigned n = i / 2;
- T result = factorial<T>(n);
+ T result = factorial<T>(n, pol);
if(ldexp(tools::max_value<T>(), -(int)n) > result)
return result * ldexp(T(1), (int)n);
}
//
// If we fall through to here then the result is infinite:
//
- return tools::overflow_error<T>(BOOST_CURRENT_FUNCTION);
+ return policy::raise_overflow_error<T>("boost::math::double_factorial<%1%>(unsigned)", 0, pol);
+}
+
+template <class T>
+inline T double_factorial(unsigned i)
+{
+ return double_factorial<T>(i, policy::policy<>());
}
namespace detail{
-template <class T>
-T rising_factorial_imp(T x, int n)
+template <class T, class Policy>
+T rising_factorial_imp(T x, int n, const Policy& pol)
{
if(x < 0)
{
@@ -108,7 +121,7 @@
n = -n;
inv = true;
}
- T result = ((n&1) ? -1 : 1) * falling_factorial(-x, n);
+ T result = ((n&1) ? -1 : 1) * falling_factorial(-x, n, pol);
if(inv)
result = 1 / result;
return result;
@@ -120,11 +133,11 @@
// tgamma_delta_ratio is alreay optimised for that
// use case:
//
- return 1 / tgamma_delta_ratio(x, static_cast<T>(n));
+ return 1 / tgamma_delta_ratio(x, static_cast<T>(n), pol);
}
-template <class T>
-inline T falling_factorial_imp(T x, unsigned n)
+template <class T, class Policy>
+inline T falling_factorial_imp(T x, unsigned n, const Policy& pol)
{
using namespace std; // ADL of std names
if(x == 0)
@@ -135,7 +148,7 @@
// For x < 0 we really have a rising factorial
// modulo a possible change of sign:
//
- return (n&1 ? -1 : 1) * rising_factorial(-x, n);
+ return (n&1 ? -1 : 1) * rising_factorial(-x, n, pol);
}
if(n == 0)
return 1;
@@ -149,12 +162,12 @@
unsigned n2 = tools::real_cast<unsigned>(floor(xp1));
if(n2 == xp1)
return 0;
- T result = tgamma_delta_ratio(xp1, -static_cast<T>(n2));
+ T result = tgamma_delta_ratio(xp1, -static_cast<T>(n2), pol);
x -= n2;
result *= x;
++n2;
if(n2 < n)
- result *= falling_factorial(x - 1, n - n2);
+ result *= falling_factorial(x - 1, n - n2, pol);
return result;
}
//
@@ -164,7 +177,7 @@
// because tgamma_delta_ratio is alreay optimised
// for that use case:
//
- return tgamma_delta_ratio(x + 1, -static_cast<T>(n));
+ return tgamma_delta_ratio(x + 1, -static_cast<T>(n), pol);
}
} // namespace detail
@@ -175,7 +188,16 @@
{
typedef typename tools::promote_args<RT>::type result_type;
return detail::falling_factorial_imp(
- static_cast<result_type>(x), n);
+ static_cast<result_type>(x), n, policy::policy<>());
+}
+
+template <class RT, class Policy>
+inline typename tools::promote_args<RT>::type
+ falling_factorial(RT x, unsigned n, const Policy& pol)
+{
+ typedef typename tools::promote_args<RT>::type result_type;
+ return detail::falling_factorial_imp(
+ static_cast<result_type>(x), n, pol);
}
template <class RT>
@@ -184,7 +206,16 @@
{
typedef typename tools::promote_args<RT>::type result_type;
return detail::rising_factorial_imp(
- static_cast<result_type>(x), n);
+ static_cast<result_type>(x), n, policy::policy<>());
+}
+
+template <class RT, class Policy>
+inline typename tools::promote_args<RT>::type
+ rising_factorial(RT x, unsigned n, const Policy& pol)
+{
+ typedef typename tools::promote_args<RT>::type result_type;
+ return detail::rising_factorial_imp(
+ static_cast<result_type>(x), n, pol);
}
} // namespace math
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/hypot.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/hypot.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/hypot.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -8,6 +8,7 @@
#include <boost/math/tools/config.hpp>
#include <boost/math/tools/precision.hpp>
+#include <boost/math/policy/error_handling.hpp>
#include <boost/math/special_functions/math_fwd.hpp>
#include <cmath>
#include <algorithm> // for swap
@@ -18,8 +19,8 @@
namespace boost{ namespace math{ namespace detail{
-template <class T>
-T hypot_imp(T x, T y)
+template <class T, class Policy>
+T hypot_imp(T x, T y, const Policy& pol)
{
//
// Normalize x and y, so that both are positive and x >= y:
@@ -37,7 +38,7 @@
if(std::numeric_limits<T>::has_infinity
&& ((x == std::numeric_limits<T>::infinity())
|| (y == std::numeric_limits<T>::infinity())))
- return std::numeric_limits<T>::infinity();
+ return policy::raise_overflow_error<T>("boost::math::hypot<%1%>(%1%,%1%)", 0, pol);
#ifdef BOOST_MSVC
#pragma warning(pop)
#endif
@@ -60,7 +61,16 @@
{
typedef typename tools::promote_args<T1, T2>::type result_type;
return detail::hypot_imp(
- static_cast<result_type>(x), static_cast<result_type>(y));
+ static_cast<result_type>(x), static_cast<result_type>(y), policy::policy<>());
+}
+
+template <class T1, class T2, class Policy>
+inline typename tools::promote_args<T1, T2>::type
+ hypot(T1 x, T2 y, const Policy& pol)
+{
+ typedef typename tools::promote_args<T1, T2>::type result_type;
+ return detail::hypot_imp(
+ static_cast<result_type>(x), static_cast<result_type>(y), pol);
}
} // namespace math
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/sin_pi.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/sin_pi.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/sin_pi.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -40,6 +40,12 @@
return invert ? -rem : rem;
}
+template <class T, class Policy>
+inline T sin_pi(T x, const Policy&)
+{
+ return boost::math::sin_pi(x);
+}
+
} // namespace math
} // namespace boost
#endif
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/sinc.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/sinc.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/sinc.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -13,6 +13,7 @@
#include <boost/math/tools/config.hpp>
#include <boost/math/tools/precision.hpp>
+#include <boost/math/policy/policy.hpp>
#include <boost/math/special_functions/math_fwd.hpp>
#include <cmath>
#include <boost/limits.hpp>
@@ -98,6 +99,13 @@
return detail::sinc_pi_imp(static_cast<result_type>(x));
}
+ template <class T, class Policy>
+ inline typename tools::promote_args<T>::type sinc_pi(T x, const Policy&)
+ {
+ typedef typename tools::promote_args<T>::type result_type;
+ return detail::sinc_pi_imp(static_cast<result_type>(x));
+ }
+
#ifdef BOOST_NO_TEMPLATE_TEMPLATES
#else /* BOOST_NO_TEMPLATE_TEMPLATES */
template<typename T, template<typename> class U>
@@ -151,6 +159,12 @@
return(result);
}
}
+
+ template<typename T, template<typename> class U, class Policy>
+ inline U<T> sinc_pi(const U<T> x, const Policy&)
+ {
+ return sinc_pi(x);
+ }
#endif /* BOOST_NO_TEMPLATE_TEMPLATES */
}
}
Modified: sandbox/math_toolkit/policy/boost/math/tools/test.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/tools/test.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/tools/test.hpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -177,9 +177,21 @@
{
const row_type& row = a[i];
value_type point;
- try{
+ try
+ {
point = test_func(row);
- }catch(const std::exception& e)
+ }
+ catch(const std::underflow_error&)
+ {
+ point = 0;
+ }
+ catch(const std::overflow_error&)
+ {
+ point = std::numeric_limits<value_type>::has_infinity ?
+ std::numeric_limits<value_type>::infinity()
+ : tools::max_value<value_type>();
+ }
+ catch(const std::exception& e)
{
std::cerr << e.what() << std::endl;
print_row(row);
Modified: sandbox/math_toolkit/policy/libs/math/test/hypot_test.cpp
==============================================================================
--- sandbox/math_toolkit/policy/libs/math/test/hypot_test.cpp (original)
+++ sandbox/math_toolkit/policy/libs/math/test/hypot_test.cpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -3,6 +3,7 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/special_functions/hypot.hpp>
Modified: sandbox/math_toolkit/policy/libs/math/test/test_bessel_i.cpp
==============================================================================
--- sandbox/math_toolkit/policy/libs/math/test/test_bessel_i.cpp (original)
+++ sandbox/math_toolkit/policy/libs/math/test/test_bessel_i.cpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -3,6 +3,7 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
#include <boost/math/concepts/real_concept.hpp>
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
Modified: sandbox/math_toolkit/policy/libs/math/test/test_bessel_k.cpp
==============================================================================
--- sandbox/math_toolkit/policy/libs/math/test/test_bessel_k.cpp (original)
+++ sandbox/math_toolkit/policy/libs/math/test/test_bessel_k.cpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -3,6 +3,7 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
#include <boost/math/concepts/real_concept.hpp>
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
Modified: sandbox/math_toolkit/policy/libs/math/test/test_bessel_y.cpp
==============================================================================
--- sandbox/math_toolkit/policy/libs/math/test/test_bessel_y.cpp (original)
+++ sandbox/math_toolkit/policy/libs/math/test/test_bessel_y.cpp 2007-07-13 08:48:57 EDT (Fri, 13 Jul 2007)
@@ -3,6 +3,7 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
#include <boost/math/concepts/real_concept.hpp>
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk