Add 128bit support

Author: mborland

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>
+#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__)
@@ -786,6 +791,39 @@ inline std::int64_t float_distance(double a, double b)
    return result;
 }
 
+#if defined(BOOST_MATH_USE_FLOAT128) && !defined(BOOST_MATH_STANDALONE)
+boost::multiprecision::int128_type float_distance(boost::multiprecision::float128 a, boost::multiprecision::float128 b)
+{
+   using std::abs;
+   constexpr boost::multiprecision::float128 tol = 2 * (std::numeric_limits<boost::multiprecision::float128>::min)();
+
+   if (abs(a) == 0 || abs(b) == 0)
+   {
+      return static_cast<boost::multiprecision::int128_type>(float_distance(a, b, policies::policy<>()));
+   }
+   else if (abs(a) < tol || abs(b) < tol)
+   {
+      return static_cast<boost::multiprecision::int128_type>(float_distance(a, b, policies::policy<>()));
+   }
+
+   static_assert(sizeof(boost::multiprecision::int128_type) == sizeof(boost::multiprecision::float128), "128 bit float is the wrong size");
+
+   boost::multiprecision::int128_type ai;
+   boost::multiprecision::int128_type bi;
+   std::memcpy(&ai, &a, sizeof(boost::multiprecision::float128));
+   std::memcpy(&ai, &a, sizeof(boost::multiprecision::float128));
+
+   boost::multiprecision::int128_type result = bi - ai;
+
+   if (ai < 0 || bi < 0)
+   {
+      result = -result;
+   }
+
+   return result;
+}
+#endif
+
 namespace detail{
 
 template <class T, class Policy>