From 08d2c34e66104891c9be328e5c17c2e7e928d656 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Sun, 21 Aug 2022 19:54:55 -0700 Subject: [PATCH] Offer 128 bit support --- .../special_functions/fast_float_distance.hpp | 81 +++++++++++++++++++ include/boost/math/special_functions/next.hpp | 5 ++ test/Jamfile.v2 | 1 + test/test_fast_float_distance.cpp | 44 ++++++++++ 4 files changed, 131 insertions(+) create mode 100644 include/boost/math/special_functions/fast_float_distance.hpp create mode 100644 test/test_fast_float_distance.cpp diff --git a/include/boost/math/special_functions/fast_float_distance.hpp b/include/boost/math/special_functions/fast_float_distance.hpp new file mode 100644 index 0000000000..ad07a1d3f0 --- /dev/null +++ b/include/boost/math/special_functions/fast_float_distance.hpp @@ -0,0 +1,81 @@ +// (C) Copyright Matt Borland 2022. +// 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_SF_FAST_FLOAT_DISTANCE +#define BOOST_MATH_SF_FAST_FLOAT_DISTANCE + +#include +#include +#include + +#if defined(BOOST_MATH_USE_FLOAT128) && !defined(BOOST_MATH_STANDALONE) +#include +#include +#define BOOST_MATH_USE_FAST_FLOAT128 +#endif + +namespace boost { namespace math { + +// https://randomascii.wordpress.com/2012/01/23/stupid-float-tricks-2/ +// https://blog.regehr.org/archives/959 +inline std::int32_t fast_float_distance(float a, float b) +{ + return boost::math::float_distance(a, b); +} + +inline std::int64_t fast_float_distance(double a, double b) +{ + return boost::math::float_distance(a, b); +} + +#ifdef BOOST_MATH_USE_FAST_FLOAT128 +boost::multiprecision::int128_type fast_float_distance(boost::multiprecision::float128_type a, boost::multiprecision::float128_type b) +{ + using std::abs; + using std::isfinite; + + constexpr boost::multiprecision::float128_type tol = BOOST_MP_QUAD_MIN; + + // 0, very small, and large magnitude distances all need special handling + if (abs(a) == 0 || abs(b) == 0) + { + return 0; + } + else if (abs(a) < tol || abs(b) < tol) + { + BOOST_MATH_THROW_EXCEPTION(std::domain_error("special handling is required for tiny distances. Please use boost::math::float_distance for a slower but safe solution")); + } + + if (!(isfinite)(a)) + { + BOOST_MATH_THROW_EXCEPTION(std::domain_error("Both arguments to fast_float_distnace must be finite")); + } + else if (!(isfinite)(b)) + { + BOOST_MATH_THROW_EXCEPTION(std::domain_error("Both arguments to fast_float_distnace must be finite")); + } + + static_assert(sizeof(boost::multiprecision::int128_type) == sizeof(boost::multiprecision::float128_type)); + + boost::multiprecision::int128_type ai; + boost::multiprecision::int128_type bi; + std::memcpy(&ai, &a, sizeof(boost::multiprecision::float128_type)); + std::memcpy(&bi, &b, sizeof(boost::multiprecision::float128_type)); + + boost::multiprecision::int128_type result = bi - ai; + + if (ai < 0 || bi < 0) + { + result = -result; + } + + return result; +} + +#endif + +}} // Namespaces + +#endif diff --git a/include/boost/math/special_functions/next.hpp b/include/boost/math/special_functions/next.hpp index 4336c07e0e..91f58f594e 100644 --- a/include/boost/math/special_functions/next.hpp +++ b/include/boost/math/special_functions/next.hpp @@ -22,6 +22,11 @@ #include #include +#if defined(BOOST_MATH_USE_FLOAT128) && !defined(BOOST_MATH_STANDALONE) +#include +#include +#define BOOST_MATH_USE_FAST_FLOAT128 +#endif #if !defined(_CRAYC) && !defined(__CUDACC__) && (!defined(__GNUC__) || (__GNUC__ > 3) || ((__GNUC__ == 3) && (__GNUC_MINOR__ > 3))) #if (defined(_M_IX86_FP) && (_M_IX86_FP >= 2)) || defined(__SSE2__) diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 17aa8484a0..ea1580676d 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -526,6 +526,7 @@ test-suite special_fun : [ run test_ldouble_simple.cpp ../../test/build//boost_unit_test_framework ] # Needs to run in release mode, as it's rather slow: [ run test_next.cpp pch ../../test/build//boost_unit_test_framework : : : release ] + [ run test_fast_float_distance.cpp ../../test/build//boost_unit_test_framework : : : release [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : "-Bstatic -lquadmath -Bdynamic" : no ] ] [ run test_next_decimal.cpp pch ../../test/build//boost_unit_test_framework : : : release ] [ run test_owens_t.cpp ../../test/build//boost_unit_test_framework ] [ run test_polygamma.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ] diff --git a/test/test_fast_float_distance.cpp b/test/test_fast_float_distance.cpp new file mode 100644 index 0000000000..b2820233dc --- /dev/null +++ b/test/test_fast_float_distance.cpp @@ -0,0 +1,44 @@ +// (C) Copyright John Maddock 2008. +// (C) Copyright Matt Borland 2022. +// 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 +#include +#include + +#include "math_unit_test.hpp" + + +template +void test_value(const T& val) +{ + using namespace boost::math; + + assert(fast_float_distance(float_next(val), val) == -1); + assert(float_next(val) > val); + assert(float_next(float_prior(val)) == val); + + assert(fast_float_distance(float_advance(val, 4), val) == -4); + assert(fast_float_distance(float_advance(val, -4), val) == 4); + if(std::numeric_limits::is_specialized && (std::numeric_limits::has_denorm == std::denorm_present)) + { + assert(fast_float_distance(float_advance(float_next(float_next(val)), 4), float_next(float_next(val))) == -4); + assert(fast_float_distance(float_advance(float_next(float_next(val)), -4), float_next(float_next(val))) == 4); + } +} + +int main(void) +{ + test_value(1.0f); + test_value(1.0); + + #ifdef BOOST_MATH_USE_FAST_FLOAT128 + test_value(boost::multiprecision::float128_type(0)); + test_value(__float128(0)); + #endif +}