Fix rounding for non-representable numbers

include/boost/math/ccmath/abs.hpp

@@ -31,7 +31,7 @@ namespace detail {
 template <typename T> 
 constexpr T abs_impl(T x) noexcept
 {
-    if (boost::math::ccmath::isnan(x))
+    if ((boost::math::ccmath::isnan)(x))
     {
         return std::numeric_limits<T>::quiet_NaN();
     }

include/boost/math/ccmath/fpclassify.hpp

@@ -10,6 +10,7 @@
 #include <limits>
 #include <type_traits>
 #include <boost/math/tools/is_constant_evaluated.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
 #include <boost/math/ccmath/abs.hpp>
 #include <boost/math/ccmath/isinf.hpp>
 #include <boost/math/ccmath/isnan.hpp>
@@ -18,19 +19,19 @@
 namespace boost::math::ccmath {
 
 template <typename T, std::enable_if_t<!std::is_integral_v<T>, bool> = true>
-inline constexpr int fpclassify(T x)
+inline constexpr int fpclassify BOOST_PREVENT_MACRO_SUBSTITUTION(T x)
 {
     if(BOOST_MATH_IS_CONSTANT_EVALUATED(x))
     {
-        return boost::math::ccmath::isnan(x) ? FP_NAN :
-               boost::math::ccmath::isinf(x) ? FP_INFINITE :
+        return (boost::math::ccmath::isnan)(x) ? FP_NAN :
+               (boost::math::ccmath::isinf)(x) ? FP_INFINITE :
                boost::math::ccmath::abs(x) == T(0) ? FP_ZERO :
                boost::math::ccmath::abs(x) > 0 && boost::math::ccmath::abs(x) < (std::numeric_limits<T>::min)() ? FP_SUBNORMAL : FP_NORMAL;
     }
     else
     {
-        using std::fpclassify;
-        return fpclassify(x);
+        using boost::math::fpclassify;
+        return (fpclassify)(x);
     }
 }
 

include/boost/math/ccmath/isinf.hpp

@@ -10,6 +10,7 @@
 #include <limits>
 #include <type_traits>
 #include <boost/math/tools/is_constant_evaluated.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
 
 #include <boost/math/tools/is_standalone.hpp>
 #ifndef BOOST_MATH_STANDALONE
@@ -22,7 +23,7 @@
 namespace boost::math::ccmath {
 
 template <typename T>
-constexpr bool isinf(T x) noexcept
+constexpr bool isinf BOOST_PREVENT_MACRO_SUBSTITUTION(T x) noexcept
 {
     if(BOOST_MATH_IS_CONSTANT_EVALUATED(x))
     {
@@ -37,15 +38,15 @@ constexpr bool isinf(T x) noexcept
     }
     else
     {
-        using std::isinf;
+        using boost::math::isinf;
 
         if constexpr (!std::is_integral_v<T>)
         {
-            return isinf(x);
+            return (isinf)(x);
         }
         else
         {
-            return isinf(static_cast<double>(x));
+            return (isinf)(static_cast<double>(x));
         }
     }
 }

include/boost/math/ccmath/isnan.hpp

@@ -9,6 +9,7 @@
 #include <cmath>
 #include <type_traits>
 #include <boost/math/tools/is_constant_evaluated.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
 
 #include <boost/math/tools/is_standalone.hpp>
 #ifndef BOOST_MATH_STANDALONE
@@ -21,23 +22,23 @@
 namespace boost::math::ccmath {
 
 template <typename T>
-inline constexpr bool isnan(T x)
+inline constexpr bool isnan BOOST_PREVENT_MACRO_SUBSTITUTION(T x)
 {
     if(BOOST_MATH_IS_CONSTANT_EVALUATED(x))
     {
         return x != x;
     }
     else
     {
-        using std::isnan;
+        using boost::math::isnan;
 
         if constexpr (!std::is_integral_v<T>)
         {
-            return isnan(x);
+            return (isnan)(x);
         }
         else
         {
-            return isnan(static_cast<double>(x));
+            return (isnan)(static_cast<double>(x));
         }
     }
 }

include/boost/math/ccmath/ldexp.hpp

@@ -45,8 +45,8 @@ inline constexpr Real ldexp(Real arg, int exp) noexcept
     if(BOOST_MATH_IS_CONSTANT_EVALUATED(arg))
     {
         return boost::math::ccmath::abs(arg) == Real(0) ? arg :
-               boost::math::ccmath::isinf(arg) ? arg :
-               boost::math::ccmath::isnan(arg) ? arg :
+               (boost::math::ccmath::isinf)(arg) ? arg :
+               (boost::math::ccmath::isnan)(arg) ? arg :
                boost::math::ccmath::detail::ldexp_impl(arg, exp);
     }
     else

include/boost/math/special_functions/round.hpp

@@ -1,4 +1,5 @@
 //  Copyright John Maddock 2007.
+//  Copyright Matt Borland 2023.
 //  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)
@@ -14,6 +15,16 @@
 #include <boost/math/policies/error_handling.hpp>
 #include <boost/math/special_functions/math_fwd.hpp>
 #include <boost/math/special_functions/fpclassify.hpp>
+#include <type_traits>
+#include <limits>
+#include <cmath>
+
+#ifndef BOOST_NO_CXX17_IF_CONSTEXPR
+#include <boost/math/ccmath/ldexp.hpp>
+#  if !defined(BOOST_MATH_NO_CONSTEXPR_DETECTION)
+#    define BOOST_MATH_HAS_CONSTEXPR_LDEXP
+#  endif
+#endif
 
 namespace boost{ namespace math{
 
@@ -91,11 +102,40 @@ inline int iround(const T& v, const Policy& pol)
    BOOST_MATH_STD_USING
    using result_type = tools::promote_args_t<T>;
 
-   T r = boost::math::round(v, pol);
-   if(r > static_cast<result_type>((std::numeric_limits<int>::max)()) || r < static_cast<result_type>((std::numeric_limits<int>::min)()))
+   result_type r = boost::math::round(v, pol);
+
+   #ifdef BOOST_MATH_HAS_CONSTEXPR_LDEXP
+   if constexpr (std::is_arithmetic_v<result_type>
+                 #ifdef BOOST_MATH_FLOAT128_TYPE
+                 && !std::is_same_v<BOOST_MATH_FLOAT128_TYPE, result_type>
+                 #endif
+                )
+   {
+      constexpr result_type max_val = boost::math::ccmath::ldexp(static_cast<result_type>(1), std::numeric_limits<int>::digits);
+
+      if (r >= max_val || r < -max_val)
+      {
+         return static_cast<int>(boost::math::policies::raise_rounding_error("boost::math::iround<%1%>(%1%)", nullptr, v, static_cast<int>(0), pol));
+      }
+   }
+   else
+   {
+      static const result_type max_val = ldexp(static_cast<result_type>(1), std::numeric_limits<int>::digits);
+
+      if (r >= max_val || r < -max_val)
+      {
+         return static_cast<int>(boost::math::policies::raise_rounding_error("boost::math::iround<%1%>(%1%)", nullptr, v, static_cast<int>(0), pol));
+      }
+   }
+   #else
+   static const result_type max_val = ldexp(static_cast<result_type>(1), std::numeric_limits<int>::digits);
+
+   if (r >= max_val || r < -max_val)
    {
-      return static_cast<int>(policies::raise_rounding_error("boost::math::iround<%1%>(%1%)", nullptr, v, 0, pol));
+      return static_cast<int>(boost::math::policies::raise_rounding_error("boost::math::iround<%1%>(%1%)", nullptr, v, static_cast<int>(0), pol));
    }
+   #endif
+
    return static_cast<int>(r);
 }
 template <class T>
@@ -109,12 +149,42 @@ inline long lround(const T& v, const Policy& pol)
 {
    BOOST_MATH_STD_USING
    using result_type = tools::promote_args_t<T>;
-   T r = boost::math::round(v, pol);
-   if(r > static_cast<result_type>((std::numeric_limits<long>::max)()) || r < static_cast<result_type>((std::numeric_limits<long>::min)()))
+
+   result_type r = boost::math::round(v, pol);
+
+   #ifdef BOOST_MATH_HAS_CONSTEXPR_LDEXP
+   if constexpr (std::is_arithmetic_v<result_type>
+                 #ifdef BOOST_MATH_FLOAT128_TYPE
+                 && !std::is_same_v<BOOST_MATH_FLOAT128_TYPE, result_type>
+                 #endif
+                )
+   {
+      constexpr result_type max_val = boost::math::ccmath::ldexp(static_cast<result_type>(1), std::numeric_limits<long>::digits);
+
+      if (r >= max_val || r < -max_val)
+      {
+         return static_cast<long>(boost::math::policies::raise_rounding_error("boost::math::lround<%1%>(%1%)", nullptr, v, static_cast<long>(0), pol));
+      }
+   }
+   else
    {
-      return static_cast<long int>(policies::raise_rounding_error("boost::math::lround<%1%>(%1%)", nullptr, v, 0L, pol));
+      static const result_type max_val = ldexp(static_cast<result_type>(1), std::numeric_limits<long>::digits);
+
+      if (r >= max_val || r < -max_val)
+      {
+         return static_cast<long>(boost::math::policies::raise_rounding_error("boost::math::lround<%1%>(%1%)", nullptr, v, static_cast<long>(0), pol));
+      }
    }
-   return static_cast<long int>(r);
+   #else
+   static const result_type max_val = ldexp(static_cast<result_type>(1), std::numeric_limits<long>::digits);
+
+   if (r >= max_val || r < -max_val)
+   {
+      return static_cast<long>(boost::math::policies::raise_rounding_error("boost::math::lround<%1%>(%1%)", nullptr, v, static_cast<long>(0), pol));
+   }
+   #endif
+
+   return static_cast<long>(r);
 }
 template <class T>
 inline long lround(const T& v)
@@ -126,14 +196,42 @@ template <class T, class Policy>
 inline long long llround(const T& v, const Policy& pol)
 {
    BOOST_MATH_STD_USING
-   using result_type = tools::promote_args_t<T>;
+   using result_type = boost::math::tools::promote_args_t<T>;
+
+   result_type r = boost::math::round(v, pol);
 
-   T r = boost::math::round(v, pol);
-   if(r > static_cast<result_type>((std::numeric_limits<long long>::max)()) ||
-      r < static_cast<result_type>((std::numeric_limits<long long>::min)()))
+   #ifdef BOOST_MATH_HAS_CONSTEXPR_LDEXP
+   if constexpr (std::is_arithmetic_v<result_type>
+                 #ifdef BOOST_MATH_FLOAT128_TYPE
+                 && !std::is_same_v<BOOST_MATH_FLOAT128_TYPE, result_type>
+                 #endif
+                )
    {
-      return static_cast<long long>(policies::raise_rounding_error("boost::math::llround<%1%>(%1%)", nullptr, v, static_cast<long long>(0), pol));
+      constexpr result_type max_val = boost::math::ccmath::ldexp(static_cast<result_type>(1), std::numeric_limits<long long>::digits);
+
+      if (r >= max_val || r < -max_val)
+      {
+         return static_cast<long long>(boost::math::policies::raise_rounding_error("boost::math::llround<%1%>(%1%)", nullptr, v, static_cast<long long>(0), pol));
+      }
    }
+   else
+   {
+      static const result_type max_val = ldexp(static_cast<result_type>(1), std::numeric_limits<long long>::digits);
+
+      if (r >= max_val || r < -max_val)
+      {
+         return static_cast<long long>(boost::math::policies::raise_rounding_error("boost::math::llround<%1%>(%1%)", nullptr, v, static_cast<long long>(0), pol));
+      }
+   }
+   #else
+   static const result_type max_val = ldexp(static_cast<result_type>(1), std::numeric_limits<long long>::digits);
+
+   if (r >= max_val || r < -max_val)
+   {
+      return static_cast<long long>(boost::math::policies::raise_rounding_error("boost::math::llround<%1%>(%1%)", nullptr, v, static_cast<long long>(0), pol));
+   }
+   #endif
+
    return static_cast<long long>(r);
 }
 template <class T>