From 9b0d249f682cacb888c5ff5dfd503bb91f382299 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Fri, 5 Apr 2019 18:54:10 +0100 Subject: [PATCH] pFq: Add rescaling if term when b crosses the origin underflows. --- .../hypergeometric_pFq_checked_series.hpp | 22 ++++++++++++++++--- 1 file changed, 19 insertions(+), 3 deletions(-) diff --git a/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp b/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp index e4a7c3c04..2ad01c351 100644 --- a/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp +++ b/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp @@ -114,9 +114,11 @@ { if (k < crossover_locations[n]) { - for(auto ai = aj.begin(); ai != aj.end(); ++ai) + for (auto ai = aj.begin(); ai != aj.end(); ++ai) + { if ((*ai < 0) && (floor(*ai) == *ai) && (*ai > crossover_locations[n])) return std::make_pair(result, abs_result); // b's will never cross the origin! + } // // b hasn't crossed the origin yet and the series may spring back into life at that point // so we need to jump forward to that term and then evaluate forwards and backwards from there: @@ -147,8 +149,22 @@ if (z < 0) s1 *= (s & 1 ? -1 : 1); term -= local_scaling; - //std::cout << "term = " << term << std::endl; - if (term > -tools::log_max_value()) + if (term <= tools::log_min_value()) + { + // rescale if we can: + int scale = itrunc(floor(term - tools::log_min_value())); + if (scale < tools::log_max_value()) + { + Real ex = exp(Real(-scale)); + if (tools::max_value() / ex > result) + { + term -= scale; + log_scale += scale; + result *= ex; + } + } + } + if (term > tools::log_min_value()) { if (term > 10) {