180 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_ASYM_HPP11#define BOOST_MATH_HYPERGEOMETRIC_ASYM_HPP12 13#include <boost/math/special_functions/gamma.hpp>14#include <boost/math/special_functions/hypergeometric_2F0.hpp>15 16#ifdef _MSC_VER17#pragma warning(push)18#pragma warning(disable:4127)19#endif20 21 namespace boost { namespace math {22 23 namespace detail {24 25 //26 // Asymptotic series based on https://dlmf.nist.gov/13.7#E127 //28 // Note that a and b must not be negative integers, in addition29 // we require z > 0 and so apply Kummer's relation for z < 0.30 //31 template <class T, class Policy>32 inline T hypergeometric_1F1_asym_large_z_series(T a, const T& b, T z, const Policy& pol, long long& log_scaling)33 {34 BOOST_MATH_STD_USING35 static const char* function = "boost::math::hypergeometric_1F1_asym_large_z_series<%1%>(%1%, %1%, %1%)";36 T prefix;37 long long e;38 int s;39 if (z < 0)40 {41 a = b - a;42 z = -z;43 prefix = 1;44 }45 else46 {47 e = z > static_cast<T>((std::numeric_limits<long long>::max)()) ? (std::numeric_limits<long long>::max)() : lltrunc(z, pol);48 log_scaling += e;49 prefix = exp(z - e);50 }51 if ((fabs(a) < 10) && (fabs(b) < 10))52 {53 prefix *= pow(z, a) * pow(z, -b) * boost::math::tgamma(b, pol) / boost::math::tgamma(a, pol);54 }55 else56 {57 T t = log(z) * (a - b);58 e = lltrunc(t, pol);59 log_scaling += e;60 prefix *= exp(t - e);61 62 t = boost::math::lgamma(b, &s, pol);63 e = lltrunc(t, pol);64 log_scaling += e;65 prefix *= s * exp(t - e);66 67 t = boost::math::lgamma(a, &s, pol);68 e = lltrunc(t, pol);69 log_scaling -= e;70 prefix /= s * exp(t - e);71 }72 //73 // Checked 2F0:74 //75 unsigned k = 0;76 T a1_poch(1 - a);77 T a2_poch(b - a);78 T z_mult(1 / z);79 T sum = 0;80 T abs_sum = 0;81 T term = 1;82 T last_term = 0;83 do84 {85 sum += term;86 last_term = term;87 abs_sum += fabs(sum);88 term *= a1_poch * a2_poch * z_mult;89 term /= ++k;90 a1_poch += 1;91 a2_poch += 1;92 if (fabs(sum) * boost::math::policies::get_epsilon<T, Policy>() > fabs(term))93 break;94 if(fabs(sum) / abs_sum < boost::math::policies::get_epsilon<T, Policy>())95 return boost::math::policies::raise_evaluation_error<T>(function, "Large-z asymptotic approximation to 1F1 has destroyed all the digits in the result due to cancellation. Current best guess is %1%", 96 prefix * sum, Policy());97 if(k > boost::math::policies::get_max_series_iterations<Policy>())98 return boost::math::policies::raise_evaluation_error<T>(function, "1F1: Unable to locate solution in a reasonable time:"99 " large-z asymptotic approximation. Current best guess is %1%", prefix * sum, Policy());100 if((k > 10) && (fabs(term) > fabs(last_term)))101 return boost::math::policies::raise_evaluation_error<T>(function, "Large-z asymptotic approximation to 1F1 is divergent. Current best guess is %1%", prefix * sum, Policy());102 } while (true);103 104 return prefix * sum;105 }106 107 108 // experimental range109 template <class T, class Policy>110 inline bool hypergeometric_1F1_asym_region(const T& a, const T& b, const T& z, const Policy&)111 {112 BOOST_MATH_STD_USING113 int half_digits = policies::digits<T, Policy>() / 2;114 bool in_region = false;115 116 if (fabs(a) < 0.001f)117 return false; // Haven't been able to make this work, why not? TODO!118 119 //120 // We use the following heuristic, if after we have had half_digits terms121 // of the 2F0 series, we require terms to be decreasing in size by a factor122 // of at least 0.7. Assuming the earlier terms were converging much faster123 // than this, then this should be enough to achieve convergence before the124 // series shoots off to infinity.125 //126 if (z > 0)127 {128 T one_minus_a = 1 - a;129 T b_minus_a = b - a;130 if (fabs((one_minus_a + half_digits) * (b_minus_a + half_digits) / (half_digits * z)) < 0.7)131 {132 in_region = true;133 //134 // double check that we are not divergent at the start if a,b < 0:135 //136 if ((one_minus_a < 0) || (b_minus_a < 0))137 {138 if (fabs(one_minus_a * b_minus_a / z) > 0.5)139 in_region = false;140 }141 }142 }143 else if (fabs((1 - (b - a) + half_digits) * (a + half_digits) / (half_digits * z)) < 0.7)144 {145 if ((floor(b - a) == (b - a)) && (b - a < 0))146 return false; // Can't have a negative integer b-a.147 in_region = true;148 //149 // double check that we are not divergent at the start if a,b < 0:150 //151 T a1 = 1 - (b - a);152 if ((a1 < 0) || (a < 0))153 {154 if (fabs(a1 * a / z) > 0.5)155 in_region = false;156 }157 }158 //159 // Check for a and b negative integers as these aren't supported by the approximation:160 //161 if (in_region)162 {163 if ((a < 0) && (floor(a) == a))164 in_region = false;165 if ((b < 0) && (floor(b) == b))166 in_region = false;167 if (fabs(z) < 40)168 in_region = false;169 }170 return in_region;171 }172 173 } } } // namespaces174 175#ifdef _MSC_VER176#pragma warning(pop)177#endif178 179#endif // BOOST_MATH_HYPERGEOMETRIC_ASYM_HPP180