From: john_at_[hidden]
Date: 2008-04-29 08:01:22
Author: johnmaddock
Date: 2008-04-29 08:01:22 EDT (Tue, 29 Apr 2008)
New Revision: 44878
URL: http://svn.boost.org/trac/boost/changeset/44878
Log:
Added first cut of nextafter family of functions.
Added:
sandbox/math_toolkit/boost/math/special_functions/next.hpp (contents, props changed)
sandbox/math_toolkit/libs/math/test/compile_test/sf_next_incl_test.cpp (contents, props changed)
sandbox/math_toolkit/libs/math/test/test_next.cpp (contents, props changed)
Text files modified:
sandbox/math_toolkit/boost/math/special_functions.hpp | 1 +
sandbox/math_toolkit/boost/math/special_functions/math_fwd.hpp | 23 +++++++++++++++++++++++
sandbox/math_toolkit/libs/math/test/Jamfile.v2 | 3 +++
sandbox/math_toolkit/libs/math/test/compile_test/instantiate.hpp | 12 ++++++++++++
4 files changed, 39 insertions(+), 0 deletions(-)
Modified: sandbox/math_toolkit/boost/math/special_functions.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions.hpp 2008-04-29 08:01:22 EDT (Tue, 29 Apr 2008)
@@ -41,6 +41,7 @@
#include <boost/math/special_functions/legendre.hpp>
#include <boost/math/special_functions/log1p.hpp>
#include <boost/math/special_functions/math_fwd.hpp>
+#include <boost/math/special_functions/next.hpp>
#include <boost/math/special_functions/powm1.hpp>
#include <boost/math/special_functions/sign.hpp>
#include <boost/math/special_functions/sin_pi.hpp>
Modified: sandbox/math_toolkit/boost/math/special_functions/math_fwd.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/math_fwd.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/math_fwd.hpp 2008-04-29 08:01:22 EDT (Tue, 29 Apr 2008)
@@ -657,6 +657,24 @@
template <class T>
typename tools::promote_args<T>::type zeta(T s);
+ // next:
+ template <class T, class Policy>
+ T nextafter(const T&, const T&, const Policy&);
+ template <class T>
+ T nextafter(const T&, const T&);
+ template <class T, class Policy>
+ T next_greater(const T&, const Policy&);
+ template <class T>
+ T next_greater(const T&);
+ template <class T, class Policy>
+ T next_less(const T&, const Policy&);
+ template <class T>
+ T next_less(const T&);
+ template <class T, class Policy>
+ T edit_distance(const T&, const T&, const Policy&);
+ template <class T>
+ T edit_distance(const T&, const T&);
+
} // namespace math
} // namespace boost
@@ -1003,6 +1021,11 @@
\
template <int N, class T>\
inline typename boost::math::tools::promote_args<T>::type pow(T v){ return boost::math::pow<N>(v, Policy()); }\
+ template <class T> T nextafter(const T& a, const T& b){ return boost::math::nextafter(a, b, Policy()); }\
+ template <class T> T next_greater(const T& a){ return boost::math::next_greater(a, Policy()); }\
+ template <class T> T next_less(const T& a){ return boost::math::next_less(a, Policy()); }\
+ template <class T> T edit_distance(const T& a, const T& b){ return boost::math::edit_distance(a, b, Policy()); }\
+
#endif // BOOST_MATH_SPECIAL_MATH_FWD_HPP
Added: sandbox/math_toolkit/boost/math/special_functions/next.hpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/boost/math/special_functions/next.hpp 2008-04-29 08:01:22 EDT (Tue, 29 Apr 2008)
@@ -0,0 +1,240 @@
+// (C) Copyright John Maddock 2008.
+// 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_SPECIAL_NEXT_HPP
+#define BOOST_MATH_SPECIAL_NEXT_HPP
+
+#ifdef _MSC_VER
+#pragma once
+#endif
+
+#include <boost/math/policies/error_handling.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/math/special_functions/sign.hpp>
+
+#ifdef BOOST_MSVC
+#include <float.h>
+#endif
+
+namespace boost{ namespace math{
+
+namespace detail{
+
+template <class T>
+inline T get_smallest_value(mpl::true_ const&)
+{
+ return std::numeric_limits<T>::denorm_min();
+}
+
+template <class T>
+inline T get_smallest_value(mpl::false_ const&)
+{
+ return tools::min_value<T>();
+}
+
+template <class T>
+inline T get_smallest_value()
+{
+ return get_smallest_value<T>(mpl::bool_<std::numeric_limits<T>::is_specialized && std::numeric_limits<T>::has_denorm>());
+}
+
+}
+
+template <class T, class Policy>
+T next_greater(const T& val, const Policy& pol)
+{
+ BOOST_MATH_STD_USING
+ int expon;
+ static const char* function = "next_greater<%1%>(%1%)";
+
+ if(!(boost::math::isfinite)(val))
+ return policies::raise_domain_error<T>(
+ function,
+ "Argument must be finite, but got %1%", val, pol);
+
+ if(val >= tools::max_value<T>())
+ return policies::raise_overflow_error<T>(function, 0, pol);
+
+ if(val == 0)
+ return detail::get_smallest_value<T>();
+
+ if(-0.5f == frexp(val, &expon))
+ --expon; // reduce exponent when val is a power of two, and negative.
+ T diff = ldexp(T(1), expon - tools::digits<T>());
+ if(diff == 0)
+ diff = detail::get_smallest_value<T>();
+ return val + diff;
+}
+
+#ifdef BOOST_MSVC
+template <class Policy>
+inline double next_greater(const double& val, const Policy& pol)
+{
+ static const char* function = "next_greater<%1%>(%1%)";
+
+ if(!(boost::math::isfinite)(val))
+ return policies::raise_domain_error<double>(
+ function,
+ "Argument must be finite, but got %1%", val, pol);
+
+ if(val >= tools::max_value<double>())
+ return policies::raise_overflow_error<double>(function, 0, pol);
+
+ return ::_nextafter(val, tools::max_value<double>());
+}
+#endif
+
+template <class T>
+inline T next_greater(const T& val)
+{
+ return next_greater(val, policies::policy<>());
+}
+
+template <class T, class Policy>
+T next_less(const T& val, const Policy& pol)
+{
+ BOOST_MATH_STD_USING
+ int expon;
+ static const char* function = "next_less<%1%>(%1%)";
+
+ if(!(boost::math::isfinite)(val))
+ return policies::raise_domain_error<T>(
+ function,
+ "Argument must be finite, but got %1%", val, pol);
+
+ if(val <= -tools::max_value<T>())
+ return -policies::raise_overflow_error<T>(function, 0, pol);
+
+ if(val == 0)
+ return -detail::get_smallest_value<T>();
+
+ T remain = frexp(val, &expon);
+ if(remain == 0.5)
+ --expon; // when val is a power of two we must reduce the exponent
+ T diff = ldexp(T(1), expon - tools::digits<T>());
+ if(diff == 0)
+ diff = detail::get_smallest_value<T>();
+ return val - diff;
+}
+
+#ifdef BOOST_MSVC
+template <class Policy>
+inline double next_less(const double& val, const Policy& pol)
+{
+ static const char* function = "next_less<%1%>(%1%)";
+
+ if(!(boost::math::isfinite)(val))
+ return policies::raise_domain_error<double>(
+ function,
+ "Argument must be finite, but got %1%", val, pol);
+
+ if(val <= -tools::max_value<double>())
+ return -policies::raise_overflow_error<double>(function, 0, pol);
+
+ return ::_nextafter(val, -tools::max_value<double>());
+}
+#endif
+
+template <class T>
+inline T next_less(const T& val)
+{
+ return next_less(val, policies::policy<>());
+}
+
+template <class T, class Policy>
+inline T nextafter(const T& val, const T& direction, const Policy& pol)
+{
+ return val < direction ? boost::math::next_greater(val, pol) : val == direction ? val : boost::math::next_less(val, pol);
+}
+
+template <class T>
+inline T nextafter(const T& val, const T& direction)
+{
+ return nextafter(val, direction, policies::policy<>());
+}
+
+template <class T, class Policy>
+T edit_distance(const T& a, const T& b, const Policy& pol)
+{
+ BOOST_MATH_STD_USING
+ //
+ // Error handling:
+ //
+ static const char* function = "edit_distance<%1%>(%1%, %1%)";
+ if(!(boost::math::isfinite)(a))
+ return policies::raise_domain_error<T>(
+ function,
+ "Argument a must be finite, but got %1%", a, pol);
+ if(!(boost::math::isfinite)(b))
+ return policies::raise_domain_error<T>(
+ function,
+ "Argument b must be finite, but got %1%", b, pol);
+ //
+ // Special cases:
+ //
+ if(a == b)
+ return 0;
+ if(a == 0)
+ return 1 + edit_distance(boost::math::sign(b) * detail::get_smallest_value<T>(), b, pol);
+ if(b == 0)
+ return 1 + edit_distance(boost::math::sign(a) * detail::get_smallest_value<T>(), a, pol);
+ if(boost::math::sign(a) != boost::math::sign(b))
+ return 2 + edit_distance(boost::math::sign(b) * detail::get_smallest_value<T>(), b, pol)
+ + edit_distance(boost::math::sign(a) * detail::get_smallest_value<T>(), a, pol);
+
+ 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 edit_distance(a, split, pol) + edit_distance(split, b, pol);
+ }
+
+ 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
+ // 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);
+ expon = tools::digits<T>() - expon;
+ //
+ // 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:
+ //
+ T mb = -b;
+ T x = a + mb;
+ T z = x - a;
+ T y = (a - (x - z)) + (mb - z);
+ if(x < 0)
+ {
+ x = -x;
+ y = -y;
+ }
+
+ T result = ldexp(x, expon) + ldexp(y, expon);
+ //
+ // Result must be an integer:
+ //
+ BOOST_ASSERT(result == floor(result));
+ return result;
+}
+
+template <class T>
+T edit_distance(const T& a, const T& b)
+{
+ return boost::math::edit_distance(a, b, policies::policy<>());
+}
+
+}} // namespaces
+
+#endif // BOOST_MATH_SPECIAL_NEXT_HPP
+
Modified: sandbox/math_toolkit/libs/math/test/Jamfile.v2
==============================================================================
--- sandbox/math_toolkit/libs/math/test/Jamfile.v2 (original)
+++ sandbox/math_toolkit/libs/math/test/Jamfile.v2 2008-04-29 08:01:22 EDT (Tue, 29 Apr 2008)
@@ -246,6 +246,7 @@
: # requirements
<define>TEST_REAL_CONCEPT
: test_negative_binomial_real_concept ;
+run test_next.cpp ;
run test_nc_chi_squared.cpp
: # command line
: # input files
@@ -444,6 +445,7 @@
compile compile_test/sf_log1p_incl_test.cpp ;
compile compile_test/sf_math_fwd_incl_test.cpp ;
compile compile_test/sf_modf_incl_test.cpp ;
+compile compile_test/sf_next_incl_test.cpp ;
compile compile_test/sf_powm1_incl_test.cpp ;
compile compile_test/sf_round_incl_test.cpp ;
compile compile_test/sf_sign_incl_test.cpp ;
@@ -472,3 +474,4 @@
compile compile_test/tools_test_inc_test.cpp ;
compile compile_test/tools_toms748_inc_test.cpp ;
+
Modified: sandbox/math_toolkit/libs/math/test/compile_test/instantiate.hpp
==============================================================================
--- sandbox/math_toolkit/libs/math/test/compile_test/instantiate.hpp (original)
+++ sandbox/math_toolkit/libs/math/test/compile_test/instantiate.hpp 2008-04-29 08:01:22 EDT (Tue, 29 Apr 2008)
@@ -246,6 +246,10 @@
boost::math::modf(v1, &ll);
#endif
boost::math::pow<2>(v1);
+ boost::math::nextafter(v1, v1);
+ boost::math::next_greater(v1);
+ boost::math::next_less(v1);
+ boost::math::edit_distance(v1, v1);
//
// All over again, with a policy this time:
//
@@ -368,6 +372,10 @@
modf(v1, &ll, pol);
#endif
boost::math::pow<2>(v1, pol);
+ boost::math::nextafter(v1, v1, pol);
+ boost::math::next_greater(v1, pol);
+ boost::math::next_less(v1, pol);
+ boost::math::edit_distance(v1, v1, pol);
//
// All over again with the versions in test::
//
@@ -483,6 +491,10 @@
test::modf(v1, &ll);
#endif
test::pow<2>(v1);
+ test::nextafter(v1, v1);
+ test::next_greater(v1);
+ test::next_less(v1);
+ test::edit_distance(v1, v1);
}
template <class RealType>
Added: sandbox/math_toolkit/libs/math/test/compile_test/sf_next_incl_test.cpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/libs/math/test/compile_test/sf_next_incl_test.cpp 2008-04-29 08:01:22 EDT (Tue, 29 Apr 2008)
@@ -0,0 +1,42 @@
+// Copyright John Maddock 2006.
+// 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)
+//
+// Basic sanity check that header <boost/math/special_functions/next.hpp>
+// #includes all the files that it needs to.
+//
+#include <boost/math/special_functions/next.hpp>
+//
+// Note this header includes no other headers, this is
+// important if this test is to be meaningful:
+//
+#include "test_compile_result.hpp"
+
+void check()
+{
+ check_result<float>(boost::math::nextafter<float>(f, f));
+ check_result<double>(boost::math::nextafter<double>(d, d));
+#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+ check_result<long double>(boost::math::nextafter<long double>(l, l));
+#endif
+
+ check_result<float>(boost::math::next_greater<float>(f));
+ check_result<double>(boost::math::next_greater<double>(d));
+#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+ check_result<long double>(boost::math::next_greater<long double>(l));
+#endif
+
+ check_result<float>(boost::math::next_less<float>(f));
+ check_result<double>(boost::math::next_less<double>(d));
+#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+ check_result<long double>(boost::math::next_less<long double>(l));
+#endif
+
+ check_result<float>(boost::math::edit_distance<float>(f, f));
+ check_result<double>(boost::math::edit_distance<double>(d, d));
+#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+ check_result<long double>(boost::math::edit_distance<long double>(l, l));
+#endif
+
+}
Added: sandbox/math_toolkit/libs/math/test/test_next.cpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/libs/math/test/test_next.cpp 2008-04-29 08:01:22 EDT (Tue, 29 Apr 2008)
@@ -0,0 +1,102 @@
+// (C) Copyright John Maddock 2008.
+// 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)
+
+#include <boost/math/concepts/real_concept.hpp>
+#include <boost/test/included/test_exec_monitor.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+#include <boost/math/special_functions/next.hpp>
+
+template <class T>
+void test_value(const T& val, const char* name)
+{
+ using namespace boost::math;
+ T upper = tools::max_value<T>();
+ T lower = -upper;
+
+ std::cout << "Testing type " << name << " with initial value " << val << std::endl;
+
+ BOOST_CHECK_EQUAL(edit_distance(next_greater(val), val), 1);
+ BOOST_CHECK(next_greater(val) > val);
+ BOOST_CHECK_EQUAL(edit_distance(next_less(val), val), 1);
+ BOOST_CHECK(next_less(val) < val);
+ BOOST_CHECK_EQUAL(edit_distance(nextafter(val, upper), val), 1);
+ BOOST_CHECK(nextafter(val, upper) > val);
+ BOOST_CHECK_EQUAL(edit_distance(nextafter(val, lower), val), 1);
+ BOOST_CHECK(nextafter(val, lower) < val);
+ BOOST_CHECK_EQUAL(edit_distance(next_greater(next_greater(val)), val), 2);
+ BOOST_CHECK_EQUAL(edit_distance(next_less(next_less(val)), val), 2);
+ BOOST_CHECK_EQUAL(edit_distance(next_less(next_greater(val)), val), 0);
+ BOOST_CHECK_EQUAL(edit_distance(next_greater(next_less(val)), val), 0);
+ BOOST_CHECK_EQUAL(next_less(next_greater(val)), val);
+ BOOST_CHECK_EQUAL(next_greater(next_less(val)), val);
+}
+
+template <class T>
+void test_values(const T& val, const char* name)
+{
+ static const T a = static_cast<T>(1.3456724e22);
+ static const T b = static_cast<T>(1.3456724e-22);
+ static const T z = 0;
+ static const T one = 1;
+ static const T two = 2;
+
+ test_value(a, name);
+ test_value(-a, name);
+ test_value(b, name);
+ test_value(-b, name);
+ test_value(boost::math::tools::epsilon<T>(), name);
+ test_value(-boost::math::tools::epsilon<T>(), name);
+ test_value(boost::math::tools::min_value<T>(), name);
+ test_value(-boost::math::tools::min_value<T>(), name);
+ test_value(z, name);
+ test_value(-z, name);
+ test_value(one, name);
+ test_value(-one, name);
+ test_value(two, name);
+ test_value(-two, name);
+ if(std::numeric_limits<T>::is_specialized && std::numeric_limits<T>::has_denorm)
+ {
+ test_value(std::numeric_limits<T>::denorm_min(), name);
+ test_value(-std::numeric_limits<T>::denorm_min(), name);
+ test_value(2 * std::numeric_limits<T>::denorm_min(), name);
+ test_value(-2 * std::numeric_limits<T>::denorm_min(), name);
+ }
+
+ static const unsigned primes[] = {
+ 11, 13, 17, 19, 23, 29,
+ 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
+ 73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
+ 127, 131, 137, 139, 149, 151, 157, 163, 167, 173,
+ 179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
+ 233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
+ 283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
+ 353, 359, 367, 373, 379, 383, 389, 397, 401, 409,
+ 419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
+ };
+
+ for(unsigned i = 0; i < sizeof(primes)/sizeof(primes[0]); ++i)
+ {
+ T v1 = val;
+ T v2 = val;
+ for(unsigned j = 0; j < primes[i]; ++j)
+ {
+ v1 = boost::math::next_greater(v1);
+ v2 = boost::math::next_less(v2);
+ }
+ BOOST_CHECK_EQUAL(boost::math::edit_distance(v1, val), primes[i]);
+ BOOST_CHECK_EQUAL(boost::math::edit_distance(v2, val), primes[i]);
+ }
+}
+
+int test_main(int, char* [])
+{
+ std::cout << boost::math::edit_distance(1.0, 0.0) << std::endl;
+ test_values(1.0f, "float");
+ test_values(1.0, "double");
+ test_values(1.0L, "long double");
+ test_values(boost::math::concepts::real_concept(0), "real_concept");
+ return 0;
+}
+