]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
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_HYPERGEOMETRIC_ASYM_HPP | |
11 | #define BOOST_MATH_HYPERGEOMETRIC_ASYM_HPP | |
12 | ||
13 | #include <boost/math/special_functions/gamma.hpp> | |
14 | #include <boost/math/special_functions/hypergeometric_2F0.hpp> | |
15 | ||
1e59de90 | 16 | #ifdef _MSC_VER |
92f5a8d4 TL |
17 | #pragma warning(push) |
18 | #pragma warning(disable:4127) | |
19 | #endif | |
20 | ||
21 | namespace boost { namespace math { | |
22 | ||
23 | namespace detail { | |
24 | ||
25 | // | |
26 | // Asymptotic series based on https://dlmf.nist.gov/13.7#E1 | |
27 | // | |
28 | // Note that a and b must not be negative integers, in addition | |
29 | // we require z > 0 and so apply Kummer's relation for z < 0. | |
30 | // | |
31 | template <class T, class Policy> | |
1e59de90 | 32 | inline T hypergeometric_1F1_asym_large_z_series(T a, const T& b, T z, const Policy& pol, long long& log_scaling) |
92f5a8d4 TL |
33 | { |
34 | BOOST_MATH_STD_USING | |
35 | static const char* function = "boost::math::hypergeometric_1F1_asym_large_z_series<%1%>(%1%, %1%, %1%)"; | |
36 | T prefix; | |
1e59de90 TL |
37 | long long e; |
38 | int s; | |
92f5a8d4 TL |
39 | if (z < 0) |
40 | { | |
41 | a = b - a; | |
42 | z = -z; | |
43 | prefix = 1; | |
44 | } | |
45 | else | |
46 | { | |
1e59de90 | 47 | e = z > static_cast<T>((std::numeric_limits<long long>::max)()) ? (std::numeric_limits<long long>::max)() : lltrunc(z, pol); |
92f5a8d4 TL |
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 | else | |
56 | { | |
57 | T t = log(z) * (a - b); | |
1e59de90 | 58 | e = lltrunc(t, pol); |
92f5a8d4 TL |
59 | log_scaling += e; |
60 | prefix *= exp(t - e); | |
61 | ||
62 | t = boost::math::lgamma(b, &s, pol); | |
1e59de90 | 63 | e = lltrunc(t, pol); |
92f5a8d4 TL |
64 | log_scaling += e; |
65 | prefix *= s * exp(t - e); | |
66 | ||
67 | t = boost::math::lgamma(a, &s, pol); | |
1e59de90 | 68 | e = lltrunc(t, pol); |
92f5a8d4 TL |
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 | do | |
84 | { | |
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 range | |
109 | 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_USING | |
113 | 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 terms | |
121 | // of the 2F0 series, we require terms to be decreasing in size by a factor | |
122 | // of at least 0.7. Assuming the earlier terms were converging much faster | |
123 | // than this, then this should be enough to achieve convergence before the | |
124 | // 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 | } } } // namespaces | |
174 | ||
1e59de90 | 175 | #ifdef _MSC_VER |
92f5a8d4 TL |
176 | #pragma warning(pop) |
177 | #endif | |
178 | ||
179 | #endif // BOOST_MATH_HYPERGEOMETRIC_ASYM_HPP |