823 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_HYPERGEOMETRIC_1F1_HPP11#define BOOST_MATH_HYPERGEOMETRIC_1F1_HPP12 13#include <boost/math/tools/config.hpp>14#include <boost/math/policies/policy.hpp>15#include <boost/math/policies/error_handling.hpp>16#include <boost/math/special_functions/detail/hypergeometric_series.hpp>17#include <boost/math/special_functions/detail/hypergeometric_asym.hpp>18#include <boost/math/special_functions/detail/hypergeometric_rational.hpp>19#include <boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp>20#include <boost/math/special_functions/detail/hypergeometric_1F1_by_ratios.hpp>21#include <boost/math/special_functions/detail/hypergeometric_pade.hpp>22#include <boost/math/special_functions/detail/hypergeometric_1F1_bessel.hpp>23#include <boost/math/special_functions/detail/hypergeometric_1F1_scaled_series.hpp>24#include <boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp>25#include <boost/math/special_functions/detail/hypergeometric_1F1_addition_theorems_on_z.hpp>26#include <boost/math/special_functions/detail/hypergeometric_1F1_large_abz.hpp>27#include <boost/math/special_functions/detail/hypergeometric_1F1_small_a_negative_b_by_ratio.hpp>28#include <boost/math/special_functions/detail/hypergeometric_1F1_negative_b_regions.hpp>29 30namespace boost { namespace math { namespace detail {31 32 // check when 1F1 series can't decay to polynom33 template <class T>34 inline bool check_hypergeometric_1F1_parameters(const T& a, const T& b)35 {36 BOOST_MATH_STD_USING37 38 if ((b <= 0) && (b == floor(b)))39 {40 if ((a >= 0) || (a < b) || (a != floor(a)))41 return false;42 }43 44 return true;45 }46 47 template <class T, class Policy>48 T hypergeometric_1F1_divergent_fallback(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling)49 {50 BOOST_MATH_STD_USING51 const char* function = "hypergeometric_1F1_divergent_fallback<%1%>(%1%,%1%,%1%)";52 //53 // We get here if either:54 // 1) We decide up front that Tricomi's method won't work, or:55 // 2) We've called Tricomi's method and it's failed.56 //57 if (b > 0)58 {59 // Commented out since recurrence seems to always be better?60#if 061 if ((z < b) && (a > -50))62 // Might as well use a recurrence in preference to z-recurrence:63 return hypergeometric_1F1_backward_recurrence_for_negative_a(a, b, z, pol, function, log_scaling);64 T z_limit = fabs((2 * a - b) / (sqrt(fabs(a))));65 int k = 1 + itrunc(z - z_limit);66 // If k is too large we destroy all the digits in the result:67 T convergence_at_50 = (b - a + 50) * k / (z * 50);68 if ((k > 0) && (k < 50) && (fabs(convergence_at_50) < 1) && (z > z_limit))69 {70 return boost::math::detail::hypergeometric_1f1_recurrence_on_z_minus_zero(a, b, T(z - k), k, pol, log_scaling);71 }72#endif73 if (z < b)74 return hypergeometric_1F1_backward_recurrence_for_negative_a(a, b, z, pol, function, log_scaling);75 else76 return hypergeometric_1F1_backwards_recursion_on_b_for_negative_a(a, b, z, pol, function, log_scaling);77 }78 else // b < 079 {80 if (a < 0)81 {82 if ((b < a) && (z < -b / 4))83 // Defensive programming: it is *almost* certain that we can never get here, proving that is hard though...84 return hypergeometric_1F1_from_function_ratio_negative_ab(a, b, z, pol, log_scaling); // LCOV_EXCL_LINE85 else86 {87 //88 // Solve (a+n)z/((b+n)n) == 1 for n, the number of iterations till the series starts to converge.89 // If this is well away from the origin then it's probably better to use the series to evaluate this.90 // Note that if sqr is negative then we have no solution, so assign an arbitrarily large value to the91 // number of iterations.92 //93 bool can_use_recursion = (z - b + 100 < boost::math::policies::get_max_series_iterations<Policy>()) && (100 - a < boost::math::policies::get_max_series_iterations<Policy>());94 T sqr = 4 * a * z + b * b - 2 * b * z + z * z;95 T iterations_to_convergence = sqr > 0 ? T(0.5f * (-sqrt(sqr) - b + z)) : T(-a - b);96 if(can_use_recursion && ((std::max)(a, b) + iterations_to_convergence > -300))97 return hypergeometric_1F1_backwards_recursion_on_b_for_negative_a(a, b, z, pol, function, log_scaling);98 //99 // When a < b and if we fall through to the series, then we get divergent behaviour when b crosses the origin100 // so ideally we would pick another method. Otherwise the terms immediately after b crosses the origin may101 // suffer catastrophic cancellation....102 //103 if((a < b) && can_use_recursion)104 return hypergeometric_1F1_backwards_recursion_on_b_for_negative_a(a, b, z, pol, function, log_scaling);105 }106 }107 else108 {109 //110 // Start by getting the domain of the recurrence relations, we get either:111 // -1 Backwards recursion is stable and the CF will converge to double precision.112 // +1 Forwards recursion is stable and the CF will converge to double precision.113 // 0 No man's land, we're not far enough away from the crossover point to get double precision from either CF.114 //115 // At higher than double precision we need to be further away from the crossover location to116 // get full converge, but it's not clear how much further - indeed at quad precision it's117 // basically impossible to ever get forwards iteration to work. Backwards seems to work118 // OK as long as a > 1 whatever the precision though.119 //120 int domain = hypergeometric_1F1_negative_b_recurrence_region(a, b, z);121 if ((domain < 0) && ((a > 1) || (boost::math::policies::digits<T, Policy>() <= 64)))122 return hypergeometric_1F1_from_function_ratio_negative_b(a, b, z, pol, log_scaling);123 else if (domain > 0)124 {125 if (boost::math::policies::digits<T, Policy>() <= 64)126 return hypergeometric_1F1_from_function_ratio_negative_b_forwards(a, b, z, pol, log_scaling);127 // LCOV_EXCL_START, what follows is multiprecision only128#ifndef BOOST_MATH_NO_EXCEPTIONS129 try130#endif131 {132 return hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling);133 }134#ifndef BOOST_MATH_NO_EXCEPTIONS135 catch (const evaluation_error&)136 {137 //138 // The series failed, try the recursions instead and hope we get at least double precision:139 //140 return hypergeometric_1F1_from_function_ratio_negative_b_forwards(a, b, z, pol, log_scaling);141 }142#endif143 // LCOV_EXCL_STOP144 }145 //146 // We could fall back to Tricomi's approximation if we're in the transition zone147 // between the above two regions. However, I've been unable to find any examples148 // where this is better than the series, and there are many cases where it leads to149 // quite grievous errors.150 /*151 else if (allow_tricomi)152 {153 T aa = a < 1 ? T(1) : a;154 if (z < fabs((2 * aa - b) / (sqrt(fabs(aa * b)))))155 return hypergeometric_1F1_AS_13_3_7_tricomi(a, b, z, pol, log_scaling);156 }157 */158 }159 }160 161 // If we get here, then we've run out of methods to try, use the checked series which will162 // raise an error if the result is garbage:163 return hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling);164 }165 166#if 0167 // Archived, not used, see comments at call site.168 template <class T>169 bool is_convergent_negative_z_series(const T& a, const T& b, const T& z, const T& b_minus_a)170 {171 BOOST_MATH_STD_USING172 //173 // Filter out some cases we don't want first:174 //175 if((b_minus_a > 0) && (b > 0))176 {177 if (a < 0)178 return false;179 }180 //181 // Generic check: we have small initial divergence and are convergent after 10 terms:182 //183 if ((fabs(z * a / b) < 2) && (fabs(z * (a + 10) / ((b + 10) * 10)) < 1))184 {185 // Double check for divergence when we cross the origin on a and b:186 if (a < 0)187 {188 T n = 3 - floor(a);189 if (fabs((a + n) * z / ((b + n) * n)) < 1)190 {191 if (b < 0)192 {193 T m = 3 - floor(b);194 if (fabs((a + m) * z / ((b + m) * m)) < 1)195 return true;196 }197 else198 return true;199 }200 }201 else if (b < 0)202 {203 T n = 3 - floor(b);204 if (fabs((a + n) * z / ((b + n) * n)) < 1)205 return true;206 }207 }208 if ((b > 0) && (a < 0))209 {210 //211 // For a and z both negative, we're OK with some initial divergence as long as212 // it occurs before we hit the origin, as to start with all the terms have the213 // same sign.214 //215 // https://www.wolframalpha.com/input/?i=solve+(a%2Bn)z+%2F+((b%2Bn)n)+%3D%3D+1+for+n216 //217 T sqr = 4 * a * z + b * b - 2 * b * z + z * z;218 T iterations_to_convergence = sqr > 0 ? T(0.5f * (-sqrt(sqr) - b + z)) : T(-a + b);219 if (iterations_to_convergence < 0)220 iterations_to_convergence = 0.5f * (sqrt(sqr) - b + z);221 if (a + iterations_to_convergence < -50)222 {223 // Need to check for divergence when we cross the origin on a:224 if (a > -1)225 return true;226 T n = 300 - floor(a);227 if(fabs((a + n) * z / ((b + n) * n)) < 1)228 return true;229 }230 }231 return false;232 }233#endif234 template <class T>235 inline T cyl_bessel_i_shrinkage_rate(const T& z)236 {237 // Approximately the ratio I_10.5(z/2) / I_9.5(z/2), this gives us an idea of how quickly238 // the Bessel terms in A&S 13.6.4 are converging:239 if (z < -160)240 return 1;241 if (z < -40)242 return 0.75f;243 if (z < -20)244 return 0.5f;245 if (z < -7)246 return 0.25f;247 if (z < -2)248 return 0.1f;249 return 0.05f;250 }251 252 template <class T>253 inline bool hypergeometric_1F1_is_13_3_6_region(const T& a, const T& b, const T& z)254 {255 BOOST_MATH_STD_USING256 if(fabs(a) == 0.5)257 return false;258 if ((z < 0) && (fabs(10 * a / b) < 1) && (fabs(a) < 50))259 {260 T shrinkage = cyl_bessel_i_shrinkage_rate(z);261 // We want the first term not too divergent, and convergence by term 10:262 if ((fabs((2 * a - 1) * (2 * a - b) / b) < 2) && (fabs(shrinkage * (2 * a + 9) * (2 * a - b + 10) / (10 * (b + 10))) < 0.75))263 return true;264 }265 return false;266 }267 268 template <class T>269 inline bool hypergeometric_1F1_need_kummer_reflection(const T& a, const T& b, const T& z)270 {271 BOOST_MATH_STD_USING272 //273 // Check to see if we should apply Kummer's relation or not:274 //275 if (z > 0)276 return false;277 if (z < -1)278 return true;279 //280 // When z is small and negative, things get more complex.281 // More often than not we do not need apply Kummer's relation and the282 // series is convergent as is, but we do need to check:283 //284 if (a > 0)285 {286 if (b > 0)287 {288 return fabs((a + 10) * z / (10 * (b + 10))) < 1; // Is the 10'th term convergent?289 }290 else291 {292 return true; // Likely to be divergent as b crosses the origin293 }294 }295 else // a < 0296 {297 if (b > 0)298 {299 return false; // Terms start off all positive and then by the time a crosses the origin we *must* be convergent.300 }301 else302 {303 return true; // Likely to be divergent as b crosses the origin, but hard to rationalise about!304 }305 }306 }307 308 309 template <class T, class Policy>310 T hypergeometric_1F1_imp(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling)311 {312 BOOST_MATH_STD_USING // exp, fabs, sqrt313 314 static const char* const function = "boost::math::hypergeometric_1F1<%1%,%1%,%1%>(%1%,%1%,%1%)";315 316 if ((z == 0) || (a == 0))317 return T(1);318 319 // undefined result:320 if (!detail::check_hypergeometric_1F1_parameters(a, b))321 return policies::raise_domain_error<T>(function, "Function is indeterminate for negative integer b = %1%.", b, pol);322 323 // other checks:324 if (a == -1)325 {326 T r = 1 - (z / b);327 if (fabs(r) < 0.5)328 r = (b - z) / b;329 return r;330 }331 332 const T b_minus_a = b - a;333 334 // 0f0 a == b case;335 if (b_minus_a == 0)336 {337 if ((a < 0) && (floor(a) == a))338 {339 // Special case, use the truncated series to match what Mathematica does.340 if ((a < -20) && (z > 0) && (z < 1))341 {342 // https://functions.wolfram.com/HypergeometricFunctions/Hypergeometric1F1/03/01/04/02/0002/343 return exp(z) * boost::math::gamma_q(1 - a, z, pol);344 }345 // https://functions.wolfram.com/HypergeometricFunctions/Hypergeometric1F1/03/01/04/02/0003/346 return hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling);347 }348 long long scale = lltrunc(z, pol);349 log_scaling += scale;350 return exp(z - scale);351 }352 // Special case for b-a = -1, we don't use for small a as it throws the digits of a away and leads to large errors:353 if ((b_minus_a == -1) && (fabs(a) > 0.5))354 {355 // for negative small integer a it is reasonable to use truncated series - polynomial356 if ((a < 0) && (a == ceil(a)) && (a > -50))357 return detail::hypergeometric_1F1_generic_series(a, b, z, pol, log_scaling, function);358 359 log_scaling = lltrunc(floor(z));360 T local_z = z - log_scaling;361 return (b + z) * exp(local_z) / b;362 }363 364 if ((a == 1) && (b == 2))365 return boost::math::expm1(z, pol) / z;366 367 if ((b - a == b) && (fabs(z / b) < policies::get_epsilon<T, Policy>()))368 return 1;369 //370 // Special case for A&S 13.3.6:371 //372 if (z < 0)373 {374 if (hypergeometric_1F1_is_13_3_6_region(a, b, z))375 {376 // a is tiny compared to b, and z < 0377 // 13.3.6 appears to be the most efficient and often the most accurate method.378 T r = boost::math::detail::hypergeometric_1F1_AS_13_3_6(b_minus_a, b, T(-z), a, pol, log_scaling);379 long long scale = lltrunc(z, pol);380 log_scaling += scale;381 return r * exp(z - scale);382 }383 if ((b < 0) && (fabs(a) < 1e-2))384 {385 //386 // This is a tricky area, potentially we have no good method at all:387 //388 if (b - ceil(b) == a)389 {390 // Fractional parts of a and b are genuinely equal, we might as well391 // apply Kummer's relation and get a truncated series:392 long long scaling = lltrunc(z);393 T r = exp(z - scaling) * detail::hypergeometric_1F1_imp<T>(b_minus_a, b, -z, pol, log_scaling);394 log_scaling += scaling;395 return r;396 }397 if ((b < -1) && (max_b_for_1F1_small_a_negative_b_by_ratio(z) < b))398 return hypergeometric_1F1_small_a_negative_b_by_ratio(a, b, z, pol, log_scaling);399 if ((b > -1) && (b < -0.5f))400 {401 // Recursion is meta-stable:402 T first = hypergeometric_1F1_imp(a, T(b + 2), z, pol);403 T second = hypergeometric_1F1_imp(a, T(b + 1), z, pol);404 return tools::apply_recurrence_relation_backward(hypergeometric_1F1_recurrence_small_b_coefficients<T>(a, b, z, 1), 1, first, second);405 }406 //407 // We've got nothing left but 13.3.6, even though it may be initially divergent:408 //409 T r = boost::math::detail::hypergeometric_1F1_AS_13_3_6(b_minus_a, b, T(-z), a, pol, log_scaling);410 long long scale = lltrunc(z, pol);411 log_scaling += scale;412 return r * exp(z - scale);413 }414 }415 //416 // Asymptotic expansion for large z417 // TODO: check region for higher precision types.418 // Use recurrence relations to move to this region when a and b are also large.419 //420 if (detail::hypergeometric_1F1_asym_region(a, b, z, pol))421 {422 long long saved_scale = log_scaling;423#ifndef BOOST_MATH_NO_EXCEPTIONS424 try425#endif426 {427 return hypergeometric_1F1_asym_large_z_series(a, b, z, pol, log_scaling);428 }429#ifndef BOOST_MATH_NO_EXCEPTIONS430 catch (const evaluation_error&)431 {432 }433#endif434 //435 // Very occasionally our convergence criteria don't quite go to full precision436 // and we have to try another method:437 //438 log_scaling = saved_scale;439 }440 441 if ((fabs(a * z / b) < 3.5) && (fabs(z * 100) < fabs(b)) && ((fabs(a) > 1e-2) || (b < -5)))442 return detail::hypergeometric_1F1_rational(a, b, z, pol);443 444 if (hypergeometric_1F1_need_kummer_reflection(a, b, z))445 {446 if (a == 1)447 return detail::hypergeometric_1F1_pade(b, z, pol);448#if 0449 //450 // Commented out: is_convergent_negative_z_series is fine so far as it goes451 // but there appear to be no cases that use it, and in extremis, we will452 // fall through to the series evaluation anyway.453 //454 if (is_convergent_negative_z_series(a, b, z, b_minus_a))455 {456 if ((boost::math::sign(b_minus_a) == boost::math::sign(b)) && ((b > 0) || (b < -200)))457 {458 // Series is close enough to convergent that we should be OK,459 // In this domain b - a ~ b and since 1F1[a, a, z] = e^z 1F1[b-a, b, -z]460 // and 1F1[a, a, -z] = e^-z the result must necessarily be somewhere near unity.461 // We have to rule out b small and negative because if b crosses the origin early462 // in the series (before we're pretty much converged) then all bets are off.463 // Note that this can go badly wrong when b and z are both large and negative,464 // in that situation the series goes in waves of large and small values which465 // may or may not cancel out. Likewise the initial part of the series may or may466 // not converge, and even if it does may or may not give a correct answer!467 // For example 1F1[-small, -1252.5, -1043.7] can loose up to ~800 digits due to468 // cancellation and is basically incalculable via this method.469 return hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling);470 }471 }472#endif473 if ((b < 0) && (floor(b) == b))474 {475 // Negative integer b, so a must be a negative integer too.476 // Kummer's transformation fails here!477 if(a > -50)478 return detail::hypergeometric_1F1_generic_series(a, b, z, pol, log_scaling, function);479 // Is there anything better than this??480 return hypergeometric_1F1_imp(a, float_next(b), z, pol, log_scaling);481 }482 else483 {484 // Let's otherwise make z positive (almost always)485 // by Kummer's transformation486 // (we also don't transform if z belongs to [-1,0])487 // Also note that Kummer's transformation fails when b is 488 // a negative integer, although this seems to be unmentioned489 // in the literature...490 long long scaling = lltrunc(z);491 T r = exp(z - scaling) * detail::hypergeometric_1F1_imp<T>(b_minus_a, b, -z, pol, log_scaling);492 log_scaling += scaling;493 return r;494 }495 }496 //497 // Check for initial divergence:498 //499 bool series_is_divergent = (a + 1) * z / (b + 1) < -1;500 if (series_is_divergent && (a < 0) && (b < 0) && (a > -1))501 series_is_divergent = false; // Best off taking the series in this situation502 //503 // If series starts off non-divergent, and becomes divergent later504 // then it's because both a and b are negative, so check for later505 // divergence as well:506 //507 if (!series_is_divergent && (a < 0) && (b < 0) && (b > a))508 {509 //510 // We need to exclude situations where we're over the initial "hump"511 // in the series terms (ie series has already converged by the time512 // b crosses the origin:513 //514 //T fa = fabs(a);515 //T fb = fabs(b);516 T convergence_point = sqrt((a - 1) * (a - b)) - a;517 if (-b < convergence_point)518 {519 T n = -floor(b);520 series_is_divergent = (a + n) * z / ((b + n) * n) < -1;521 }522 }523 else if (!series_is_divergent && (b < 0) && (a > 0))524 {525 // Series almost always become divergent as b crosses the origin:526 series_is_divergent = true;527 }528 if (series_is_divergent && (b < -1) && (b > -5) && (a > b))529 series_is_divergent = false; // don't bother with divergence, series will be OK530 531 //532 // Test for alternating series due to negative a,533 // in particular, see if the series is initially divergent534 // If so use the recurrence relation on a:535 //536 if (series_is_divergent)537 {538 if((a < 0) && (floor(a) == a) && (-a < policies::get_max_series_iterations<Policy>()))539 // This works amazingly well for negative integer a:540 return hypergeometric_1F1_backward_recurrence_for_negative_a(a, b, z, pol, function, log_scaling);541 //542 // In what follows we have to set limits on how large z can be otherwise543 // the Bessel series become large and divergent and all the digits cancel out.544 // The criteria are distinctly empiracle rather than based on a firm analysis545 // of the terms in the series.546 //547 if (b > 0)548 {549 T z_limit = fabs((2 * a - b) / (sqrt(fabs(a))));550 if ((z < z_limit) && hypergeometric_1F1_is_tricomi_viable_positive_b(a, b, z))551 return detail::hypergeometric_1F1_AS_13_3_7_tricomi(a, b, z, pol, log_scaling);552 }553 else // b < 0554 {555 if (a < 0)556 {557 T z_limit = fabs((2 * a - b) / (sqrt(fabs(a))));558 //559 // I hate these hard limits, but they're about the best we can do to try and avoid560 // Bessel function internal failures: these will be caught and handled561 // but up the expense of this function call:562 //563 if (((z < z_limit) || (a > -500)) && ((b > -500) || (b - 2 * a > 0)) && (z < -a))564 {565 //566 // Outside this domain we will probably get better accuracy from the recursive methods.567 //568 if(!(((a < b) && (z > -b)) || (z > z_limit)))569 return detail::hypergeometric_1F1_AS_13_3_7_tricomi(a, b, z, pol, log_scaling);570 //571 // When b and z are both very small, we get large errors from the recurrence methods572 // in the fallbacks. Tricomi seems to work well here, as does direct series evaluation573 // at least some of the time. Picking the right method is not easy, and sometimes this574 // is much worse than the fallback. Overall though, it's a reasonable choice that keeps575 // the very worst errors under control.576 //577 if(b > -1)578 return detail::hypergeometric_1F1_AS_13_3_7_tricomi(a, b, z, pol, log_scaling);579 }580 }581 //582 // We previously used Tricomi here, but it appears to be worse than583 // the recurrence-based algorithms in hypergeometric_1F1_divergent_fallback.584 /*585 else586 {587 T aa = a < 1 ? T(1) : a;588 if (z < fabs((2 * aa - b) / (sqrt(fabs(aa * b)))))589 return detail::hypergeometric_1F1_AS_13_3_7_tricomi(a, b, z, pol, log_scaling);590 }*/591 }592 593 return hypergeometric_1F1_divergent_fallback(a, b, z, pol, log_scaling);594 }595 596 if (hypergeometric_1F1_is_13_3_6_region(b_minus_a, b, T(-z)))597 {598 // b_minus_a is tiny compared to b, and -z < 0599 // 13.3.6 appears to be the most efficient and often the most accurate method.600 return boost::math::detail::hypergeometric_1F1_AS_13_3_6(a, b, z, b_minus_a, pol, log_scaling);601 }602#if 0603 if ((a > 0) && (b > 0) && (a * z / b > 2))604 {605 //606 // Series is initially divergent and slow to converge, see if applying607 // Kummer's relation can improve things:608 //609 if (is_convergent_negative_z_series(b_minus_a, b, T(-z), b_minus_a))610 {611 long long scaling = lltrunc(z);612 T r = exp(z - scaling) * detail::hypergeometric_1F1_checked_series_impl(b_minus_a, b, T(-z), pol, log_scaling);613 log_scaling += scaling;614 return r;615 }616 617 }618#endif619 if ((a > 0) && (b > 0) && (a * z > 50))620 return detail::hypergeometric_1F1_large_abz(a, b, z, pol, log_scaling);621 622 if (b < 0)623 return detail::hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling);624 625 return detail::hypergeometric_1F1_generic_series(a, b, z, pol, log_scaling, function);626 }627 628 template <class T, class Policy>629 inline T hypergeometric_1F1_imp(const T& a, const T& b, const T& z, const Policy& pol)630 {631 BOOST_MATH_STD_USING // exp, fabs, sqrt632 long long log_scaling = 0;633 T result = hypergeometric_1F1_imp(a, b, z, pol, log_scaling);634 //635 // Actual result will be result * e^log_scaling.636 //637 static const thread_local long long max_scaling = lltrunc(boost::math::tools::log_max_value<T>()) - 2;638 static const thread_local T max_scale_factor = exp(T(max_scaling));639 640 while (log_scaling > max_scaling)641 {642 result *= max_scale_factor;643 log_scaling -= max_scaling;644 }645 while (log_scaling < -max_scaling)646 {647 result /= max_scale_factor;648 log_scaling += max_scaling;649 }650 if (log_scaling)651 result *= exp(T(log_scaling));652 return result;653 }654 655 template <class T, class Policy>656 inline T log_hypergeometric_1F1_imp(const T& a, const T& b, const T& z, int* sign, const Policy& pol)657 {658 BOOST_MATH_STD_USING // exp, fabs, sqrt659 long long log_scaling = 0;660 T result = hypergeometric_1F1_imp(a, b, z, pol, log_scaling);661 if (sign)662 *sign = result < 0 ? -1 : 1;663 result = log(fabs(result)) + log_scaling;664 return result;665 }666 667 template <class T, class Policy>668 inline T hypergeometric_1F1_regularized_imp(const T& a, const T& b, const T& z, const Policy& pol)669 {670 BOOST_MATH_STD_USING // exp, fabs, sqrt671 long long log_scaling = 0;672 T result = hypergeometric_1F1_imp(a, b, z, pol, log_scaling);673 //674 // Actual result will be result * e^log_scaling / tgamma(b).675 //676 int result_sign = 1;677 T scale = log_scaling - boost::math::lgamma(b, &result_sign, pol);678 679 static const thread_local T max_scaling = boost::math::tools::log_max_value<T>() - 2;680 static const thread_local T max_scale_factor = exp(max_scaling);681 682 while (scale > max_scaling)683 {684 if((fabs(result) > 1) && (fabs(tools::max_value<T>()) / result <= max_scale_factor))685 return policies::raise_overflow_error<T>("hypergeometric_1F1_regularized", nullptr, pol);686 // This is *probably* unreachable:687 // LCOV_EXCL_START688 result *= max_scale_factor;689 scale -= max_scaling;690 // LCOV_EXCL_STOP691 }692 while (scale < -max_scaling)693 {694 result /= max_scale_factor;695 scale += max_scaling;696 }697 if (scale != 0)698 {699 scale = exp(scale);700 if ((scale > 1) && (fabs(result) > 1) && (fabs(tools::max_value<T>() / result) <= scale))701 return policies::raise_overflow_error<T>("hypergeometric_1F1_regularized", nullptr, pol);702 result *= scale;703 }704 return result * result_sign;705 }706 707} // namespace detail708 709template <class T1, class T2, class T3, class Policy>710inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_1F1(T1 a, T2 b, T3 z, const Policy& /* pol */)711{712 BOOST_FPU_EXCEPTION_GUARD713 typedef typename tools::promote_args<T1, T2, T3>::type result_type;714 typedef typename policies::evaluation<result_type, Policy>::type value_type;715 typedef typename policies::normalise<716 Policy,717 policies::promote_float<false>,718 policies::promote_double<false>,719 policies::discrete_quantile<>,720 policies::assert_undefined<> >::type forwarding_policy;721 return policies::checked_narrowing_cast<result_type, Policy>(722 detail::hypergeometric_1F1_imp<value_type>(723 static_cast<value_type>(a),724 static_cast<value_type>(b),725 static_cast<value_type>(z),726 forwarding_policy()),727 "boost::math::hypergeometric_1F1<%1%>(%1%,%1%,%1%)");728}729 730template <class T1, class T2, class T3>731inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_1F1(T1 a, T2 b, T3 z)732{733 return hypergeometric_1F1(a, b, z, policies::policy<>());734}735 736template <class T1, class T2, class T3, class Policy>737inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_1F1_regularized(T1 a, T2 b, T3 z, const Policy& /* pol */)738{739 BOOST_FPU_EXCEPTION_GUARD740 typedef typename tools::promote_args<T1, T2, T3>::type result_type;741 typedef typename policies::evaluation<result_type, Policy>::type value_type;742 typedef typename policies::normalise<743 Policy,744 policies::promote_float<false>,745 policies::promote_double<false>,746 policies::discrete_quantile<>,747 policies::assert_undefined<> >::type forwarding_policy;748 return policies::checked_narrowing_cast<result_type, Policy>(749 detail::hypergeometric_1F1_regularized_imp<value_type>(750 static_cast<value_type>(a),751 static_cast<value_type>(b),752 static_cast<value_type>(z),753 forwarding_policy()),754 "boost::math::hypergeometric_1F1<%1%>(%1%,%1%,%1%)");755}756 757template <class T1, class T2, class T3>758inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_1F1_regularized(T1 a, T2 b, T3 z)759{760 return hypergeometric_1F1_regularized(a, b, z, policies::policy<>());761}762 763template <class T1, class T2, class T3, class Policy>764inline typename tools::promote_args<T1, T2, T3>::type log_hypergeometric_1F1(T1 a, T2 b, T3 z, const Policy& /* pol */)765{766 BOOST_FPU_EXCEPTION_GUARD767 typedef typename tools::promote_args<T1, T2, T3>::type result_type;768 typedef typename policies::evaluation<result_type, Policy>::type value_type;769 typedef typename policies::normalise<770 Policy,771 policies::promote_float<false>,772 policies::promote_double<false>,773 policies::discrete_quantile<>,774 policies::assert_undefined<> >::type forwarding_policy;775 return policies::checked_narrowing_cast<result_type, Policy>(776 detail::log_hypergeometric_1F1_imp<value_type>(777 static_cast<value_type>(a),778 static_cast<value_type>(b),779 static_cast<value_type>(z),780 0,781 forwarding_policy()),782 "boost::math::hypergeometric_1F1<%1%>(%1%,%1%,%1%)");783}784 785template <class T1, class T2, class T3>786inline typename tools::promote_args<T1, T2, T3>::type log_hypergeometric_1F1(T1 a, T2 b, T3 z)787{788 return log_hypergeometric_1F1(a, b, z, policies::policy<>());789}790 791template <class T1, class T2, class T3, class Policy>792inline typename tools::promote_args<T1, T2, T3>::type log_hypergeometric_1F1(T1 a, T2 b, T3 z, int* sign, const Policy& /* pol */)793{794 BOOST_FPU_EXCEPTION_GUARD795 typedef typename tools::promote_args<T1, T2, T3>::type result_type;796 typedef typename policies::evaluation<result_type, Policy>::type value_type;797 typedef typename policies::normalise<798 Policy,799 policies::promote_float<false>,800 policies::promote_double<false>,801 policies::discrete_quantile<>,802 policies::assert_undefined<> >::type forwarding_policy;803 return policies::checked_narrowing_cast<result_type, Policy>(804 detail::log_hypergeometric_1F1_imp<value_type>(805 static_cast<value_type>(a),806 static_cast<value_type>(b),807 static_cast<value_type>(z),808 sign,809 forwarding_policy()),810 "boost::math::hypergeometric_1F1<%1%>(%1%,%1%,%1%)");811}812 813template <class T1, class T2, class T3>814inline typename tools::promote_args<T1, T2, T3>::type log_hypergeometric_1F1(T1 a, T2 b, T3 z, int* sign)815{816 return log_hypergeometric_1F1(a, b, z, sign, policies::policy<>());817}818 819 820 } } // namespace boost::math821 822#endif // BOOST_MATH_HYPERGEOMETRIC_HPP823