435 lines · plain
1///////////////////////////////////////////////////////////////////////////////2// Copyright 2014 Anton Bikineev3// Copyright 2014 Christopher Kormanyos4// Copyright 2014 John Maddock5// Copyright 2014 Paul Bristow6// Distributed under the Boost7// Software License, Version 1.0. (See accompanying file8// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)9//10#ifndef BOOST_MATH_DETAIL_HYPERGEOMETRIC_SERIES_HPP11#define BOOST_MATH_DETAIL_HYPERGEOMETRIC_SERIES_HPP12 13#include <cmath>14#include <cstdint>15#include <boost/math/tools/series.hpp>16#include <boost/math/special_functions/gamma.hpp>17#include <boost/math/special_functions/trunc.hpp>18#include <boost/math/policies/error_handling.hpp>19 20 namespace boost { namespace math { namespace detail {21 22 // primary template for term of Taylor series23 template <class T, unsigned p, unsigned q>24 struct hypergeometric_pFq_generic_series_term;25 26 // partial specialization for 0F127 template <class T>28 struct hypergeometric_pFq_generic_series_term<T, 0u, 1u>29 {30 typedef T result_type;31 32 hypergeometric_pFq_generic_series_term(const T& b, const T& z)33 : n(0), term(1), b(b), z(z)34 {35 }36 37 T operator()()38 {39 BOOST_MATH_STD_USING40 const T r = term;41 term *= ((1 / ((b + n) * (n + 1))) * z);42 ++n;43 return r;44 }45 46 private:47 unsigned n;48 T term;49 const T b, z;50 };51 52 // partial specialization for 1F053 template <class T>54 struct hypergeometric_pFq_generic_series_term<T, 1u, 0u>55 {56 typedef T result_type;57 58 hypergeometric_pFq_generic_series_term(const T& a, const T& z)59 : n(0), term(1), a(a), z(z)60 {61 }62 63 T operator()()64 {65 BOOST_MATH_STD_USING66 const T r = term;67 term *= (((a + n) / (n + 1)) * z);68 ++n;69 return r;70 }71 72 private:73 unsigned n;74 T term;75 const T a, z;76 };77 78 // partial specialization for 1F179 template <class T>80 struct hypergeometric_pFq_generic_series_term<T, 1u, 1u>81 {82 typedef T result_type;83 84 hypergeometric_pFq_generic_series_term(const T& a, const T& b, const T& z)85 : n(0), term(1), a(a), b(b), z(z)86 {87 }88 89 T operator()()90 {91 BOOST_MATH_STD_USING92 const T r = term;93 term *= (((a + n) / ((b + n) * (n + 1))) * z);94 ++n;95 return r;96 }97 98 private:99 unsigned n;100 T term;101 const T a, b, z;102 };103 104 // partial specialization for 1F2105 template <class T>106 struct hypergeometric_pFq_generic_series_term<T, 1u, 2u>107 {108 typedef T result_type;109 110 hypergeometric_pFq_generic_series_term(const T& a, const T& b1, const T& b2, const T& z)111 : n(0), term(1), a(a), b1(b1), b2(b2), z(z)112 {113 }114 115 T operator()()116 {117 BOOST_MATH_STD_USING118 const T r = term;119 term *= (((a + n) / ((b1 + n) * (b2 + n) * (n + 1))) * z);120 ++n;121 return r;122 }123 124 private:125 unsigned n;126 T term;127 const T a, b1, b2, z;128 };129 130 // partial specialization for 2F0131 template <class T>132 struct hypergeometric_pFq_generic_series_term<T, 2u, 0u>133 {134 typedef T result_type;135 136 hypergeometric_pFq_generic_series_term(const T& a1, const T& a2, const T& z)137 : n(0), term(1), a1(a1), a2(a2), z(z)138 {139 }140 141 T operator()()142 {143 BOOST_MATH_STD_USING144 const T r = term;145 term *= (((a1 + n) * (a2 + n) / (n + 1)) * z);146 ++n;147 return r;148 }149 150 private:151 unsigned n;152 T term;153 const T a1, a2, z;154 };155 156 // partial specialization for 2F1157 template <class T>158 struct hypergeometric_pFq_generic_series_term<T, 2u, 1u>159 {160 typedef T result_type;161 162 hypergeometric_pFq_generic_series_term(const T& a1, const T& a2, const T& b, const T& z)163 : n(0), term(1), a1(a1), a2(a2), b(b), z(z)164 {165 }166 167 T operator()()168 {169 BOOST_MATH_STD_USING170 const T r = term;171 term *= (((a1 + n) * (a2 + n) / ((b + n) * (n + 1))) * z);172 ++n;173 return r;174 }175 176 private:177 unsigned n;178 T term;179 const T a1, a2, b, z;180 };181 182 // we don't need to define extra check and make a polinom from183 // series, when p(i) and q(i) are negative integers and p(i) >= q(i)184 // as described in functions.wolfram.alpha, because we always185 // stop summation when result (in this case numerator) is zero.186 template <class T, unsigned p, unsigned q, class Policy>187 inline T sum_pFq_series(detail::hypergeometric_pFq_generic_series_term<T, p, q>& term, const Policy& pol)188 {189 BOOST_MATH_STD_USING190 std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();191 192 const T result = boost::math::tools::sum_series(term, boost::math::policies::get_epsilon<T, Policy>(), max_iter);193 194 policies::check_series_iterations<T>("boost::math::hypergeometric_pFq_generic_series<%1%>(%1%,%1%,%1%)", max_iter, pol);195 return result;196 }197 198 template <class T, class Policy>199 inline T hypergeometric_0F1_generic_series(const T& b, const T& z, const Policy& pol)200 {201 detail::hypergeometric_pFq_generic_series_term<T, 0u, 1u> s(b, z);202 return detail::sum_pFq_series(s, pol);203 }204 205 template <class T, class Policy>206 inline T hypergeometric_1F0_generic_series(const T& a, const T& z, const Policy& pol)207 {208 detail::hypergeometric_pFq_generic_series_term<T, 1u, 0u> s(a, z);209 return detail::sum_pFq_series(s, pol);210 }211 212 template <class T, class Policy>213 inline T log_pochhammer(T z, unsigned n, const Policy pol, int* s = nullptr)214 {215 BOOST_MATH_STD_USING216#if 0217 if (z < 0)218 {219 if (n < -z)220 {221 if(s)222 *s = (n & 1 ? -1 : 1);223 return log_pochhammer(T(-z + (1 - (int)n)), n, pol);224 }225 else226 {227 int cross = itrunc(ceil(-z));228 return log_pochhammer(T(-z + (1 - cross)), cross, pol, s) + log_pochhammer(T(cross + z), n - cross, pol);229 }230 }231 else232#endif233 {234 if (z + n < 0)235 {236 T r = log_pochhammer(T(-z - n + 1), n, pol, s);237 if (s)238 *s *= (n & 1 ? -1 : 1);239 return r;240 }241 int s1, s2;242 auto r = static_cast<T>(boost::math::lgamma(T(z + n), &s1, pol) - boost::math::lgamma(z, &s2, pol));243 if(s)244 *s = s1 * s2;245 return r;246 }247 }248 249 template <class T, class Policy>250 inline T hypergeometric_1F1_generic_series(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling, const char* function)251 {252 BOOST_MATH_STD_USING253 T sum(0), term(1), upper_limit(sqrt(boost::math::tools::max_value<T>())), diff;254 T lower_limit(1 / upper_limit);255 unsigned n = 0;256 long long log_scaling_factor = lltrunc(boost::math::tools::log_max_value<T>()) - 2;257 T scaling_factor = exp(T(log_scaling_factor));258 T term_m1 = 0;259 long long local_scaling = 0;260 //261 // When a is very small, then (a+n)/n => 1 faster than262 // z / (b+n) => 1, as a result the series starts off263 // converging, then at some unspecified time very gradually264 // starts to diverge, potentially resulting in some very large265 // values being missed. As a result we need a check for small266 // a in the convergence criteria. Note that this issue occurs267 // even when all the terms are positive.268 //269 bool small_a = fabs(a) < 0.25;270 271 unsigned summit_location = 0;272 bool have_minima = false;273 T sq = 4 * a * z + b * b - 2 * b * z + z * z;274 if (sq >= 0)275 {276 T t = (-sqrt(sq) - b + z) / 2;277 if (t > 1) // Don't worry about a minima between 0 and 1.278 have_minima = true;279 t = (sqrt(sq) - b + z) / 2;280 if (t > 0)281 summit_location = itrunc(t);282 }283 284 if (summit_location > boost::math::policies::get_max_series_iterations<Policy>() / 4)285 {286 //287 // Skip forward to the location of the largest term in the series and288 // evaluate outwards from there:289 //290 int s1, s2;291 term = log_pochhammer(a, summit_location, pol, &s1) + summit_location * log(z) - log_pochhammer(b, summit_location, pol, &s2) - lgamma(T(summit_location + 1), pol);292 //std::cout << term << " " << log_pochhammer(boost::multiprecision::mpfr_float(a), summit_location, pol, &s1) + summit_location * log(boost::multiprecision::mpfr_float(z)) - log_pochhammer(boost::multiprecision::mpfr_float(b), summit_location, pol, &s2) - lgamma(boost::multiprecision::mpfr_float(summit_location + 1), pol) << std::endl;293 local_scaling = lltrunc(term);294 log_scaling += local_scaling;295 term = s1 * s2 * exp(term - local_scaling);296 //std::cout << term << " " << exp(log_pochhammer(boost::multiprecision::mpfr_float(a), summit_location, pol, &s1) + summit_location * log(boost::multiprecision::mpfr_float(z)) - log_pochhammer(boost::multiprecision::mpfr_float(b), summit_location, pol, &s2) - lgamma(boost::multiprecision::mpfr_float(summit_location + 1), pol) - local_scaling) << std::endl;297 n = summit_location;298 }299 else300 summit_location = 0;301 302 T saved_term = term;303 long long saved_scale = local_scaling;304 305 do306 {307 sum += term;308 //std::cout << n << " " << term * exp(boost::multiprecision::mpfr_float(local_scaling)) << " " << rising_factorial(boost::multiprecision::mpfr_float(a), n) * pow(boost::multiprecision::mpfr_float(z), n) / (rising_factorial(boost::multiprecision::mpfr_float(b), n) * factorial<boost::multiprecision::mpfr_float>(n)) << std::endl;309 if (fabs(sum) >= upper_limit)310 {311 sum /= scaling_factor;312 term /= scaling_factor;313 log_scaling += log_scaling_factor;314 local_scaling += log_scaling_factor;315 }316 if (fabs(sum) < lower_limit)317 {318 sum *= scaling_factor;319 term *= scaling_factor;320 log_scaling -= log_scaling_factor;321 local_scaling -= log_scaling_factor;322 }323 term_m1 = term;324 term *= (((a + n) / ((b + n) * (n + 1))) * z);325 if (n - summit_location > boost::math::policies::get_max_series_iterations<Policy>())326 return boost::math::policies::raise_evaluation_error(function, "Series did not converge, best value is %1%", sum, pol);327 ++n;328 diff = fabs(term / sum);329 } while ((diff > boost::math::policies::get_epsilon<T, Policy>()) || (fabs(term_m1) < fabs(term)) || (small_a && n < 10));330 331 //332 // See if we need to go backwards as well:333 //334 if (summit_location)335 {336 //337 // Backup state:338 //339 term = saved_term * exp(T(local_scaling - saved_scale));340 n = summit_location;341 term *= (b + (n - 1)) * n / ((a + (n - 1)) * z);342 --n;343 344 do345 {346 sum += term;347 //std::cout << n << " " << term * exp(boost::multiprecision::mpfr_float(local_scaling)) << " " << rising_factorial(boost::multiprecision::mpfr_float(a), n) * pow(boost::multiprecision::mpfr_float(z), n) / (rising_factorial(boost::multiprecision::mpfr_float(b), n) * factorial<boost::multiprecision::mpfr_float>(n)) << std::endl;348 if (n == 0)349 break;350 if (fabs(sum) >= upper_limit)351 {352 sum /= scaling_factor;353 term /= scaling_factor;354 log_scaling += log_scaling_factor;355 local_scaling += log_scaling_factor;356 }357 if (fabs(sum) < lower_limit)358 {359 sum *= scaling_factor;360 term *= scaling_factor;361 log_scaling -= log_scaling_factor;362 local_scaling -= log_scaling_factor;363 }364 term_m1 = term;365 term *= (b + (n - 1)) * n / ((a + (n - 1)) * z);366 if (summit_location - n > boost::math::policies::get_max_series_iterations<Policy>())367 return boost::math::policies::raise_evaluation_error(function, "Series did not converge, best value is %1%", sum, pol);368 --n;369 diff = fabs(term / sum);370 } while ((diff > boost::math::policies::get_epsilon<T, Policy>()) || (fabs(term_m1) < fabs(term)));371 }372 373 if (have_minima && n && summit_location)374 {375 //376 // There are a few terms starting at n == 0 which377 // haven't been accounted for yet...378 //379 unsigned backstop = n;380 n = 0;381 term = exp(T(-local_scaling));382 do383 {384 sum += term;385 //std::cout << n << " " << term << " " << sum << std::endl;386 if (fabs(sum) >= upper_limit)387 {388 sum /= scaling_factor;389 term /= scaling_factor;390 log_scaling += log_scaling_factor;391 }392 if (fabs(sum) < lower_limit)393 {394 sum *= scaling_factor;395 term *= scaling_factor;396 log_scaling -= log_scaling_factor;397 }398 //term_m1 = term;399 term *= (((a + n) / ((b + n) * (n + 1))) * z);400 if (n > boost::math::policies::get_max_series_iterations<Policy>())401 return boost::math::policies::raise_evaluation_error(function, "Series did not converge, best value is %1%", sum, pol);402 if (++n == backstop)403 break; // we've caught up with ourselves.404 diff = fabs(term / sum);405 } while ((diff > boost::math::policies::get_epsilon<T, Policy>())/* || (fabs(term_m1) < fabs(term))*/);406 }407 //std::cout << sum << std::endl;408 return sum;409 }410 411 template <class T, class Policy>412 inline T hypergeometric_1F2_generic_series(const T& a, const T& b1, const T& b2, const T& z, const Policy& pol)413 {414 detail::hypergeometric_pFq_generic_series_term<T, 1u, 2u> s(a, b1, b2, z);415 return detail::sum_pFq_series(s, pol);416 }417 418 template <class T, class Policy>419 inline T hypergeometric_2F0_generic_series(const T& a1, const T& a2, const T& z, const Policy& pol)420 {421 detail::hypergeometric_pFq_generic_series_term<T, 2u, 0u> s(a1, a2, z);422 return detail::sum_pFq_series(s, pol);423 }424 425 template <class T, class Policy>426 inline T hypergeometric_2F1_generic_series(const T& a1, const T& a2, const T& b, const T& z, const Policy& pol)427 {428 detail::hypergeometric_pFq_generic_series_term<T, 2u, 1u> s(a1, a2, b, z);429 return detail::sum_pFq_series(s, pol);430 }431 432 } } } // namespaces433 434#endif // BOOST_MATH_DETAIL_HYPERGEOMETRIC_SERIES_HPP435