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)
7 #ifndef BOOST_HYPERGEOMETRIC_PFQ_SERIES_HPP_
8 #define BOOST_HYPERGEOMETRIC_PFQ_SERIES_HPP_
10 #ifndef BOOST_MATH_PFQ_MAX_B_TERMS
11 # define BOOST_MATH_PFQ_MAX_B_TERMS 5
16 #include <boost/math/special_functions/detail/hypergeometric_series.hpp>
18 namespace boost { namespace math { namespace detail {
20 template <class Seq, class Real>
21 unsigned set_crossover_locations(const Seq& aj, const Seq& bj, const Real& z, unsigned int* crossover_locations)
26 if(aj.size() == 1 && bj.size() == 1)
29 // For 1F1 we can work out where the peaks in the series occur,
30 // which is to say when:
32 // (a + k)z / (k(b + k)) == +-1
34 // Then we are at either a maxima or a minima in the series, and the
35 // last such point must be a maxima since the series is globally convergent.
36 // Potentially then we are solving 2 quadratic equations and have up to 4
37 // solutions, any solutions which are complex or negative are discarded,
38 // leaving us with 4 options:
40 // 0 solutions: The series is directly convergent.
41 // 1 solution : The series diverges to a maxima before converging.
42 // 2 solutions: The series is initially convergent, followed by divergence to a maxima before final convergence.
43 // 3 solutions: The series diverges to a maxima, converges to a minima before diverging again to a second maxima before final convergence.
44 // 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.
46 // The first 2 situations are adequately handled by direct series evaluation, while the 2,3 and 4 solutions are not.
50 Real sq = 4 * a * z + b * b - 2 * b * z + z * z;
53 Real t = (-sqrt(sq) - b + z) / 2;
56 crossover_locations[N_terms] = itrunc(t);
59 t = (sqrt(sq) - b + z) / 2;
62 crossover_locations[N_terms] = itrunc(t);
66 sq = -4 * a * z + b * b + 2 * b * z + z * z;
69 Real t = (-sqrt(sq) - b - z) / 2;
72 crossover_locations[N_terms] = itrunc(t);
75 t = (sqrt(sq) - b - z) / 2;
78 crossover_locations[N_terms] = itrunc(t);
82 std::sort(crossover_locations, crossover_locations + N_terms, std::less<Real>());
84 // Now we need to discard every other terms, as these are the minima:
92 crossover_locations[0] = crossover_locations[1];
96 crossover_locations[1] = crossover_locations[2];
100 crossover_locations[0] = crossover_locations[1];
101 crossover_locations[1] = crossover_locations[3];
109 for (auto bi = bj.begin(); bi != bj.end(); ++bi, ++n)
111 crossover_locations[n] = *bi >= 0 ? 0 : itrunc(-*bi) + 1;
113 std::sort(crossover_locations, crossover_locations + bj.size(), std::less<Real>());
114 N_terms = (unsigned)bj.size();
119 template <class Seq, class Real, class Policy, class Terminal>
120 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)
127 Real tol = boost::math::policies::get_epsilon<Real, Policy>();
128 std::uintmax_t k = 0;
129 Real upper_limit(sqrt(boost::math::tools::max_value<Real>())), diff;
130 Real lower_limit(1 / upper_limit);
131 long long log_scaling_factor = lltrunc(boost::math::tools::log_max_value<Real>()) - 2;
132 Real scaling_factor = exp(Real(log_scaling_factor));
134 long long local_scaling = 0;
136 if ((aj.size() == 1) && (bj.size() == 0))
140 if ((z > 0) && (floor(*aj.begin()) != *aj.begin()))
142 Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p == 1 and q == 0 and |z| > 1, result is imaginary", z, pol);
143 return std::make_pair(r, r);
145 std::pair<Real, Real> r = hypergeometric_pFq_checked_series_impl(aj, bj, Real(1 / z), pol, termination, log_scale);
146 Real mul = pow(-z, -*aj.begin());
153 if (aj.size() > bj.size())
155 if (aj.size() == bj.size() + 1)
159 Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p == q+1 and |z| > 1, series does not converge", z, pol);
160 return std::make_pair(r, r);
165 for (auto i = bj.begin(); i != bj.end(); ++i)
167 for (auto i = aj.begin(); i != aj.end(); ++i)
169 if ((z == 1) && (s <= 0))
171 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);
172 return std::make_pair(r, r);
174 if ((z == -1) && (s <= -1))
176 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);
177 return std::make_pair(r, r);
183 Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p > q+1, series does not converge", z, pol);
184 return std::make_pair(r, r);
188 while (!termination(k))
190 for (auto ai = aj.begin(); ai != aj.end(); ++ai)
196 // There is a negative integer in the aj's:
197 return std::make_pair(result, abs_result);
199 for (auto bi = bj.begin(); bi != bj.end(); ++bi)
203 // The series is undefined:
204 result = boost::math::policies::raise_domain_error("boost::math::hypergeometric_pFq<%1%>", "One of the b values was the negative integer %1%", *bi, pol);
205 return std::make_pair(result, result);
212 //std::cout << k << " " << *bj.begin() + k << " " << result << " " << term << /*" " << term_at_k(*aj.begin(), *bj.begin(), z, k, pol) <<*/ std::endl;
214 abs_result += abs(term);
215 //std::cout << "k = " << k << " term = " << term * exp(log_scale) << " result = " << result * exp(log_scale) << " abs_result = " << abs_result * exp(log_scale) << std::endl;
220 if (fabs(abs_result) >= upper_limit)
222 abs_result /= scaling_factor;
223 result /= scaling_factor;
224 term /= scaling_factor;
225 log_scale += log_scaling_factor;
226 local_scaling += log_scaling_factor;
228 if (fabs(abs_result) < lower_limit)
230 abs_result *= scaling_factor;
231 result *= scaling_factor;
232 term *= scaling_factor;
233 log_scale -= log_scaling_factor;
234 local_scaling -= log_scaling_factor;
237 if ((abs(result * tol) > abs(term)) && (abs(term0) > abs(term)))
239 if (abs_result * tol > abs(result))
241 // We have no correct bits in the result... just give up!
242 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);
243 return std::make_pair(result, result);
247 //std::cout << "result = " << result << std::endl;
248 //std::cout << "local_scaling = " << local_scaling << std::endl;
249 //std::cout << "Norm result = " << std::setprecision(35) << boost::multiprecision::mpfr_float_50(result) * exp(boost::multiprecision::mpfr_float_50(local_scaling)) << std::endl;
251 // We have to be careful when one of the b's crosses the origin:
253 if(bj.size() > BOOST_MATH_PFQ_MAX_B_TERMS)
254 policies::raise_domain_error<Real>("boost::math::hypergeometric_pFq<%1%>(Seq, Seq, %1%)",
255 "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%.",
256 Real(bj.size()), pol);
258 unsigned crossover_locations[BOOST_MATH_PFQ_MAX_B_TERMS];
260 unsigned N_crossovers = set_crossover_locations(aj, bj, z, crossover_locations);
262 bool terminate = false; // Set to true if one of the a's passes through the origin and terminates the series.
264 for (unsigned n = 0; n < N_crossovers; ++n)
266 if (k < crossover_locations[n])
268 for (auto ai = aj.begin(); ai != aj.end(); ++ai)
270 if ((*ai < 0) && (floor(*ai) == *ai) && (*ai > crossover_locations[n]))
271 return std::make_pair(result, abs_result); // b's will never cross the origin!
276 Real loop_result = 0;
277 Real loop_abs_result = 0;
278 long long loop_scale = 0;
280 // loop_error_scale will be used to increase the size of the error
281 // estimate (absolute sum), based on the errors inherent in calculating
282 // the pochhammer symbols.
284 Real loop_error_scale = 0;
285 //boost::multiprecision::mpfi_float err_est = 0;
287 // b hasn't crossed the origin yet and the series may spring back into life at that point
288 // so we need to jump forward to that term and then evaluate forwards and backwards from there:
290 unsigned s = crossover_locations[n];
291 std::uintmax_t backstop = k;
292 long long s1(1), s2(1);
294 for (auto ai = aj.begin(); ai != aj.end(); ++ai)
296 if ((floor(*ai) == *ai) && (*ai < 0) && (-*ai <= s))
298 // One of the a terms has passed through zero and terminated the series:
305 Real p = log_pochhammer(*ai, s, pol, &ls);
308 loop_error_scale = (std::max)(p, loop_error_scale);
309 //err_est += boost::multiprecision::mpfi_float(p);
312 //std::cout << "term = " << term << std::endl;
315 for (auto bi = bj.begin(); bi != bj.end(); ++bi)
318 Real p = log_pochhammer(*bi, s, pol, &ls);
321 loop_error_scale = (std::max)(p, loop_error_scale);
322 //err_est -= boost::multiprecision::mpfi_float(p);
324 //std::cout << "term = " << term << std::endl;
325 Real p = lgamma(Real(s + 1), pol);
327 loop_error_scale = (std::max)(p, loop_error_scale);
328 //err_est -= boost::multiprecision::mpfi_float(p);
329 p = s * log(fabs(z));
331 loop_error_scale = (std::max)(p, loop_error_scale);
332 //err_est += boost::multiprecision::mpfi_float(p);
333 //err_est = exp(err_est);
334 //std::cout << err_est << std::endl;
336 // Convert loop_error scale to the absolute error
337 // in term after exp is applied:
339 loop_error_scale *= tools::epsilon<Real>();
341 // Convert to relative error after exp:
343 loop_error_scale = fabs(expm1(loop_error_scale, pol));
345 // Convert to multiplier for the error term:
347 loop_error_scale /= tools::epsilon<Real>();
350 s1 *= (s & 1 ? -1 : 1);
352 if (term <= tools::log_min_value<Real>())
354 // rescale if we can:
355 long long scale = lltrunc(floor(term - tools::log_min_value<Real>()) - 2);
361 int scale = itrunc(floor(term));
365 //std::cout << "term = " << term << std::endl;
366 term = s1 * s2 * exp(term);
367 //std::cout << "term = " << term << std::endl;
368 //std::cout << "loop_scale = " << loop_scale << std::endl;
371 long long saved_loop_scale = loop_scale;
372 bool terms_are_growing = true;
373 bool trivial_small_series_check = false;
377 loop_abs_result += fabs(term);
378 //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;
379 if (fabs(loop_result) >= upper_limit)
381 loop_result /= scaling_factor;
382 loop_abs_result /= scaling_factor;
383 term /= scaling_factor;
384 loop_scale += log_scaling_factor;
386 if (fabs(loop_result) < lower_limit)
388 loop_result *= scaling_factor;
389 loop_abs_result *= scaling_factor;
390 term *= scaling_factor;
391 loop_scale -= log_scaling_factor;
394 for (auto ai = aj.begin(); ai != aj.end(); ++ai)
400 // There is a negative integer in the aj's:
401 return std::make_pair(result, abs_result);
403 for (auto bi = bj.begin(); bi != bj.end(); ++bi)
407 // The series is undefined:
408 result = boost::math::policies::raise_domain_error("boost::math::hypergeometric_pFq<%1%>", "One of the b values was the negative integer %1%", *bi, pol);
409 return std::make_pair(result, result);
416 diff = fabs(term / loop_result);
417 terms_are_growing = fabs(term) > fabs(term_m1);
418 if (!trivial_small_series_check && !terms_are_growing)
421 // Now that we have started to converge, check to see if the value of
422 // this local sum is trivially small compared to the result. If so
423 // abort this part of the series.
425 trivial_small_series_check = true;
427 if (loop_scale > local_scaling)
429 long long rescale = local_scaling - loop_scale;
430 if (rescale < tools::log_min_value<Real>())
431 d = 1; // arbitrary value, we want to keep going
433 d = fabs(term / (result * exp(Real(rescale))));
437 long long rescale = loop_scale - local_scaling;
438 if (rescale < tools::log_min_value<Real>())
439 d = 0; // terminate this loop
441 d = fabs(term * exp(Real(rescale)) / result);
443 if (d < boost::math::policies::get_epsilon<Real, Policy>())
446 } while (!termination(k - s) && ((diff > boost::math::policies::get_epsilon<Real, Policy>()) || terms_are_growing));
448 //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;
450 // We now need to combine the results of the first series summation with whatever
451 // local results we have now. First though, rescale abs_result by loop_error_scale
452 // to factor in the error in the pochhammer terms at the start of this block:
454 std::uintmax_t next_backstop = k;
455 loop_abs_result += loop_error_scale * fabs(loop_result);
456 if (loop_scale > local_scaling)
459 // Need to shrink previous result:
461 long long rescale = local_scaling - loop_scale;
462 local_scaling = loop_scale;
463 log_scale -= rescale;
464 Real ex = exp(Real(rescale));
467 result += loop_result;
468 abs_result += loop_abs_result;
470 else if (local_scaling > loop_scale)
473 // Need to shrink local result:
475 long long rescale = loop_scale - local_scaling;
476 Real ex = exp(Real(rescale));
478 loop_abs_result *= ex;
479 result += loop_result;
480 abs_result += loop_abs_result;
484 result += loop_result;
485 abs_result += loop_abs_result;
488 // Now go backwards as well:
494 loop_scale = saved_loop_scale;
495 trivial_small_series_check = false;
502 for (auto ai = aj.begin(); ai != aj.end(); ++ai)
506 for (auto bi = bj.begin(); bi != bj.end(); ++bi)
510 // The series is undefined:
511 result = boost::math::policies::raise_domain_error("boost::math::hypergeometric_pFq<%1%>", "One of the b values was the negative integer %1%", *bi, pol);
512 return std::make_pair(result, result);
518 loop_abs_result += fabs(term);
520 if (!trivial_small_series_check && (fabs(term) < fabs(term_m1)))
523 // Now that we have started to converge, check to see if the value of
524 // this local sum is trivially small compared to the result. If so
525 // abort this part of the series.
527 trivial_small_series_check = true;
529 if (loop_scale > local_scaling)
531 long long rescale = local_scaling - loop_scale;
532 if (rescale < tools::log_min_value<Real>())
535 d = fabs(term / (result * exp(Real(rescale))));
539 long long rescale = loop_scale - local_scaling;
540 if (rescale < tools::log_min_value<Real>())
541 d = 0; // stop, underflow
543 d = fabs(term * exp(Real(rescale)) / result);
545 if (d < boost::math::policies::get_epsilon<Real, Policy>())
549 //std::cout << "k = " << k << " result = " << result << " abs_result = " << abs_result << std::endl;
550 if (fabs(loop_result) >= upper_limit)
552 loop_result /= scaling_factor;
553 loop_abs_result /= scaling_factor;
554 term /= scaling_factor;
555 loop_scale += log_scaling_factor;
557 if (fabs(loop_result) < lower_limit)
559 loop_result *= scaling_factor;
560 loop_abs_result *= scaling_factor;
561 term *= scaling_factor;
562 loop_scale -= log_scaling_factor;
564 diff = fabs(term / loop_result);
565 } while (!termination(s - k) && ((diff > boost::math::policies::get_epsilon<Real, Policy>()) || (fabs(term) > fabs(term_m1))));
567 //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;
569 // We now need to combine the results of the first series summation with whatever
570 // local results we have now. First though, rescale abs_result by loop_error_scale
571 // to factor in the error in the pochhammer terms at the start of this block:
573 loop_abs_result += loop_error_scale * fabs(loop_result);
575 if (loop_scale > local_scaling)
578 // Need to shrink previous result:
580 long long rescale = local_scaling - loop_scale;
581 local_scaling = loop_scale;
582 log_scale -= rescale;
583 Real ex = exp(Real(rescale));
586 result += loop_result;
587 abs_result += loop_abs_result;
589 else if (local_scaling > loop_scale)
592 // Need to shrink local result:
594 long long rescale = loop_scale - local_scaling;
595 Real ex = exp(Real(rescale));
597 loop_abs_result *= ex;
598 result += loop_result;
599 abs_result += loop_abs_result;
603 result += loop_result;
604 abs_result += loop_abs_result;
607 // Reset k to the largest k we reached
613 return std::make_pair(result, abs_result);
616 struct iteration_terminator
618 iteration_terminator(std::uintmax_t i) : m(i) {}
620 bool operator()(std::uintmax_t v) const { return v >= m; }
625 template <class Seq, class Real, class Policy>
626 Real hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol, long long& log_scale)
629 iteration_terminator term(boost::math::policies::get_max_series_iterations<Policy>());
630 std::pair<Real, Real> result = hypergeometric_pFq_checked_series_impl(aj, bj, z, pol, term, log_scale);
632 // Check to see how many digits we've lost, if it's more than half, raise an evaluation error -
633 // this is an entirely arbitrary cut off, but not unreasonable.
635 if (result.second * sqrt(boost::math::policies::get_epsilon<Real, Policy>()) > abs(result.first))
637 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);
642 template <class Real, class Policy>
643 inline Real hypergeometric_1F1_checked_series_impl(const Real& a, const Real& b, const Real& z, const Policy& pol, long long& log_scale)
645 std::array<Real, 1> aj = { a };
646 std::array<Real, 1> bj = { b };
647 return hypergeometric_pFq_checked_series_impl(aj, bj, z, pol, log_scale);
652 #endif // BOOST_HYPERGEOMETRIC_PFQ_SERIES_HPP_