special_functions/detail/hypergeometric_pFq_checked_series.hpp

0001 ///////////////////////////////////////////////////////////////////////////////
0002 //  Copyright 2018 John Maddock
0003 //  Distributed under the Boost
0004 //  Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_HYPERGEOMETRIC_PFQ_SERIES_HPP_
0008 #define BOOST_HYPERGEOMETRIC_PFQ_SERIES_HPP_
0009
0010 #ifndef BOOST_MATH_PFQ_MAX_B_TERMS
0011 #  define BOOST_MATH_PFQ_MAX_B_TERMS 5
0012 #endif
0013
0014 #include <array>
0015 #include <cstdint>
0016 #include <boost/math/special_functions/gamma.hpp>
0017 #include <boost/math/special_functions/expm1.hpp>
0018 #include <boost/math/special_functions/detail/hypergeometric_series.hpp>
0019
0020   namespace boost { namespace math { namespace detail {
0021
0022      template <class Seq, class Real>
0023      unsigned set_crossover_locations(const Seq& aj, const Seq& bj, const Real& z, unsigned int* crossover_locations)
0024      {
0025         BOOST_MATH_STD_USING
0026         unsigned N_terms = 0;
0027
0028         if(aj.size() == 1 && bj.size() == 1)
0029         {
0030            //
0031            // For 1F1 we can work out where the peaks in the series occur,
0032            //  which is to say when:
0033            //
0034            // (a + k)z / (k(b + k)) == +-1
0035            //
0036            // Then we are at either a maxima or a minima in the series, and the
0037            // last such point must be a maxima since the series is globally convergent.
0038            // Potentially then we are solving 2 quadratic equations and have up to 4
0039            // solutions, any solutions which are complex or negative are discarded,
0040            // leaving us with 4 options:
0041            //
0042            // 0 solutions: The series is directly convergent.
0043            // 1 solution : The series diverges to a maxima before converging.
0044            // 2 solutions: The series is initially convergent, followed by divergence to a maxima before final convergence.
0045            // 3 solutions: The series diverges to a maxima, converges to a minima before diverging again to a second maxima before final convergence.
0046            // 4 solutions: The series converges to a minima before diverging to a maxima, converging to a minima, diverging to a second maxima and then converging.
0047            //
0048            // The first 2 situations are adequately handled by direct series evaluation, while the 2,3 and 4 solutions are not.
0049            //
0050            Real a = *aj.begin();
0051            Real b = *bj.begin();
0052            Real sq = 4 * a * z + b * b - 2 * b * z + z * z;
0053            if (sq >= 0)
0054            {
0055               Real t = (-sqrt(sq) - b + z) / 2;
0056               if (t >= 0)
0057               {
0058                  crossover_locations[N_terms] = itrunc(t);
0059                  ++N_terms;
0060               }
0061               t = (sqrt(sq) - b + z) / 2;
0062               if (t >= 0)
0063               {
0064                  crossover_locations[N_terms] = itrunc(t);
0065                  ++N_terms;
0066               }
0067            }
0068            sq = -4 * a * z + b * b + 2 * b * z + z * z;
0069            if (sq >= 0)
0070            {
0071               Real t = (-sqrt(sq) - b - z) / 2;
0072               if (t >= 0)
0073               {
0074                  crossover_locations[N_terms] = itrunc(t);
0075                  ++N_terms;
0076               }
0077               t = (sqrt(sq) - b - z) / 2;
0078               if (t >= 0)
0079               {
0080                  crossover_locations[N_terms] = itrunc(t);
0081                  ++N_terms;
0082               }
0083            }
0084            std::sort(crossover_locations, crossover_locations + N_terms, std::less<Real>());
0085            //
0086            // Now we need to discard every other terms, as these are the minima:
0087            //
0088            switch (N_terms)
0089            {
0090            case 0:
0091            case 1:
0092               break;
0093            case 2:
0094               crossover_locations[0] = crossover_locations[1];
0095               --N_terms;
0096               break;
0097            case 3:
0098               crossover_locations[1] = crossover_locations[2];
0099               --N_terms;
0100               break;
0101            case 4:
0102               crossover_locations[0] = crossover_locations[1];
0103               crossover_locations[1] = crossover_locations[3];
0104               N_terms -= 2;
0105               break;
0106            }
0107         }
0108         else
0109         {
0110            unsigned n = 0;
0111            for (auto bi = bj.begin(); bi != bj.end(); ++bi, ++n)
0112            {
0113               crossover_locations[n] = *bi >= 0 ? 0 : itrunc(-*bi) + 1;
0114            }
0115            std::sort(crossover_locations, crossover_locations + bj.size(), std::less<Real>());
0116            N_terms = (unsigned)bj.size();
0117         }
0118         return N_terms;
0119      }
0120
0121      template <class Seq, class Real, class Policy, class Terminal>
0122      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)
0123      {
0124         BOOST_MATH_STD_USING
0125         Real result = 1;
0126         Real abs_result = 1;
0127         Real term = 1;
0128         Real term0 = 0;
0129         Real tol = boost::math::policies::get_epsilon<Real, Policy>();
0130         std::uintmax_t k = 0;
0131         Real upper_limit(sqrt(boost::math::tools::max_value<Real>())), diff;
0132         Real lower_limit(1 / upper_limit);
0133         long long log_scaling_factor = lltrunc(boost::math::tools::log_max_value<Real>()) - 2;
0134         Real scaling_factor = exp(Real(log_scaling_factor));
0135         Real term_m1;
0136         long long local_scaling = 0;
0137         bool have_no_correct_bits = false;
0138
0139         if ((aj.size() == 1) && (bj.size() == 0))
0140         {
0141            if (fabs(z) > 1)
0142            {
0143               if ((z > 0) && (floor(*aj.begin()) != *aj.begin()))
0144               {
0145                  Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p == 1 and q == 0 and |z| > 1, result is imaginary", z, pol);
0146                  return std::make_pair(r, r);
0147               }
0148               std::pair<Real, Real> r = hypergeometric_pFq_checked_series_impl(aj, bj, Real(1 / z), pol, termination, log_scale);
0149
0150               #if (defined(__GNUC__) && __GNUC__ == 13)
0151               Real mul = pow(-z, Real(-*aj.begin()));
0152               #else
0153               Real mul = pow(-z, -*aj.begin());
0154               #endif
0155
0156               r.first *= mul;
0157               r.second *= mul;
0158               return r;
0159            }
0160         }
0161
0162         if (aj.size() > bj.size())
0163         {
0164            if (aj.size() == bj.size() + 1)
0165            {
0166               if (fabs(z) > 1)
0167               {
0168                  Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p == q+1 and |z| > 1, series does not converge", z, pol);
0169                  return std::make_pair(r, r);
0170               }
0171               if (fabs(z) == 1)
0172               {
0173                  Real s = 0;
0174                  for (auto i = bj.begin(); i != bj.end(); ++i)
0175                     s += *i;
0176                  for (auto i = aj.begin(); i != aj.end(); ++i)
0177                     s -= *i;
0178                  if ((z == 1) && (s <= 0))
0179                  {
0180                     Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p == q+1 and |z| == 1, in a situation where the series does not converge", z, pol);
0181                     return std::make_pair(r, r);
0182                  }
0183                  if ((z == -1) && (s <= -1))
0184                  {
0185                     Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p == q+1 and |z| == 1, in a situation where the series does not converge", z, pol);
0186                     return std::make_pair(r, r);
0187                  }
0188               }
0189            }
0190            else
0191            {
0192               Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p > q+1, series does not converge", z, pol);
0193               return std::make_pair(r, r);
0194            }
0195         }
0196
0197         while (!termination(k))
0198         {
0199            for (auto ai = aj.begin(); ai != aj.end(); ++ai)
0200            {
0201               term *= *ai + k;
0202            }
0203            if (term == 0)
0204            {
0205               // There is a negative integer in the aj's:
0206               return std::make_pair(result, abs_result);
0207            }
0208            for (auto bi = bj.begin(); bi != bj.end(); ++bi)
0209            {
0210               if (*bi + k == 0)
0211               {
0212                  // The series is undefined:
0213                  result = boost::math::policies::raise_domain_error("boost::math::hypergeometric_pFq<%1%>", "One of the b values was the negative integer %1%", *bi, pol);
0214                  return std::make_pair(result, result);
0215               }
0216               term /= *bi + k;
0217            }
0218            term *= z;
0219            ++k;
0220            term /= k;
0221            //std::cout << k << " " << *bj.begin() + k << " " << result << " " << term << /*" " << term_at_k(*aj.begin(), *bj.begin(), z, k, pol) <<*/ std::endl;
0222            result += term;
0223            abs_result += abs(term);
0224            //std::cout << "k = " << k << " term = " << term * exp(log_scale) << " result = " << result * exp(log_scale) << " abs_result = " << abs_result * exp(log_scale) << std::endl;
0225
0226            //
0227            // Rescaling:
0228            //
0229            if (fabs(abs_result) >= upper_limit)
0230            {
0231               abs_result /= scaling_factor;
0232               result /= scaling_factor;
0233               term /= scaling_factor;
0234               log_scale += log_scaling_factor;
0235               local_scaling += log_scaling_factor;
0236            }
0237            if (fabs(abs_result) < lower_limit)
0238            {
0239               abs_result *= scaling_factor;
0240               result *= scaling_factor;
0241               term *= scaling_factor;
0242               log_scale -= log_scaling_factor;
0243               local_scaling -= log_scaling_factor;
0244            }
0245
0246            if ((abs(result * tol) > abs(term)) && (abs(term0) > abs(term)))
0247               break;
0248            if (abs_result * tol > abs(result))
0249            {
0250               // Check if result is so small compared to abs_resuslt that there are no longer any
0251               // correct bits... we require two consecutive passes here before aborting to
0252               // avoid false positives when result transiently drops to near zero then rebounds.
0253               if (have_no_correct_bits)
0254               {
0255                  // We have no correct bits in the result... just give up!
0256                  result = boost::math::policies::raise_evaluation_error("boost::math::hypergeometric_pFq<%1%>", "Cancellation is so severe that no bits in the result are correct, last result was %1%", Real(result * exp(Real(log_scale))), pol);
0257                  return std::make_pair(result, result);
0258               }
0259               else
0260                  have_no_correct_bits = true;
0261            }
0262            else
0263               have_no_correct_bits = false;
0264            term0 = term;
0265         }
0266         //std::cout << "result = " << result << std::endl;
0267         //std::cout << "local_scaling = " << local_scaling << std::endl;
0268         //std::cout << "Norm result = " << std::setprecision(35) << boost::multiprecision::mpfr_float_50(result) * exp(boost::multiprecision::mpfr_float_50(local_scaling)) << std::endl;
0269         //
0270         // We have to be careful when one of the b's crosses the origin:
0271         //
0272         if(bj.size() > BOOST_MATH_PFQ_MAX_B_TERMS)
0273            policies::raise_domain_error<Real>("boost::math::hypergeometric_pFq<%1%>(Seq, Seq, %1%)",
0274               "The number of b terms must be less than the value of BOOST_MATH_PFQ_MAX_B_TERMS (" BOOST_MATH_STRINGIZE(BOOST_MATH_PFQ_MAX_B_TERMS)  "), but got %1%.",
0275               Real(bj.size()), pol);
0276
0277         unsigned crossover_locations[BOOST_MATH_PFQ_MAX_B_TERMS];
0278
0279         unsigned N_crossovers = set_crossover_locations(aj, bj, z, crossover_locations);
0280
0281         bool terminate = false;   // Set to true if one of the a's passes through the origin and terminates the series.
0282
0283         for (unsigned n = 0; n < N_crossovers; ++n)
0284         {
0285            if (k < crossover_locations[n])
0286            {
0287               for (auto ai = aj.begin(); ai != aj.end(); ++ai)
0288               {
0289                  if ((*ai < 0) && (floor(*ai) == *ai) && (*ai > static_cast<decltype(*ai)>(crossover_locations[n])))
0290                     return std::make_pair(result, abs_result);  // b's will never cross the origin!
0291               }
0292               //
0293               // local results:
0294               //
0295               Real loop_result = 0;
0296               Real loop_abs_result = 0;
0297               long long loop_scale = 0;
0298               //
0299               // loop_error_scale will be used to increase the size of the error
0300               // estimate (absolute sum), based on the errors inherent in calculating
0301               // the pochhammer symbols.
0302               //
0303               Real loop_error_scale = 0;
0304               //boost::multiprecision::mpfi_float err_est = 0;
0305               //
0306               // b hasn't crossed the origin yet and the series may spring back into life at that point
0307               // so we need to jump forward to that term and then evaluate forwards and backwards from there:
0308               //
0309               unsigned s = crossover_locations[n];
0310               std::uintmax_t backstop = k;
0311               long long s1(1), s2(1);
0312               term = 0;
0313               for (auto ai = aj.begin(); ai != aj.end(); ++ai)
0314               {
0315                  if ((floor(*ai) == *ai) && (*ai < 0) && (-*ai <= static_cast<decltype(*ai)>(s)))
0316                  {
0317                     // One of the a terms has passed through zero and terminated the series:
0318                     terminate = true;
0319                     break;
0320                  }
0321                  else
0322                  {
0323                     int ls = 1;
0324                     Real p = log_pochhammer(*ai, s, pol, &ls);
0325                     s1 *= ls;
0326                     term += p;
0327                     loop_error_scale = (std::max)(p, loop_error_scale);
0328                     //err_est += boost::multiprecision::mpfi_float(p);
0329                  }
0330               }
0331               //std::cout << "term = " << term << std::endl;
0332               if (terminate)
0333                  break;
0334               for (auto bi = bj.begin(); bi != bj.end(); ++bi)
0335               {
0336                  int ls = 1;
0337                  Real p = log_pochhammer(*bi, s, pol, &ls);
0338                  s2 *= ls;
0339                  term -= p;
0340                  loop_error_scale = (std::max)(p, loop_error_scale);
0341                  //err_est -= boost::multiprecision::mpfi_float(p);
0342               }
0343               //std::cout << "term = " << term << std::endl;
0344               Real p = lgamma(Real(s + 1), pol);
0345               term -= p;
0346               loop_error_scale = (std::max)(p, loop_error_scale);
0347               //err_est -= boost::multiprecision::mpfi_float(p);
0348               p = s * log(fabs(z));
0349               term += p;
0350               loop_error_scale = (std::max)(p, loop_error_scale);
0351               //err_est += boost::multiprecision::mpfi_float(p);
0352               //err_est = exp(err_est);
0353               //std::cout << err_est << std::endl;
0354               //
0355               // Convert loop_error scale to the absolute error
0356               // in term after exp is applied:
0357               //
0358               loop_error_scale *= tools::epsilon<Real>();
0359               //
0360               // Convert to relative error after exp:
0361               //
0362               loop_error_scale = fabs(expm1(loop_error_scale, pol));
0363               //
0364               // Convert to multiplier for the error term:
0365               //
0366               loop_error_scale /= tools::epsilon<Real>();
0367
0368               if (z < 0)
0369                  s1 *= (s & 1 ? -1 : 1);
0370
0371               if (term <= tools::log_min_value<Real>())
0372               {
0373                  // rescale if we can:
0374                  long long scale = lltrunc(floor(term - tools::log_min_value<Real>()) - 2);
0375                  term -= scale;
0376                  loop_scale += scale;
0377               }
0378                if (term > 10)
0379                {
0380                   int scale = itrunc(floor(term));
0381                   term -= scale;
0382                   loop_scale += scale;
0383                }
0384                //std::cout << "term = " << term << std::endl;
0385                term = s1 * s2 * exp(term);
0386                //std::cout << "term = " << term << std::endl;
0387                //std::cout << "loop_scale = " << loop_scale << std::endl;
0388                k = s;
0389                term0 = term;
0390                long long saved_loop_scale = loop_scale;
0391                bool terms_are_growing = true;
0392                bool trivial_small_series_check = false;
0393                do
0394                {
0395                   loop_result += term;
0396                   loop_abs_result += fabs(term);
0397                   //std::cout << "k = " << k << " term = " << term * exp(loop_scale) << " result = " << loop_result * exp(loop_scale) << " abs_result = " << loop_abs_result * exp(loop_scale) << std::endl;
0398                   if (fabs(loop_result) >= upper_limit)
0399                   {
0400                      loop_result /= scaling_factor;
0401                      loop_abs_result /= scaling_factor;
0402                      term /= scaling_factor;
0403                      loop_scale += log_scaling_factor;
0404                   }
0405                   if (fabs(loop_result) < lower_limit)
0406                   {
0407                      loop_result *= scaling_factor;
0408                      loop_abs_result *= scaling_factor;
0409                      term *= scaling_factor;
0410                      loop_scale -= log_scaling_factor;
0411                   }
0412                   term_m1 = term;
0413                   for (auto ai = aj.begin(); ai != aj.end(); ++ai)
0414                   {
0415                      term *= *ai + k;
0416                   }
0417                   if (term == 0)
0418                   {
0419                      // There is a negative integer in the aj's:
0420                      return std::make_pair(result, abs_result);
0421                   }
0422                   for (auto bi = bj.begin(); bi != bj.end(); ++bi)
0423                   {
0424                      if (*bi + k == 0)
0425                      {
0426                         // The series is undefined:
0427                         result = boost::math::policies::raise_domain_error("boost::math::hypergeometric_pFq<%1%>", "One of the b values was the negative integer %1%", *bi, pol);
0428                         return std::make_pair(result, result);
0429                      }
0430                      term /= *bi + k;
0431                   }
0432                   term *= z / (k + 1);
0433
0434                   ++k;
0435                   diff = fabs(term / loop_result);
0436                   terms_are_growing = fabs(term) > fabs(term_m1);
0437                   if (!trivial_small_series_check && !terms_are_growing)
0438                   {
0439                      //
0440                      // Now that we have started to converge, check to see if the value of
0441                      // this local sum is trivially small compared to the result.  If so
0442                      // abort this part of the series.
0443                      //
0444                      trivial_small_series_check = true;
0445                      Real d;
0446                      if (loop_scale > local_scaling)
0447                      {
0448                         long long rescale = local_scaling - loop_scale;
0449                         if (rescale < tools::log_min_value<Real>())
0450                            d = 1;  // arbitrary value, we want to keep going
0451                         else
0452                            d = fabs(term / (result * exp(Real(rescale))));
0453                      }
0454                      else
0455                      {
0456                         long long rescale = loop_scale - local_scaling;
0457                         if (rescale < tools::log_min_value<Real>())
0458                            d = 0;  // terminate this loop
0459                         else
0460                            d = fabs(term * exp(Real(rescale)) / result);
0461                      }
0462                      if (d < boost::math::policies::get_epsilon<Real, Policy>())
0463                         break;
0464                   }
0465                } while (!termination(k - s) && ((diff > boost::math::policies::get_epsilon<Real, Policy>()) || terms_are_growing));
0466
0467                //std::cout << "Norm loop result = " << std::setprecision(35) << boost::multiprecision::mpfr_float_50(loop_result)* exp(boost::multiprecision::mpfr_float_50(loop_scale)) << std::endl;
0468                //
0469                // We now need to combine the results of the first series summation with whatever
0470                // local results we have now.  First though, rescale abs_result by loop_error_scale
0471                // to factor in the error in the pochhammer terms at the start of this block:
0472                //
0473                std::uintmax_t next_backstop = k;
0474                loop_abs_result += loop_error_scale * fabs(loop_result);
0475                if (loop_scale > local_scaling)
0476                {
0477                   //
0478                   // Need to shrink previous result:
0479                   //
0480                   long long rescale = local_scaling - loop_scale;
0481                   local_scaling = loop_scale;
0482                   log_scale -= rescale;
0483                   Real ex = exp(Real(rescale));
0484                   result *= ex;
0485                   abs_result *= ex;
0486                   result += loop_result;
0487                   abs_result += loop_abs_result;
0488                }
0489                else if (local_scaling > loop_scale)
0490                {
0491                   //
0492                   // Need to shrink local result:
0493                   //
0494                   long long rescale = loop_scale - local_scaling;
0495                   Real ex = exp(Real(rescale));
0496                   loop_result *= ex;
0497                   loop_abs_result *= ex;
0498                   result += loop_result;
0499                   abs_result += loop_abs_result;
0500                }
0501                else
0502                {
0503                   result += loop_result;
0504                   abs_result += loop_abs_result;
0505                }
0506                //
0507                // Now go backwards as well:
0508                //
0509                k = s;
0510                term = term0;
0511                loop_result = 0;
0512                loop_abs_result = 0;
0513                loop_scale = saved_loop_scale;
0514                trivial_small_series_check = false;
0515                do
0516                {
0517                   --k;
0518                   if (k == backstop)
0519                      break;
0520                   term_m1 = term;
0521                   for (auto ai = aj.begin(); ai != aj.end(); ++ai)
0522                   {
0523                      term /= *ai + k;
0524                   }
0525                   for (auto bi = bj.begin(); bi != bj.end(); ++bi)
0526                   {
0527                      if (*bi + k == 0)
0528                      {
0529                         // The series is undefined:
0530                         result = boost::math::policies::raise_domain_error("boost::math::hypergeometric_pFq<%1%>", "One of the b values was the negative integer %1%", *bi, pol);
0531                         return std::make_pair(result, result);
0532                      }
0533                      term *= *bi + k;
0534                   }
0535                   term *= (k + 1) / z;
0536                   loop_result += term;
0537                   loop_abs_result += fabs(term);
0538
0539                   if (!trivial_small_series_check && (fabs(term) < fabs(term_m1)))
0540                   {
0541                      //
0542                      // Now that we have started to converge, check to see if the value of
0543                      // this local sum is trivially small compared to the result.  If so
0544                      // abort this part of the series.
0545                      //
0546                      trivial_small_series_check = true;
0547                      Real d;
0548                      if (loop_scale > local_scaling)
0549                      {
0550                         long long rescale = local_scaling - loop_scale;
0551                         if (rescale < tools::log_min_value<Real>())
0552                            d = 1;  // keep going
0553                         else
0554                            d = fabs(term / (result * exp(Real(rescale))));
0555                      }
0556                      else
0557                      {
0558                         long long rescale = loop_scale - local_scaling;
0559                         if (rescale < tools::log_min_value<Real>())
0560                            d = 0;  // stop, underflow
0561                         else
0562                            d = fabs(term * exp(Real(rescale)) / result);
0563                      }
0564                      if (d < boost::math::policies::get_epsilon<Real, Policy>())
0565                         break;
0566                   }
0567
0568                   //std::cout << "k = " << k << " result = " << result << " abs_result = " << abs_result << std::endl;
0569                   if (fabs(loop_result) >= upper_limit)
0570                   {
0571                      loop_result /= scaling_factor;
0572                      loop_abs_result /= scaling_factor;
0573                      term /= scaling_factor;
0574                      loop_scale += log_scaling_factor;
0575                   }
0576                   if (fabs(loop_result) < lower_limit)
0577                   {
0578                      loop_result *= scaling_factor;
0579                      loop_abs_result *= scaling_factor;
0580                      term *= scaling_factor;
0581                      loop_scale -= log_scaling_factor;
0582                   }
0583                   diff = fabs(term / loop_result);
0584                } while (!termination(s - k) && ((diff > boost::math::policies::get_epsilon<Real, Policy>()) || (fabs(term) > fabs(term_m1))));
0585
0586                //std::cout << "Norm loop result = " << std::setprecision(35) << boost::multiprecision::mpfr_float_50(loop_result)* exp(boost::multiprecision::mpfr_float_50(loop_scale)) << std::endl;
0587                //
0588                // We now need to combine the results of the first series summation with whatever
0589                // local results we have now.  First though, rescale abs_result by loop_error_scale
0590                // to factor in the error in the pochhammer terms at the start of this block:
0591                //
0592                loop_abs_result += loop_error_scale * fabs(loop_result);
0593                //
0594                if (loop_scale > local_scaling)
0595                {
0596                   //
0597                   // Need to shrink previous result:
0598                   //
0599                   long long rescale = local_scaling - loop_scale;
0600                   local_scaling = loop_scale;
0601                   log_scale -= rescale;
0602                   Real ex = exp(Real(rescale));
0603                   result *= ex;
0604                   abs_result *= ex;
0605                   result += loop_result;
0606                   abs_result += loop_abs_result;
0607                }
0608                else if (local_scaling > loop_scale)
0609                {
0610                   //
0611                   // Need to shrink local result:
0612                   //
0613                   long long rescale = loop_scale - local_scaling;
0614                   Real ex = exp(Real(rescale));
0615                   loop_result *= ex;
0616                   loop_abs_result *= ex;
0617                   result += loop_result;
0618                   abs_result += loop_abs_result;
0619                }
0620                else
0621                {
0622                   result += loop_result;
0623                   abs_result += loop_abs_result;
0624                }
0625                //
0626                // Reset k to the largest k we reached
0627                //
0628                k = next_backstop;
0629            }
0630         }
0631
0632         return std::make_pair(result, abs_result);
0633      }
0634
0635      struct iteration_terminator
0636      {
0637         iteration_terminator(std::uintmax_t i) : m(i) {}
0638
0639         bool operator()(std::uintmax_t v) const { return v >= m; }
0640
0641         std::uintmax_t m;
0642      };
0643
0644      template <class Seq, class Real, class Policy>
0645      Real hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol, long long& log_scale)
0646      {
0647         BOOST_MATH_STD_USING
0648         iteration_terminator term(boost::math::policies::get_max_series_iterations<Policy>());
0649         std::pair<Real, Real> result = hypergeometric_pFq_checked_series_impl(aj, bj, z, pol, term, log_scale);
0650         //
0651         // Check to see how many digits we've lost, if it's more than half, raise an evaluation error -
0652         // this is an entirely arbitrary cut off, but not unreasonable.
0653         //
0654         if (result.second * sqrt(boost::math::policies::get_epsilon<Real, Policy>()) > abs(result.first))
0655         {
0656            return boost::math::policies::raise_evaluation_error("boost::math::hypergeometric_pFq<%1%>", "Cancellation is so severe that fewer than half the bits in the result are correct, last result was %1%", Real(result.first * exp(Real(log_scale))), pol);
0657         }
0658         return result.first;
0659      }
0660
0661      template <class Real, class Policy>
0662      inline Real hypergeometric_1F1_checked_series_impl(const Real& a, const Real& b, const Real& z, const Policy& pol, long long& log_scale)
0663      {
0664         std::array<Real, 1> aj = { a };
0665         std::array<Real, 1> bj = { b };
0666         return hypergeometric_pFq_checked_series_impl(aj, bj, z, pol, log_scale);
0667      }
0668
0669   } } } // namespaces
0670
0671 #endif // BOOST_HYPERGEOMETRIC_PFQ_SERIES_HPP_