]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
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_MATH_HYPERGEOMETRIC_1F1_ADDITION_THEOREMS_ON_Z_HPP | |
9 | #define BOOST_MATH_HYPERGEOMETRIC_1F1_ADDITION_THEOREMS_ON_Z_HPP | |
10 | ||
11 | #include <boost/math/tools/series.hpp> | |
12 | ||
13 | // | |
14 | // This file implements the addition theorems for 1F1 on z, specifically | |
15 | // each function returns 1F1[a, b, z + k] for some integer k - there's | |
16 | // no particular reason why k needs to be an integer, but no reason why | |
17 | // it shouldn't be either. | |
18 | // | |
19 | // The functions are named hypergeometric_1f1_recurrence_on_z_[plus|minus|zero]_[plus|minus|zero] | |
20 | // where a "plus" indicates forward recurrence, minus backwards recurrence, and zero no recurrence. | |
21 | // So for example hypergeometric_1f1_recurrence_on_z_zero_plus uses forward recurrence on b and | |
22 | // hypergeometric_1f1_recurrence_on_z_minus_minus uses backwards recurrence on both a and b. | |
23 | // | |
24 | // See https://dlmf.nist.gov/13.13 | |
25 | // | |
26 | ||
27 | namespace boost { namespace math { namespace detail { | |
28 | ||
29 | // | |
30 | // This works moderately well for a < 0, but has some very strange behaviour with | |
31 | // strings of values of the same sign followed by a sign switch then another | |
32 | // series all the same sign and so on.... doesn't converge smoothly either | |
33 | // but rises and falls in wave-like behaviour.... very slow to converge... | |
34 | // | |
35 | template <class T, class Policy> | |
36 | struct hypergeometric_1f1_recurrence_on_z_minus_zero_series | |
37 | { | |
38 | typedef T result_type; | |
39 | ||
40 | hypergeometric_1f1_recurrence_on_z_minus_zero_series(const T& a, const T& b, const T& z, int k_, const Policy& pol) | |
41 | : term(1), b_minus_a_plus_n(b - a), a_(a), b_(b), z_(z), n(0), k(k_) | |
42 | { | |
43 | BOOST_MATH_STD_USING | |
44 | int scale1(0), scale2(0); | |
45 | M = boost::math::detail::hypergeometric_1F1_imp(a, b, z, pol, scale1); | |
46 | M_next = boost::math::detail::hypergeometric_1F1_imp(T(a - 1), b, z, pol, scale2); | |
47 | if (scale1 != scale2) | |
48 | M_next *= exp(scale2 - scale1); | |
49 | if (M > 1e10f) | |
50 | { | |
51 | // rescale: | |
52 | int rescale = itrunc(log(fabs(M))); | |
53 | M *= exp(T(-rescale)); | |
54 | M_next *= exp(T(-rescale)); | |
55 | scale1 += rescale; | |
56 | } | |
57 | scaling = scale1; | |
58 | } | |
59 | T operator()() | |
60 | { | |
61 | T result = term * M; | |
62 | term *= b_minus_a_plus_n * k / ((z_ + k) * ++n); | |
63 | b_minus_a_plus_n += 1; | |
64 | T M2 = -((2 * (a_ - n) - b_ + z_) * M_next - (a_ - n) * M) / (b_ - (a_ - n)); | |
65 | M = M_next; | |
66 | M_next = M2; | |
67 | ||
68 | return result; | |
69 | } | |
70 | int scale()const { return scaling; } | |
71 | private: | |
72 | T term, b_minus_a_plus_n, M, M_next, a_, b_, z_; | |
73 | int n, k, scaling; | |
74 | }; | |
75 | ||
76 | template <class T, class Policy> | |
77 | T hypergeometric_1f1_recurrence_on_z_minus_zero(const T& a, const T& b, const T& z, int k, const Policy& pol, int& log_scaling) | |
78 | { | |
79 | BOOST_MATH_STD_USING | |
80 | BOOST_ASSERT((z + k) / z > 0.5f); | |
81 | hypergeometric_1f1_recurrence_on_z_minus_zero_series<T, Policy> s(a, b, z, k, pol); | |
82 | boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); | |
83 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); | |
84 | log_scaling += s.scale(); | |
85 | boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_recurrence_on_z_plus_plus<%1%>(%1%,%1%,%1%)", max_iter, pol); | |
86 | return result * exp(T(k)) * pow(z / (z + k), b - a); | |
87 | } | |
88 | ||
89 | #if 0 | |
90 | ||
91 | // | |
92 | // These are commented out as they are currently unused, but may find use in the future: | |
93 | // | |
94 | ||
95 | template <class T, class Policy> | |
96 | struct hypergeometric_1f1_recurrence_on_z_plus_plus_series | |
97 | { | |
98 | typedef T result_type; | |
99 | ||
100 | hypergeometric_1f1_recurrence_on_z_plus_plus_series(const T& a, const T& b, const T& z, int k_, const Policy& pol) | |
101 | : term(1), a_plus_n(a), b_plus_n(b), z_(z), n(0), k(k_) | |
102 | { | |
103 | M = boost::math::detail::hypergeometric_1F1_imp(a, b, z, pol); | |
104 | M_next = boost::math::detail::hypergeometric_1F1_imp(a + 1, b + 1, z, pol); | |
105 | } | |
106 | T operator()() | |
107 | { | |
108 | T result = term * M; | |
109 | term *= a_plus_n * k / (b_plus_n * ++n); | |
110 | a_plus_n += 1; | |
111 | b_plus_n += 1; | |
112 | // The a_plus_n == 0 case below isn't actually correct, but doesn't matter as that term will be zero | |
f67539c2 | 113 | // anyway, we just need to not divide by zero and end up with a NaN in the result. |
92f5a8d4 TL |
114 | T M2 = (a_plus_n == -1) ? 1 : (a_plus_n == 0) ? 0 : (M_next * b_plus_n * (1 - b_plus_n + z_) + b_plus_n * (b_plus_n - 1) * M) / (a_plus_n * z_); |
115 | M = M_next; | |
116 | M_next = M2; | |
117 | return result; | |
118 | } | |
119 | T term, a_plus_n, b_plus_n, M, M_next, z_; | |
120 | int n, k; | |
121 | }; | |
122 | ||
123 | template <class T, class Policy> | |
124 | T hypergeometric_1f1_recurrence_on_z_plus_plus(const T& a, const T& b, const T& z, int k, const Policy& pol) | |
125 | { | |
126 | hypergeometric_1f1_recurrence_on_z_plus_plus_series<T, Policy> s(a, b, z, k, pol); | |
127 | boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); | |
128 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); | |
129 | boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_recurrence_on_z_plus_plus<%1%>(%1%,%1%,%1%)", max_iter, pol); | |
130 | return result; | |
131 | } | |
132 | ||
133 | template <class T, class Policy> | |
134 | struct hypergeometric_1f1_recurrence_on_z_zero_minus_series | |
135 | { | |
136 | typedef T result_type; | |
137 | ||
138 | hypergeometric_1f1_recurrence_on_z_zero_minus_series(const T& a, const T& b, const T& z, int k_, const Policy& pol) | |
139 | : term(1), b_pochhammer(1 - b), x_k_power(-k_ / z), b_minus_n(b), a_(a), z_(z), b_(b), n(0), k(k_) | |
140 | { | |
141 | M = boost::math::detail::hypergeometric_1F1_imp(a, b, z, pol); | |
142 | M_next = boost::math::detail::hypergeometric_1F1_imp(a, b - 1, z, pol); | |
143 | } | |
144 | T operator()() | |
145 | { | |
146 | BOOST_MATH_STD_USING | |
147 | T result = term * M; | |
148 | term *= b_pochhammer * x_k_power / ++n; | |
149 | b_pochhammer += 1; | |
150 | b_minus_n -= 1; | |
151 | T M2 = (M_next * b_minus_n * (1 - b_minus_n - z_) + z_ * (b_minus_n - a_) * M) / (-b_minus_n * (b_minus_n - 1)); | |
152 | M = M_next; | |
153 | M_next = M2; | |
154 | return result; | |
155 | } | |
156 | T term, b_pochhammer, x_k_power, M, M_next, b_minus_n, a_, z_, b_; | |
157 | int n, k; | |
158 | }; | |
159 | ||
160 | template <class T, class Policy> | |
161 | T hypergeometric_1f1_recurrence_on_z_zero_minus(const T& a, const T& b, const T& z, int k, const Policy& pol) | |
162 | { | |
163 | BOOST_MATH_STD_USING | |
164 | BOOST_ASSERT(abs(k) < fabs(z)); | |
165 | hypergeometric_1f1_recurrence_on_z_zero_minus_series<T, Policy> s(a, b, z, k, pol); | |
166 | boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); | |
167 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); | |
168 | boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_recurrence_on_z_plus_plus<%1%>(%1%,%1%,%1%)", max_iter, pol); | |
169 | return result * pow((z + k) / z, 1 - b); | |
170 | } | |
171 | ||
172 | template <class T, class Policy> | |
173 | struct hypergeometric_1f1_recurrence_on_z_plus_zero_series | |
174 | { | |
175 | typedef T result_type; | |
176 | ||
177 | hypergeometric_1f1_recurrence_on_z_plus_zero_series(const T& a, const T& b, const T& z, int k_, const Policy& pol) | |
178 | : term(1), a_pochhammer(a), z_plus_k(z + k_), b_(b), a_(a), z_(z), n(0), k(k_) | |
179 | { | |
180 | M = boost::math::detail::hypergeometric_1F1_imp(a, b, z, pol); | |
181 | M_next = boost::math::detail::hypergeometric_1F1_imp(a + 1, b, z, pol); | |
182 | } | |
183 | T operator()() | |
184 | { | |
185 | T result = term * M; | |
186 | term *= a_pochhammer * k / (++n * z_plus_k); | |
187 | a_pochhammer += 1; | |
188 | T M2 = (a_pochhammer == -1) ? 1 : (a_pochhammer == 0) ? 0 : (M_next * (2 * a_pochhammer - b_ + z_) + (b_ - a_pochhammer) * M) / a_pochhammer; | |
189 | M = M_next; | |
190 | M_next = M2; | |
191 | ||
192 | return result; | |
193 | } | |
194 | T term, a_pochhammer, z_plus_k, M, M_next, b_minus_n, a_, b_, z_; | |
195 | int n, k; | |
196 | }; | |
197 | ||
198 | template <class T, class Policy> | |
199 | T hypergeometric_1f1_recurrence_on_z_plus_zero(const T& a, const T& b, const T& z, int k, const Policy& pol) | |
200 | { | |
201 | BOOST_MATH_STD_USING | |
202 | BOOST_ASSERT(k / z > -0.5f); | |
203 | //BOOST_ASSERT(floor(a) != a || a > 0); | |
204 | hypergeometric_1f1_recurrence_on_z_plus_zero_series<T, Policy> s(a, b, z, k, pol); | |
205 | boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); | |
206 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); | |
207 | boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_recurrence_on_z_plus_plus<%1%>(%1%,%1%,%1%)", max_iter, pol); | |
208 | return result * pow(z / (z + k), a); | |
209 | } | |
210 | ||
211 | template <class T, class Policy> | |
212 | struct hypergeometric_1f1_recurrence_on_z_zero_plus_series | |
213 | { | |
214 | typedef T result_type; | |
215 | ||
216 | hypergeometric_1f1_recurrence_on_z_zero_plus_series(const T& a, const T& b, const T& z, int k_, const Policy& pol) | |
217 | : term(1), b_minus_a_plus_n(b - a), b_plus_n(b), a_(a), z_(z), n(0), k(k_) | |
218 | { | |
219 | M = boost::math::detail::hypergeometric_1F1_imp(a, b, z, pol); | |
220 | M_next = boost::math::detail::hypergeometric_1F1_imp(a, b + 1, z, pol); | |
221 | } | |
222 | T operator()() | |
223 | { | |
224 | T result = term * M; | |
225 | term *= b_minus_a_plus_n * -k / (b_plus_n * ++n); | |
226 | b_minus_a_plus_n += 1; | |
227 | b_plus_n += 1; | |
228 | T M2 = (b_plus_n * (b_plus_n - 1) * M + b_plus_n * (1 - b_plus_n - z_) * M_next) / (-z_ * b_minus_a_plus_n); | |
229 | M = M_next; | |
230 | M_next = M2; | |
231 | ||
232 | return result; | |
233 | } | |
234 | T term, b_minus_a_plus_n, M, M_next, b_minus_n, a_, b_plus_n, z_; | |
235 | int n, k; | |
236 | }; | |
237 | ||
238 | template <class T, class Policy> | |
239 | T hypergeometric_1f1_recurrence_on_z_zero_plus(const T& a, const T& b, const T& z, int k, const Policy& pol) | |
240 | { | |
241 | BOOST_MATH_STD_USING | |
242 | hypergeometric_1f1_recurrence_on_z_zero_plus_series<T, Policy> s(a, b, z, k, pol); | |
243 | boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); | |
244 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); | |
245 | boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_recurrence_on_z_plus_plus<%1%>(%1%,%1%,%1%)", max_iter, pol); | |
246 | return result * exp(T(k)); | |
247 | } | |
248 | // | |
249 | // I'm unable to find any situation where this series isn't divergent and therefore | |
250 | // is probably quite useless: | |
251 | // | |
252 | template <class T, class Policy> | |
253 | struct hypergeometric_1f1_recurrence_on_z_minus_minus_series | |
254 | { | |
255 | typedef T result_type; | |
256 | ||
257 | hypergeometric_1f1_recurrence_on_z_minus_minus_series(const T& a, const T& b, const T& z, int k_, const Policy& pol) | |
258 | : term(1), one_minus_b_plus_n(1 - b), a_(a), b_(b), z_(z), n(0), k(k_) | |
259 | { | |
260 | M = boost::math::detail::hypergeometric_1F1_imp(a, b, z, pol); | |
261 | M_next = boost::math::detail::hypergeometric_1F1_imp(a - 1, b - 1, z, pol); | |
262 | } | |
263 | T operator()() | |
264 | { | |
265 | T result = term * M; | |
266 | term *= one_minus_b_plus_n * k / (z_ * ++n); | |
267 | one_minus_b_plus_n += 1; | |
268 | T M2 = -((b_ - n) * (1 - b_ + n + z_) * M_next - (a_ - n) * z_ * M) / ((b_ - n) * (b_ - n - 1)); | |
269 | M = M_next; | |
270 | M_next = M2; | |
271 | ||
272 | return result; | |
273 | } | |
274 | T term, one_minus_b_plus_n, M, M_next, a_, b_, z_; | |
275 | int n, k; | |
276 | }; | |
277 | ||
278 | template <class T, class Policy> | |
279 | T hypergeometric_1f1_recurrence_on_z_minus_minus(const T& a, const T& b, const T& z, int k, const Policy& pol) | |
280 | { | |
281 | BOOST_MATH_STD_USING | |
282 | hypergeometric_1f1_recurrence_on_z_minus_minus_series<T, Policy> s(a, b, z, k, pol); | |
283 | boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); | |
284 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); | |
285 | boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_recurrence_on_z_plus_plus<%1%>(%1%,%1%,%1%)", max_iter, pol); | |
286 | return result * exp(T(k)) * pow((z + k) / z, 1 - b); | |
287 | } | |
288 | #endif | |
289 | ||
290 | } } } // namespaces | |
291 | ||
292 | #endif // BOOST_MATH_HYPERGEOMETRIC_1F1_ADDITION_THEOREMS_ON_Z_HPP |