Skip to content

Commit

Permalink
Fix inverse_discrete_quantile for large guess
Browse files Browse the repository at this point in the history
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);
```
  • Loading branch information
Yuhta committed Jul 31, 2023
1 parent 281f491 commit 071d338
Show file tree
Hide file tree
Showing 2 changed files with 27 additions and 10 deletions.
24 changes: 14 additions & 10 deletions include/boost/math/distributions/detail/inv_discrete_quantile.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
Expand Down
13 changes: 13 additions & 0 deletions test/test_binomial.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down

0 comments on commit 071d338

Please sign in to comment.