|
|
Boost-Commit : |
From: johnmaddock_at_[hidden]
Date: 2007-06-20 14:03:09
Author: johnmaddock
Date: 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
New Revision: 7112
URL: http://svn.boost.org/trac/boost/changeset/7112
Log:
Converted gamma.hpp to policy based implementation.
Added:
sandbox/math_toolkit/policy/boost/math/special_functions/detail/lgamma_small.hpp
Text files modified:
sandbox/math_toolkit/policy/boost/math/concepts/real_concept.hpp | 17
sandbox/math_toolkit/policy/boost/math/policy/error_handling.hpp | 21
sandbox/math_toolkit/policy/boost/math/policy/policy.hpp | 68 ++
sandbox/math_toolkit/policy/boost/math/special_functions/beta.hpp | 6
sandbox/math_toolkit/policy/boost/math/special_functions/detail/erf_inv.hpp | 12
sandbox/math_toolkit/policy/boost/math/special_functions/detail/gamma_inva.hpp | 38 +
sandbox/math_toolkit/policy/boost/math/special_functions/detail/igamma_inverse.hpp | 94 ++-
sandbox/math_toolkit/policy/boost/math/special_functions/detail/igamma_large.hpp | 41
sandbox/math_toolkit/policy/boost/math/special_functions/detail/t_distribution_inv.hpp | 4
sandbox/math_toolkit/policy/boost/math/special_functions/erf.hpp | 11
sandbox/math_toolkit/policy/boost/math/special_functions/expm1.hpp | 86 +-
sandbox/math_toolkit/policy/boost/math/special_functions/gamma.hpp | 989 +++++++++++----------------------------
sandbox/math_toolkit/policy/boost/math/special_functions/lanczos.hpp | 87 ++
sandbox/math_toolkit/policy/boost/math/special_functions/log1p.hpp | 153 ++++--
sandbox/math_toolkit/policy/boost/math/special_functions/math_fwd.hpp | 14
sandbox/math_toolkit/policy/boost/math/special_functions/powm1.hpp | 16
sandbox/math_toolkit/policy/boost/math/special_functions/sqrt1pm1.hpp | 12
sandbox/math_toolkit/policy/boost/math/tools/promotion.hpp | 10
sandbox/math_toolkit/policy/libs/math/test/test_igamma.cpp | 10
sandbox/math_toolkit/policy/libs/math/test/test_igamma_inv.cpp | 8
sandbox/math_toolkit/policy/libs/math/test/test_igamma_inva.cpp | 6
sandbox/math_toolkit/policy/libs/math/test/test_tgamma_ratio.cpp | 9
22 files changed, 802 insertions(+), 910 deletions(-)
Modified: sandbox/math_toolkit/policy/boost/math/concepts/real_concept.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/concepts/real_concept.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/concepts/real_concept.hpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -24,6 +24,7 @@
#include <boost/limits.hpp>
#include <boost/math/tools/real_cast.hpp>
#include <boost/math/tools/precision.hpp>
+#include <boost/math/policy/policy.hpp>
#include <ostream>
#include <istream>
@@ -336,6 +337,22 @@
} // namespace tools
+namespace policy{
+
+template <>
+inline int digits<concepts::real_concept, policy<> >(BOOST_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::real_concept))
+{ // Assume number of significand bits is same as long double,
+ // unless std::numeric_limits<T>::is_specialized to provide digits.
+ return tools::digits<long double>();
+ // Note that if numeric_limits real concept is NOT specialized to provide digits10
+ // (or max_digits10) then the default precision of 6 decimal digits will be used
+ // by Boost test (giving misleading error messages like
+ // "difference between {9.79796} and {9.79796} exceeds 5.42101e-19%"
+ // and by Boost lexical cast and serialization causing loss of accuracy.
+}
+
+}
+
} // namespace math
} // namespace boost
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-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -419,7 +419,7 @@
using namespace std;
if(fabs(val) > tools::max_value<R>())
{
- *result = raise_overflow_error<R>(function, 0, pol);
+ *result = static_cast<R>(raise_overflow_error<R>(function, 0, pol));
return true;
}
return false;
@@ -429,7 +429,7 @@
{
if((val != 0) && (static_cast<R>(val) == 0))
{
- *result = raise_underflow_error<R>(function, 0, pol);
+ *result = static_cast<R>(raise_underflow_error<R>(function, 0, pol));
return true;
}
return false;
@@ -440,18 +440,18 @@
using namespace std;
if((fabs(val) < static_cast<T>(tools::min_value<R>())) && (static_cast<R>(val) != 0))
{
- *result = raise_denorm_error<R>(function, 0, val, pol);
+ *result = static_cast<R>(raise_denorm_error<R>(function, 0, static_cast<R>(val), pol));
return true;
}
return false;
}
template <class R, class T>
-inline bool check_overflow(T val, R* result, const overflow_error<ignore_error>&){ return false; }
+inline bool check_overflow(T val, R* result, const char* function, const overflow_error<ignore_error>&){ return false; }
template <class R, class T>
-inline bool check_underflow(T val, R* result, const overflow_error<ignore_error>&){ return false; }
+inline bool check_underflow(T val, R* result, const char* function, const underflow_error<ignore_error>&){ return false; }
template <class R, class T>
-inline bool check_denormflow(T val, R* result, const overflow_error<ignore_error>&){ return false; }
+inline bool check_denorm(T val, R* result, const char* function, const denorm_error<ignore_error>&){ return false; }
}
@@ -475,6 +475,15 @@
return static_cast<R>(val);
}
+template <class Policy>
+inline void check_series_iterations(const char* function, boost::uintmax_t max_iter, const Policy& pol)
+{
+ if(max_iter >= BOOST_MATH_MAX_ITER)
+ raise_evaluation_error<boost::uintmax_t>(
+ function,
+ "Series evaluation exceeded %1% iterations, giving up now.", max_iter, pol);
+}
+
} //namespace policy
}} // namespaces boost/math
Modified: sandbox/math_toolkit/policy/boost/math/policy/policy.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/policy/policy.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/policy/policy.hpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -15,6 +15,7 @@
#include <boost/mpl/size.hpp>
#include <boost/type_traits/is_same.hpp>
#include <boost/static_assert.hpp>
+#include <boost/math/tools/config.hpp>
namespace boost{ namespace math{ namespace policy{
@@ -46,8 +47,12 @@
#define BOOST_MATH_PROMOTE_FLOAT_POLICY true
#endif
#ifndef BOOST_MATH_PROMOTE_DOUBLE_POLICY
+#ifdef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false
+#else
#define BOOST_MATH_PROMOTE_DOUBLE_POLICY true
#endif
+#endif
#ifndef BOOST_MATH_DISCRETE_QUANTILE_POLICY
#define BOOST_MATH_DISCRETE_QUANTILE_POLICY integer_outside
#endif
@@ -603,6 +608,69 @@
typedef typename mpl::if_<typename Policy::double_promote_type, long double, double>::type type;
};
+template <class Real, class Policy>
+struct precision
+{
+ typedef typename Policy::precision_type precision_type;
+ typedef typename mpl::if_c<
+ ((std::numeric_limits<Real>::is_specialized == 0) || (std::numeric_limits<Real>::digits == 0)),
+ // Possibly unknown precision:
+ precision_type,
+ typename mpl::if_c<
+#ifndef __BORLANDC__
+ ((::std::numeric_limits<Real>::digits <= precision_type::value)
+ || (Policy::precision_type::value <= 0)),
+#else
+ ((::std::numeric_limits<Real>::digits <= ::boost::math::policy::precision::precision_type::value)
+ || (::boost::math::policy::precision::precision_type::value <= 0)),
+#endif
+ // Default case, full precision for RealType:
+ digits2< ::std::numeric_limits<Real>::digits>,
+ // User customised precision:
+ precision_type
+ >::type
+ >::type type;
+};
+
+template <class T, class Policy>
+inline int digits(BOOST_EXPLICIT_TEMPLATE_TYPE(T))
+{
+#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
+ BOOST_STATIC_ASSERT( ::std::numeric_limits<T>::is_specialized);
+#else
+ BOOST_ASSERT(::std::numeric_limits<T>::is_specialized);
+#endif
+ typedef typename precision<T, Policy>::type p_t;
+ return p_t::value;
+}
+
+namespace detail{
+
+template <class A1,
+ class A2,
+ class A3,
+ class A4,
+ class A5,
+ class A6,
+ class A7,
+ class A8,
+ class A9,
+ class A10,
+ class A11>
+char test_is_policy(const policy<A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11>*);
+double test_is_policy(...);
+
+template <class P>
+struct is_policy_imp
+{
+ BOOST_STATIC_CONSTANT(bool, value = (sizeof(test_is_policy(static_cast<P*>(0))) == 1));
+};
+
+}
+
+template <class P>
+struct is_policy : public mpl::bool_< ::boost::math::policy::detail::is_policy_imp<P>::value> {};
+
}}} // namespaces
#endif // BOOST_MATH_POLICY_HPP
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/beta.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/beta.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/beta.hpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -711,17 +711,17 @@
u = -t * lx;
// and from from 9.2:
T prefix;
- T h = regularised_gamma_prefix(b, u, l);
+ T h = regularised_gamma_prefix(b, u, policy::policy<>(), l);
if(h <= tools::min_value<T>())
return s0;
if(normalised)
{
- prefix = h / tgamma_delta_ratio_imp(a, b, l);
+ prefix = h / tgamma_delta_ratio_imp(a, b, policy::policy<>());
prefix /= pow(t, b);
}
else
{
- prefix = full_igamma_prefix(b, u) / pow(t, b);
+ prefix = full_igamma_prefix(b, u, policy::policy<>()) / pow(t, b);
}
prefix *= mult;
//
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/detail/erf_inv.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/detail/erf_inv.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/detail/erf_inv.hpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -428,6 +428,18 @@
detail::erf_inv_imp(static_cast<eval_type>(p), static_cast<eval_type>(q), static_cast<tag_type const*>(0)), BOOST_CURRENT_FUNCTION);
}
+template <class T, class Policy>
+inline typename tools::promote_args<T>::type erfc_inv(T z, const Policy&)
+{
+ return erfc_inv(z);
+}
+
+template <class T, class Policy>
+inline typename tools::promote_args<T>::type erf_inv(T z, const Policy&)
+{
+ return erf_inv(z);
+}
+
} // namespace math
} // namespace boost
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/detail/gamma_inva.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/detail/gamma_inva.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/detail/gamma_inva.hpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -18,21 +18,21 @@
namespace boost{ namespace math{ namespace detail{
-template <class T>
+template <class T, class Policy>
struct gamma_inva_t
{
gamma_inva_t(T z_, T p_, bool invert_) : z(z_), p(p_), invert(invert_) {}
T operator()(T a)
{
- return invert ? p - boost::math::gamma_q(a, z) : boost::math::gamma_p(a, z) - p;
+ return invert ? p - boost::math::gamma_q(a, z, Policy()) : boost::math::gamma_p(a, z, Policy()) - p;
}
private:
T z, p;
bool invert;
};
-template <class T>
-T inverse_poisson_cornish_fisher(T lambda, T p, T q)
+template <class T, class Policy>
+T inverse_poisson_cornish_fisher(T lambda, T p, T q, const Policy& pol)
{
using namespace std;
// mean:
@@ -44,7 +44,7 @@
// kurtosis:
// T k = 1/lambda;
// Get the inverse of a std normal distribution:
- T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p) * constants::root_two<T>();
+ T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>();
// Set the sign:
if(p < 0.5)
x = -x;
@@ -62,8 +62,8 @@
return w > tools::min_value<T>() ? w : tools::min_value<T>();
}
-template <class T>
-T gamma_inva_imp(const T& z, const T& p, const T& q)
+template <class T, class Policy>
+T gamma_inva_imp(const T& z, const T& p, const T& q, const Policy& pol)
{
using namespace std; // for ADL of std lib math functions
//
@@ -81,11 +81,11 @@
// Function object, this is the functor whose root
// we have to solve:
//
- gamma_inva_t<T> f(z, (p < q) ? p : q, (p < q) ? false : true);
+ gamma_inva_t<T, Policy> f(z, (p < q) ? p : q, (p < q) ? false : true);
//
// Tolerance: full precision.
//
- tools::eps_tolerance<T> tol(tools::digits<T>());
+ tools::eps_tolerance<T> tol(policy::digits<T, Policy>());
//
// Now figure out a starting guess for what a may be,
// we'll start out with a value that'll put p or q
@@ -105,7 +105,7 @@
// when z is very small. Also set our step-factor according
// to how accurate we think the result is likely to be:
//
- guess = 1 + inverse_poisson_cornish_fisher(z, q, p);
+ guess = 1 + inverse_poisson_cornish_fisher(z, q, p, pol);
if(z > 5)
{
if(z > 1000)
@@ -140,22 +140,34 @@
//
std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, false, tol, max_iter);
if(max_iter >= 200)
- tools::logic_error<T>(BOOST_CURRENT_FUNCTION, "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first);
+ policy::raise_evaluation_error<T>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first, pol);
return (r.first + r.second) / 2;
}
} // namespace detail
+template <class T, class Policy>
+inline T gamma_p_inva(T x, T p, const Policy& pol)
+{
+ return detail::gamma_inva_imp(x, p, 1 - p, pol);
+}
+
+template <class T, class Policy>
+inline T gamma_q_inva(T x, T q, const Policy& pol)
+{
+ return detail::gamma_inva_imp(x, 1 - q, q, pol);
+}
+
template <class T>
inline T gamma_p_inva(T x, T p)
{
- return detail::gamma_inva_imp(x, p, 1 - p);
+ return detail::gamma_inva_imp(x, p, 1 - p, policy::policy<>());
}
template <class T>
inline T gamma_q_inva(T x, T q)
{
- return detail::gamma_inva_imp(x, 1 - q, q);
+ return detail::gamma_inva_imp(x, 1 - q, q, policy::policy<>());
}
} // namespace math
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/detail/igamma_inverse.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/detail/igamma_inverse.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/detail/igamma_inverse.hpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -71,8 +71,8 @@
return sum;
}
-template <class T>
-inline T didonato_FN(T p, T a, T x, unsigned N, T tolerance)
+template <class T, class Policy>
+inline T didonato_FN(T p, T a, T x, unsigned N, T tolerance, const Policy& pol)
{
//
// Computation of the Incomplete Gamma Function Ratios and their Inverse
@@ -83,12 +83,12 @@
// See equation 34.
//
using namespace std;
- T u = log(p) + boost::math::lgamma(a + 1);
+ T u = log(p) + boost::math::lgamma(a + 1, pol);
return exp((u + x - log(didonato_SN(a, x, N, tolerance))) / a);
}
-template <class T>
-T estimate_inverse_gamma(T a, T p, T q)
+template <class T, class Policy>
+T estimate_inverse_gamma(T a, T p, T q, const Policy& pol)
{
//
// In order to understand what's going on here, you will
@@ -107,7 +107,7 @@
result = -log(q);
else if(a < 1)
{
- T g = boost::math::tgamma(a);
+ T g = boost::math::tgamma(a, pol);
T b = q * g;
if((b > 0.6) || ((b >= 0.45) && (a >= 0.3)))
{
@@ -205,7 +205,7 @@
else
{
T D = (std::max)(T(2), a * (a - 1));
- T lg = boost::math::lgamma(a);
+ T lg = boost::math::lgamma(a, pol);
T lb = log(q) + lg;
if(lb < -D * 2.3)
{
@@ -243,10 +243,10 @@
else
{
// DiDonato and Morris Eq 35:
- T z = didonato_FN(p, a, w, 0, T(0));
- z = didonato_FN(p, a, z, 2, T(0));
- z = didonato_FN(p, a, z, 2, T(0));
- z = didonato_FN(p, a, z, 3, T(0));
+ T z = didonato_FN(p, a, w, 0, T(0), pol);
+ z = didonato_FN(p, a, z, 2, T(0), pol);
+ z = didonato_FN(p, a, z, 2, T(0), pol);
+ z = didonato_FN(p, a, z, 3, T(0), pol);
if((z <= 0.01 * (a + 1)) || (z > 0.7 * (a + 1)))
{
@@ -255,8 +255,8 @@
else
{
// DiDonato and Morris Eq 36:
- T zb = didonato_FN(p, a, z, 100, T(1e-4));
- T u = log(p) + boost::math::lgamma(a + 1);
+ T zb = didonato_FN(p, a, z, 100, T(1e-4), pol);
+ T u = log(p) + boost::math::lgamma(a + 1, pol);
result = zb * (1 - (a * log(zb) - zb - u + log(didonato_SN(a, z, 100, T(1e-4)))) / (a - zb));
}
}
@@ -264,7 +264,7 @@
return result;
}
-template <class T>
+template <class T, class Policy>
struct gamma_p_inverse_func
{
gamma_p_inverse_func(T a_, T p_, bool inv) : a(a_), p(p_), invert(inv)
@@ -302,8 +302,8 @@
f = static_cast<T>(boost::math::detail::gamma_incomplete_imp(
static_cast<value_type>(a),
static_cast<value_type>(x),
- true, invert, evaluation_type(),
- &ft));
+ true, invert,
+ Policy(), &ft));
f1 = static_cast<T>(ft);
T f2;
T div = (a - x - 1) / x;
@@ -331,20 +331,22 @@
bool invert;
};
-template <class T>
-T gamma_p_inv_imp(T a, T p)
+template <class T, class Policy>
+T gamma_p_inv_imp(T a, T p, const Policy& pol)
{
using namespace std; // ADL of std functions.
+ static const char* function = "boost::math::gamma_p_inv<%1%>(%1%, %1%)";
+
if(a <= 0)
- tools::domain_error<T>(BOOST_CURRENT_FUNCTION, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a);
+ policy::raise_domain_error<T>(function, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a, pol);
if((p < 0) || (p > 1))
- tools::domain_error<T>(BOOST_CURRENT_FUNCTION, "Probabilty must be in the range [0,1] in the incomplete gamma function inverse (got p=%1%).", p);
+ policy::raise_domain_error<T>(function, "Probabilty must be in the range [0,1] in the incomplete gamma function inverse (got p=%1%).", p, pol);
if(p == 1)
return tools::max_value<T>();
if(p == 0)
return 0;
- T guess = detail::estimate_inverse_gamma(a, p, 1 - p);
+ T guess = detail::estimate_inverse_gamma(a, p, 1 - p, pol);
T lower = tools::min_value<T>();
if(guess <= lower)
guess = tools::min_value<T>();
@@ -355,36 +357,38 @@
// precision to prevent premature termination of the root-finding routine.
//
unsigned digits = (tools::digits<T>() * 2) / 3;
- if((a < 0.125) && (fabs(gamma_p_derivative(a, guess)) > 1 / sqrt(tools::epsilon<T>())))
+ if((a < 0.125) && (fabs(gamma_p_derivative(a, guess, pol)) > 1 / sqrt(tools::epsilon<T>())))
digits = tools::digits<T>() - 2;
//
// Go ahead and iterate:
//
guess = tools::halley_iterate(
- detail::gamma_p_inverse_func<T>(a, p, false),
+ detail::gamma_p_inverse_func<T, Policy>(a, p, false),
guess,
lower,
tools::max_value<T>(),
digits);
if(guess == lower)
- guess = tools::underflow_error<T>(BOOST_CURRENT_FUNCTION, "Expected result known to be non-zero, but is smaller than the smallest available number.");
+ guess = policy::raise_underflow_error<T>(function, "Expected result known to be non-zero, but is smaller than the smallest available number.", pol);
return guess;
}
-template <class T>
-T gamma_q_inv_imp(T a, T q)
+template <class T, class Policy>
+T gamma_q_inv_imp(T a, T q, const Policy& pol)
{
using namespace std; // ADL of std functions.
+ static const char* function = "boost::math::gamma_q_inv<%1%>(%1%, %1%)";
+
if(a <= 0)
- tools::domain_error<T>(BOOST_CURRENT_FUNCTION, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a);
+ policy::raise_domain_error<T>(function, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a, pol);
if((q < 0) || (q > 1))
- tools::domain_error<T>(BOOST_CURRENT_FUNCTION, "Probabilty must be in the range [0,1] in the incomplete gamma function inverse (got q=%1%).", q);
+ policy::raise_domain_error<T>(function, "Probabilty must be in the range [0,1] in the incomplete gamma function inverse (got q=%1%).", q, pol);
if(q == 0)
return tools::max_value<T>();
if(q == 1)
return 0;
- T guess = detail::estimate_inverse_gamma(a, 1 - q, q);
+ T guess = detail::estimate_inverse_gamma(a, 1 - q, q, pol);
T lower = tools::min_value<T>();
if(guess <= lower)
guess = tools::min_value<T>();
@@ -395,42 +399,56 @@
// precision to prevent premature termination of the root-finding routine.
//
unsigned digits = (tools::digits<T>() * 2) / 3;
- if((a < 0.125) && (fabs(gamma_p_derivative(a, guess)) > 1 / sqrt(tools::epsilon<T>())))
+ if((a < 0.125) && (fabs(gamma_p_derivative(a, guess, pol)) > 1 / sqrt(tools::epsilon<T>())))
digits = tools::digits<T>();
//
// Go ahead and iterate:
//
guess = tools::halley_iterate(
- detail::gamma_p_inverse_func<T>(a, q, true),
+ detail::gamma_p_inverse_func<T, Policy>(a, q, true),
guess,
lower,
tools::max_value<T>(),
digits);
if(guess == lower)
- guess = tools::underflow_error<T>(BOOST_CURRENT_FUNCTION, "Expected result known to be non-zero, but is smaller than the smallest available number.");
+ guess = policy::raise_underflow_error<T>(function, "Expected result known to be non-zero, but is smaller than the smallest available number.", pol);
return guess;
}
} // namespace detail
-template <class T1, class T2>
+template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type
- gamma_p_inv(T1 a, T2 p)
+ gamma_p_inv(T1 a, T2 p, const Policy& pol)
{
typedef typename tools::promote_args<T1, T2>::type result_type;
return detail::gamma_p_inv_imp(
static_cast<result_type>(a),
- static_cast<result_type>(p));
+ static_cast<result_type>(p), pol);
}
-template <class T1, class T2>
+template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type
- gamma_q_inv(T1 a, T2 p)
+ gamma_q_inv(T1 a, T2 p, const Policy& pol)
{
typedef typename tools::promote_args<T1, T2>::type result_type;
return detail::gamma_q_inv_imp(
static_cast<result_type>(a),
- static_cast<result_type>(p));
+ static_cast<result_type>(p), pol);
+}
+
+template <class T1, class T2>
+inline typename tools::promote_args<T1, T2>::type
+ gamma_p_inv(T1 a, T2 p)
+{
+ return gamma_p_inv(a, p, policy::policy<>());
+}
+
+template <class T1, class T2>
+inline typename tools::promote_args<T1, T2>::type
+ gamma_q_inv(T1 a, T2 p)
+{
+ return gamma_q_inv(a, p, policy::policy<>());
}
} // namespace math
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/detail/igamma_large.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/detail/igamma_large.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/detail/igamma_large.hpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -48,16 +48,11 @@
namespace boost{ namespace math{ namespace detail{
//
-// Forward declaration, computes log(1+x) - x:
-//
-template <class T>
-T log1pmx(T x);
-//
// This version will never be called (at runtime), it's a stub used
// when T is unsuitable to be passed to these routines:
//
-template <class T>
-inline T igamma_temme_large(T, T, mpl::int_<0> const *)
+template <class T, class Policy>
+inline T igamma_temme_large(T, T, const Policy& pol, mpl::int_<0> const *)
{
// stub function, should never actually be called
BOOST_ASSERT(0);
@@ -67,12 +62,12 @@
// This version is accurate for up to 64-bit mantissa's,
// (80-bit long double, or 10^-20).
//
-template <class T>
-T igamma_temme_large(T a, T x, mpl::int_<64> const *)
+template <class T, class Policy>
+T igamma_temme_large(T a, T x, const Policy& pol, mpl::int_<64> const *)
{
using namespace std; // ADL of std functions
T sigma = (x - a) / a;
- T phi = -log1pmx(sigma);
+ T phi = -log1pmx(sigma, pol);
T y = a * phi;
T z = sqrt(2 * phi);
if(x < a)
@@ -264,7 +259,7 @@
if(x < a)
result = -result;
- result += boost::math::erfc(sqrt(y)) / 2;
+ result += boost::math::erfc(sqrt(y), pol) / 2;
return result;
}
@@ -272,12 +267,12 @@
// This one is accurate for 53-bit mantissa's
// (IEEE double precision or 10^-17).
//
-template <class T>
-T igamma_temme_large(T a, T x, mpl::int_<53> const *)
+template <class T, class Policy>
+T igamma_temme_large(T a, T x, const Policy& pol, mpl::int_<53> const *)
{
using namespace std; // ADL of std functions
T sigma = (x - a) / a;
- T phi = -log1pmx(sigma);
+ T phi = -log1pmx(sigma, pol);
T y = a * phi;
T z = sqrt(2 * phi);
if(x < a)
@@ -406,7 +401,7 @@
if(x < a)
result = -result;
- result += boost::math::erfc(sqrt(y)) / 2;
+ result += boost::math::erfc(sqrt(y), pol) / 2;
return result;
}
@@ -414,12 +409,12 @@
// This one is accurate for 24-bit mantissa's
// (IEEE float precision, or 10^-8)
//
-template <class T>
-T igamma_temme_large(T a, T x, mpl::int_<24> const *)
+template <class T, class Policy>
+T igamma_temme_large(T a, T x, const Policy& pol, mpl::int_<24> const *)
{
using namespace std; // ADL of std functions
T sigma = (x - a) / a;
- T phi = -log1pmx(sigma);
+ T phi = -log1pmx(sigma, pol);
T y = a * phi;
T z = sqrt(2 * phi);
if(x < a)
@@ -459,7 +454,7 @@
if(x < a)
result = -result;
- result += boost::math::erfc(sqrt(y)) / 2;
+ result += boost::math::erfc(sqrt(y), pol) / 2;
return result;
}
@@ -470,12 +465,12 @@
// It's use for a < 200 is not recomended, that would
// require many more terms in the polynomials.
//
-template <class T>
-T igamma_temme_large(T a, T x, mpl::int_<113> const *)
+template <class T, class Policy>
+T igamma_temme_large(T a, T x, const Policy& pol, mpl::int_<113> const *)
{
using namespace std; // ADL of std functions
T sigma = (x - a) / a;
- T phi = -log1pmx(sigma);
+ T phi = -log1pmx(sigma, pol);
T y = a * phi;
T z = sqrt(2 * phi);
if(x < a)
@@ -756,7 +751,7 @@
if(x < a)
result = -result;
- result += boost::math::erfc(sqrt(y)) / 2;
+ result += boost::math::erfc(sqrt(y), pol) / 2;
return result;
}
Added: sandbox/math_toolkit/policy/boost/math/special_functions/detail/lgamma_small.hpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/detail/lgamma_small.hpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -0,0 +1,496 @@
+// (C) Copyright John Maddock 2006.
+// 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_SPECIAL_FUNCTIONS_DETAIL_LGAMMA_SMALL
+#define BOOST_MATH_SPECIAL_FUNCTIONS_DETAIL_LGAMMA_SMALL
+
+namespace boost{ namespace math{ namespace detail{
+
+//
+// lgamma for small arguments:
+//
+template <class T, class Policy>
+T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<64>&, const Policy& /* l */)
+{
+ // This version uses rational approximations for small
+ // values of z accurate enough for 64-bit mantissas
+ // (80-bit long doubles), works well for 53-bit doubles as well.
+ // L is only used to select the Lanczos function.
+
+ using namespace std; // for ADL of std names
+ T result = 0;
+ if(z < tools::epsilon<T>())
+ {
+ result = -log(z);
+ }
+ else if((zm1 == 0) || (zm2 == 0))
+ {
+ // nothing to do, result is zero....
+ }
+ else if(z > 2)
+ {
+ //
+ // Begin by performing argument reduction until
+ // z is in [2,3):
+ //
+ if(z >= 3)
+ {
+ do
+ {
+ z -= 1;
+ zm2 -= 1;
+ result += log(z);
+ }while(z >= 3);
+ // Update zm2, we need it below:
+ zm2 = z - 2;
+ }
+
+ //
+ // Use the following form:
+ //
+ // lgamma(z) = (z-2)(z+1)(Y + R(z-2))
+ //
+ // where R(z-2) is a rational approximation optimised for
+ // low absolute error - as long as it's absolute error
+ // is small compared to the constant Y - then any rounding
+ // error in it's computation will get wiped out.
+ //
+ // R(z-2) has the following properties:
+ //
+ // At double: Max error found: 4.231e-18
+ // At long double: Max error found: 1.987e-21
+ // Maximum Deviation Found (approximation error): 5.900e-24
+ //
+ static const T P[] = {
+ static_cast<T>(-0.180355685678449379109e-1L),
+ static_cast<T>(0.25126649619989678683e-1L),
+ static_cast<T>(0.494103151567532234274e-1L),
+ static_cast<T>(0.172491608709613993966e-1L),
+ static_cast<T>(-0.259453563205438108893e-3L),
+ static_cast<T>(-0.541009869215204396339e-3L),
+ static_cast<T>(-0.324588649825948492091e-4L)
+ };
+ static const T Q[] = {
+ static_cast<T>(0.1e1),
+ static_cast<T>(0.196202987197795200688e1L),
+ static_cast<T>(0.148019669424231326694e1L),
+ static_cast<T>(0.541391432071720958364e0L),
+ static_cast<T>(0.988504251128010129477e-1L),
+ static_cast<T>(0.82130967464889339326e-2L),
+ static_cast<T>(0.224936291922115757597e-3L),
+ static_cast<T>(-0.223352763208617092964e-6L)
+ };
+
+ static const float Y = 0.158963680267333984375e0f;
+
+ T r = zm2 * (z + 1);
+ T R = tools::evaluate_polynomial(P, zm2);
+ R /= tools::evaluate_polynomial(Q, zm2);
+
+ result += r * Y + r * R;
+ }
+ else
+ {
+ //
+ // If z is less than 1 use recurrance to shift to
+ // z in the interval [1,2]:
+ //
+ if(z < 1)
+ {
+ result += -log(z);
+ zm2 = zm1;
+ zm1 = z;
+ z += 1;
+ }
+ //
+ // Two approximations, on for z in [1,1.5] and
+ // one for z in [1.5,2]:
+ //
+ if(z <= 1.5)
+ {
+ //
+ // Use the following form:
+ //
+ // lgamma(z) = (z-1)(z-2)(Y + R(z-1))
+ //
+ // where R(z-1) is a rational approximation optimised for
+ // low absolute error - as long as it's absolute error
+ // is small compared to the constant Y - then any rounding
+ // error in it's computation will get wiped out.
+ //
+ // R(z-1) has the following properties:
+ //
+ // At double precision: Max error found: 1.230011e-17
+ // At 80-bit long double precision: Max error found: 5.631355e-21
+ // Maximum Deviation Found: 3.139e-021
+ // Expected Error Term: 3.139e-021
+
+ //
+ static const float Y = 0.52815341949462890625f;
+
+ static const T P[] = {
+ static_cast<T>(0.490622454069039543534e-1L),
+ static_cast<T>(-0.969117530159521214579e-1L),
+ static_cast<T>(-0.414983358359495381969e0L),
+ static_cast<T>(-0.406567124211938417342e0L),
+ static_cast<T>(-0.158413586390692192217e0L),
+ static_cast<T>(-0.240149820648571559892e-1L),
+ static_cast<T>(-0.100346687696279557415e-2L)
+ };
+ static const T Q[] = {
+ static_cast<T>(0.1e1L),
+ static_cast<T>(0.302349829846463038743e1L),
+ static_cast<T>(0.348739585360723852576e1L),
+ static_cast<T>(0.191415588274426679201e1L),
+ static_cast<T>(0.507137738614363510846e0L),
+ static_cast<T>(0.577039722690451849648e-1L),
+ static_cast<T>(0.195768102601107189171e-2L)
+ };
+
+ T r = tools::evaluate_polynomial(P, zm1) / tools::evaluate_polynomial(Q, zm1);
+ T prefix = zm1 * zm2;
+
+ result += prefix * Y + prefix * r;
+ }
+ else
+ {
+ //
+ // Use the following form:
+ //
+ // lgamma(z) = (2-z)(1-z)(Y + R(2-z))
+ //
+ // where R(2-z) is a rational approximation optimised for
+ // low absolute error - as long as it's absolute error
+ // is small compared to the constant Y - then any rounding
+ // error in it's computation will get wiped out.
+ //
+ // R(2-z) has the following properties:
+ //
+ // At double precision, max error found: 1.797565e-17
+ // At 80-bit long double precision, max error found: 9.306419e-21
+ // Maximum Deviation Found: 2.151e-021
+ // Expected Error Term: 2.150e-021
+ //
+ static const float Y = 0.452017307281494140625f;
+
+ static const T P[] = {
+ static_cast<T>(-0.292329721830270012337e-1L),
+ static_cast<T>(0.144216267757192309184e0L),
+ static_cast<T>(-0.142440390738631274135e0L),
+ static_cast<T>(0.542809694055053558157e-1L),
+ static_cast<T>(-0.850535976868336437746e-2L),
+ static_cast<T>(0.431171342679297331241e-3L)
+ };
+ static const T Q[] = {
+ static_cast<T>(0.1e1),
+ static_cast<T>(-0.150169356054485044494e1L),
+ static_cast<T>(0.846973248876495016101e0L),
+ static_cast<T>(-0.220095151814995745555e0L),
+ static_cast<T>(0.25582797155975869989e-1L),
+ static_cast<T>(-0.100666795539143372762e-2L),
+ static_cast<T>(-0.827193521891290553639e-6L)
+ };
+ T r = zm2 * zm1;
+ T R = tools::evaluate_polynomial(P, -zm2) / tools::evaluate_polynomial(Q, -zm2);
+
+ result += r * Y + r * R;
+ }
+ }
+ return result;
+}
+template <class T, class Policy>
+T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<113>&, const Policy& /* l */)
+{
+ //
+ // This version uses rational approximations for small
+ // values of z accurate enough for 113-bit mantissas
+ // (128-bit long doubles).
+ //
+ using namespace std; // for ADL of std names
+ T result = 0;
+ if(z < tools::epsilon<T>())
+ {
+ result = -log(z);
+ }
+ else if((zm1 == 0) || (zm2 == 0))
+ {
+ // nothing to do, result is zero....
+ }
+ else if(z > 2)
+ {
+ //
+ // Begin by performing argument reduction until
+ // z is in [2,3):
+ //
+ if(z >= 3)
+ {
+ do
+ {
+ z -= 1;
+ result += log(z);
+ }while(z >= 3);
+ zm2 = z - 2;
+ }
+
+ //
+ // Use the following form:
+ //
+ // lgamma(z) = (z-2)(z+1)(Y + R(z-2))
+ //
+ // where R(z-2) is a rational approximation optimised for
+ // low absolute error - as long as it's absolute error
+ // is small compared to the constant Y - then any rounding
+ // error in it's computation will get wiped out.
+ //
+ // Maximum Deviation Found (approximation error) 3.73e-37
+
+ static const T P[] = {
+ -0.018035568567844937910504030027467476655L,
+ 0.013841458273109517271750705401202404195L,
+ 0.062031842739486600078866923383017722399L,
+ 0.052518418329052161202007865149435256093L,
+ 0.01881718142472784129191838493267755758L,
+ 0.0025104830367021839316463675028524702846L,
+ -0.00021043176101831873281848891452678568311L,
+ -0.00010249622350908722793327719494037981166L,
+ -0.11381479670982006841716879074288176994e-4L,
+ -0.49999811718089980992888533630523892389e-6L,
+ -0.70529798686542184668416911331718963364e-8L
+ };
+ static const T Q[] = {
+ 1L,
+ 2.5877485070422317542808137697939233685L,
+ 2.8797959228352591788629602533153837126L,
+ 1.8030885955284082026405495275461180977L,
+ 0.69774331297747390169238306148355428436L,
+ 0.17261566063277623942044077039756583802L,
+ 0.02729301254544230229429621192443000121L,
+ 0.0026776425891195270663133581960016620433L,
+ 0.00015244249160486584591370355730402168106L,
+ 0.43997034032479866020546814475414346627e-5L,
+ 0.46295080708455613044541885534408170934e-7L,
+ -0.93326638207459533682980757982834180952e-11L,
+ 0.42316456553164995177177407325292867513e-13L
+ };
+
+ T R = tools::evaluate_polynomial(P, zm2);
+ R /= tools::evaluate_polynomial(Q, zm2);
+
+ static const float Y = 0.158963680267333984375F;
+
+ T r = zm2 * (z + 1);
+
+ result += r * Y + r * R;
+ }
+ else
+ {
+ //
+ // If z is less than 1 use recurrance to shift to
+ // z in the interval [1,2]:
+ //
+ if(z < 1)
+ {
+ result += -log(z);
+ zm2 = zm1;
+ zm1 = z;
+ z += 1;
+ }
+ //
+ // Three approximations, on for z in [1,1.35], [1.35,1.625] and [1.625,1]
+ //
+ if(z <= 1.35)
+ {
+ //
+ // Use the following form:
+ //
+ // lgamma(z) = (z-1)(z-2)(Y + R(z-1))
+ //
+ // where R(z-1) is a rational approximation optimised for
+ // low absolute error - as long as it's absolute error
+ // is small compared to the constant Y - then any rounding
+ // error in it's computation will get wiped out.
+ //
+ // R(z-1) has the following properties:
+ //
+ // Maximum Deviation Found (approximation error) 1.659e-36
+ // Expected Error Term (theoretical error) 1.343e-36
+ // Max error found at 128-bit long double precision 1.007e-35
+ //
+ static const float Y = 0.54076099395751953125f;
+
+ static const T P[] = {
+ 0.036454670944013329356512090082402429697L,
+ -0.066235835556476033710068679907798799959L,
+ -0.67492399795577182387312206593595565371L,
+ -1.4345555263962411429855341651960000166L,
+ -1.4894319559821365820516771951249649563L,
+ -0.87210277668067964629483299712322411566L,
+ -0.29602090537771744401524080430529369136L,
+ -0.0561832587517836908929331992218879676L,
+ -0.0053236785487328044334381502530383140443L,
+ -0.00018629360291358130461736386077971890789L,
+ -0.10164985672213178500790406939467614498e-6L,
+ 0.13680157145361387405588201461036338274e-8L
+ };
+ static const T Q[] = {
+ 1,
+ 4.9106336261005990534095838574132225599L,
+ 10.258804800866438510889341082793078432L,
+ 11.88588976846826108836629960537466889L,
+ 8.3455000546999704314454891036700998428L,
+ 3.6428823682421746343233362007194282703L,
+ 0.97465989807254572142266753052776132252L,
+ 0.15121052897097822172763084966793352524L,
+ 0.012017363555383555123769849654484594893L,
+ 0.0003583032812720649835431669893011257277L
+ };
+
+ T r = tools::evaluate_polynomial(P, zm1) / tools::evaluate_polynomial(Q, zm1);
+ T prefix = zm1 * zm2;
+
+ result += prefix * Y + prefix * r;
+ }
+ else if(z <= 1.625)
+ {
+ //
+ // Use the following form:
+ //
+ // lgamma(z) = (2-z)(1-z)(Y + R(2-z))
+ //
+ // where R(2-z) is a rational approximation optimised for
+ // low absolute error - as long as it's absolute error
+ // is small compared to the constant Y - then any rounding
+ // error in it's computation will get wiped out.
+ //
+ // R(2-z) has the following properties:
+ //
+ // Max error found at 128-bit long double precision 9.634e-36
+ // Maximum Deviation Found (approximation error) 1.538e-37
+ // Expected Error Term (theoretical error) 2.350e-38
+ //
+ static const float Y = 0.483787059783935546875f;
+
+ static const T P[] = {
+ -0.017977422421608624353488126610933005432L,
+ 0.18484528905298309555089509029244135703L,
+ -0.40401251514859546989565001431430884082L,
+ 0.40277179799147356461954182877921388182L,
+ -0.21993421441282936476709677700477598816L,
+ 0.069595742223850248095697771331107571011L,
+ -0.012681481427699686635516772923547347328L,
+ 0.0012489322866834830413292771335113136034L,
+ -0.57058739515423112045108068834668269608e-4L,
+ 0.8207548771933585614380644961342925976e-6L
+ };
+ static const T Q[] = {
+ 1,
+ -2.9629552288944259229543137757200262073L,
+ 3.7118380799042118987185957298964772755L,
+ -2.5569815272165399297600586376727357187L,
+ 1.0546764918220835097855665680632153367L,
+ -0.26574021300894401276478730940980810831L,
+ 0.03996289731752081380552901986471233462L,
+ -0.0033398680924544836817826046380586480873L,
+ 0.00013288854760548251757651556792598235735L,
+ -0.17194794958274081373243161848194745111e-5L
+ };
+ T r = zm2 * zm1;
+ T R = tools::evaluate_polynomial(P, 0.625 - zm1) / tools::evaluate_polynomial(Q, 0.625 - zm1);
+
+ result += r * Y + r * R;
+ }
+ else
+ {
+ //
+ // Same form as above.
+ //
+ // Max error found (at 128-bit long double precision) 1.831e-35
+ // Maximum Deviation Found (approximation error) 8.588e-36
+ // Expected Error Term (theoretical error) 1.458e-36
+ //
+ static const float Y = 0.443811893463134765625f;
+
+ static const T P[] = {
+ -0.021027558364667626231512090082402429494L,
+ 0.15128811104498736604523586803722368377L,
+ -0.26249631480066246699388544451126410278L,
+ 0.21148748610533489823742352180628489742L,
+ -0.093964130697489071999873506148104370633L,
+ 0.024292059227009051652542804957550866827L,
+ -0.0036284453226534839926304745756906117066L,
+ 0.0002939230129315195346843036254392485984L,
+ -0.11088589183158123733132268042570710338e-4L,
+ 0.13240510580220763969511741896361984162e-6L
+ };
+ static const T Q[] = {
+ 1,
+ -2.4240003754444040525462170802796471996L,
+ 2.4868383476933178722203278602342786002L,
+ -1.4047068395206343375520721509193698547L,
+ 0.47583809087867443858344765659065773369L,
+ -0.09865724264554556400463655444270700132L,
+ 0.012238223514176587501074150988445109735L,
+ -0.00084625068418239194670614419707491797097L,
+ 0.2796574430456237061420839429225710602e-4L,
+ -0.30202973883316730694433702165188835331e-6L
+ };
+ // (2 - x) * (1 - x) * (c + R(2 - x))
+ T r = zm2 * zm1;
+ T R = tools::evaluate_polynomial(P, -zm2) / tools::evaluate_polynomial(Q, -zm2);
+
+ result += r * Y + r * R;
+ }
+ }
+ return result;
+}
+template <class T, class Policy>
+T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<0>&, const Policy& pol)
+{
+ typedef typename lanczos::lanczos<T, Policy>::type L;
+ //
+ // No rational approximations are available because either
+ // T has no numeric_limits support (so we can't tell how
+ // many digits it has), or T has more digits than we know
+ // what to do with.... we do have a Lanczos approximation
+ // though, and that can be used to keep errors under control.
+ //
+ using namespace std; // for ADL of std names
+ T result = 0;
+ if(z < tools::epsilon<T>())
+ {
+ result = -log(z);
+ }
+ else if(z < 0.5)
+ {
+ // taking the log of tgamma reduces the error, no danger of overflow here:
+ result = log(gamma_imp(z, pol));
+ }
+ else if(z >= 3)
+ {
+ // taking the log of tgamma reduces the error, no danger of overflow here:
+ result = log(gamma_imp(z, pol));
+ }
+ else if(z >= 1.5)
+ {
+ // special case near 2:
+ T dz = zm2;
+ result = dz * log((z + L::g() - T(0.5)) / boost::math::constants::e<T>());
+ result += boost::math::log1p(dz / (L::g() + T(1.5)), pol) * T(1.5);
+ result += boost::math::log1p(L::lanczos_sum_near_2(dz), pol);
+ }
+ else
+ {
+ // special case near 1:
+ T dz = zm1;
+ result = dz * log((z + L::g() - T(0.5)) / boost::math::constants::e<T>());
+ result += boost::math::log1p(dz / (L::g() + T(0.5)), pol) / 2;
+ result += boost::math::log1p(L::lanczos_sum_near_1(dz), pol);
+ }
+ return result;
+}
+
+}}} // namespaces
+
+#endif // BOOST_MATH_SPECIAL_FUNCTIONS_DETAIL_LGAMMA_SMALL
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/detail/t_distribution_inv.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/detail/t_distribution_inv.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/detail/t_distribution_inv.hpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -72,7 +72,7 @@
using namespace std;
// Tail series expansion, see section 6 of Shaw's paper.
// w is calculated using Eq 60:
- T w = detail::tgamma_delta_ratio_imp(df / 2, constants::half<T>(), l)
+ T w = detail::tgamma_delta_ratio_imp(df / 2, constants::half<T>(), policy::policy<>())
* sqrt(df * constants::pi<T>()) * v;
// define some variables:
T np2 = df + 2;
@@ -126,7 +126,7 @@
//
// Start with Eq 56 of Shaw:
//
- T v = detail::tgamma_delta_ratio_imp(df / 2, constants::half<T>(), l)
+ T v = detail::tgamma_delta_ratio_imp(df / 2, constants::half<T>(), policy::policy<>())
* sqrt(df * constants::pi<T>()) * (u - constants::half<T>());
//
// Workspace for the polynomial coefficients:
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/erf.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/erf.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/erf.hpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -113,13 +113,13 @@
result = z * exp(-x);
result /= sqrt(boost::math::constants::pi<T>());
if(result != 0)
- result *= 2 * detail::lower_gamma_series(T(0.5f), x, boost::math::tools::digits<T>());
+ result *= 2 * detail::lower_gamma_series(T(0.5f), x, policy::policy<>());
}
else if(x < 1.1f)
{
// Compute Q:
invert = !invert;
- result = tgamma_small_upper_part(T(0.5f), x, l);
+ result = tgamma_small_upper_part(T(0.5f), x, policy::policy<>());
result /= sqrt(boost::math::constants::pi<T>());
}
else
@@ -774,6 +774,13 @@
tag_type()), BOOST_CURRENT_FUNCTION);
}
+template <class T, class Policy>
+inline typename tools::promote_args<T>::type erfc(T z, const Policy& pol)
+{
+ // Temporary Hack:
+ return erfc(z);
+}
+
} // namespace math
} // namespace boost
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/expm1.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/expm1.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/expm1.hpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -12,9 +12,10 @@
#include <boost/math/tools/config.hpp>
#include <boost/math/tools/series.hpp>
#include <boost/math/tools/precision.hpp>
-#include <boost/math/tools/evaluation_type.hpp>
+#include <boost/math/policy/error_handling.hpp>
#include <boost/math/tools/rational.hpp>
#include <boost/math/special_functions/math_fwd.hpp>
+#include <boost/mpl/less_equal.hpp>
#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
# include <boost/static_assert.hpp>
@@ -64,8 +65,8 @@
//
// This version uses a Taylor series expansion for 0.5 > |x| > epsilon.
//
-template <class T>
-T expm1_imp(T x, const mpl::int_<0>&)
+template <class T, class Policy>
+T expm1_imp(T x, const mpl::int_<0>&, const Policy& pol)
{
using namespace std;
@@ -77,17 +78,17 @@
detail::expm1_series<T> s(x);
boost::uintmax_t max_iter = BOOST_MATH_MAX_ITER;
#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
- T result = tools::sum_series(s, tools::digits<T>() + 2, max_iter);
+ T result = tools::sum_series(s, policy::digits<T, Policy>(), max_iter);
#else
T zero = 0;
- T result = tools::sum_series(s, tools::digits<T>() + 2, max_iter, zero);
+ T result = tools::sum_series(s, policy::digits<T, Policy>(), max_iter, zero);
#endif
- tools::check_series_iterations(BOOST_CURRENT_FUNCTION, max_iter);
+ policy::check_series_iterations(BOOST_CURRENT_FUNCTION, max_iter, pol);
return result;
}
-template <class T>
-T expm1_imp(T x, const mpl::int_<53>&)
+template <class T, class P>
+T expm1_imp(T x, const mpl::int_<53>&, const P&)
{
using namespace std;
@@ -105,8 +106,8 @@
return result;
}
-template <class T>
-T expm1_imp(T x, const mpl::int_<64>&)
+template <class T, class P>
+T expm1_imp(T x, const mpl::int_<64>&, const P&)
{
using namespace std;
@@ -140,8 +141,8 @@
return result;
}
-template <class T>
-T expm1_imp(T x, const mpl::int_<113>&)
+template <class T, class P>
+T expm1_imp(T x, const mpl::int_<113>&, const P&)
{
using namespace std;
@@ -184,23 +185,24 @@
} // namespace detail
-template <class T>
-inline typename tools::promote_args<T>::type expm1(T x)
+template <class T, class Policy>
+inline typename tools::promote_args<T>::type expm1(T x, const Policy& pol)
{
typedef typename tools::promote_args<T>::type result_type;
- typedef typename tools::evaluation<result_type>::type value_type;
+ typedef typename policy::evaluation<result_type, Policy>::type value_type;
+ typedef typename policy::precision<result_type, Policy>::type precision_type;
typedef typename mpl::if_c<
::std::numeric_limits<result_type>::is_specialized == 0,
mpl::int_<0>, // no numeric_limits, use generic solution
- typename mpl::if_c<
- ::std::numeric_limits<result_type>::digits <= 53,
+ typename mpl::if_<
+ typename mpl::less_equal<precision_type, mpl::int_<53> >::type,
mpl::int_<53>, // double
- typename mpl::if_c<
- ::std::numeric_limits<result_type>::digits <= 64,
+ typename mpl::if_<
+ typename mpl::less_equal<precision_type, mpl::int_<64> >::type,
mpl::int_<64>, // 80-bit long double
- typename mpl::if_c<
- ::std::numeric_limits<result_type>::digits <= 113,
+ typename mpl::if_<
+ typename mpl::less_equal<precision_type, mpl::int_<113> >::type,
mpl::int_<113>, // 128-bit long double
mpl::int_<0> // too many bits, use generic version.
>::type
@@ -210,7 +212,30 @@
return tools::checked_narrowing_cast<result_type>(detail::expm1_imp(
static_cast<value_type>(x),
- tag_type()), BOOST_CURRENT_FUNCTION);
+ tag_type(), pol), BOOST_CURRENT_FUNCTION);
+}
+
+#ifdef expm1
+# ifndef BOOST_HAS_expm1
+# define BOOST_HAS_expm1
+# endif
+# undef expm1
+#endif
+
+#ifdef BOOST_HAS_EXPM1
+# if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901)
+inline float expm1(float x, const policy::policy<>&){ return ::expm1f(x); }
+inline long double expm1(long double x, const policy::policy<>&){ return ::expm1l(x); }
+#else
+inline float expm1(float x, const policy::policy<>&){ return ::expm1(x); }
+#endif
+inline double expm1(double x, const policy::policy<>&){ return ::expm1(x); }
+#endif
+
+template <class T>
+inline typename tools::promote_args<T>::type expm1(T x)
+{
+ return expm1(x, policy::policy<>());
}
#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x564))
@@ -228,23 +253,6 @@
}
#endif
-#ifdef expm1
-# ifndef BOOST_HAS_expm1
-# define BOOST_HAS_expm1
-# endif
-# undef expm1
-#endif
-
-#ifdef BOOST_HAS_EXPM1
-# if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901)
-inline float expm1(float x){ return ::expm1f(x); }
-inline long double expm1(long double x){ return ::expm1l(x); }
-#else
-inline float expm1(float x){ return ::expm1(x); }
-#endif
-inline double expm1(double x){ return ::expm1(x); }
-#endif
-
} // namespace math
} // namespace boost
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/gamma.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/gamma.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/gamma.hpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -33,8 +33,11 @@
#include <boost/math/special_functions/fpclassify.hpp>
#include <boost/math/special_functions/detail/igamma_large.hpp>
#include <boost/math/special_functions/detail/unchecked_factorial.hpp>
+#include <boost/math/special_functions/detail/lgamma_small.hpp>
#include <boost/type_traits/is_convertible.hpp>
#include <boost/assert.hpp>
+#include <boost/mpl/greater.hpp>
+#include <boost/mpl/equal_to.hpp>
#include <cmath>
#include <algorithm>
@@ -111,8 +114,8 @@
//
// tgamma(z), with Lanczos support:
//
-template <class T, class L>
-T gamma_imp(T z, const L& l)
+template <class T, class Policy, class L>
+T gamma_imp(T z, const Policy& pol, const L& l)
{
using namespace std;
@@ -126,19 +129,20 @@
b = true;
}
#endif
+ static const char* function = "boost::math::tgamma<%1%>(%1%)";
if((z <= 0) && (floor(z) == z))
- return tools::pole_error<T>(BOOST_CURRENT_FUNCTION, "Evaluation of tgamma at a negative integer %1%.", z);
+ return policy::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
if(z <= -20)
{
- result = gamma_imp(-z, l) * sinpx(z);
+ result = gamma_imp(-z, pol, l) * sinpx(z);
if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
- return tools::overflow_error<T>(BOOST_CURRENT_FUNCTION, "Result of tgamma is too large to represent.");
+ return policy::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
result = -boost::math::constants::pi<T>() / result;
if(result == 0)
- return tools::underflow_error<T>(BOOST_CURRENT_FUNCTION, "Result of tgamma is too small to represent.");
+ return policy::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
if(boost::math::fpclassify(result) == FP_SUBNORMAL)
- return tools::denorm_error<T>(result, BOOST_CURRENT_FUNCTION, "Result of tgamma is denormalized.");
+ return policy::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);
return result;
}
@@ -160,11 +164,11 @@
// we're going to overflow unless this is done with care:
T zgh = (z + static_cast<T>(L::g()) - boost::math::constants::half<T>());
if(log(zgh) * z / 2 > tools::log_max_value<T>())
- return tools::overflow_error<T>(BOOST_CURRENT_FUNCTION, "Result of tgamma is too large to represent.");
+ return policy::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
T hp = pow(zgh, (z / 2) - T(0.25));
result *= hp / exp(zgh);
if(tools::max_value<T>() / hp < result)
- return tools::overflow_error<T>(BOOST_CURRENT_FUNCTION, "Result of tgamma is too large to represent.");
+ return policy::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
result *= hp;
}
else
@@ -176,491 +180,10 @@
return result;
}
//
-// lgamma for small arguments:
-//
-template <class T, class L>
-T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<64>&, const L& /* l */)
-{
- // This version uses rational approximations for small
- // values of z accurate enough for 64-bit mantissas
- // (80-bit long doubles), works well for 53-bit doubles as well.
- // L is only used to select the Lanczos function.
-
- using namespace std; // for ADL of std names
- T result = 0;
- if(z < tools::epsilon<T>())
- {
- result = -log(z);
- }
- else if((zm1 == 0) || (zm2 == 0))
- {
- // nothing to do, result is zero....
- }
- else if(z > 2)
- {
- //
- // Begin by performing argument reduction until
- // z is in [2,3):
- //
- if(z >= 3)
- {
- do
- {
- z -= 1;
- zm2 -= 1;
- result += log(z);
- }while(z >= 3);
- // Update zm2, we need it below:
- zm2 = z - 2;
- }
-
- //
- // Use the following form:
- //
- // lgamma(z) = (z-2)(z+1)(Y + R(z-2))
- //
- // where R(z-2) is a rational approximation optimised for
- // low absolute error - as long as it's absolute error
- // is small compared to the constant Y - then any rounding
- // error in it's computation will get wiped out.
- //
- // R(z-2) has the following properties:
- //
- // At double: Max error found: 4.231e-18
- // At long double: Max error found: 1.987e-21
- // Maximum Deviation Found (approximation error): 5.900e-24
- //
- static const T P[] = {
- static_cast<T>(-0.180355685678449379109e-1L),
- static_cast<T>(0.25126649619989678683e-1L),
- static_cast<T>(0.494103151567532234274e-1L),
- static_cast<T>(0.172491608709613993966e-1L),
- static_cast<T>(-0.259453563205438108893e-3L),
- static_cast<T>(-0.541009869215204396339e-3L),
- static_cast<T>(-0.324588649825948492091e-4L)
- };
- static const T Q[] = {
- static_cast<T>(0.1e1),
- static_cast<T>(0.196202987197795200688e1L),
- static_cast<T>(0.148019669424231326694e1L),
- static_cast<T>(0.541391432071720958364e0L),
- static_cast<T>(0.988504251128010129477e-1L),
- static_cast<T>(0.82130967464889339326e-2L),
- static_cast<T>(0.224936291922115757597e-3L),
- static_cast<T>(-0.223352763208617092964e-6L)
- };
-
- static const float Y = 0.158963680267333984375e0f;
-
- T r = zm2 * (z + 1);
- T R = tools::evaluate_polynomial(P, zm2);
- R /= tools::evaluate_polynomial(Q, zm2);
-
- result += r * Y + r * R;
- }
- else
- {
- //
- // If z is less than 1 use recurrance to shift to
- // z in the interval [1,2]:
- //
- if(z < 1)
- {
- result += -log(z);
- zm2 = zm1;
- zm1 = z;
- z += 1;
- }
- //
- // Two approximations, on for z in [1,1.5] and
- // one for z in [1.5,2]:
- //
- if(z <= 1.5)
- {
- //
- // Use the following form:
- //
- // lgamma(z) = (z-1)(z-2)(Y + R(z-1))
- //
- // where R(z-1) is a rational approximation optimised for
- // low absolute error - as long as it's absolute error
- // is small compared to the constant Y - then any rounding
- // error in it's computation will get wiped out.
- //
- // R(z-1) has the following properties:
- //
- // At double precision: Max error found: 1.230011e-17
- // At 80-bit long double precision: Max error found: 5.631355e-21
- // Maximum Deviation Found: 3.139e-021
- // Expected Error Term: 3.139e-021
-
- //
- static const float Y = 0.52815341949462890625f;
-
- static const T P[] = {
- static_cast<T>(0.490622454069039543534e-1L),
- static_cast<T>(-0.969117530159521214579e-1L),
- static_cast<T>(-0.414983358359495381969e0L),
- static_cast<T>(-0.406567124211938417342e0L),
- static_cast<T>(-0.158413586390692192217e0L),
- static_cast<T>(-0.240149820648571559892e-1L),
- static_cast<T>(-0.100346687696279557415e-2L)
- };
- static const T Q[] = {
- static_cast<T>(0.1e1L),
- static_cast<T>(0.302349829846463038743e1L),
- static_cast<T>(0.348739585360723852576e1L),
- static_cast<T>(0.191415588274426679201e1L),
- static_cast<T>(0.507137738614363510846e0L),
- static_cast<T>(0.577039722690451849648e-1L),
- static_cast<T>(0.195768102601107189171e-2L)
- };
-
- T r = tools::evaluate_polynomial(P, zm1) / tools::evaluate_polynomial(Q, zm1);
- T prefix = zm1 * zm2;
-
- result += prefix * Y + prefix * r;
- }
- else
- {
- //
- // Use the following form:
- //
- // lgamma(z) = (2-z)(1-z)(Y + R(2-z))
- //
- // where R(2-z) is a rational approximation optimised for
- // low absolute error - as long as it's absolute error
- // is small compared to the constant Y - then any rounding
- // error in it's computation will get wiped out.
- //
- // R(2-z) has the following properties:
- //
- // At double precision, max error found: 1.797565e-17
- // At 80-bit long double precision, max error found: 9.306419e-21
- // Maximum Deviation Found: 2.151e-021
- // Expected Error Term: 2.150e-021
- //
- static const float Y = 0.452017307281494140625f;
-
- static const T P[] = {
- static_cast<T>(-0.292329721830270012337e-1L),
- static_cast<T>(0.144216267757192309184e0L),
- static_cast<T>(-0.142440390738631274135e0L),
- static_cast<T>(0.542809694055053558157e-1L),
- static_cast<T>(-0.850535976868336437746e-2L),
- static_cast<T>(0.431171342679297331241e-3L)
- };
- static const T Q[] = {
- static_cast<T>(0.1e1),
- static_cast<T>(-0.150169356054485044494e1L),
- static_cast<T>(0.846973248876495016101e0L),
- static_cast<T>(-0.220095151814995745555e0L),
- static_cast<T>(0.25582797155975869989e-1L),
- static_cast<T>(-0.100666795539143372762e-2L),
- static_cast<T>(-0.827193521891290553639e-6L)
- };
- T r = zm2 * zm1;
- T R = tools::evaluate_polynomial(P, -zm2) / tools::evaluate_polynomial(Q, -zm2);
-
- result += r * Y + r * R;
- }
- }
- return result;
-}
-template <class T, class L>
-T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<113>&, const L& /* l */)
-{
- //
- // This version uses rational approximations for small
- // values of z accurate enough for 113-bit mantissas
- // (128-bit long doubles).
- //
- using namespace std; // for ADL of std names
- T result = 0;
- if(z < tools::epsilon<T>())
- {
- result = -log(z);
- }
- else if((zm1 == 0) || (zm2 == 0))
- {
- // nothing to do, result is zero....
- }
- else if(z > 2)
- {
- //
- // Begin by performing argument reduction until
- // z is in [2,3):
- //
- if(z >= 3)
- {
- do
- {
- z -= 1;
- result += log(z);
- }while(z >= 3);
- zm2 = z - 2;
- }
-
- //
- // Use the following form:
- //
- // lgamma(z) = (z-2)(z+1)(Y + R(z-2))
- //
- // where R(z-2) is a rational approximation optimised for
- // low absolute error - as long as it's absolute error
- // is small compared to the constant Y - then any rounding
- // error in it's computation will get wiped out.
- //
- // Maximum Deviation Found (approximation error) 3.73e-37
-
- static const T P[] = {
- -0.018035568567844937910504030027467476655L,
- 0.013841458273109517271750705401202404195L,
- 0.062031842739486600078866923383017722399L,
- 0.052518418329052161202007865149435256093L,
- 0.01881718142472784129191838493267755758L,
- 0.0025104830367021839316463675028524702846L,
- -0.00021043176101831873281848891452678568311L,
- -0.00010249622350908722793327719494037981166L,
- -0.11381479670982006841716879074288176994e-4L,
- -0.49999811718089980992888533630523892389e-6L,
- -0.70529798686542184668416911331718963364e-8L
- };
- static const T Q[] = {
- 1L,
- 2.5877485070422317542808137697939233685L,
- 2.8797959228352591788629602533153837126L,
- 1.8030885955284082026405495275461180977L,
- 0.69774331297747390169238306148355428436L,
- 0.17261566063277623942044077039756583802L,
- 0.02729301254544230229429621192443000121L,
- 0.0026776425891195270663133581960016620433L,
- 0.00015244249160486584591370355730402168106L,
- 0.43997034032479866020546814475414346627e-5L,
- 0.46295080708455613044541885534408170934e-7L,
- -0.93326638207459533682980757982834180952e-11L,
- 0.42316456553164995177177407325292867513e-13L
- };
-
- T R = tools::evaluate_polynomial(P, zm2);
- R /= tools::evaluate_polynomial(Q, zm2);
-
- static const float Y = 0.158963680267333984375F;
-
- T r = zm2 * (z + 1);
-
- result += r * Y + r * R;
- }
- else
- {
- //
- // If z is less than 1 use recurrance to shift to
- // z in the interval [1,2]:
- //
- if(z < 1)
- {
- result += -log(z);
- zm2 = zm1;
- zm1 = z;
- z += 1;
- }
- //
- // Three approximations, on for z in [1,1.35], [1.35,1.625] and [1.625,1]
- //
- if(z <= 1.35)
- {
- //
- // Use the following form:
- //
- // lgamma(z) = (z-1)(z-2)(Y + R(z-1))
- //
- // where R(z-1) is a rational approximation optimised for
- // low absolute error - as long as it's absolute error
- // is small compared to the constant Y - then any rounding
- // error in it's computation will get wiped out.
- //
- // R(z-1) has the following properties:
- //
- // Maximum Deviation Found (approximation error) 1.659e-36
- // Expected Error Term (theoretical error) 1.343e-36
- // Max error found at 128-bit long double precision 1.007e-35
- //
- static const float Y = 0.54076099395751953125f;
-
- static const T P[] = {
- 0.036454670944013329356512090082402429697L,
- -0.066235835556476033710068679907798799959L,
- -0.67492399795577182387312206593595565371L,
- -1.4345555263962411429855341651960000166L,
- -1.4894319559821365820516771951249649563L,
- -0.87210277668067964629483299712322411566L,
- -0.29602090537771744401524080430529369136L,
- -0.0561832587517836908929331992218879676L,
- -0.0053236785487328044334381502530383140443L,
- -0.00018629360291358130461736386077971890789L,
- -0.10164985672213178500790406939467614498e-6L,
- 0.13680157145361387405588201461036338274e-8L
- };
- static const T Q[] = {
- 1,
- 4.9106336261005990534095838574132225599L,
- 10.258804800866438510889341082793078432L,
- 11.88588976846826108836629960537466889L,
- 8.3455000546999704314454891036700998428L,
- 3.6428823682421746343233362007194282703L,
- 0.97465989807254572142266753052776132252L,
- 0.15121052897097822172763084966793352524L,
- 0.012017363555383555123769849654484594893L,
- 0.0003583032812720649835431669893011257277L
- };
-
- T r = tools::evaluate_polynomial(P, zm1) / tools::evaluate_polynomial(Q, zm1);
- T prefix = zm1 * zm2;
-
- result += prefix * Y + prefix * r;
- }
- else if(z <= 1.625)
- {
- //
- // Use the following form:
- //
- // lgamma(z) = (2-z)(1-z)(Y + R(2-z))
- //
- // where R(2-z) is a rational approximation optimised for
- // low absolute error - as long as it's absolute error
- // is small compared to the constant Y - then any rounding
- // error in it's computation will get wiped out.
- //
- // R(2-z) has the following properties:
- //
- // Max error found at 128-bit long double precision 9.634e-36
- // Maximum Deviation Found (approximation error) 1.538e-37
- // Expected Error Term (theoretical error) 2.350e-38
- //
- static const float Y = 0.483787059783935546875f;
-
- static const T P[] = {
- -0.017977422421608624353488126610933005432L,
- 0.18484528905298309555089509029244135703L,
- -0.40401251514859546989565001431430884082L,
- 0.40277179799147356461954182877921388182L,
- -0.21993421441282936476709677700477598816L,
- 0.069595742223850248095697771331107571011L,
- -0.012681481427699686635516772923547347328L,
- 0.0012489322866834830413292771335113136034L,
- -0.57058739515423112045108068834668269608e-4L,
- 0.8207548771933585614380644961342925976e-6L
- };
- static const T Q[] = {
- 1,
- -2.9629552288944259229543137757200262073L,
- 3.7118380799042118987185957298964772755L,
- -2.5569815272165399297600586376727357187L,
- 1.0546764918220835097855665680632153367L,
- -0.26574021300894401276478730940980810831L,
- 0.03996289731752081380552901986471233462L,
- -0.0033398680924544836817826046380586480873L,
- 0.00013288854760548251757651556792598235735L,
- -0.17194794958274081373243161848194745111e-5L
- };
- T r = zm2 * zm1;
- T R = tools::evaluate_polynomial(P, 0.625 - zm1) / tools::evaluate_polynomial(Q, 0.625 - zm1);
-
- result += r * Y + r * R;
- }
- else
- {
- //
- // Same form as above.
- //
- // Max error found (at 128-bit long double precision) 1.831e-35
- // Maximum Deviation Found (approximation error) 8.588e-36
- // Expected Error Term (theoretical error) 1.458e-36
- //
- static const float Y = 0.443811893463134765625f;
-
- static const T P[] = {
- -0.021027558364667626231512090082402429494L,
- 0.15128811104498736604523586803722368377L,
- -0.26249631480066246699388544451126410278L,
- 0.21148748610533489823742352180628489742L,
- -0.093964130697489071999873506148104370633L,
- 0.024292059227009051652542804957550866827L,
- -0.0036284453226534839926304745756906117066L,
- 0.0002939230129315195346843036254392485984L,
- -0.11088589183158123733132268042570710338e-4L,
- 0.13240510580220763969511741896361984162e-6L
- };
- static const T Q[] = {
- 1,
- -2.4240003754444040525462170802796471996L,
- 2.4868383476933178722203278602342786002L,
- -1.4047068395206343375520721509193698547L,
- 0.47583809087867443858344765659065773369L,
- -0.09865724264554556400463655444270700132L,
- 0.012238223514176587501074150988445109735L,
- -0.00084625068418239194670614419707491797097L,
- 0.2796574430456237061420839429225710602e-4L,
- -0.30202973883316730694433702165188835331e-6L
- };
- // (2 - x) * (1 - x) * (c + R(2 - x))
- T r = zm2 * zm1;
- T R = tools::evaluate_polynomial(P, -zm2) / tools::evaluate_polynomial(Q, -zm2);
-
- result += r * Y + r * R;
- }
- }
- return result;
-}
-template <class T, class L>
-T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<0>&, const L& l)
-{
- //
- // No rational approximations are available because either
- // T has no numeric_limits support (so we can't tell how
- // many digits it has), or T has more digits than we know
- // what to do with.... we do have a Lanczos approximation
- // though, and that can be used to keep errors under control.
- //
- using namespace std; // for ADL of std names
- T result = 0;
- if(z < tools::epsilon<T>())
- {
- result = -log(z);
- }
- else if(z < 0.5)
- {
- // taking the log of tgamma reduces the error, no danger of overflow here:
- result = log(gamma_imp(z, l));
- }
- else if(z >= 3)
- {
- // taking the log of tgamma reduces the error, no danger of overflow here:
- result = log(gamma_imp(z, l));
- }
- else if(z >= 1.5)
- {
- // special case near 2:
- T dz = zm2;
- result = dz * log((z + L::g() - T(0.5)) / boost::math::constants::e<T>());
- result += boost::math::log1p(dz / (L::g() + T(1.5))) * T(1.5);
- result += boost::math::log1p(L::lanczos_sum_near_2(dz));
- }
- else
- {
- // special case near 1:
- T dz = zm1;
- result = dz * log((z + L::g() - T(0.5)) / boost::math::constants::e<T>());
- result += boost::math::log1p(dz / (L::g() + T(0.5))) / 2;
- result += boost::math::log1p(L::lanczos_sum_near_1(dz));
- }
- return result;
-}
-//
// lgamma(z) with Lanczos support:
//
-template <class T, class L>
-T lgamma_imp(T z, const L& l, int* sign = 0)
+template <class T, class Policy, class L>
+T lgamma_imp(T z, const Policy& pol, const L& l, int* sign = 0)
{
#ifdef BOOST_MATH_INSTRUMENT
static bool b = false;
@@ -673,13 +196,15 @@
using namespace std;
+ static const char* function = "boost::math::lgamma<%1%>(%1%)";
+
T result = 0;
int sresult = 1;
if(z <= 0)
{
// reflection formula:
if(floor(z) == z)
- return tools::pole_error<T>(BOOST_CURRENT_FUNCTION, "Evaluation of lgamma at a negative integer %1%.", z);
+ return policy::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol);
T t = sinpx(z);
z = -z;
@@ -691,23 +216,24 @@
{
sresult = -sresult;
}
- result = log(boost::math::constants::pi<T>()) - lgamma_imp(z, l) - log(t);
+ result = log(boost::math::constants::pi<T>()) - lgamma_imp(z, pol, l) - log(t);
}
else if(z < 15)
{
- typedef typename mpl::if_c<
- (std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::digits <= 64)),
+ typedef typename policy::precision<T, Policy>::type precision_type;
+ typedef typename mpl::if_<
+ mpl::less_equal<precision_type, mpl::int_<64> >,
mpl::int_<64>,
- typename mpl::if_c<
- (std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::digits <= 113)),
+ typename mpl::if_<
+ mpl::less_equal<precision_type, mpl::int_<113> >,
mpl::int_<113>, mpl::int_<0> >::type
>::type tag_type;
- result = lgamma_small_imp(z, z - 1, z - 2, tag_type(), l);
+ result = lgamma_small_imp(z, z - 1, z - 2, tag_type(), pol);
}
else if((z >= 3) && (z < 100))
{
// taking the log of tgamma reduces the error, no danger of overflow here:
- result = log(gamma_imp(z, l));
+ result = log(gamma_imp(z, pol, l));
}
else
{
@@ -777,21 +303,22 @@
}
};
-template <class T>
-inline T lower_gamma_series(T a, T z, int bits)
+template <class T, class Policy>
+inline T lower_gamma_series(T a, T z, const Policy& pol)
{
// Multiply result by ((z^a) * (e^-z) / a) to get the full
// lower incomplete integral. Then divide by tgamma(a)
// to get the normalised value.
lower_incomplete_gamma_series<T> s(a, z);
boost::uintmax_t max_iter = BOOST_MATH_MAX_ITER;
+ int bits = policy::digits<T, Policy>();
#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
T zero = 0;
T result = boost::math::tools::sum_series(s, bits, max_iter, zero);
#else
T result = boost::math::tools::sum_series(s, bits, max_iter);
#endif
- tools::check_series_iterations(BOOST_CURRENT_FUNCTION, max_iter);
+ policy::check_series_iterations("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol);
return result;
}
@@ -799,22 +326,23 @@
// Fully generic tgamma and lgamma use the incomplete partial
// sums added together:
//
-template <class T>
-T gamma_imp(T z, const lanczos::undefined_lanczos& l)
+template <class T, class Policy>
+T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l)
{
+ static const char* function = "boost::math::tgamma<%1%>(%1%)";
using namespace std;
if((z <= 0) && (floor(z) == z))
- return tools::pole_error<T>(BOOST_CURRENT_FUNCTION, "Evaluation of tgamma at a negative integer %1%.", z);
+ return policy::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
if(z <= -20)
{
- T result = gamma_imp(-z, l) * sinpx(z);
+ T result = gamma_imp(-z, pol, l) * sinpx(z);
if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
- return tools::overflow_error<T>(BOOST_CURRENT_FUNCTION, "Result of tgamma is too large to represent.");
+ return policy::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
result = -boost::math::constants::pi<T>() / result;
if(result == 0)
- return tools::underflow_error<T>(BOOST_CURRENT_FUNCTION, "Result of tgamma is too small to represent.");
+ return policy::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
if(boost::math::fpclassify(result) == FP_SUBNORMAL)
- return tools::denorm_error<T>(result, BOOST_CURRENT_FUNCTION, "Result of tgamma is denormalized.");
+ return policy::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);
return result;
}
//
@@ -835,26 +363,27 @@
else
{
prefix = prefix * pow(z / boost::math::constants::e<T>(), z);
- T sum = detail::lower_gamma_series(z, z, ::boost::math::tools::digits<T>()) / z;
- sum += detail::upper_gamma_fraction(z, z, ::boost::math::tools::digits<T>());
+ T sum = detail::lower_gamma_series(z, z, pol) / z;
+ sum += detail::upper_gamma_fraction(z, z, ::boost::math::policy::digits<T, Policy>());
if(fabs(tools::max_value<T>() / prefix) < fabs(sum))
- return tools::overflow_error<T>(BOOST_CURRENT_FUNCTION, "Result of tgamma is too large to represent.");
+ return policy::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
return sum * prefix;
}
return prefix;
}
-template <class T>
-T lgamma_imp(T z, const lanczos::undefined_lanczos&, int*sign)
+template <class T, class Policy>
+T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l, int*sign)
{
using namespace std;
+ static const char* function = "boost::math::lgamma<%1%>(%1%)";
T result = 0;
int sresult = 1;
if(z <= 0)
{
if(floor(z) == z)
- return tools::pole_error<T>(BOOST_CURRENT_FUNCTION, "Evaluation of tgamma at a negative integer %1%.", z);
+ return policy::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
T t = detail::sinpx(z);
z = -z;
if(t < 0)
@@ -865,14 +394,14 @@
{
sresult = -sresult;
}
- result = log(boost::math::constants::pi<T>()) - lgamma(z) - log(t);
+ result = log(boost::math::constants::pi<T>()) - lgamma_imp(z, pol, l, 0) - log(t);
}
else if((z != 1) && (z != 2))
{
T limit = (std::max)(z+1, T(10));
T prefix = z * log(limit) - limit;
- T sum = detail::lower_gamma_series(z, limit, ::boost::math::tools::digits<T>()) / z;
- sum += detail::upper_gamma_fraction(z, limit, ::boost::math::tools::digits<T>());
+ T sum = detail::lower_gamma_series(z, limit, pol) / z;
+ sum += detail::upper_gamma_fraction(z, limit, ::boost::math::policy::digits<T, Policy>());
result = log(sum) + prefix;
}
if(sign)
@@ -883,24 +412,32 @@
// This helper calculates tgamma(dz+1)-1 without cancellation errors,
// used by the upper incomplete gamma with z < 1:
//
-template <class T, class Tag, class L>
-T tgammap1m1_imp(T dz, Tag const& tag, const L& l)
+template <class T, class Policy, class L>
+T tgammap1m1_imp(T dz, Policy const& pol, const L& l)
{
using namespace std;
+ typedef typename mpl::if_<
+ mpl::less_equal<typename policy::precision<T,Policy>::type, mpl::int_<64> >,
+ mpl::int_<64>,
+ typename mpl::if_<
+ mpl::less_equal<typename policy::precision<T,Policy>::type, mpl::int_<113> >,
+ mpl::int_<113>, mpl::int_<0> >::type
+ >::type tag_type;
+
T result;
if(dz < 0)
{
if(dz < -0.5)
{
// Best method is simply to subtract 1 from tgamma:
- result = boost::math::tgamma(1+dz) - 1;
+ result = boost::math::tgamma(1+dz, pol) - 1;
}
else
{
// Use expm1 on lgamma:
- result = boost::math::expm1(-boost::math::log1p(dz)
- + lgamma_small_imp(dz+2, dz + 1, dz, tag, l));
+ result = boost::math::expm1(-boost::math::log1p(dz, pol)
+ + lgamma_small_imp(dz+2, dz + 1, dz, tag_type(), pol));
}
}
else
@@ -908,20 +445,20 @@
if(dz < 2)
{
// Use expm1 on lgamma:
- result = boost::math::expm1(lgamma_small_imp(dz+1, dz, dz-1, tag, l));
+ result = boost::math::expm1(lgamma_small_imp(dz+1, dz, dz-1, tag_type(), pol), pol);
}
else
{
// Best method is simply to subtract 1 from tgamma:
- result = boost::math::tgamma(1+dz) - 1;
+ result = boost::math::tgamma(1+dz, pol) - 1;
}
}
return result;
}
-template <class T, class Tag>
-inline T tgammap1m1_imp(T dz, Tag const& /*tag*/,
+template <class T, class Policy>
+inline T tgammap1m1_imp(T dz, Policy const& pol,
const ::boost::math::lanczos::undefined_lanczos& l)
{
using namespace std; // ADL of std names
@@ -930,7 +467,7 @@
// algebra isn't easy for the general case....
// Start by subracting 1 from tgamma:
//
- T result = gamma_imp(1 + dz, l) - 1;
+ T result = gamma_imp(1 + dz, pol, l) - 1;
//
// Test the level of cancellation error observed: we loose one bit
// for each power of 2 the result is less than 1. If we would get
@@ -939,25 +476,11 @@
//
if((dz > -0.5) && (dz < 2) && (ldexp(1.0, boost::math::tools::digits<T>()) * fabs(result) < 1e34))
{
- result = tgammap1m1_imp(dz, mpl::int_<113>(), boost::math::lanczos::lanczos24m113());
+ result = tgammap1m1_imp(dz, pol, boost::math::lanczos::lanczos24m113());
}
return result;
}
-template <class T, class L>
-inline T tgammap1m1_imp(T dz, const L& l)
-{
- // Simple forwarding function.
- typedef typename mpl::if_c<
- (std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::digits <= 64)),
- mpl::int_<64>,
- typename mpl::if_c<
- (std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::digits <= 113)),
- mpl::int_<113>, mpl::int_<0> >::type
- >::type tag_type;
- return tgammap1m1_imp(dz, tag_type(), l);
-}
-
//
// Series representation for upper fraction when z is small:
//
@@ -985,8 +508,8 @@
// calculate power term prefix (z^a)(e^-z) used in the non-normalised
// incomplete gammas:
//
-template <class T>
-T full_igamma_prefix(T a, T z)
+template <class T, class Policy>
+T full_igamma_prefix(T a, T z, const Policy& pol)
{
using namespace std;
@@ -1028,34 +551,16 @@
// rather than before it...
//
if(boost::math::fpclassify(prefix) == FP_INFINITE)
- tools::overflow_error<T>(BOOST_CURRENT_FUNCTION, "Result of incomplete gamma function is too large to represent.");
+ policy::raise_overflow_error<T>("boost::math::detail::full_igamma_prefix<%1%>(%1%, %1%)", "Result of incomplete gamma function is too large to represent.", pol);
return prefix;
}
//
-// Helper to compute log(1+x)-x:
-//
-template <class T>
-inline T log1pmx(T x)
-{
- boost::math::detail::log1p_series<T> s(x);
- s();
- boost::uintmax_t max_iter = BOOST_MATH_MAX_ITER;
-#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
- T zero = 0;
- T result = boost::math::tools::sum_series(s, tools::digits<T>() + 2, max_iter, zero);
-#else
- T result = boost::math::tools::sum_series(s, tools::digits<T>() + 2, max_iter);
-#endif
- tools::check_series_iterations(BOOST_CURRENT_FUNCTION, max_iter);
- return result;
-}
-//
// Compute (z^a)(e^-z)/tgamma(a)
// most if the error occurs in this function:
//
-template <class T, class L>
-T regularised_gamma_prefix(T a, T z, const L& l)
+template <class T, class Policy, class L>
+T regularised_gamma_prefix(T a, T z, const Policy& pol, const L& l)
{
using namespace std;
T agh = a + static_cast<T>(L::g()) - T(0.5);
@@ -1076,19 +581,19 @@
if(z <= tools::log_min_value<T>())
{
// Oh dear, have to use logs, should be free of cancellation errors though:
- return exp(a * log(z) - z - lgamma_imp(a, l));
+ return exp(a * log(z) - z - lgamma_imp(a, pol, l));
}
else
{
// direct calculation, no danger of overflow as gamma(a) < 1/a
// for small a.
- return pow(z, a) * exp(-z) / gamma_imp(a, l);
+ return pow(z, a) * exp(-z) / gamma_imp(a, pol, l);
}
}
else if((fabs(d*d*a) <= 100) && (a > 150))
{
// special case for large a and a ~ z.
- prefix = a * log1pmx(d) + z * static_cast<T>(0.5 - L::g()) / agh;
+ prefix = a * log1pmx(d, pol) + z * static_cast<T>(0.5 - L::g()) / agh;
prefix = exp(prefix);
}
else
@@ -1136,14 +641,14 @@
//
// And again, without Lanczos support:
//
-template <class T>
-T regularised_gamma_prefix(T a, T z, const lanczos::undefined_lanczos&)
+template <class T, class Policy>
+T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined_lanczos&)
{
using namespace std;
T limit = (std::max)(T(10), a);
- T sum = detail::lower_gamma_series(a, limit, ::boost::math::tools::digits<T>()) / a;
- sum += detail::upper_gamma_fraction(a, limit, ::boost::math::tools::digits<T>());
+ T sum = detail::lower_gamma_series(a, limit, pol) / a;
+ sum += detail::upper_gamma_fraction(a, limit, ::boost::math::policy::digits<T, Policy>());
if(a < 10)
{
@@ -1184,24 +689,24 @@
//
// Upper gamma fraction for very small a:
//
-template <class T, class L>
-inline T tgamma_small_upper_part(T a, T x, const L& l)
+template <class T, class Policy>
+inline T tgamma_small_upper_part(T a, T x, const Policy& pol)
{
using namespace std; // ADL of std functions.
//
// Compute the full upper fraction (Q) when a is very small:
//
- T result = tgammap1m1_imp(a, l) - powm1(x, a);
+ T result = tgamma1pm1(a, pol) - powm1(x, a, pol);
result /= a;
detail::small_gamma2_series<T> s(a, x);
boost::uintmax_t max_iter = BOOST_MATH_MAX_ITER;
#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
T zero = 0;
- result -= pow(x, a) * tools::sum_series(s, boost::math::tools::digits<T>(), max_iter, zero);
+ result -= pow(x, a) * tools::sum_series(s, boost::math::policy::digits<T, Policy>(), max_iter, zero);
#else
- result -= pow(x, a) * tools::sum_series(s, boost::math::tools::digits<T>(), max_iter);
+ result -= pow(x, a) * tools::sum_series(s, boost::math::policy::digits<T, Policy>(), max_iter);
#endif
- tools::check_series_iterations(BOOST_CURRENT_FUNCTION, max_iter);
+ policy::check_series_iterations("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol);
return result;
}
//
@@ -1230,14 +735,14 @@
//
// Upper gamma fraction for half integer a:
//
-template <class T>
-T finite_half_gamma_q(T a, T x, T* p_derivative)
+template <class T, class Policy>
+T finite_half_gamma_q(T a, T x, T* p_derivative, const Policy& pol)
{
//
// Calculates normalised Q when a is a half-integer:
//
using namespace std;
- T e = boost::math::erfc(sqrt(x));
+ T e = boost::math::erfc(sqrt(x), pol);
if((e != 0) && (a > 1))
{
T term = exp(-x) / sqrt(constants::pi<T>() * x);
@@ -1267,9 +772,9 @@
//
// Main incomplete gamma entry point, handles all four incomplete gamma's:
//
-template <class T, class L>
+template <class T, class Policy>
T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
- const L& l, T* p_derivative)
+ const Policy& pol, T* p_derivative)
{
if(a <= 0)
tools::domain_error<T>(BOOST_CURRENT_FUNCTION, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a);
@@ -1278,6 +783,8 @@
using namespace std;
+ typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
+
T result;
BOOST_ASSERT((p_derivative == 0) || (normalised == true));
@@ -1292,20 +799,20 @@
invert = !invert;
result = finite_gamma_q(a, x);
if(normalised == false)
- result *= gamma_imp(a, l);
+ result *= boost::math::tgamma(a, pol);
// TODO: calculate derivative inside sum:
if(p_derivative)
- *p_derivative = regularised_gamma_prefix(a, x, l);
+ *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
}
else if(is_half_int && is_small_a && (x > 0.2))
{
// calculate Q via finite sum for half integer a:
invert = !invert;
- result = finite_half_gamma_q(a, x, p_derivative);
+ result = finite_half_gamma_q(a, x, p_derivative, pol);
if(normalised == false)
- result *= gamma_imp(a, l);
+ result *= boost::math::tgamma(a, pol);
if(p_derivative && (*p_derivative == 0))
- *p_derivative = regularised_gamma_prefix(a, x, l);
+ *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
}
else if(x < 0.5)
{
@@ -1315,21 +822,21 @@
if(-0.4 / log(x) < a)
{
// Compute P:
- result = normalised ? regularised_gamma_prefix(a, x, l) : full_igamma_prefix(a, x);
+ result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
if(p_derivative)
*p_derivative = result;
if(result != 0)
- result *= detail::lower_gamma_series(a, x, boost::math::tools::digits<T>()) / a;
+ result *= detail::lower_gamma_series(a, x, pol) / a;
}
else
{
// Compute Q:
invert = !invert;
- result = tgamma_small_upper_part(a, x, l);
+ result = tgamma_small_upper_part(a, x, pol);
if(normalised)
- result /= gamma_imp(a, l);
+ result /= boost::math::tgamma(a, pol);
if(p_derivative)
- *p_derivative = regularised_gamma_prefix(a, x, l);
+ *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
}
}
else if(x < 1.1)
@@ -1340,21 +847,21 @@
if(x * 1.1f < a)
{
// Compute P:
- result = normalised ? regularised_gamma_prefix(a, x, l) : full_igamma_prefix(a, x);
+ result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
if(p_derivative)
*p_derivative = result;
if(result != 0)
- result *= detail::lower_gamma_series(a, x, boost::math::tools::digits<T>()) / a;
+ result *= detail::lower_gamma_series(a, x, pol) / a;
}
else
{
// Compute Q:
invert = !invert;
- result = tgamma_small_upper_part(a, x, l);
+ result = tgamma_small_upper_part(a, x, pol);
if(normalised)
- result /= gamma_imp(a, l);
+ result /= boost::math::tgamma(a, pol);
if(p_derivative)
- *p_derivative = regularised_gamma_prefix(a, x, l);
+ *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
}
}
else
@@ -1400,27 +907,28 @@
// these expansions, in that case we'll be calling
// an empty function.
//
- 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 <= 53),
+ typename mpl::if_<
+ mpl::less_equal<precision_type, mpl::int_<53> >,
mpl::int_<53>,
- typename mpl::if_c<
- (std::numeric_limits<T>::digits <= 64),
+ typename mpl::if_<
+ mpl::less_equal<precision_type, mpl::int_<64> >,
mpl::int_<64>,
mpl::int_<113>
>::type
>::type
>::type tag_type;
- result = igamma_temme_large(a, x, static_cast<tag_type const*>(0));
+ result = igamma_temme_large(a, x, pol, static_cast<tag_type const*>(0));
if(x >= a)
invert = !invert;
if(p_derivative)
- *p_derivative = regularised_gamma_prefix(a, x, l);
+ *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
}
else
{
@@ -1429,28 +937,28 @@
//
// Changeover here occurs at P ~ Q ~ 0.5
//
- result = normalised ? regularised_gamma_prefix(a, x, l) : full_igamma_prefix(a, x);
+ result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
if(p_derivative)
*p_derivative = result;
if(x < a)
{
// Compute P:
if(result != 0)
- result *= detail::lower_gamma_series(a, x, boost::math::tools::digits<T>()) / a;
+ result *= detail::lower_gamma_series(a, x, pol) / a;
}
else
{
// Compute Q:
invert = !invert;
if(result != 0)
- result *= upper_gamma_fraction(a, x, tools::digits<T>());
+ result *= upper_gamma_fraction(a, x, policy::digits<T, Policy>());
}
}
}
if(invert)
{
- T gam = normalised ? 1 : gamma_imp(a, l);
+ T gam = normalised ? 1 : boost::math::tgamma(a, pol);
result = gam - result;
}
if(p_derivative)
@@ -1470,24 +978,18 @@
return result;
}
-template <class T, class L>
-inline T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, const L& l)
-{
- return gamma_incomplete_imp(a, x, normalised, invert, l, static_cast<T*>(0));
-}
-
//
// Ratios of two gamma functions:
//
-template <class T, class L>
-T tgamma_delta_ratio_imp_lanczos(T z, T delta, const L&)
+template <class T, class Policy, class L>
+T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const L&)
{
using namespace std;
T zgh = z + L::g() - constants::half<T>();
T result;
if(fabs(delta) < 10)
{
- result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh));
+ result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));
}
else
{
@@ -1500,8 +1002,8 @@
//
// And again without Lanczos support this time:
//
-template <class T>
-T tgamma_delta_ratio_imp_lanczos(T z, T delta, const lanczos::undefined_lanczos&)
+template <class T, class Policy>
+T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const lanczos::undefined_lanczos&)
{
using namespace std;
//
@@ -1520,32 +1022,32 @@
}
if(delta < 10)
{
- prefix *= exp(-z * boost::math::log1p(delta / z));
+ prefix *= exp(-z * boost::math::log1p(delta / z, pol));
}
else
{
prefix *= pow(z / zd, z);
}
prefix *= pow(constants::e<T>() / zd, delta);
- T sum = detail::lower_gamma_series(z, z, ::boost::math::tools::digits<T>()) / z;
- sum += detail::upper_gamma_fraction(z, z, ::boost::math::tools::digits<T>());
- T sumd = detail::lower_gamma_series(zd, zd, ::boost::math::tools::digits<T>()) / zd;
- sumd += detail::upper_gamma_fraction(zd, zd, ::boost::math::tools::digits<T>());
+ T sum = detail::lower_gamma_series(z, z, pol) / z;
+ sum += detail::upper_gamma_fraction(z, z, ::boost::math::policy::digits<T, Policy>());
+ T sumd = detail::lower_gamma_series(zd, zd, pol) / zd;
+ sumd += detail::upper_gamma_fraction(zd, zd, ::boost::math::policy::digits<T, Policy>());
sum /= sumd;
if(fabs(tools::max_value<T>() / prefix) < fabs(sum))
- return tools::overflow_error<T>(BOOST_CURRENT_FUNCTION, "Result of tgamma is too large to represent.");
+ return policy::raise_overflow_error<T>("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Result of tgamma is too large to represent.", pol);
return sum * prefix;
}
-template <class T, class L>
-T tgamma_delta_ratio_imp(T z, T delta, const L& l)
+template <class T, class Policy>
+T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol)
{
using namespace std;
if(z <= 0)
- tools::domain_error<T>(BOOST_CURRENT_FUNCTION, "Gamma function ratios only implemented for positive arguments (got a=%1%).", z);
+ policy::raise_domain_error<T>("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got a=%1%).", z, pol);
if(z+delta <= 0)
- tools::domain_error<T>(BOOST_CURRENT_FUNCTION, "Gamma function ratios only implemented for positive arguments (got b=%1%).", z+delta);
+ policy::raise_domain_error<T>("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got b=%1%).", z+delta, pol);
if(floor(delta) == delta)
{
@@ -1590,35 +1092,37 @@
}
}
}
- return tgamma_delta_ratio_imp_lanczos(z, delta, l);
+ typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
+ return tgamma_delta_ratio_imp_lanczos(z, delta, pol, lanczos_type());
}
-template <class T, class L>
-T gamma_p_derivative_imp(T a, T x, const L& l)
+template <class T, class Policy>
+T gamma_p_derivative_imp(T a, T x, const Policy& pol)
{
//
// Usual error checks first:
//
if(a <= 0)
- tools::domain_error<T>(BOOST_CURRENT_FUNCTION, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a);
+ policy::raise_domain_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
if(x < 0)
- tools::domain_error<T>(BOOST_CURRENT_FUNCTION, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x);
+ policy::raise_domain_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
//
// Now special cases:
//
if(x == 0)
{
return (a > 1) ? 0 :
- (a == 1) ? 1 : tools::overflow_error<T>(BOOST_CURRENT_FUNCTION);
+ (a == 1) ? 1 : policy::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
}
//
// Normal case:
//
- T f1 = detail::regularised_gamma_prefix(a, x, l);
+ typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
+ T f1 = detail::regularised_gamma_prefix(a, x, pol, lanczos_type());
if((x < 1) && (tools::max_value<T>() * x < f1))
{
// overflow:
- return tools::overflow_error<T>(BOOST_CURRENT_FUNCTION);
+ return policy::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
}
f1 /= x;
@@ -1626,55 +1130,96 @@
return f1;
}
+template <class T, class Policy>
+inline typename tools::promote_args<T>::type
+ tgamma(T z, const Policy& pol, const mpl::true_)
+{
+ BOOST_FPU_EXCEPTION_GUARD
+ typedef typename tools::promote_args<T>::type result_type;
+ typedef typename policy::evaluation<result_type, Policy>::type value_type;
+ typedef typename lanczos::lanczos<result_type, Policy>::type evaluation_type;
+ return policy::checked_narrowing_cast<result_type, Policy>(detail::gamma_imp(static_cast<value_type>(z), pol, evaluation_type()), "boost::math::tgamma<%1%>(%1%)");
+}
+
+template <class T1, class T2, class Policy>
+inline typename tools::promote_args<T1, T2>::type
+ tgamma(T1 a, T2 z, const Policy& pol, const mpl::false_)
+{
+ BOOST_FPU_EXCEPTION_GUARD
+ typedef typename tools::promote_args<T1, T2>::type result_type;
+ typedef typename policy::evaluation<result_type, Policy>::type value_type;
+ typedef typename lanczos::lanczos<result_type, Policy>::type evaluation_type;
+ return policy::checked_narrowing_cast<result_type, Policy>(
+ detail::gamma_incomplete_imp(static_cast<value_type>(a),
+ static_cast<value_type>(z), false, true,
+ pol, static_cast<value_type*>(0)), "boost::math::tgamma<%1%>(%1%, %1%)");
+}
+
+template <class T1, class T2>
+inline typename tools::promote_args<T1, T2>::type
+ tgamma(T1 a, T2 z, const mpl::false_ tag)
+{
+ return tgamma(a, z, policy::policy<>(), tag);
+}
+
} // namespace detail
template <class T>
inline typename tools::promote_args<T>::type
tgamma(T z)
{
+ return tgamma(z, policy::policy<>());
+}
+
+template <class T, class Policy>
+inline typename tools::promote_args<T>::type
+ lgamma(T z, int* sign, const Policy& pol)
+{
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args<T>::type result_type;
- typedef typename lanczos::lanczos_traits<result_type>::value_type value_type;
- typedef typename lanczos::lanczos_traits<result_type>::evaluation_type evaluation_type;
- return tools::checked_narrowing_cast<result_type>(detail::gamma_imp(static_cast<value_type>(z), evaluation_type()), BOOST_CURRENT_FUNCTION);
+ typedef typename policy::evaluation<result_type, Policy>::type value_type;
+ typedef typename lanczos::lanczos<result_type, Policy>::type evaluation_type;
+ return policy::checked_narrowing_cast<result_type, Policy>(detail::lgamma_imp(static_cast<value_type>(z), pol, evaluation_type(), sign), "boost::math::lgamma<%1%>(%1%)");
}
template <class T>
inline typename tools::promote_args<T>::type
lgamma(T z, int* sign)
{
- BOOST_FPU_EXCEPTION_GUARD
- typedef typename tools::promote_args<T>::type result_type;
- typedef typename lanczos::lanczos_traits<result_type>::value_type value_type;
- typedef typename lanczos::lanczos_traits<result_type>::evaluation_type evaluation_type;
- return tools::checked_narrowing_cast<result_type>(detail::lgamma_imp(static_cast<value_type>(z), evaluation_type(), sign), BOOST_CURRENT_FUNCTION);
+ return lgamma(z, sign, policy::policy<>());
+}
+
+template <class T, class Policy>
+inline typename tools::promote_args<T>::type
+ lgamma(T x, const Policy& pol)
+{
+ return ::boost::math::lgamma(x, 0, pol);
}
template <class T>
inline typename tools::promote_args<T>::type
lgamma(T x)
{
- BOOST_FPU_EXCEPTION_GUARD
- return ::boost::math::lgamma(x, 0);
+ return ::boost::math::lgamma(x, 0, policy::policy<>());
}
-template <class T>
+template <class T, class Policy>
inline typename tools::promote_args<T>::type
- tgamma1pm1(T z)
+ tgamma1pm1(T z, const Policy& pol)
{
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args<T>::type result_type;
- typedef typename lanczos::lanczos_traits<result_type>::value_type value_type;
- typedef typename lanczos::lanczos_traits<result_type>::evaluation_type evaluation_type;
- typedef typename mpl::if_c<
- (std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::digits <= 64)),
- mpl::int_<64>,
- typename mpl::if_c<
- (std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::digits <= 113)),
- mpl::int_<113>, mpl::int_<0> >::type
- >::type tag_type;
+ typedef typename policy::evaluation<result_type, Policy>::type value_type;
+ typedef typename lanczos::lanczos<result_type, Policy>::type evaluation_type;
- return tools::checked_narrowing_cast<typename remove_cv<T>::type>(detail::tgammap1m1_imp(static_cast<value_type>(z), tag_type(), evaluation_type()), BOOST_CURRENT_FUNCTION);
+ return policy::checked_narrowing_cast<typename remove_cv<result_type>::type, Policy>(detail::tgammap1m1_imp(static_cast<value_type>(z), pol, evaluation_type()), "boost::math::tgamma1pm1<%!%>(%1%)");
+}
+
+template <class T>
+inline typename tools::promote_args<T>::type
+ tgamma1pm1(T z)
+{
+ return tgamma1pm1(z, policy::policy<>());
}
//
@@ -1684,95 +1229,131 @@
inline typename tools::promote_args<T1, T2>::type
tgamma(T1 a, T2 z)
{
- BOOST_FPU_EXCEPTION_GUARD
- typedef typename tools::promote_args<T1, T2>::type result_type;
- typedef typename lanczos::lanczos_traits<result_type>::value_type value_type;
- typedef typename lanczos::lanczos_traits<result_type>::evaluation_type evaluation_type;
- return tools::checked_narrowing_cast<result_type>(
- detail::gamma_incomplete_imp(static_cast<value_type>(a),
- static_cast<value_type>(z), false, true,
- evaluation_type()), BOOST_CURRENT_FUNCTION);
+ //
+ // Type T2 could be a policy object, or a value, select the
+ // right overload based on T2:
+ //
+ typedef typename policy::is_policy<T2>::type maybe_policy;
+ return detail::tgamma(a, z, maybe_policy());
+}
+template <class T1, class T2, class Policy>
+inline typename tools::promote_args<T1, T2>::type
+ tgamma(T1 a, T2 z, const Policy& pol)
+{
+ return detail::tgamma(a, z, pol, mpl::false_());
}
//
// Full lower incomplete gamma:
//
-template <class T1, class T2>
+template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type
- tgamma_lower(T1 a, T2 z)
+ tgamma_lower(T1 a, T2 z, const Policy& pol)
{
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args<T1, T2>::type result_type;
- typedef typename lanczos::lanczos_traits<result_type>::value_type value_type;
- typedef typename lanczos::lanczos_traits<result_type>::evaluation_type evaluation_type;
- return tools::checked_narrowing_cast<result_type>(
+ typedef typename policy::evaluation<result_type, Policy>::type value_type;
+ typedef typename lanczos::lanczos<result_type, Policy>::type evaluation_type;
+ return policy::checked_narrowing_cast<result_type, Policy>(
detail::gamma_incomplete_imp(static_cast<value_type>(a),
static_cast<value_type>(z), false, false,
- evaluation_type()), BOOST_CURRENT_FUNCTION);
+ pol, static_cast<value_type*>(0)), "tgamma_lower<%1%>(%1%, %1%)");
+}
+template <class T1, class T2>
+inline typename tools::promote_args<T1, T2>::type
+ tgamma_lower(T1 a, T2 z)
+{
+ return tgamma_lower(a, z, policy::policy<>());
}
//
// Regularised upper incomplete gamma:
//
-template <class T1, class T2>
+template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type
- gamma_q(T1 a, T2 z)
+ gamma_q(T1 a, T2 z, const Policy& pol)
{
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args<T1, T2>::type result_type;
- typedef typename lanczos::lanczos_traits<result_type>::value_type value_type;
- typedef typename lanczos::lanczos_traits<result_type>::evaluation_type evaluation_type;
- return tools::checked_narrowing_cast<result_type>(
+ typedef typename policy::evaluation<result_type, Policy>::type value_type;
+ typedef typename lanczos::lanczos<result_type, Policy>::type evaluation_type;
+ return policy::checked_narrowing_cast<result_type, Policy>(
detail::gamma_incomplete_imp(static_cast<value_type>(a),
static_cast<value_type>(z), true, true,
- evaluation_type()), BOOST_CURRENT_FUNCTION);
+ pol, static_cast<value_type*>(0)), "gamma_q<%1%>(%1%, %1%)");
+}
+template <class T1, class T2>
+inline typename tools::promote_args<T1, T2>::type
+ gamma_q(T1 a, T2 z)
+{
+ return gamma_q(a, z, policy::policy<>());
}
//
// Regularised lower incomplete gamma:
//
-template <class T1, class T2>
+template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type
- gamma_p(T1 a, T2 z)
+ gamma_p(T1 a, T2 z, const Policy& pol)
{
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args<T1, T2>::type result_type;
- typedef typename lanczos::lanczos_traits<result_type>::value_type value_type;
- typedef typename lanczos::lanczos_traits<result_type>::evaluation_type evaluation_type;
- return tools::checked_narrowing_cast<result_type>(
+ typedef typename policy::evaluation<result_type, Policy>::type value_type;
+ typedef typename lanczos::lanczos<result_type, Policy>::type evaluation_type;
+ return policy::checked_narrowing_cast<result_type, Policy>(
detail::gamma_incomplete_imp(static_cast<value_type>(a),
static_cast<value_type>(z), true, false,
- evaluation_type()), BOOST_CURRENT_FUNCTION);
+ pol, static_cast<value_type*>(0)), "gamma_p<%1%>(%1%, %1%)");
+}
+template <class T1, class T2>
+inline typename tools::promote_args<T1, T2>::type
+ gamma_p(T1 a, T2 z)
+{
+ return gamma_p(a, z, policy::policy<>());
}
// ratios of gamma functions:
+template <class T1, class T2, class Policy>
+inline typename tools::promote_args<T1, T2>::type
+ tgamma_delta_ratio(T1 z, T2 delta, const Policy& pol)
+{
+ BOOST_FPU_EXCEPTION_GUARD
+ typedef typename tools::promote_args<T1, T2>::type result_type;
+ typedef typename policy::evaluation<result_type, Policy>::type value_type;
+ return policy::checked_narrowing_cast<result_type, Policy>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(z), static_cast<value_type>(delta), pol), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
+}
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type
tgamma_delta_ratio(T1 z, T2 delta)
{
- BOOST_FPU_EXCEPTION_GUARD
+ return tgamma_delta_ratio(z, delta, policy::policy<>());
+}
+template <class T1, class T2, class Policy>
+inline typename tools::promote_args<T1, T2>::type
+ tgamma_ratio(T1 a, T2 b, const Policy& pol)
+{
typedef typename tools::promote_args<T1, T2>::type result_type;
- typedef typename lanczos::lanczos_traits<result_type>::value_type value_type;
- typedef typename lanczos::lanczos_traits<result_type>::evaluation_type evaluation_type;
- return tools::checked_narrowing_cast<result_type>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(z), static_cast<value_type>(delta), evaluation_type()), BOOST_CURRENT_FUNCTION);
+ typedef typename policy::evaluation<result_type, Policy>::type value_type;
+ return policy::checked_narrowing_cast<result_type, Policy>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(b) - static_cast<value_type>(a), pol), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
}
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type
tgamma_ratio(T1 a, T2 b)
{
+ return tgamma_ratio(a, b, policy::policy<>());
+}
+
+template <class T1, class T2, class Policy>
+inline typename tools::promote_args<T1, T2>::type
+ gamma_p_derivative(T1 a, T2 x, const Policy& pol)
+{
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args<T1, T2>::type result_type;
- typedef typename lanczos::lanczos_traits<result_type>::value_type value_type;
- typedef typename lanczos::lanczos_traits<result_type>::evaluation_type evaluation_type;
- return tools::checked_narrowing_cast<result_type>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(b - a), evaluation_type()), BOOST_CURRENT_FUNCTION);
+ typedef typename policy::evaluation<result_type, Policy>::type value_type;
+ return policy::checked_narrowing_cast<result_type, Policy>(detail::gamma_p_derivative_imp(static_cast<value_type>(a), static_cast<value_type>(x), pol), "boost::math::gamma_p_derivative<%1%>(%1%, %1%)");
}
-
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type
gamma_p_derivative(T1 a, T2 x)
{
- BOOST_FPU_EXCEPTION_GUARD
- typedef typename tools::promote_args<T1, T2>::type result_type;
- typedef typename lanczos::lanczos_traits<result_type>::value_type value_type;
- typedef typename lanczos::lanczos_traits<result_type>::evaluation_type evaluation_type;
- return tools::checked_narrowing_cast<result_type>(detail::gamma_p_derivative_imp(static_cast<value_type>(a), static_cast<value_type>(x), evaluation_type()), BOOST_CURRENT_FUNCTION);
+ return gamma_p_derivative(a, x, policy::policy<>());
}
} // namespace math
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/lanczos.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/lanczos.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/lanczos.hpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -12,6 +12,8 @@
#include <boost/cstdint.hpp>
#include <boost/math/tools/evaluation_type.hpp>
#include <boost/math/tools/rational.hpp>
+#include <boost/math/policy/policy.hpp>
+#include <boost/mpl/less_equal.hpp>
#include <limits.h>
@@ -32,8 +34,12 @@
// Max experimental error (with arbitary precision arithmetic) 9.516e-12
// Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
//
-struct lanczos6
+struct lanczos6 : public mpl::int_<35>
{
+ //
+ // Produces slightly better than float precision when evaluated at
+ // double precision:
+ //
template <class T>
static T lanczos_sum(const T& z)
{
@@ -124,8 +130,12 @@
// Max experimental error (with arbitary precision arithmetic) 2.16676e-19
// Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
//
-struct lanczos11
+struct lanczos11 : public mpl::int_<60>
{
+ //
+ // Produces slightly better than double precision when evaluated at
+ // extended-double precision:
+ //
template <class T>
static T lanczos_sum(const T& z)
{
@@ -246,8 +256,12 @@
// Max experimental error (with arbitary precision arithmetic) 9.2213e-23
// Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
//
-struct lanczos13
+struct lanczos13 : public mpl::int_<72>
{
+ //
+ // Produces slightly better than extended-double precision when evaluated at
+ // higher precision:
+ //
template <class T>
static T lanczos_sum(const T& z)
{
@@ -380,8 +394,12 @@
// Max experimental error (with arbitary precision arithmetic) 2.9524e-38
// Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
//
-struct lanczos22
+struct lanczos22 : public mpl::int_<120>
{
+ //
+ // Produces slightly better than 128-bit long-double precision when
+ // evaluated at higher precision:
+ //
template <class T>
static T lanczos_sum(const T& z)
{
@@ -568,8 +586,11 @@
// Max experimental error (with arbitary precision arithmetic) 8.111667e-8
// Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
//
-struct lanczos6m24
+struct lanczos6m24 : public mpl::int_<24>
{
+ //
+ // Use for float precision, when evaluated as a float:
+ //
template <class T>
static T lanczos_sum(const T& z)
{
@@ -660,8 +681,11 @@
// Max experimental error (with arbitary precision arithmetic) 1.196214e-17
// Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
//
-struct lanczos13m53
+struct lanczos13m53 : public mpl::int_<53>
{
+ //
+ // Use for double precision, when evaluated as a double:
+ //
template <class T>
static T lanczos_sum(const T& z)
{
@@ -794,8 +818,11 @@
// Max experimental error (with arbitary precision arithmetic) 2.7699e-26
// Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
//
-struct lanczos17m64
+struct lanczos17m64 : public mpl::int_<64>
{
+ //
+ // Use for extended-double precision, when evaluated as an extended-double:
+ //
template <class T>
static T lanczos_sum(const T& z)
{
@@ -952,8 +979,11 @@
// Max experimental error (with arbitary precision arithmetic) 1.0541e-38
// Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
//
-struct lanczos24m113
+struct lanczos24m113 : public mpl::int_<113>
{
+ //
+ // Use for long-double precision, when evaluated as an long-double:
+ //
template <class T>
static T lanczos_sum(const T& z)
{
@@ -1151,9 +1181,7 @@
//
// placeholder for no lanczos info available:
//
-struct undefined_lanczos
-{
-};
+struct undefined_lanczos : public mpl::int_<INT_MAX - 1> { };
//
// lanczos_traits, specialise this to point tgamma etc to
@@ -1183,17 +1211,17 @@
>::type
>::type evaluation_type;
};
-
+#if 0
#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
-#define BOOST_MATH_FLT_DIGITS std::numeric_limits<float>::digits
-#define BOOST_MATH_DBL_DIGITS std::numeric_limits<double>::digits
-#define BOOST_MATH_LDBL_DIGITS std::numeric_limits<long double>::digits
+#define BOOST_MATH_FLT_DIGITS ::std::numeric_limits<float>::digits
+#define BOOST_MATH_DBL_DIGITS ::std::numeric_limits<double>::digits
+#define BOOST_MATH_LDBL_DIGITS ::std::numeric_limits<long double>::digits
#else
#define BOOST_MATH_FLT_DIGITS FLT_MANT_DIG
#define BOOST_MATH_DBL_DIGITS DBL_MANT_DIG
#define BOOST_MATH_LDBL_DIGITS LDBL_MANT_DIG
#endif
-
+#endif
template<>
struct lanczos_traits<float>
{
@@ -1256,6 +1284,33 @@
>::type evaluation_type;
};
+typedef mpl::list<
+ lanczos6m24,
+ lanczos6,
+ lanczos13m53,
+ lanczos13,
+ lanczos17m64,
+ lanczos24m113,
+ lanczos22,
+ undefined_lanczos> lanczos_list;
+
+template <class Real, class Policy>
+struct lanczos
+{
+ typedef typename mpl::if_<
+ typename mpl::less_equal<
+ typename policy::precision<Real, Policy>::type,
+ mpl::int_<0>
+ >::type,
+ mpl::int_<INT_MAX - 2>,
+ typename policy::precision<Real, Policy>::type
+ >::type target_precision;
+
+ typedef typename mpl::deref<typename mpl::find_if<
+ lanczos_list,
+ mpl::less_equal<target_precision, mpl::_1> >::type>::type type;
+};
+
} // namespace lanczos
} // namespace math
} // namespace boost
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/log1p.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/log1p.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/log1p.hpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -11,7 +11,7 @@
#include <boost/limits.hpp>
#include <boost/math/tools/config.hpp>
#include <boost/math/tools/series.hpp>
-#include <boost/math/tools/precision.hpp>
+#include <boost/math/policy/error_handling.hpp>
#include <boost/math/special_functions/math_fwd.hpp>
#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
@@ -63,22 +63,24 @@
// require up to std::numeric_limits<T>::digits+1 terms to be calculated.
// It would be much more efficient to use the equivalence:
// log(1+x) == (log(1+x) * x) / ((1-x) - 1)
-// Unfortunately optimizing compilers make such a mess of this, that it performs
-// no better than log(1+x): which is to say not very well at all.
+// Unfortunately many optimizing compilers make such a mess of this, that
+// it performs no better than log(1+x): which is to say not very well at all.
//
-template <class T>
-typename tools::promote_args<T>::type log1p(T x)
+template <class T, class Policy>
+typename tools::promote_args<T>::type log1p(T x, const Policy& pol)
{ // The function returns the natural logarithm of 1 + x.
// A domain error occurs if x < -1. TODO should there be a check?
typedef typename tools::promote_args<T>::type result_type;
using namespace std;
+ static const char* function = "boost::math::log1p<%1%>(%1%)";
+
if(x < -1)
- return tools::domain_error<T>(
- BOOST_CURRENT_FUNCTION, "log1p(x) requires x > -1, but got x = %1%.", x);
+ return policy::raise_domain_error<T>(
+ function, "log1p(x) requires x > -1, but got x = %1%.", x, pol);
if(x == -1)
- return -tools::overflow_error<T>(
- BOOST_CURRENT_FUNCTION);
+ return -policy::raise_overflow_error<T>(
+ function, 0, pol);
result_type a = abs(result_type(x));
if(a > result_type(0.5L))
@@ -90,12 +92,12 @@
detail::log1p_series<result_type> s(x);
boost::uintmax_t max_iter = BOOST_MATH_MAX_ITER;
#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
- result_type result = tools::sum_series(s, tools::digits<T>() + 2, max_iter);
+ result_type result = tools::sum_series(s, policy::digits<result_type, Policy>(), max_iter);
#else
result_type zero = 0;
- result_type result = tools::sum_series(s, tools::digits<T>() + 2, max_iter, zero);
+ result_type result = tools::sum_series(s, policy::digits<result_type, Policy>(), max_iter, zero);
#endif
- tools::check_series_iterations(BOOST_CURRENT_FUNCTION, max_iter);
+ policy::check_series_iterations(function, max_iter, pol);
return result;
}
@@ -126,50 +128,50 @@
# if (defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901)) \
|| defined(linux) || defined(__linux) || defined(__linux__) \
|| defined(__hpux)
-template <>
-inline float log1p<float>(float x)
+template <class Policy>
+inline float log1p(float x, const Policy& pol)
{
if(x < -1)
- return tools::domain_error<float>(
- BOOST_CURRENT_FUNCTION, "log1p(x) requires x > -1, but got x = %1%.", x);
+ return policy::raise_domain_error<float>(
+ "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol);
if(x == -1)
- return -tools::overflow_error<float>(
- BOOST_CURRENT_FUNCTION);
+ return -policy::raise_overflow_error<float>(
+ "log1p<%1%>(%1%)", 0, pol);
return ::log1pf(x);
}
-template <>
-inline long double log1p<long double>(long double x)
+template <class Policy>
+inline long double log1p(long double x, const Policy& pol)
{
if(x < -1)
- return tools::domain_error<long double>(
- BOOST_CURRENT_FUNCTION, "log1p(x) requires x > -1, but got x = %1%.", x);
+ return policy::raise_domain_error<long double>(
+ "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol);
if(x == -1)
- return -tools::overflow_error<long double>(
- BOOST_CURRENT_FUNCTION);
+ return -policy::raise_overflow_error<long double>(
+ "log1p<%1%>(%1%), 0, pol");
return ::log1pl(x);
}
#else
-template <>
-inline float log1p<float>(float x)
+template <class Policy>
+inline float log1p(float x, const Policy& pol)
{
if(x < -1)
- return tools::domain_error<float>(
- BOOST_CURRENT_FUNCTION, "log1p(x) requires x > -1, but got x = %1%.", x);
+ return policy::raise_domain_error<float>(
+ "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol);
if(x == -1)
- return -tools::overflow_error<float>(
- BOOST_CURRENT_FUNCTION);
+ return -policy::raise_overflow_error<float>(
+ "log1p<%1%>(%1%)", 0, pol);
return ::log1p(x);
}
#endif
-template <>
-inline double log1p<double>(double x)
+template <class Policy>
+inline double log1p(double x, const Policy& pol)
{
if(x < -1)
- return tools::domain_error<double>(
- BOOST_CURRENT_FUNCTION, "log1p(x) requires x > -1, but got x = %1%.", x);
+ return policy::raise_domain_error<double>(
+ "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol);
if(x == -1)
- return -tools::overflow_error<double>(
- BOOST_CURRENT_FUNCTION);
+ return -policy::raise_overflow_error<double>(
+ "log1p<%1%>(%1%)", 0, pol);
return ::log1p(x);
}
#elif defined(_MSC_VER) && (BOOST_MSVC >= 1400)
@@ -178,35 +180,35 @@
// that your compilers optimizer won't mess this code up!!
// Currently tested with VC8 and Intel 9.1.
//
-template <>
-inline double log1p<double>(double x)
+template <class Policy>
+inline double log1p(double x, const Policy& pol)
{
if(x < -1)
- return tools::domain_error<double>(
- BOOST_CURRENT_FUNCTION, "log1p(x) requires x > -1, but got x = %1%.", x);
+ return policy::raise_domain_error<double>(
+ "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol);
if(x == -1)
- return -tools::overflow_error<double>(
- BOOST_CURRENT_FUNCTION);
+ return -policy::raise_overflow_error<double>(
+ "log1p<%1%>(%1%)", 0, pol);
double u = 1+x;
if(u == 1.0)
return x;
else
return log(u)*(x/(u-1.0));
}
-template <>
-inline float log1p<float>(float x)
+template <class Policy>
+inline float log1p(float x, const Policy& pol)
{
- return static_cast<float>(boost::math::log1p<double>(x));
+ return static_cast<float>(boost::math::log1p(static_cast<double>(x), pol));
}
-template <>
-inline long double log1p<long double>(long double x)
+template <class Policy>
+inline long double log1p(long double x, const Policy& pol)
{
if(x < -1)
- return tools::domain_error<long double>(
- BOOST_CURRENT_FUNCTION, "log1p(x) requires x > -1, but got x = %1%.", x);
+ return policy::raise_domain_error<long double>(
+ "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol);
if(x == -1)
- return -tools::overflow_error<long double>(
- BOOST_CURRENT_FUNCTION);
+ return -policy::raise_overflow_error<long double>(
+ "log1p<%1%>(%1%)", 0, pol);
long double u = 1+x;
if(u == 1.0)
return x;
@@ -215,6 +217,55 @@
}
#endif
+template <class T>
+inline typename tools::promote_args<T>::type log1p(T x)
+{
+ return boost::math::log1p(x, policy::policy<>());
+}
+//
+// Compute log(1+x)-x:
+//
+template <class T, class Policy>
+inline typename tools::promote_args<T>::type
+ log1pmx(T x, const Policy& pol)
+{
+ typedef typename tools::promote_args<T>::type result_type;
+ using namespace std;
+ static const char* function = "boost::math::log1pmx<%1%>(%1%)";
+
+ if(x < -1)
+ return policy::raise_domain_error<T>(
+ function, "log1pmx(x) requires x > -1, but got x = %1%.", x, pol);
+ if(x == -1)
+ return -policy::raise_overflow_error<T>(
+ function, 0, pol);
+
+ result_type a = abs(result_type(x));
+ if(a > result_type(0.95L))
+ return log(1 + result_type(x)) - result_type(x);
+ // Note that without numeric_limits specialisation support,
+ // epsilon just returns zero, and our "optimisation" will always fail:
+ if(a < tools::epsilon<result_type>())
+ return -x * x / 2;
+ boost::math::detail::log1p_series<T> s(x);
+ s();
+ boost::uintmax_t max_iter = BOOST_MATH_MAX_ITER;
+#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
+ T zero = 0;
+ T result = boost::math::tools::sum_series(s, policy::digits<T, Policy>(), max_iter, zero);
+#else
+ T result = boost::math::tools::sum_series(s, policy::digits<T, Policy>(), max_iter);
+#endif
+ tools::check_series_iterations(function, max_iter);
+ return result;
+}
+
+template <class T>
+inline T log1pmx(T x)
+{
+ return log1pmx(x, policy::policy<>());
+}
+
} // namespace math
} // namespace boost
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/math_fwd.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/math_fwd.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/math_fwd.hpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -71,15 +71,21 @@
// erf & erfc error functions.
template <class RT> // Error function.
typename tools::promote_args<RT>::type erf(RT z);
+ template <class RT, class Policy> // Error function.
+ typename tools::promote_args<RT>::type erf(RT z, const Policy&);
template <class RT>// Error function complement.
typename tools::promote_args<RT>::type erfc(RT z);
+ template <class RT, class Policy>// Error function complement.
+ typename tools::promote_args<RT>::type erfc(RT z, const Policy&);
template <class RT>// Error function inverse.
typename tools::promote_args<RT>::type erf_inv(RT z);
+ template <class RT, class Policy>// Error function inverse.
+ typename tools::promote_args<RT>::type erf_inv(RT z, const Policy& pol);
- template <class RT>// Error function complement inverse.
- typename tools::promote_args<RT>::type erfc_inv(RT z);
+ template <class RT, class Policy>// Error function complement inverse.
+ typename tools::promote_args<RT>::type erfc_inv(RT z, const Policy& pol);
// Polynomials:
template <class T1, class T2, class T3>
@@ -245,6 +251,10 @@
template <class T>
typename tools::promote_args<T>::type log1p(T);
+ // log1pmx is log(x + 1) - x
+ template <class T>
+ typename tools::promote_args<T>::type log1pmx(T);
+
// Exp (x) minus 1 functions.
template <class T>
typename tools::promote_args<T>::type expm1(T);
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/powm1.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/powm1.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/powm1.hpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -13,8 +13,8 @@
namespace boost{ namespace math{ namespace detail{
-template <class T>
-inline T powm1_imp(const T a, const T z)
+template <class T, class Policy>
+inline T powm1_imp(const T a, const T z, const Policy& pol)
{
using namespace std;
@@ -22,7 +22,7 @@
{
T p = log(a) * z;
if(fabs(p) < 2)
- return boost::math::expm1(p);
+ return boost::math::expm1(p, pol);
// otherwise fall though:
}
return pow(a, z) - 1;
@@ -35,7 +35,15 @@
powm1(const T1 a, const T2 z)
{
typedef typename tools::promote_args<T1, T2>::type result_type;
- return detail::powm1_imp(static_cast<result_type>(a), static_cast<result_type>(z));
+ return detail::powm1_imp(static_cast<result_type>(a), static_cast<result_type>(z), policy::policy<>());
+}
+
+template <class T1, class T2, class Policy>
+inline typename tools::promote_args<T1, T2>::type
+ powm1(const T1 a, const T2 z, const Policy& pol)
+{
+ typedef typename tools::promote_args<T1, T2>::type result_type;
+ return detail::powm1_imp(static_cast<result_type>(a), static_cast<result_type>(z), pol);
}
} // namespace math
Modified: sandbox/math_toolkit/policy/boost/math/special_functions/sqrt1pm1.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/special_functions/sqrt1pm1.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/special_functions/sqrt1pm1.hpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -16,15 +16,21 @@
namespace boost{ namespace math{
-template <class T>
-inline typename tools::promote_args<T>::type sqrt1pm1(const T& val)
+template <class T, class Policy>
+inline typename tools::promote_args<T>::type sqrt1pm1(const T& val, const Policy& pol)
{
typedef typename tools::promote_args<T>::type result_type;
using namespace std;
if(fabs(result_type(val)) > 0.75)
return sqrt(1 + result_type(val)) - 1;
- return boost::math::expm1(boost::math::log1p(val) / 2);
+ return boost::math::expm1(boost::math::log1p(val, pol) / 2, pol);
+}
+
+template <class T>
+inline typename tools::promote_args<T>::type sqrt1pm1(const T& val)
+{
+ return sqrt1pm1(val, policy::policy<>());
}
} // namespace math
Modified: sandbox/math_toolkit/policy/boost/math/tools/promotion.hpp
==============================================================================
--- sandbox/math_toolkit/policy/boost/math/tools/promotion.hpp (original)
+++ sandbox/math_toolkit/policy/boost/math/tools/promotion.hpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -26,6 +26,8 @@
#include <boost/type_traits/remove_cv.hpp>// for boost::remove_cv
// Boost Template meta programming:
#include <boost/mpl/if.hpp> // for boost::mpl::if_c.
+#include <boost/mpl/and.hpp> // for boost::mpl::if_c.
+#include <boost/mpl/or.hpp> // for boost::mpl::if_c.
namespace boost
{
@@ -66,11 +68,11 @@
typedef typename promote_arg<T1>::type T1P; // T1 perhaps promoted.
typedef typename promote_arg<T2>::type T2P; // T2 perhaps promoted.
- typedef typename mpl::if_c<
- ::boost::is_floating_point<T1P>::value && ::boost::is_floating_point<T2P>::value, // both T1P and T2P are floating-point?
- typename mpl::if_c< ::boost::is_same<long double, T1P>::value || ::boost::is_same<long double, T2P>::value, // either long double?
+ typedef typename mpl::if_<
+ typename mpl::and_<is_floating_point<T1P>, is_floating_point<T2P> >::type, // both T1P and T2P are floating-point?
+ typename mpl::if_< typename mpl::or_<is_same<long double, T1P>, is_same<long double, T2P> >::type, // either long double?
long double, // then result type is long double.
- typename mpl::if_c< ::boost::is_same<double, T1P>::value || ::boost::is_same<double, T2P>::value, // either double?
+ typename mpl::if_< typename mpl::or_<is_same<double, T1P>, is_same<double, T2P> >::type, // either double?
double, // result type is double.
float // else result type is float.
>::type
Modified: sandbox/math_toolkit/policy/libs/math/test/test_igamma.cpp
==============================================================================
--- sandbox/math_toolkit/policy/libs/math/test/test_igamma.cpp (original)
+++ sandbox/math_toolkit/policy/libs/math/test/test_igamma.cpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -3,6 +3,8 @@
// 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/math/special_functions/gamma.hpp>
#include <boost/test/included/test_exec_monitor.hpp>
@@ -12,8 +14,10 @@
#include <boost/math/constants/constants.hpp>
#include <boost/type_traits/is_floating_point.hpp>
#include <boost/array.hpp>
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
#include <boost/lambda/lambda.hpp>
#include <boost/lambda/bind.hpp>
+#endif
#include "test_gamma_hooks.hpp"
#include "handle_test_result.hpp"
@@ -228,6 +232,7 @@
template <class T>
void do_test_gamma_2(const T& data, const char* type_name, const char* test_name)
{
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
typedef typename T::value_type row_type;
typedef typename row_type::value_type value_type;
@@ -308,6 +313,7 @@
}
#endif
std::cout << std::endl;
+#endif
}
template <class T>
@@ -390,16 +396,20 @@
test_spots(0.0);
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
test_spots(0.0L);
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
test_spots(boost::math::concepts::real_concept(0.1));
#endif
+#endif
test_gamma(0.1F, "float");
test_gamma(0.1, "double");
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
test_gamma(0.1L, "long double");
#ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
test_gamma(boost::math::concepts::real_concept(0.1), "real_concept");
#endif
+#endif
#else
std::cout << "<note>The long double tests have been disabled on this platform "
"either because the long double overloads of the usual math functions are "
Modified: sandbox/math_toolkit/policy/libs/math/test/test_igamma_inv.cpp
==============================================================================
--- sandbox/math_toolkit/policy/libs/math/test/test_igamma_inv.cpp (original)
+++ sandbox/math_toolkit/policy/libs/math/test/test_igamma_inv.cpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -12,8 +12,10 @@
#include <boost/math/constants/constants.hpp>
#include <boost/type_traits/is_floating_point.hpp>
#include <boost/array.hpp>
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
#include <boost/lambda/lambda.hpp>
#include <boost/lambda/bind.hpp>
+#endif
#include "test_gamma_hooks.hpp"
#include "handle_test_result.hpp"
@@ -254,6 +256,7 @@
template <class T>
void do_test_gamma_inv(const T& data, const char* type_name, const char* test_name)
{
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
typedef typename T::value_type row_type;
typedef typename row_type::value_type value_type;
@@ -284,6 +287,7 @@
bind(funcp, ret<value_type>(_1[0]), ret<value_type>(_1[1])),
ret<value_type>(_1[3]));
handle_test_result(result, data[result.worst()], result.worst(), type_name, "boost::math::gamma_q_inv", test_name);
+#endif
}
template <class T>
@@ -365,16 +369,20 @@
test_spots(0.0, "double");
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
test_spots(0.0L, "long double");
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
test_spots(boost::math::concepts::real_concept(0.1), "real_concept");
#endif
+#endif
test_gamma(0.1F, "float");
test_gamma(0.1, "double");
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
test_gamma(0.1L, "long double");
#ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
test_gamma(boost::math::concepts::real_concept(0.1), "real_concept");
#endif
+#endif
#else
std::cout << "<note>The long double tests have been disabled on this platform "
"either because the long double overloads of the usual math functions are "
Modified: sandbox/math_toolkit/policy/libs/math/test/test_igamma_inva.cpp
==============================================================================
--- sandbox/math_toolkit/policy/libs/math/test/test_igamma_inva.cpp (original)
+++ sandbox/math_toolkit/policy/libs/math/test/test_igamma_inva.cpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -12,8 +12,10 @@
#include <boost/math/constants/constants.hpp>
#include <boost/type_traits/is_floating_point.hpp>
#include <boost/array.hpp>
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
#include <boost/lambda/lambda.hpp>
#include <boost/lambda/bind.hpp>
+#endif
#include "handle_test_result.hpp"
@@ -180,6 +182,7 @@
template <class T>
void do_test_gamma_inva(const T& data, const char* type_name, const char* test_name)
{
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
typedef typename T::value_type row_type;
typedef typename row_type::value_type value_type;
@@ -210,6 +213,7 @@
bind(funcp, ret<value_type>(_1[0]), ret<value_type>(_1[1])),
ret<value_type>(_1[3]));
handle_test_result(result, data[result.worst()], result.worst(), type_name, "boost::math::gamma_q_inva", test_name);
+#endif
}
template <class T>
@@ -250,8 +254,10 @@
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
test_gamma(0.1L, "long double");
#ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
test_gamma(boost::math::concepts::real_concept(0.1), "real_concept");
#endif
+#endif
#else
std::cout << "<note>The long double tests have been disabled on this platform "
"either because the long double overloads of the usual math functions are "
Modified: sandbox/math_toolkit/policy/libs/math/test/test_tgamma_ratio.cpp
==============================================================================
--- sandbox/math_toolkit/policy/libs/math/test/test_tgamma_ratio.cpp (original)
+++ sandbox/math_toolkit/policy/libs/math/test/test_tgamma_ratio.cpp 2007-06-20 14:03:05 EDT (Wed, 20 Jun 2007)
@@ -2,6 +2,9 @@
// 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)
+
+#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>
@@ -9,8 +12,10 @@
#include <boost/math/tools/stats.hpp>
#include <boost/math/tools/test.hpp>
#include <boost/array.hpp>
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
#include <boost/lambda/lambda.hpp>
#include <boost/lambda/bind.hpp>
+#endif
#include "handle_test_result.hpp"
@@ -112,6 +117,7 @@
template <class T>
void do_test_tgamma_delta_ratio(const T& data, const char* type_name, const char* test_name)
{
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
typedef typename T::value_type row_type;
typedef typename row_type::value_type value_type;
@@ -140,6 +146,7 @@
-boost::lambda::ret<value_type>(boost::lambda::_1[1])),
boost::lambda::ret<value_type>(boost::lambda::_1[3]));
handle_test_result(result, data[result.worst()], result.worst(), type_name, "boost::math::tgamma_delta_ratio(a -delta)", test_name);
+#endif
}
template <class T>
@@ -201,8 +208,10 @@
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
test_tgamma_ratio(0.1L, "long double");
#ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
test_tgamma_ratio(boost::math::concepts::real_concept(0.1), "real_concept");
#endif
+#endif
#else
std::cout << "<note>The long double tests have been disabled on this platform "
"either because the long double overloads of the usual math functions are "
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