196 lines · plain
1 2///////////////////////////////////////////////////////////////////////////////3// Copyright 2018 John Maddock4// Distributed under the Boost5// Software License, Version 1.0. (See accompanying file6// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)7 8#ifndef BOOST_MATH_HYPERGEOMETRIC_PFQ_HPP9#define BOOST_MATH_HYPERGEOMETRIC_PFQ_HPP10 11#include <boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp>12#include <boost/math/tools/throw_exception.hpp>13#include <chrono>14#include <initializer_list>15 16namespace boost {17 namespace math {18 19 namespace detail {20 21 struct pFq_termination_exception : public std::runtime_error22 {23 pFq_termination_exception(const char* p) : std::runtime_error(p) {}24 };25 26 struct timed_iteration_terminator27 {28 timed_iteration_terminator(std::uintmax_t i, double t) : max_iter(i), max_time(t), start_time(std::chrono::system_clock::now()) {}29 30 bool operator()(std::uintmax_t iter)const31 {32 if (iter > max_iter)33 BOOST_MATH_THROW_EXCEPTION(boost::math::detail::pFq_termination_exception("pFq exceeded maximum permitted iterations."));34 if (std::chrono::duration<double>(std::chrono::system_clock::now() - start_time).count() > max_time)35 BOOST_MATH_THROW_EXCEPTION(boost::math::detail::pFq_termination_exception("pFq exceeded maximum permitted evaluation time."));36 return false;37 }38 39 std::uintmax_t max_iter;40 double max_time;41 std::chrono::system_clock::time_point start_time;42 };43 44 }45 46 template <class Seq, class Real, class Policy>47 inline typename tools::promote_args<Real, typename Seq::value_type>::type hypergeometric_pFq(const Seq& aj, const Seq& bj, const Real& z, Real* p_abs_error, const Policy& pol)48 {49 typedef typename tools::promote_args<Real, typename Seq::value_type>::type result_type;50 typedef typename policies::evaluation<result_type, Policy>::type value_type;51 typedef typename policies::normalise<52 Policy,53 policies::promote_float<false>,54 policies::promote_double<false>,55 policies::discrete_quantile<>,56 policies::assert_undefined<> >::type forwarding_policy;57 58 BOOST_MATH_STD_USING59 60 long long scale = 0;61 std::pair<value_type, value_type> r = boost::math::detail::hypergeometric_pFq_checked_series_impl(aj, bj, value_type(z), pol, boost::math::detail::iteration_terminator(boost::math::policies::get_max_series_iterations<forwarding_policy>()), scale);62 r.first *= exp(Real(scale));63 r.second *= exp(Real(scale));64 if (p_abs_error)65 *p_abs_error = static_cast<Real>(r.second) * boost::math::tools::epsilon<Real>();66 return policies::checked_narrowing_cast<result_type, Policy>(r.first, "boost::math::hypergeometric_pFq<%1%>(%1%,%1%,%1%)");67 }68 69 template <class Seq, class Real>70 inline typename tools::promote_args<Real, typename Seq::value_type>::type hypergeometric_pFq(const Seq& aj, const Seq& bj, const Real& z, Real* p_abs_error = 0)71 {72 return hypergeometric_pFq(aj, bj, z, p_abs_error, boost::math::policies::policy<>());73 }74 75 template <class R, class Real, class Policy>76 inline typename tools::promote_args<Real, R>::type hypergeometric_pFq(const std::initializer_list<R>& aj, const std::initializer_list<R>& bj, const Real& z, Real* p_abs_error, const Policy& pol)77 {78 return hypergeometric_pFq<std::initializer_list<R>, Real, Policy>(aj, bj, z, p_abs_error, pol);79 }80 81 template <class R, class Real>82 inline typename tools::promote_args<Real, R>::type hypergeometric_pFq(const std::initializer_list<R>& aj, const std::initializer_list<R>& bj, const Real& z, Real* p_abs_error = nullptr)83 {84 return hypergeometric_pFq<std::initializer_list<R>, Real>(aj, bj, z, p_abs_error);85 }86 87#ifndef BOOST_MATH_NO_EXCEPTIONS88 template <class T>89 struct scoped_precision90 {91 scoped_precision(unsigned p)92 {93 old_p = T::default_precision();94 T::default_precision(p);95 }96 ~scoped_precision()97 {98 T::default_precision(old_p);99 }100 unsigned old_p;101 };102 103 template <class Seq, class Real, class Policy>104 Real hypergeometric_pFq_precision(const Seq& aj, const Seq& bj, Real z, unsigned digits10, double timeout, const Policy& pol)105 {106 unsigned current_precision = digits10 + 5;107 108 for (auto ai = aj.begin(); ai != aj.end(); ++ai)109 {110 current_precision = (std::max)(current_precision, ai->precision());111 }112 for (auto bi = bj.begin(); bi != bj.end(); ++bi)113 {114 current_precision = (std::max)(current_precision, bi->precision());115 }116 current_precision = (std::max)(current_precision, z.precision());117 118 Real r, norm;119 std::vector<Real> aa(aj), bb(bj);120 do121 {122 scoped_precision<Real> p(current_precision);123 for (auto ai = aa.begin(); ai != aa.end(); ++ai)124 ai->precision(current_precision);125 for (auto bi = bb.begin(); bi != bb.end(); ++bi)126 bi->precision(current_precision);127 z.precision(current_precision);128 try129 {130 long long scale = 0;131 std::pair<Real, Real> rp = boost::math::detail::hypergeometric_pFq_checked_series_impl(aa, bb, z, pol, boost::math::detail::timed_iteration_terminator(boost::math::policies::get_max_series_iterations<Policy>(), timeout), scale);132 rp.first *= exp(Real(scale));133 rp.second *= exp(Real(scale));134 135 r = rp.first;136 norm = rp.second;137 138 unsigned cancellation;139 try {140 cancellation = itrunc(log10(abs(norm / r)));141 }142 catch (const boost::math::rounding_error&)143 {144 // Happens when r is near enough zero:145 cancellation = UINT_MAX;146 }147 if (cancellation >= current_precision - 1)148 {149 current_precision *= 2;150 continue;151 }152 unsigned precision_obtained = current_precision - 1 - cancellation;153 if (precision_obtained < digits10)154 {155 current_precision += digits10 - precision_obtained + 5;156 }157 else158 break;159 }160 catch (const boost::math::evaluation_error&)161 {162 current_precision *= 2;163 }164 catch (const detail::pFq_termination_exception& e)165 {166 //167 // Either we have exhausted the number of series iterations, or the timeout.168 // Either way we quit now.169 throw boost::math::evaluation_error(e.what());170 }171 } while (true);172 173 return r;174 }175 template <class Seq, class Real>176 Real hypergeometric_pFq_precision(const Seq& aj, const Seq& bj, const Real& z, unsigned digits10, double timeout = 0.5)177 {178 return hypergeometric_pFq_precision(aj, bj, z, digits10, timeout, boost::math::policies::policy<>());179 }180 181 template <class Real, class Policy>182 Real hypergeometric_pFq_precision(const std::initializer_list<Real>& aj, const std::initializer_list<Real>& bj, const Real& z, unsigned digits10, double timeout, const Policy& pol)183 {184 return hypergeometric_pFq_precision< std::initializer_list<Real>, Real>(aj, bj, z, digits10, timeout, pol);185 }186 template <class Real>187 Real hypergeometric_pFq_precision(const std::initializer_list<Real>& aj, const std::initializer_list<Real>& bj, const Real& z, unsigned digits10, double timeout = 0.5)188 {189 return hypergeometric_pFq_precision< std::initializer_list<Real>, Real>(aj, bj, z, digits10, timeout, boost::math::policies::policy<>());190 }191#endif192 }193} // namespaces194 195#endif // BOOST_MATH_BESSEL_ITERATORS_HPP196