From: john_at_[hidden]
Date: 2008-01-05 05:10:03
Author: johnmaddock
Date: 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
New Revision: 42474
URL: http://svn.boost.org/trac/boost/changeset/42474
Log:
Added initial support for interval arithmetic and mpfr extended-precision support.
Added:
sandbox/interval_math_toolkit/boost/math/concepts/real_type_concept.hpp (contents, props changed)
sandbox/interval_math_toolkit/libs/math/test/interval_concept_check.cpp (contents, props changed)
sandbox/interval_math_toolkit/libs/math/test/mpfr_concept_check.cpp (contents, props changed)
sandbox/interval_math_toolkit/libs/math/test/ntl_concept_check.cpp (contents, props changed)
Text files modified:
sandbox/interval_math_toolkit/boost/math/bindings/mpfr.hpp | 1
sandbox/interval_math_toolkit/boost/math/concepts/distributions.hpp | 2
sandbox/interval_math_toolkit/boost/math/constants/constants.hpp | 16 +++
sandbox/interval_math_toolkit/boost/math/distributions/binomial.hpp | 6
sandbox/interval_math_toolkit/boost/math/distributions/cauchy.hpp | 2
sandbox/interval_math_toolkit/boost/math/distributions/chi_squared.hpp | 6
sandbox/interval_math_toolkit/boost/math/distributions/detail/inv_discrete_quantile.hpp | 9 +
sandbox/interval_math_toolkit/boost/math/distributions/exponential.hpp | 4
sandbox/interval_math_toolkit/boost/math/distributions/extreme_value.hpp | 4
sandbox/interval_math_toolkit/boost/math/distributions/fisher_f.hpp | 16 ++--
sandbox/interval_math_toolkit/boost/math/distributions/lognormal.hpp | 2
sandbox/interval_math_toolkit/boost/math/distributions/negative_binomial.hpp | 10 +-
sandbox/interval_math_toolkit/boost/math/distributions/normal.hpp | 6
sandbox/interval_math_toolkit/boost/math/distributions/poisson.hpp | 8 +-
sandbox/interval_math_toolkit/boost/math/distributions/rayleigh.hpp | 4
sandbox/interval_math_toolkit/boost/math/distributions/triangular.hpp | 6
sandbox/interval_math_toolkit/boost/math/distributions/weibull.hpp | 24 +++---
sandbox/interval_math_toolkit/boost/math/policies/error_handling.hpp | 20 ++++-
sandbox/interval_math_toolkit/boost/math/policies/policy.hpp | 13 +-
sandbox/interval_math_toolkit/boost/math/special_functions/bessel.hpp | 12 +-
sandbox/interval_math_toolkit/boost/math/special_functions/beta.hpp | 61 ++++++++-------
sandbox/interval_math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp | 10 +-
sandbox/interval_math_toolkit/boost/math/special_functions/detail/bessel_jy.hpp | 27 +++---
sandbox/interval_math_toolkit/boost/math/special_functions/detail/bessel_jy_asym.hpp | 12 ++-
sandbox/interval_math_toolkit/boost/math/special_functions/detail/gamma_inva.hpp | 10 +-
sandbox/interval_math_toolkit/boost/math/special_functions/detail/ibeta_inv_ab.hpp | 27 +++---
sandbox/interval_math_toolkit/boost/math/special_functions/detail/ibeta_inverse.hpp | 19 ++--
sandbox/interval_math_toolkit/boost/math/special_functions/detail/igamma_inverse.hpp | 11 +-
sandbox/interval_math_toolkit/boost/math/special_functions/detail/lgamma_small.hpp | 4
sandbox/interval_math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp | 93 +++++++++++++++--------
sandbox/interval_math_toolkit/boost/math/special_functions/detail/unchecked_factorial.hpp | 19 ++++
sandbox/interval_math_toolkit/boost/math/special_functions/ellint_1.hpp | 6
sandbox/interval_math_toolkit/boost/math/special_functions/ellint_2.hpp | 4
sandbox/interval_math_toolkit/boost/math/special_functions/ellint_3.hpp | 8 +-
sandbox/interval_math_toolkit/boost/math/special_functions/ellint_rj.hpp | 2
sandbox/interval_math_toolkit/boost/math/special_functions/erf.hpp | 6
sandbox/interval_math_toolkit/boost/math/special_functions/expint.hpp | 4
sandbox/interval_math_toolkit/boost/math/special_functions/expm1.hpp | 2
sandbox/interval_math_toolkit/boost/math/special_functions/factorials.hpp | 14 +-
sandbox/interval_math_toolkit/boost/math/special_functions/gamma.hpp | 60 +++++++++------
sandbox/interval_math_toolkit/boost/math/special_functions/legendre.hpp | 2
sandbox/interval_math_toolkit/boost/math/special_functions/log1p.hpp | 4
sandbox/interval_math_toolkit/boost/math/special_functions/sin_pi.hpp | 2
sandbox/interval_math_toolkit/boost/math/special_functions/spherical_harmonic.hpp | 4
sandbox/interval_math_toolkit/boost/math/special_functions/sqrt1pm1.hpp | 2
sandbox/interval_math_toolkit/boost/math/special_functions/zeta.hpp | 4
sandbox/interval_math_toolkit/boost/math/tools/config.hpp | 153 +++++++++++++++++++++++++++++++++++++++
sandbox/interval_math_toolkit/boost/math/tools/fraction.hpp | 8 +-
sandbox/interval_math_toolkit/boost/math/tools/minima.hpp | 4
sandbox/interval_math_toolkit/boost/math/tools/rational.hpp | 12 +-
sandbox/interval_math_toolkit/boost/math/tools/real_cast.hpp | 7 +
sandbox/interval_math_toolkit/boost/math/tools/roots.hpp | 151 +++++++++++++++++++-------------------
sandbox/interval_math_toolkit/boost/math/tools/series.hpp | 12 +-
sandbox/interval_math_toolkit/boost/math/tools/test.hpp | 115 +++++++++++++++++++++++++++++
sandbox/interval_math_toolkit/boost/math/tools/toms748_solve.hpp | 16 ++-
sandbox/interval_math_toolkit/libs/math/test/handle_test_result.hpp | 47 ++++++++---
sandbox/interval_math_toolkit/libs/math/test/log1p_expm1_test.cpp | 11 ++
sandbox/interval_math_toolkit/libs/math/test/powm1_sqrtp1m1_test.cpp | 9 ++
sandbox/interval_math_toolkit/libs/math/test/test_beta.cpp | 9 ++
sandbox/interval_math_toolkit/libs/math/test/test_binomial_coeff.cpp | 9 ++
sandbox/interval_math_toolkit/libs/math/test/test_carlson.cpp | 16 +++
sandbox/interval_math_toolkit/libs/math/test/test_cbrt.cpp | 9 ++
sandbox/interval_math_toolkit/libs/math/test/test_digamma.cpp | 9 ++
sandbox/interval_math_toolkit/libs/math/test/test_ellint_1.cpp | 11 ++
sandbox/interval_math_toolkit/libs/math/test/test_ellint_2.cpp | 10 ++
sandbox/interval_math_toolkit/libs/math/test/test_ellint_3.cpp | 10 ++
sandbox/interval_math_toolkit/libs/math/test/test_erf.cpp | 9 ++
sandbox/interval_math_toolkit/libs/math/test/test_expint.cpp | 9 ++
sandbox/interval_math_toolkit/libs/math/test/test_gamma.cpp | 9 ++
sandbox/interval_math_toolkit/libs/math/test/test_hermite.cpp | 11 ++
sandbox/interval_math_toolkit/libs/math/test/test_ibeta.cpp | 9 ++
sandbox/interval_math_toolkit/libs/math/test/test_ibeta_inv.cpp | 12 +++
72 files changed, 904 insertions(+), 371 deletions(-)
Modified: sandbox/interval_math_toolkit/boost/math/bindings/mpfr.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/bindings/mpfr.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/bindings/mpfr.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -705,7 +705,6 @@
return value;
}
-
} // namespace detail
}}
Modified: sandbox/interval_math_toolkit/boost/math/concepts/distributions.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/concepts/distributions.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/concepts/distributions.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -162,6 +162,7 @@
v = quantile(complement(dist, d));
v = hazard(dist, d);
v = chf(dist, d);
+#ifndef TEST_MPFR
long double ld = 1;
v = cdf(dist, ld);
v = cdf(complement(dist, ld));
@@ -170,6 +171,7 @@
v = quantile(complement(dist, ld));
v = hazard(dist, ld);
v = chf(dist, ld);
+#endif
int i = 1;
v = cdf(dist, i);
v = cdf(complement(dist, i));
Added: sandbox/interval_math_toolkit/boost/math/concepts/real_type_concept.hpp
==============================================================================
--- (empty file)
+++ sandbox/interval_math_toolkit/boost/math/concepts/real_type_concept.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -0,0 +1,111 @@
+// Copyright John Maddock 2007-8.
+// 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_REAL_TYPE_CONCEPT_HPP
+#define BOOST_MATH_REAL_TYPE_CONCEPT_HPP
+
+#ifdef BOOST_MSVC
+#pragma warning(push)
+#pragma warning(disable: 4100)
+#pragma warning(disable: 4510)
+#pragma warning(disable: 4610)
+#endif
+#include <boost/concept_check.hpp>
+#ifdef BOOST_MSVC
+#pragma warning(pop)
+#endif
+#include <boost/math/tools/config.hpp>
+#include <boost/math/tools/precision.hpp>
+
+
+namespace boost{ namespace math{ namespace concepts{
+
+template <class RealType>
+struct RealTypeConcept
+{
+ template <class Other>
+ void check_binary_ops(Other o)
+ {
+ using namespace boost::numeric::interval_math_compare;
+ RealType r(o);
+ r = o;
+ r -= o;
+ r += o;
+ r *= o;
+ r /= o;
+ r = r - o;
+ r = o - r;
+ r = r + o;
+ r = o + r;
+ r = o * r;
+ r = r * o;
+ r = r / o;
+ r = o / r;
+ bool b;
+ b = r == o;
+ b = o == r;
+ b = r != o;
+ b = o != r;
+ b = r <= o;
+ b = o <= r;
+ b = r >= o;
+ b = o >= r;
+ b = r < o;
+ b = o < r;
+ b = r > o;
+ b = o > r;
+ }
+
+ void constraints()
+ {
+ BOOST_MATH_STD_USING
+ using namespace boost::numeric::interval_math_compare;
+
+ RealType r;
+ check_binary_ops(r);
+ check_binary_ops(0.5f);
+ check_binary_ops(0.5);
+ //check_binary_ops(0.5L);
+ check_binary_ops(1);
+ //check_binary_ops(1u);
+ check_binary_ops(1L);
+ //check_binary_ops(1uL);
+#ifndef BOOST_HAS_LONG_LONG
+ check_binary_ops(1LL);
+#endif
+ RealType r2 = +r;
+ r2 = -r;
+
+ r2 = fabs(r);
+ r2 = abs(r);
+ r2 = ceil(r);
+ r2 = floor(r);
+ r2 = exp(r);
+ r2 = pow(r, r2);
+ r2 = sqrt(r);
+ r2 = log(r);
+ r2 = cos(r);
+ r2 = sin(r);
+ r2 = tan(r);
+ r2 = asin(r);
+ r2 = acos(r);
+ r2 = atan(r);
+ int i;
+ r2 = ldexp(r, i);
+ r2 = frexp(r, &i);
+ i = boost::math::tools::digits<RealType>();
+ r2 = boost::math::tools::max_value<RealType>();
+ r2 = boost::math::tools::min_value<RealType>();
+ r2 = boost::math::tools::log_max_value<RealType>();
+ r2 = boost::math::tools::log_min_value<RealType>();
+ r2 = boost::math::tools::epsilon<RealType>();
+ }
+}; // struct DistributionConcept
+
+
+}}} // namespaces
+
+#endif
+
Modified: sandbox/interval_math_toolkit/boost/math/constants/constants.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/constants/constants.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/constants/constants.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -35,12 +35,24 @@
// (This is necessary because you can't use a numeric constant
// since even a long double might not have enough digits).
+// TODO: this does not create propper intervals!
#define BOOST_DEFINE_MATH_CONSTANT(name, x, y, exp)\
+ template <class T> T name(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T));\
+ namespace detail{\
+ template <class T> inline T name(boost::mpl::false_& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))\
+ {\
+ static const T result = ::boost::lexical_cast<T>(BOOST_STRINGIZE(BOOST_JOIN(BOOST_JOIN(x, y), BOOST_JOIN(e, exp))));\
+ return result;\
+ }\
+ template <class T> inline T name(boost::mpl::true_& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))\
+ {\
+ return boost::math::constants::name<typename T::base_type>();\
+ }\
+ }\
template <class T> inline T name(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T))\
{\
- static const T result = ::boost::lexical_cast<T>(BOOST_STRINGIZE(BOOST_JOIN(BOOST_JOIN(x, y), BOOST_JOIN(e, exp))));\
- return result;\
+ return detail::name<T>(boost::numeric::is_interval<T>());\
}\
template <> inline float name<float>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(float))\
{ return BOOST_JOIN(BOOST_JOIN(x, BOOST_JOIN(e, exp)), F); }\
Modified: sandbox/interval_math_toolkit/boost/math/distributions/binomial.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/distributions/binomial.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/distributions/binomial.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -173,7 +173,7 @@
// kurtosis:
// T k = (1 - 6 * sf * (1 - sf) ) / (n * sf * (1 - sf));
// Get the inverse of a std normal distribution:
- T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>();
+ T x = boost::math::erfc_inv(p > q ? static_cast<T>(2 * q) : static_cast<T>(2 * p), pol) * constants::root_two<T>();
// Set the sign:
if(p < 0.5)
x = -x;
@@ -649,13 +649,13 @@
template <class RealType, class Policy>
inline RealType quantile(const binomial_distribution<RealType, Policy>& dist, const RealType& p)
{
- return binomial_detail::quantile_imp(dist, p, 1-p);
+ return binomial_detail::quantile_imp(dist, p, static_cast<RealType>(1-p));
} // quantile
template <class RealType, class Policy>
RealType quantile(const complemented2_type<binomial_distribution<RealType, Policy>, RealType>& c)
{
- return binomial_detail::quantile_imp(c.dist, 1-c.param, c.param);
+ return binomial_detail::quantile_imp(c.dist, static_cast<RealType>(1-c.param), c.param);
} // quantile
template <class RealType, class Policy>
Modified: sandbox/interval_math_toolkit/boost/math/distributions/cauchy.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/distributions/cauchy.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/distributions/cauchy.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -140,7 +140,7 @@
return location;
}
result = -scale / tan(constants::pi<RealType>() * P);
- return complement ? location - result : location + result;
+ return complement ? static_cast<RealType>(location - result) : static_cast<RealType>(location + result);
} // quantile
} // namespace detail
Modified: sandbox/interval_math_toolkit/boost/math/distributions/chi_squared.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/distributions/chi_squared.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/distributions/chi_squared.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -107,7 +107,7 @@
}
}
- return gamma_p_derivative(degrees_of_freedom / 2, chi_square / 2, Policy()) / 2;
+ return gamma_p_derivative(static_cast<RealType>(degrees_of_freedom / 2), static_cast<RealType>(chi_square / 2), Policy()) / 2;
} // pdf
template <class RealType, class Policy>
@@ -128,7 +128,7 @@
function, "Chi Square parameter was %1%, but must be > 0 !", chi_square, Policy());
}
- return boost::math::gamma_p(degrees_of_freedom / 2, chi_square / 2, Policy());
+ return boost::math::gamma_p(static_cast<RealType>(degrees_of_freedom / 2), static_cast<RealType>(chi_square / 2), Policy());
} // cdf
template <class RealType, class Policy>
@@ -165,7 +165,7 @@
function, "Chi Square parameter was %1%, but must be > 0 !", chi_square, Policy());
}
- return boost::math::gamma_q(degrees_of_freedom / 2, chi_square / 2, Policy());
+ return boost::math::gamma_q(static_cast<RealType>(degrees_of_freedom / 2), static_cast<RealType>(chi_square / 2), Policy());
}
template <class RealType, class Policy>
Modified: sandbox/interval_math_toolkit/boost/math/distributions/detail/inv_discrete_quantile.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/distributions/detail/inv_discrete_quantile.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/distributions/detail/inv_discrete_quantile.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -127,8 +127,9 @@
}
else
{
+ using std::max;
b = a;
- a = (std::max)(b - 1, value_type(0));
+ a = (max)(static_cast<value_type>(b - 1), value_type(0));
if(a < min_bound)
a = min_bound;
fa = f(a);
@@ -158,7 +159,8 @@
}
else
{
- b = (std::max)(a - adder, value_type(0));
+ using std::max;
+ b = max BOOST_PREVENT_MACRO_SUBSTITUTION(static_cast<value_type>(a - adder), value_type(0));
if(b < min_bound)
b = min_bound;
}
@@ -182,7 +184,8 @@
}
else
{
- b = (std::max)(a - adder, value_type(0));
+ using std::max;
+ b = max BOOST_PREVENT_MACRO_SUBSTITUTION(static_cast<value_type>(a - adder), value_type(0));
if(b < min_bound)
b = min_bound;
}
Modified: sandbox/interval_math_toolkit/boost/math/distributions/exponential.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/distributions/exponential.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/distributions/exponential.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -124,7 +124,7 @@
return result;
if(0 == detail::verify_exp_x(function, x, &result, Policy()))
return result;
- result = -boost::math::expm1(-x * lambda, Policy());
+ result = -boost::math::expm1(static_cast<RealType>(-x * lambda), Policy());
return result;
} // cdf
@@ -148,7 +148,7 @@
if(p == 1)
return policies::raise_overflow_error<RealType>(function, 0, Policy());
- result = -boost::math::log1p(-p, Policy()) / lambda;
+ result = -boost::math::log1p(static_cast<RealType>(-p), Policy()) / lambda;
return result;
} // quantile
Modified: sandbox/interval_math_toolkit/boost/math/distributions/extreme_value.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/distributions/extreme_value.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/distributions/extreme_value.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -153,7 +153,7 @@
if(0 == detail::verify_scale_b("boost::math::cdf(const extreme_value_distribution<%1%>&, %1%)", b, &result, Policy()))
return result;
- result = -boost::math::expm1(-exp((a-c.param)/b), Policy());
+ result = -boost::math::expm1(static_cast<RealType>(-exp((a-c.param)/b)), Policy());
return result;
}
@@ -179,7 +179,7 @@
if(q == 1)
return -policies::raise_overflow_error<RealType>(function, 0, Policy());
- result = a - log(-boost::math::log1p(-q, Policy())) * b;
+ result = a - log(-boost::math::log1p(static_cast<RealType>(-q), Policy())) * b;
return result;
}
Modified: sandbox/interval_math_toolkit/boost/math/distributions/fisher_f.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/distributions/fisher_f.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/distributions/fisher_f.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -115,13 +115,13 @@
if(v1x > df2)
{
result = (df2 * df1) / ((df2 + v1x) * (df2 + v1x));
- result *= ibeta_derivative(df2 / 2, df1 / 2, df2 / (df2 + v1x), Policy());
+ result *= ibeta_derivative(static_cast<RealType>(df2 / 2), static_cast<RealType>(df1 / 2), static_cast<RealType>(df2 / (df2 + v1x)), Policy());
}
else
{
result = df2 + df1 * x;
result = (result * df1 - x * df1 * df1) / (result * result);
- result *= ibeta_derivative(df1 / 2, df2 / 2, v1x / (df2 + v1x), Policy());
+ result *= ibeta_derivative(static_cast<RealType>(df1 / 2), static_cast<RealType>(df2 / 2), static_cast<RealType>(v1x / (df2 + v1x)), Policy());
}
return result;
} // pdf
@@ -157,8 +157,8 @@
// to switch things around so we're passing 1-z instead.
//
return v1x > df2
- ? boost::math::ibetac(df2 / 2, df1 / 2, df2 / (df2 + v1x), Policy())
- : boost::math::ibeta(df1 / 2, df2 / 2, v1x / (df2 + v1x), Policy());
+ ? boost::math::ibetac(static_cast<RealType>(df2 / 2), static_cast<RealType>(df1 / 2), static_cast<RealType>(df2 / (df2 + v1x)), Policy())
+ : boost::math::ibeta(static_cast<RealType>(df1 / 2), static_cast<RealType>(df2 / 2), static_cast<RealType>(v1x / (df2 + v1x)), Policy());
} // cdf
template <class RealType, class Policy>
@@ -179,7 +179,7 @@
RealType x, y;
- x = boost::math::ibeta_inv(df1 / 2, df2 / 2, p, &y, Policy());
+ x = boost::math::ibeta_inv(static_cast<RealType>(df1 / 2), static_cast<RealType>(df2 / 2), p, &y, Policy());
return df2 * x / (df1 * y);
} // quantile
@@ -216,8 +216,8 @@
// to switch things around so we're passing 1-z instead.
//
return v1x > df2
- ? boost::math::ibeta(df2 / 2, df1 / 2, df2 / (df2 + v1x), Policy())
- : boost::math::ibetac(df1 / 2, df2 / 2, v1x / (df2 + v1x), Policy());
+ ? boost::math::ibeta(static_cast<RealType>(df2 / 2), static_cast<RealType>(df1 / 2), static_cast<RealType>(df2 / (df2 + v1x)), Policy())
+ : boost::math::ibetac(static_cast<RealType>(df1 / 2), static_cast<RealType>(df2 / 2), static_cast<RealType>(v1x / (df2 + v1x)), Policy());
}
template <class RealType, class Policy>
@@ -239,7 +239,7 @@
RealType x, y;
- x = boost::math::ibetac_inv(df1 / 2, df2 / 2, p, &y, Policy());
+ x = boost::math::ibetac_inv(static_cast<RealType>(df1 / 2), static_cast<RealType>(df2 / 2), p, &y, Policy());
return df2 * x / (df1 * y);
}
Modified: sandbox/interval_math_toolkit/boost/math/distributions/lognormal.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/distributions/lognormal.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/distributions/lognormal.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -221,7 +221,7 @@
if(0 == detail::check_scale("boost::math::variance(const lognormal_distribution<%1%>&)", sigma, &result, Policy()))
return result;
- return boost::math::expm1(sigma * sigma, Policy()) * exp(2 * mu + sigma * sigma);
+ return boost::math::expm1(static_cast<RealType>(sigma * sigma), Policy()) * exp(2 * mu + sigma * sigma);
}
template <class RealType, class Policy>
Modified: sandbox/interval_math_toolkit/boost/math/distributions/negative_binomial.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/distributions/negative_binomial.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/distributions/negative_binomial.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -431,6 +431,7 @@
// BUT Cephes/CodeCogs says: finds argument p (0 to 1) such that cdf(k, n, p) = y
static const char* function = "boost::math::quantile(const negative_binomial_distribution<%1%>&, %1%)";
BOOST_MATH_STD_USING // ADL of std functions.
+ using std::min;
RealType p = dist.success_fraction();
RealType r = dist.successes();
@@ -468,14 +469,14 @@
RealType guess = 0;
RealType factor = 5;
if(r * r * r * P * p > 0.005)
- guess = detail::inverse_negative_binomial_cornish_fisher(r, p, 1-p, P, 1-P, Policy());
+ guess = detail::inverse_negative_binomial_cornish_fisher(r, p, static_cast<RealType>(1-p), P, static_cast<RealType>(1-P), Policy());
if(guess < 10)
{
//
// Cornish-Fisher Negative binomial approximation not accurate in this area:
//
- guess = (std::min)(r * 2, RealType(10));
+ guess = min BOOST_PREVENT_MACRO_SUBSTITUTION(static_cast<RealType>(r * 2), RealType(10));
}
else
factor = (1-P < sqrt(tools::epsilon<RealType>())) ? 2 : (guess < 20 ? 1.2f : 1.1f);
@@ -503,6 +504,7 @@
// complement of the probability Q = 1 - P.
static const char* function = "boost::math::quantile(const negative_binomial_distribution<%1%>&, %1%)";
BOOST_MATH_STD_USING
+ using std::min;
// Error checks:
RealType Q = c.param;
@@ -545,14 +547,14 @@
RealType guess = 0;
RealType factor = 5;
if(r * r * r * (1-Q) * p > 0.005)
- guess = detail::inverse_negative_binomial_cornish_fisher(r, p, 1-p, 1-Q, Q, Policy());
+ guess = detail::inverse_negative_binomial_cornish_fisher(r, p, static_cast<RealType>(1-p), static_cast<RealType>(1-Q), Q, Policy());
if(guess < 10)
{
//
// Cornish-Fisher Negative binomial approximation not accurate in this area:
//
- guess = (std::min)(r * 2, RealType(10));
+ guess = min BOOST_PREVENT_MACRO_SUBSTITUTION(static_cast<RealType>(r * 2), RealType(10));
}
else
factor = (Q < sqrt(tools::epsilon<RealType>())) ? 2 : (guess < 20 ? 1.2f : 1.1f);
Modified: sandbox/interval_math_toolkit/boost/math/distributions/normal.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/distributions/normal.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/distributions/normal.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -166,7 +166,7 @@
return result;
}
RealType diff = (x - mean) / (sd * constants::root_two<RealType>());
- result = boost::math::erfc(-diff, Policy()) / 2;
+ result = boost::math::erfc(static_cast<RealType>(-diff), Policy()) / 2;
return result;
} // cdf
@@ -187,7 +187,7 @@
if(false == detail::check_probability(function, p, &result, Policy()))
return result;
- result= boost::math::erfc_inv(2 * p, Policy());
+ result= boost::math::erfc_inv(static_cast<RealType>(2 * p), Policy());
result = -result;
result *= sd * constants::root_two<RealType>();
result += mean;
@@ -247,7 +247,7 @@
RealType q = c.param;
if(false == detail::check_probability(function, q, &result, Policy()))
return result;
- result = boost::math::erfc_inv(2 * q, Policy());
+ result = boost::math::erfc_inv(static_cast<RealType>(2 * q), Policy());
result *= sd * constants::root_two<RealType>();
result += mean;
return result;
Modified: sandbox/interval_math_toolkit/boost/math/distributions/poisson.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/distributions/poisson.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/distributions/poisson.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -347,7 +347,7 @@
// (e ^ -mean * mean ^ k) / k!
// == exp(log(e ^ -mean) + log (mean ^ k) - lgamma(k+1))
// exp( -mean + log(mean) * k - lgamma(k+1))
- return exp(-mean + log(mean) * k - boost::math::lgamma(k+1, Policy()));
+ return exp(-mean + log(mean) * k - boost::math::lgamma(static_cast<RealType>(k+1), Policy()));
// return gamma_p_derivative(k+1, mean); // equivalent & also passes tests.
}
} // pdf
@@ -444,7 +444,7 @@
}
if (k == 0)
{ // Avoid repeated checks on k and mean in gamma_p.
- return -boost::math::expm1(-mean, Policy());
+ return -boost::math::expm1(static_cast<RealType>(-mean), Policy());
}
// Unlike un-complemented cdf (sum from 0 to k),
// can't use finite sum from k+1 to infinity for small integral k,
@@ -493,7 +493,7 @@
if(z < 1)
guess = z;
else
- guess = boost::math::detail::inverse_poisson_cornish_fisher(z, p, 1-p, Policy());
+ guess = boost::math::detail::inverse_poisson_cornish_fisher(z, p, static_cast<RealType>(1-p), Policy());
if(z > 5)
{
if(z > 1000)
@@ -561,7 +561,7 @@
if(z < 1)
guess = z;
else
- guess = boost::math::detail::inverse_poisson_cornish_fisher(z, 1-q, q, Policy());
+ guess = boost::math::detail::inverse_poisson_cornish_fisher(z, static_cast<RealType>(1-q), q, Policy());
if(z > 5)
{
if(z > 1000)
Modified: sandbox/interval_math_toolkit/boost/math/distributions/rayleigh.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/distributions/rayleigh.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/distributions/rayleigh.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -129,7 +129,7 @@
{
return result;
}
- result = -boost::math::expm1(-x * x / ( 2 * sigma * sigma), Policy());
+ result = -boost::math::expm1(static_cast<RealType>(-x * x / ( 2 * sigma * sigma)), Policy());
return result;
} // cdf
@@ -154,7 +154,7 @@
{
return policies::raise_overflow_error<RealType>(function, 0, Policy());
}
- result = sqrt(-2 * sigma * sigma * boost::math::log1p(-p, Policy()));
+ result = sqrt(-2 * sigma * sigma * boost::math::log1p(static_cast<RealType>(-p), Policy()));
return result;
} // quantile
Modified: sandbox/interval_math_toolkit/boost/math/distributions/triangular.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/distributions/triangular.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/distributions/triangular.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -212,11 +212,11 @@
}
if (x == lower)
{ // (mode - lower) == 0 which would lead to divide by zero!
- return (mode == lower) ? 2 / (upper - lower) : 0;
+ return (mode == lower) ? RealType(2) / (upper - lower) : RealType(0);
}
else if (x == upper)
{
- return (mode == upper) ? 2 / (upper - lower) : 0;
+ return (mode == upper) ? RealType(2) / (upper - lower) : RealType(0);
}
else if (x <= mode)
{
@@ -468,7 +468,7 @@
RealType mode = dist.mode();
RealType upper = dist.upper();
RealType result; // of checks.
- if(false == detail::check_triangular(function,lower, mode, upper, &result, Policy()))
+ if(false == boost::math::detail::check_triangular(function,lower, mode, upper, &result, Policy()))
{
return result;
}
Modified: sandbox/interval_math_toolkit/boost/math/distributions/weibull.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/distributions/weibull.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/distributions/weibull.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -158,7 +158,7 @@
if(false == detail::check_weibull_x(function, x, &result, Policy()))
return result;
- result = -boost::math::expm1(-pow(x / scale, shape), Policy());
+ result = -boost::math::expm1(static_cast<RealType>(-pow(x / scale, shape)), Policy());
return result;
}
@@ -182,7 +182,7 @@
if(p == 1)
return policies::raise_overflow_error<RealType>(function, 0, Policy());
- result = scale * pow(-boost::math::log1p(-p, Policy()), 1 / shape);
+ result = scale * pow(-boost::math::log1p(static_cast<RealType>(-p), Policy()), 1 / shape);
return result;
}
@@ -247,7 +247,7 @@
if(false == detail::check_weibull(function, scale, shape, &result, Policy()))
return result;
- result = scale * boost::math::tgamma(1 + 1 / shape, Policy());
+ result = scale * boost::math::tgamma(static_cast<RealType>(1 + 1 / shape), Policy());
return result;
}
@@ -264,9 +264,9 @@
{
return result;
}
- result = boost::math::tgamma(1 + 1 / shape, Policy());
+ result = boost::math::tgamma(static_cast<RealType>(1 + 1 / shape), Policy());
result *= -result;
- result += boost::math::tgamma(1 + 2 / shape, Policy());
+ result += boost::math::tgamma(static_cast<RealType>(1 + 2 / shape), Policy());
result *= scale * scale;
return result;
}
@@ -327,9 +327,9 @@
}
RealType g1, g2, g3, d;
- g1 = boost::math::tgamma(1 + 1 / shape, Policy());
- g2 = boost::math::tgamma(1 + 2 / shape, Policy());
- g3 = boost::math::tgamma(1 + 3 / shape, Policy());
+ g1 = boost::math::tgamma(static_cast<RealType>(1 + 1 / shape), Policy());
+ g2 = boost::math::tgamma(static_cast<RealType>(1 + 2 / shape), Policy());
+ g3 = boost::math::tgamma(static_cast<RealType>(1 + 3 / shape), Policy());
d = pow(g2 - g1 * g1, RealType(1.5));
result = (2 * g1 * g1 * g1 - 3 * g1 * g2 + g3) / d;
@@ -352,10 +352,10 @@
RealType g1, g2, g3, g4, d, g1_2, g1_4;
- g1 = boost::math::tgamma(1 + 1 / shape, Policy());
- g2 = boost::math::tgamma(1 + 2 / shape, Policy());
- g3 = boost::math::tgamma(1 + 3 / shape, Policy());
- g4 = boost::math::tgamma(1 + 4 / shape, Policy());
+ g1 = boost::math::tgamma(static_cast<RealType>(1 + 1 / shape), Policy());
+ g2 = boost::math::tgamma(static_cast<RealType>(1 + 2 / shape), Policy());
+ g3 = boost::math::tgamma(static_cast<RealType>(1 + 3 / shape), Policy());
+ g4 = boost::math::tgamma(static_cast<RealType>(1 + 4 / shape), Policy());
g1_2 = g1 * g1;
g1_4 = g1_2 * g1_2;
d = g2 - g1_2;
Modified: sandbox/interval_math_toolkit/boost/math/policies/error_handling.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/policies/error_handling.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/policies/error_handling.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -18,6 +18,7 @@
#include <boost/math/policies/policy.hpp>
#include <boost/math/tools/precision.hpp>
#include <boost/cstdint.hpp>
+#include <boost/tr1/tuple.hpp>
#ifdef BOOST_MSVC
# pragma warning(push) // Quiet warnings in boost/format.hpp
# pragma warning(disable: 4996) // _SCL_SECURE_NO_DEPRECATE
@@ -55,6 +56,17 @@
template <class T>
T user_evaluation_error(const char* function, const char* message, const T& val);
+template <class T>
+inline T make_printable_type(const T& v)
+{ return v; }
+
+template <class T, class Policy>
+inline std::tr1::tuple<T, T>
+make_printable_type(const boost::numeric::interval<T, Policy>& inter)
+{
+ return std::tr1::tuple<T, T>(inter.lower(), inter.upper());
+}
+
namespace detail
{
@@ -89,7 +101,7 @@
msg += message;
int prec = 2 + (boost::math::policies::digits<T, boost::math::policies::policy<> >() * 30103UL) / 100000UL;
- msg = (boost::format(msg) % boost::io::group(std::setprecision(prec), val)).str();
+ msg = (boost::format(msg) % boost::io::group(std::setprecision(prec), make_printable_type(val))).str();
E e(msg);
boost::throw_exception(e);
@@ -429,7 +441,7 @@
inline bool check_overflow(T val, R* result, const char* function, const Policy& pol)
{
BOOST_MATH_STD_USING
- if(fabs(val) > tools::max_value<R>())
+ if(boost::math::tools::maybe_greater(fabs(val), static_cast<T>(tools::max_value<R>())))
{
*result = static_cast<R>(boost::math::policies::detail::raise_overflow_error<R>(function, 0, pol));
return true;
@@ -439,7 +451,7 @@
template <class R, class T, class Policy>
inline bool check_underflow(T val, R* result, const char* function, const Policy& pol)
{
- if((val != 0) && (static_cast<R>(val) == 0))
+ if(!boost::math::tools::maybe_equal(val, T(0)) && boost::math::tools::maybe_equal(static_cast<R>(val), R(0)))
{
*result = static_cast<R>(boost::math::policies::detail::raise_underflow_error<R>(function, 0, pol));
return true;
@@ -450,7 +462,7 @@
inline bool check_denorm(T val, R* result, const char* function, const Policy& pol)
{
BOOST_MATH_STD_USING
- if((fabs(val) < static_cast<T>(tools::min_value<R>())) && (static_cast<R>(val) != 0))
+ if(boost::math::tools::maybe_less(fabs(val), static_cast<T>(tools::min_value<R>())) && !boost::math::tools::maybe_equal(static_cast<R>(val), R(0)))
{
*result = static_cast<R>(boost::math::policies::detail::raise_denorm_error<R>(function, 0, static_cast<R>(val), pol));
return true;
Modified: sandbox/interval_math_toolkit/boost/math/policies/policy.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/policies/policy.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/policies/policy.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -730,25 +730,26 @@
template <class Real, class Policy>
struct precision
{
+ typedef typename boost::numeric::interval_base_type<Real>::type base_type;
#ifndef __BORLANDC__
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)),
+ ((::std::numeric_limits<base_type>::is_specialized == 0) || (::std::numeric_limits<base_type>::digits == 0)),
// Possibly unknown precision:
precision_type,
typename mpl::if_c<
- ((::std::numeric_limits<Real>::digits <= precision_type::value)
+ ((::std::numeric_limits<base_type>::digits <= precision_type::value)
|| (Policy::precision_type::value <= 0)),
// Default case, full precision for RealType:
- digits2< ::std::numeric_limits<Real>::digits>,
+ digits2< ::std::numeric_limits<base_type>::digits>,
// User customised precision:
precision_type
>::type
>::type type;
#else
typedef typename Policy::precision_type precision_type;
- typedef mpl::int_< ::std::numeric_limits<Real>::digits> digits_t;
- typedef mpl::bool_< ::std::numeric_limits<Real>::is_specialized> spec_t;
+ typedef mpl::int_< ::std::numeric_limits<base_type>::digits> digits_t;
+ typedef mpl::bool_< ::std::numeric_limits<base_type>::is_specialized> spec_t;
typedef typename mpl::if_<
mpl::or_<mpl::equal_to<spec_t, mpl::false_>, mpl::equal_to<digits_t, mpl::int_<0> > >,
// Possibly unknown precision:
@@ -756,7 +757,7 @@
typename mpl::if_<
mpl::or_<mpl::less_equal<digits_t, precision_type>, mpl::less_equal<precision_type, mpl::int_<0> > >,
// Default case, full precision for RealType:
- digits2< ::std::numeric_limits<Real>::digits>,
+ digits2< ::std::numeric_limits<base_type>::digits>,
// User customised precision:
precision_type
>::type
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/bessel.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/bessel.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/bessel.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -41,7 +41,7 @@
{
BOOST_MATH_STD_USING
mult = x / 2;
- term = pow(mult, v) / boost::math::tgamma(v+1, Policy());
+ term = pow(mult, v) / boost::math::tgamma(static_cast<T>(v+1), Policy());
mult *= -mult;
}
T operator()()
@@ -68,7 +68,7 @@
{
BOOST_MATH_STD_USING
mult = x / 2;
- term = pow(mult, T(v)) / boost::math::tgamma(v+1+T(0.5f), Policy());
+ term = pow(mult, T(v)) / boost::math::tgamma(static_cast<T>(v+1+T(0.5f)), Policy());
mult *= -mult;
}
T operator()()
@@ -126,7 +126,7 @@
// better have integer v:
if(floor(v) == v)
{
- T r = cyl_bessel_j_imp(v, -x, t, pol);
+ T r = cyl_bessel_j_imp(v, static_cast<T>(-x), t, pol);
if(tools::real_cast<int>(v) & 1)
r = -r;
return r;
@@ -207,7 +207,7 @@
// Default case is just a naive evaluation of the definition:
//
return sqrt(constants::pi<T>() / (2 * x))
- * cyl_bessel_j_imp(T(n)+T(0.5f), x, bessel_no_int_tag(), pol);
+ * cyl_bessel_j_imp(static_cast<T>(T(n)+T(0.5f)), x, bessel_no_int_tag(), pol);
}
template <class T, class Policy>
@@ -225,7 +225,7 @@
// better have integer v:
if(floor(v) == v)
{
- T r = cyl_bessel_i_imp(v, -x, pol);
+ T r = cyl_bessel_i_imp(v, static_cast<T>(-x), pol);
if(tools::real_cast<int>(v) & 1)
r = -r;
return r;
@@ -378,7 +378,7 @@
if(x < 2 * tools::min_value<T>())
return -policies::raise_overflow_error<T>(function, 0, pol);
- T result = cyl_neumann_imp(T(v)+0.5f, x, bessel_no_int_tag(), pol);
+ T result = cyl_neumann_imp(static_cast<T>(T(v)+0.5f), x, bessel_no_int_tag(), pol);
T tx = sqrt(constants::pi<T>() / (2 * x));
if((tx > 1) && (tools::max_value<T>() / tx < result))
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/beta.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/beta.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/beta.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -89,7 +89,7 @@
{
// Special case where the base of the power term is close to 1
// compute (1+x)^y instead:
- result *= exp(ambh * boost::math::log1p(-b / cgh, pol));
+ result *= exp(ambh * boost::math::log1p(static_cast<T>(-b / cgh), pol));
}
else
{
@@ -154,9 +154,10 @@
std::swap(a, b);
// set integration limits:
- T la = (std::max)(T(10), a);
- T lb = (std::max)(T(10), b);
- T lc = (std::max)(T(10), a+b);
+ using std::max;
+ T la = max BOOST_PREVENT_MACRO_SUBSTITUTION(T(10), a);
+ T lb = max BOOST_PREVENT_MACRO_SUBSTITUTION(T(10), b);
+ T lc = max BOOST_PREVENT_MACRO_SUBSTITUTION(T(10), static_cast<T>(a+b));
// calculate the fraction parts:
T sa = detail::lower_gamma_series(a, la, pol) / a;
@@ -200,6 +201,8 @@
const Policy& pol)
{
BOOST_MATH_STD_USING
+ using std::min;
+ using std::max;
if(!normalised)
{
@@ -221,11 +224,11 @@
// l1 and l2 are the base of the exponents minus one:
T l1 = (x * b - y * agh) / agh;
T l2 = (y * a - x * bgh) / bgh;
- if(((std::min)(fabs(l1), fabs(l2)) < 0.2))
+ if((min BOOST_PREVENT_MACRO_SUBSTITUTION(fabs(l1), fabs(l2)) < 0.2))
{
// when the base of the exponent is very near 1 we get really
// gross errors unless extra care is taken:
- if((l1 * l2 > 0) || ((std::min)(a, b) < 1))
+ if((l1 * l2 > 0) || (min BOOST_PREVENT_MACRO_SUBSTITUTION(a, b) < 1))
{
//
// This first branch handles the simple cases where either:
@@ -249,7 +252,7 @@
else
result *= pow((y * cgh) / bgh, b);
}
- else if((std::max)(fabs(l1), fabs(l2)) < 0.5)
+ else if(max BOOST_PREVENT_MACRO_SUBSTITUTION(fabs(l1), fabs(l2)) < 0.5)
{
//
// Both exponents are near one and both the exponents are
@@ -274,14 +277,14 @@
T ratio = b / a;
if((small_a && (ratio * l2 < 0.1)) || (!small_a && (l1 / ratio > 0.1)))
{
- T l3 = boost::math::expm1(ratio * boost::math::log1p(l2, pol), pol);
+ T l3 = boost::math::expm1(static_cast<T>(ratio * boost::math::log1p(l2, pol)), pol);
l3 = l1 + l3 + l3 * l1;
l3 = a * boost::math::log1p(l3, pol);
result *= exp(l3);
}
else
{
- T l3 = boost::math::expm1(boost::math::log1p(l1, pol) / ratio, pol);
+ T l3 = boost::math::expm1(static_cast<T>(boost::math::log1p(l1, pol) / ratio), pol);
l3 = l2 + l3 + l3 * l2;
l3 = b * boost::math::log1p(l3, pol);
result *= exp(l3);
@@ -460,7 +463,7 @@
T cgh = c + L::g() - T(0.5);
result = L::lanczos_sum_expG_scaled(c) / (L::lanczos_sum_expG_scaled(a) * L::lanczos_sum_expG_scaled(b));
if(a * b < bgh * 10)
- result *= exp((b - 0.5f) * boost::math::log1p(a / bgh, pol));
+ result *= exp((b - 0.5f) * boost::math::log1p(static_cast<T>(a / bgh), pol));
else
result *= pow(cgh / bgh, b - 0.5f);
result *= pow(x * cgh / agh, a);
@@ -717,7 +720,7 @@
T t = a + bm1 / 2;
T lx, u;
if(y < 0.35)
- lx = boost::math::log1p(-y, pol);
+ lx = boost::math::log1p(static_cast<T>(-y), pol);
else
lx = log(x);
u = -t * lx;
@@ -821,7 +824,7 @@
BOOST_MATH_STD_USING // ADL of std names
T result = pow(x, n);
T term = result;
- for(unsigned i = tools::real_cast<unsigned>(n - 1); i > k; --i)
+ for(unsigned i = tools::real_cast<unsigned>(static_cast<T>(n - 1)); i > k; --i)
{
term *= ((i + 1) * y) / ((n - i) * x) ;
result += term;
@@ -843,6 +846,8 @@
static const char* function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)";
typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
BOOST_MATH_STD_USING // for ADL of std math functions.
+ using std::max;
+ using std::min;
bool invert = inv;
T fract;
@@ -873,7 +878,7 @@
{
if(p_derivative)
{
- *p_derivative = (a == 1) ? 1 : (a < 1) ? tools::max_value<T>() / 2 : tools::min_value<T>() * 2;
+ *p_derivative = (a == 1) ? 1 : (a < 1) ? static_cast<T>(tools::max_value<T>() / 2) : static_cast<T>(tools::min_value<T>() * 2);
}
return (invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) : 0);
}
@@ -881,12 +886,12 @@
{
if(p_derivative)
{
- *p_derivative = (b == 1) ? 1 : (b < 1) ? tools::max_value<T>() / 2 : tools::min_value<T>() * 2;
+ *p_derivative = (b == 1) ? 1 : (b < 1) ? static_cast<T>(tools::max_value<T>() / 2) : static_cast<T>(tools::min_value<T>() * 2);
}
return (invert == 0 ? (normalised ? 1 : boost::math::beta(a, b, pol)) : 0);
}
- if((std::min)(a, b) <= 1)
+ if(min BOOST_PREVENT_MACRO_SUBSTITUTION(a, b) <= 1)
{
if(x > 0.5)
{
@@ -894,10 +899,10 @@
std::swap(x, y);
invert = !invert;
}
- if((std::max)(a, b) <= 1)
+ if(max BOOST_PREVENT_MACRO_SUBSTITUTION(a, b) <= 1)
{
// Both a,b < 1:
- if((a >= (std::min)(T(0.2), b)) || (pow(x, a) <= 0.9))
+ if((a >= min BOOST_PREVENT_MACRO_SUBSTITUTION(T(0.2), b)) || (pow(x, a) <= 0.9))
{
if(!invert)
fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
@@ -930,7 +935,7 @@
T prefix;
if(!normalised)
{
- prefix = rising_factorial_ratio(a+b, a, 20);
+ prefix = rising_factorial_ratio(static_cast<T>(a+b), a, 20);
}
else
{
@@ -938,12 +943,12 @@
}
fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative);
if(!invert)
- fract = beta_small_b_large_a_series(a + 20, b, x, y, fract, prefix, pol, normalised);
+ fract = beta_small_b_large_a_series(static_cast<T>(a + 20), b, x, y, fract, prefix, pol, normalised);
else
{
fract -= (normalised ? 1 : boost::math::beta(a, b, pol));
invert = false;
- fract = -beta_small_b_large_a_series(a + 20, b, x, y, fract, prefix, pol, normalised);
+ fract = -beta_small_b_large_a_series(static_cast<T>(a + 20), b, x, y, fract, prefix, pol, normalised);
}
}
}
@@ -996,7 +1001,7 @@
T prefix;
if(!normalised)
{
- prefix = rising_factorial_ratio(a+b, a, 20);
+ prefix = rising_factorial_ratio(static_cast<T>(a+b), a, 20);
}
else
{
@@ -1004,12 +1009,12 @@
}
fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative);
if(!invert)
- fract = beta_small_b_large_a_series(a + 20, b, x, y, fract, prefix, pol, normalised);
+ fract = beta_small_b_large_a_series(static_cast<T>(a + 20), b, x, y, fract, prefix, pol, normalised);
else
{
fract -= (normalised ? 1 : boost::math::beta(a, b, pol));
invert = false;
- fract = -beta_small_b_large_a_series(a + 20, b, x, y, fract, prefix, pol, normalised);
+ fract = -beta_small_b_large_a_series(static_cast<T>(a + 20), b, x, y, fract, prefix, pol, normalised);
}
}
}
@@ -1059,14 +1064,14 @@
else if(a > 15)
{
// sidestep so we can use the series representation:
- int n = static_cast<int>(boost::math::tools::real_cast<long double>(floor(b)));
+ int n = tools::real_cast<int>(static_cast<T>(floor(b)));
if(n == b)
--n;
T bbar = b - n;
T prefix;
if(!normalised)
{
- prefix = rising_factorial_ratio(a+bbar, bbar, n);
+ prefix = rising_factorial_ratio(static_cast<T>(a+bbar), bbar, n);
}
else
{
@@ -1081,7 +1086,7 @@
// the formula here for the non-normalised case is tricky to figure
// out (for me!!), and requires two pochhammer calculations rather
// than one, so leave it for now....
- int n = static_cast<int>(boost::math::tools::real_cast<long double>(floor(b)));
+ int n = tools::real_cast<int>(static_cast<T>(floor(b)));
T bbar = b - n;
if(bbar <= 0)
{
@@ -1093,7 +1098,7 @@
if(invert)
fract -= (normalised ? 1 : boost::math::beta(a, b, pol));
//fract = ibeta_series(a+20, bbar, x, fract, l, normalised, p_derivative, y);
- fract = beta_small_b_large_a_series(a+20, bbar, x, y, fract, T(1), pol, normalised);
+ fract = beta_small_b_large_a_series(static_cast<T>(a+20), bbar, x, y, fract, T(1), pol, normalised);
if(invert)
{
fract = -fract;
@@ -1166,7 +1171,7 @@
// Now the regular cases:
//
typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
- T f1 = ibeta_power_terms(a, b, x, 1 - x, lanczos_type(), true, pol);
+ T f1 = ibeta_power_terms(a, b, x, static_cast<T>(1 - x), lanczos_type(), true, pol);
T y = (1 - x) * x;
if(f1 == 0)
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -42,17 +42,17 @@
BOOST_ASSERT(abs(v) <= 0.5f);
T gp = boost::math::tgamma1pm1(v, pol);
- T gm = boost::math::tgamma1pm1(-v, pol);
+ T gm = boost::math::tgamma1pm1(static_cast<T>(-v), pol);
a = log(x / 2);
b = exp(v * a);
sigma = -a * v;
c = abs(v) < tools::epsilon<T>() ?
- 1 : boost::math::sin_pi(v) / (v * pi<T>());
+ static_cast<T>(1) : static_cast<T>(boost::math::sin_pi(v) / (v * pi<T>()));
d = abs(sigma) < tools::epsilon<T>() ?
- 1 : sinh(sigma) / sigma;
+ static_cast<T>(1) : static_cast<T>(sinh(sigma) / sigma);
gamma1 = abs(v) < tools::epsilon<T>() ?
- -euler<T>() : (0.5f / v) * (gp - gm) * c;
+ static_cast<T>(-euler<T>()) : static_cast<T>((0.5f / v) * (gp - gm) * c);
gamma2 = (2 + gp + gm) * c / 2;
// initial values
@@ -234,7 +234,7 @@
v = -v; // v is non-negative from here
kind |= need_k;
}
- n = tools::real_cast<unsigned>(v + 0.5f);
+ n = tools::real_cast<unsigned>(static_cast<T>(v + 0.5f));
u = v - n; // -1/2 <= u < 1/2
BOOST_MATH_INSTRUMENT_VARIABLE(n);
BOOST_MATH_INSTRUMENT_VARIABLE(u);
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/detail/bessel_jy.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/detail/bessel_jy.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/detail/bessel_jy.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -46,21 +46,21 @@
BOOST_ASSERT(fabs(v) <= 0.5f); // precondition for using this routine
T gp = boost::math::tgamma1pm1(v, pol);
- T gm = boost::math::tgamma1pm1(-v, pol);
+ T gm = boost::math::tgamma1pm1(static_cast<T>(-v), pol);
T spv = boost::math::sin_pi(v, pol);
- T spv2 = boost::math::sin_pi(v/2, pol);
+ T spv2 = boost::math::sin_pi(static_cast<T>(v/2), pol);
T xp = pow(x/2, v);
a = log(x / 2);
sigma = -a * v;
d = abs(sigma) < tools::epsilon<T>() ?
T(1) : sinh(sigma) / sigma;
- e = abs(v) < tools::epsilon<T>() ? v*pi<T>()*pi<T>() / 2
- : 2 * spv2 * spv2 / v;
+ e = abs(v) < tools::epsilon<T>() ? static_cast<T>(v*pi<T>()*pi<T>() / 2)
+ : static_cast<T>(2 * spv2 * spv2 / v);
- T g1 = (v == 0) ? -euler<T>() : (gp - gm) / ((1 + gp) * (1 + gm) * 2 * v);
+ T g1 = (v == 0) ? static_cast<T>(-euler<T>()) : static_cast<T>((gp - gm) / ((1 + gp) * (1 + gm) * 2 * v));
T g2 = (2 + gp + gm) / ((1 + gp) * (1 + gm) * 2);
- T vspv = (fabs(v) < tools::epsilon<T>()) ? 1/constants::pi<T>() : v / spv;
+ T vspv = (fabs(v) < tools::epsilon<T>()) ? static_cast<T>(1/constants::pi<T>()) : static_cast<T>(v / spv);
f = (g1 * cosh(sigma) - g2 * a * d) * 2 * vspv;
p = vspv / (xp * (1 + gm));
@@ -118,7 +118,7 @@
tolerance = 2 * tools::epsilon<T>();
tiny = sqrt(tools::min_value<T>());
C = f = tiny; // b0 = 0, replace with tiny
- D = 0.0L;
+ D = 0;
for (k = 1; k < policies::get_max_series_iterations<Policy>() * 100; k++)
{
a = -1;
@@ -131,7 +131,7 @@
delta = C * D;
f *= delta;
if (D < 0) { s = -s; }
- if (abs(delta - 1.0L) < tolerance)
+ if (abs(delta - 1) < tolerance)
{ break; }
}
policies::check_series_iterations("boost::math::bessel_jy<%1%>(%1%,%1%) in CF1_jy", k / 100, pol);
@@ -158,7 +158,7 @@
typedef typename complex_trait<T>::type complex_type;
complex_type C, D, f, a, b, delta, one(1);
- T tiny, zero(0.0L);
+ T tiny, zero(0);
unsigned long k;
// |x| >= |v|, CF2_jy converges rapidly
@@ -169,7 +169,7 @@
// Lentz, Applied Optics, vol 15, 668 (1976)
T tolerance = 2 * tools::epsilon<T>();
tiny = sqrt(tools::min_value<T>());
- C = f = complex_type(-0.5f/x, 1.0L);
+ C = f = complex_type(static_cast<T>(-0.5f/x), 1);
D = 0;
for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)
{
@@ -225,7 +225,7 @@
v = -v; // v is non-negative from here
kind = need_j|need_y; // need both for reflection formula
}
- n = real_cast<unsigned>(v + 0.5L);
+ n = real_cast<unsigned>(static_cast<T>(v + 0.5f));
u = v - n; // -1/2 <= u < 1/2
if (x == 0)
@@ -279,7 +279,8 @@
lim = asymptotic_bessel_y_limit<T>(tag_type());
break;
default:
- lim = (std::max)(
+ using std::max;
+ lim = max BOOST_PREVENT_MACRO_SUBSTITUTION(
asymptotic_bessel_j_limit<T>(v, tag_type()),
asymptotic_bessel_y_limit<T>(tag_type()));
break;
@@ -289,7 +290,7 @@
if(kind&need_y)
{
Yu = asymptotic_bessel_y_large_x_2(u, x);
- Yu1 = asymptotic_bessel_y_large_x_2(u + 1, x);
+ Yu1 = asymptotic_bessel_y_large_x_2(static_cast<T>(u + 1), x);
}
else
Yu = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/detail/bessel_jy_asym.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/detail/bessel_jy_asym.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/detail/bessel_jy_asym.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -208,28 +208,32 @@
{
// default case:
BOOST_MATH_STD_USING
- T v2 = (std::max)(T(3), v * v);
+ using std::max;
+ T v2 = max BOOST_PREVENT_MACRO_SUBSTITUTION(T(3), static_cast<T>(v * v));
return v2 / pow(100 * tools::epsilon<T>() / T(2e-5f), T(0.17f));
}
template <class T>
inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<53>&)
{
// double case:
- T v2 = (std::max)(T(3), v * v);
+ using std::max;
+ T v2 = max BOOST_PREVENT_MACRO_SUBSTITUTION(T(3), v * v);
return v2 * 33 /*73*/;
}
template <class T>
inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<64>&)
{
// 80-bit extended-double case:
- T v2 = (std::max)(T(3), v * v);
+ using std::max;
+ T v2 = max BOOST_PREVENT_MACRO_SUBSTITUTION(T(3), v * v);
return v2 * 121 /*266*/;
}
template <class T>
inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<113>&)
{
// 128-bit long double case:
- T v2 = (std::max)(T(3), v * v);
+ using std::max;
+ T v2 = max BOOST_PREVENT_MACRO_SUBSTITUTION(T(3), v * v);
return v2 * 39154 /*85700*/;
}
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/detail/gamma_inva.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/detail/gamma_inva.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/detail/gamma_inva.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -48,7 +48,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, pol) * constants::root_two<T>();
+ T x = boost::math::erfc_inv(p > q ? static_cast<T>(2 * q) : static_cast<T>(2 * p), pol) * constants::root_two<T>();
// Set the sign:
if(p < 0.5)
x = -x;
@@ -153,25 +153,25 @@
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);
+ return detail::gamma_inva_imp(x, p, static_cast<T>(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);
+ return detail::gamma_inva_imp(x, static_cast<T>(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, policies::policy<>());
+ return detail::gamma_inva_imp(x, p, static_cast<T>(1 - p), policies::policy<>());
}
template <class T>
inline T gamma_q_inva(T x, T q)
{
- return detail::gamma_inva_imp(x, 1 - q, q, policies::policy<>());
+ return detail::gamma_inva_imp(x, static_cast<T>(1 - q), q, policies::policy<>());
}
} // namespace math
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/detail/ibeta_inv_ab.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/detail/ibeta_inv_ab.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/detail/ibeta_inv_ab.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -51,7 +51,7 @@
// kurtosis:
T k = (6 - sf * (5+sfc)) / (n * (sfc));
// Get the inverse of a std normal distribution:
- T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>();
+ T x = boost::math::erfc_inv(p > q ? static_cast<T>(2 * q) : static_cast<T>(2 * p), pol) * constants::root_two<T>();
// Set the sign:
if(p < 0.5)
x = -x;
@@ -74,6 +74,7 @@
T ibeta_inv_ab_imp(const T& b, const T& z, const T& p, const T& q, bool swap_ab, const Policy& pol)
{
BOOST_MATH_STD_USING // for ADL of std lib math functions
+ using std::min;
//
// Special cases first:
//
@@ -120,11 +121,11 @@
//
if((p < q) != swap_ab)
{
- guess = (std::min)(b * 2, T(1));
+ guess = min BOOST_PREVENT_MACRO_SUBSTITUTION(static_cast<T>(b * 2), T(1));
}
else
{
- guess = (std::min)(b / 2, T(1));
+ guess = min BOOST_PREVENT_MACRO_SUBSTITUTION(static_cast<T>(b / 2), T(1));
}
}
if(n * n * n * u * sf > 0.005)
@@ -137,11 +138,11 @@
//
if((p < q) != swap_ab)
{
- guess = (std::min)(b * 2, T(10));
+ guess = min BOOST_PREVENT_MACRO_SUBSTITUTION(static_cast<T>(b * 2), T(10));
}
else
{
- guess = (std::min)(b / 2, T(10));
+ guess = min BOOST_PREVENT_MACRO_SUBSTITUTION(static_cast<T>(b / 2), T(10));
}
}
else
@@ -162,49 +163,49 @@
template <class T, class Policy>
inline T ibeta_inva(T b, T x, T p, const Policy& pol)
{
- return detail::ibeta_inv_ab_imp(b, x, p, 1 - p, false, pol);
+ return detail::ibeta_inv_ab_imp(b, x, p, static_cast<T>(1 - p), false, pol);
}
template <class T, class Policy>
inline T ibetac_inva(T b, T x, T q, const Policy& pol)
{
- return detail::ibeta_inv_ab_imp(b, x, 1 - q, q, false, pol);
+ return detail::ibeta_inv_ab_imp(b, x, static_cast<T>(1 - q), q, false, pol);
}
template <class T, class Policy>
inline T ibeta_invb(T b, T x, T p, const Policy& pol)
{
- return detail::ibeta_inv_ab_imp(b, x, p, 1 - p, true, pol);
+ return detail::ibeta_inv_ab_imp(b, x, p, static_cast<T>(1 - p), true, pol);
}
template <class T, class Policy>
inline T ibetac_invb(T b, T x, T q, const Policy& pol)
{
- return detail::ibeta_inv_ab_imp(b, x, 1 - q, q, true, pol);
+ return detail::ibeta_inv_ab_imp(b, x, static_cast<T>(1 - q), q, true, pol);
}
template <class T>
inline T ibeta_inva(T b, T x, T p)
{
- return detail::ibeta_inv_ab_imp(b, x, p, 1 - p, false, policies::policy<>());
+ return detail::ibeta_inv_ab_imp(b, x, p, static_cast<T>(1 - p), false, policies::policy<>());
}
template <class T>
inline T ibetac_inva(T b, T x, T q)
{
- return detail::ibeta_inv_ab_imp(b, x, 1 - q, q, false, policies::policy<>());
+ return detail::ibeta_inv_ab_imp(b, x, static_cast<T>(1 - q), q, false, policies::policy<>());
}
template <class T>
inline T ibeta_invb(T b, T x, T p)
{
- return detail::ibeta_inv_ab_imp(b, x, p, 1 - p, true, policies::policy<>());
+ return detail::ibeta_inv_ab_imp(b, x, p, static_cast<T>(1 - p), true, policies::policy<>());
}
template <class T>
inline T ibetac_invb(T b, T x, T q)
{
- return detail::ibeta_inv_ab_imp(b, x, 1 - q, q, true, policies::policy<>());
+ return detail::ibeta_inv_ab_imp(b, x, static_cast<T>(1 - q), q, true, policies::policy<>());
}
} // namespace math
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/detail/ibeta_inverse.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/detail/ibeta_inverse.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/detail/ibeta_inverse.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -66,7 +66,7 @@
// get the first approximation for eta from the inverse
// error function (Eq: 2.9 and 2.10).
//
- T eta0 = boost::math::erfc_inv(2 * z, pol);
+ T eta0 = boost::math::erfc_inv(static_cast<T>(2 * z), pol);
eta0 /= -sqrt(a / 2);
T terms[4] = { eta0 };
@@ -106,7 +106,7 @@
//
// Bring them together to get a final estimate for eta:
//
- T eta = tools::evaluate_polynomial(terms, 1/a, 4);
+ T eta = tools::evaluate_polynomial(terms, static_cast<T>(1/a), 4);
//
// now we need to convert eta to x, by solving the appropriate
// quadratic equation:
@@ -143,7 +143,7 @@
// Get first estimate for eta, see Eq 3.9 and 3.10,
// but note there is a typo in Eq 3.10:
//
- T eta0 = boost::math::erfc_inv(2 * z, pol);
+ T eta0 = boost::math::erfc_inv(static_cast<T>(2 * z), pol);
eta0 /= -sqrt(r / 2);
T s = sin(theta);
@@ -207,7 +207,7 @@
// Bring the correction terms together to evaluate eta,
// this is the last equation on page 151:
//
- T eta = tools::evaluate_polynomial(terms, 1/r, 4);
+ T eta = tools::evaluate_polynomial(terms, static_cast<T>(1/r), 4);
//
// Now that we have eta we need to back solve for x,
// we seek the value of x that gives eta in Eq 3.2.
@@ -453,6 +453,8 @@
T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
{
BOOST_MATH_STD_USING // For ADL of math functions.
+ using std::max;
+ using std::min;
//
// The flag invert is set to true if we swap a for b and p for q,
@@ -521,8 +523,8 @@
std::swap(p, q);
invert = !invert;
}
- T minv = (std::min)(a, b);
- T maxv = (std::max)(a, b);
+ T minv = min BOOST_PREVENT_MACRO_SUBSTITUTION(a, b);
+ T maxv = max BOOST_PREVENT_MACRO_SUBSTITUTION(a, b);
if((sqrt(minv) > (maxv - minv)) && (minv > 5))
{
//
@@ -660,7 +662,7 @@
//
T lx = log(p * a * boost::math::beta(a, b, pol)) / a;
x = exp(lx);
- y = x < 0.9 ? 1 - x : -boost::math::expm1(lx, pol);
+ y = x < 0.9 ? static_cast<T>(1 - x) : static_cast<T>(-boost::math::expm1(lx, pol));
if((b < a) && (x < 0.2))
{
@@ -734,7 +736,6 @@
x = 1 - y;
}
}
-
//
// Now we have a guess for x (and for y) we can set things up for
// iteration. If x > 0.5 it pays to swap things round:
@@ -804,7 +805,7 @@
//BOOST_ASSERT(x != upper);
//BOOST_ASSERT((x != lower) || (x == boost::math::tools::min_value<T>()) || (x == boost::math::tools::epsilon<T>()));
//
- // Tidy up, if we "lower" was too high then zero is the best answer we have:
+ // Tidy up, if "lower" was too high then zero is the best answer we have:
//
if(x == lower)
x = 0;
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/detail/igamma_inverse.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/detail/igamma_inverse.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/detail/igamma_inverse.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -88,7 +88,7 @@
// See equation 34.
//
BOOST_MATH_STD_USING
- T u = log(p) + boost::math::lgamma(a + 1, pol);
+ T u = log(p) + boost::math::lgamma(static_cast<T>(a + 1), pol);
return exp((u + x - log(didonato_SN(a, x, N, tolerance))) / a);
}
@@ -105,6 +105,7 @@
// December 1986, Pages 377-393.
//
BOOST_MATH_STD_USING
+ using std::max;
T result;
@@ -209,7 +210,7 @@
}
else
{
- T D = (std::max)(T(2), a * (a - 1));
+ T D = max BOOST_PREVENT_MACRO_SUBSTITUTION(T(2), static_cast<T>(a * (a - 1)));
T lg = boost::math::lgamma(a, pol);
T lb = log(q) + lg;
if(lb < -D * 2.3)
@@ -261,7 +262,7 @@
{
// DiDonato and Morris Eq 36:
T zb = didonato_FN(p, a, z, 100, T(1e-4), pol);
- T u = log(p) + boost::math::lgamma(a + 1, pol);
+ T u = log(p) + boost::math::lgamma(static_cast<T>(a + 1), pol);
result = zb * (1 - (a * log(zb) - zb - u + log(didonato_SN(a, z, 100, T(1e-4)))) / (a - zb));
}
}
@@ -357,7 +358,7 @@
return tools::max_value<T>();
if(p == 0)
return 0;
- T guess = detail::find_inverse_gamma(a, p, 1 - p, pol);
+ T guess = detail::find_inverse_gamma(a, p, static_cast<T>(1 - p), pol);
T lower = tools::min_value<T>();
if(guess <= lower)
guess = tools::min_value<T>();
@@ -399,7 +400,7 @@
return tools::max_value<T>();
if(q == 1)
return 0;
- T guess = detail::find_inverse_gamma(a, 1 - q, q, pol);
+ T guess = detail::find_inverse_gamma(a, static_cast<T>(1 - q), q, pol);
T lower = tools::min_value<T>();
if(guess <= lower)
guess = tools::min_value<T>();
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/detail/lgamma_small.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/detail/lgamma_small.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/detail/lgamma_small.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -492,7 +492,7 @@
// 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(static_cast<T>(dz / (L::g() + T(1.5))), pol) * T(1.5);
result += boost::math::log1p(L::lanczos_sum_near_2(dz), pol);
}
else
@@ -500,7 +500,7 @@
// 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(static_cast<T>(dz / (L::g() + T(0.5))), pol) / 2;
result += boost::math::log1p(L::lanczos_sum_near_1(dz), pol);
}
return result;
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -30,7 +30,7 @@
T a, b, c, d, q, x, y;
if (ndf > 1e20f)
- return -boost::math::erfc_inv(2 * u, pol) * constants::root_two<T>();
+ return -boost::math::erfc_inv(static_cast<T>(2 * u), pol) * constants::root_two<T>();
a = 1 / (ndf - 0.5f);
b = 48 / (a * a);
@@ -43,14 +43,14 @@
//
// Asymptotic inverse expansion about normal:
//
- x = -boost::math::erfc_inv(2 * u, pol) * constants::root_two<T>();
+ x = -boost::math::erfc_inv(static_cast<T>(2 * u), pol) * constants::root_two<T>();
y = x * x;
if (ndf < 5)
c += 0.3f * (ndf - 4.5f) * (x + 0.6f);
c += (((0.05f * d * x - 5) * x - 7) * x - 2) * x + b;
y = (((((0.4f * y + 6.3f) * y + 36) * y + 94.5f) / c - y - 3) / b + 1) * x;
- y = boost::math::expm1(a * y * y, pol);
+ y = boost::math::expm1(static_cast<T>(a * y * y), pol);
}
else
{
@@ -141,35 +141,60 @@
// Figure out what the coefficients are, note these depend
// only on the degrees of freedom (Eq 57 of Shaw):
//
- c[2] = T(1) / 6 + T(1) / (6 * df);
+ c[2] = 0.16666666666666666667 + 0.16666666666666666667 / df;
T in = 1 / df;
- c[3] = (((T(1) / 120) * in) + (T(1) / 15)) * in + (T(7) / 120);
- c[4] = ((((T(1) / 5040) * in + (T(1) / 560)) * in + (T(3) / 112)) * in + T(127) / 5040);
- c[5] = ((((T(1) / 362880) * in + (T(17) / 45360)) * in + (T(67) / 60480)) * in + (T(479) / 45360)) * in + (T(4369) / 362880);
- c[6] = ((((((T(1) / 39916800) * in + (T(2503) / 39916800)) * in + (T(11867) / 19958400)) * in + (T(1285) / 798336)) * in + (T(153161) / 39916800)) * in + (T(34807) / 5702400));
- c[7] = (((((((T(1) / 6227020800LL) * in + (T(37) / 2402400)) * in +
- (T(339929) / 2075673600LL)) * in + (T(67217) / 97297200)) * in +
- (T(870341) / 691891200LL)) * in + (T(70691) / 64864800LL)) * in +
- (T(20036983LL) / 6227020800LL));
- c[8] = (((((((T(1) / 1307674368000LL) * in + (T(1042243LL) / 261534873600LL)) * in +
- (T(21470159) / 435891456000LL)) * in + (T(326228899LL) / 1307674368000LL)) * in +
- (T(843620579) / 1307674368000LL)) * in + (T(332346031LL) / 435891456000LL)) * in +
- (T(43847599) / 1307674368000LL)) * in + (T(2280356863LL) / 1307674368000LL);
- c[9] = (((((((((T(1) / 355687428096000LL)) * in + (T(24262727LL) / 22230464256000LL)) * in +
- (T(123706507LL) / 8083805184000LL)) * in + (T(404003599LL) / 4446092851200LL)) * in +
- (T(51811946317LL) / 177843714048000LL)) * in + (T(91423417LL) / 177843714048LL)) * in +
- (T(32285445833LL) / 88921857024000LL)) * in + (T(531839683LL) / 1710035712000LL)) * in +
- (T(49020204823LL) / 50812489728000LL);
- c[10] = (((((((((T(1) / 121645100408832000LL) * in +
- (T(4222378423LL) / 13516122267648000LL)) * in +
- (T(49573465457LL) / 10137091700736000LL)) * in +
- (T(176126809LL) / 5304600576000LL)) * in +
- (T(44978231873LL) / 355687428096000LL)) * in +
- (T(5816850595639LL) / 20274183401472000LL)) * in +
- (T(73989712601LL) / 206879422464000LL)) * in +
- (T(26591354017LL) / 259925428224000LL)) * in +
- (T(14979648446341LL) / 40548366802944000LL)) * in +
- (T(65967241200001LL) / 121645100408832000LL);
+ c[3] = (0.0083333333333333333333 * in
+ + 0.066666666666666666667) * in
+ + 0.058333333333333333333;
+ c[4] = ((0.00019841269841269841270 * in
+ + 0.0017857142857142857143) * in
+ + 0.026785714285714285714) * in
+ + 0.025198412698412698413;
+ c[5] = (((2.7557319223985890653e10-6 * in
+ + 0.00037477954144620811287) * in
+ - 0.0011078042328042328042) * in
+ + 0.010559964726631393298) * in
+ + 0.012039792768959435626;
+ c[6] = ((((2.5052108385441718775e-8 * in
+ - 0.000062705427288760622094) * in
+ + 0.00059458674042007375341) * in
+ - 0.0016095979637646304313) * in
+ + 0.0061039211560044893378) * in
+ + 0.0038370059724226390893;
+ c[7] = (((((1.6059043836821614599e-10 * in
+ + 0.000015401265401265401265) * in
+ - 0.00016376804137220803887) * in
+ + 0.00069084207973096861986) * in
+ - 0.0012579159844784844785) * in
+ + 0.0010898206731540064873) * in
+ + 0.0032177478835464946576;
+ c[8] = ((((((7.6471637318198164759e-13 * in
+ - 3.9851014346715404916e-6) * in
+ + 0.000049255746366361445727) * in
+ - 0.00024947258047043099953) * in
+ + 0.00064513046951456342991) * in
+ - 0.00076245135440323932387) * in
+ + 0.000033530976880017885309) * in
+ + 0.0017438262298340009980;
+ c[9] = (((((((2.8114572543455207632e-15 * in
+ + 1.0914179173496789432e-6) * in
+ - 0.000015303004486655377567) * in
+ + 0.000090867107935219902229) * in
+ - 0.00029133414466938067350) * in
+ + 0.00051406605788341121363) * in
+ - 0.00036307660358786885787) * in
+ - 0.00031101086326318780412) * in
+ + 0.00096472747321388644237;
+ c[10] = ((((((((8.2206352466243297170e-18 * in
+ - 3.1239569599829868045e-7) * in
+ + 4.8903045291975346210e-6) * in
+ - 0.000033202652391372058698) * in
+ + 0.00012645437628698076975) * in
+ - 0.00028690924218514613987) * in
+ + 0.00035764655430568632777) * in
+ - 0.00010230378073700412687) * in
+ - 0.00036942667800009661203) * in
+ + 0.00054229262813129686486;
//
// The result is then a polynomial in v (see Eq 56 of Shaw):
//
@@ -258,7 +283,7 @@
//
T a = 4 * (u - u * u);//1 - 4 * (u - 0.5f) * (u - 0.5f);
T b = boost::math::cbrt(a);
- static const T c = 0.85498797333834849467655443627193L;
+ static const T c = 0.85498797333834849467655443627193;
T p = 6 * (1 + c * (1 / b - 1));
T p0;
do{
@@ -368,7 +393,7 @@
// where we use Shaw's tail series.
// The crossover point is roughly exponential in -df:
//
- T crossover = ldexp(1.0f, tools::real_cast<int>(df / -0.654f));
+ T crossover = ldexp(1.0f, tools::real_cast<int>(static_cast<T>(df / -0.654f)));
if(u > crossover)
{
result = boost::math::detail::inverse_students_t_hill(df, u, pol);
@@ -385,7 +410,7 @@
template <class T, class Policy>
inline T find_ibeta_inv_from_t_dist(T a, T p, T q, T* py, const Policy& pol)
{
- T u = (p > q) ? 0.5f - q / 2 : p / 2;
+ T u = (p > q) ? static_cast<T>(0.5f - q / 2) : static_cast<T>(p / 2);
T v = 1 - u; // u < 0.5 so no cancellation error
T df = a * 2;
T t = boost::math::detail::inverse_students_t(df, u, v, pol);
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/detail/unchecked_factorial.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/detail/unchecked_factorial.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/detail/unchecked_factorial.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -279,8 +279,16 @@
value = ::boost::math::max_factorial<long double>::value);
};
+namespace detail{
+
template <class T>
-inline T unchecked_factorial(unsigned i BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T))
+inline T unchecked_factorial_imp(const mpl::true_&, unsigned i BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T))
+{
+ return boost::math::unchecked_factorial<typename T::base_type>(i);
+}
+
+template <class T>
+inline T unchecked_factorial_imp(const mpl::false_&, unsigned i BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T))
{
static const boost::array<T, 101> factorials = {{
boost::lexical_cast<T>("1"),
@@ -389,6 +397,15 @@
return factorials[i];
}
+} // namespace detail
+
+template <class T>
+inline T unchecked_factorial(unsigned i BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T))
+{
+ typedef typename boost::numeric::is_interval<T>::type tag_type;
+ return detail::unchecked_factorial_imp<T>(tag_type(), i);
+}
+
template <class T>
struct max_factorial
{
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/ellint_1.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/ellint_1.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/ellint_1.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -88,9 +88,9 @@
// so rewritten to use fmod instead:
//
BOOST_MATH_INSTRUMENT_CODE("pi/2 = " << constants::pi<T>() / 2);
- T rphi = boost::math::tools::fmod_workaround(phi, constants::pi<T>() / 2);
+ T rphi = boost::math::tools::fmod_workaround(phi, static_cast<T>(constants::pi<T>() / 2));
BOOST_MATH_INSTRUMENT_VARIABLE(rphi);
- T m = 2 * (phi - rphi) / constants::pi<T>();
+ T m = floor(2 * (phi - rphi) / constants::pi<T>() + 0.5f);
BOOST_MATH_INSTRUMENT_VARIABLE(m);
int s = 1;
if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5)
@@ -104,7 +104,7 @@
T cosp = cos(rphi);
BOOST_MATH_INSTRUMENT_VARIABLE(sinp);
BOOST_MATH_INSTRUMENT_VARIABLE(cosp);
- result = s * sinp * ellint_rf_imp(cosp * cosp, 1 - k * k * sinp * sinp, T(1), pol);
+ result = s * sinp * ellint_rf_imp(static_cast<T>(cosp * cosp), static_cast<T>(1 - k * k * sinp * sinp), T(1), pol);
BOOST_MATH_INSTRUMENT_VARIABLE(result);
if(m != 0)
{
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/ellint_2.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/ellint_2.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/ellint_2.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -74,8 +74,8 @@
// but that fails if T has more digits than a long long,
// so rewritten to use fmod instead:
//
- T rphi = boost::math::tools::fmod_workaround(phi, constants::pi<T>() / 2);
- T m = 2 * (phi - rphi) / constants::pi<T>();
+ T rphi = boost::math::tools::fmod_workaround(phi, static_cast<T>(constants::pi<T>() / 2));
+ T m = floor(2 * (phi - rphi) / constants::pi<T>() + 0.5f);
int s = 1;
if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5)
{
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/ellint_3.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/ellint_3.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/ellint_3.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -61,7 +61,7 @@
T sphi = sin(fabs(phi));
- if(v > 1 / (sphi * sphi))
+ if(tools::maybe_greater(v, static_cast<T>(1 / (sphi * sphi))))
{
// Complex result is a domain error:
return policies::raise_domain_error<T>(function,
@@ -102,7 +102,7 @@
// v > 1:
T vcr = sqrt(-vc);
T arg = vcr * tan(phi);
- return (boost::math::log1p(arg, pol) - boost::math::log1p(-arg, pol)) / (2 * vcr);
+ return (boost::math::log1p(arg, pol) - boost::math::log1p(static_cast<T>(-arg), pol)) / (2 * vcr);
}
}
@@ -185,8 +185,8 @@
}
else
{
- T rphi = boost::math::tools::fmod_workaround(fabs(phi), constants::pi<T>() / 2);
- T m = 2 * (fabs(phi) - rphi) / constants::pi<T>();
+ T rphi = boost::math::tools::fmod_workaround(fabs(phi), static_cast<T>(constants::pi<T>() / 2));
+ T m = floor(2 * (fabs(phi) - rphi) / constants::pi<T>() + 0.5f);
int sign = 1;
if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5)
{
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/ellint_rj.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/ellint_rj.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/ellint_rj.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -94,7 +94,7 @@
value = boost::math::ellint_rj(x, y, z, p, pol);
value *= pmy;
value -= 3 * boost::math::ellint_rf(x, y, z, pol);
- value += 3 * sqrt((x * y * z) / (x * z + p * q)) * boost::math::ellint_rc(x * z + p * q, p * q, pol);
+ value += 3 * sqrt((x * y * z) / (x * z + p * q)) * boost::math::ellint_rc(static_cast<T>(x * z + p * q), static_cast<T>(p * q), pol);
value /= (y + q);
return value;
}
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/erf.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/erf.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/erf.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -122,9 +122,9 @@
if(z < 0)
{
if(!invert)
- return -erf_imp(-z, invert, pol, t);
+ return -erf_imp(static_cast<T>(-z), invert, pol, t);
else
- return 1 + erf_imp(-z, false, pol, t);
+ return 1 + erf_imp(static_cast<T>(-z), false, pol, t);
}
T result;
@@ -142,7 +142,7 @@
if(x < 0.6)
{
// Compute P:
- result = z * exp(-x);
+ result = z * exp(static_cast<T>(-x));
result /= sqrt(boost::math::constants::pi<T>());
if(result != 0)
result *= 2 * detail::lower_gamma_series(T(0.5f), x, pol);
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/expint.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/expint.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/expint.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -436,7 +436,7 @@
if(z < 0)
return policies::raise_domain_error<T>(function, "Function requires z >= 0 but got %1%.", z, pol);
if(z == 0)
- return n == 1 ? policies::raise_overflow_error<T>(function, 0, pol) : 1 / (static_cast<T>(n - 1));
+ return n == 1 ? policies::raise_overflow_error<T>(function, 0, pol) : static_cast<T>(1 / (static_cast<T>(n - 1)));
T result;
@@ -504,7 +504,7 @@
{
static const char* function = "boost::math::expint<%1%>(%1%)";
if(z < 0)
- return -expint_imp(1, -z, pol, tag);
+ return -expint_imp(1, static_cast<T>(-z), pol, tag);
if(z == 0)
return -policies::raise_overflow_error<T>(function, 0, pol);
return expint_i_as_series(z, pol);
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/expm1.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/expm1.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/expm1.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -75,7 +75,7 @@
BOOST_MATH_STD_USING
T a = fabs(x);
- if(a > T(0.5L))
+ if(a > T(0.5f))
return exp(x) - T(1);
if(a < tools::epsilon<T>())
return x;
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/factorials.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/factorials.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/factorials.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -79,7 +79,7 @@
// Fallthrough: i is too large to use table lookup, try the
// gamma function instead.
//
- T result = boost::math::tgamma(static_cast<T>(i) / 2 + 1, pol) / sqrt(constants::pi<T>());
+ T result = boost::math::tgamma(static_cast<T>(static_cast<T>(i) / 2 + 1), pol) / sqrt(constants::pi<T>());
if(ldexp(tools::max_value<T>(), -static_cast<int>(i+1) / 2) > result)
return ceil(result * ldexp(T(1), (i+1) / 2) - 0.5f);
}
@@ -125,7 +125,7 @@
n = -n;
inv = true;
}
- T result = ((n&1) ? -1 : 1) * falling_factorial(-x, n, pol);
+ T result = ((n&1) ? -1 : 1) * falling_factorial(static_cast<T>(-x), n, pol);
if(inv)
result = 1 / result;
return result;
@@ -152,7 +152,7 @@
// For x < 0 we really have a rising factorial
// modulo a possible change of sign:
//
- return (n&1 ? -1 : 1) * rising_factorial(-x, n, pol);
+ return (n&1 ? -1 : 1) * rising_factorial(static_cast<T>(-x), n, pol);
}
if(n == 0)
return 1;
@@ -163,15 +163,15 @@
// handle it, split the product up into three parts:
//
T xp1 = x + 1;
- unsigned n2 = tools::real_cast<unsigned>(floor(xp1));
+ unsigned n2 = tools::real_cast<unsigned>(static_cast<T>(floor(xp1)));
if(n2 == xp1)
return 0;
- T result = boost::math::tgamma_delta_ratio(xp1, -static_cast<T>(n2), pol);
+ T result = boost::math::tgamma_delta_ratio(xp1, static_cast<T>(-static_cast<T>(n2)), pol);
x -= n2;
result *= x;
++n2;
if(n2 < n)
- result *= falling_factorial(x - 1, n - n2, pol);
+ result *= falling_factorial(static_cast<T>(x - 1), n - n2, pol);
return result;
}
//
@@ -181,7 +181,7 @@
// because tgamma_delta_ratio is alreay optimised
// for that use case:
//
- return boost::math::tgamma_delta_ratio(x + 1, -static_cast<T>(n), pol);
+ return boost::math::tgamma_delta_ratio(static_cast<T>(x + 1), static_cast<T>(-static_cast<T>(n)), pol);
}
} // namespace detail
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/gamma.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/gamma.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/gamma.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -146,7 +146,7 @@
return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
if(z <= -20)
{
- result = gamma_imp(-z, pol, l) * sinpx(z);
+ result = gamma_imp(static_cast<T>(-z), pol, l) * sinpx(z);
if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
result = -boost::math::constants::pi<T>() / result;
@@ -233,13 +233,17 @@
{
typedef typename policies::precision<T, Policy>::type precision_type;
typedef typename mpl::if_<
- mpl::less_equal<precision_type, mpl::int_<64> >,
- mpl::int_<64>,
- 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(), pol, l);
+ mpl::less_equal<precision_type, mpl::int_<0> >,
+ mpl::int_<0>,
+ mpl::if_<
+ mpl::less_equal<precision_type, mpl::int_<64> >,
+ mpl::int_<64>,
+ typename mpl::if_<
+ mpl::less_equal<precision_type, mpl::int_<113> >,
+ mpl::int_<113>, mpl::int_<0> >::type
+ >::type
+ >::type tag_type;
+ result = lgamma_small_imp(z, static_cast<T>(z - 1), static_cast<T>(z - 2), tag_type(), pol, l);
}
else if((z >= 3) && (z < 100))
{
@@ -346,7 +350,7 @@
return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
if(z <= -20)
{
- T result = gamma_imp(-z, pol, l) * sinpx(z);
+ T result = gamma_imp(static_cast<T>(-z), pol, l) * sinpx(z);
if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
result = -boost::math::constants::pi<T>() / result;
@@ -414,7 +418,8 @@
}
else if((z != 1) && (z != 2))
{
- T limit = (std::max)(z+1, T(10));
+ using std::max;
+ T limit = max BOOST_PREVENT_MACRO_SUBSTITUTION(static_cast<T>(z+1), T(10));
T prefix = z * log(limit) - limit;
T sum = detail::lower_gamma_series(z, limit, pol) / z;
sum += detail::upper_gamma_fraction(z, limit, ::boost::math::policies::digits<T, Policy>());
@@ -456,14 +461,14 @@
if(dz < -0.5)
{
// Best method is simply to subtract 1 from tgamma:
- result = boost::math::tgamma(1+dz, pol) - 1;
+ result = boost::math::tgamma(static_cast<T>(1+dz), pol) - 1;
BOOST_MATH_INSTRUMENT_CODE(result);
}
else
{
// Use expm1 on lgamma:
- result = boost::math::expm1(-boost::math::log1p(dz, pol)
- + lgamma_small_imp(dz+2, dz + 1, dz, tag_type(), pol, l));
+ result = boost::math::expm1(static_cast<T>(-boost::math::log1p(dz, pol)
+ + lgamma_small_imp(static_cast<T>(dz+2), static_cast<T>(dz + 1), dz, tag_type(), pol, l)));
BOOST_MATH_INSTRUMENT_CODE(result);
}
}
@@ -472,13 +477,13 @@
if(dz < 2)
{
// Use expm1 on lgamma:
- result = boost::math::expm1(lgamma_small_imp(dz+1, dz, dz-1, tag_type(), pol, l), pol);
+ result = boost::math::expm1(lgamma_small_imp(static_cast<T>(dz+1), dz, static_cast<T>(dz-1), tag_type(), pol, l), pol);
BOOST_MATH_INSTRUMENT_CODE(result);
}
else
{
// Best method is simply to subtract 1 from tgamma:
- result = boost::math::tgamma(1+dz, pol) - 1;
+ result = boost::math::tgamma(static_cast<T>(1+dz), pol) - 1;
BOOST_MATH_INSTRUMENT_CODE(result);
}
}
@@ -496,7 +501,7 @@
// algebra isn't easy for the general case....
// Start by subracting 1 from tgamma:
//
- T result = gamma_imp(1 + dz, pol, l) - 1;
+ T result = gamma_imp(static_cast<T>(1 + dz), pol, l) - 1;
BOOST_MATH_INSTRUMENT_CODE(result);
//
// Test the level of cancellation error observed: we loose one bit
@@ -597,6 +602,9 @@
T regularised_gamma_prefix(T a, T z, const Policy& pol, const L& l)
{
BOOST_MATH_STD_USING
+ using std::max;
+ using std::min;
+
T agh = a + static_cast<T>(L::g()) - T(0.5);
T prefix;
T d = ((z - a) - static_cast<T>(L::g()) + T(0.5)) / agh;
@@ -639,16 +647,16 @@
//
T alz = a * log(z / agh);
T amz = a - z;
- if(((std::min)(alz, amz) <= tools::log_min_value<T>()) || ((std::max)(alz, amz) >= tools::log_max_value<T>()))
+ if((min BOOST_PREVENT_MACRO_SUBSTITUTION(alz, amz) <= tools::log_min_value<T>()) || (max BOOST_PREVENT_MACRO_SUBSTITUTION(alz, amz) >= tools::log_max_value<T>()))
{
T amza = amz / a;
- if(((std::min)(alz, amz)/2 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/2 < tools::log_max_value<T>()))
+ if((min BOOST_PREVENT_MACRO_SUBSTITUTION(alz, amz)/2 > tools::log_min_value<T>()) && (max BOOST_PREVENT_MACRO_SUBSTITUTION(alz, amz)/2 < tools::log_max_value<T>()))
{
// compute square root of the result and then square it:
T sq = pow(z / agh, a / 2) * exp(amz / 2);
prefix = sq * sq;
}
- else if(((std::min)(alz, amz)/4 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/4 < tools::log_max_value<T>()) && (z > a))
+ else if((min BOOST_PREVENT_MACRO_SUBSTITUTION(alz, amz)/4 > tools::log_min_value<T>()) && (max BOOST_PREVENT_MACRO_SUBSTITUTION(alz, amz)/4 < tools::log_max_value<T>()) && (z > a))
{
// compute the 4th root of the result then square it twice:
T sq = pow(z / agh, a / 4) * exp(amz / 4);
@@ -679,8 +687,10 @@
T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined_lanczos&)
{
BOOST_MATH_STD_USING
+ using std::max;
+ using std::min;
- T limit = (std::max)(T(10), a);
+ T limit = max BOOST_PREVENT_MACRO_SUBSTITUTION(T(10), a);
T sum = detail::lower_gamma_series(a, limit, pol) / a;
sum += detail::upper_gamma_fraction(a, limit, ::boost::math::policies::digits<T, Policy>());
@@ -701,7 +711,7 @@
T amz = a - z;
T alzoa = a * log(zoa);
T prefix;
- if(((std::min)(alzoa, amz) <= tools::log_min_value<T>()) || ((std::max)(alzoa, amz) >= tools::log_max_value<T>()))
+ if((min BOOST_PREVENT_MACRO_SUBSTITUTION(alzoa, amz) <= tools::log_min_value<T>()) || (max BOOST_PREVENT_MACRO_SUBSTITUTION(alzoa, amz) >= tools::log_max_value<T>()))
{
T amza = amz / a;
if((amza <= tools::log_min_value<T>()) || (amza >= tools::log_max_value<T>()))
@@ -777,7 +787,7 @@
// Calculates normalised Q when a is a half-integer:
//
BOOST_MATH_STD_USING
- T e = boost::math::erfc(sqrt(x), pol);
+ T e = boost::math::erfc(static_cast<T>(sqrt(x)), pol);
if((e != 0) && (a > 1))
{
T term = exp(-x) / sqrt(constants::pi<T>() * x);
@@ -1025,7 +1035,7 @@
T result;
if(fabs(delta) < 10)
{
- result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));
+ result = exp((constants::half<T>() - z) * boost::math::log1p(static_cast<T>(delta / zgh), pol));
}
else
{
@@ -1095,7 +1105,7 @@
//
if((z <= max_factorial<T>::value) && (z + delta <= max_factorial<T>::value))
{
- return unchecked_factorial<T>(tools::real_cast<unsigned>(z) - 1) / unchecked_factorial<T>(tools::real_cast<unsigned>(z + delta) - 1);
+ return unchecked_factorial<T>(tools::real_cast<unsigned>(z) - 1) / unchecked_factorial<T>(tools::real_cast<unsigned>(static_cast<T>(z + delta)) - 1);
}
}
if(fabs(delta) < 20)
@@ -1426,7 +1436,7 @@
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
- return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(b) - static_cast<value_type>(a), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
+ return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(static_cast<value_type>(b) - static_cast<value_type>(a)), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
}
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/legendre.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/legendre.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/legendre.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -45,7 +45,7 @@
if(second)
{
// A solution of the second kind (Q):
- p0 = (boost::math::log1p(x, pol) - boost::math::log1p(-x, pol)) / 2;
+ p0 = (boost::math::log1p(x, pol) - boost::math::log1p(static_cast<T>(-x), pol)) / 2;
p1 = x * p0 - 1;
}
else
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/log1p.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/log1p.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/log1p.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -88,7 +88,7 @@
function, 0, pol);
result_type a = abs(result_type(x));
- if(a > result_type(0.5L))
+ if(a > result_type(0.5f))
return log(1 + result_type(x));
// Note that without numeric_limits specialisation support,
// epsilon just returns zero, and our "optimisation" will always fail:
@@ -248,7 +248,7 @@
function, 0, pol);
result_type a = abs(result_type(x));
- if(a > result_type(0.95L))
+ if(a > result_type(0.95f))
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:
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/sin_pi.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/sin_pi.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/sin_pi.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -22,7 +22,7 @@
{
BOOST_MATH_STD_USING // ADL of std names
if(x < 0)
- return -sin_pi(-x);
+ return -sin_pi(static_cast<T>(-x));
// sin of pi*x:
bool invert;
if(x < 0.5)
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/spherical_harmonic.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/spherical_harmonic.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/spherical_harmonic.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -61,7 +61,7 @@
if(m&1)
{
// Check phase if theta is outside [0, PI]:
- T mod = boost::math::tools::fmod_workaround(theta, 2 * constants::pi<T>());
+ T mod = boost::math::tools::fmod_workaround(theta, static_cast<T>(2 * constants::pi<T>()));
if(mod < 0)
mod += 2 * constants::pi<T>();
if(mod > constants::pi<T>())
@@ -88,7 +88,7 @@
if(m&1)
{
// Check phase if theta is outside [0, PI]:
- T mod = boost::math::tools::fmod_workaround(theta, 2 * constants::pi<T>());
+ T mod = boost::math::tools::fmod_workaround(theta, static_cast<T>(2 * constants::pi<T>()));
if(mod < 0)
mod += 2 * constants::pi<T>();
if(mod > constants::pi<T>())
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/sqrt1pm1.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/sqrt1pm1.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/sqrt1pm1.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -28,7 +28,7 @@
if(fabs(result_type(val)) > 0.75)
return sqrt(1 + result_type(val)) - 1;
- return boost::math::expm1(boost::math::log1p(val, pol) / 2, pol);
+ return boost::math::expm1(static_cast<T>(boost::math::log1p(val, pol) / 2), pol);
}
template <class T>
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/zeta.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/zeta.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/zeta.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -836,7 +836,7 @@
result = 0;
else
{
- result = boost::math::sin_pi(0.5f * sc, pol)
+ result = boost::math::sin_pi(static_cast<T>(0.5f * sc), pol)
* 2 * pow(2 * constants::pi<T>(), -s)
* boost::math::tgamma(s, pol)
* zeta_imp(s, sc, pol, tag);
@@ -884,7 +884,7 @@
return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::zeta_imp(
static_cast<value_type>(s),
- 1 - static_cast<value_type>(s),
+ static_cast<value_type>(1 - static_cast<value_type>(s)),
forwarding_policy(),
tag_type()), "boost::math::zeta<%1%>(%1%)");
}
Modified: sandbox/interval_math_toolkit/boost/math/tools/config.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/tools/config.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/tools/config.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -12,7 +12,9 @@
#include <boost/cstdint.hpp> // for boost::uintmax_t
#include <boost/config.hpp>
+#include <boost/mpl/bool.hpp>
#include <boost/detail/workaround.hpp>
+#include <boost/tr1/tuple.hpp>
#include <algorithm> // for min and max
#include <cmath>
#include <climits>
@@ -159,6 +161,63 @@
# define BOOST_MATH_CONTROL_FP
#endif
//
+// Forward declare interval type so we can handle it as a special
+// case where required.
+//
+namespace boost{ namespace numeric{
+
+template <class T, class Policy>
+class interval;
+
+template <class T>
+struct is_interval : public boost::mpl::false_{};
+
+template <class T, class Policy>
+struct is_interval<interval<T, Policy> > : public boost::mpl::true_{};
+
+template <class T>
+struct interval_base_type
+{
+ typedef T type;
+};
+
+template <class T, class Policy>
+struct interval_base_type<interval<T, Policy> >
+{
+ typedef T type;
+};
+
+// Namespace forward declarations:
+namespace interval_math_compare{}
+namespace interval_lib{ namespace compare{ namespace possible{} } }
+
+template <class T, class Policy>
+inline std::tr1::tuple<T,T> make_printable(const boost::numeric::interval<T, Policy>& v)
+{
+ return std::make_pair(v.lower(), v.upper());
+}
+
+template <class T>
+inline const T& make_printable(const T& v)
+{
+ return v;
+}
+
+#ifdef BOOST_MATH_INSTRUMENT
+//
+// We need a stream operator for interval types, this may well conflict
+// with a user provided definition:
+//
+template <class T, class Policy>
+std::ostream& operator << (std::ostream& os, const interval<T, Policy>& val)
+{
+ os << "{ " << val.lower() << ", " << val.upper() << " }";
+ return os;
+}
+#endif
+
+}} // namespaces
+//
// Helper macro for using statements:
//
#define BOOST_MATH_STD_USING \
@@ -184,7 +243,8 @@
using std::ceil;\
using std::floor;\
using std::log10;\
- using std::sqrt;
+ using std::sqrt;\
+ using namespace boost::numeric::interval_math_compare;
namespace boost{ namespace math{
@@ -194,14 +254,101 @@
template <class T>
inline T max BOOST_PREVENT_MACRO_SUBSTITUTION(T a, T b, T c)
{
- return (std::max)((std::max)(a, b), c);
+ using std::max;
+ return max BOOST_PREVENT_MACRO_SUBSTITUTION(max BOOST_PREVENT_MACRO_SUBSTITUTION(a, b), c);
}
template <class T>
inline T max BOOST_PREVENT_MACRO_SUBSTITUTION(T a, T b, T c, T d)
{
- return (std::max)((std::max)(a, b), (std::max)(c, d));
+ using std::max;
+ return max BOOST_PREVENT_MACRO_SUBSTITUTION(max BOOST_PREVENT_MACRO_SUBSTITUTION(a, b), max BOOST_PREVENT_MACRO_SUBSTITUTION(c, d));
+}
+//
+// Interval aware "soft" comparisons:
+//
+template <class T>
+inline bool maybe_equal(const T& a, const T& b)
+{
+ return a == b;
+}
+template <class T>
+inline bool maybe_less(const T& a, const T& b)
+{
+ return a < b;
+}
+template <class T>
+inline bool maybe_less_equal(const T& a, const T& b)
+{
+ return a <= b;
+}
+template <class T>
+inline bool maybe_greater(const T& a, const T& b)
+{
+ return a > b;
+}
+template <class T>
+inline bool maybe_greater_equal(const T& a, const T& b)
+{
+ return a >= b;
+}
+
+template <class T, class Policy>
+inline bool maybe_equal(const boost::numeric::interval<T, Policy>& a, const boost::numeric::interval<T, Policy>& b)
+{
+ using namespace boost::numeric::interval_lib::compare::possible;
+ return a == b;
+}
+template <class T, class Policy>
+inline bool maybe_less(const boost::numeric::interval<T, Policy>& a, const boost::numeric::interval<T, Policy>& b)
+{
+ using namespace boost::numeric::interval_lib::compare::possible;
+ return a < b;
+}
+template <class T, class Policy>
+inline bool maybe_less_equal(const boost::numeric::interval<T, Policy>& a, const boost::numeric::interval<T, Policy>& b)
+{
+ using namespace boost::numeric::interval_lib::compare::possible;
+ return a <= b;
}
+template <class T, class Policy>
+inline bool maybe_greater(const boost::numeric::interval<T, Policy>& a, const boost::numeric::interval<T, Policy>& b)
+{
+ using namespace boost::numeric::interval_lib::compare::possible;
+ return a > b;
+}
+template <class T, class Policy>
+inline bool maybe_greater_equal(const boost::numeric::interval<T, Policy>& a, const boost::numeric::interval<T, Policy>& b)
+{
+ using namespace boost::numeric::interval_lib::compare::possible;
+ return a >= b;
+}
+
+//
+// Normalisation functions:
+//
+template <class T>
+inline void normalize_median(T&){}
+template <class T, class Policy>
+inline void normalize_median(boost::numeric::interval<T, Policy>& v)
+{
+ v = median(v);
+}
+template <class T>
+inline void normalize_up(T&){}
+template <class T, class Policy>
+inline void normalize_up(boost::numeric::interval<T, Policy>& v)
+{
+ v = v.upper();
+}
+template <class T>
+inline void normalize_down(T&){}
+template <class T, class Policy>
+inline void normalize_down(boost::numeric::interval<T, Policy>& v)
+{
+ v = v.lower();
+}
+
} // namespace tools
}} // namespace boost namespace math
Modified: sandbox/interval_math_toolkit/boost/math/tools/fraction.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/tools/fraction.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/tools/fraction.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -116,7 +116,7 @@
D = 1/D;
delta = C*D;
f = f * delta;
- }while(fabs(delta - 1) > factor);
+ }while(!maybe_less_equal(fabs(delta - 1), factor));
return f;
}
@@ -155,7 +155,7 @@
D = 1/D;
delta = C*D;
f = f * delta;
- }while((fabs(delta - 1) > factor) && --counter);
+ }while(!maybe_less_equal(fabs(delta - 1), factor) && --counter);
max_terms = max_terms - counter;
@@ -209,7 +209,7 @@
D = 1/D;
delta = C*D;
f = f * delta;
- }while(fabs(delta - 1) > factor);
+ }while(!maybe_less_equal(fabs(delta - 1), factor));
return a0/f;
}
@@ -249,7 +249,7 @@
D = 1/D;
delta = C*D;
f = f * delta;
- }while((fabs(delta - 1) > factor) && --counter);
+ }while(!maybe_less_equal(fabs(delta - 1), factor) && --counter);
max_terms = max_terms - counter;
Modified: sandbox/interval_math_toolkit/boost/math/tools/minima.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/tools/minima.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/tools/minima.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -23,7 +23,9 @@
std::pair<T, T> brent_find_minima(F f, T min, T max, int bits, boost::uintmax_t& max_iter)
{
BOOST_MATH_STD_USING
- bits = (std::min)(policies::digits<T, policies::policy<> >() / 2, bits);
+ using std::min;
+
+ bits = min BOOST_PREVENT_MACRO_SUBSTITUTION(policies::digits<T, policies::policy<> >() / 2, bits);
T tolerance = static_cast<T>(ldexp(1.0, 1-bits));
T x; // minima so far
T w; // second best point
Modified: sandbox/interval_math_toolkit/boost/math/tools/rational.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/tools/rational.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/tools/rational.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -98,19 +98,19 @@
template <class T, class U>
inline U evaluate_even_polynomial(const T* poly, U z, std::size_t count)
{
- return evaluate_polynomial(poly, z*z, count);
+ return evaluate_polynomial(poly, static_cast<U>(z*z), count);
}
template <std::size_t N, class T, class V>
inline V evaluate_even_polynomial(const T(&a)[N], const V& z)
{
- return evaluate_polynomial(a, z*z);
+ return evaluate_polynomial(a, static_cast<V>(z*z));
}
template <std::size_t N, class T, class V>
inline V evaluate_even_polynomial(const boost::array<T,N>& a, const V& z)
{
- return evaluate_polynomial(a, z*z);
+ return evaluate_polynomial(a, static_cast<V>(z*z));
}
//
// Odd polynomials come next:
@@ -118,21 +118,21 @@
template <class T, class U>
inline U evaluate_odd_polynomial(const T* poly, U z, std::size_t count)
{
- return poly[0] + z * evaluate_polynomial(poly+1, z*z, count-1);
+ return poly[0] + z * evaluate_polynomial(poly+1, static_cast<U>(z*z), count-1);
}
template <std::size_t N, class T, class V>
inline V evaluate_odd_polynomial(const T(&a)[N], const V& z)
{
typedef mpl::int_<N-1> tag_type;
- return a[0] + z * detail::evaluate_polynomial_c_imp(static_cast<const T*>(a) + 1, z*z, static_cast<tag_type const*>(0));
+ return a[0] + z * detail::evaluate_polynomial_c_imp(static_cast<const T*>(a) + 1, static_cast<V>(z*z), static_cast<tag_type const*>(0));
}
template <std::size_t N, class T, class V>
inline V evaluate_odd_polynomial(const boost::array<T,N>& a, const V& z)
{
typedef mpl::int_<N-1> tag_type;
- return a[0] + z * detail::evaluate_polynomial_c_imp(static_cast<const T*>(a.data()) + 1, z*z, static_cast<tag_type const*>(0));
+ return a[0] + z * detail::evaluate_polynomial_c_imp(static_cast<const T*>(a.data()) + 1, static_cast<V>(z*z), static_cast<tag_type const*>(0));
}
template <class T, class U, class V>
Modified: sandbox/interval_math_toolkit/boost/math/tools/real_cast.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/tools/real_cast.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/tools/real_cast.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -10,6 +10,8 @@
#pragma once
#endif
+#include <boost/math/tools/config.hpp>
+
namespace boost{ namespace math
{
namespace tools
@@ -19,6 +21,11 @@
{
return static_cast<To>(t);
}
+ template <class To, class T, class Policy>
+ inline To real_cast(boost::numeric::interval<T, Policy> t)
+ {
+ return real_cast<To>(norm(t));
+ }
} // namespace tools
} // namespace math
} // namespace boost
Modified: sandbox/interval_math_toolkit/boost/math/tools/roots.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/tools/roots.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/tools/roots.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -44,20 +44,20 @@
T& delta,
T& result,
T& guess,
- const T& min,
- const T& max)
+ const T& vmin,
+ const T& vmax)
{
if(last_f0 == 0)
{
// this must be the first iteration, pretend that we had a
- // previous one at either min or max:
- if(result == min)
+ // previous one at either vmin or vmax:
+ if(result == vmin)
{
- guess = max;
+ guess = vmax;
}
else
{
- guess = min;
+ guess = vmin;
}
last_f0 = std::tr1::get<0>(f(guess));
delta = guess - result;
@@ -67,11 +67,11 @@
// we've crossed over so move in opposite direction to last step:
if(delta < 0)
{
- delta = (result - min) / 2;
+ delta = (result - vmin) / 2;
}
else
{
- delta = (result - max) / 2;
+ delta = (result - vmax) / 2;
}
}
else
@@ -79,11 +79,11 @@
// move in same direction as last step:
if(delta < 0)
{
- delta = (result - max) / 2;
+ delta = (result - vmax) / 2;
}
else
{
- delta = (result - min) / 2;
+ delta = (result - vmin) / 2;
}
}
}
@@ -91,28 +91,28 @@
} // namespace
template <class F, class T, class Tol, class Policy>
-std::pair<T, T> bisect(F f, T min, T max, Tol tol, boost::uintmax_t& max_iter, const Policy& pol)
+std::pair<T, T> bisect(F f, T vmin, T vmax, Tol tol, boost::uintmax_t& max_iter, const Policy& pol)
{
- T fmin = f(min);
- T fmax = f(max);
+ T fmin = f(vmin);
+ T fmax = f(vmax);
if(fmin == 0)
- return std::make_pair(min, min);
+ return std::make_pair(vmin, vmin);
if(fmax == 0)
- return std::make_pair(max, max);
+ return std::make_pair(vmax, vmax);
//
// Error checking:
//
static const char* function = "boost::math::tools::bisect<%1%>";
- if(min >= max)
+ if(vmin >= vmax)
{
policies::raise_evaluation_error(function,
- "Arguments in wrong order in boost::math::tools::bisect (first arg=%1%)", min, pol);
+ "Arguments in wrong order in boost::math::tools::bisect (first arg=%1%)", vmin, pol);
}
if(fmin * fmax >= 0)
{
policies::raise_evaluation_error(function,
- "No change of sign in boost::math::tools::bisect, either there is no root to find, or there are multiple roots in the interval (f(min) = %1%).", fmin, pol);
+ "No change of sign in boost::math::tools::bisect, either there is no root to find, or there are multiple roots in the interval (f(vmin) = %1%).", fmin, pol);
}
//
@@ -124,25 +124,25 @@
else
count -= 3;
- while(count && (0 == tol(min, max)))
+ while(count && (0 == tol(vmin, vmax)))
{
- T mid = (min + max) / 2;
+ T mid = (vmin + vmax) / 2;
T fmid = f(mid);
- if((mid == max) || (mid == min))
+ if((mid == vmax) || (mid == vmin))
break;
if(fmid == 0)
{
- min = max = mid;
+ vmin = vmax = mid;
break;
}
else if(sign(fmid) * sign(fmin) < 0)
{
- max = mid;
+ vmax = mid;
fmax = fmid;
}
else
{
- min = mid;
+ vmin = mid;
fmin = fmid;
}
--count;
@@ -161,24 +161,24 @@
}
#endif
- return std::make_pair(min, max);
+ return std::make_pair(vmin, vmax);
}
template <class F, class T, class Tol>
-inline std::pair<T, T> bisect(F f, T min, T max, Tol tol, boost::uintmax_t& max_iter)
+inline std::pair<T, T> bisect(F f, T vmin, T vmax, Tol tol, boost::uintmax_t& max_iter)
{
- return bisect(f, min, max, tol, max_iter, policies::policy<>());
+ return bisect(f, vmin, vmax, tol, max_iter, policies::policy<>());
}
template <class F, class T, class Tol>
-inline std::pair<T, T> bisect(F f, T min, T max, Tol tol)
+inline std::pair<T, T> bisect(F f, T vmin, T vmax, Tol tol)
{
boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
- return bisect(f, min, max, tol, m, policies::policy<>());
+ return bisect(f, vmin, vmax, tol, m, policies::policy<>());
}
template <class F, class T>
-T newton_raphson_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter)
+T newton_raphson_iterate(F f, T guess, T vmin, T vmax, int digits, boost::uintmax_t& max_iter)
{
BOOST_MATH_STD_USING
@@ -205,7 +205,7 @@
#ifdef BOOST_MATH_INSTRUMENT
std::cout << "Newton iteration, zero derivative found" << std::endl;
#endif
- detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
+ detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, vmin, vmax);
}
else
{
@@ -217,29 +217,29 @@
if(fabs(delta * 2) > fabs(delta2))
{
// last two steps haven't converged, try bisection:
- delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
+ delta = (delta > 0) ? (result - vmin) / 2 : (result - vmax) / 2;
}
guess = result;
result -= delta;
- if(result <= min)
+ if(result <= vmin)
{
- delta = 0.5F * (guess - min);
+ delta = 0.5F * (guess - vmin);
result = guess - delta;
- if((result == min) || (result == max))
+ if((result == vmin) || (result == vmax))
break;
}
- else if(result >= max)
+ else if(result >= vmax)
{
- delta = 0.5F * (guess - max);
+ delta = 0.5F * (guess - vmax);
result = guess - delta;
- if((result == min) || (result == max))
+ if((result == vmin) || (result == vmax))
break;
}
// update brackets:
if(delta > 0)
- max = guess;
+ vmax = guess;
else
- min = guess;
+ vmin = guess;
}while(--count && (fabs(result * factor) < fabs(delta)));
max_iter -= count;
@@ -259,22 +259,23 @@
}
template <class F, class T>
-inline T newton_raphson_iterate(F f, T guess, T min, T max, int digits)
+inline T newton_raphson_iterate(F f, T guess, T vmin, T vmax, int digits)
{
boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
- return newton_raphson_iterate(f, guess, min, max, digits, m);
+ return newton_raphson_iterate(f, guess, vmin, vmax, digits, m);
}
template <class F, class T>
-T halley_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter)
+T halley_iterate(F f, T guess, T vmin, T vmax, int digits, boost::uintmax_t& max_iter)
{
BOOST_MATH_STD_USING
+ using std::max;
T f0(0), f1, f2;
T result = guess;
T factor = static_cast<T>(ldexp(1.0, 1 - digits));
- T delta = (std::max)(10000000 * guess, T(10000000)); // arbitarily large delta
+ T delta = max BOOST_PREVENT_MACRO_SUBSTITUTION(static_cast<T>(10000000 * guess), T(10000000)); // arbitarily large delta
T last_f0 = 0;
T delta1 = delta;
T delta2 = delta;
@@ -300,7 +301,7 @@
#ifdef BOOST_MATH_INSTRUMENT
std::cout << "Halley iteration, zero derivative found" << std::endl;
#endif
- detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
+ detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, vmin, vmax);
}
else
{
@@ -331,7 +332,7 @@
if((convergence > 0.8) && (convergence < 2))
{
// last two steps haven't converged, try bisection:
- delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
+ delta = (delta > 0) ? (result - vmin) / 2 : (result - vmax) / 2;
// reset delta2 so that this branch will *not* be taken on the
// next iteration:
delta2 = delta * 3;
@@ -340,53 +341,53 @@
result -= delta;
// check for out of bounds step:
- if(result < min)
+ if(result < vmin)
{
- T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min))) ? 1000 : result / min;
+ T diff = ((fabs(vmin) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(vmin))) ? static_cast<T>(1000) : static_cast<T>(result / vmin);
if(fabs(diff) < 1)
diff = 1 / diff;
if(!out_of_bounds_sentry && (diff > 0) && (diff < 3))
{
// Only a small out of bounds step, lets assume that the result
- // is probably approximately at min:
- delta = 0.99f * (guess - min);
+ // is probably approximately at vmin:
+ delta = 0.99f * (guess - vmin);
result = guess - delta;
out_of_bounds_sentry = true; // only take this branch once!
}
else
{
- delta = (guess - min) / 2;
+ delta = (guess - vmin) / 2;
result = guess - delta;
- if((result == min) || (result == max))
+ if((result == vmin) || (result == vmax))
break;
}
}
- else if(result > max)
+ else if(result > vmax)
{
- T diff = ((fabs(max) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(max))) ? 1000 : result / max;
+ T diff = ((fabs(vmax) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(vmax))) ? static_cast<T>(1000) : static_cast<T>(result / vmax);
if(fabs(diff) < 1)
diff = 1 / diff;
if(!out_of_bounds_sentry && (diff > 0) && (diff < 3))
{
// Only a small out of bounds step, lets assume that the result
- // is probably approximately at min:
- delta = 0.99f * (guess - max);
+ // is probably approximately at vmin:
+ delta = 0.99f * (guess - vmax);
result = guess - delta;
out_of_bounds_sentry = true; // only take this branch once!
}
else
{
- delta = (guess - max) / 2;
+ delta = (guess - vmax) / 2;
result = guess - delta;
- if((result == min) || (result == max))
+ if((result == vmin) || (result == vmax))
break;
}
}
// update brackets:
if(delta > 0)
- max = guess;
+ vmax = guess;
else
- min = guess;
+ vmin = guess;
}while(--count && (fabs(result * factor) < fabs(delta)));
max_iter -= count;
@@ -406,14 +407,14 @@
}
template <class F, class T>
-inline T halley_iterate(F f, T guess, T min, T max, int digits)
+inline T halley_iterate(F f, T guess, T vmin, T vmax, int digits)
{
boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
- return halley_iterate(f, guess, min, max, digits, m);
+ return halley_iterate(f, guess, vmin, vmax, digits, m);
}
template <class F, class T>
-T schroeder_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter)
+T schroeder_iterate(F f, T guess, T vmin, T vmax, int digits, boost::uintmax_t& max_iter)
{
BOOST_MATH_STD_USING
@@ -444,7 +445,7 @@
#ifdef BOOST_MATH_INSTRUMENT
std::cout << "Halley iteration, zero derivative found" << std::endl;
#endif
- detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
+ detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, vmin, vmax);
}
else
{
@@ -462,32 +463,32 @@
if(fabs(delta * 2) > fabs(delta2))
{
// last two steps haven't converged, try bisection:
- delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
+ delta = (delta > 0) ? (result - vmin) / 2 : (result - vmax) / 2;
}
guess = result;
result -= delta;
#ifdef BOOST_MATH_INSTRUMENT
std::cout << "Halley iteration, delta = " << delta << std::endl;
#endif
- if(result <= min)
+ if(result <= vmin)
{
- delta = 0.5F * (guess - min);
+ delta = 0.5F * (guess - vmin);
result = guess - delta;
- if((result == min) || (result == max))
+ if((result == vmin) || (result == vmax))
break;
}
- else if(result >= max)
+ else if(result >= vmax)
{
- delta = 0.5F * (guess - max);
+ delta = 0.5F * (guess - vmax);
result = guess - delta;
- if((result == min) || (result == max))
+ if((result == vmin) || (result == vmax))
break;
}
// update brackets:
if(delta > 0)
- max = guess;
+ vmax = guess;
else
- min = guess;
+ vmin = guess;
}while(--count && (fabs(result * factor) < fabs(delta)));
max_iter -= count;
@@ -507,10 +508,10 @@
}
template <class F, class T>
-inline T schroeder_iterate(F f, T guess, T min, T max, int digits)
+inline T schroeder_iterate(F f, T guess, T vmin, T vmax, int digits)
{
boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
- return schroeder_iterate(f, guess, min, max, digits, m);
+ return schroeder_iterate(f, guess, vmin, vmax, digits, m);
}
} // namespace tools
Modified: sandbox/interval_math_toolkit/boost/math/tools/series.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/tools/series.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/tools/series.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -33,7 +33,7 @@
next_term = func();
result += next_term;
}
- while(fabs(result) < fabs(factor * next_term));
+ while(!maybe_greater_equal(fabs(result), fabs(factor * next_term)));
return result;
}
@@ -53,7 +53,7 @@
next_term = func();
result += next_term;
}
- while((fabs(result) < fabs(factor * next_term)) && --counter);
+ while(!maybe_greater_equal(fabs(result), fabs(factor * next_term)) && --counter);
// set max_terms to the actual number of terms of the series evaluated:
max_terms = max_terms - counter;
@@ -75,7 +75,7 @@
next_term = func();
result += next_term;
}
- while(fabs(result) < fabs(factor * next_term));
+ while(!maybe_greater_equal(fabs(result), fabs(factor * next_term)));
return result;
}
@@ -96,7 +96,7 @@
next_term = func();
result += next_term;
}
- while((fabs(result) < fabs(factor * next_term)) && --counter);
+ while(!maybe_greater_equal(fabs(result), fabs(factor * next_term)) && --counter);
// set max_terms to the actual number of terms of the series evaluated:
max_terms = max_terms - counter;
@@ -135,7 +135,7 @@
carry -= y;
result = t;
}
- while(fabs(result) < fabs(factor * next_term));
+ while(!maybe_greater_equal(fabs(result), fabs(factor * next_term)));
return result;
}
@@ -160,7 +160,7 @@
carry -= y;
result = t;
}
- while((fabs(result) < fabs(factor * next_term)) && --counter);
+ while(!maybe_greater_equal(fabs(result), fabs(factor * next_term)) && --counter);
// set max_terms to the actual number of terms of the series evaluated:
max_terms = max_terms - counter;
Modified: sandbox/interval_math_toolkit/boost/math/tools/test.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/tools/test.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/tools/test.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -14,6 +14,7 @@
#include <boost/math/tools/stats.hpp>
#include <boost/math/special_functions/fpclassify.hpp>
#include <boost/test/test_tools.hpp>
+#include <boost/mpl/bool.hpp>
#include <stdexcept>
namespace boost{ namespace math{ namespace tools{
@@ -158,7 +159,7 @@
{
if(i)
std::cout << ", ";
- std::cout << row[i];
+ std::cout << boost::numeric::make_printable(row[i]);
}
std::cout << std::endl;
}
@@ -171,7 +172,7 @@
// expressions.
//
template <class A, class F1, class F2>
-test_result<typename calculate_result_type<A>::value_type> test(const A& a, F1 test_func, F2 expect_func)
+test_result<typename calculate_result_type<A>::value_type> test(const A& a, F1 test_func, F2 expect_func, const mpl::false_&)
{
typedef typename A::value_type row_type;
typedef typename row_type::value_type value_type;
@@ -238,6 +239,116 @@
return result;
}
+template <class A, class F1, class F2>
+test_result<typename calculate_result_type<A>::value_type> test(const A& a, F1 test_func, F2 expect_func, const mpl::true_&)
+{
+ typedef typename A::value_type row_type;
+ typedef typename row_type::value_type value_type;
+
+ test_result<value_type> result;
+
+ for(unsigned i = 0; i < a.size(); ++i)
+ {
+ const row_type& row = a[i];
+ value_type point;
+ try
+ {
+ point = test_func(row);
+ }
+ catch(const std::underflow_error&)
+ {
+ point = 0;
+ }
+ catch(const std::overflow_error&)
+ {
+ point = std::numeric_limits<value_type>::has_infinity ?
+ std::numeric_limits<value_type>::infinity()
+ : tools::max_value<value_type>();
+ }
+ catch(const std::exception& e)
+ {
+ std::cerr << e.what() << std::endl;
+ print_row(row);
+ BOOST_ERROR("Unexpected exception.");
+ // so we don't get further errors:
+ point = expect_func(row);
+ }
+ value_type expected = expect_func(row);
+ value_type err = (median(point) != 0)
+ ? width(point) / median(point)
+ : width(point);
+
+#ifdef BOOST_INSTRUMENT
+ if(err != 0)
+ {
+ std::cout << std::setprecision(35);
+ std::cout << row[0] << " " << err;
+ std::cout << " (" << err.lower() / tools::epsilon<value_type>() << "eps)";
+ std::cout << std::endl;
+ }
+#endif
+ if(!(boost::math::isfinite)(point) && (boost::math::isfinite)(expected))
+ {
+ std::cout << std::setprecision(35);
+ std::cout << "CAUTION: Found non-finite result, when a finite value was expected at entry " << i << "\n";
+ std::cout << "Found: " << boost::numeric::make_printable(point) << " Expected " << boost::numeric::make_printable(expected) << " Error: " << boost::numeric::make_printable(err) << std::endl;
+ print_row(row);
+ BOOST_ERROR("Unexpected non-finite result");
+ }
+ if(err > 0.5)
+ {
+ std::cout << std::setprecision(35);
+ std::cout << "CAUTION: Gross error found at entry " << i << ".\n";
+ std::cout << "Found: " << boost::numeric::make_printable(point) << " Expected " << boost::numeric::make_printable(expected) << " Error: " << boost::numeric::make_printable(err) << std::endl;
+ print_row(row);
+ BOOST_ERROR("Gross error");
+ }
+ //
+ // Check whether the interval includes the result:
+ //
+ if(!overlap(point, expected))
+ {
+ //
+ // Since the transcendental functions do *not* guarentee
+ // to bracket the true value, we only raise an error if
+ // we're more than a few epsilon from the true value:
+ //
+ bool have_error = false;
+ int max_eps = 3;
+ typename value_type::base_type err = 0;
+ if(point.upper() < expected.lower())
+ err = (relative_error(point.upper(), expected.lower()) / epsilon<typename value_type::base_type>());
+ else if(point.lower() > expected.upper())
+ err = (relative_error(point.lower(), expected.upper()) / epsilon<typename value_type::base_type>());
+ have_error = max_eps < err;
+
+ if(have_error)
+ {
+ std::cout << std::setprecision(35);
+ std::cout << "Interval and expected result do not overlap.\n";
+ std::cout << "Found: " << boost::numeric::make_printable(point) << " Expected " << boost::numeric::make_printable(expected) << " Error: " << boost::numeric::make_printable(err) << std::endl;
+ print_row(row);
+ BOOST_ERROR("Non Overlapping Interval Error");
+ }
+ }
+
+ result.add(err);
+ if(maybe_greater_equal((result.max)(), err))
+ result.set_worst(i);
+ }
+ return result;
+}
+
+template <class A, class F1, class F2>
+inline test_result<typename calculate_result_type<A>::value_type> test(const A& a, F1 test_func, F2 expect_func)
+{
+ typedef typename A::value_type row_type;
+ typedef typename row_type::value_type value_type;
+ typedef typename boost::numeric::is_interval<value_type>::type tag_type;
+
+ return test(a, test_func, expect_func, tag_type());
+}
+
} // namespace tools
} // namespace math
} // namespace boost
Modified: sandbox/interval_math_toolkit/boost/math/tools/toms748_solve.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/tools/toms748_solve.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/tools/toms748_solve.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -26,12 +26,14 @@
eps_tolerance(unsigned bits)
{
BOOST_MATH_STD_USING
- eps = (std::max)(T(ldexp(1.0F, 1-bits)), 2 * tools::epsilon<T>());
+ using std::max;
+ eps = max BOOST_PREVENT_MACRO_SUBSTITUTION(T(ldexp(1.0F, 1-bits)), T(2 * tools::epsilon<T>()));
}
bool operator()(const T& a, const T& b)
{
BOOST_MATH_STD_USING
- return (fabs(a - b) / (std::min)(fabs(a), fabs(b))) <= eps;
+ using std::min;
+ return (fabs(a - b) / min BOOST_PREVENT_MACRO_SUBSTITUTION(fabs(a), fabs(b))) <= eps;
}
private:
T eps;
@@ -193,9 +195,9 @@
//
// Start by obtaining the coefficients of the quadratic polynomial:
//
- T B = safe_div(fb - fa, b - a, tools::max_value<T>());
- T A = safe_div(fd - fb, d - b, tools::max_value<T>());
- A = safe_div(A - B, d - a, T(0));
+ T B = safe_div(static_cast<T>(fb - fa), static_cast<T>(b - a), tools::max_value<T>());
+ T A = safe_div(static_cast<T>(fd - fb), static_cast<T>(d - b), tools::max_value<T>());
+ A = safe_div(static_cast<T>(A - B), static_cast<T>(d - a), T(0));
if(a == 0)
{
@@ -220,7 +222,7 @@
for(unsigned i = 1; i <= count; ++i)
{
//c -= safe_div(B * c, (B + A * (2 * c - a - b)), 1 + c - a);
- c -= safe_div(fa+(B+A*(c-b))*(c-a), (B + A * (2 * c - a - b)), 1 + c - a);
+ c -= safe_div(static_cast<T>(fa+(B+A*(c-b))*(c-a)), static_cast<T>((B + A * (2 * c - a - b))), static_cast<T>(1 + c - a));
}
if((c <= a) || (c >= b))
{
@@ -436,7 +438,7 @@
//
e = d;
fe = fd;
- detail::bracket(f, a, b, a + (b - a) / 2, fa, fb, d, fd);
+ detail::bracket(f, a, b, static_cast<T>(a + (b - a) / 2), fa, fb, d, fd);
--count;
BOOST_MATH_INSTRUMENT_CODE("Not converging: Taking a bisection!!!!");
BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
Modified: sandbox/interval_math_toolkit/libs/math/test/handle_test_result.hpp
==============================================================================
--- sandbox/interval_math_toolkit/libs/math/test/handle_test_result.hpp (original)
+++ sandbox/interval_math_toolkit/libs/math/test/handle_test_result.hpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -111,7 +111,7 @@
const char* test_name,
const char* group_name)
{
- using namespace std; // To aid selection of the right pow.
+ BOOST_MATH_STD_USING
T eps = boost::math::tools::epsilon<T>();
std::cout << std::setprecision(4);
@@ -120,8 +120,8 @@
//
// Begin by printing the main tag line with the results:
//
- std::cout << test_name << "<" << type_name << "> Max = " << max_error_found
- << " RMS Mean=" << mean_error_found;
+ std::cout << test_name << "<" << type_name << "> Max = " << boost::numeric::make_printable(max_error_found)
+ << " RMS Mean=" << boost::numeric::make_printable(mean_error_found);
//
// If the max error is non-zero, give the row of the table that
// produced the worst error:
@@ -137,7 +137,7 @@
#if defined(__SGI_STL_PORT)
std::cout << boost::math::tools::real_cast<double>(worst[i]);
#else
- std::cout << worst[i];
+ std::cout << boost::numeric::make_printable(worst[i]);
#endif
}
std::cout << " }";
@@ -146,16 +146,19 @@
//
// Now verify that the results are within our expected bounds:
//
- std::pair<boost::uintmax_t, boost::uintmax_t> const& bounds = get_max_errors(type_name, test_name, group_name);
- if(bounds.first < max_error_found)
+ if(!boost::numeric::is_interval<T>::value)
{
- std::cerr << "Peak error greater than expected value of " << bounds.first << std::endl;
- BOOST_CHECK(bounds.first >= max_error_found);
- }
- if(bounds.second < mean_error_found)
- {
- std::cerr << "Mean error greater than expected value of " << bounds.second << std::endl;
- BOOST_CHECK(bounds.second >= mean_error_found);
+ std::pair<boost::uintmax_t, boost::uintmax_t> const& bounds = get_max_errors(type_name, test_name, group_name);
+ if(bounds.first < max_error_found)
+ {
+ std::cerr << "Peak error greater than expected value of " << bounds.first << std::endl;
+ BOOST_CHECK(bounds.first >= max_error_found);
+ }
+ if(bounds.second < mean_error_found)
+ {
+ std::cerr << "Mean error greater than expected value of " << bounds.second << std::endl;
+ BOOST_CHECK(bounds.second >= mean_error_found);
+ }
}
std::cout << std::endl;
}
@@ -194,5 +197,23 @@
std::cout << std::endl;
}
+#ifdef TEST_INTERVAL
+
+#include <boost/math/bindings/interval.hpp>
+
+typedef boost::numeric::interval<double,
+ boost::numeric::interval_lib::policies<
+ boost::numeric::interval_lib::save_state<
+ boost::numeric::interval_lib::rounded_transc_nearest<
+ double,
+ boost::numeric::interval_lib::rounded_arith_opp<double>
+ >
+ >,
+ boost::numeric::interval_lib::checking_base<double>
+ >
+> interval_type;
+
+#endif
+
#endif // BOOST_MATH_HANDLE_TEST_RESULT
Added: sandbox/interval_math_toolkit/libs/math/test/interval_concept_check.cpp
==============================================================================
--- (empty file)
+++ sandbox/interval_math_toolkit/libs/math/test/interval_concept_check.cpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -0,0 +1,37 @@
+// Copyright John Maddock 2007-8.
+// 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)
+
+//
+// This tests two things: that boost::interval meets our
+// conceptual requirements, and that we can instantiate
+// all our distributions and special functions on this type.
+//
+#define BOOST_MATH_INSTRUMENT
+
+#define BOOST_MATH_ASSERT_UNDEFINED_POLICY false
+
+#include <boost/math/bindings/interval.hpp>
+#include <boost/math/concepts/real_type_concept.hpp>
+#include "compile_test/instantiate.hpp"
+
+
+typedef boost::numeric::interval<double,
+ boost::numeric::interval_lib::policies<
+ boost::numeric::interval_lib::save_state<
+ boost::numeric::interval_lib::rounded_transc_opp<double> >,
+ boost::numeric::interval_lib::checking_base<double> > > interval_type;
+
+
+void foo()
+{
+ instantiate(interval_type());
+}
+
+int main()
+{
+ BOOST_CONCEPT_ASSERT((boost::math::concepts::RealTypeConcept<double>));
+ BOOST_CONCEPT_ASSERT((boost::math::concepts::RealTypeConcept<interval_type>));
+}
+
Modified: sandbox/interval_math_toolkit/libs/math/test/log1p_expm1_test.cpp
==============================================================================
--- sandbox/interval_math_toolkit/libs/math/test/log1p_expm1_test.cpp (original)
+++ sandbox/interval_math_toolkit/libs/math/test/log1p_expm1_test.cpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -5,7 +5,12 @@
#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
+#ifdef TEST_INTERVAL
+#include <boost/math/bindings/interval.hpp>
+#else
#include <boost/math/concepts/real_concept.hpp>
+#endif
+
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/special_functions/log1p.hpp>
@@ -104,6 +109,7 @@
do_test(log1p_expm1_data, type_name, "expm1 and log1p");
+#ifndef TEST_INTERVAL
//
// C99 Appendix F special cases:
static const T zero = 0;
@@ -117,11 +123,13 @@
BOOST_CHECK_EQUAL(boost::math::expm1(-std::numeric_limits<T>::infinity()), m_one);
BOOST_CHECK_EQUAL(boost::math::expm1(std::numeric_limits<T>::infinity()), std::numeric_limits<T>::infinity());
}
+#endif
}
int test_main(int, char* [])
{
+#ifndef TEST_INTERVAL
expected_results();
BOOST_MATH_CONTROL_FP;
test(float(0), "float");
@@ -142,5 +150,8 @@
"not available at all, or because they are too inaccurate for these tests "
"to pass.</note>" << std::cout;
#endif
+#else // TEST_INTERVAL
+ test((interval_type)(0), "boost::interval<double>");
+#endif // TEST_INTERVAL
return 0;
}
Added: sandbox/interval_math_toolkit/libs/math/test/mpfr_concept_check.cpp
==============================================================================
--- (empty file)
+++ sandbox/interval_math_toolkit/libs/math/test/mpfr_concept_check.cpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -0,0 +1,30 @@
+// Copyright John Maddock 2007-8.
+// 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)
+
+//
+// This tests two things: that mpfr_class meets our
+// conceptual requirements, and that we can instantiate
+// all our distributions and special functions on this type.
+//
+#define BOOST_MATH_ASSERT_UNDEFINED_POLICY false
+#define TEST_MPFR
+
+#pragma warning(disable:4800)
+
+#include <boost/math/bindings/mpfr.hpp>
+#include <boost/math/concepts/real_type_concept.hpp>
+#include "compile_test/instantiate.hpp"
+
+void foo()
+{
+ //instantiate(mpfr_class());
+}
+
+int main()
+{
+ BOOST_CONCEPT_ASSERT((mpfr_class));
+ BOOST_CONCEPT_ASSERT((mpfr_class));
+}
+
Added: sandbox/interval_math_toolkit/libs/math/test/ntl_concept_check.cpp
==============================================================================
--- (empty file)
+++ sandbox/interval_math_toolkit/libs/math/test/ntl_concept_check.cpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -0,0 +1,27 @@
+// Copyright John Maddock 2007-8.
+// 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)
+
+//
+// This tests two things: that boost::math::ntl::RR meets our
+// conceptual requirements, and that we can instantiate
+// all our distributions and special functions on this type.
+//
+#define BOOST_MATH_ASSERT_UNDEFINED_POLICY false
+
+#include <boost/math/bindings/rr.hpp>
+#include <boost/math/concepts/real_type_concept.hpp>
+#include "compile_test/instantiate.hpp"
+
+void foo()
+{
+ instantiate(boost::math::ntl::RR());
+}
+
+int main()
+{
+ BOOST_CONCEPT_ASSERT((boost::math::ntl::RR));
+ BOOST_CONCEPT_ASSERT((boost::math::ntl::RR));
+}
+
Modified: sandbox/interval_math_toolkit/libs/math/test/powm1_sqrtp1m1_test.cpp
==============================================================================
--- sandbox/interval_math_toolkit/libs/math/test/powm1_sqrtp1m1_test.cpp (original)
+++ sandbox/interval_math_toolkit/libs/math/test/powm1_sqrtp1m1_test.cpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -3,7 +3,12 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#ifdef TEST_INTERVAL
+#include <boost/math/bindings/interval.hpp>
+#else
#include <boost/math/concepts/real_concept.hpp>
+#endif
+
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/special_functions/sqrt1pm1.hpp>
@@ -1637,6 +1642,7 @@
int test_main(int, char* [])
{
+#ifndef TEST_INTERVAL
expected_results();
BOOST_MATH_CONTROL_FP;
test_powm1_sqrtp1m1(1.0F, "float");
@@ -1652,6 +1658,9 @@
"not available at all, or because they are too inaccurate for these tests "
"to pass.</note>" << std::cout;
#endif
+#else // TEST_INTERVAL
+ test_powm1_sqrtp1m1(interval_type(), "boost::numeric::interval<double>");
+#endif // TEST_INTERVAL
return 0;
}
Modified: sandbox/interval_math_toolkit/libs/math/test/test_beta.cpp
==============================================================================
--- sandbox/interval_math_toolkit/libs/math/test/test_beta.cpp (original)
+++ sandbox/interval_math_toolkit/libs/math/test/test_beta.cpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -5,7 +5,12 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#ifdef TEST_INTERVAL
+#include <boost/math/bindings/interval.hpp>
+#else
#include <boost/math/concepts/real_concept.hpp>
+#endif
+
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/special_functions/beta.hpp>
@@ -186,6 +191,7 @@
int test_main(int, char* [])
{
+#ifndef TEST_INTERVAL
expected_results();
BOOST_MATH_CONTROL_FP;
test_spots(0.0F);
@@ -217,6 +223,9 @@
"not available at all, or because they are too inaccurate for these tests "
"to pass.</note>" << std::cout;
#endif
+#else // TEST_INTERVAL
+ test_beta(interval_type(), "boost::numeric::interval<double>");
+#endif // TEST_INTERVAL
return 0;
}
Modified: sandbox/interval_math_toolkit/libs/math/test/test_binomial_coeff.cpp
==============================================================================
--- sandbox/interval_math_toolkit/libs/math/test/test_binomial_coeff.cpp (original)
+++ sandbox/interval_math_toolkit/libs/math/test/test_binomial_coeff.cpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -3,7 +3,12 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#ifdef TEST_INTERVAL
+#include <boost/math/bindings/interval.hpp>
+#else
#include <boost/math/concepts/real_concept.hpp>
+#endif
+
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/special_functions/binomial.hpp>
@@ -163,6 +168,7 @@
int test_main(int, char* [])
{
+#ifndef TEST_INTERVAL
expected_results();
BOOST_MATH_CONTROL_FP;
@@ -186,6 +192,9 @@
"not available at all, or because they are too inaccurate for these tests "
"to pass.</note>" << std::cout;
#endif
+#else // TEST_INTERVAL
+ test_binomial(interval_type(), "boost:math::interval<double>");
+#endif // TEST_INTERVAL
return 0;
}
Modified: sandbox/interval_math_toolkit/libs/math/test/test_carlson.cpp
==============================================================================
--- sandbox/interval_math_toolkit/libs/math/test/test_carlson.cpp (original)
+++ sandbox/interval_math_toolkit/libs/math/test/test_carlson.cpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -5,7 +5,12 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#ifdef TEST_INTERVAL
+#include <boost/math/bindings/interval.hpp>
+#else
#include <boost/math/concepts/real_concept.hpp>
+#endif
+
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/special_functions/ellint_rf.hpp>
@@ -195,7 +200,7 @@
std::cout << "Testing: " << test << std::endl;
- value_type (*fp)(value_type, value_type, value_type) = boost::math::ellint_rd;
+ value_type (*fp)(value_type, value_type, value_type) = boost::math::ellint_rd;
boost::math::tools::test_result<value_type> result;
result = boost::math::tools::test(
@@ -214,6 +219,8 @@
{
using namespace boost::math;
using namespace std;
+
+#ifndef TEST_INTERVAL
// Spot values from Numerical Computation of Real or Complex
// Elliptic Integrals, B. C. Carlson: http://arxiv.org/abs/math.CA/9409227
// RF:
@@ -281,7 +288,7 @@
s2 = ellint_rj(x, y, z, z);
BOOST_CHECK_CLOSE_FRACTION(s1, s2, eps40);
}
-
+#endif // TEST_INTERVAL
//
// Now random spot values:
//
@@ -304,6 +311,7 @@
int test_main(int, char* [])
{
+#ifndef TEST_INTERVAL
expected_results();
BOOST_MATH_CONTROL_FP;
@@ -322,7 +330,9 @@
"not available at all, or because they are too inaccurate for these tests "
"to pass.</note>" << std::cout;
#endif
-
+#else // TEST_INTERVAL
+ test_spots(interval_type(0), "boost::numeric::interval<double>");
+#endif // TEST_INTERVAL
return 0;
}
Modified: sandbox/interval_math_toolkit/libs/math/test/test_cbrt.cpp
==============================================================================
--- sandbox/interval_math_toolkit/libs/math/test/test_cbrt.cpp (original)
+++ sandbox/interval_math_toolkit/libs/math/test/test_cbrt.cpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -3,7 +3,12 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#ifdef TEST_INTERVAL
+#include <boost/math/bindings/interval.hpp>
+#else
#include <boost/math/concepts/real_concept.hpp>
+#endif
+
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/tools/stats.hpp>
@@ -115,6 +120,7 @@
int test_main(int, char* [])
{
+#ifndef TEST_INTERVAL
expected_results();
BOOST_MATH_CONTROL_FP;
test_cbrt(0.1F, "float");
@@ -125,6 +131,9 @@
test_cbrt(boost::math::concepts::real_concept(0.1), "real_concept");
#endif
#endif
+#else // TEST_INTERVAL
+ test_cbrt(interval_type(0.1), "boost::numeric::interval<double>");
+#endif // TEST_INTERVAL
return 0;
}
Modified: sandbox/interval_math_toolkit/libs/math/test/test_digamma.cpp
==============================================================================
--- sandbox/interval_math_toolkit/libs/math/test/test_digamma.cpp (original)
+++ sandbox/interval_math_toolkit/libs/math/test/test_digamma.cpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -3,7 +3,12 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#ifdef TEST_INTERVAL
+#include <boost/math/bindings/interval.hpp>
+#else
#include <boost/math/concepts/real_concept.hpp>
+#endif
+
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/special_functions/digamma.hpp>
@@ -149,6 +154,7 @@
int test_main(int, char* [])
{
+#ifndef TEST_INTERVAL
BOOST_MATH_CONTROL_FP;
test_spots(0.0F, "float");
test_spots(0.0, "double");
@@ -172,6 +178,9 @@
"not available at all, or because they are too inaccurate for these tests "
"to pass.</note>" << std::cout;
#endif
+#else // TEST_INTERVAL
+ test_digamma(interval_type(0.1), "boost::numeric::interval<double>");
+#endif
return 0;
}
Modified: sandbox/interval_math_toolkit/libs/math/test/test_ellint_1.cpp
==============================================================================
--- sandbox/interval_math_toolkit/libs/math/test/test_ellint_1.cpp (original)
+++ sandbox/interval_math_toolkit/libs/math/test/test_ellint_1.cpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -11,7 +11,12 @@
// Constants are too big for float case, but this doesn't matter for test.
#endif
+#ifdef TEST_INTERVAL
+#include <boost/math/bindings/interval.hpp>
+#else
#include <boost/math/concepts/real_concept.hpp>
+#endif
+
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/special_functions/ellint_1.hpp>
@@ -192,6 +197,8 @@
int test_main(int, char* [])
{
expected_results();
+
+#ifndef TEST_INTERVAL
BOOST_MATH_CONTROL_FP;
test_spots(0.0F, "float");
@@ -207,6 +214,8 @@
"not available at all, or because they are too inaccurate for these tests "
"to pass.</note>" << std::cout;
#endif
-
+#else // TEST_INTERVAL
+ test_spots(interval_type(0), "boost::numeric::interval<T>");
+#endif // TEST_INTERVAL
return 0;
}
Modified: sandbox/interval_math_toolkit/libs/math/test/test_ellint_2.cpp
==============================================================================
--- sandbox/interval_math_toolkit/libs/math/test/test_ellint_2.cpp (original)
+++ sandbox/interval_math_toolkit/libs/math/test/test_ellint_2.cpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -11,7 +11,12 @@
// Constants are too big for float case, but this doesn't matter for test.
#endif
+#ifdef TEST_INTERVAL
+#include <boost/math/bindings/interval.hpp>
+#else
#include <boost/math/concepts/real_concept.hpp>
+#endif
+
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/special_functions/ellint_2.hpp>
@@ -183,6 +188,7 @@
int test_main(int, char* [])
{
expected_results();
+#ifndef TEST_INTERVAL
BOOST_MATH_CONTROL_FP;
#ifndef BOOST_MATH_BUGGY_LARGE_FLOAT_CONSTANTS
test_spots(0.0F, "float");
@@ -199,6 +205,8 @@
"not available at all, or because they are too inaccurate for these tests "
"to pass.</note>" << std::cout;
#endif
-
+#else // TEST_INTERVAL
+ test_spots(interval_type(0), "boost::numeric::interval<double>");
+#endif // TEST_INTERVAL
return 0;
}
Modified: sandbox/interval_math_toolkit/libs/math/test/test_ellint_3.cpp
==============================================================================
--- sandbox/interval_math_toolkit/libs/math/test/test_ellint_3.cpp (original)
+++ sandbox/interval_math_toolkit/libs/math/test/test_ellint_3.cpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -11,7 +11,12 @@
// Constants are too big for float case, but this doesn't matter for test.
#endif
+#ifdef TEST_INTERVAL
+#include <boost/math/bindings/interval.hpp>
+#else
#include <boost/math/concepts/real_concept.hpp>
+#endif
+
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/special_functions/ellint_3.hpp>
@@ -220,6 +225,7 @@
int test_main(int, char* [])
{
expected_results();
+#ifndef TEST_INTERVAL
BOOST_MATH_CONTROL_FP;
test_spots(0.0F, "float");
test_spots(0.0, "double");
@@ -234,6 +240,8 @@
"not available at all, or because they are too inaccurate for these tests "
"to pass.</note>" << std::cout;
#endif
-
+#else // TEST_INTERVAL
+ test_spots(interval_type(0.0), "boost::numeric::interval<double>");
+#endif // TEST_INTERVAL
return 0;
}
Modified: sandbox/interval_math_toolkit/libs/math/test/test_erf.cpp
==============================================================================
--- sandbox/interval_math_toolkit/libs/math/test/test_erf.cpp (original)
+++ sandbox/interval_math_toolkit/libs/math/test/test_erf.cpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -4,7 +4,12 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#ifdef TEST_INTERVAL
+#include <boost/math/bindings/interval.hpp>
+#else
#include <boost/math/concepts/real_concept.hpp>
+#endif
+
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/special_functions/erf.hpp>
@@ -317,6 +322,7 @@
int test_main(int, char* [])
{
+#ifndef TEST_INTERVAL
BOOST_MATH_CONTROL_FP;
test_spots(0.0F, "float");
test_spots(0.0, "double");
@@ -344,6 +350,9 @@
"not available at all, or because they are too inaccurate for these tests "
"to pass.</note>" << std::cout;
#endif
+#else // TEST_INTERVAL
+ test_erf(interval_type(0), "boost::numeric::interval<double>");
+#endif // TEST_INTERVAL
return 0;
}
Modified: sandbox/interval_math_toolkit/libs/math/test/test_expint.cpp
==============================================================================
--- sandbox/interval_math_toolkit/libs/math/test/test_expint.cpp (original)
+++ sandbox/interval_math_toolkit/libs/math/test/test_expint.cpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -3,7 +3,12 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#ifdef TEST_INTERVAL
+#include <boost/math/bindings/interval.hpp>
+#else
#include <boost/math/concepts/real_concept.hpp>
+#endif
+
#include <boost/math/special_functions/expint.hpp>
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
@@ -293,6 +298,7 @@
int test_main(int, char* [])
{
expected_results();
+#ifndef TEST_INTERVAL
BOOST_MATH_CONTROL_FP;
boost::math::expint(114.7);
@@ -320,6 +326,9 @@
"not available at all, or because they are too inaccurate for these tests "
"to pass.</note>" << std::cout;
#endif
+#else // TEST_INTERVAL
+ test_expint(interval_type(), "boost::numeric::interval<double>");
+#endif // TEST_INTERVAL
return 0;
}
Modified: sandbox/interval_math_toolkit/libs/math/test/test_gamma.cpp
==============================================================================
--- sandbox/interval_math_toolkit/libs/math/test/test_gamma.cpp (original)
+++ sandbox/interval_math_toolkit/libs/math/test/test_gamma.cpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -5,7 +5,12 @@
#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
+#ifdef TEST_INTERVAL
+#include <boost/math/bindings/interval.hpp>
+#else
#include <boost/math/concepts/real_concept.hpp>
+#endif
+
#include <boost/math/special_functions/gamma.hpp>
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
@@ -450,6 +455,7 @@
int test_main(int, char* [])
{
expected_results();
+#ifndef TEST_INTERVAL
BOOST_MATH_CONTROL_FP;
#ifndef BOOST_MATH_BUGGY_LARGE_FLOAT_CONSTANTS
@@ -478,6 +484,9 @@
"not available at all, or because they are too inaccurate for these tests "
"to pass.</note>" << std::cout;
#endif
+#else // TEST_INTERVAL
+ test_gamma(interval_type(), "boost::numeric::interval<double>");
+#endif // TEST_INTERVAL
return 0;
}
Modified: sandbox/interval_math_toolkit/libs/math/test/test_hermite.cpp
==============================================================================
--- sandbox/interval_math_toolkit/libs/math/test/test_hermite.cpp (original)
+++ sandbox/interval_math_toolkit/libs/math/test/test_hermite.cpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -12,7 +12,12 @@
// Constants are too big for float case, but this doesn't matter for test.
#endif
+#ifdef TEST_INTERVAL
+#include <boost/math/bindings/interval.hpp>
+#else
#include <boost/math/concepts/real_concept.hpp>
+#endif
+
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/special_functions/hermite.hpp>
@@ -168,8 +173,7 @@
{
BOOST_MATH_CONTROL_FP;
- boost::math::hermite(51, 915.0);
-
+#ifndef TEST_INTERVAL
#ifndef BOOST_MATH_BUGGY_LARGE_FLOAT_CONSTANTS
test_spots(0.0F, "float");
#endif
@@ -196,6 +200,9 @@
"not available at all, or because they are too inaccurate for these tests "
"to pass.</note>" << std::cout;
#endif
+#else // TEST_INTERVAL
+ test_hermite(interval_type(), "boost::numeric::interval<double>");
+#endif // TEST_INTERVAL
return 0;
}
Modified: sandbox/interval_math_toolkit/libs/math/test/test_ibeta.cpp
==============================================================================
--- sandbox/interval_math_toolkit/libs/math/test/test_ibeta.cpp (original)
+++ sandbox/interval_math_toolkit/libs/math/test/test_ibeta.cpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -3,7 +3,12 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#ifdef TEST_INTERVAL
+#include <boost/math/bindings/interval.hpp>
+#else
#include <boost/math/concepts/real_concept.hpp>
+#endif
+
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/special_functions/beta.hpp>
@@ -515,6 +520,7 @@
int test_main(int, char* [])
{
expected_results();
+#ifndef TEST_INTERVAL
BOOST_MATH_CONTROL_FP;
#ifdef TEST_GSL
gsl_set_error_handler_off();
@@ -557,6 +563,9 @@
"not available at all, or because they are too inaccurate for these tests "
"to pass.</note>" << std::cout;
#endif
+#else // TEST_INTERVAL
+ test_beta(interval_type(), "boost::numeric::interval<double>");
+#endif // TEST_INTERVAL
return 0;
}
Modified: sandbox/interval_math_toolkit/libs/math/test/test_ibeta_inv.cpp
==============================================================================
--- sandbox/interval_math_toolkit/libs/math/test/test_ibeta_inv.cpp (original)
+++ sandbox/interval_math_toolkit/libs/math/test/test_ibeta_inv.cpp 2008-01-05 05:09:55 EST (Sat, 05 Jan 2008)
@@ -3,7 +3,12 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#ifdef TEST_INTERVAL
+#include <boost/math/bindings/interval.hpp>
+#else
#include <boost/math/concepts/real_concept.hpp>
+#endif
+
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/special_functions/beta.hpp>
@@ -270,6 +275,7 @@
// The contents are as follows, each row of data contains
// five items, input value a, input value b, integration limits x, beta(a, b, x) and ibeta(a, b, x):
//
+#ifndef TEST_INTERVAL
# include "ibeta_small_data.ipp"
test_inverses(ibeta_small_data);
@@ -282,6 +288,8 @@
test_inverses(ibeta_large_data);
+#endif // TEST_INTERVAL
+
# include "ibeta_inv_data.ipp"
test_inverses2(ibeta_inv_data, name, "Inverse incomplete beta");
@@ -324,6 +332,7 @@
{
BOOST_MATH_CONTROL_FP;
expected_results();
+#ifndef TEST_INTERVAL
#ifdef TEST_GSL
gsl_set_error_handler_off();
#endif
@@ -363,6 +372,9 @@
"not available at all, or because they are too inaccurate for these tests "
"to pass.</note>" << std::cout;
#endif
+#else // TEST_INTERVAL
+ test_beta(interval_type(), "boost::numeric::interval<double>");
+#endif // TEST_INTERVAL
return 0;
}