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_HYPERGEOMETRIC_1F1_LARGE_ABZ_HPP_ 9 #define BOOST_HYPERGEOMETRIC_1F1_LARGE_ABZ_HPP_ 10 11 #include <boost/math/special_functions/detail/hypergeometric_1F1_bessel.hpp> 12 #include <boost/math/special_functions/detail/hypergeometric_series.hpp> 13 #include <boost/math/special_functions/gamma.hpp> 14 15 16 namespace boost { namespace math { namespace detail { 17 18 template <class T> is_negative_integer(const T & x)19 inline bool is_negative_integer(const T& x) 20 { 21 using std::floor; 22 return (x <= 0) && (floor(x) == x); 23 } 24 25 26 template <class T, class Policy> 27 struct hypergeometric_1F1_igamma_series 28 { 29 enum{ cache_size = 64 }; 30 31 typedef T result_type; hypergeometric_1F1_igamma_seriesboost::math::detail::hypergeometric_1F1_igamma_series32 hypergeometric_1F1_igamma_series(const T& alpha, const T& delta, const T& x, const Policy& pol) 33 : delta_poch(-delta), alpha_poch(alpha), x(x), k(0), cache_offset(0), pol(pol) 34 { 35 BOOST_MATH_STD_USING 36 T log_term = log(x) * -alpha; 37 log_scaling = itrunc(log_term - 3 - boost::math::tools::log_min_value<T>() / 50); 38 term = exp(log_term - log_scaling); 39 refill_cache(); 40 } operator ()boost::math::detail::hypergeometric_1F1_igamma_series41 T operator()() 42 { 43 if (k - cache_offset >= cache_size) 44 { 45 cache_offset += cache_size; 46 refill_cache(); 47 } 48 T result = term * gamma_cache[k - cache_offset]; 49 term *= delta_poch * alpha_poch / (++k * x); 50 delta_poch += 1; 51 alpha_poch += 1; 52 return result; 53 } refill_cacheboost::math::detail::hypergeometric_1F1_igamma_series54 void refill_cache() 55 { 56 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; 57 58 gamma_cache[cache_size - 1] = boost::math::gamma_p(alpha_poch + ((int)cache_size - 1), x, pol); 59 for (int i = cache_size - 1; i > 0; --i) 60 { 61 gamma_cache[i - 1] = gamma_cache[i] >= 1 ? T(1) : T(gamma_cache[i] + regularised_gamma_prefix(T(alpha_poch + (i - 1)), x, pol, lanczos_type()) / (alpha_poch + (i - 1))); 62 } 63 } 64 T delta_poch, alpha_poch, x, term; 65 T gamma_cache[cache_size]; 66 int k; 67 int log_scaling; 68 int cache_offset; 69 Policy pol; 70 }; 71 72 template <class T, class Policy> 73 T hypergeometric_1F1_igamma(const T& a, const T& b, const T& x, const T& b_minus_a, const Policy& pol, int& log_scaling) 74 { 75 BOOST_MATH_STD_USING 76 if (b_minus_a == 0) 77 { 78 // special case: M(a,a,z) == exp(z) 79 int scale = itrunc(x, pol); 80 log_scaling += scale; 81 return exp(x - scale); 82 } 83 hypergeometric_1F1_igamma_series<T, Policy> s(b_minus_a, a - 1, x, pol); 84 log_scaling += s.log_scaling; 85 boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); 86 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); 87 boost::math::policies::check_series_iterations<T>("boost::math::tgamma<%1%>(%1%,%1%)", max_iter, pol); 88 T log_prefix = x + boost::math::lgamma(b, pol) - boost::math::lgamma(a, pol); 89 int scale = itrunc(log_prefix); 90 log_scaling += scale; 91 return result * exp(log_prefix - scale); 92 } 93 94 template <class T, class Policy> 95 T hypergeometric_1F1_shift_on_a(T h, const T& a_local, const T& b_local, const T& x, int a_shift, const Policy& pol, int& log_scaling) 96 { 97 BOOST_MATH_STD_USING 98 T a = a_local + a_shift; 99 if (a_shift == 0) 100 return h; 101 else if (a_shift > 0) 102 { 103 // 104 // Forward recursion on a is stable as long as 2a-b+z > 0. 105 // If 2a-b+z < 0 then backwards recursion is stable even though 106 // the function may be strictly increasing with a. Potentially 107 // we may need to split the recurrnce in 2 sections - one using 108 // forward recursion, and one backwards. 109 // 110 // We will get the next seed value from the ratio 111 // on b as that's stable and quick to compute. 112 // 113 114 T crossover_a = (b_local - x) / 2; 115 int crossover_shift = itrunc(crossover_a - a_local); 116 117 if (crossover_shift > 1) 118 { 119 // 120 // Forwards recursion will start off unstable, but may switch to the stable direction later. 121 // Start in the middle and go in both directions: 122 // 123 if (crossover_shift > a_shift) 124 crossover_shift = a_shift; 125 crossover_a = a_local + crossover_shift; 126 boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(crossover_a, b_local, x); 127 boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); 128 T b_ratio = boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter); 129 boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol); 130 // 131 // Convert to a ratio: 132 // (1+a-b)M(a, b, z) - aM(a+1, b, z) + (b-1)M(a, b-1, z) = 0 133 // 134 // hence: M(a+1,b,z) = ((1+a-b) / a) M(a,b,z) + ((b-1) / a) M(a,b,z)/b_ratio 135 // 136 T first = 1; 137 T second = ((1 + crossover_a - b_local) / crossover_a) + ((b_local - 1) / crossover_a) / b_ratio; 138 // 139 // Recurse down to a_local, compare values and re-nomralise first and second: 140 // 141 boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> a_coef(crossover_a, b_local, x); 142 int backwards_scale = 0; 143 T comparitor = boost::math::tools::apply_recurrence_relation_backward(a_coef, crossover_shift, second, first, &backwards_scale); 144 log_scaling -= backwards_scale; 145 if ((h < 1) && (tools::max_value<T>() * h > comparitor)) 146 { 147 // Need to rescale! 148 int scale = itrunc(log(h), pol) + 1; 149 h *= exp(T(-scale)); 150 log_scaling += scale; 151 } 152 comparitor /= h; 153 first /= comparitor; 154 second /= comparitor; 155 // 156 // Now we can recurse forwards for the rest of the range: 157 // 158 if (crossover_shift < a_shift) 159 { 160 boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> a_coef_2(crossover_a + 1, b_local, x); 161 h = boost::math::tools::apply_recurrence_relation_forward(a_coef_2, a_shift - crossover_shift - 1, first, second, &log_scaling); 162 } 163 else 164 h = first; 165 } 166 else 167 { 168 // 169 // Regular case where forwards iteration is stable right from the start: 170 // 171 boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(a_local, b_local, x); 172 boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); 173 T b_ratio = boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter); 174 boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol); 175 // 176 // Convert to a ratio: 177 // (1+a-b)M(a, b, z) - aM(a+1, b, z) + (b-1)M(a, b-1, z) = 0 178 // 179 // hence: M(a+1,b,z) = ((1+a-b) / a) M(a,b,z) + ((b-1) / a) M(a,b,z)/b_ratio 180 // 181 T second = ((1 + a_local - b_local) / a_local) * h + ((b_local - 1) / a_local) * h / b_ratio; 182 boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> a_coef(a_local + 1, b_local, x); 183 h = boost::math::tools::apply_recurrence_relation_forward(a_coef, --a_shift, h, second, &log_scaling); 184 } 185 } 186 else 187 { 188 // 189 // We've calculated h for a larger value of a than we want, and need to recurse down. 190 // However, only forward iteration is stable, so calculate the ratio, compare values, 191 // and normalise. Note that we calculate the ratio on b and convert to a since the 192 // direction is the minimal solution for N->+INF. 193 // 194 // IMPORTANT: this is only currently called for a > b and therefore forwards iteration 195 // is the only stable direction as we will only iterate down until a ~ b, but we 196 // will check this with an assert: 197 // 198 BOOST_ASSERT(2 * a - b_local + x > 0); 199 boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(a, b_local, x); 200 boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); 201 T b_ratio = boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter); 202 boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol); 203 // 204 // Convert to a ratio: 205 // (1+a-b)M(a, b, z) - aM(a+1, b, z) + (b-1)M(a, b-1, z) = 0 206 // 207 // hence: M(a+1,b,z) = (1+a-b) / a M(a,b,z) + (b-1) / a M(a,b,z)/ (M(a,b,z)/M(a,b-1,z)) 208 // 209 T first = 1; // arbitrary value; 210 T second = ((1 + a - b_local) / a) + ((b_local - 1) / a) * (1 / b_ratio); 211 212 if (a_shift == -1) 213 h = h / second; 214 else 215 { 216 boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> a_coef(a + 1, b_local, x); 217 T comparitor = boost::math::tools::apply_recurrence_relation_forward(a_coef, -(a_shift + 1), first, second); 218 if (boost::math::tools::min_value<T>() * comparitor > h) 219 { 220 // Ooops, need to rescale h: 221 int rescale = itrunc(log(fabs(h))); 222 T scale = exp(T(-rescale)); 223 h *= scale; 224 log_scaling += rescale; 225 } 226 h = h / comparitor; 227 } 228 } 229 return h; 230 } 231 232 template <class T, class Policy> 233 T hypergeometric_1F1_shift_on_b(T h, const T& a, const T& b_local, const T& x, int b_shift, const Policy& pol, int& log_scaling) 234 { 235 BOOST_MATH_STD_USING 236 237 T b = b_local + b_shift; 238 if (b_shift == 0) 239 return h; 240 else if (b_shift > 0) 241 { 242 // 243 // We get here for b_shift > 0 when b > z. We can't use forward recursion on b - it's unstable, 244 // so grab the ratio and work backwards to b - b_shift and normalise. 245 // 246 boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(a, b, x); 247 boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); 248 249 T first = 1; // arbitrary value; 250 T second = 1 / boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter); 251 boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol); 252 if (b_shift == 1) 253 h = h / second; 254 else 255 { 256 // 257 // Reset coefficients and recurse: 258 // 259 boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef_2(a, b - 1, x); 260 int local_scale = 0; 261 T comparitor = boost::math::tools::apply_recurrence_relation_backward(b_coef_2, --b_shift, first, second, &local_scale); 262 log_scaling -= local_scale; 263 if (boost::math::tools::min_value<T>() * comparitor > h) 264 { 265 // Ooops, need to rescale h: 266 int rescale = itrunc(log(fabs(h))); 267 T scale = exp(T(-rescale)); 268 h *= scale; 269 log_scaling += rescale; 270 } 271 h = h / comparitor; 272 } 273 } 274 else 275 { 276 T second; 277 if (a == b_local) 278 { 279 // recurrence is trivial for a == b and method of ratios fails as the c-term goes to zero: 280 second = -b_local * (1 - b_local - x) * h / (b_local * (b_local - 1)); 281 } 282 else 283 { 284 BOOST_ASSERT(!is_negative_integer(b - a)); 285 boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(a, b_local, x); 286 boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); 287 second = h / boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter); 288 boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol); 289 } 290 if (b_shift == -1) 291 h = second; 292 else 293 { 294 boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef_2(a, b_local - 1, x); 295 h = boost::math::tools::apply_recurrence_relation_backward(b_coef_2, -(++b_shift), h, second, &log_scaling); 296 } 297 } 298 return h; 299 } 300 301 302 template <class T, class Policy> 303 T hypergeometric_1F1_large_igamma(const T& a, const T& b, const T& x, const T& b_minus_a, const Policy& pol, int& log_scaling) 304 { 305 BOOST_MATH_STD_USING 306 // 307 // We need a < b < z in order to ensure there's at least a chance of convergence, 308 // we can use recurrence relations to shift forwards on a+b or just a to achieve this, 309 // for decent accuracy, try to keep 2b - 1 < a < 2b < z 310 // 311 int b_shift = b * 2 < x ? 0 : itrunc(b - x / 2); 312 int a_shift = a > b - b_shift ? -itrunc(b - b_shift - a - 1) : -itrunc(b - b_shift - a); 313 314 if (a_shift < 0) 315 { 316 // Might as well do all the shifting on b as scale a downwards: 317 b_shift -= a_shift; 318 a_shift = 0; 319 } 320 321 T a_local = a - a_shift; 322 T b_local = b - b_shift; 323 T b_minus_a_local = (a_shift == 0) && (b_shift == 0) ? b_minus_a : b_local - a_local; 324 int local_scaling = 0; 325 T h = hypergeometric_1F1_igamma(a_local, b_local, x, b_minus_a_local, pol, local_scaling); 326 log_scaling += local_scaling; 327 328 // 329 // Apply shifts on a and b as required: 330 // 331 h = hypergeometric_1F1_shift_on_a(h, a_local, b_local, x, a_shift, pol, log_scaling); 332 h = hypergeometric_1F1_shift_on_b(h, a, b_local, x, b_shift, pol, log_scaling); 333 334 return h; 335 } 336 337 template <class T, class Policy> 338 T hypergeometric_1F1_large_series(const T& a, const T& b, const T& z, const Policy& pol, int& log_scaling) 339 { 340 BOOST_MATH_STD_USING 341 // 342 // We make a small, and b > z: 343 // 344 int a_shift(0), b_shift(0); 345 if (a * z > b) 346 { 347 a_shift = itrunc(a) - 5; 348 b_shift = b < z ? itrunc(b - z - 1) : 0; 349 } 350 // 351 // If a_shift is trivially small, there's really not much point in losing 352 // accuracy to save a couple of iterations: 353 // 354 if (a_shift < 5) 355 a_shift = 0; 356 T a_local = a - a_shift; 357 T b_local = b - b_shift; 358 T h = boost::math::detail::hypergeometric_1F1_generic_series(a_local, b_local, z, pol, log_scaling, "hypergeometric_1F1_large_series<%1%>(a,b,z)"); 359 // 360 // Apply shifts on a and b as required: 361 // 362 if (a_shift && (a_local == 0)) 363 { 364 // 365 // Shifting on a via method of ratios in hypergeometric_1F1_shift_on_a fails when 366 // a_local == 0. However, the value of h calculated was trivial (unity), so 367 // calculate a second 1F1 for a == 1 and recurse as normal: 368 // 369 int scale = 0; 370 T h2 = boost::math::detail::hypergeometric_1F1_generic_series(T(a_local + 1), b_local, z, pol, scale, "hypergeometric_1F1_large_series<%1%>(a,b,z)"); 371 if (scale != log_scaling) 372 { 373 h2 *= exp(T(scale - log_scaling)); 374 } 375 boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> coef(a_local + 1, b_local, z); 376 h = boost::math::tools::apply_recurrence_relation_forward(coef, a_shift - 1, h, h2, &log_scaling); 377 h = hypergeometric_1F1_shift_on_b(h, a, b_local, z, b_shift, pol, log_scaling); 378 } 379 else 380 { 381 h = hypergeometric_1F1_shift_on_a(h, a_local, b_local, z, a_shift, pol, log_scaling); 382 h = hypergeometric_1F1_shift_on_b(h, a, b_local, z, b_shift, pol, log_scaling); 383 } 384 return h; 385 } 386 387 template <class T, class Policy> 388 T hypergeometric_1F1_large_13_3_6_series(const T& a, const T& b, const T& z, const Policy& pol, int& log_scaling) 389 { 390 BOOST_MATH_STD_USING 391 // 392 // A&S 13.3.6 is good only when a ~ b, but isn't too fussy on the size of z. 393 // So shift b to match a (b shifting seems to be more stable via method of ratios). 394 // 395 int b_shift = itrunc(b - a); 396 T b_local = b - b_shift; 397 T h = boost::math::detail::hypergeometric_1F1_AS_13_3_6(a, b_local, z, T(b_local - a), pol, log_scaling); 398 return hypergeometric_1F1_shift_on_b(h, a, b_local, z, b_shift, pol, log_scaling); 399 } 400 401 template <class T, class Policy> 402 T hypergeometric_1F1_large_abz(const T& a, const T& b, const T& z, const Policy& pol, int& log_scaling) 403 { 404 BOOST_MATH_STD_USING 405 // 406 // This is the selection logic to pick the "best" method. 407 // We have a,b,z >> 0 and need to comute the approximate cost of each method 408 // and then select whichever wins out. 409 // 410 enum method 411 { 412 method_series = 0, 413 method_shifted_series, 414 method_gamma, 415 method_bessel 416 }; 417 // 418 // Cost of direct series, is the approx number of terms required until we hit convergence: 419 // 420 T current_cost = (sqrt(16 * z * (3 * a + z) + 9 * b * b - 24 * b * z) - 3 * b + 4 * z) / 6; 421 method current_method = method_series; 422 // 423 // Cost of shifted series, is the number of recurrences required to move to a zone where 424 // the series is convergent right from the start. 425 // Note that recurrence relations fail for very small b, and too many recurrences on a 426 // will completely destroy all our digits. 427 // Also note that the method fails when b-a is a negative integer unless b is already 428 // larger than z and thus does not need shifting. 429 // 430 T cost = a + ((b < z) ? T(z - b) : T(0)); 431 if((b > 1) && (cost < current_cost) && ((b > z) || !is_negative_integer(b-a))) 432 { 433 current_method = method_shifted_series; 434 current_cost = cost; 435 } 436 // 437 // Cost for gamma function method is the number of recurrences required to move it 438 // into a convergent zone, note that recurrence relations fail for very small b. 439 // Also add on a fudge factor to account for the fact that this method is both 440 // more expensive to compute (requires gamma functions), and less accurate than the 441 // methods above: 442 // 443 T b_shift = fabs(b * 2 < z ? T(0) : T(b - z / 2)); 444 T a_shift = fabs(a > b - b_shift ? T(-(b - b_shift - a - 1)) : T(-(b - b_shift - a))); 445 cost = 1000 + b_shift + a_shift; 446 if((b > 1) && (cost <= current_cost)) 447 { 448 current_method = method_gamma; 449 current_cost = cost; 450 } 451 // 452 // Cost for bessel approximation is the number of recurrences required to make a ~ b, 453 // Note that recurrence relations fail for very small b. We also have issue with large 454 // z: either overflow/numeric instability or else the series goes divergent. We seem to be 455 // OK for z smaller than log_max_value<Quad> though, maybe we can stretch a little further 456 // but that's not clear... 457 // Also need to add on a fudge factor to the cost to account for the fact that we need 458 // to calculate the Bessel functions, this is not quite as high as the gamma function 459 // method above as this is generally more accurate and so prefered if the methods are close: 460 // 461 cost = 50 + fabs(b - a); 462 if((b > 1) && (cost <= current_cost) && (z < tools::log_max_value<T>()) && (z < 11356) && (b - a != 0.5f)) 463 { 464 current_method = method_bessel; 465 current_cost = cost; 466 } 467 468 switch (current_method) 469 { 470 case method_series: 471 return detail::hypergeometric_1F1_generic_series(a, b, z, pol, log_scaling, "hypergeometric_1f1_large_abz<%1%>(a,b,z)"); 472 case method_shifted_series: 473 return detail::hypergeometric_1F1_large_series(a, b, z, pol, log_scaling); 474 case method_gamma: 475 return detail::hypergeometric_1F1_large_igamma(a, b, z, T(b - a), pol, log_scaling); 476 case method_bessel: 477 return detail::hypergeometric_1F1_large_13_3_6_series(a, b, z, pol, log_scaling); 478 } 479 return 0; // We don't get here! 480 } 481 482 } } } // namespaces 483 484 #endif // BOOST_HYPERGEOMETRIC_1F1_LARGE_ABZ_HPP_ 485