|
|
Boost-Commit : |
From: john_at_[hidden]
Date: 2008-01-23 05:41:08
Author: johnmaddock
Date: 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
New Revision: 42920
URL: http://svn.boost.org/trac/boost/changeset/42920
Log:
Added error handling to the rounding functions.
Added better error handling to the non-central chi squared, and updated the tests.
Added:
sandbox/math_toolkit/libs/math/test/compile_test/dist_nc_chi_squ_incl_test.cpp (contents, props changed)
Text files modified:
sandbox/math_toolkit/boost/math/concepts/real_concept.hpp | 36 +++++++++++-----------
sandbox/math_toolkit/boost/math/concepts/std_real_concept.hpp | 36 +++++++++++-----------
sandbox/math_toolkit/boost/math/distributions/non_central_chi_squared.hpp | 39 ++++++++++++++++++++----
sandbox/math_toolkit/boost/math/distributions/poisson.hpp | 2
sandbox/math_toolkit/boost/math/policies/error_handling.hpp | 64 ++++++++++++++++++++++++++++++++++++++++
sandbox/math_toolkit/boost/math/policies/policy.hpp | 9 +++++
sandbox/math_toolkit/boost/math/special_functions/bessel.hpp | 10 +++---
sandbox/math_toolkit/boost/math/special_functions/beta.hpp | 4 +-
sandbox/math_toolkit/boost/math/special_functions/cos_pi.hpp | 12 +++---
sandbox/math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp | 2
sandbox/math_toolkit/boost/math/special_functions/detail/bessel_jy.hpp | 2
sandbox/math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp | 4 +-
sandbox/math_toolkit/boost/math/special_functions/factorials.hpp | 2
sandbox/math_toolkit/boost/math/special_functions/gamma.hpp | 6 +-
sandbox/math_toolkit/boost/math/special_functions/modf.hpp | 36 +++++++++++++++++-----
sandbox/math_toolkit/boost/math/special_functions/round.hpp | 43 +++++++++++++++++++++++---
sandbox/math_toolkit/boost/math/special_functions/sin_pi.hpp | 12 +++---
sandbox/math_toolkit/boost/math/special_functions/trunc.hpp | 43 +++++++++++++++++++++++---
sandbox/math_toolkit/libs/math/test/Jamfile.v2 | 1
sandbox/math_toolkit/libs/math/test/test_nc_chi_squared.cpp | 44 +++++++++++++-------------
sandbox/math_toolkit/libs/math/test/test_round.cpp | 61 ++++++++++++++++++++++++++++++++++++++
21 files changed, 357 insertions(+), 111 deletions(-)
Modified: sandbox/math_toolkit/boost/math/concepts/real_concept.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/concepts/real_concept.hpp (original)
+++ sandbox/math_toolkit/boost/math/concepts/real_concept.hpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -313,46 +313,46 @@
//
// Conversion and truncation routines:
//
-template <>
-inline int iround<concepts::real_concept>(const concepts::real_concept& v)
+template <class Policy>
+inline int iround(const concepts::real_concept& v, const Policy& pol)
{
- return iround(v.value());
+ return iround(v.value(), pol);
}
-template <>
-inline long lround<concepts::real_concept>(const concepts::real_concept& v)
+template <class Policy>
+inline long lround(const concepts::real_concept& v, const Policy& pol)
{
- return lround(v.value());
+ return lround(v.value(), pol);
}
#ifdef BOOST_HAS_LONG_LONG
-template <>
-inline long long llround<concepts::real_concept>(const concepts::real_concept& v)
+template <class Policy>
+inline long long llround(const concepts::real_concept& v, const Policy& pol)
{
- return llround(v.value());
+ return llround(v.value(), pol);
}
#endif
-template <>
-inline int itrunc<concepts::real_concept>(const concepts::real_concept& v)
+template <class Policy>
+inline int itrunc(const concepts::real_concept& v, const Policy& pol)
{
- return itrunc(v.value());
+ return itrunc(v.value(), pol);
}
-template <>
-inline long ltrunc<concepts::real_concept>(const concepts::real_concept& v)
+template <class Policy>
+inline long ltrunc(const concepts::real_concept& v, const Policy& pol)
{
- return ltrunc(v.value());
+ return ltrunc(v.value(), pol);
}
#ifdef BOOST_HAS_LONG_LONG
-template <>
-inline long long lltrunc<concepts::real_concept>(const concepts::real_concept& v)
+template <class Policy>
+inline long long lltrunc(const concepts::real_concept& v, const Policy& pol)
{
- return lltrunc(v.value());
+ return lltrunc(v.value(), pol);
}
#endif
Modified: sandbox/math_toolkit/boost/math/concepts/std_real_concept.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/concepts/std_real_concept.hpp (original)
+++ sandbox/math_toolkit/boost/math/concepts/std_real_concept.hpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -282,46 +282,46 @@
//
// Conversion and truncation routines:
//
-template <>
-inline int iround<concepts::std_real_concept>(const concepts::std_real_concept& v)
+template <class Policy>
+inline int iround(const concepts::std_real_concept& v, const Policy& pol)
{
- return iround(v.value());
+ return iround(v.value(), pol);
}
-template <>
-inline long lround<concepts::std_real_concept>(const concepts::std_real_concept& v)
+template <class Policy>
+inline long lround(const concepts::std_real_concept& v, const Policy& pol)
{
- return lround(v.value());
+ return lround(v.value(), pol);
}
#ifdef BOOST_HAS_LONG_LONG
-template <>
-inline long long llround<concepts::std_real_concept>(const concepts::std_real_concept& v)
+template <class Policy>
+inline long long llround(const concepts::std_real_concept& v, const Policy& pol)
{
- return llround(v.value());
+ return llround(v.value(), pol);
}
#endif
-template <>
-inline int itrunc<concepts::std_real_concept>(const concepts::std_real_concept& v)
+template <class Policy>
+inline int itrunc(const concepts::std_real_concept& v, const Policy& pol)
{
- return itrunc(v.value());
+ return itrunc(v.value(), pol);
}
-template <>
-inline long ltrunc<concepts::std_real_concept>(const concepts::std_real_concept& v)
+template <class Policy>
+inline long ltrunc(const concepts::std_real_concept& v, const Policy& pol)
{
- return ltrunc(v.value());
+ return ltrunc(v.value(), pol);
}
#ifdef BOOST_HAS_LONG_LONG
-template <>
-inline long long lltrunc<concepts::std_real_concept>(const concepts::std_real_concept& v)
+template <class Policy>
+inline long long lltrunc(const concepts::std_real_concept& v, const Policy& pol)
{
- return lltrunc(v.value());
+ return lltrunc(v.value(), pol);
}
#endif
Modified: sandbox/math_toolkit/boost/math/distributions/non_central_chi_squared.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/non_central_chi_squared.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/non_central_chi_squared.hpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -70,7 +70,7 @@
// k is chosen as the peek of the Poisson weights, which
// will occur *before* the largest term.
//
- int k = iround(lambda);
+ int k = iround(lambda, pol);
// Forwards and backwards Poisson weights:
T poisf = boost::math::gamma_p_derivative(1 + k, lambda, pol);
T poisb = poisf * k / lambda;
@@ -87,7 +87,8 @@
// stable direction for the gamma function
// recurrences:
//
- for(int i = k; static_cast<boost::uintmax_t>(i) < max_iter; ++i)
+ int i;
+ for(i = k; static_cast<boost::uintmax_t>(i-k) < max_iter; ++i)
{
T term = poisf * gamf;
sum += term;
@@ -97,6 +98,11 @@
if((sum == 0) || (term / sum < errtol))
break;
}
+ //Error check:
+ if(static_cast<boost::uintmax_t>(i-k) >= max_iter)
+ policies::raise_evaluation_error(
+ "cdf(non_central_chi_squared_distribution<%1%>, %1%)",
+ "Series did not converge, closest value was %1%", sum, pol);
//
// Now backwards iteration: the gamma
// function recurrences are unstable in this
@@ -153,7 +159,8 @@
boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
T errtol = ldexp(1.0, -boost::math::policies::digits<T, Policy>());
- for(int i = 1; static_cast<boost::uintmax_t>(i) < max_iter; ++i)
+ int i;
+ for(i = 1; static_cast<boost::uintmax_t>(i) < max_iter; ++i)
{
tk = tk * x / (f + 2 * i);
uk = uk * lambda / i;
@@ -163,6 +170,11 @@
if(term / sum < errtol)
break;
}
+ //Error check:
+ if(static_cast<boost::uintmax_t>(i) >= max_iter)
+ policies::raise_evaluation_error(
+ "cdf(non_central_chi_squared_distribution<%1%>, %1%)",
+ "Series did not converge, closest value was %1%", sum, pol);
return sum;
}
@@ -185,7 +197,7 @@
//
BOOST_MATH_STD_USING
// Special case:
- if(x == 0)
+ if(y == 0)
return 0;
boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
T errtol = ldexp(1.0, -boost::math::policies::digits<T, Policy>());
@@ -200,7 +212,7 @@
// function, which ocurrs *after* the largest term in the
// sum.
//
- int k = iround(del);
+ int k = iround(del, pol);
T a = n / 2 + k;
// Central chi squared term for forward iteration:
T gamkf = boost::math::gamma_p(a, x, pol);
@@ -256,6 +268,12 @@
++i;
}while((errorf / sum > errtol) && (static_cast<boost::uintmax_t>(i) < max_iter));
+ //Error check:
+ if(static_cast<boost::uintmax_t>(i) >= max_iter)
+ policies::raise_evaluation_error(
+ "cdf(non_central_chi_squared_distribution<%1%>, %1%)",
+ "Series did not converge, closest value was %1%", sum, pol);
+
return sum;
}
@@ -511,14 +529,21 @@
v = pdf(dist, lower_bound);
}while(maxval < v);
- boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
+ boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
- return tools::brent_find_minima(
+ RealType result = tools::brent_find_minima(
pdf_minimizer<RealType, Policy>(dist),
lower_bound,
upper_bound,
policies::digits<RealType, Policy>(),
max_iter).first;
+ if(max_iter == policies::get_max_root_iterations<Policy>())
+ {
+ return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
+ " either there is no answer to the mode of the non central chi squared distribution"
+ " or the answer is infinite. Current best guess is %1%", result, Policy());
+ }
+ return result;
}
template <class RealType, class Policy>
Modified: sandbox/math_toolkit/boost/math/distributions/poisson.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/poisson.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/poisson.hpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -341,7 +341,7 @@
&& k < max_factorial<RealType>::value)
{ // k is small enough (for float 34, double 170 ...) to use factorial(k).
return exp(-mean) * pow(mean, k) /
- unchecked_factorial<RealType>(itrunc(floork));
+ unchecked_factorial<RealType>(itrunc(floork, Policy()));
}
else
{ // Need to use log(factorial(k)) = lgamma(k+1)
Modified: sandbox/math_toolkit/boost/math/policies/error_handling.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/policies/error_handling.hpp (original)
+++ sandbox/math_toolkit/boost/math/policies/error_handling.hpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -37,6 +37,12 @@
evaluation_error(const std::string& s) : std::runtime_error(s){}
};
+class rounding_error : public std::runtime_error
+{
+public:
+ rounding_error(const std::string& s) : std::runtime_error(s){}
+};
+
namespace policies{
//
// Forward declarations of user error handlers,
@@ -54,6 +60,8 @@
T user_denorm_error(const char* function, const char* message, const T& val);
template <class T>
T user_evaluation_error(const char* function, const char* message, const T& val);
+template <class T>
+T user_rounding_error(const char* function, const char* message, const T& val);
namespace detail
{
@@ -362,6 +370,53 @@
return user_evaluation_error(function, message, val);
}
+template <class T>
+inline T raise_rounding_error(
+ const char* function,
+ const char* message,
+ const T& val,
+ const ::boost::math::policies::rounding_error< ::boost::math::policies::throw_on_error>&)
+{
+ raise_error<boost::math::rounding_error, T>(function, message, val);
+ // we never get here:
+ return T(0);
+}
+
+template <class T>
+inline T raise_rounding_error(
+ const char* ,
+ const char* ,
+ const T& val,
+ const ::boost::math::policies::rounding_error< ::boost::math::policies::ignore_error>&)
+{
+ // This may or may not do the right thing, but the user asked for the error
+ // to be ignored so here we go anyway:
+ return val;
+}
+
+template <class T>
+inline T raise_rounding_error(
+ const char* ,
+ const char* ,
+ const T& val,
+ const ::boost::math::policies::rounding_error< ::boost::math::policies::errno_on_error>&)
+{
+ errno = ERANGE;
+ // This may or may not do the right thing, but the user asked for the error
+ // to be silent so here we go anyway:
+ return val;
+}
+
+template <class T>
+inline T raise_rounding_error(
+ const char* function,
+ const char* message,
+ const T& val,
+ const ::boost::math::policies::rounding_error< ::boost::math::policies::user_error>&)
+{
+ return user_rounding_error(function, message, val);
+}
+
} // namespace detail
template <class T, class Policy>
@@ -419,6 +474,15 @@
val, policy_type());
}
+template <class T, class Policy>
+inline T raise_rounding_error(const char* function, const char* message, const T& val, const Policy&)
+{
+ typedef typename Policy::rounding_error_type policy_type;
+ return detail::raise_rounding_error(
+ function, message ? message : "Value %1% can not be represented in the target integer type.",
+ val, policy_type());
+}
+
//
// checked_narrowing_cast:
//
Modified: sandbox/math_toolkit/boost/math/policies/policy.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/policies/policy.hpp (original)
+++ sandbox/math_toolkit/boost/math/policies/policy.hpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -54,6 +54,9 @@
#ifndef BOOST_MATH_EVALUATION_ERROR_POLICY
#define BOOST_MATH_EVALUATION_ERROR_POLICY throw_on_error
#endif
+#ifndef BOOST_MATH_ROUNDING_ERROR_POLICY
+#define BOOST_MATH_ROUNDING_ERROR_POLICY throw_on_error
+#endif
#ifndef BOOST_MATH_UNDERFLOW_ERROR_POLICY
#define BOOST_MATH_UNDERFLOW_ERROR_POLICY ignore_error
#endif
@@ -178,6 +181,7 @@
BOOST_MATH_META_INT(error_policy_type, underflow_error, BOOST_MATH_UNDERFLOW_ERROR_POLICY)
BOOST_MATH_META_INT(error_policy_type, denorm_error, BOOST_MATH_DENORM_ERROR_POLICY)
BOOST_MATH_META_INT(error_policy_type, evaluation_error, BOOST_MATH_EVALUATION_ERROR_POLICY)
+BOOST_MATH_META_INT(error_policy_type, rounding_error, BOOST_MATH_ROUNDING_ERROR_POLICY)
//
// Policy types for internal promotion:
@@ -398,6 +402,7 @@
typedef typename detail::find_arg<arg_list, is_underflow_error<mpl::_1>, underflow_error<> >::type underflow_error_type;
typedef typename detail::find_arg<arg_list, is_denorm_error<mpl::_1>, denorm_error<> >::type denorm_error_type;
typedef typename detail::find_arg<arg_list, is_evaluation_error<mpl::_1>, evaluation_error<> >::type evaluation_error_type;
+ typedef typename detail::find_arg<arg_list, is_rounding_error<mpl::_1>, rounding_error<> >::type rounding_error_type;
private:
//
// Now work out the precision:
@@ -440,6 +445,7 @@
typedef underflow_error<> underflow_error_type;
typedef denorm_error<> denorm_error_type;
typedef evaluation_error<> evaluation_error_type;
+ typedef rounding_error<> rounding_error_type;
#if BOOST_MATH_DIGITS10_POLICY == 0
typedef digits2<> precision_type;
#else
@@ -463,6 +469,7 @@
typedef underflow_error<> underflow_error_type;
typedef denorm_error<> denorm_error_type;
typedef evaluation_error<> evaluation_error_type;
+ typedef rounding_error<> rounding_error_type;
#if BOOST_MATH_DIGITS10_POLICY == 0
typedef digits2<> precision_type;
#else
@@ -500,6 +507,7 @@
typedef typename detail::find_arg<arg_list, is_underflow_error<mpl::_1>, typename Policy::underflow_error_type >::type underflow_error_type;
typedef typename detail::find_arg<arg_list, is_denorm_error<mpl::_1>, typename Policy::denorm_error_type >::type denorm_error_type;
typedef typename detail::find_arg<arg_list, is_evaluation_error<mpl::_1>, typename Policy::evaluation_error_type >::type evaluation_error_type;
+ typedef typename detail::find_arg<arg_list, is_rounding_error<mpl::_1>, typename Policy::rounding_error_type >::type rounding_error_type;
//
// Now work out the precision:
//
@@ -534,6 +542,7 @@
underflow_error_type,
denorm_error_type,
evaluation_error_type,
+ rounding_error_type,
precision_type,
promote_float_type,
promote_double_type,
Modified: sandbox/math_toolkit/boost/math/special_functions/bessel.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/bessel.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/bessel.hpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -129,7 +129,7 @@
if(floor(v) == v)
{
T r = cyl_bessel_j_imp(v, -x, t, pol);
- if(iround(v) & 1)
+ if(iround(v, pol) & 1)
r = -r;
return r;
}
@@ -165,7 +165,7 @@
if(fabs(x) > asymptotic_bessel_j_limit<T>(v, tag_type()))
return asymptotic_bessel_j_large_x_2(v, x);
else
- return bessel_jn(iround(v), x, pol);
+ return bessel_jn(iround(v, pol), x, pol);
}
return cyl_bessel_j_imp(v, x, bessel_no_int_tag(), pol);
}
@@ -228,7 +228,7 @@
if(floor(v) == v)
{
T r = cyl_bessel_i_imp(v, -x, pol);
- if(iround(v) & 1)
+ if(iround(v, pol) & 1)
r = -r;
return r;
}
@@ -337,12 +337,12 @@
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);
- if((v < 0) && (itrunc(v) & 1))
+ if((v < 0) && (itrunc(v, pol) & 1))
r = -r;
return r;
}
else
- return bessel_yn(itrunc(v), x, pol);
+ return bessel_yn(itrunc(v, pol), x, pol);
}
return cyl_neumann_imp<T>(v, x, bessel_no_int_tag(), pol);
}
Modified: sandbox/math_toolkit/boost/math/special_functions/beta.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/beta.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/beta.hpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -1060,7 +1060,7 @@
else if(a > 15)
{
// sidestep so we can use the series representation:
- int n = itrunc(floor(b));
+ int n = itrunc(floor(b), pol);
if(n == b)
--n;
T bbar = b - n;
@@ -1082,7 +1082,7 @@
// the formula here for the non-normalised case is tricky to figure
// out (for me!!), and requires two pochhammer calculations rather
// than one, so leave it for now....
- int n = itrunc(floor(b));
+ int n = itrunc(floor(b), pol);
T bbar = b - n;
if(bbar <= 0)
{
Modified: sandbox/math_toolkit/boost/math/special_functions/cos_pi.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/cos_pi.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/cos_pi.hpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -17,8 +17,8 @@
namespace boost{ namespace math{
-template <class T>
-T cos_pi(T x)
+template <class T, class Policy>
+T cos_pi(T x, const Policy& pol)
{
BOOST_MATH_STD_USING // ADL of std names
// cos of pi*x:
@@ -31,7 +31,7 @@
}
T rem = floor(x);
- if(itrunc(rem) & 1)
+ if(itrunc(rem, pol) & 1)
invert = !invert;
rem = x - rem;
if(rem > 0.5f)
@@ -46,10 +46,10 @@
return invert ? -rem : rem;
}
-template <class T, class Policy>
-inline T cos_pi(T x, const Policy&)
+template <class T>
+inline T cos_pi(T x)
{
- return cos_pi(x);
+ return cos_pi(x, policies::policy<>());
}
} // namespace math
Modified: sandbox/math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -234,7 +234,7 @@
v = -v; // v is non-negative from here
kind |= need_k;
}
- n = iround(v);
+ n = iround(v, pol);
u = v - n; // -1/2 <= u < 1/2
BOOST_MATH_INSTRUMENT_VARIABLE(n);
BOOST_MATH_INSTRUMENT_VARIABLE(u);
Modified: sandbox/math_toolkit/boost/math/special_functions/detail/bessel_jy.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/bessel_jy.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/bessel_jy.hpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -225,7 +225,7 @@
v = -v; // v is non-negative from here
kind = need_j|need_y; // need both for reflection formula
}
- n = iround(v);
+ n = iround(v, pol);
u = v - n; // -1/2 <= u < 1/2
if (x == 0)
Modified: sandbox/math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -205,7 +205,7 @@
//
T tolerance = ldexp(1.0f, (2 * policies::digits<T, Policy>()) / 3);
- switch(itrunc(df))
+ switch(itrunc(df, Policy()))
{
case 1:
{
@@ -369,7 +369,7 @@
// where we use Shaw's tail series.
// The crossover point is roughly exponential in -df:
//
- T crossover = ldexp(1.0f, iround(df / -0.654f));
+ T crossover = ldexp(1.0f, iround(df / -0.654f, pol));
if(u > crossover)
{
result = boost::math::detail::inverse_students_t_hill(df, u, pol);
Modified: sandbox/math_toolkit/boost/math/special_functions/factorials.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/factorials.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/factorials.hpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -163,7 +163,7 @@
// handle it, split the product up into three parts:
//
T xp1 = x + 1;
- unsigned n2 = itrunc(floor(xp1));
+ unsigned n2 = itrunc(floor(xp1), pol);
if(n2 == xp1)
return 0;
T result = boost::math::tgamma_delta_ratio(xp1, -static_cast<T>(n2), pol);
Modified: sandbox/math_toolkit/boost/math/special_functions/gamma.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/gamma.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/gamma.hpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -165,7 +165,7 @@
}
if((floor(z) == z) && (z < max_factorial<T>::value))
{
- result *= unchecked_factorial<T>(itrunc(z) - 1);
+ result *= unchecked_factorial<T>(itrunc(z, pol) - 1);
}
else
{
@@ -370,7 +370,7 @@
BOOST_MATH_INSTRUMENT_CODE(prefix);
if((floor(z) == z) && (z < max_factorial<T>::value))
{
- prefix *= unchecked_factorial<T>(itrunc(z) - 1);
+ prefix *= unchecked_factorial<T>(itrunc(z, pol) - 1);
}
else
{
@@ -1095,7 +1095,7 @@
//
if((z <= max_factorial<T>::value) && (z + delta <= max_factorial<T>::value))
{
- return unchecked_factorial<T>((unsigned)itrunc(z) - 1) / unchecked_factorial<T>((unsigned)itrunc(z + delta) - 1);
+ return unchecked_factorial<T>((unsigned)itrunc(z, pol) - 1) / unchecked_factorial<T>((unsigned)itrunc(z + delta) - 1);
}
}
if(fabs(delta) < 20)
Modified: sandbox/math_toolkit/boost/math/special_functions/modf.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/modf.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/modf.hpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -15,33 +15,53 @@
namespace boost{ namespace math{
+template <class T, class Policy>
+inline T modf(const T& v, T* ipart, const Policy& pol)
+{
+ *ipart = trunc(v, pol);
+ return v - *ipart;
+}
template <class T>
inline T modf(const T& v, T* ipart)
{
- *ipart = trunc(v);
- return v - *ipart;
+ return modf(v, ipart, policies::policy<>());
}
+template <class T, class Policy>
+inline T modf(const T& v, int* ipart, const Policy& pol)
+{
+ *ipart = itrunc(v, pol);
+ return v - *ipart;
+}
template <class T>
inline T modf(const T& v, int* ipart)
{
- *ipart = itrunc(v);
- return v - *ipart;
+ return modf(v, ipart, policies::policy<>());
}
+template <class T, class Policy>
+inline T modf(const T& v, long* ipart, const Policy& pol)
+{
+ *ipart = ltrunc(v, pol);
+ return v - *ipart;
+}
template <class T>
inline T modf(const T& v, long* ipart)
{
- *ipart = ltrunc(v);
- return v - *ipart;
+ return modf(v, ipart, policies::policy<>());
}
#ifdef BOOST_HAS_LONG_LONG
+template <class T, class Policy>
+inline T modf(const T& v, long long* ipart, const Policy& pol)
+{
+ *ipart = lltrunc(v, pol);
+ return v - *ipart;
+}
template <class T>
inline T modf(const T& v, long long* ipart)
{
- *ipart = lltrunc(v);
- return v - *ipart;
+ return modf(v, ipart, policies::policy<>());
}
#endif
Modified: sandbox/math_toolkit/boost/math/special_functions/round.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/round.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/round.hpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -11,15 +11,24 @@
#endif
#include <boost/math/tools/config.hpp>
+#include <boost/math/policies/error_handling.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
namespace boost{ namespace math{
-template <class T>
-inline T round(const T& v)
+template <class T, class Policy>
+inline T round(const T& v, const Policy& pol)
{
BOOST_MATH_STD_USING
+ if(!(boost::math::isfinite)(v))
+ return policies::raise_rounding_error("boost::math::round<%1%>(%1%)", 0, v, pol);
return floor(v + 0.5f);
}
+template <class T>
+inline T round(const T& v)
+{
+ return round(v, policies::policy<>());
+}
//
// The following functions will not compile unless T has an
// implicit convertion to the integer types. For user-defined
@@ -29,24 +38,48 @@
// namespace as the UDT: these will then be found via argument
// dependent lookup. See our concept archetypes for examples.
//
+template <class T, class Policy>
+inline int iround(const T& v, const Policy& pol)
+{
+ T r = boost::math::round(v, pol);
+ if(fabs(r) > (std::numeric_limits<int>::max)())
+ return static_cast<int>(policies::raise_rounding_error("boost::math::iround<%1%>(%1%)", 0, v, pol));
+ return static_cast<int>(r);
+}
template <class T>
inline int iround(const T& v)
{
- return static_cast<int>(boost::math::round(v));
+ return iround(v, policies::policy<>());
}
+template <class T, class Policy>
+inline long lround(const T& v, const Policy& pol)
+{
+ T r = boost::math::round(v, pol);
+ if(fabs(r) > (std::numeric_limits<long>::max)())
+ return static_cast<long int>(policies::raise_rounding_error("boost::math::lround<%1%>(%1%)", 0, v, pol));
+ return static_cast<long int>(r);
+}
template <class T>
inline long lround(const T& v)
{
- return static_cast<long int>(boost::math::round(v));
+ return lround(v, policies::policy<>());
}
#ifdef BOOST_HAS_LONG_LONG
+template <class T, class Policy>
+inline long long llround(const T& v, const Policy& pol)
+{
+ T r = boost::math::round(v, pol);
+ if(fabs(r) > (std::numeric_limits<long long>::max)())
+ return static_cast<long long>(policies::raise_rounding_error("boost::math::llround<%1%>(%1%)", 0, v, pol));
+ return static_cast<long long>(r);
+}
template <class T>
inline long long llround(const T& v)
{
- return static_cast<long long>(boost::math::round(v));
+ return llround(v, policies::policy<>());
}
#endif
Modified: sandbox/math_toolkit/boost/math/special_functions/sin_pi.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/sin_pi.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/sin_pi.hpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -17,8 +17,8 @@
namespace boost{ namespace math{
-template <class T>
-T sin_pi(T x)
+template <class T, class Policy>
+T sin_pi(T x, const Policy& pol)
{
BOOST_MATH_STD_USING // ADL of std names
if(x < 0)
@@ -36,7 +36,7 @@
invert = false;
T rem = floor(x);
- if(itrunc(rem) & 1)
+ if(itrunc(rem, pol) & 1)
invert = !invert;
rem = x - rem;
if(rem > 0.5f)
@@ -48,10 +48,10 @@
return invert ? -rem : rem;
}
-template <class T, class Policy>
-inline T sin_pi(T x, const Policy&)
+template <class T>
+inline T sin_pi(T x)
{
- return boost::math::sin_pi(x);
+ return boost::math::sin_pi(x, policies::policy<>());
}
} // namespace math
Modified: sandbox/math_toolkit/boost/math/special_functions/trunc.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/trunc.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/trunc.hpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -11,15 +11,24 @@
#endif
#include <boost/math/tools/config.hpp>
+#include <boost/math/policies/error_handling.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
namespace boost{ namespace math{
-template <class T>
-inline T trunc(const T& v)
+template <class T, class Policy>
+inline T trunc(const T& v, const Policy& pol)
{
BOOST_MATH_STD_USING
+ if(!(boost::math::isfinite)(v))
+ return policies::raise_rounding_error("boost::math::trunc<%1%>(%1%)", 0, v, pol);
return (v >= 0) ? floor(v) : ceil(v);
}
+template <class T>
+inline T trunc(const T& v)
+{
+ return trunc(v, policies::policy<>());
+}
//
// The following functions will not compile unless T has an
// implicit convertion to the integer types. For user-defined
@@ -29,24 +38,48 @@
// namespace as the UDT: these will then be found via argument
// dependent lookup. See our concept archetypes for examples.
//
+template <class T, class Policy>
+inline int itrunc(const T& v, const Policy& pol)
+{
+ T r = boost::math::trunc(v, pol);
+ if(fabs(r) > (std::numeric_limits<int>::max)())
+ return static_cast<int>(policies::raise_rounding_error("boost::math::itrunc<%1%>(%1%)", 0, v, pol));
+ return static_cast<int>(r);
+}
template <class T>
inline int itrunc(const T& v)
{
- return static_cast<int>(boost::math::trunc(v));
+ return itrunc(v, policies::policy<>());
}
+template <class T, class Policy>
+inline long ltrunc(const T& v, const Policy& pol)
+{
+ T r = boost::math::trunc(v, pol);
+ if(fabs(r) > (std::numeric_limits<long>::max)())
+ return static_cast<long>(policies::raise_rounding_error("boost::math::ltrunc<%1%>(%1%)", 0, v, pol));
+ return static_cast<long>(r);
+}
template <class T>
inline long ltrunc(const T& v)
{
- return static_cast<long int>(boost::math::trunc(v));
+ return ltrunc(v, policies::policy<>());
}
#ifdef BOOST_HAS_LONG_LONG
+template <class T, class Policy>
+inline long long lltrunc(const T& v, const Policy& pol)
+{
+ T r = boost::math::trunc(v, pol);
+ if(fabs(r) > (std::numeric_limits<long long>::max)())
+ return static_cast<long long>(policies::raise_rounding_error("boost::math::lltrunc<%1%>(%1%)", 0, v, pol));
+ return static_cast<long long>(r);
+}
template <class T>
inline long long lltrunc(const T& v)
{
- return static_cast<long long>(boost::math::trunc(v));
+ return lltrunc(v, policies::policy<>());
}
#endif
Modified: sandbox/math_toolkit/libs/math/test/Jamfile.v2
==============================================================================
--- sandbox/math_toolkit/libs/math/test/Jamfile.v2 (original)
+++ sandbox/math_toolkit/libs/math/test/Jamfile.v2 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -349,6 +349,7 @@
compile compile_test/dist_gamma_incl_test.cpp ;
compile compile_test/dist_lognormal_incl_test.cpp ;
compile compile_test/dist_neg_binom_incl_test.cpp ;
+compile compile_test/dist_nc_chi_squ_incl_test.cpp ;
compile compile_test/dist_normal_incl_test.cpp ;
compile compile_test/dist_poisson_incl_test.cpp ;
compile compile_test/dist_students_t_incl_test.cpp ;
Added: sandbox/math_toolkit/libs/math/test/compile_test/dist_nc_chi_squ_incl_test.cpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/libs/math/test/compile_test/dist_nc_chi_squ_incl_test.cpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -0,0 +1,25 @@
+// Copyright John Maddock 2008.
+// 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)
+//
+// Basic sanity check that header <boost/math/distributions/non_central_chi_squared.hpp>
+// #includes all the files that it needs to.
+//
+#include <boost/math/distributions/non_central_chi_squared.hpp>
+//
+// Note this header includes no other headers, this is
+// important if this test is to be meaningful:
+//
+#include "test_compile_result.hpp"
+
+void check()
+{
+ TEST_DIST_FUNC(non_central_chi_squared)
+}
+
+template class boost::math::non_central_chi_squared_distribution<float, boost::math::policies::policy<> >;
+template class boost::math::non_central_chi_squared_distribution<double, boost::math::policies::policy<> >;
+#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+template class boost::math::non_central_chi_squared_distribution<long double, boost::math::policies::policy<> >;
+#endif
Modified: sandbox/math_toolkit/libs/math/test/test_nc_chi_squared.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/test/test_nc_chi_squared.cpp (original)
+++ sandbox/math_toolkit/libs/math/test/test_nc_chi_squared.cpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -438,30 +438,30 @@
value_type pt = data[i][2];
BOOST_CHECK_CLOSE_EX(pt, p, precision, i);
}
- //
- // Sanity check mode as well, note this may well overflow
- // since we don't know how to compute PDF's (and hence the mode)
- // for large values of the parameters, in addition accuracy of
- // the mode is at *best* the square root of the accuracy of the PDF:
- //
- try
- {
- value_type m = mode(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]));
- value_type p = pdf(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), m);
- BOOST_CHECK_EX(pdf(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), m * (1 + sqrt(precision))) <= p, i);
- BOOST_CHECK_EX(pdf(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), m * (1 - sqrt(precision))) <= p, i);
- }
- catch(const std::overflow_error&)
- {
- }
- //
- // Sanity check degrees-of-freedom finder, don't bother at float
- // precision though as there's not enough data in the probability
- // values to get back to the correct degrees of freedom or
- // non-cenrality parameter:
- //
if(boost::math::tools::digits<value_type>() > 50)
{
+ //
+ // Sanity check mode, note this may well overflow
+ // since we don't know how to compute PDF's (and hence the mode)
+ // for large values of the parameters, in addition accuracy of
+ // the mode is at *best* the square root of the accuracy of the PDF:
+ //
+ try
+ {
+ value_type m = mode(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]));
+ value_type p = pdf(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), m);
+ BOOST_CHECK_EX(pdf(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), m * (1 + sqrt(precision) * 10)) <= p, i);
+ BOOST_CHECK_EX(pdf(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), m * (1 - sqrt(precision)) * 10) <= p, i);
+ }
+ catch(const std::overflow_error&)
+ {
+ }
+ //
+ // Sanity check degrees-of-freedom finder, don't bother at float
+ // precision though as there's not enough data in the probability
+ // values to get back to the correct degrees of freedom or
+ // non-cenrality parameter:
+ //
try{
if((data[i][3] < 0.99) && (data[i][3] != 0))
{
Modified: sandbox/math_toolkit/libs/math/test/test_round.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/test/test_round.cpp (original)
+++ sandbox/math_toolkit/libs/math/test/test_round.cpp 2008-01-23 05:41:06 EST (Wed, 23 Jan 2008)
@@ -152,6 +152,67 @@
}
#endif
}
+ //
+ // Finish off by testing the error handlers:
+ //
+ BOOST_CHECK_THROW(boost::math::iround(static_cast<T>(1e20)), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::iround(static_cast<T>(-1e20)), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::lround(static_cast<T>(1e20)), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::lround(static_cast<T>(-1e20)), boost::math::rounding_error);
+#ifdef BOOST_HAS_LONG_LONG
+ BOOST_CHECK_THROW(boost::math::llround(static_cast<T>(1e20)), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::llround(static_cast<T>(-1e20)), boost::math::rounding_error);
+#endif
+ if(std::numeric_limits<T>::has_infinity)
+ {
+ BOOST_CHECK_THROW(boost::math::round(std::numeric_limits<T>::infinity()), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::iround(std::numeric_limits<T>::infinity()), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::iround(-std::numeric_limits<T>::infinity()), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::lround(std::numeric_limits<T>::infinity()), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::lround(-std::numeric_limits<T>::infinity()), boost::math::rounding_error);
+ #ifdef BOOST_HAS_LONG_LONG
+ BOOST_CHECK_THROW(boost::math::llround(std::numeric_limits<T>::infinity()), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::llround(-std::numeric_limits<T>::infinity()), boost::math::rounding_error);
+ #endif
+ }
+ if(std::numeric_limits<T>::has_quiet_NaN)
+ {
+ BOOST_CHECK_THROW(boost::math::round(std::numeric_limits<T>::quiet_NaN()), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::iround(std::numeric_limits<T>::quiet_NaN()), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::lround(std::numeric_limits<T>::quiet_NaN()), boost::math::rounding_error);
+ #ifdef BOOST_HAS_LONG_LONG
+ BOOST_CHECK_THROW(boost::math::llround(std::numeric_limits<T>::quiet_NaN()), boost::math::rounding_error);
+ #endif
+ }
+ BOOST_CHECK_THROW(boost::math::itrunc(static_cast<T>(1e20)), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::itrunc(static_cast<T>(-1e20)), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::ltrunc(static_cast<T>(1e20)), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::ltrunc(static_cast<T>(-1e20)), boost::math::rounding_error);
+#ifdef BOOST_HAS_LONG_LONG
+ BOOST_CHECK_THROW(boost::math::lltrunc(static_cast<T>(1e20)), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::lltrunc(static_cast<T>(-1e20)), boost::math::rounding_error);
+#endif
+ if(std::numeric_limits<T>::has_infinity)
+ {
+ BOOST_CHECK_THROW(boost::math::trunc(std::numeric_limits<T>::infinity()), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::itrunc(std::numeric_limits<T>::infinity()), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::itrunc(-std::numeric_limits<T>::infinity()), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::ltrunc(std::numeric_limits<T>::infinity()), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::ltrunc(-std::numeric_limits<T>::infinity()), boost::math::rounding_error);
+ #ifdef BOOST_HAS_LONG_LONG
+ BOOST_CHECK_THROW(boost::math::lltrunc(std::numeric_limits<T>::infinity()), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::lltrunc(-std::numeric_limits<T>::infinity()), boost::math::rounding_error);
+ #endif
+ }
+ if(std::numeric_limits<T>::has_quiet_NaN)
+ {
+ BOOST_CHECK_THROW(boost::math::trunc(std::numeric_limits<T>::quiet_NaN()), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::itrunc(std::numeric_limits<T>::quiet_NaN()), boost::math::rounding_error);
+ BOOST_CHECK_THROW(boost::math::ltrunc(std::numeric_limits<T>::quiet_NaN()), boost::math::rounding_error);
+ #ifdef BOOST_HAS_LONG_LONG
+ BOOST_CHECK_THROW(boost::math::lltrunc(std::numeric_limits<T>::quiet_NaN()), boost::math::rounding_error);
+ #endif
+ }
}
int test_main(int, char* [])
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