123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671 |
- #ifndef BOOST_HYPERGEOMETRIC_PFQ_SERIES_HPP_
- #define BOOST_HYPERGEOMETRIC_PFQ_SERIES_HPP_
- #ifndef BOOST_MATH_PFQ_MAX_B_TERMS
- # define BOOST_MATH_PFQ_MAX_B_TERMS 5
- #endif
- #include <array>
- #include <cstdint>
- #include <boost/math/special_functions/gamma.hpp>
- #include <boost/math/special_functions/expm1.hpp>
- #include <boost/math/special_functions/detail/hypergeometric_series.hpp>
- namespace boost { namespace math { namespace detail {
- template <class Seq, class Real>
- unsigned set_crossover_locations(const Seq& aj, const Seq& bj, const Real& z, unsigned int* crossover_locations)
- {
- BOOST_MATH_STD_USING
- unsigned N_terms = 0;
- if(aj.size() == 1 && bj.size() == 1)
- {
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- Real a = *aj.begin();
- Real b = *bj.begin();
- Real sq = 4 * a * z + b * b - 2 * b * z + z * z;
- if (sq >= 0)
- {
- Real t = (-sqrt(sq) - b + z) / 2;
- if (t >= 0)
- {
- crossover_locations[N_terms] = itrunc(t);
- ++N_terms;
- }
- t = (sqrt(sq) - b + z) / 2;
- if (t >= 0)
- {
- crossover_locations[N_terms] = itrunc(t);
- ++N_terms;
- }
- }
- sq = -4 * a * z + b * b + 2 * b * z + z * z;
- if (sq >= 0)
- {
- Real t = (-sqrt(sq) - b - z) / 2;
- if (t >= 0)
- {
- crossover_locations[N_terms] = itrunc(t);
- ++N_terms;
- }
- t = (sqrt(sq) - b - z) / 2;
- if (t >= 0)
- {
- crossover_locations[N_terms] = itrunc(t);
- ++N_terms;
- }
- }
- std::sort(crossover_locations, crossover_locations + N_terms, std::less<Real>());
-
-
-
- switch (N_terms)
- {
- case 0:
- case 1:
- break;
- case 2:
- crossover_locations[0] = crossover_locations[1];
- --N_terms;
- break;
- case 3:
- crossover_locations[1] = crossover_locations[2];
- --N_terms;
- break;
- case 4:
- crossover_locations[0] = crossover_locations[1];
- crossover_locations[1] = crossover_locations[3];
- N_terms -= 2;
- break;
- }
- }
- else
- {
- unsigned n = 0;
- for (auto bi = bj.begin(); bi != bj.end(); ++bi, ++n)
- {
- crossover_locations[n] = *bi >= 0 ? 0 : itrunc(-*bi) + 1;
- }
- std::sort(crossover_locations, crossover_locations + bj.size(), std::less<Real>());
- N_terms = (unsigned)bj.size();
- }
- return N_terms;
- }
- template <class Seq, class Real, class Policy, class Terminal>
- 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)
- {
- BOOST_MATH_STD_USING
- Real result = 1;
- Real abs_result = 1;
- Real term = 1;
- Real term0 = 0;
- Real tol = boost::math::policies::get_epsilon<Real, Policy>();
- std::uintmax_t k = 0;
- Real upper_limit(sqrt(boost::math::tools::max_value<Real>())), diff;
- Real lower_limit(1 / upper_limit);
- long long log_scaling_factor = lltrunc(boost::math::tools::log_max_value<Real>()) - 2;
- Real scaling_factor = exp(Real(log_scaling_factor));
- Real term_m1;
- long long local_scaling = 0;
- bool have_no_correct_bits = false;
- if ((aj.size() == 1) && (bj.size() == 0))
- {
- if (fabs(z) > 1)
- {
- if ((z > 0) && (floor(*aj.begin()) != *aj.begin()))
- {
- Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p == 1 and q == 0 and |z| > 1, result is imaginary", z, pol);
- return std::make_pair(r, r);
- }
- std::pair<Real, Real> r = hypergeometric_pFq_checked_series_impl(aj, bj, Real(1 / z), pol, termination, log_scale);
-
- #if (defined(__GNUC__) && __GNUC__ == 13)
- Real mul = pow(-z, Real(-*aj.begin()));
- #else
- Real mul = pow(-z, -*aj.begin());
- #endif
-
- r.first *= mul;
- r.second *= mul;
- return r;
- }
- }
- if (aj.size() > bj.size())
- {
- if (aj.size() == bj.size() + 1)
- {
- if (fabs(z) > 1)
- {
- Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p == q+1 and |z| > 1, series does not converge", z, pol);
- return std::make_pair(r, r);
- }
- if (fabs(z) == 1)
- {
- Real s = 0;
- for (auto i = bj.begin(); i != bj.end(); ++i)
- s += *i;
- for (auto i = aj.begin(); i != aj.end(); ++i)
- s -= *i;
- if ((z == 1) && (s <= 0))
- {
- 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);
- return std::make_pair(r, r);
- }
- if ((z == -1) && (s <= -1))
- {
- 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);
- return std::make_pair(r, r);
- }
- }
- }
- else
- {
- Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p > q+1, series does not converge", z, pol);
- return std::make_pair(r, r);
- }
- }
- while (!termination(k))
- {
- for (auto ai = aj.begin(); ai != aj.end(); ++ai)
- {
- term *= *ai + k;
- }
- if (term == 0)
- {
-
- return std::make_pair(result, abs_result);
- }
- for (auto bi = bj.begin(); bi != bj.end(); ++bi)
- {
- if (*bi + k == 0)
- {
-
- result = boost::math::policies::raise_domain_error("boost::math::hypergeometric_pFq<%1%>", "One of the b values was the negative integer %1%", *bi, pol);
- return std::make_pair(result, result);
- }
- term /= *bi + k;
- }
- term *= z;
- ++k;
- term /= k;
-
- result += term;
- abs_result += abs(term);
-
-
-
-
- if (fabs(abs_result) >= upper_limit)
- {
- abs_result /= scaling_factor;
- result /= scaling_factor;
- term /= scaling_factor;
- log_scale += log_scaling_factor;
- local_scaling += log_scaling_factor;
- }
- if (fabs(abs_result) < lower_limit)
- {
- abs_result *= scaling_factor;
- result *= scaling_factor;
- term *= scaling_factor;
- log_scale -= log_scaling_factor;
- local_scaling -= log_scaling_factor;
- }
- if ((abs(result * tol) > abs(term)) && (abs(term0) > abs(term)))
- break;
- if (abs_result * tol > abs(result))
- {
-
-
-
- if (have_no_correct_bits)
- {
-
- 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);
- return std::make_pair(result, result);
- }
- else
- have_no_correct_bits = true;
- }
- else
- have_no_correct_bits = false;
- term0 = term;
- }
-
-
-
-
-
-
- if(bj.size() > BOOST_MATH_PFQ_MAX_B_TERMS)
- policies::raise_domain_error<Real>("boost::math::hypergeometric_pFq<%1%>(Seq, Seq, %1%)",
- "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%.",
- Real(bj.size()), pol);
- unsigned crossover_locations[BOOST_MATH_PFQ_MAX_B_TERMS];
- unsigned N_crossovers = set_crossover_locations(aj, bj, z, crossover_locations);
- bool terminate = false;
- for (unsigned n = 0; n < N_crossovers; ++n)
- {
- if (k < crossover_locations[n])
- {
- for (auto ai = aj.begin(); ai != aj.end(); ++ai)
- {
- if ((*ai < 0) && (floor(*ai) == *ai) && (*ai > static_cast<decltype(*ai)>(crossover_locations[n])))
- return std::make_pair(result, abs_result);
- }
-
-
-
- Real loop_result = 0;
- Real loop_abs_result = 0;
- long long loop_scale = 0;
-
-
-
-
-
- Real loop_error_scale = 0;
-
-
-
-
-
- unsigned s = crossover_locations[n];
- std::uintmax_t backstop = k;
- long long s1(1), s2(1);
- term = 0;
- for (auto ai = aj.begin(); ai != aj.end(); ++ai)
- {
- if ((floor(*ai) == *ai) && (*ai < 0) && (-*ai <= static_cast<decltype(*ai)>(s)))
- {
-
- terminate = true;
- break;
- }
- else
- {
- int ls = 1;
- Real p = log_pochhammer(*ai, s, pol, &ls);
- s1 *= ls;
- term += p;
- loop_error_scale = (std::max)(p, loop_error_scale);
-
- }
- }
-
- if (terminate)
- break;
- for (auto bi = bj.begin(); bi != bj.end(); ++bi)
- {
- int ls = 1;
- Real p = log_pochhammer(*bi, s, pol, &ls);
- s2 *= ls;
- term -= p;
- loop_error_scale = (std::max)(p, loop_error_scale);
-
- }
-
- Real p = lgamma(Real(s + 1), pol);
- term -= p;
- loop_error_scale = (std::max)(p, loop_error_scale);
-
- p = s * log(fabs(z));
- term += p;
- loop_error_scale = (std::max)(p, loop_error_scale);
-
-
-
-
-
-
-
- loop_error_scale *= tools::epsilon<Real>();
-
-
-
- loop_error_scale = fabs(expm1(loop_error_scale, pol));
-
-
-
- loop_error_scale /= tools::epsilon<Real>();
- if (z < 0)
- s1 *= (s & 1 ? -1 : 1);
- if (term <= tools::log_min_value<Real>())
- {
-
- long long scale = lltrunc(floor(term - tools::log_min_value<Real>()) - 2);
- term -= scale;
- loop_scale += scale;
- }
- if (term > 10)
- {
- int scale = itrunc(floor(term));
- term -= scale;
- loop_scale += scale;
- }
-
- term = s1 * s2 * exp(term);
-
-
- k = s;
- term0 = term;
- long long saved_loop_scale = loop_scale;
- bool terms_are_growing = true;
- bool trivial_small_series_check = false;
- do
- {
- loop_result += term;
- loop_abs_result += fabs(term);
-
- if (fabs(loop_result) >= upper_limit)
- {
- loop_result /= scaling_factor;
- loop_abs_result /= scaling_factor;
- term /= scaling_factor;
- loop_scale += log_scaling_factor;
- }
- if (fabs(loop_result) < lower_limit)
- {
- loop_result *= scaling_factor;
- loop_abs_result *= scaling_factor;
- term *= scaling_factor;
- loop_scale -= log_scaling_factor;
- }
- term_m1 = term;
- for (auto ai = aj.begin(); ai != aj.end(); ++ai)
- {
- term *= *ai + k;
- }
- if (term == 0)
- {
-
- return std::make_pair(result, abs_result);
- }
- for (auto bi = bj.begin(); bi != bj.end(); ++bi)
- {
- if (*bi + k == 0)
- {
-
- result = boost::math::policies::raise_domain_error("boost::math::hypergeometric_pFq<%1%>", "One of the b values was the negative integer %1%", *bi, pol);
- return std::make_pair(result, result);
- }
- term /= *bi + k;
- }
- term *= z / (k + 1);
- ++k;
- diff = fabs(term / loop_result);
- terms_are_growing = fabs(term) > fabs(term_m1);
- if (!trivial_small_series_check && !terms_are_growing)
- {
-
-
-
-
-
- trivial_small_series_check = true;
- Real d;
- if (loop_scale > local_scaling)
- {
- long long rescale = local_scaling - loop_scale;
- if (rescale < tools::log_min_value<Real>())
- d = 1;
- else
- d = fabs(term / (result * exp(Real(rescale))));
- }
- else
- {
- long long rescale = loop_scale - local_scaling;
- if (rescale < tools::log_min_value<Real>())
- d = 0;
- else
- d = fabs(term * exp(Real(rescale)) / result);
- }
- if (d < boost::math::policies::get_epsilon<Real, Policy>())
- break;
- }
- } while (!termination(k - s) && ((diff > boost::math::policies::get_epsilon<Real, Policy>()) || terms_are_growing));
-
-
-
-
-
-
- std::uintmax_t next_backstop = k;
- loop_abs_result += loop_error_scale * fabs(loop_result);
- if (loop_scale > local_scaling)
- {
-
-
-
- long long rescale = local_scaling - loop_scale;
- local_scaling = loop_scale;
- log_scale -= rescale;
- Real ex = exp(Real(rescale));
- result *= ex;
- abs_result *= ex;
- result += loop_result;
- abs_result += loop_abs_result;
- }
- else if (local_scaling > loop_scale)
- {
-
-
-
- long long rescale = loop_scale - local_scaling;
- Real ex = exp(Real(rescale));
- loop_result *= ex;
- loop_abs_result *= ex;
- result += loop_result;
- abs_result += loop_abs_result;
- }
- else
- {
- result += loop_result;
- abs_result += loop_abs_result;
- }
-
-
-
- k = s;
- term = term0;
- loop_result = 0;
- loop_abs_result = 0;
- loop_scale = saved_loop_scale;
- trivial_small_series_check = false;
- do
- {
- --k;
- if (k == backstop)
- break;
- term_m1 = term;
- for (auto ai = aj.begin(); ai != aj.end(); ++ai)
- {
- term /= *ai + k;
- }
- for (auto bi = bj.begin(); bi != bj.end(); ++bi)
- {
- if (*bi + k == 0)
- {
-
- result = boost::math::policies::raise_domain_error("boost::math::hypergeometric_pFq<%1%>", "One of the b values was the negative integer %1%", *bi, pol);
- return std::make_pair(result, result);
- }
- term *= *bi + k;
- }
- term *= (k + 1) / z;
- loop_result += term;
- loop_abs_result += fabs(term);
- if (!trivial_small_series_check && (fabs(term) < fabs(term_m1)))
- {
-
-
-
-
-
- trivial_small_series_check = true;
- Real d;
- if (loop_scale > local_scaling)
- {
- long long rescale = local_scaling - loop_scale;
- if (rescale < tools::log_min_value<Real>())
- d = 1;
- else
- d = fabs(term / (result * exp(Real(rescale))));
- }
- else
- {
- long long rescale = loop_scale - local_scaling;
- if (rescale < tools::log_min_value<Real>())
- d = 0;
- else
- d = fabs(term * exp(Real(rescale)) / result);
- }
- if (d < boost::math::policies::get_epsilon<Real, Policy>())
- break;
- }
-
- if (fabs(loop_result) >= upper_limit)
- {
- loop_result /= scaling_factor;
- loop_abs_result /= scaling_factor;
- term /= scaling_factor;
- loop_scale += log_scaling_factor;
- }
- if (fabs(loop_result) < lower_limit)
- {
- loop_result *= scaling_factor;
- loop_abs_result *= scaling_factor;
- term *= scaling_factor;
- loop_scale -= log_scaling_factor;
- }
- diff = fabs(term / loop_result);
- } while (!termination(s - k) && ((diff > boost::math::policies::get_epsilon<Real, Policy>()) || (fabs(term) > fabs(term_m1))));
-
-
-
-
-
-
- loop_abs_result += loop_error_scale * fabs(loop_result);
-
- if (loop_scale > local_scaling)
- {
-
-
-
- long long rescale = local_scaling - loop_scale;
- local_scaling = loop_scale;
- log_scale -= rescale;
- Real ex = exp(Real(rescale));
- result *= ex;
- abs_result *= ex;
- result += loop_result;
- abs_result += loop_abs_result;
- }
- else if (local_scaling > loop_scale)
- {
-
-
-
- long long rescale = loop_scale - local_scaling;
- Real ex = exp(Real(rescale));
- loop_result *= ex;
- loop_abs_result *= ex;
- result += loop_result;
- abs_result += loop_abs_result;
- }
- else
- {
- result += loop_result;
- abs_result += loop_abs_result;
- }
-
-
-
- k = next_backstop;
- }
- }
- return std::make_pair(result, abs_result);
- }
- struct iteration_terminator
- {
- iteration_terminator(std::uintmax_t i) : m(i) {}
- bool operator()(std::uintmax_t v) const { return v >= m; }
- std::uintmax_t m;
- };
- template <class Seq, class Real, class Policy>
- Real hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol, long long& log_scale)
- {
- BOOST_MATH_STD_USING
- iteration_terminator term(boost::math::policies::get_max_series_iterations<Policy>());
- std::pair<Real, Real> result = hypergeometric_pFq_checked_series_impl(aj, bj, z, pol, term, log_scale);
-
-
-
-
- if (result.second * sqrt(boost::math::policies::get_epsilon<Real, Policy>()) > abs(result.first))
- {
- 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);
- }
- return result.first;
- }
- template <class Real, class Policy>
- inline Real hypergeometric_1F1_checked_series_impl(const Real& a, const Real& b, const Real& z, const Policy& pol, long long& log_scale)
- {
- std::array<Real, 1> aj = { a };
- std::array<Real, 1> bj = { b };
- return hypergeometric_pFq_checked_series_impl(aj, bj, z, pol, log_scale);
- }
- } } }
- #endif
|