|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r56598 - in trunk/boost/math: policies special_functions special_functions/detail tools
From: john_at_[hidden]
Date: 2009-10-05 13:36:59
Author: johnmaddock
Date: 2009-10-05 13:36:56 EDT (Mon, 05 Oct 2009)
New Revision: 56598
URL: http://svn.boost.org/trac/boost/changeset/56598
Log:
Add a check on iterations when Halley iterating.
Correct igamma inverse initial calculations, and adjust digits needed for convergence.
Fix erf calculation limits.
Adjust non-convergence behaviour in Halley iteration.
For issue #3408.
Text files modified:
trunk/boost/math/policies/error_handling.hpp | 9 ++++
trunk/boost/math/special_functions/cbrt.hpp | 7 ++
trunk/boost/math/special_functions/detail/erf_inv.hpp | 6 +
trunk/boost/math/special_functions/detail/ibeta_inverse.hpp | 4 +
trunk/boost/math/special_functions/detail/igamma_inverse.hpp | 84 +++++++++++++++++++++++++++++++--------
trunk/boost/math/special_functions/erf.hpp | 24 ++++++----
trunk/boost/math/tools/roots.hpp | 2
7 files changed, 103 insertions(+), 33 deletions(-)
Modified: trunk/boost/math/policies/error_handling.hpp
==============================================================================
--- trunk/boost/math/policies/error_handling.hpp (original)
+++ trunk/boost/math/policies/error_handling.hpp 2009-10-05 13:36:56 EDT (Mon, 05 Oct 2009)
@@ -632,6 +632,15 @@
"Series evaluation exceeded %1% iterations, giving up now.", max_iter, pol);
}
+template <class Policy>
+inline void check_root_iterations(const char* function, boost::uintmax_t max_iter, const Policy& pol)
+{
+ if(max_iter >= policies::get_max_root_iterations<Policy>())
+ raise_evaluation_error<boost::uintmax_t>(
+ function,
+ "Root finding evaluation exceeded %1% iterations, giving up now.", max_iter, pol);
+}
+
} //namespace policies
#ifdef BOOST_MSVC
Modified: trunk/boost/math/special_functions/cbrt.hpp
==============================================================================
--- trunk/boost/math/special_functions/cbrt.hpp (original)
+++ trunk/boost/math/special_functions/cbrt.hpp 2009-10-05 13:36:56 EDT (Mon, 05 Oct 2009)
@@ -32,7 +32,7 @@
};
template <class T, class Policy>
-T cbrt_imp(T z, const Policy&)
+T cbrt_imp(T z, const Policy& pol)
{
BOOST_MATH_STD_USING
int i_exp, sign(1);
@@ -49,7 +49,10 @@
T max = static_cast<T>(ldexp(2.0, i_exp/3));
T guess = static_cast<T>(ldexp(1.0, i_exp/3));
int digits = (policies::digits<T, Policy>()) / 2;
- return sign * tools::halley_iterate(detail::cbrt_functor<T>(z), guess, min, max, digits);
+ boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
+ guess = sign * tools::halley_iterate(detail::cbrt_functor<T>(z), guess, min, max, digits, max_iter);
+ policies::check_root_iterations("boost::math::cbrt<%1%>", max_iter, pol);
+ return guess;
}
} // namespace detail
Modified: trunk/boost/math/special_functions/detail/erf_inv.hpp
==============================================================================
--- trunk/boost/math/special_functions/detail/erf_inv.hpp (original)
+++ trunk/boost/math/special_functions/detail/erf_inv.hpp 2009-10-05 13:36:56 EDT (Mon, 05 Oct 2009)
@@ -304,14 +304,16 @@
//
if(policies::digits<T, Policy>() > 64)
{
+ boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
if(p <= 0.5)
{
- result = tools::halley_iterate(detail::erf_roots<typename remove_cv<T>::type, Policy>(p, 1), guess, static_cast<T>(0), tools::max_value<T>(), (policies::digits<T, Policy>() * 2) / 3);
+ result = tools::halley_iterate(detail::erf_roots<typename remove_cv<T>::type, Policy>(p, 1), guess, static_cast<T>(0), tools::max_value<T>(), (policies::digits<T, Policy>() * 2) / 3, max_iter);
}
else
{
- result = tools::halley_iterate(detail::erf_roots<typename remove_cv<T>::type, Policy>(q, -1), guess, static_cast<T>(0), tools::max_value<T>(), (policies::digits<T, Policy>() * 2) / 3);
+ result = tools::halley_iterate(detail::erf_roots<typename remove_cv<T>::type, Policy>(q, -1), guess, static_cast<T>(0), tools::max_value<T>(), (policies::digits<T, Policy>() * 2) / 3, max_iter);
}
+ policies::check_root_iterations("boost::math::erf_inv<%1%>", max_iter, pol);
}
else
{
Modified: trunk/boost/math/special_functions/detail/ibeta_inverse.hpp
==============================================================================
--- trunk/boost/math/special_functions/detail/ibeta_inverse.hpp (original)
+++ trunk/boost/math/special_functions/detail/ibeta_inverse.hpp 2009-10-05 13:36:56 EDT (Mon, 05 Oct 2009)
@@ -795,8 +795,10 @@
// Now iterate, we can use either p or q as the target here
// depending on which is smaller:
//
+ boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
x = boost::math::tools::halley_iterate(
- boost::math::detail::ibeta_roots<T, Policy>(a, b, (p < q ? p : q), (p < q ? false : true)), x, lower, upper, digits);
+ boost::math::detail::ibeta_roots<T, Policy>(a, b, (p < q ? p : q), (p < q ? false : true)), x, lower, upper, digits, max_iter);
+ policies::check_root_iterations("boost::math::ibeta<%1%>(%1%, %1%, %1%)", max_iter, pol);
//
// We don't really want these asserts here, but they are useful for sanity
// checking that we have the limits right, uncomment if you suspect bugs *only*.
Modified: trunk/boost/math/special_functions/detail/igamma_inverse.hpp
==============================================================================
--- trunk/boost/math/special_functions/detail/igamma_inverse.hpp (original)
+++ trunk/boost/math/special_functions/detail/igamma_inverse.hpp 2009-10-05 13:36:56 EDT (Mon, 05 Oct 2009)
@@ -93,7 +93,7 @@
}
template <class T, class Policy>
-T find_inverse_gamma(T a, T p, T q, const Policy& pol)
+T find_inverse_gamma(T a, T p, T q, const Policy& pol, bool* p_has_10_digits)
{
//
// In order to understand what's going on here, you will
@@ -107,6 +107,7 @@
BOOST_MATH_STD_USING
T result;
+ *p_has_10_digits = false;
if(a == 1)
{
@@ -191,6 +192,8 @@
T y_4 = y_2 * y_2;
result = y + c1 + (c2 / y) + (c3 / y_2) + (c4 / y_3) + (c5 / y_4);
BOOST_MATH_INSTRUMENT_VARIABLE(result);
+ if(b < 1e-28f)
+ *p_has_10_digits = true;
}
}
else
@@ -218,6 +221,7 @@
if((a >= 500) && (fabs(1 - w / a) < 1e-6))
{
result = w;
+ *p_has_10_digits = true;
BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
else if (p > 0.5)
@@ -269,25 +273,37 @@
}
else
{
- // DiDonato and Morris Eq 35:
- T z = didonato_FN(p, a, w, 0, T(0), pol);
- z = didonato_FN(p, a, z, 2, T(0), pol);
- z = didonato_FN(p, a, z, 2, T(0), pol);
- z = didonato_FN(p, a, z, 3, T(0), pol);
-
- BOOST_MATH_INSTRUMENT_VARIABLE(z);
+ T z = w;
+ T ap1 = a + 1;
+ if(w < 0.15f * ap1)
+ {
+ // DiDonato and Morris Eq 35:
+ T v = log(p) + boost::math::lgamma(ap1, pol);
+ T s = 1;
+ z = exp((v + w) / a);
+ s = boost::math::log1p(z / ap1 * (1 + z / (a + 2)));
+ z = exp((v + z - s) / a);
+ z = exp((v + z - s) / a);
+ s = boost::math::log1p(z / ap1 * (1 + z / (a + 2) * (1 + z / (a + 3))));
+ z = exp((v + z - s) / a);
+ BOOST_MATH_INSTRUMENT_VARIABLE(z);
+ }
- if((z <= 0.01 * (a + 1)) || (z > 0.7 * (a + 1)))
+ if((z <= 0.01 * ap1) || (z > 0.7 * ap1))
{
result = z;
+ if(z <= 0.002 * ap1)
+ *p_has_10_digits = true;
BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
else
{
// 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);
- result = zb * (1 - (a * log(zb) - zb - u + log(didonato_SN(a, z, 100, T(1e-4)))) / (a - zb));
+ T ls = log(didonato_SN(a, z, 100, T(1e-4)));
+ T v = log(p) + boost::math::lgamma(ap1, pol);
+ z = exp((v + z - ls) / a);
+ result = z * (1 - (a * log(z) - z - v + ls) / (a - z));
+
BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
}
@@ -386,7 +402,10 @@
return tools::max_value<T>();
if(p == 0)
return 0;
- T guess = detail::find_inverse_gamma<T>(a, p, 1 - p, pol);
+ bool has_10_digits;
+ T guess = detail::find_inverse_gamma<T>(a, p, 1 - p, pol, &has_10_digits);
+ if((policies::digits<T, Policy>() <= 36) && has_10_digits)
+ return guess;
T lower = tools::min_value<T>();
if(guess <= lower)
guess = tools::min_value<T>();
@@ -397,18 +416,31 @@
// large convergence is slow, so we'll bump it up to full
// precision to prevent premature termination of the root-finding routine.
//
- unsigned digits = (policies::digits<T, Policy>() * 2) / 3;
+ unsigned digits = policies::digits<T, Policy>();
+ if(digits < 30)
+ {
+ digits *= 2;
+ digits /= 3;
+ }
+ else
+ {
+ digits /= 2;
+ digits -= 1;
+ }
if((a < 0.125) && (fabs(gamma_p_derivative(a, guess, pol)) > 1 / sqrt(tools::epsilon<T>())))
digits = policies::digits<T, Policy>() - 2;
//
// Go ahead and iterate:
//
+ boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
guess = tools::halley_iterate(
detail::gamma_p_inverse_func<T, Policy>(a, p, false),
guess,
lower,
tools::max_value<T>(),
- digits);
+ digits,
+ max_iter);
+ policies::check_root_iterations(function, max_iter, pol);
BOOST_MATH_INSTRUMENT_VARIABLE(guess);
if(guess == lower)
guess = policies::raise_underflow_error<T>(function, "Expected result known to be non-zero, but is smaller than the smallest available number.", pol);
@@ -430,7 +462,10 @@
return tools::max_value<T>();
if(q == 1)
return 0;
- T guess = detail::find_inverse_gamma<T>(a, 1 - q, q, pol);
+ bool has_10_digits;
+ T guess = detail::find_inverse_gamma<T>(a, 1 - q, q, pol, &has_10_digits);
+ if((policies::digits<T, Policy>() <= 36) && has_10_digits)
+ return guess;
T lower = tools::min_value<T>();
if(guess <= lower)
guess = tools::min_value<T>();
@@ -440,18 +475,31 @@
// large convergence is slow, so we'll bump it up to full
// precision to prevent premature termination of the root-finding routine.
//
- unsigned digits = (policies::digits<T, Policy>() * 2) / 3;
+ unsigned digits = policies::digits<T, Policy>();
+ if(digits < 30)
+ {
+ digits *= 2;
+ digits /= 3;
+ }
+ else
+ {
+ digits /= 2;
+ digits -= 1;
+ }
if((a < 0.125) && (fabs(gamma_p_derivative(a, guess, pol)) > 1 / sqrt(tools::epsilon<T>())))
digits = policies::digits<T, Policy>();
//
// Go ahead and iterate:
//
+ boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
guess = tools::halley_iterate(
detail::gamma_p_inverse_func<T, Policy>(a, q, true),
guess,
lower,
tools::max_value<T>(),
- digits);
+ digits,
+ max_iter);
+ policies::check_root_iterations(function, max_iter, pol);
if(guess == lower)
guess = policies::raise_underflow_error<T>(function, "Expected result known to be non-zero, but is smaller than the smallest available number.", pol);
return guess;
Modified: trunk/boost/math/special_functions/erf.hpp
==============================================================================
--- trunk/boost/math/special_functions/erf.hpp (original)
+++ trunk/boost/math/special_functions/erf.hpp 2009-10-05 13:36:56 EDT (Mon, 05 Oct 2009)
@@ -197,13 +197,16 @@
//
// We're going to calculate erf:
//
- if(z == 0)
+ if(z < 1e-10)
{
- result = T(0);
- }
- else if(z < 1e-10)
- {
- result = static_cast<T>(z * 1.125f + z * 0.003379167095512573896158903121545171688L);
+ if(z == 0)
+ {
+ result = T(0);
+ }
+ else
+ {
+ result = static_cast<T>(z * 1.125f + z * 0.003379167095512573896158903121545171688L);
+ }
}
else
{
@@ -227,10 +230,11 @@
0.00858571925074406212772L,
0.000370900071787748000569L,
};
- result = z * (Y + tools::evaluate_polynomial(P, z * z) / tools::evaluate_polynomial(Q, z * z));
+ T zz = z * z;
+ result = z * (Y + tools::evaluate_polynomial(P, zz) / tools::evaluate_polynomial(Q, zz));
}
}
- else if((z < 14) || ((z < 28) && invert))
+ else if(invert ? (z < 28) : (z < 5.8f))
{
//
// We'll be calculating erfc:
@@ -425,7 +429,7 @@
result = z * (Y + tools::evaluate_polynomial(P, z * z) / tools::evaluate_polynomial(Q, z * z));
}
}
- else if((z < 110) || ((z < 110) && invert)) // TODO FIXME!!!
+ else if(invert ? (z < 110) : (z < 6.4f))
{
//
// We'll be calculating erfc:
@@ -634,7 +638,7 @@
result = z * (Y + tools::evaluate_polynomial(P, z * z) / tools::evaluate_polynomial(Q, z * z));
}
}
- else if((z < 9) || ((z < 110) && invert))
+ else if(invert ? (z < 110) : (z < 8.65f))
{
//
// We'll be calculating erfc:
Modified: trunk/boost/math/tools/roots.hpp
==============================================================================
--- trunk/boost/math/tools/roots.hpp (original)
+++ trunk/boost/math/tools/roots.hpp 2009-10-05 13:36:56 EDT (Mon, 05 Oct 2009)
@@ -346,6 +346,8 @@
{
// last two steps haven't converged, try bisection:
delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
+ if(fabs(delta) > result)
+ delta = sign(delta) * result; // protect against huge jumps!
// reset delta2 so that this branch will *not* be taken on the
// next iteration:
delta2 = delta * 3;
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk