Skip to content

Commit f695d39

Browse files
jzmaddockjzmaddock
authored
Merge pull request #1406 from boostorg/issue1404
pFq: prevent spurious overflow.
2 parents 6dea961 + 45fe999 commit f695d39

3 files changed

Lines changed: 68 additions & 16 deletions

File tree

‎include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp‎

Lines changed: 29 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,7 @@
2323
unsigned set_crossover_locations(const Seq& aj, const Seq& bj, const Real& z, unsigned int* crossover_locations)
2424
{
2525
BOOST_MATH_STD_USING
26+
using nothrow_policy = typename boost::math::policies::normalise<boost::math::policies::policy<>, boost::math::policies::rounding_error<boost::math::policies::ignore_error>>::type;
2627
unsigned N_terms = 0;
2728

2829
if(aj.size() == 1 && bj.size() == 1)
@@ -55,13 +56,13 @@
5556
Real t = (-sqrt(sq) - b + z) / 2;
5657
if (t >= 0)
5758
{
58-
crossover_locations[N_terms] = itrunc(t);
59+
crossover_locations[N_terms] = itrunc(t, nothrow_policy());
5960
++N_terms;
6061
}
6162
t = (sqrt(sq) - b + z) / 2;
6263
if (t >= 0)
6364
{
64-
crossover_locations[N_terms] = itrunc(t);
65+
crossover_locations[N_terms] = itrunc(t, nothrow_policy());
6566
++N_terms;
6667
}
6768
}
@@ -71,13 +72,13 @@
7172
Real t = (-sqrt(sq) - b - z) / 2;
7273
if (t >= 0)
7374
{
74-
crossover_locations[N_terms] = itrunc(t);
75+
crossover_locations[N_terms] = itrunc(t, nothrow_policy());
7576
++N_terms;
7677
}
7778
t = (sqrt(sq) - b - z) / 2;
7879
if (t >= 0)
7980
{
80-
crossover_locations[N_terms] = itrunc(t);
81+
crossover_locations[N_terms] = itrunc(t, nothrow_policy());
8182
++N_terms;
8283
}
8384
}
@@ -110,7 +111,7 @@
110111
unsigned n = 0;
111112
for (auto bi = bj.begin(); bi != bj.end(); ++bi, ++n)
112113
{
113-
crossover_locations[n] = *bi >= 0 ? 0 : itrunc(-*bi) + 1;
114+
crossover_locations[n] = *bi >= 0 ? 0 : itrunc(-*bi, nothrow_policy()) + 1;
114115
}
115116
std::sort(crossover_locations, crossover_locations + bj.size(), std::less<Real>());
116117
N_terms = (unsigned)bj.size();
@@ -122,15 +123,27 @@
122123
std::pair<Real, Real> hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol, const Terminal& termination, long long& log_scale)
123124
{
124125
BOOST_MATH_STD_USING
126+
using nothrow_policy = typename boost::math::policies::normalise<Policy, boost::math::policies::rounding_error<boost::math::policies::ignore_error>>::type;
125127
Real result = 1;
126128
Real abs_result = 1;
127129
Real term = 1;
128130
Real term0 = 0;
129131
Real tol = boost::math::policies::get_epsilon<Real, Policy>();
130132
std::uintmax_t k = 0;
131133
Real upper_limit(sqrt(boost::math::tools::max_value<Real>())), diff;
134+
if ((tools::max_value<Real>() / fabs(z) < upper_limit))
135+
{
136+
upper_limit = tools::max_value<Real>() / fabs(z);
137+
}
138+
for (auto pa = aj.begin(); pa != aj.end(); ++pa)
139+
{
140+
if (tools::max_value<Real>() / fabs(*pa) < upper_limit)
141+
{
142+
upper_limit = tools::max_value<Real>() / fabs(*pa);
143+
}
144+
}
132145
Real lower_limit(1 / upper_limit);
133-
long long log_scaling_factor = lltrunc(boost::math::tools::log_max_value<Real>()) - 2;
146+
long long log_scaling_factor = lltrunc(boost::math::tools::log_max_value<Real>(), nothrow_policy()) - 2;
134147
Real scaling_factor = exp(Real(log_scaling_factor));
135148
Real term_m1;
136149
long long local_scaling = 0;
@@ -234,15 +247,18 @@
234247
log_scale += log_scaling_factor;
235248
local_scaling += log_scaling_factor;
236249
}
237-
if (fabs(abs_result) < lower_limit)
250+
/*
251+
* rescaling to avoid underflow is pointless here, given that the first term is 1.
252+
*
253+
else if ((fabs(abs_result) < lower_limit) && (fabs(abs_result) * scaling_factor < upper_limit))
238254
{
239255
abs_result *= scaling_factor;
240256
result *= scaling_factor;
241257
term *= scaling_factor;
242258
log_scale -= log_scaling_factor;
243259
local_scaling -= log_scaling_factor;
244260
}
245-
261+
*/
246262
if ((abs(result * tol) > abs(term)) && (abs(term0) > abs(term)))
247263
break;
248264
if (abs_result * tol > abs(result))
@@ -321,7 +337,7 @@
321337
else
322338
{
323339
int ls = 1;
324-
Real p = log_pochhammer(*ai, s, pol, &ls);
340+
Real p = log_pochhammer(static_cast<Real>(*ai), s, pol, &ls);
325341
s1 *= ls;
326342
term += p;
327343
loop_error_scale = (std::max)(p, loop_error_scale);
@@ -334,7 +350,7 @@
334350
for (auto bi = bj.begin(); bi != bj.end(); ++bi)
335351
{
336352
int ls = 1;
337-
Real p = log_pochhammer(*bi, s, pol, &ls);
353+
Real p = log_pochhammer(static_cast<Real>(*bi), s, pol, &ls);
338354
s2 *= ls;
339355
term -= p;
340356
loop_error_scale = (std::max)(p, loop_error_scale);
@@ -371,13 +387,13 @@
371387
if (term <= tools::log_min_value<Real>())
372388
{
373389
// rescale if we can:
374-
long long scale = lltrunc(floor(term - tools::log_min_value<Real>()) - 2);
390+
long long scale = lltrunc(floor(term - tools::log_min_value<Real>()) - 2, nothrow_policy());
375391
term -= scale;
376392
loop_scale += scale;
377393
}
378394
if (term > 10)
379395
{
380-
int scale = itrunc(floor(term));
396+
long long scale = lltrunc(floor(term), nothrow_policy());
381397
term -= scale;
382398
loop_scale += scale;
383399
}
@@ -402,7 +418,7 @@
402418
term /= scaling_factor;
403419
loop_scale += log_scaling_factor;
404420
}
405-
if (fabs(loop_result) < lower_limit)
421+
else if ((fabs(loop_result) < lower_limit) && (fabs(loop_result) * scaling_factor < upper_limit))
406422
{
407423
loop_result *= scaling_factor;
408424
loop_abs_result *= scaling_factor;

‎include/boost/math/special_functions/hypergeometric_pFq.hpp‎

Lines changed: 13 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -58,12 +58,22 @@ namespace boost {
5858
BOOST_MATH_STD_USING
5959

6060
long long scale = 0;
61+
static const char* function = "boost::math::hypergeometric_pFq<%1%>(%1%,%1%,%1%)";
6162
std::pair<value_type, value_type> r = boost::math::detail::hypergeometric_pFq_checked_series_impl(aj, bj, value_type(z), pol, boost::math::detail::iteration_terminator(boost::math::policies::get_max_series_iterations<forwarding_policy>()), scale);
62-
r.first *= exp(Real(scale));
63-
r.second *= exp(Real(scale));
63+
//
64+
// Overflow check:
65+
//
66+
if (static_cast<Real>(scale) > tools::log_max_value<Real>())
67+
return (r.first < 0 ? -1 : 1) * policies::raise_overflow_error<Real, Policy>(function, nullptr, pol);
68+
Real mul = exp(Real(scale));
69+
if(fabs(r.first) > 1)
70+
if(tools::max_value<Real>() / fabs(r.first) < mul)
71+
return (r.first < 0 ? -1 : 1) * policies::raise_overflow_error<Real, Policy>(function, nullptr, pol);
72+
r.first *= mul;
73+
r.second *= mul;
6474
if (p_abs_error)
6575
*p_abs_error = static_cast<Real>(r.second) * boost::math::tools::epsilon<Real>();
66-
return policies::checked_narrowing_cast<result_type, Policy>(r.first, "boost::math::hypergeometric_pFq<%1%>(%1%,%1%,%1%)");
76+
return policies::checked_narrowing_cast<result_type, Policy>(r.first, function);
6777
}
6878

6979
template <class Seq, class Real>

‎test/test_pFq.hpp‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -298,6 +298,31 @@ void test_spots_2F2(T, const char*)
298298
}
299299
}
300300

301+
template <class T>
302+
void test_special_cases(T, const char*)
303+
{
304+
typedef boost::math::policies::policy<boost::math::policies::overflow_error<boost::math::policies::ignore_error>, boost::math::policies::promote_double<false> > nothrow_policy;
305+
typedef boost::math::policies::policy<boost::math::policies::overflow_error<boost::math::policies::throw_on_error>, boost::math::policies::promote_double<false> > throw_policy;
306+
T error_rate;
307+
BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ 2 }, { 3 }, static_cast<T>(1e10), &error_rate, throw_policy()), std::overflow_error);
308+
BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ static_cast<T>(2) }, { static_cast<T>(3) }, static_cast<T>(1e10), &error_rate, throw_policy()), std::overflow_error);
309+
BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ 2 }, { 3 }, boost::math::tools::max_value<T>(), &error_rate, throw_policy()), std::overflow_error);
310+
BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ static_cast<T>(2) }, { static_cast<T>(3) }, boost::math::tools::max_value<T>(), &error_rate, throw_policy()), std::overflow_error);
311+
BOOST_MATH_IF_CONSTEXPR(std::numeric_limits<T>::is_specialized)
312+
{
313+
// Can't figure out how to make this work with real_concept (we get an evaluation_error):
314+
BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ boost::math::tools::max_value<T>() }, { static_cast<T>(3) }, static_cast<T>(0.5f), &error_rate, throw_policy()), std::overflow_error);
315+
}
316+
BOOST_MATH_IF_CONSTEXPR(std::numeric_limits<T>::has_infinity)
317+
{
318+
BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ 2 }, { 3 }, static_cast<T>(1e10), &error_rate, nothrow_policy()), std::numeric_limits<T>::infinity());
319+
BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ static_cast<T>(2) }, { static_cast<T>(3) }, static_cast<T>(1e10), &error_rate, nothrow_policy()), std::numeric_limits<T>::infinity());
320+
BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ 2 }, { 3 }, boost::math::tools::max_value<T>(), &error_rate, nothrow_policy()), std::numeric_limits<T>::infinity());
321+
BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ static_cast<T>(2) }, { static_cast<T>(3) }, boost::math::tools::max_value<T>(), &error_rate, nothrow_policy()), std::numeric_limits<T>::infinity());
322+
BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ boost::math::tools::max_value<T>() }, { static_cast<T>(3) }, static_cast<T>(0.5f), &error_rate, nothrow_policy()), std::numeric_limits<T>::infinity());
323+
}
324+
}
325+
301326
template <class T>
302327
void test_spots(T z, const char* type_name)
303328
{
@@ -309,4 +334,5 @@ void test_spots(T z, const char* type_name)
309334
test_spots_1F2(z, type_name);
310335
test_spots_2F2(z, type_name);
311336
test_spots_2F1(z, type_name);
337+
test_special_cases(z, type_name);
312338
}

0 commit comments

Comments
 (0)