From d7e23fd2a413bea8ac2296a7f795dbc512aba59b Mon Sep 17 00:00:00 2001 From: John Maddock Date: Tue, 29 Apr 2008 12:01:22 +0000 Subject: [PATCH] Added first cut of nextafter family of functions. [SVN r44878] --- include/boost/math/special_functions.hpp | 1 + .../boost/math/special_functions/math_fwd.hpp | 23 ++ include/boost/math/special_functions/next.hpp | 240 ++++++++++++++++++ test/Jamfile.v2 | 3 + test/compile_test/instantiate.hpp | 12 + test/compile_test/sf_next_incl_test.cpp | 42 +++ test/test_next.cpp | 102 ++++++++ 7 files changed, 423 insertions(+) create mode 100644 include/boost/math/special_functions/next.hpp create mode 100644 test/compile_test/sf_next_incl_test.cpp create mode 100644 test/test_next.cpp diff --git a/include/boost/math/special_functions.hpp b/include/boost/math/special_functions.hpp index 44ea43b8c..bfb26a77b 100644 --- a/include/boost/math/special_functions.hpp +++ b/include/boost/math/special_functions.hpp @@ -41,6 +41,7 @@ #include #include #include +#include #include #include #include diff --git a/include/boost/math/special_functions/math_fwd.hpp b/include/boost/math/special_functions/math_fwd.hpp index 947b182b3..d95057e0e 100644 --- a/include/boost/math/special_functions/math_fwd.hpp +++ b/include/boost/math/special_functions/math_fwd.hpp @@ -657,6 +657,24 @@ namespace boost template typename tools::promote_args::type zeta(T s); + // next: + template + T nextafter(const T&, const T&, const Policy&); + template + T nextafter(const T&, const T&); + template + T next_greater(const T&, const Policy&); + template + T next_greater(const T&); + template + T next_less(const T&, const Policy&); + template + T next_less(const T&); + template + T edit_distance(const T&, const T&, const Policy&); + template + T edit_distance(const T&, const T&); + } // namespace math } // namespace boost @@ -1003,6 +1021,11 @@ namespace boost \ template \ inline typename boost::math::tools::promote_args::type pow(T v){ return boost::math::pow(v, Policy()); }\ + template T nextafter(const T& a, const T& b){ return boost::math::nextafter(a, b, Policy()); }\ + template T next_greater(const T& a){ return boost::math::next_greater(a, Policy()); }\ + template T next_less(const T& a){ return boost::math::next_less(a, Policy()); }\ + template T edit_distance(const T& a, const T& b){ return boost::math::edit_distance(a, b, Policy()); }\ + #endif // BOOST_MATH_SPECIAL_MATH_FWD_HPP diff --git a/include/boost/math/special_functions/next.hpp b/include/boost/math/special_functions/next.hpp new file mode 100644 index 000000000..87b4aac3c --- /dev/null +++ b/include/boost/math/special_functions/next.hpp @@ -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 +#include +#include + +#ifdef BOOST_MSVC +#include +#endif + +namespace boost{ namespace math{ + +namespace detail{ + +template +inline T get_smallest_value(mpl::true_ const&) +{ + return std::numeric_limits::denorm_min(); +} + +template +inline T get_smallest_value(mpl::false_ const&) +{ + return tools::min_value(); +} + +template +inline T get_smallest_value() +{ + return get_smallest_value(mpl::bool_::is_specialized && std::numeric_limits::has_denorm>()); +} + +} + +template +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( + function, + "Argument must be finite, but got %1%", val, pol); + + if(val >= tools::max_value()) + return policies::raise_overflow_error(function, 0, pol); + + if(val == 0) + return detail::get_smallest_value(); + + 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()); + if(diff == 0) + diff = detail::get_smallest_value(); + return val + diff; +} + +#ifdef BOOST_MSVC +template +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( + function, + "Argument must be finite, but got %1%", val, pol); + + if(val >= tools::max_value()) + return policies::raise_overflow_error(function, 0, pol); + + return ::_nextafter(val, tools::max_value()); +} +#endif + +template +inline T next_greater(const T& val) +{ + return next_greater(val, policies::policy<>()); +} + +template +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( + function, + "Argument must be finite, but got %1%", val, pol); + + if(val <= -tools::max_value()) + return -policies::raise_overflow_error(function, 0, pol); + + if(val == 0) + return -detail::get_smallest_value(); + + 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()); + if(diff == 0) + diff = detail::get_smallest_value(); + return val - diff; +} + +#ifdef BOOST_MSVC +template +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( + function, + "Argument must be finite, but got %1%", val, pol); + + if(val <= -tools::max_value()) + return -policies::raise_overflow_error(function, 0, pol); + + return ::_nextafter(val, -tools::max_value()); +} +#endif + +template +inline T next_less(const T& val) +{ + return next_less(val, policies::policy<>()); +} + +template +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 +inline T nextafter(const T& val, const T& direction) +{ + return nextafter(val, direction, policies::policy<>()); +} + +template +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( + function, + "Argument a must be finite, but got %1%", a, pol); + if(!(boost::math::isfinite)(b)) + return policies::raise_domain_error( + 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(), b, pol); + if(b == 0) + return 1 + edit_distance(boost::math::sign(a) * detail::get_smallest_value(), a, pol); + if(boost::math::sign(a) != boost::math::sign(b)) + return 2 + edit_distance(boost::math::sign(b) * detail::get_smallest_value(), b, pol) + + edit_distance(boost::math::sign(a) * detail::get_smallest_value(), a, pol); + + if((std::min)(fabs(a), fabs(b)) / (std::max)(fabs(a), fabs(b)) < 2 * tools::epsilon()) + { + bool biga = fabs(a) > fabs(b); + T split = ldexp(biga ? b : a, tools::digits() - 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() + // significant bits in the representation: + // + frexp((boost::math::fpclassify(mv) == FP_SUBNORMAL) ? tools::min_value() : mv, &expon); + expon = tools::digits() - 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 +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 + diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index e4470e805..88421f05b 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -246,6 +246,7 @@ run test_negative_binomial.cpp : # requirements 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_legendre_incl_test.cpp ; 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_data_inc_test.cpp ; compile compile_test/tools_test_inc_test.cpp ; compile compile_test/tools_toms748_inc_test.cpp ; + diff --git a/test/compile_test/instantiate.hpp b/test/compile_test/instantiate.hpp index 7fad7b4a0..f676d3e5a 100644 --- a/test/compile_test/instantiate.hpp +++ b/test/compile_test/instantiate.hpp @@ -246,6 +246,10 @@ void instantiate(RealType) 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 @@ void instantiate(RealType) 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 @@ void instantiate(RealType) 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 diff --git a/test/compile_test/sf_next_incl_test.cpp b/test/compile_test/sf_next_incl_test.cpp new file mode 100644 index 000000000..feac333ee --- /dev/null +++ b/test/compile_test/sf_next_incl_test.cpp @@ -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 +// #includes all the files that it needs to. +// +#include +// +// 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(boost::math::nextafter(f, f)); + check_result(boost::math::nextafter(d, d)); +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + check_result(boost::math::nextafter(l, l)); +#endif + + check_result(boost::math::next_greater(f)); + check_result(boost::math::next_greater(d)); +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + check_result(boost::math::next_greater(l)); +#endif + + check_result(boost::math::next_less(f)); + check_result(boost::math::next_less(d)); +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + check_result(boost::math::next_less(l)); +#endif + + check_result(boost::math::edit_distance(f, f)); + check_result(boost::math::edit_distance(d, d)); +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + check_result(boost::math::edit_distance(l, l)); +#endif + +} diff --git a/test/test_next.cpp b/test/test_next.cpp new file mode 100644 index 000000000..131b38f71 --- /dev/null +++ b/test/test_next.cpp @@ -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 +#include +#include +#include + +template +void test_value(const T& val, const char* name) +{ + using namespace boost::math; + T upper = tools::max_value(); + 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 +void test_values(const T& val, const char* name) +{ + static const T a = static_cast(1.3456724e22); + static const T b = static_cast(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(), name); + test_value(-boost::math::tools::epsilon(), name); + test_value(boost::math::tools::min_value(), name); + test_value(-boost::math::tools::min_value(), 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::is_specialized && std::numeric_limits::has_denorm) + { + test_value(std::numeric_limits::denorm_min(), name); + test_value(-std::numeric_limits::denorm_min(), name); + test_value(2 * std::numeric_limits::denorm_min(), name); + test_value(-2 * std::numeric_limits::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; +} +