From 071d338765847a15eb924574b65a5fa2b8e66037 Mon Sep 17 00:00:00 2001
From: Jimmy Lu <gongchuo.lu@gmail.com>
Date: Sun, 30 Jul 2023 22:27:50 -0400
Subject: [PATCH] Fix inverse_discrete_quantile for large guess

If `guess` passed to `inverse_discrete_quantile` cannot be represented
as floating point number, it is possible that `guess + 1` or `guess - 1`
does not change the value at all and we are stuck in infinite loop
inside `round_to_floor` or `round_to_ceil`.  Fix this by
increase/decrease more than 1 in these cases.

Example code to reproduce this:
```c++
boost::math::binomial_distribution<> dist(9079765771874083840, 0.561815);
boost::math::quantile(dist, 0.0365346);
```
---
 .../detail/inv_discrete_quantile.hpp          | 24 +++++++++++--------
 test/test_binomial.cpp                        | 13 ++++++++++
 2 files changed, 27 insertions(+), 10 deletions(-)

diff --git a/include/boost/math/distributions/detail/inv_discrete_quantile.hpp b/include/boost/math/distributions/detail/inv_discrete_quantile.hpp
index 25216e88f6..444f2e8eeb 100644
--- a/include/boost/math/distributions/detail/inv_discrete_quantile.hpp
+++ b/include/boost/math/distributions/detail/inv_discrete_quantile.hpp
@@ -302,15 +302,17 @@ inline typename Dist::value_type round_to_floor(const Dist& d, typename Dist::va
    //
    while(result != 0)
    {
-      cc = result - 1;
+      cc = result;
+      for (typename Dist::value_type d = 1; cc == result; d += 1)
+      {
+         cc = result - d;
+      }
       if(cc < support(d).first)
          break;
       pp = c ? cdf(complement(d, cc)) : cdf(d, cc);
-      if(pp == p)
-         result = cc;
-      else if(c ? pp > p : pp < p)
+      if(c ? pp > p : pp < p)
          break;
-      result -= 1;
+      result = cc;
    }
 
    return result;
@@ -336,15 +338,17 @@ inline typename Dist::value_type round_to_ceil(const Dist& d, typename Dist::val
    //
    while(true)
    {
-      cc = result + 1;
+      cc = result;
+      for (typename Dist::value_type d = 1; cc == result; d += 1)
+      {
+         cc = result + d;
+      }
       if(cc > support(d).second)
          break;
       pp = c ? cdf(complement(d, cc)) : cdf(d, cc);
-      if(pp == p)
-         result = cc;
-      else if(c ? pp < p : pp > p)
+      if(c ? pp < p : pp > p)
          break;
-      result += 1;
+      result = cc;
    }
 
    return result;
diff --git a/test/test_binomial.cpp b/test/test_binomial.cpp
index fa50da7f5a..ced7bc2edd 100644
--- a/test/test_binomial.cpp
+++ b/test/test_binomial.cpp
@@ -716,6 +716,19 @@ void test_spots(RealType T)
 
    check_out_of_range<boost::math::binomial_distribution<RealType> >(1, 1); // (All) valid constructor parameter values.
 
+  {
+     binomial_distribution<RealType> dist(9079765771874083840, 0.561815);
+#ifdef TEST_FLOAT
+     RealType expected = 5101149012694663168;
+#endif
+#ifdef TEST_DOUBLE
+     RealType expected = 5101148604445670400;
+#endif
+#ifdef TEST_LDOUBLE
+     RealType expected = 5101148604445669693;
+#endif
+     BOOST_CHECK_EQUAL(quantile(dist, 0.0365346), expected);
+  }
 
 } // template <class RealType>void test_spots(RealType)