boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Added first cut of nextafter family of functions.

Browse Source 
[SVN r44878]
This commit is contained in:
John Maddock
2008-04-29 12:01:22 +00:00
parent 3ccbd30ab3
commit d7e23fd2a4
7 changed files with 423 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
include/boost/math/special_functions.hpp
View File

@@ -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>

23
include/boost/math/special_functions/math_fwd.hpp
View File

@@ -657,6 +657,24 @@ namespace boost
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 @@ namespace boost
\
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

240
include/boost/math/special_functions/next.hpp Normal file
View File

@@ -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

3
test/Jamfile.v2
View File

@@ -246,6 +246,7 @@ run test_negative_binomial.cpp
: # 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_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 ;

12
test/compile_test/instantiate.hpp
View File

@@ -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 <class RealType>

42
test/compile_test/sf_next_incl_test.cpp Normal file
View File

@@ -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
}

102
test/test_next.cpp Normal file
View File

@@ -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;
}