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
14 #include <boost/array.hpp>
15 #include <boost/math/special_functions/detail/hypergeometric_series.hpp>
17 namespace boost { namespace math { namespace detail {
19 template <class Seq, class Real>
20 unsigned set_crossover_locations(const Seq& aj, const Seq& bj, const Real& z, unsigned int* crossover_locations)
25 if(aj.size() == 1 && bj.size() == 1)
28 // For 1F1 we can work out where the peaks in the series occur,
29 // which is to say when:
31 // (a + k)z / (k(b + k)) == +-1
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:
39 // 0 solutions: The series is directly convergent.
40 // 1 solution : The series diverges to a maxima before converging.
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.
45 // The first 2 situations are adequately handled by direct series evaluation, while the 2,3 and 4 solutions are not.
49 Real sq = 4 * a * z + b * b - 2 * b * z + z * z;
52 Real t = (-sqrt(sq) - b + z) / 2;
55 crossover_locations[N_terms] = itrunc(t);
58 t = (sqrt(sq) - b + z) / 2;
61 crossover_locations[N_terms] = itrunc(t);
65 sq = -4 * a * z + b * b + 2 * b * z + z * z;
68 Real t = (-sqrt(sq) - b - z) / 2;
71 crossover_locations[N_terms] = itrunc(t);
74 t = (sqrt(sq) - b - z) / 2;
77 crossover_locations[N_terms] = itrunc(t);
81 std::sort(crossover_locations, crossover_locations + N_terms, std::less<Real>());
83 // Now we need to discard every other terms, as these are the minima:
91 crossover_locations[0] = crossover_locations[1];
95 crossover_locations[1] = crossover_locations[2];
99 crossover_locations[0] = crossover_locations[1];
100 crossover_locations[1] = crossover_locations[3];
108 for (auto bi = bj.begin(); bi != bj.end(); ++bi, ++n)
110 crossover_locations[n] = *bi >= 0 ? 0 : itrunc(-*bi) + 1;
112 std::sort(crossover_locations, crossover_locations + bj.size(), std::less<Real>());
113 N_terms = (unsigned)bj.size();
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)
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));
133 int local_scaling = 0;
135 if ((aj.size() == 1) && (bj.size() == 0))
139 if ((z > 0) && (floor(*aj.begin()) != *aj.begin()))
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);
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());
152 if (aj.size() > bj.size())
154 if (aj.size() == bj.size() + 1)
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);
164 for (auto i = bj.begin(); i != bj.end(); ++i)
166 for (auto i = aj.begin(); i != aj.end(); ++i)
168 if ((z == 1) && (s <= 0))
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);
173 if ((z == -1) && (s <= -1))
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);
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);
187 while (!termination(k))
189 for (auto ai = aj.begin(); ai != aj.end(); ++ai)
195 // There is a negative integer in the aj's:
196 return std::make_pair(result, abs_result);
198 for (auto bi = bj.begin(); bi != bj.end(); ++bi)
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);
211 //std::cout << k << " " << *bj.begin() + k << " " << result << " " << term << /*" " << term_at_k(*aj.begin(), *bj.begin(), z, k, pol) <<*/ std::endl;
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;
219 if (fabs(abs_result) >= upper_limit)
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;
227 if (fabs(abs_result) < lower_limit)
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;
236 if ((abs(result * tol) > abs(term)) && (abs(term0) > abs(term)))
238 if (abs_result * tol > abs(result))
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);
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;
250 // We have to be careful when one of the b's crosses the origin:
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);
257 unsigned crossover_locations[BOOST_MATH_PFQ_MAX_B_TERMS];
259 unsigned N_crossovers = set_crossover_locations(aj, bj, z, crossover_locations);
261 bool terminate = false; // Set to true if one of the a's passes through the origin and terminates the series.
263 for (unsigned n = 0; n < N_crossovers; ++n)
265 if (k < crossover_locations[n])
267 for (auto ai = aj.begin(); ai != aj.end(); ++ai)
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!
275 Real loop_result = 0;
276 Real loop_abs_result = 0;
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.
283 Real loop_error_scale = 0;
284 //boost::multiprecision::mpfi_float err_est = 0;
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:
289 unsigned s = crossover_locations[n];
290 boost::uintmax_t backstop = k;
293 for (auto ai = aj.begin(); ai != aj.end(); ++ai)
295 if ((floor(*ai) == *ai) && (*ai < 0) && (-*ai <= s))
297 // One of the a terms has passed through zero and terminated the series:
304 Real p = log_pochhammer(*ai, s, pol, &ls);
307 loop_error_scale = (std::max)(p, loop_error_scale);
308 //err_est += boost::multiprecision::mpfi_float(p);
311 //std::cout << "term = " << term << std::endl;
314 for (auto bi = bj.begin(); bi != bj.end(); ++bi)
317 Real p = log_pochhammer(*bi, s, pol, &ls);
320 loop_error_scale = (std::max)(p, loop_error_scale);
321 //err_est -= boost::multiprecision::mpfi_float(p);
323 //std::cout << "term = " << term << std::endl;
324 Real p = lgamma(Real(s + 1), pol);
326 loop_error_scale = (std::max)(p, loop_error_scale);
327 //err_est -= boost::multiprecision::mpfi_float(p);
328 p = s * log(fabs(z));
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;
335 // Convert loop_error scale to the absolute error
336 // in term after exp is applied:
338 loop_error_scale *= tools::epsilon<Real>();
340 // Convert to relative error after exp:
342 loop_error_scale = fabs(expm1(loop_error_scale));
344 // Convert to multiplier for the error term:
346 loop_error_scale /= tools::epsilon<Real>();
349 s1 *= (s & 1 ? -1 : 1);
351 if (term <= tools::log_min_value<Real>())
353 // rescale if we can:
354 int scale = itrunc(floor(term - tools::log_min_value<Real>()) - 2);
360 int scale = itrunc(floor(term));
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;
370 int saved_loop_scale = loop_scale;
371 bool terms_are_growing = true;
372 bool trivial_small_series_check = false;
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)
380 loop_result /= scaling_factor;
381 loop_abs_result /= scaling_factor;
382 term /= scaling_factor;
383 loop_scale += log_scaling_factor;
385 if (fabs(loop_result) < lower_limit)
387 loop_result *= scaling_factor;
388 loop_abs_result *= scaling_factor;
389 term *= scaling_factor;
390 loop_scale -= log_scaling_factor;
393 for (auto ai = aj.begin(); ai != aj.end(); ++ai)
399 // There is a negative integer in the aj's:
400 return std::make_pair(result, abs_result);
402 for (auto bi = bj.begin(); bi != bj.end(); ++bi)
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);
415 diff = fabs(term / loop_result);
416 terms_are_growing = fabs(term) > fabs(term_m1);
417 if (!trivial_small_series_check && !terms_are_growing)
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.
424 trivial_small_series_check = true;
426 if (loop_scale > local_scaling)
428 int rescale = local_scaling - loop_scale;
429 if (rescale < tools::log_min_value<Real>())
430 d = 1; // arbitrary value, we want to keep going
432 d = fabs(term / (result * exp(Real(rescale))));
436 int rescale = loop_scale - local_scaling;
437 if (rescale < tools::log_min_value<Real>())
438 d = 0; // terminate this loop
440 d = fabs(term * exp(Real(rescale)) / result);
442 if (d < boost::math::policies::get_epsilon<Real, Policy>())
445 } while (!termination(k - s) && ((diff > boost::math::policies::get_epsilon<Real, Policy>()) || terms_are_growing));
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;
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:
453 boost::uintmax_t next_backstop = k;
454 loop_abs_result += loop_error_scale * fabs(loop_result);
455 if (loop_scale > local_scaling)
458 // Need to shrink previous result:
460 int rescale = local_scaling - loop_scale;
461 local_scaling = loop_scale;
462 log_scale -= rescale;
463 Real ex = exp(Real(rescale));
466 result += loop_result;
467 abs_result += loop_abs_result;
469 else if (local_scaling > loop_scale)
472 // Need to shrink local result:
474 int rescale = loop_scale - local_scaling;
475 Real ex = exp(Real(rescale));
477 loop_abs_result *= ex;
478 result += loop_result;
479 abs_result += loop_abs_result;
483 result += loop_result;
484 abs_result += loop_abs_result;
487 // Now go backwards as well:
493 loop_scale = saved_loop_scale;
494 trivial_small_series_check = false;
501 for (auto ai = aj.begin(); ai != aj.end(); ++ai)
505 for (auto bi = bj.begin(); bi != bj.end(); ++bi)
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);
517 loop_abs_result += fabs(term);
519 if (!trivial_small_series_check && (fabs(term) < fabs(term_m1)))
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.
526 trivial_small_series_check = true;
528 if (loop_scale > local_scaling)
530 int rescale = local_scaling - loop_scale;
531 if (rescale < tools::log_min_value<Real>())
534 d = fabs(term / (result * exp(Real(rescale))));
538 int rescale = loop_scale - local_scaling;
539 if (rescale < tools::log_min_value<Real>())
540 d = 0; // stop, underflow
542 d = fabs(term * exp(Real(rescale)) / result);
544 if (d < boost::math::policies::get_epsilon<Real, Policy>())
548 //std::cout << "k = " << k << " result = " << result << " abs_result = " << abs_result << std::endl;
549 if (fabs(loop_result) >= upper_limit)
551 loop_result /= scaling_factor;
552 loop_abs_result /= scaling_factor;
553 term /= scaling_factor;
554 loop_scale += log_scaling_factor;
556 if (fabs(loop_result) < lower_limit)
558 loop_result *= scaling_factor;
559 loop_abs_result *= scaling_factor;
560 term *= scaling_factor;
561 loop_scale -= log_scaling_factor;
563 diff = fabs(term / loop_result);
564 } while (!termination(s - k) && ((diff > boost::math::policies::get_epsilon<Real, Policy>()) || (fabs(term) > fabs(term_m1))));
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;
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:
572 loop_abs_result += loop_error_scale * fabs(loop_result);
574 if (loop_scale > local_scaling)
577 // Need to shrink previous result:
579 int rescale = local_scaling - loop_scale;
580 local_scaling = loop_scale;
581 log_scale -= rescale;
582 Real ex = exp(Real(rescale));
585 result += loop_result;
586 abs_result += loop_abs_result;
588 else if (local_scaling > loop_scale)
591 // Need to shrink local result:
593 int rescale = loop_scale - local_scaling;
594 Real ex = exp(Real(rescale));
596 loop_abs_result *= ex;
597 result += loop_result;
598 abs_result += loop_abs_result;
602 result += loop_result;
603 abs_result += loop_abs_result;
606 // Reset k to the largest k we reached
612 return std::make_pair(result, abs_result);
615 struct iteration_terminator
617 iteration_terminator(boost::uintmax_t i) : m(i) {}
619 bool operator()(boost::uintmax_t v) const { return v >= m; }
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)
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);
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.
634 if (result.second * sqrt(boost::math::policies::get_epsilon<Real, Policy>()) > abs(result.first))
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);
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)
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);
651 #endif // BOOST_HYPERGEOMETRIC_PFQ_SERIES_HPP_