Offer 128 bit support

Author: mborland

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 <boost/math/special_functions/next.hpp>
+#include <boost/math/tools/throw_exception.hpp>
+#include <stdexcept>
+
+#if defined(BOOST_MATH_USE_FLOAT128) && !defined(BOOST_MATH_STANDALONE)
+#include <boost/multiprecision/float128.hpp>
+#include <boost/multiprecision/detail/standalone_config.hpp>
+#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

include/boost/math/special_functions/next.hpp

@@ -22,6 +22,11 @@
 #include <cstdint>
 #include <cstring>
 
+#if defined(BOOST_MATH_USE_FLOAT128) && !defined(BOOST_MATH_STANDALONE)
+#include <boost/multiprecision/float128.hpp>
+#include <boost/multiprecision/detail/standalone_config.hpp>
+#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__)

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" : <linkflags>"-Bstatic -lquadmath -Bdynamic" : <build>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  ]

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 <boost/math/concepts/real_concept.hpp>
+#include <boost/test/tools/floating_point_comparison.hpp>
+#include <boost/math/special_functions/next.hpp>
+#include <boost/math/special_functions/ulp.hpp>
+#include <boost/math/special_functions/fast_float_distance.hpp>
+#include <boost/multiprecision/float128.hpp>
+
+#include "math_unit_test.hpp"
+
+
+template <class T>
+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<T>::is_specialized && (std::numeric_limits<T>::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
+}