Skip to content

Commit

Permalink
Fix inverse_discrete_quantile for large guess (#1007)
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 authored Aug 24, 2023
1 parent 3d8e1d5 commit 74bb4bc
Show file tree
Hide file tree
Showing 2 changed files with 18 additions and 10 deletions.
16 changes: 6 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,13 @@ inline typename Dist::value_type round_to_floor(const Dist& d, typename Dist::va
//
while(result != 0)
{
cc = result - 1;
cc = floor(float_prior(result));
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 +334,13 @@ inline typename Dist::value_type round_to_ceil(const Dist& d, typename Dist::val
//
while(true)
{
cc = result + 1;
cc = ceil(float_next(result));
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
12 changes: 12 additions & 0 deletions test/test_binomial.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -716,6 +716,18 @@ void test_spots(RealType T)

check_out_of_range<boost::math::binomial_distribution<RealType> >(1, 1); // (All) valid constructor parameter values.

// TODO: Generic ibeta_power_terms has accuracy issue when long
// double is not precise enough, causing overflow in this case.
if(!std::is_same<RealType, real_concept>::value ||
sizeof(boost::math::concepts::real_concept_base_type) > 8)
{
using namespace boost::math::policies;
typedef policy<discrete_quantile<integer_round_outwards> > Policy;
binomial_distribution<RealType, Policy> dist(9079765771874083840, 0.561815);
// Accuracy is not too important here; the main purpose is to
// make sure it is not stuck.
BOOST_CHECK_CLOSE(quantile(dist, 0.0365346), 5101148604445670400, 1e12);
}

} // template <class RealType>void test_spots(RealType)

Expand Down

0 comments on commit 74bb4bc

Please sign in to comment.