|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r49828 - in trunk/boost/math/special_functions: . detail
From: john_at_[hidden]
Date: 2008-11-18 14:47:51
Author: johnmaddock
Date: 2008-11-18 14:47:50 EST (Tue, 18 Nov 2008)
New Revision: 49828
URL: http://svn.boost.org/trac/boost/changeset/49828
Log:
Fixed bug in elliptic integral periodic-reduction code.
Added some more debugging/tracing aids.
Text files modified:
trunk/boost/math/special_functions/bessel.hpp | 23 +++++++++++-
trunk/boost/math/special_functions/beta.hpp | 74 ++++++++++++++++++++++++++++++++++++++++
trunk/boost/math/special_functions/detail/igamma_inverse.hpp | 31 ++++++++++++++++
trunk/boost/math/special_functions/ellint_1.hpp | 2
trunk/boost/math/special_functions/ellint_2.hpp | 2
trunk/boost/math/special_functions/ellint_3.hpp | 2
6 files changed, 129 insertions(+), 5 deletions(-)
Modified: trunk/boost/math/special_functions/bessel.hpp
==============================================================================
--- trunk/boost/math/special_functions/bessel.hpp (original)
+++ trunk/boost/math/special_functions/bessel.hpp 2008-11-18 14:47:50 EST (Tue, 18 Nov 2008)
@@ -307,6 +307,10 @@
inline T cyl_neumann_imp(T v, T x, const bessel_no_int_tag&, const Policy& pol)
{
static const char* function = "boost::math::cyl_neumann<%1%>(%1%,%1%)";
+
+ BOOST_MATH_INSTRUMENT_VARIABLE(v);
+ BOOST_MATH_INSTRUMENT_VARIABLE(x);
+
if(x <= 0)
{
return (v == 0) && (x == 0) ?
@@ -332,6 +336,10 @@
{
BOOST_MATH_STD_USING
typedef typename bessel_asymptotic_tag<T, Policy>::type tag_type;
+
+ BOOST_MATH_INSTRUMENT_VARIABLE(v);
+ BOOST_MATH_INSTRUMENT_VARIABLE(x);
+
if(floor(v) == v)
{
if((fabs(x) > asymptotic_bessel_y_limit<T>(tag_type())) && (fabs(x) > 5 * abs(v)))
@@ -339,12 +347,19 @@
T r = asymptotic_bessel_y_large_x_2(static_cast<T>(abs(v)), x);
if((v < 0) && (itrunc(v, pol) & 1))
r = -r;
+ BOOST_MATH_INSTRUMENT_VARIABLE(r);
return r;
}
else
- return bessel_yn(itrunc(v, pol), x, pol);
+ {
+ T r = bessel_yn(itrunc(v, pol), x, pol);
+ BOOST_MATH_INSTRUMENT_VARIABLE(r);
+ return r;
+ }
}
- return cyl_neumann_imp<T>(v, x, bessel_no_int_tag(), pol);
+ T r = cyl_neumann_imp<T>(v, x, bessel_no_int_tag(), pol);
+ BOOST_MATH_INSTRUMENT_VARIABLE(r);
+ return r;
}
template <class T, class Policy>
@@ -352,6 +367,10 @@
{
BOOST_MATH_STD_USING
typedef typename bessel_asymptotic_tag<T, Policy>::type tag_type;
+
+ BOOST_MATH_INSTRUMENT_VARIABLE(v);
+ BOOST_MATH_INSTRUMENT_VARIABLE(x);
+
if((fabs(x) > asymptotic_bessel_y_limit<T>(tag_type())) && (fabs(x) > 5 * abs(v)))
{
T r = asymptotic_bessel_y_large_x_2(static_cast<T>(abs(v)), x);
Modified: trunk/boost/math/special_functions/beta.hpp
==============================================================================
--- trunk/boost/math/special_functions/beta.hpp (original)
+++ trunk/boost/math/special_functions/beta.hpp 2008-11-18 14:47:50 EST (Tue, 18 Nov 2008)
@@ -242,13 +242,25 @@
// since one of the power terms will evaluate to a number close to 1.
//
if(fabs(l1) < 0.1)
+ {
result *= exp(a * boost::math::log1p(l1, pol));
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
+ }
else
+ {
result *= pow((x * cgh) / agh, a);
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
+ }
if(fabs(l2) < 0.1)
+ {
result *= exp(b * boost::math::log1p(l2, pol));
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
+ }
else
+ {
result *= pow((y * cgh) / bgh, b);
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
+ }
}
else if((std::max)(fabs(l1), fabs(l2)) < 0.5)
{
@@ -279,6 +291,7 @@
l3 = l1 + l3 + l3 * l1;
l3 = a * boost::math::log1p(l3, pol);
result *= exp(l3);
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
else
{
@@ -286,6 +299,7 @@
l3 = l2 + l3 + l3 * l2;
l3 = b * boost::math::log1p(l3, pol);
result *= exp(l3);
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
}
else if(fabs(l1) < fabs(l2))
@@ -294,6 +308,7 @@
T l = a * boost::math::log1p(l1, pol)
+ b * log((y * cgh) / bgh);
result *= exp(l);
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
else
{
@@ -301,6 +316,7 @@
T l = b * boost::math::log1p(l2, pol)
+ a * log((x * cgh) / agh);
result *= exp(l);
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
}
else
@@ -321,11 +337,13 @@
result *= pow(pow(b2, b/a) * b1, a);
else
result *= pow(pow(b1, a/b) * b2, b);
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
else
{
// finally the normal case:
result *= pow(b1, a) * pow(b2, b);
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
}
// combine with the leftover terms from the Lanczos approximation:
@@ -333,6 +351,8 @@
result *= sqrt(agh / cgh);
result *= prefix;
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
+
return result;
}
//
@@ -624,6 +644,9 @@
T ibeta_a_step(T a, T b, T x, T y, int k, const Policy& pol, bool normalised, T* p_derivative)
{
typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
+
+ BOOST_MATH_INSTRUMENT_VARIABLE(k);
+
T prefix = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol);
if(p_derivative)
{
@@ -662,6 +685,7 @@
// This is only called with small k, for large k
// it is grossly inefficient, do not use outside it's
// intended purpose!!!
+ BOOST_MATH_INSTRUMENT_VARIABLE(k);
if(k == 0)
return 1;
T result = 1;
@@ -845,6 +869,12 @@
typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
BOOST_MATH_STD_USING // for ADL of std math functions.
+ BOOST_MATH_INSTRUMENT_VARIABLE(a);
+ BOOST_MATH_INSTRUMENT_VARIABLE(b);
+ BOOST_MATH_INSTRUMENT_VARIABLE(x);
+ BOOST_MATH_INSTRUMENT_VARIABLE(inv);
+ BOOST_MATH_INSTRUMENT_VARIABLE(normalised);
+
bool invert = inv;
T fract;
T y = 1 - x;
@@ -894,6 +924,7 @@
std::swap(a, b);
std::swap(x, y);
invert = !invert;
+ BOOST_MATH_INSTRUMENT_VARIABLE(invert);
}
if((std::max)(a, b) <= 1)
{
@@ -901,12 +932,16 @@
if((a >= (std::min)(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);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
+ }
else
{
fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
invert = false;
fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
}
}
else
@@ -917,12 +952,16 @@
if(y >= 0.3)
{
if(!invert)
+ {
fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
+ }
else
{
fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
invert = false;
fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
}
}
else
@@ -939,12 +978,16 @@
}
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);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
+ }
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);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
}
}
}
@@ -955,12 +998,16 @@
if((b <= 1) || ((x < 0.1) && (pow(b * x, a) <= 0.7)))
{
if(!invert)
+ {
fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
+ }
else
{
fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
invert = false;
fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
}
}
else
@@ -972,23 +1019,31 @@
if(y >= 0.3)
{
if(!invert)
+ {
fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
+ }
else
{
fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
invert = false;
fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
}
}
else if(a >= 15)
{
if(!invert)
+ {
fract = beta_small_b_large_a_series(a, b, x, y, T(0), T(1), pol, normalised);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
+ }
else
{
fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
invert = false;
fract = -beta_small_b_large_a_series(a, b, x, y, fract, T(1), pol, normalised);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
}
}
else
@@ -1004,13 +1059,18 @@
prefix = 1;
}
fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
if(!invert)
+ {
fract = beta_small_b_large_a_series(a + 20, b, x, y, fract, prefix, pol, normalised);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
+ }
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);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
}
}
}
@@ -1033,6 +1093,7 @@
std::swap(a, b);
std::swap(x, y);
invert = !invert;
+ BOOST_MATH_INSTRUMENT_VARIABLE(invert);
}
if(b < 40)
@@ -1045,16 +1106,21 @@
fract = binomial_ccdf(n, k, x, y);
if(!normalised)
fract *= boost::math::beta(a, b, pol);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
}
else if(b * x <= 0.7)
{
if(!invert)
+ {
fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
+ }
else
{
fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
invert = false;
fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
}
}
else if(a > 15)
@@ -1076,6 +1142,7 @@
fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(0));
fract = beta_small_b_large_a_series(a, bbar, x, y, fract, T(1), pol, normalised);
fract /= prefix;
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
}
else if(normalised)
{
@@ -1100,12 +1167,19 @@
fract = -fract;
invert = false;
}
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
}
else
+ {
fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
+ }
}
else
+ {
fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
+ BOOST_MATH_INSTRUMENT_VARIABLE(fract);
+ }
}
if(p_derivative)
{
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 2008-11-18 14:47:50 EST (Tue, 18 Nov 2008)
@@ -109,11 +109,16 @@
T result;
if(a == 1)
+ {
result = -log(q);
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
+ }
else if(a < 1)
{
T g = boost::math::tgamma(a, pol);
T b = q * g;
+ BOOST_MATH_INSTRUMENT_VARIABLE(g);
+ BOOST_MATH_INSTRUMENT_VARIABLE(b);
if((b > 0.6) || ((b >= 0.45) && (a >= 0.3)))
{
// DiDonato & Morris Eq 21:
@@ -127,12 +132,15 @@
if((b * q > 1e-8) && (q > 1e-5))
{
u = pow(p * g * a, 1 / a);
+ BOOST_MATH_INSTRUMENT_VARIABLE(u);
}
else
{
u = exp((-q / a) - constants::euler<T>());
+ BOOST_MATH_INSTRUMENT_VARIABLE(u);
}
result = u / (1 - (u / (a + 1)));
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
else if((a < 0.3) && (b >= 0.35))
{
@@ -140,6 +148,7 @@
T t = exp(-constants::euler<T>() - b);
T u = t * exp(t);
result = t * exp(u);
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
else if((b > 0.15) || (a >= 0.3))
{
@@ -147,6 +156,7 @@
T y = -log(b);
T u = y - (1 - a) * log(y);
result = y - (1 - a) * log(u) - log(1 + (1 - a) / (1 + u));
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
else if (b > 0.1)
{
@@ -154,6 +164,7 @@
T y = -log(b);
T u = y - (1 - a) * log(y);
result = y - (1 - a) * log(u) - log((u * u + 2 * (3 - a) * u + (2 - a) * (3 - a)) / (u * u + (5 - a) * u + 2));
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
else
{
@@ -179,6 +190,7 @@
T y_3 = y_2 * y;
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);
}
}
else
@@ -186,26 +198,34 @@
// DiDonato and Morris Eq 31:
T s = find_inverse_s(p, q);
+ BOOST_MATH_INSTRUMENT_VARIABLE(s);
+
T s_2 = s * s;
T s_3 = s_2 * s;
T s_4 = s_2 * s_2;
T s_5 = s_4 * s;
T ra = sqrt(a);
+ BOOST_MATH_INSTRUMENT_VARIABLE(ra);
+
T w = a + s * ra + (s * s -1) / 3;
w += (s_3 - 7 * s) / (36 * ra);
w -= (3 * s_4 + 7 * s_2 - 16) / (810 * a);
w += (9 * s_5 + 256 * s_3 - 433 * s) / (38880 * a * ra);
+ BOOST_MATH_INSTRUMENT_VARIABLE(w);
+
if((a >= 500) && (fabs(1 - w / a) < 1e-6))
{
result = w;
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
else if (p > 0.5)
{
if(w < 3 * a)
{
result = w;
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
else
{
@@ -236,12 +256,14 @@
T y_3 = y_2 * y;
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);
}
else
{
// DiDonato and Morris Eq 33:
T u = -lb + (a - 1) * log(w) - log(1 + (1 - a) / (1 + w));
result = -lb + (a - 1) * log(u) - log(1 + (1 - a) / (1 + u));
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
}
}
@@ -253,9 +275,12 @@
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);
+
if((z <= 0.01 * (a + 1)) || (z > 0.7 * (a + 1)))
{
result = z;
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
else
{
@@ -263,6 +288,7 @@
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));
+ BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
}
}
@@ -349,6 +375,9 @@
static const char* function = "boost::math::gamma_p_inv<%1%>(%1%, %1%)";
+ BOOST_MATH_INSTRUMENT_VARIABLE(a);
+ BOOST_MATH_INSTRUMENT_VARIABLE(p);
+
if(a <= 0)
policies::raise_domain_error<T>(function, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a, pol);
if((p < 0) || (p > 1))
@@ -361,6 +390,7 @@
T lower = tools::min_value<T>();
if(guess <= lower)
guess = tools::min_value<T>();
+ BOOST_MATH_INSTRUMENT_VARIABLE(guess);
//
// Work out how many digits to converge to, normally this is
// 2/3 of the digits in T, but if the first derivative is very
@@ -379,6 +409,7 @@
lower,
tools::max_value<T>(),
digits);
+ 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);
return guess;
Modified: trunk/boost/math/special_functions/ellint_1.hpp
==============================================================================
--- trunk/boost/math/special_functions/ellint_1.hpp (original)
+++ trunk/boost/math/special_functions/ellint_1.hpp 2008-11-18 14:47:50 EST (Tue, 18 Nov 2008)
@@ -90,7 +90,7 @@
BOOST_MATH_INSTRUMENT_CODE("pi/2 = " << constants::pi<T>() / 2);
T rphi = boost::math::tools::fmod_workaround(phi, constants::pi<T>() / 2);
BOOST_MATH_INSTRUMENT_VARIABLE(rphi);
- T m = 2 * (phi - rphi) / constants::pi<T>();
+ T m = floor((2 * phi) / constants::pi<T>());
BOOST_MATH_INSTRUMENT_VARIABLE(m);
int s = 1;
if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5)
Modified: trunk/boost/math/special_functions/ellint_2.hpp
==============================================================================
--- trunk/boost/math/special_functions/ellint_2.hpp (original)
+++ trunk/boost/math/special_functions/ellint_2.hpp 2008-11-18 14:47:50 EST (Tue, 18 Nov 2008)
@@ -75,7 +75,7 @@
// 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 m = floor((2 * phi) / constants::pi<T>());
int s = 1;
if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5)
{
Modified: trunk/boost/math/special_functions/ellint_3.hpp
==============================================================================
--- trunk/boost/math/special_functions/ellint_3.hpp (original)
+++ trunk/boost/math/special_functions/ellint_3.hpp 2008-11-18 14:47:50 EST (Tue, 18 Nov 2008)
@@ -183,7 +183,7 @@
else
{
T rphi = boost::math::tools::fmod_workaround(fabs(phi), constants::pi<T>() / 2);
- T m = 2 * (fabs(phi) - rphi) / constants::pi<T>();
+ T m = floor((2 * fabs(phi)) / constants::pi<T>());
int sign = 1;
if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5)
{
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