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