Boost logo

Boost-Commit :

From: john_at_[hidden]
Date: 2008-05-13 12:43:47


Author: johnmaddock
Date: 2008-05-13 12:43:46 EDT (Tue, 13 May 2008)
New Revision: 45332
URL: http://svn.boost.org/trac/boost/changeset/45332

Log:
Fixed float_next and added extra tests to detect the bug.
Fixed expected assoc_legendre results.
Text files modified:
   sandbox/math_toolkit/boost/math/special_functions/next.hpp | 41 ++++++++++++++++++++++-----------------
   sandbox/math_toolkit/libs/math/test/test_next.cpp | 2
   sandbox/math_toolkit/libs/math/test/test_tr1.cpp | 20 +++++++++---------
   3 files changed, 34 insertions(+), 29 deletions(-)

Modified: sandbox/math_toolkit/boost/math/special_functions/next.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/next.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/next.hpp 2008-05-13 12:43:46 EDT (Tue, 13 May 2008)
@@ -185,34 +185,40 @@
    if(boost::math::sign(a) != boost::math::sign(b))
       return 2 + fabs(float_distance(boost::math::sign(b) * detail::get_smallest_value<T>(), b, pol))
          + fabs(float_distance(boost::math::sign(a) * detail::get_smallest_value<T>(), a, pol));
+ //
+ // By the time we get here, both a and b must have the same sign, we want
+ // b > a and both postive for the following logic:
+ //
+ if(a < 0)
+ return float_distance(-b, -a);
 
- if((std::min)(fabs(a), fabs(b)) / (std::max)(fabs(a), fabs(b)) < 2 * tools::epsilon<T>())
- {
- bool biga = fabs(a) > fabs(b);
- T split = ldexp(biga ? b : a, tools::digits<T>() - 2);
- return fabs(float_distance(a, split, pol) + float_distance(split, b, pol));
- }
+ BOOST_ASSERT(a >= 0);
+ BOOST_ASSERT(b >= a);
 
    BOOST_MATH_STD_USING
    int expon;
    //
- // We're going to left shift the result by the exponent of the
- // smaller of the two values (irrespective of sign):
- //
- T mv = (std::min)(fabs(a), fabs(b));
- //
- // Note that if mv is a denorm then the usual formula fails
+ // Note that if a is a denorm then the usual formula fails
    // because we actually have fewer than tools::digits<T>()
    // significant bits in the representation:
    //
- frexp((boost::math::fpclassify(mv) == FP_SUBNORMAL) ? tools::min_value<T>() : mv, &expon);
+ frexp((boost::math::fpclassify(a) == FP_SUBNORMAL) ? tools::min_value<T>() : a, &expon);
+ T upper = ldexp(T(1), expon);
+ T result = 0;
    expon = tools::digits<T>() - expon;
    //
+ // If b is greater than upper, then we *must* split the calculation
+ // as the size of the ULP changes with each order of magnitude change:
+ //
+ if(b > upper)
+ {
+ result = float_distance(upper, b);
+ }
+ //
    // Use compensated double-double addition to avoid rounding
- // errors in the subtraction, note this will still fail if
- // the two values differ by many orders of magnitute:
+ // errors in the subtraction:
    //
- T mb = -b;
+ T mb = -(std::min)(upper, b);
    T x = a + mb;
    T z = x - a;
    T y = (a - (x - z)) + (mb - z);
@@ -221,8 +227,7 @@
       x = -x;
       y = -y;
    }
-
- T result = ldexp(x, expon) + ldexp(y, expon);
+ result += ldexp(x, expon) + ldexp(y, expon);
    //
    // Result must be an integer:
    //

Modified: sandbox/math_toolkit/libs/math/test/test_next.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/test/test_next.cpp (original)
+++ sandbox/math_toolkit/libs/math/test/test_next.cpp 2008-05-13 12:43:46 EDT (Tue, 13 May 2008)
@@ -31,6 +31,7 @@
    BOOST_CHECK(nextafter(val, lower) < val);
    BOOST_CHECK_EQUAL(float_distance(float_next(float_next(val)), val), -2);
    BOOST_CHECK_EQUAL(float_distance(float_prior(float_prior(val)), val), 2);
+ BOOST_CHECK_EQUAL(float_distance(float_prior(float_prior(val)), float_next(float_next(val))), 4);
    BOOST_CHECK_EQUAL(float_distance(float_prior(float_next(val)), val), 0);
    BOOST_CHECK_EQUAL(float_distance(float_next(float_prior(val)), val), 0);
    BOOST_CHECK_EQUAL(float_prior(float_next(val)), val);
@@ -96,7 +97,6 @@
 
 int test_main(int, char* [])
 {
- std::cout << boost::math::float_distance(1.0, 0.0) << std::endl;
    test_values(1.0f, "float");
    test_values(1.0, "double");
    test_values(1.0L, "long double");

Modified: sandbox/math_toolkit/libs/math/test/test_tr1.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/test/test_tr1.cpp (original)
+++ sandbox/math_toolkit/libs/math/test/test_tr1.cpp 2008-05-13 12:43:46 EDT (Tue, 13 May 2008)
@@ -192,9 +192,9 @@
    BOOST_CHECK_CLOSE_FRACTION(tr1::laguerref(50, static_cast<float>(4.5L)), static_cast<float>(-0.7795068145562651416494321484050019245248L), eps * 100);
 
    BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendref(4, 2, static_cast<float>(0.5L)), static_cast<float>(4.218750000000000000000000000000000000000L), eps * 100);
- BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendref(7, 5, static_cast<float>(0.5L)), static_cast<float>(-5696.789530152175143607977274672800795328L), eps * 100);
+ BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendref(7, 5, static_cast<float>(0.5L)), static_cast<float>(5696.789530152175143607977274672800795328L), eps * 100);
    BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendref(4, 2, static_cast<float>(-0.5L)), static_cast<float>(4.218750000000000000000000000000000000000L), eps * 100);
- BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendref(7, 5, static_cast<float>(-0.5L)), static_cast<float>(-5696.789530152175143607977274672800795328L), eps * 100);
+ BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendref(7, 5, static_cast<float>(-0.5L)), static_cast<float>(5696.789530152175143607977274672800795328L), eps * 100);
 
    BOOST_CHECK_CLOSE_FRACTION(tr1::legendref(1, static_cast<float>(0.5L)), static_cast<float>(0.5L), eps * 100);
    BOOST_CHECK_CLOSE_FRACTION(tr1::legendref(4, static_cast<float>(0.5L)), static_cast<float>(-0.2890625000000000000000000000000000000000L), eps * 100);
@@ -434,9 +434,9 @@
    BOOST_CHECK_CLOSE_FRACTION(tr1::laguerre(50, static_cast<float>(4.5L)), static_cast<float>(-0.7795068145562651416494321484050019245248L), eps * 100);
 
    BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendre(4, 2, static_cast<float>(0.5L)), static_cast<float>(4.218750000000000000000000000000000000000L), eps * 100);
- BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendre(7, 5, static_cast<float>(0.5L)), static_cast<float>(-5696.789530152175143607977274672800795328L), eps * 100);
+ BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendre(7, 5, static_cast<float>(0.5L)), static_cast<float>(5696.789530152175143607977274672800795328L), eps * 100);
    BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendre(4, 2, static_cast<float>(-0.5L)), static_cast<float>(4.218750000000000000000000000000000000000L), eps * 100);
- BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendre(7, 5, static_cast<float>(-0.5L)), static_cast<float>(-5696.789530152175143607977274672800795328L), eps * 100);
+ BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendre(7, 5, static_cast<float>(-0.5L)), static_cast<float>(5696.789530152175143607977274672800795328L), eps * 100);
 
    BOOST_CHECK_CLOSE_FRACTION(tr1::legendre(1, static_cast<float>(0.5L)), static_cast<float>(0.5L), eps * 100);
    BOOST_CHECK_CLOSE_FRACTION(tr1::legendre(4, static_cast<float>(0.5L)), static_cast<float>(-0.2890625000000000000000000000000000000000L), eps * 100);
@@ -750,9 +750,9 @@
    BOOST_CHECK_CLOSE_FRACTION(tr1::laguerre(50, static_cast<double>(4.5L)), static_cast<double>(-0.7795068145562651416494321484050019245248L), eps * 100);
 
    BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendre(4, 2, static_cast<double>(0.5L)), static_cast<double>(4.218750000000000000000000000000000000000L), eps * 100);
- BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendre(7, 5, static_cast<double>(0.5L)), static_cast<double>(-5696.789530152175143607977274672800795328L), eps * 100);
+ BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendre(7, 5, static_cast<double>(0.5L)), static_cast<double>(5696.789530152175143607977274672800795328L), eps * 100);
    BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendre(4, 2, static_cast<double>(-0.5L)), static_cast<double>(4.218750000000000000000000000000000000000L), eps * 100);
- BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendre(7, 5, static_cast<double>(-0.5L)), static_cast<double>(-5696.789530152175143607977274672800795328L), eps * 100);
+ BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendre(7, 5, static_cast<double>(-0.5L)), static_cast<double>(5696.789530152175143607977274672800795328L), eps * 100);
 
    BOOST_CHECK_CLOSE_FRACTION(tr1::legendre(1, static_cast<double>(0.5L)), static_cast<double>(0.5L), eps * 100);
    BOOST_CHECK_CLOSE_FRACTION(tr1::legendre(4, static_cast<double>(0.5L)), static_cast<double>(-0.2890625000000000000000000000000000000000L), eps * 100);
@@ -1140,9 +1140,9 @@
    BOOST_CHECK_CLOSE_FRACTION(tr1::laguerrel(50L, static_cast<long double>(4.5L)), static_cast<long double>(-0.7795068145562651416494321484050019245248L), eps * 100L);
 
    BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendrel(4L, 2L, static_cast<long double>(0.5L)), static_cast<long double>(4.218750000000000000000000000000000000000L), eps * 100L);
- BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendrel(7L, 5L, static_cast<long double>(0.5L)), static_cast<long double>(-5696.789530152175143607977274672800795328L), eps * 100L);
+ BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendrel(7L, 5L, static_cast<long double>(0.5L)), static_cast<long double>(5696.789530152175143607977274672800795328L), eps * 100L);
    BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendrel(4L, 2L, static_cast<long double>(-0.5L)), static_cast<long double>(4.218750000000000000000000000000000000000L), eps * 100L);
- BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendrel(7L, 5L, static_cast<long double>(-0.5L)), static_cast<long double>(-5696.789530152175143607977274672800795328L), eps * 100L);
+ BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendrel(7L, 5L, static_cast<long double>(-0.5L)), static_cast<long double>(5696.789530152175143607977274672800795328L), eps * 100L);
 
    BOOST_CHECK_CLOSE_FRACTION(tr1::legendrel(1L, static_cast<long double>(0.5L)), static_cast<long double>(0.5L), eps * 100L);
    BOOST_CHECK_CLOSE_FRACTION(tr1::legendrel(4L, static_cast<long double>(0.5L)), static_cast<long double>(-0.2890625000000000000000000000000000000000L), eps * 100L);
@@ -1382,9 +1382,9 @@
    BOOST_CHECK_CLOSE_FRACTION(tr1::laguerre(50L, static_cast<long double>(4.5L)), static_cast<long double>(-0.7795068145562651416494321484050019245248L), eps * 100L);
 
    BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendre(4L, 2L, static_cast<long double>(0.5L)), static_cast<long double>(4.218750000000000000000000000000000000000L), eps * 100L);
- BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendre(7L, 5L, static_cast<long double>(0.5L)), static_cast<long double>(-5696.789530152175143607977274672800795328L), eps * 100L);
+ BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendre(7L, 5L, static_cast<long double>(0.5L)), static_cast<long double>(5696.789530152175143607977274672800795328L), eps * 100L);
    BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendre(4L, 2L, static_cast<long double>(-0.5L)), static_cast<long double>(4.218750000000000000000000000000000000000L), eps * 100L);
- BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendre(7L, 5L, static_cast<long double>(-0.5L)), static_cast<long double>(-5696.789530152175143607977274672800795328L), eps * 100L);
+ BOOST_CHECK_CLOSE_FRACTION(tr1::assoc_legendre(7L, 5L, static_cast<long double>(-0.5L)), static_cast<long double>(5696.789530152175143607977274672800795328L), eps * 100L);
 
    BOOST_CHECK_CLOSE_FRACTION(tr1::legendre(1L, static_cast<long double>(0.5L)), static_cast<long double>(0.5L), eps * 100L);
    BOOST_CHECK_CLOSE_FRACTION(tr1::legendre(4L, static_cast<long double>(0.5L)), static_cast<long double>(-0.2890625000000000000000000000000000000000L), eps * 100L);


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