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) 9 // 10 #ifndef BOOST_MATH_DETAIL_HYPERGEOMETRIC_SERIES_HPP 11 #define BOOST_MATH_DETAIL_HYPERGEOMETRIC_SERIES_HPP 12 13 #include <boost/math/tools/series.hpp> 14 #include <boost/math/special_functions/trunc.hpp> 15 #include <boost/math/policies/error_handling.hpp> 16 17 namespace boost { namespace math { namespace detail { 18 19 // primary template for term of Taylor series 20 template <class T, unsigned p, unsigned q> 21 struct hypergeometric_pFq_generic_series_term; 22 23 // partial specialization for 0F1 24 template <class T> 25 struct hypergeometric_pFq_generic_series_term<T, 0u, 1u> 26 { 27 typedef T result_type; 28 hypergeometric_pFq_generic_series_termboost::math::detail::hypergeometric_pFq_generic_series_term29 hypergeometric_pFq_generic_series_term(const T& b, const T& z) 30 : n(0), term(1), b(b), z(z) 31 { 32 } 33 operator ()boost::math::detail::hypergeometric_pFq_generic_series_term34 T operator()() 35 { 36 BOOST_MATH_STD_USING 37 const T r = term; 38 term *= ((1 / ((b + n) * (n + 1))) * z); 39 ++n; 40 return r; 41 } 42 43 private: 44 unsigned n; 45 T term; 46 const T b, z; 47 }; 48 49 // partial specialization for 1F0 50 template <class T> 51 struct hypergeometric_pFq_generic_series_term<T, 1u, 0u> 52 { 53 typedef T result_type; 54 hypergeometric_pFq_generic_series_termboost::math::detail::hypergeometric_pFq_generic_series_term55 hypergeometric_pFq_generic_series_term(const T& a, const T& z) 56 : n(0), term(1), a(a), z(z) 57 { 58 } 59 operator ()boost::math::detail::hypergeometric_pFq_generic_series_term60 T operator()() 61 { 62 BOOST_MATH_STD_USING 63 const T r = term; 64 term *= (((a + n) / (n + 1)) * z); 65 ++n; 66 return r; 67 } 68 69 private: 70 unsigned n; 71 T term; 72 const T a, z; 73 }; 74 75 // partial specialization for 1F1 76 template <class T> 77 struct hypergeometric_pFq_generic_series_term<T, 1u, 1u> 78 { 79 typedef T result_type; 80 hypergeometric_pFq_generic_series_termboost::math::detail::hypergeometric_pFq_generic_series_term81 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) 83 { 84 } 85 operator ()boost::math::detail::hypergeometric_pFq_generic_series_term86 T operator()() 87 { 88 BOOST_MATH_STD_USING 89 const T r = term; 90 term *= (((a + n) / ((b + n) * (n + 1))) * z); 91 ++n; 92 return r; 93 } 94 95 private: 96 unsigned n; 97 T term; 98 const T a, b, z; 99 }; 100 101 // partial specialization for 1F2 102 template <class T> 103 struct hypergeometric_pFq_generic_series_term<T, 1u, 2u> 104 { 105 typedef T result_type; 106 hypergeometric_pFq_generic_series_termboost::math::detail::hypergeometric_pFq_generic_series_term107 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) 109 { 110 } 111 operator ()boost::math::detail::hypergeometric_pFq_generic_series_term112 T operator()() 113 { 114 BOOST_MATH_STD_USING 115 const T r = term; 116 term *= (((a + n) / ((b1 + n) * (b2 + n) * (n + 1))) * z); 117 ++n; 118 return r; 119 } 120 121 private: 122 unsigned n; 123 T term; 124 const T a, b1, b2, z; 125 }; 126 127 // partial specialization for 2F0 128 template <class T> 129 struct hypergeometric_pFq_generic_series_term<T, 2u, 0u> 130 { 131 typedef T result_type; 132 hypergeometric_pFq_generic_series_termboost::math::detail::hypergeometric_pFq_generic_series_term133 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) 135 { 136 } 137 operator ()boost::math::detail::hypergeometric_pFq_generic_series_term138 T operator()() 139 { 140 BOOST_MATH_STD_USING 141 const T r = term; 142 term *= (((a1 + n) * (a2 + n) / (n + 1)) * z); 143 ++n; 144 return r; 145 } 146 147 private: 148 unsigned n; 149 T term; 150 const T a1, a2, z; 151 }; 152 153 // partial specialization for 2F1 154 template <class T> 155 struct hypergeometric_pFq_generic_series_term<T, 2u, 1u> 156 { 157 typedef T result_type; 158 hypergeometric_pFq_generic_series_termboost::math::detail::hypergeometric_pFq_generic_series_term159 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) 161 { 162 } 163 operator ()boost::math::detail::hypergeometric_pFq_generic_series_term164 T operator()() 165 { 166 BOOST_MATH_STD_USING 167 const T r = term; 168 term *= (((a1 + n) * (a2 + n) / ((b + n) * (n + 1))) * z); 169 ++n; 170 return r; 171 } 172 173 private: 174 unsigned n; 175 T term; 176 const T a1, a2, b, z; 177 }; 178 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> sum_pFq_series(detail::hypergeometric_pFq_generic_series_term<T,p,q> & term,const Policy & pol)184 inline T sum_pFq_series(detail::hypergeometric_pFq_generic_series_term<T, p, q>& term, const Policy& pol) 185 { 186 BOOST_MATH_STD_USING 187 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); 188 #if BOOST_WORKAROUND(BOOST_BORLANDC, BOOST_TESTED_AT(0x582)) 189 const T zero = 0; 190 const T result = boost::math::tools::sum_series(term, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero); 191 #else 192 const T result = boost::math::tools::sum_series(term, boost::math::policies::get_epsilon<T, Policy>(), max_iter); 193 #endif 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> hypergeometric_0F1_generic_series(const T & b,const T & z,const Policy & pol)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> hypergeometric_1F0_generic_series(const T & a,const T & z,const Policy & pol)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> log_pochhammer(T z,unsigned n,const Policy pol,int * s=0)213 inline T log_pochhammer(T z, unsigned n, const Policy pol, int* s = 0) 214 { 215 BOOST_MATH_STD_USING 216 #if 0 217 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 else 226 { 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 else 232 #endif 233 { 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 T r = 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> hypergeometric_1F1_generic_series(const T & a,const T & b,const T & z,const Policy & pol,int & log_scaling,const char * function)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) 251 { 252 BOOST_MATH_STD_USING 253 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 int log_scaling_factor = itrunc(boost::math::tools::log_max_value<T>()) - 2; 257 T scaling_factor = exp(T(log_scaling_factor)); 258 T term_m1 = 0; 259 int local_scaling = 0; 260 // 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. 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 and 288 // 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 = 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; 297 n = summit_location; 298 } 299 else 300 summit_location = 0; 301 302 T saved_term = term; 303 int saved_scale = local_scaling; 304 305 do 306 { 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 do 345 { 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 which 377 // haven't been accounted for yet... 378 // 379 unsigned backstop = n; 380 n = 0; 381 term = exp(T(-local_scaling)); 382 do 383 { 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> hypergeometric_1F2_generic_series(const T & a,const T & b1,const T & b2,const T & z,const Policy & pol)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> hypergeometric_2F0_generic_series(const T & a1,const T & a2,const T & z,const Policy & pol)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> hypergeometric_2F1_generic_series(const T & a1,const T & a2,const T & b,const T & z,const Policy & pol)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 } } } // namespaces 433 434 #endif // BOOST_MATH_DETAIL_HYPERGEOMETRIC_SERIES_HPP 435