|
|
Boost-Commit : |
From: johnmaddock_at_[hidden]
Date: 2007-07-25 14:05:47
Author: johnmaddock
Date: 2007-07-25 14:05:46 EDT (Wed, 25 Jul 2007)
New Revision: 7538
URL: http://svn.boost.org/trac/boost/changeset/7538
Log:
Mostly fixed real_concept failures on HP-UX.
Text files modified:
sandbox/math_toolkit/boost/math/special_functions/detail/lgamma_small.hpp | 29 ++++++++++++++++++++---------
sandbox/math_toolkit/boost/math/special_functions/gamma.hpp | 33 ++++++++++++++++++++++++++-------
2 files changed, 46 insertions(+), 16 deletions(-)
Modified: sandbox/math_toolkit/boost/math/special_functions/detail/lgamma_small.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/lgamma_small.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/lgamma_small.hpp 2007-07-25 14:05:46 EDT (Wed, 25 Jul 2007)
@@ -11,8 +11,8 @@
//
// lgamma for small arguments:
//
-template <class T, class Policy>
-T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<64>&, const Policy& /* l */)
+template <class T, class Policy, class L>
+T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<64>&, const Policy& /* l */, const L&)
{
// This version uses rational approximations for small
// values of z accurate enough for 64-bit mantissas
@@ -200,8 +200,8 @@
}
return result;
}
-template <class T, class Policy>
-T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<113>&, const Policy& /* l */)
+template <class T, class Policy, class L>
+T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<113>&, const Policy& /* l */, const L&)
{
//
// This version uses rational approximations for small
@@ -213,6 +213,7 @@
if(z < tools::epsilon<T>())
{
result = -log(z);
+ BOOST_MATH_INSTRUMENT_CODE(result);
}
else if((zm1 == 0) || (zm2 == 0))
{
@@ -233,6 +234,9 @@
}while(z >= 3);
zm2 = z - 2;
}
+ BOOST_MATH_INSTRUMENT_CODE(zm2);
+ BOOST_MATH_INSTRUMENT_CODE(z);
+ BOOST_MATH_INSTRUMENT_CODE(result);
//
// Use the following form:
@@ -283,6 +287,7 @@
T r = zm2 * (z + 1);
result += r * Y + r * R;
+ BOOST_MATH_INSTRUMENT_CODE(result);
}
else
{
@@ -297,6 +302,9 @@
zm1 = z;
z += 1;
}
+ BOOST_MATH_INSTRUMENT_CODE(result);
+ BOOST_MATH_INSTRUMENT_CODE(z);
+ BOOST_MATH_INSTRUMENT_CODE(zm2);
//
// Three approximations, on for z in [1,1.35], [1.35,1.625] and [1.625,1]
//
@@ -351,6 +359,7 @@
T prefix = zm1 * zm2;
result += prefix * Y + prefix * r;
+ BOOST_MATH_INSTRUMENT_CODE(result);
}
else if(z <= 1.625)
{
@@ -400,6 +409,7 @@
T R = tools::evaluate_polynomial(P, 0.625 - zm1) / tools::evaluate_polynomial(Q, 0.625 - zm1);
result += r * Y + r * R;
+ BOOST_MATH_INSTRUMENT_CODE(result);
}
else
{
@@ -441,14 +451,15 @@
T R = tools::evaluate_polynomial(P, -zm2) / tools::evaluate_polynomial(Q, -zm2);
result += r * Y + r * R;
+ BOOST_MATH_INSTRUMENT_CODE(result);
}
}
+ BOOST_MATH_INSTRUMENT_CODE(result);
return result;
}
-template <class T, class Policy>
-T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<0>&, const Policy& pol)
+template <class T, class Policy, class L>
+T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<0>&, const Policy& pol, const L&)
{
- typedef typename lanczos::lanczos<T, Policy>::type L;
//
// No rational approximations are available because either
// T has no numeric_limits support (so we can't tell how
@@ -465,12 +476,12 @@
else if(z < 0.5)
{
// taking the log of tgamma reduces the error, no danger of overflow here:
- result = log(gamma_imp(z, pol));
+ result = log(gamma_imp(z, pol, L()));
}
else if(z >= 3)
{
// taking the log of tgamma reduces the error, no danger of overflow here:
- result = log(gamma_imp(z, pol));
+ result = log(gamma_imp(z, pol, L()));
}
else if(z >= 1.5)
{
Modified: sandbox/math_toolkit/boost/math/special_functions/gamma.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/gamma.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/gamma.hpp 2007-07-25 14:05:46 EDT (Wed, 25 Jul 2007)
@@ -228,7 +228,7 @@
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);
+ result = lgamma_small_imp(z, z - 1, z - 2, tag_type(), pol, l);
}
else if((z >= 3) && (z < 100))
{
@@ -356,6 +356,7 @@
prefix /= z;
z += 1;
}
+ BOOST_MATH_INSTRUMENT_CODE(prefix);
if((floor(z) == z) && (z < max_factorial<T>::value))
{
prefix *= unchecked_factorial<T>(tools::real_cast<unsigned>(z) - 1);
@@ -363,10 +364,14 @@
else
{
prefix = prefix * pow(z / boost::math::constants::e<T>(), z);
+ BOOST_MATH_INSTRUMENT_CODE(prefix);
T sum = detail::lower_gamma_series(z, z, pol) / z;
+ BOOST_MATH_INSTRUMENT_CODE(sum);
sum += detail::upper_gamma_fraction(z, z, ::boost::math::policy::digits<T, Policy>());
+ BOOST_MATH_INSTRUMENT_CODE(sum);
if(fabs(tools::max_value<T>() / prefix) < fabs(sum))
return policy::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
+ BOOST_MATH_INSTRUMENT_CODE((sum * prefix));
return sum * prefix;
}
return prefix;
@@ -417,12 +422,17 @@
{
using namespace std;
+ typedef typename policy::precision<T,Policy>::type precision_type;
+
typedef typename mpl::if_<
- mpl::less_equal<typename policy::precision<T,Policy>::type, mpl::int_<64> >,
- mpl::int_<64>,
+ mpl::or_<
+ mpl::less_equal<precision_type, mpl::int_<0> >,
+ mpl::greater<precision_type, mpl::int_<113> >
+ >,
+ mpl::int_<0>,
typename mpl::if_<
- mpl::less_equal<typename policy::precision<T,Policy>::type, mpl::int_<113> >,
- mpl::int_<113>, mpl::int_<0> >::type
+ mpl::less_equal<precision_type, mpl::int_<64> >,
+ mpl::int_<64>, mpl::int_<113> >::type
>::type tag_type;
T result;
@@ -432,12 +442,14 @@
{
// Best method is simply to subtract 1 from tgamma:
result = boost::math::tgamma(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));
+ + lgamma_small_imp(dz+2, dz + 1, dz, tag_type(), pol, l));
+ BOOST_MATH_INSTRUMENT_CODE(result);
}
}
else
@@ -445,12 +457,14 @@
if(dz < 2)
{
// Use expm1 on lgamma:
- result = boost::math::expm1(lgamma_small_imp(dz+1, dz, dz-1, tag_type(), pol), pol);
+ result = boost::math::expm1(lgamma_small_imp(dz+1, dz, 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;
+ BOOST_MATH_INSTRUMENT_CODE(result);
}
}
@@ -468,15 +482,20 @@
// Start by subracting 1 from tgamma:
//
T result = gamma_imp(1 + dz, pol, l) - 1;
+ BOOST_MATH_INSTRUMENT_CODE(result);
//
// Test the level of cancellation error observed: we loose one bit
// for each power of 2 the result is less than 1. If we would get
// more bits from our most precise lgamma rational approximation,
// then use that instead:
//
+ BOOST_MATH_INSTRUMENT_CODE((dz > -0.5));
+ BOOST_MATH_INSTRUMENT_CODE((dz < 2));
+ BOOST_MATH_INSTRUMENT_CODE((ldexp(1.0, boost::math::policy::digits<T, Policy>()) * fabs(result) < 1e34));
if((dz > -0.5) && (dz < 2) && (ldexp(1.0, boost::math::policy::digits<T, Policy>()) * fabs(result) < 1e34))
{
result = tgammap1m1_imp(dz, pol, boost::math::lanczos::lanczos24m113());
+ BOOST_MATH_INSTRUMENT_CODE(result);
}
return result;
}
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