1 ///////////////////////////////////////////////////////////////////////////////
2 // Copyright 2014 Anton Bikineev
3 // Copyright 2014 Christopher Kormanyos
4 // Copyright 2014 John Maddock
5 // Copyright 2014 Paul Bristow
6 // Distributed under the Boost
7 // Software License, Version 1.0. (See accompanying file
8 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
10 #ifndef BOOST_MATH_DETAIL_HYPERGEOMETRIC_SERIES_HPP
11 #define BOOST_MATH_DETAIL_HYPERGEOMETRIC_SERIES_HPP
13 #include <boost/math/tools/series.hpp>
14 #include <boost/math/special_functions/trunc.hpp>
15 #include <boost/math/policies/error_handling.hpp>
17 namespace boost { namespace math { namespace detail {
19 // primary template for term of Taylor series
20 template <class T, unsigned p, unsigned q>
21 struct hypergeometric_pFq_generic_series_term;
23 // partial specialization for 0F1
25 struct hypergeometric_pFq_generic_series_term<T, 0u, 1u>
27 typedef T result_type;
29 hypergeometric_pFq_generic_series_term(const T& b, const T& z)
30 : n(0), term(1), b(b), z(z)
38 term *= ((1 / ((b + n) * (n + 1))) * z);
49 // partial specialization for 1F0
51 struct hypergeometric_pFq_generic_series_term<T, 1u, 0u>
53 typedef T result_type;
55 hypergeometric_pFq_generic_series_term(const T& a, const T& z)
56 : n(0), term(1), a(a), z(z)
64 term *= (((a + n) / (n + 1)) * z);
75 // partial specialization for 1F1
77 struct hypergeometric_pFq_generic_series_term<T, 1u, 1u>
79 typedef T result_type;
81 hypergeometric_pFq_generic_series_term(const T& a, const T& b, const T& z)
82 : n(0), term(1), a(a), b(b), z(z)
90 term *= (((a + n) / ((b + n) * (n + 1))) * z);
101 // partial specialization for 1F2
103 struct hypergeometric_pFq_generic_series_term<T, 1u, 2u>
105 typedef T result_type;
107 hypergeometric_pFq_generic_series_term(const T& a, const T& b1, const T& b2, const T& z)
108 : n(0), term(1), a(a), b1(b1), b2(b2), z(z)
116 term *= (((a + n) / ((b1 + n) * (b2 + n) * (n + 1))) * z);
124 const T a, b1, b2, z;
127 // partial specialization for 2F0
129 struct hypergeometric_pFq_generic_series_term<T, 2u, 0u>
131 typedef T result_type;
133 hypergeometric_pFq_generic_series_term(const T& a1, const T& a2, const T& z)
134 : n(0), term(1), a1(a1), a2(a2), z(z)
142 term *= (((a1 + n) * (a2 + n) / (n + 1)) * z);
153 // partial specialization for 2F1
155 struct hypergeometric_pFq_generic_series_term<T, 2u, 1u>
157 typedef T result_type;
159 hypergeometric_pFq_generic_series_term(const T& a1, const T& a2, const T& b, const T& z)
160 : n(0), term(1), a1(a1), a2(a2), b(b), z(z)
168 term *= (((a1 + n) * (a2 + n) / ((b + n) * (n + 1))) * z);
176 const T a1, a2, b, z;
179 // we don't need to define extra check and make a polinom from
180 // series, when p(i) and q(i) are negative integers and p(i) >= q(i)
181 // as described in functions.wolfram.alpha, because we always
182 // stop summation when result (in this case numerator) is zero.
183 template <class T, unsigned p, unsigned q, class Policy>
184 inline T sum_pFq_series(detail::hypergeometric_pFq_generic_series_term<T, p, q>& term, const Policy& pol)
187 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
188 #if BOOST_WORKAROUND(BOOST_BORLANDC, BOOST_TESTED_AT(0x582))
190 const T result = boost::math::tools::sum_series(term, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
192 const T result = boost::math::tools::sum_series(term, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
194 policies::check_series_iterations<T>("boost::math::hypergeometric_pFq_generic_series<%1%>(%1%,%1%,%1%)", max_iter, pol);
198 template <class T, class Policy>
199 inline T hypergeometric_0F1_generic_series(const T& b, const T& z, const Policy& pol)
201 detail::hypergeometric_pFq_generic_series_term<T, 0u, 1u> s(b, z);
202 return detail::sum_pFq_series(s, pol);
205 template <class T, class Policy>
206 inline T hypergeometric_1F0_generic_series(const T& a, const T& z, const Policy& pol)
208 detail::hypergeometric_pFq_generic_series_term<T, 1u, 0u> s(a, z);
209 return detail::sum_pFq_series(s, pol);
212 template <class T, class Policy>
213 inline T log_pochhammer(T z, unsigned n, const Policy pol, int* s = 0)
222 *s = (n & 1 ? -1 : 1);
223 return log_pochhammer(T(-z + (1 - (int)n)), n, pol);
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);
236 T r = log_pochhammer(T(-z - n + 1), n, pol, s);
238 *s *= (n & 1 ? -1 : 1);
242 T r = boost::math::lgamma(T(z + n), &s1, pol) - boost::math::lgamma(z, &s2, pol);
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, int& log_scaling, const char* function)
253 T sum(0), term(1), upper_limit(sqrt(boost::math::tools::max_value<T>())), diff;
254 T lower_limit(1 / upper_limit);
256 int log_scaling_factor = itrunc(boost::math::tools::log_max_value<T>()) - 2;
257 T scaling_factor = exp(T(log_scaling_factor));
259 int local_scaling = 0;
261 // When a is very small, then (a+n)/n => 1 faster than
262 // z / (b+n) => 1, as a result the series starts off
263 // converging, then at some unspecified time very gradually
264 // starts to diverge, potentially resulting in some very large
265 // values being missed. As a result we need a check for small
266 // a in the convergence criteria. Note that this issue occurs
267 // even when all the terms are positive.
269 bool small_a = fabs(a) < 0.25;
271 unsigned summit_location = 0;
272 bool have_minima = false;
273 T sq = 4 * a * z + b * b - 2 * b * z + z * z;
276 T t = (-sqrt(sq) - b + z) / 2;
277 if (t > 1) // Don't worry about a minima between 0 and 1.
279 t = (sqrt(sq) - b + z) / 2;
281 summit_location = itrunc(t);
284 if (summit_location > boost::math::policies::get_max_series_iterations<Policy>() / 4)
287 // Skip forward to the location of the largest term in the series and
288 // evaluate outwards from there:
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 = itrunc(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;
303 int saved_scale = local_scaling;
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)
311 sum /= scaling_factor;
312 term /= scaling_factor;
313 log_scaling += log_scaling_factor;
314 local_scaling += log_scaling_factor;
316 if (fabs(sum) < lower_limit)
318 sum *= scaling_factor;
319 term *= scaling_factor;
320 log_scaling -= log_scaling_factor;
321 local_scaling -= log_scaling_factor;
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);
328 diff = fabs(term / sum);
329 } while ((diff > boost::math::policies::get_epsilon<T, Policy>()) || (fabs(term_m1) < fabs(term)) || (small_a && n < 10));
332 // See if we need to go backwards as well:
339 term = saved_term * exp(T(local_scaling - saved_scale));
341 term *= (b + (n - 1)) * n / ((a + (n - 1)) * z);
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;
350 if (fabs(sum) >= upper_limit)
352 sum /= scaling_factor;
353 term /= scaling_factor;
354 log_scaling += log_scaling_factor;
355 local_scaling += log_scaling_factor;
357 if (fabs(sum) < lower_limit)
359 sum *= scaling_factor;
360 term *= scaling_factor;
361 log_scaling -= log_scaling_factor;
362 local_scaling -= log_scaling_factor;
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);
369 diff = fabs(term / sum);
370 } while ((diff > boost::math::policies::get_epsilon<T, Policy>()) || (fabs(term_m1) < fabs(term)));
373 if (have_minima && n && summit_location)
376 // There are a few terms starting at n == 0 which
377 // haven't been accounted for yet...
379 unsigned backstop = n;
381 term = exp(T(-local_scaling));
385 //std::cout << n << " " << term << " " << sum << std::endl;
386 if (fabs(sum) >= upper_limit)
388 sum /= scaling_factor;
389 term /= scaling_factor;
390 log_scaling += log_scaling_factor;
392 if (fabs(sum) < lower_limit)
394 sum *= scaling_factor;
395 term *= scaling_factor;
396 log_scaling -= log_scaling_factor;
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);
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))*/);
407 //std::cout << sum << std::endl;
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)
414 detail::hypergeometric_pFq_generic_series_term<T, 1u, 2u> s(a, b1, b2, z);
415 return detail::sum_pFq_series(s, pol);
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)
421 detail::hypergeometric_pFq_generic_series_term<T, 2u, 0u> s(a1, a2, z);
422 return detail::sum_pFq_series(s, pol);
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)
428 detail::hypergeometric_pFq_generic_series_term<T, 2u, 1u> s(a1, a2, b, z);
429 return detail::sum_pFq_series(s, pol);
434 #endif // BOOST_MATH_DETAIL_HYPERGEOMETRIC_SERIES_HPP