boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp

   1 ///////////////////////////////////////////////////////////////////////////////
   2 //  Copyright 2018 John Maddock
   3 //  Distributed under the Boost
   4 //  Software License, Version 1.0. (See accompanying file
   5 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 #ifndef BOOST_HYPERGEOMETRIC_PFQ_SERIES_HPP_
   8 #define BOOST_HYPERGEOMETRIC_PFQ_SERIES_HPP_
   9
  10 #ifndef BOOST_MATH_PFQ_MAX_B_TERMS
  11 #  define BOOST_MATH_PFQ_MAX_B_TERMS 5
  12 #endif
  13
  14 #include <boost/array.hpp>
  15 #include <boost/math/special_functions/detail/hypergeometric_series.hpp>
  16
  17   namespace boost { namespace math { namespace detail {
  18
  19      template <class Seq, class Real>
  20      unsigned set_crossover_locations(const Seq& aj, const Seq& bj, const Real& z, unsigned int* crossover_locations)
  21      {
  22         BOOST_MATH_STD_USING
  23         unsigned N_terms = 0;
  24
  25         if(aj.size() == 1 && bj.size() == 1)
  26         {
  27            //
  28            // For 1F1 we can work out where the peaks in the series occur,
  29            //  which is to say when:
  30            //
  31            // (a + k)z / (k(b + k)) == +-1
  32            //
  33            // Then we are at either a maxima or a minima in the series, and the
  34            // last such point must be a maxima since the series is globally convergent.
  35            // Potentially then we are solving 2 quadratic equations and have up to 4
  36            // solutions, any solutions which are complex or negative are discarded,
  37            // leaving us with 4 options:
  38            //
  39            // 0 solutions: The series is directly convergent.
  40            // 1 solution : The series diverges to a maxima before coverging.
  41            // 2 solutions: The series is initially convergent, followed by divergence to a maxima before final convergence.
  42            // 3 solutions: The series diverges to a maxima, converges to a minima before diverging again to a second maxima before final convergence.
  43            // 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.
  44            //
  45            // The first 2 situations are adequately handled by direct series evalution, while the 2,3 and 4 solutions are not.
  46            //
  47            Real a = *aj.begin();
  48            Real b = *bj.begin();
  49            Real sq = 4 * a * z + b * b - 2 * b * z + z * z;
  50            if (sq >= 0)
  51            {
  52               Real t = (-sqrt(sq) - b + z) / 2;
  53               if (t >= 0)
  54               {
  55                  crossover_locations[N_terms] = itrunc(t);
  56                  ++N_terms;
  57               }
  58               t = (sqrt(sq) - b + z) / 2;
  59               if (t >= 0)
  60               {
  61                  crossover_locations[N_terms] = itrunc(t);
  62                  ++N_terms;
  63               }
  64            }
  65            sq = -4 * a * z + b * b + 2 * b * z + z * z;
  66            if (sq >= 0)
  67            {
  68               Real t = (-sqrt(sq) - b - z) / 2;
  69               if (t >= 0)
  70               {
  71                  crossover_locations[N_terms] = itrunc(t);
  72                  ++N_terms;
  73               }
  74               t = (sqrt(sq) - b - z) / 2;
  75               if (t >= 0)
  76               {
  77                  crossover_locations[N_terms] = itrunc(t);
  78                  ++N_terms;
  79               }
  80            }
  81            std::sort(crossover_locations, crossover_locations + N_terms, std::less<Real>());
  82            //
  83            // Now we need to discard every other terms, as these are the minima:
  84            //
  85            switch (N_terms)
  86            {
  87            case 0:
  88            case 1:
  89               break;
  90            case 2:
  91               crossover_locations[0] = crossover_locations[1];
  92               --N_terms;
  93               break;
  94            case 3:
  95               crossover_locations[1] = crossover_locations[2];
  96               --N_terms;
  97               break;
  98            case 4:
  99               crossover_locations[0] = crossover_locations[1];
 100               crossover_locations[1] = crossover_locations[3];
 101               N_terms -= 2;
 102               break;
 103            }
 104         }
 105         else
 106         {
 107            unsigned n = 0;
 108            for (auto bi = bj.begin(); bi != bj.end(); ++bi, ++n)
 109            {
 110               crossover_locations[n] = *bi >= 0 ? 0 : itrunc(-*bi) + 1;
 111            }
 112            std::sort(crossover_locations, crossover_locations + bj.size(), std::less<Real>());
 113            N_terms = (unsigned)bj.size();
 114         }
 115         return N_terms;
 116      }
 117
 118      template <class Seq, class Real, class Policy, class Terminal>
 119      std::pair<Real, Real> hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol, const Terminal& termination, int& log_scale)
 120      {
 121         BOOST_MATH_STD_USING
 122         Real result = 1;
 123         Real abs_result = 1;
 124         Real term = 1;
 125         Real term0 = 0;
 126         Real tol = boost::math::policies::get_epsilon<Real, Policy>();
 127         boost::uintmax_t k = 0;
 128         Real upper_limit(sqrt(boost::math::tools::max_value<Real>())), diff;
 129         Real lower_limit(1 / upper_limit);
 130         int log_scaling_factor = itrunc(boost::math::tools::log_max_value<Real>()) - 2;
 131         Real scaling_factor = exp(Real(log_scaling_factor));
 132         Real term_m1;
 133         int local_scaling = 0;
 134
 135         if ((aj.size() == 1) && (bj.size() == 0))
 136         {
 137            if (fabs(z) > 1)
 138            {
 139               if ((z > 0) && (floor(*aj.begin()) != *aj.begin()))
 140               {
 141                  Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p == 1 and q == 0 and |z| > 1, result is imaginary", z, pol);
 142                  return std::make_pair(r, r);
 143               }
 144               std::pair<Real, Real> r = hypergeometric_pFq_checked_series_impl(aj, bj, Real(1 / z), pol, termination, log_scale);
 145               Real mul = pow(-z, -*aj.begin());
 146               r.first *= mul;
 147               r.second *= mul;
 148               return r;
 149            }
 150         }
 151
 152         if (aj.size() > bj.size())
 153         {
 154            if (aj.size() == bj.size() + 1)
 155            {
 156               if (fabs(z) > 1)
 157               {
 158                  Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p == q+1 and |z| > 1, series does not converge", z, pol);
 159                  return std::make_pair(r, r);
 160               }
 161               if (fabs(z) == 1)
 162               {
 163                  Real s = 0;
 164                  for (auto i = bj.begin(); i != bj.end(); ++i)
 165                     s += *i;
 166                  for (auto i = aj.begin(); i != aj.end(); ++i)
 167                     s -= *i;
 168                  if ((z == 1) && (s <= 0))
 169                  {
 170                     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);
 171                     return std::make_pair(r, r);
 172                  }
 173                  if ((z == -1) && (s <= -1))
 174                  {
 175                     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);
 176                     return std::make_pair(r, r);
 177                  }
 178               }
 179            }
 180            else
 181            {
 182               Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p > q+1, series does not converge", z, pol);
 183               return std::make_pair(r, r);
 184            }
 185         }
 186
 187         while (!termination(k))
 188         {
 189            for (auto ai = aj.begin(); ai != aj.end(); ++ai)
 190            {
 191               term *= *ai + k;
 192            }
 193            if (term == 0)
 194            {
 195               // There is a negative integer in the aj's:
 196               return std::make_pair(result, abs_result);
 197            }
 198            for (auto bi = bj.begin(); bi != bj.end(); ++bi)
 199            {
 200               if (*bi + k == 0)
 201               {
 202                  // The series is undefined:
 203                  result = boost::math::policies::raise_domain_error("boost::math::hypergeometric_pFq<%1%>", "One of the b values was the negative integer %1%", *bi, pol);
 204                  return std::make_pair(result, result);
 205               }
 206               term /= *bi + k;
 207            }
 208            term *= z;
 209            ++k;
 210            term /= k;
 211            //std::cout << k << " " << *bj.begin() + k << " " << result << " " << term << /*" " << term_at_k(*aj.begin(), *bj.begin(), z, k, pol) <<*/ std::endl;
 212            result += term;
 213            abs_result += abs(term);
 214            //std::cout << "k = " << k << " term = " << term * exp(log_scale) << " result = " << result * exp(log_scale) << " abs_result = " << abs_result * exp(log_scale) << std::endl;
 215
 216            //
 217            // Rescaling:
 218            //
 219            if (fabs(abs_result) >= upper_limit)
 220            {
 221               abs_result /= scaling_factor;
 222               result /= scaling_factor;
 223               term /= scaling_factor;
 224               log_scale += log_scaling_factor;
 225               local_scaling += log_scaling_factor;
 226            }
 227            if (fabs(abs_result) < lower_limit)
 228            {
 229               abs_result *= scaling_factor;
 230               result *= scaling_factor;
 231               term *= scaling_factor;
 232               log_scale -= log_scaling_factor;
 233               local_scaling -= log_scaling_factor;
 234            }
 235
 236            if ((abs(result * tol) > abs(term)) && (abs(term0) > abs(term)))
 237               break;
 238            if (abs_result * tol > abs(result))
 239            {
 240               // We have no correct bits in the result... just give up!
 241               result = boost::math::policies::raise_evaluation_error("boost::math::hypergeometric_pFq<%1%>", "Cancellation is so severe that no bits in the reuslt are correct, last result was %1%", Real(result * exp(Real(log_scale))), pol);
 242               return std::make_pair(result, result);
 243            }
 244            term0 = term;
 245         }
 246         //std::cout << "result = " << result << std::endl;
 247         //std::cout << "local_scaling = " << local_scaling << std::endl;
 248         //std::cout << "Norm result = " << std::setprecision(35) << boost::multiprecision::mpfr_float_50(result) * exp(boost::multiprecision::mpfr_float_50(local_scaling)) << std::endl;
 249         //
 250         // We have to be careful when one of the b's crosses the origin:
 251         //
 252         if(bj.size() > BOOST_MATH_PFQ_MAX_B_TERMS)
 253            policies::raise_domain_error<Real>("boost::math::hypergeometric_pFq<%1%>(Seq, Seq, %1%)",
 254               "The number of b terms must be less than the value of BOOST_MATH_PFQ_MAX_B_TERMS (" BOOST_STRINGIZE(BOOST_MATH_PFQ_MAX_B_TERMS)  "), but got %1%.",
 255               Real(bj.size()), pol);
 256
 257         unsigned crossover_locations[BOOST_MATH_PFQ_MAX_B_TERMS];
 258
 259         unsigned N_crossovers = set_crossover_locations(aj, bj, z, crossover_locations);
 260
 261         bool terminate = false;   // Set to true if one of the a's passes through the origin and terminates the series.
 262
 263         for (unsigned n = 0; n < N_crossovers; ++n)
 264         {
 265            if (k < crossover_locations[n])
 266            {
 267               for (auto ai = aj.begin(); ai != aj.end(); ++ai)
 268               {
 269                  if ((*ai < 0) && (floor(*ai) == *ai) && (*ai > crossover_locations[n]))
 270                     return std::make_pair(result, abs_result);  // b's will never cross the origin!
 271               }
 272               //
 273               // local results:
 274               //
 275               Real loop_result = 0;
 276               Real loop_abs_result = 0;
 277               int loop_scale = 0;
 278               //
 279               // loop_error_scale will be used to increase the size of the error
 280               // estimate (absolute sum), based on the errors inherent in calculating
 281               // the pochhammer symbols.
 282               //
 283               Real loop_error_scale = 0;
 284               //boost::multiprecision::mpfi_float err_est = 0;
 285               //
 286               // b hasn't crossed the origin yet and the series may spring back into life at that point
 287               // so we need to jump forward to that term and then evaluate forwards and backwards from there:
 288               //
 289               unsigned s = crossover_locations[n];
 290               boost::uintmax_t backstop = k;
 291               int s1(1), s2(1);
 292               term = 0;
 293               for (auto ai = aj.begin(); ai != aj.end(); ++ai)
 294               {
 295                  if ((floor(*ai) == *ai) && (*ai < 0) && (-*ai <= s))
 296                  {
 297                     // One of the a terms has passed through zero and terminated the series:
 298                     terminate = true;
 299                     break;
 300                  }
 301                  else
 302                  {
 303                     int ls = 1;
 304                     Real p = log_pochhammer(*ai, s, pol, &ls);
 305                     s1 *= ls;
 306                     term += p;
 307                     loop_error_scale = (std::max)(p, loop_error_scale);
 308                     //err_est += boost::multiprecision::mpfi_float(p);
 309                  }
 310               }
 311               //std::cout << "term = " << term << std::endl;
 312               if (terminate)
 313                  break;
 314               for (auto bi = bj.begin(); bi != bj.end(); ++bi)
 315               {
 316                  int ls = 1;
 317                  Real p = log_pochhammer(*bi, s, pol, &ls);
 318                  s2 *= ls;
 319                  term -= p;
 320                  loop_error_scale = (std::max)(p, loop_error_scale);
 321                  //err_est -= boost::multiprecision::mpfi_float(p);
 322               }
 323               //std::cout << "term = " << term << std::endl;
 324               Real p = lgamma(Real(s + 1), pol);
 325               term -= p;
 326               loop_error_scale = (std::max)(p, loop_error_scale);
 327               //err_est -= boost::multiprecision::mpfi_float(p);
 328               p = s * log(fabs(z));
 329               term += p;
 330               loop_error_scale = (std::max)(p, loop_error_scale);
 331               //err_est += boost::multiprecision::mpfi_float(p);
 332               //err_est = exp(err_est);
 333               //std::cout << err_est << std::endl;
 334               //
 335               // Convert loop_error scale to the absolute error
 336               // in term after exp is applied:
 337               //
 338               loop_error_scale *= tools::epsilon<Real>();
 339               //
 340               // Convert to relative error after exp:
 341               //
 342               loop_error_scale = fabs(expm1(loop_error_scale));
 343               //
 344               // Convert to multiplier for the error term:
 345               //
 346               loop_error_scale /= tools::epsilon<Real>();
 347
 348               if (z < 0)
 349                  s1 *= (s & 1 ? -1 : 1);
 350
 351               if (term <= tools::log_min_value<Real>())
 352               {
 353                  // rescale if we can:
 354                  int scale = itrunc(floor(term - tools::log_min_value<Real>()) - 2);
 355                  term -= scale;
 356                  loop_scale += scale;
 357               }
 358                if (term > 10)
 359                {
 360                   int scale = itrunc(floor(term));
 361                   term -= scale;
 362                   loop_scale += scale;
 363                }
 364                //std::cout << "term = " << term << std::endl;
 365                term = s1 * s2 * exp(term);
 366                //std::cout << "term = " << term << std::endl;
 367                //std::cout << "loop_scale = " << loop_scale << std::endl;
 368                k = s;
 369                term0 = term;
 370                int saved_loop_scale = loop_scale;
 371                bool terms_are_growing = true;
 372                bool trivial_small_series_check = false;
 373                do
 374                {
 375                   loop_result += term;
 376                   loop_abs_result += fabs(term);
 377                   //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;
 378                   if (fabs(loop_result) >= upper_limit)
 379                   {
 380                      loop_result /= scaling_factor;
 381                      loop_abs_result /= scaling_factor;
 382                      term /= scaling_factor;
 383                      loop_scale += log_scaling_factor;
 384                   }
 385                   if (fabs(loop_result) < lower_limit)
 386                   {
 387                      loop_result *= scaling_factor;
 388                      loop_abs_result *= scaling_factor;
 389                      term *= scaling_factor;
 390                      loop_scale -= log_scaling_factor;
 391                   }
 392                   term_m1 = term;
 393                   for (auto ai = aj.begin(); ai != aj.end(); ++ai)
 394                   {
 395                      term *= *ai + k;
 396                   }
 397                   if (term == 0)
 398                   {
 399                      // There is a negative integer in the aj's:
 400                      return std::make_pair(result, abs_result);
 401                   }
 402                   for (auto bi = bj.begin(); bi != bj.end(); ++bi)
 403                   {
 404                      if (*bi + k == 0)
 405                      {
 406                         // The series is undefined:
 407                         result = boost::math::policies::raise_domain_error("boost::math::hypergeometric_pFq<%1%>", "One of the b values was the negative integer %1%", *bi, pol);
 408                         return std::make_pair(result, result);
 409                      }
 410                      term /= *bi + k;
 411                   }
 412                   term *= z / (k + 1);
 413
 414                   ++k;
 415                   diff = fabs(term / loop_result);
 416                   terms_are_growing = fabs(term) > fabs(term_m1);
 417                   if (!trivial_small_series_check && !terms_are_growing)
 418                   {
 419                      //
 420                      // Now that we have started to converge, check to see if the value of
 421                      // this local sum is trivially small compared to the result.  If so
 422                      // abort this part of the series.
 423                      //
 424                      trivial_small_series_check = true;
 425                      Real d;
 426                      if (loop_scale > local_scaling)
 427                      {
 428                         int rescale = local_scaling - loop_scale;
 429                         if (rescale < tools::log_min_value<Real>())
 430                            d = 1;  // arbtrary value, we want to keep going
 431                         else
 432                            d = fabs(term / (result * exp(Real(rescale))));
 433                      }
 434                      else
 435                      {
 436                         int rescale = loop_scale - local_scaling;
 437                         if (rescale < tools::log_min_value<Real>())
 438                            d = 0;  // terminate this loop
 439                         else
 440                            d = fabs(term * exp(Real(rescale)) / result);
 441                      }
 442                      if (d < boost::math::policies::get_epsilon<Real, Policy>())
 443                         break;
 444                   }
 445                } while (!termination(k - s) && ((diff > boost::math::policies::get_epsilon<Real, Policy>()) || terms_are_growing));
 446
 447                //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;
 448                //
 449                // We now need to combine the results of the first series summation with whatever
 450                // local results we have now.  First though, rescale abs_result by loop_error_scale
 451                // to factor in the error in the pochhammer terms at the start of this block:
 452                //
 453                boost::uintmax_t next_backstop = k;
 454                loop_abs_result += loop_error_scale * fabs(loop_result);
 455                if (loop_scale > local_scaling)
 456                {
 457                   //
 458                   // Need to shrink previous result:
 459                   //
 460                   int rescale = local_scaling - loop_scale;
 461                   local_scaling = loop_scale;
 462                   log_scale -= rescale;
 463                   Real ex = exp(Real(rescale));
 464                   result *= ex;
 465                   abs_result *= ex;
 466                   result += loop_result;
 467                   abs_result += loop_abs_result;
 468                }
 469                else if (local_scaling > loop_scale)
 470                {
 471                   //
 472                   // Need to shrink local result:
 473                   //
 474                   int rescale = loop_scale - local_scaling;
 475                   Real ex = exp(Real(rescale));
 476                   loop_result *= ex;
 477                   loop_abs_result *= ex;
 478                   result += loop_result;
 479                   abs_result += loop_abs_result;
 480                }
 481                else
 482                {
 483                   result += loop_result;
 484                   abs_result += loop_abs_result;
 485                }
 486                //
 487                // Now go backwards as well:
 488                //
 489                k = s;
 490                term = term0;
 491                loop_result = 0;
 492                loop_abs_result = 0;
 493                loop_scale = saved_loop_scale;
 494                trivial_small_series_check = false;
 495                do
 496                {
 497                   --k;
 498                   if (k == backstop)
 499                      break;
 500                   term_m1 = term;
 501                   for (auto ai = aj.begin(); ai != aj.end(); ++ai)
 502                   {
 503                      term /= *ai + k;
 504                   }
 505                   for (auto bi = bj.begin(); bi != bj.end(); ++bi)
 506                   {
 507                      if (*bi + k == 0)
 508                      {
 509                         // The series is undefined:
 510                         result = boost::math::policies::raise_domain_error("boost::math::hypergeometric_pFq<%1%>", "One of the b values was the negative integer %1%", *bi, pol);
 511                         return std::make_pair(result, result);
 512                      }
 513                      term *= *bi + k;
 514                   }
 515                   term *= (k + 1) / z;
 516                   loop_result += term;
 517                   loop_abs_result += fabs(term);
 518
 519                   if (!trivial_small_series_check && (fabs(term) < fabs(term_m1)))
 520                   {
 521                      //
 522                      // Now that we have started to converge, check to see if the value of
 523                      // this local sum is trivially small compared to the result.  If so
 524                      // abort this part of the series.
 525                      //
 526                      trivial_small_series_check = true;
 527                      Real d;
 528                      if (loop_scale > local_scaling)
 529                      {
 530                         int rescale = local_scaling - loop_scale;
 531                         if (rescale < tools::log_min_value<Real>())
 532                            d = 1;  // keep going
 533                         else
 534                            d = fabs(term / (result * exp(Real(rescale))));
 535                      }
 536                      else
 537                      {
 538                         int rescale = loop_scale - local_scaling;
 539                         if (rescale < tools::log_min_value<Real>())
 540                            d = 0;  // stop, underflow
 541                         else
 542                            d = fabs(term * exp(Real(rescale)) / result);
 543                      }
 544                      if (d < boost::math::policies::get_epsilon<Real, Policy>())
 545                         break;
 546                   }
 547
 548                   //std::cout << "k = " << k << " result = " << result << " abs_result = " << abs_result << std::endl;
 549                   if (fabs(loop_result) >= upper_limit)
 550                   {
 551                      loop_result /= scaling_factor;
 552                      loop_abs_result /= scaling_factor;
 553                      term /= scaling_factor;
 554                      loop_scale += log_scaling_factor;
 555                   }
 556                   if (fabs(loop_result) < lower_limit)
 557                   {
 558                      loop_result *= scaling_factor;
 559                      loop_abs_result *= scaling_factor;
 560                      term *= scaling_factor;
 561                      loop_scale -= log_scaling_factor;
 562                   }
 563                   diff = fabs(term / loop_result);
 564                } while (!termination(s - k) && ((diff > boost::math::policies::get_epsilon<Real, Policy>()) || (fabs(term) > fabs(term_m1))));
 565
 566                //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;
 567                //
 568                // We now need to combine the results of the first series summation with whatever
 569                // local results we have now.  First though, rescale abs_result by loop_error_scale
 570                // to factor in the error in the pochhammer terms at the start of this block:
 571                //
 572                loop_abs_result += loop_error_scale * fabs(loop_result);
 573                //
 574                if (loop_scale > local_scaling)
 575                {
 576                   //
 577                   // Need to shrink previous result:
 578                   //
 579                   int rescale = local_scaling - loop_scale;
 580                   local_scaling = loop_scale;
 581                   log_scale -= rescale;
 582                   Real ex = exp(Real(rescale));
 583                   result *= ex;
 584                   abs_result *= ex;
 585                   result += loop_result;
 586                   abs_result += loop_abs_result;
 587                }
 588                else if (local_scaling > loop_scale)
 589                {
 590                   //
 591                   // Need to shrink local result:
 592                   //
 593                   int rescale = loop_scale - local_scaling;
 594                   Real ex = exp(Real(rescale));
 595                   loop_result *= ex;
 596                   loop_abs_result *= ex;
 597                   result += loop_result;
 598                   abs_result += loop_abs_result;
 599                }
 600                else
 601                {
 602                   result += loop_result;
 603                   abs_result += loop_abs_result;
 604                }
 605                //
 606                // Reset k to the largest k we reached
 607                //
 608                k = next_backstop;
 609            }
 610         }
 611
 612         return std::make_pair(result, abs_result);
 613      }
 614
 615      struct iteration_terminator
 616      {
 617         iteration_terminator(boost::uintmax_t i) : m(i) {}
 618
 619         bool operator()(boost::uintmax_t v) const { return v >= m; }
 620
 621         boost::uintmax_t m;
 622      };
 623
 624      template <class Seq, class Real, class Policy>
 625      Real hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol, int& log_scale)
 626      {
 627         BOOST_MATH_STD_USING
 628         iteration_terminator term(boost::math::policies::get_max_series_iterations<Policy>());
 629         std::pair<Real, Real> result = hypergeometric_pFq_checked_series_impl(aj, bj, z, pol, term, log_scale);
 630         //
 631         // Check to see how many digits we've lost, if it's more than half, raise an evaluation error -
 632         // this is an entirely arbitrary cut off, but not unreasonable.
 633         //
 634         if (result.second * sqrt(boost::math::policies::get_epsilon<Real, Policy>()) > abs(result.first))
 635         {
 636            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);
 637         }
 638         return result.first;
 639      }
 640
 641      template <class Real, class Policy>
 642      inline Real hypergeometric_1F1_checked_series_impl(const Real& a, const Real& b, const Real& z, const Policy& pol, int& log_scale)
 643      {
 644         boost::array<Real, 1> aj = { a };
 645         boost::array<Real, 1> bj = { b };
 646         return hypergeometric_pFq_checked_series_impl(aj, bj, z, pol, log_scale);
 647      }
 648
 649   } } } // namespaces
 650
 651 #endif // BOOST_HYPERGEOMETRIC_PFQ_SERIES_HPP_