Skip to content

Commit 74bb4bc

Browse files
authored
Fix inverse_discrete_quantile for large guess (#1007)
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); ```
1 parent 3d8e1d5 commit 74bb4bc

2 files changed

Lines changed: 18 additions & 10 deletions

File tree

include/boost/math/distributions/detail/inv_discrete_quantile.hpp

Lines changed: 6 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -302,15 +302,13 @@ inline typename Dist::value_type round_to_floor(const Dist& d, typename Dist::va
302302
//
303303
while(result != 0)
304304
{
305-
cc = result - 1;
305+
cc = floor(float_prior(result));
306306
if(cc < support(d).first)
307307
break;
308308
pp = c ? cdf(complement(d, cc)) : cdf(d, cc);
309-
if(pp == p)
310-
result = cc;
311-
else if(c ? pp > p : pp < p)
309+
if(c ? pp > p : pp < p)
312310
break;
313-
result -= 1;
311+
result = cc;
314312
}
315313

316314
return result;
@@ -336,15 +334,13 @@ inline typename Dist::value_type round_to_ceil(const Dist& d, typename Dist::val
336334
//
337335
while(true)
338336
{
339-
cc = result + 1;
337+
cc = ceil(float_next(result));
340338
if(cc > support(d).second)
341339
break;
342340
pp = c ? cdf(complement(d, cc)) : cdf(d, cc);
343-
if(pp == p)
344-
result = cc;
345-
else if(c ? pp < p : pp > p)
341+
if(c ? pp < p : pp > p)
346342
break;
347-
result += 1;
343+
result = cc;
348344
}
349345

350346
return result;

test/test_binomial.cpp

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -716,6 +716,18 @@ void test_spots(RealType T)
716716

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

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

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

0 commit comments

Comments
 (0)