From 73b8ffc13c90ab95fdd51f2bbc0bdd15a82ed21d Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Thu, 16 Mar 2023 19:26:23 -0700
Subject: [PATCH] Fix llround for non-representable numbers

---
 .../boost/math/special_functions/round.hpp    | 44 +++++++++++++++++--
 test/Jamfile.v2                               |  1 +
 test/git_issue_430.cpp                        | 17 +++++++
 3 files changed, 59 insertions(+), 3 deletions(-)
 create mode 100644 test/git_issue_430.cpp

diff --git a/include/boost/math/special_functions/round.hpp b/include/boost/math/special_functions/round.hpp
index d46bd4480..79afcb699 100644
--- a/include/boost/math/special_functions/round.hpp
+++ b/include/boost/math/special_functions/round.hpp
@@ -14,11 +14,40 @@
 #include <boost/math/policies/error_handling.hpp>
 #include <boost/math/special_functions/math_fwd.hpp>
 #include <boost/math/special_functions/fpclassify.hpp>
+#include <limits>
+#include <cmath>
 
 namespace boost{ namespace math{
 
 namespace detail{
 
+// https://stackoverflow.com/questions/8905246/how-to-check-if-float-can-be-exactly-represented-as-an-integer/17822304#17822304
+template <typename T, typename ResultType>
+inline ResultType float_to_int(T x)
+{
+   BOOST_MATH_STD_USING
+
+   T y = floor(x);
+   if (y < static_cast<T>(0.0L))
+   {
+      return 0;
+   }
+   
+   if (y > ldexp(static_cast<T>(1.0L), sizeof(ResultType) * CHAR_BIT))
+   {
+      return (std::numeric_limits<ResultType>::max)();
+   }
+
+   return y;
+}
+
+template <typename T, typename TargetType>
+inline bool is_representable(T x)
+{
+   BOOST_MATH_STD_USING
+   return (floor(x) == x && x >= static_cast<T>(0.0L) && x < ldexp(1.0, sizeof(TargetType) * CHAR_BIT));
+}
+
 template <class T, class Policy>
 inline tools::promote_args_t<T> round(const T& v, const Policy& pol, const std::false_type&)
 {
@@ -129,12 +158,21 @@ inline long long llround(const T& v, const Policy& pol)
    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 long>::max)()) ||
-      r < static_cast<result_type>((std::numeric_limits<long long>::min)()))
+   long long return_val = boost::math::detail::float_to_int<result_type, long long>(r);
+   bool representable = boost::math::detail::is_representable<result_type, long long>(r);
+
+   if (r < static_cast<result_type>((std::numeric_limits<long long>::min)()) ||
+       (return_val == LLONG_MAX && !representable))
    {
       return static_cast<long long>(policies::raise_rounding_error("boost::math::llround<%1%>(%1%)", nullptr, v, static_cast<long long>(0), pol));
    }
-   return static_cast<long long>(r);
+
+   if (r < 0)
+   {
+      return_val = static_cast<long long>(r);
+   }
+
+   return return_val;
 }
 template <class T>
 inline long long llround(const T& v)
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index d412ca10e..9a3ea8ae2 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -890,6 +890,7 @@ test-suite distribution_tests :
    [ run complex_test.cpp ../../test/build//boost_unit_test_framework  ]
 
    [ compile test_dist_deduction_guides.cpp : [ requires cpp_deduction_guides cpp_variadic_templates ] ]
+   [ run git_issue_430.cpp ../../test/build//boost_unit_test_framework ]
    [ run git_issue_800.cpp ../../test/build//boost_unit_test_framework  ]
    [ run git_issue_845.cpp ../../test/build//boost_unit_test_framework  ]
    [ run scipy_issue_14901.cpp ../../test/build//boost_unit_test_framework ]
diff --git a/test/git_issue_430.cpp b/test/git_issue_430.cpp
new file mode 100644
index 000000000..b23ffc587
--- /dev/null
+++ b/test/git_issue_430.cpp
@@ -0,0 +1,17 @@
+// 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)
+
+#include <boost/math/special_functions/round.hpp>
+#include <cstdint>
+#include "math_unit_test.hpp"
+
+double x = 9223372036854775807.0;  // can't be represented as double, will have a different value at runtime. 
+int main()
+{
+    int64_t result = boost::math::llround(x);
+    CHECK_EQUAL(result, INT64_C(9223372036854775807));
+
+    return boost::math::test::report_errors();
+}