]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
1 | |
2 | /////////////////////////////////////////////////////////////////////////////// | |
3 | // Copyright 2014 Anton Bikineev | |
4 | // Copyright 2014 Christopher Kormanyos | |
5 | // Copyright 2014 John Maddock | |
6 | // Copyright 2014 Paul Bristow | |
7 | // Distributed under the Boost | |
8 | // Software License, Version 1.0. (See accompanying file | |
9 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
10 | // | |
11 | #ifndef BOOST_MATH_HYPERGEOMETRIC_1F1_BESSEL_HPP | |
12 | #define BOOST_MATH_HYPERGEOMETRIC_1F1_BESSEL_HPP | |
13 | ||
14 | #include <boost/math/tools/series.hpp> | |
15 | #include <boost/math/special_functions/bessel.hpp> | |
16 | #include <boost/math/special_functions/laguerre.hpp> | |
17 | #include <boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp> | |
18 | #include <boost/math/special_functions/bessel_iterators.hpp> | |
19 | ||
20 | ||
21 | namespace boost { namespace math { namespace detail { | |
22 | ||
23 | template <class T, class Policy> | |
24 | T hypergeometric_1F1_divergent_fallback(const T& a, const T& b, const T& z, const Policy& pol, int& log_scaling); | |
25 | ||
26 | template <class T> | |
27 | bool hypergeometric_1F1_is_tricomi_viable_positive_b(const T& a, const T& b, const T& z) | |
28 | { | |
29 | BOOST_MATH_STD_USING | |
30 | if ((z < b) && (a > -50)) | |
31 | return false; // might as well fall through to recursion | |
32 | if (b <= 100) | |
33 | return true; | |
34 | // Even though we're in a reasonable domain for Tricomi's approximation, | |
35 | // the arguments to the Bessel functions may be so large that we can't | |
36 | // actually evaluate them: | |
37 | T x = sqrt(fabs(2 * z * b - 4 * a * z)); | |
38 | T v = b - 1; | |
39 | return log(boost::math::constants::e<T>() * x / (2 * v)) * v > tools::log_min_value<T>(); | |
40 | } | |
41 | ||
42 | // | |
43 | // Returns an arbitrarily small value compared to "target" for use as a seed | |
44 | // value for Bessel recurrences. Note that we'd better not make it too small | |
45 | // or underflow may occur resulting in either one of the terms in the | |
46 | // recurrence being zero, or else the result being zero. Using 1/epsilon | |
47 | // as a safety factor ensures that if we do underflow to zero, all of the digits | |
48 | // will have been cancelled out anyway: | |
49 | // | |
50 | template <class T> | |
51 | T arbitrary_small_value(const T& target) | |
52 | { | |
53 | using std::fabs; | |
54 | return (tools::min_value<T>() / tools::epsilon<T>()) * (fabs(target) > 1 ? target : 1); | |
55 | } | |
56 | ||
57 | ||
58 | template <class T, class Policy> | |
59 | struct hypergeometric_1F1_AS_13_3_7_tricomi_series | |
60 | { | |
61 | typedef T result_type; | |
62 | ||
63 | enum { cache_size = 64 }; | |
64 | ||
65 | hypergeometric_1F1_AS_13_3_7_tricomi_series(const T& a, const T& b, const T& z, const Policy& pol_) | |
66 | : A_minus_2(1), A_minus_1(0), A(b / 2), mult(z / 2), term(1), b_minus_1_plus_n(b - 1), | |
67 | bessel_arg((b / 2 - a) * z), | |
68 | two_a_minus_b(2 * a - b), pol(pol_), n(2) | |
69 | { | |
70 | BOOST_MATH_STD_USING | |
71 | term /= pow(fabs(bessel_arg), b_minus_1_plus_n / 2); | |
72 | mult /= sqrt(fabs(bessel_arg)); | |
73 | bessel_cache[cache_size - 1] = bessel_arg > 0 ? boost::math::cyl_bessel_j(b_minus_1_plus_n - 1, 2 * sqrt(bessel_arg), pol) : boost::math::cyl_bessel_i(b_minus_1_plus_n - 1, 2 * sqrt(-bessel_arg), pol); | |
74 | if (fabs(bessel_cache[cache_size - 1]) < tools::min_value<T>() / tools::epsilon<T>()) | |
75 | { | |
76 | // We get very limited precision due to rapid denormalisation/underflow of the Bessel values, raise an exception and try something else: | |
77 | policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Underflow in Bessel functions", bessel_cache[cache_size - 1], pol); | |
78 | } | |
79 | if ((term * bessel_cache[cache_size - 1] < tools::min_value<T>() / (tools::epsilon<T>() * tools::epsilon<T>())) || !(boost::math::isfinite)(term) || (!std::numeric_limits<T>::has_infinity && (fabs(term) > tools::max_value<T>()))) | |
80 | { | |
81 | term = -log(fabs(bessel_arg)) * b_minus_1_plus_n / 2; | |
82 | log_scale = itrunc(term); | |
83 | term -= log_scale; | |
84 | term = exp(term); | |
85 | } | |
86 | else | |
87 | log_scale = 0; | |
88 | #ifndef BOOST_NO_CXX17_IF_CONSTEXPR | |
89 | if constexpr (std::numeric_limits<T>::has_infinity) | |
90 | { | |
91 | if (!(boost::math::isfinite)(bessel_cache[cache_size - 1])) | |
92 | policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Expected finite Bessel function result but got %1%", bessel_cache[cache_size - 1], pol); | |
93 | } | |
94 | else | |
95 | if ((boost::math::isnan)(bessel_cache[cache_size - 1]) || (fabs(bessel_cache[cache_size - 1]) >= tools::max_value<T>())) | |
96 | policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Expected finite Bessel function result but got %1%", bessel_cache[cache_size - 1], pol); | |
97 | #else | |
98 | if ((std::numeric_limits<T>::has_infinity && !(boost::math::isfinite)(bessel_cache[cache_size - 1])) | |
99 | || (!std::numeric_limits<T>::has_infinity && ((boost::math::isnan)(bessel_cache[cache_size - 1]) || (fabs(bessel_cache[cache_size - 1]) >= tools::max_value<T>())))) | |
100 | policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Expected finite Bessel function result but got %1%", bessel_cache[cache_size - 1], pol); | |
101 | #endif | |
102 | cache_offset = -cache_size; | |
103 | refill_cache(); | |
104 | } | |
105 | T operator()() | |
106 | { | |
107 | // | |
108 | // We return the n-2 term, and do 2 terms at once as every other term can be | |
109 | // very small (or zero) when b == 2a: | |
110 | // | |
111 | BOOST_MATH_STD_USING | |
112 | if(n - 2 - cache_offset >= cache_size) | |
113 | refill_cache(); | |
114 | T result = A_minus_2 * term * bessel_cache[n - 2 - cache_offset]; | |
115 | term *= mult; | |
116 | ++n; | |
117 | T A_next = ((b_minus_1_plus_n + 2) * A_minus_1 + two_a_minus_b * A_minus_2) / n; | |
118 | b_minus_1_plus_n += 1; | |
119 | A_minus_2 = A_minus_1; | |
120 | A_minus_1 = A; | |
121 | A = A_next; | |
122 | ||
123 | if (A_minus_2 != 0) | |
124 | { | |
125 | if (n - 2 - cache_offset >= cache_size) | |
126 | refill_cache(); | |
127 | result += A_minus_2 * term * bessel_cache[n - 2 - cache_offset]; | |
128 | } | |
129 | term *= mult; | |
130 | ++n; | |
131 | A_next = ((b_minus_1_plus_n + 2) * A_minus_1 + two_a_minus_b * A_minus_2) / n; | |
132 | b_minus_1_plus_n += 1; | |
133 | A_minus_2 = A_minus_1; | |
134 | A_minus_1 = A; | |
135 | A = A_next; | |
136 | ||
137 | return result; | |
138 | } | |
139 | ||
140 | int scale()const | |
141 | { | |
142 | return log_scale; | |
143 | } | |
144 | ||
145 | private: | |
146 | T A_minus_2, A_minus_1, A, mult, term, b_minus_1_plus_n, bessel_arg, two_a_minus_b; | |
147 | std::array<T, cache_size> bessel_cache; | |
148 | const Policy& pol; | |
149 | int n, log_scale, cache_offset; | |
150 | ||
151 | hypergeometric_1F1_AS_13_3_7_tricomi_series operator=(const hypergeometric_1F1_AS_13_3_7_tricomi_series&); | |
152 | ||
153 | void refill_cache() | |
154 | { | |
155 | BOOST_MATH_STD_USING | |
156 | // | |
157 | // We don't calculate a new bessel I/J value: instead start our iterator off | |
158 | // with an arbitrary small value, then when we get back to the last value in the previous cache | |
159 | // calculate the ratio and use it to renormalise all the new values. This is more efficient, but | |
160 | // also avoids problems with J_v(x) or I_v(x) underflowing to zero. | |
161 | // | |
162 | cache_offset += cache_size; | |
163 | T last_value = bessel_cache.back(); | |
164 | T ratio; | |
165 | if (bessel_arg > 0) | |
166 | { | |
167 | // | |
168 | // We will be calculating Bessel J. | |
169 | // We need a different recurrence strategy for positive and negative orders: | |
170 | // | |
171 | if (b_minus_1_plus_n > 0) | |
172 | { | |
173 | bessel_j_backwards_iterator<T> i(b_minus_1_plus_n + (int)cache_size - 1, 2 * sqrt(bessel_arg), arbitrary_small_value(last_value)); | |
174 | ||
175 | for (int j = cache_size - 1; j >= 0; --j, ++i) | |
176 | { | |
177 | bessel_cache[j] = *i; | |
178 | // | |
179 | // Depending on the value of bessel_arg, the values stored in the cache can grow so | |
180 | // large as to overflow, if that looks likely then we need to rescale all the | |
181 | // existing terms (most of which will then underflow to zero). In this situation | |
182 | // it's likely that our series will only need 1 or 2 terms of the series but we | |
183 | // can't be sure of that: | |
184 | // | |
185 | if ((j < cache_size - 2) && (tools::max_value<T>() / fabs(64 * bessel_cache[j] / bessel_cache[j + 1]) < fabs(bessel_cache[j]))) | |
186 | { | |
187 | T rescale = pow(fabs(bessel_cache[j] / bessel_cache[j + 1]), j + 1) * 2; | |
188 | if (!((boost::math::isfinite)(rescale))) | |
189 | rescale = tools::max_value<T>(); | |
190 | for (int k = j; k < cache_size; ++k) | |
191 | bessel_cache[k] /= rescale; | |
192 | bessel_j_backwards_iterator<T> ti(b_minus_1_plus_n + j, 2 * sqrt(bessel_arg), bessel_cache[j + 1], bessel_cache[j]); | |
193 | i = ti; | |
194 | } | |
195 | } | |
196 | ratio = last_value / *i; | |
197 | } | |
198 | else | |
199 | { | |
200 | // | |
201 | // Negative order is difficult: the J_v(x) recurrence relations are unstable | |
202 | // *in both directions* for v < 0, except as v -> -INF when we have | |
203 | // J_-v(x) ~= -sin(pi.v)Y_v(x). | |
204 | // For small v what we can do is compute every other Bessel function and | |
205 | // then fill in the gaps using the recurrence relation. This *is* stable | |
206 | // provided that v is not so negative that the above approximation holds. | |
207 | // | |
208 | bessel_cache[1] = cyl_bessel_j(b_minus_1_plus_n + 1, 2 * sqrt(bessel_arg), pol); | |
209 | bessel_cache[0] = (last_value + bessel_cache[1]) / (b_minus_1_plus_n / sqrt(bessel_arg)); | |
210 | int pos = 2; | |
211 | while ((pos < cache_size - 1) && (b_minus_1_plus_n + pos < 0)) | |
212 | { | |
213 | bessel_cache[pos + 1] = cyl_bessel_j(b_minus_1_plus_n + pos + 1, 2 * sqrt(bessel_arg), pol); | |
214 | bessel_cache[pos] = (bessel_cache[pos-1] + bessel_cache[pos+1]) / ((b_minus_1_plus_n + pos) / sqrt(bessel_arg)); | |
215 | pos += 2; | |
216 | } | |
217 | if (pos < cache_size) | |
218 | { | |
219 | // | |
220 | // We have crossed over into the region where backward recursion is the stable direction | |
221 | // start from arbitrary value and recurse down to "pos" and normalise: | |
222 | // | |
223 | bessel_j_backwards_iterator<T> i2(b_minus_1_plus_n + (int)cache_size - 1, 2 * sqrt(bessel_arg), arbitrary_small_value(bessel_cache[pos-1])); | |
224 | for (int loc = cache_size - 1; loc >= pos; --loc) | |
225 | bessel_cache[loc] = *i2++; | |
226 | ratio = bessel_cache[pos - 1] / *i2; | |
227 | // | |
228 | // Sanity check, if we normalised to an unusually small value then it was likely | |
229 | // to be near a root and the calculated ratio is garbage, if so perform one | |
230 | // more J_v(x) evaluation at position and normalise again: | |
231 | // | |
232 | if (fabs(bessel_cache[pos] * ratio / bessel_cache[pos - 1]) > 5) | |
233 | ratio = cyl_bessel_j(b_minus_1_plus_n + pos, 2 * sqrt(bessel_arg), pol) / bessel_cache[pos]; | |
234 | while (pos < cache_size) | |
235 | bessel_cache[pos++] *= ratio; | |
236 | } | |
237 | ratio = 1; | |
238 | } | |
239 | } | |
240 | else | |
241 | { | |
242 | // | |
243 | // Bessel I. | |
244 | // We need a different recurrence strategy for positive and negative orders: | |
245 | // | |
246 | if (b_minus_1_plus_n > 0) | |
247 | { | |
248 | bessel_i_backwards_iterator<T> i(b_minus_1_plus_n + (int)cache_size - 1, 2 * sqrt(-bessel_arg), arbitrary_small_value(last_value)); | |
249 | ||
250 | for (int j = cache_size - 1; j >= 0; --j, ++i) | |
251 | { | |
252 | bessel_cache[j] = *i; | |
253 | // | |
254 | // Depending on the value of bessel_arg, the values stored in the cache can grow so | |
255 | // large as to overflow, if that looks likely then we need to rescale all the | |
256 | // existing terms (most of which will then underflow to zero). In this situation | |
257 | // it's likely that our series will only need 1 or 2 terms of the series but we | |
258 | // can't be sure of that: | |
259 | // | |
260 | if ((j < cache_size - 2) && (tools::max_value<T>() / fabs(64 * bessel_cache[j] / bessel_cache[j + 1]) < fabs(bessel_cache[j]))) | |
261 | { | |
262 | T rescale = pow(fabs(bessel_cache[j] / bessel_cache[j + 1]), j + 1) * 2; | |
263 | if (!((boost::math::isfinite)(rescale))) | |
264 | rescale = tools::max_value<T>(); | |
265 | for (int k = j; k < cache_size; ++k) | |
266 | bessel_cache[k] /= rescale; | |
267 | i = bessel_i_backwards_iterator<T>(b_minus_1_plus_n + j, 2 * sqrt(-bessel_arg), bessel_cache[j + 1], bessel_cache[j]); | |
268 | } | |
269 | } | |
270 | ratio = last_value / *i; | |
271 | } | |
272 | else | |
273 | { | |
274 | // | |
275 | // Forwards iteration is stable: | |
276 | // | |
277 | bessel_i_forwards_iterator<T> i(b_minus_1_plus_n, 2 * sqrt(-bessel_arg)); | |
278 | int pos = 0; | |
279 | while ((pos < cache_size) && (b_minus_1_plus_n + pos < 0.5)) | |
280 | { | |
281 | bessel_cache[pos++] = *i++; | |
282 | } | |
283 | if (pos < cache_size) | |
284 | { | |
285 | // | |
286 | // We have crossed over into the region where backward recursion is the stable direction | |
287 | // start from arbitrary value and recurse down to "pos" and normalise: | |
288 | // | |
289 | bessel_i_backwards_iterator<T> i2(b_minus_1_plus_n + (int)cache_size - 1, 2 * sqrt(-bessel_arg), arbitrary_small_value(last_value)); | |
290 | for (int loc = cache_size - 1; loc >= pos; --loc) | |
291 | bessel_cache[loc] = *i2++; | |
292 | ratio = bessel_cache[pos - 1] / *i2; | |
293 | while (pos < cache_size) | |
294 | bessel_cache[pos++] *= ratio; | |
295 | } | |
296 | ratio = 1; | |
297 | } | |
298 | } | |
299 | if(ratio != 1) | |
300 | for (auto j = bessel_cache.begin(); j != bessel_cache.end(); ++j) | |
301 | *j *= ratio; | |
302 | // | |
303 | // Very occationally our normalisation fails because the normalisztion value | |
304 | // is sitting right on top of a root (or very close to it). When that happens | |
305 | // best to calculate a fresh Bessel evaluation and normalise again. | |
306 | // | |
307 | if (fabs(bessel_cache[0] / last_value) > 5) | |
308 | { | |
309 | ratio = (bessel_arg < 0 ? cyl_bessel_i(b_minus_1_plus_n, 2 * sqrt(-bessel_arg), pol) : cyl_bessel_j(b_minus_1_plus_n, 2 * sqrt(bessel_arg), pol)) / bessel_cache[0]; | |
310 | if (ratio != 1) | |
311 | for (auto j = bessel_cache.begin(); j != bessel_cache.end(); ++j) | |
312 | *j *= ratio; | |
313 | } | |
314 | } | |
315 | }; | |
316 | ||
317 | template <class T, class Policy> | |
318 | T hypergeometric_1F1_AS_13_3_7_tricomi(const T& a, const T& b, const T& z, const Policy& pol, int& log_scale) | |
319 | { | |
320 | BOOST_MATH_STD_USING | |
321 | // | |
322 | // Works for a < 0, b < 0, z > 0. | |
323 | // | |
324 | // For convergence we require A * term to be converging otherwise we get | |
325 | // a divergent alternating series. It's actually really hard to analyse this | |
326 | // and the best purely heuristic policy we've found is | |
327 | // z < fabs((2 * a - b) / (sqrt(fabs(a)))) ; b > 0 or: | |
328 | // z < fabs((2 * a - b) / (sqrt(fabs(ab)))) ; b < 0 | |
329 | // | |
330 | T prefix(0); | |
331 | int prefix_sgn(0); | |
332 | bool use_logs = false; | |
333 | int scale = 0; | |
334 | // | |
335 | // We can actually support the b == 2a case within here, but we haven't | |
336 | // as we appear never to get here in practice. Which means this get out | |
337 | // clause is a bit of defensive programming.... | |
338 | // | |
339 | if(b == 2 * a) | |
340 | return hypergeometric_1F1_divergent_fallback(a, b, z, pol, log_scale); | |
341 | ||
342 | try | |
343 | { | |
344 | prefix = boost::math::tgamma(b, pol); | |
345 | prefix *= exp(z / 2); | |
346 | } | |
347 | catch (const std::runtime_error&) | |
348 | { | |
349 | use_logs = true; | |
350 | } | |
351 | if (use_logs || (prefix == 0) || !(boost::math::isfinite)(prefix) || (!std::numeric_limits<T>::has_infinity && (fabs(prefix) >= tools::max_value<T>()))) | |
352 | { | |
353 | use_logs = true; | |
354 | prefix = boost::math::lgamma(b, &prefix_sgn, pol) + z / 2; | |
355 | scale = itrunc(prefix); | |
356 | log_scale += scale; | |
357 | prefix -= scale; | |
358 | } | |
359 | T result(0); | |
360 | boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); | |
361 | bool retry = false; | |
362 | int series_scale = 0; | |
363 | try | |
364 | { | |
365 | hypergeometric_1F1_AS_13_3_7_tricomi_series<T, Policy> s(a, b, z, pol); | |
366 | series_scale = s.scale(); | |
367 | log_scale += s.scale(); | |
368 | try | |
369 | { | |
370 | T norm = 0; | |
371 | result = 0; | |
372 | if((a < 0) && (b < 0)) | |
373 | result = boost::math::tools::checked_sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, result, norm); | |
374 | else | |
375 | result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, result); | |
376 | if (!(boost::math::isfinite)(result) || (result == 0) || (!std::numeric_limits<T>::has_infinity && (fabs(result) >= tools::max_value<T>()))) | |
377 | retry = true; | |
378 | if (norm / fabs(result) > 1 / boost::math::tools::root_epsilon<T>()) | |
379 | retry = true; // fatal cancellation | |
380 | } | |
381 | catch (const std::overflow_error&) | |
382 | { | |
383 | retry = true; | |
384 | } | |
385 | catch (const boost::math::evaluation_error&) | |
386 | { | |
387 | retry = true; | |
388 | } | |
389 | } | |
390 | catch (const std::overflow_error&) | |
391 | { | |
392 | log_scale -= scale; | |
393 | return hypergeometric_1F1_divergent_fallback(a, b, z, pol, log_scale); | |
394 | } | |
395 | catch (const boost::math::evaluation_error&) | |
396 | { | |
397 | log_scale -= scale; | |
398 | return hypergeometric_1F1_divergent_fallback(a, b, z, pol, log_scale); | |
399 | } | |
400 | if (retry) | |
401 | { | |
402 | log_scale -= scale; | |
403 | log_scale -= series_scale; | |
404 | return hypergeometric_1F1_divergent_fallback(a, b, z, pol, log_scale); | |
405 | } | |
406 | boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_AS_13_3_7<%1%>(%1%,%1%,%1%)", max_iter, pol); | |
407 | if (use_logs) | |
408 | { | |
409 | int sgn = boost::math::sign(result); | |
410 | prefix += log(fabs(result)); | |
411 | result = sgn * prefix_sgn * exp(prefix); | |
412 | } | |
413 | else | |
414 | { | |
415 | if ((fabs(result) > 1) && (fabs(prefix) > 1) && (tools::max_value<T>() / fabs(result) < fabs(prefix))) | |
416 | { | |
417 | // Overflow: | |
418 | scale = itrunc(tools::log_max_value<T>()) - 10; | |
419 | log_scale += scale; | |
420 | result /= exp(T(scale)); | |
421 | } | |
422 | result *= prefix; | |
423 | } | |
424 | return result; | |
425 | } | |
426 | ||
427 | ||
428 | template <class T> | |
429 | struct cyl_bessel_i_large_x_sum | |
430 | { | |
431 | typedef T result_type; | |
432 | ||
433 | cyl_bessel_i_large_x_sum(const T& v, const T& x) : v(v), z(x), term(1), k(0) {} | |
434 | ||
435 | T operator()() | |
436 | { | |
437 | T result = term; | |
438 | ++k; | |
439 | term *= -(4 * v * v - (2 * k - 1) * (2 * k - 1)) / (8 * k * z); | |
440 | return result; | |
441 | } | |
442 | T v, z, term; | |
443 | int k; | |
444 | }; | |
445 | ||
446 | template <class T, class Policy> | |
447 | T cyl_bessel_i_large_x_scaled(const T& v, const T& x, int& log_scaling, const Policy& pol) | |
448 | { | |
449 | BOOST_MATH_STD_USING | |
450 | cyl_bessel_i_large_x_sum<T> s(v, x); | |
451 | boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); | |
452 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); | |
453 | boost::math::policies::check_series_iterations<T>("boost::math::cyl_bessel_i_large_x<%1%>(%1%,%1%)", max_iter, pol); | |
454 | int scale = boost::math::itrunc(x); | |
455 | log_scaling += scale; | |
456 | return result * exp(x - scale) / sqrt(boost::math::constants::two_pi<T>() * x); | |
457 | } | |
458 | ||
459 | ||
460 | ||
461 | template <class T, class Policy> | |
462 | struct hypergeometric_1F1_AS_13_3_6_series | |
463 | { | |
464 | typedef T result_type; | |
465 | ||
466 | enum { cache_size = 64 }; | |
467 | // | |
468 | // This series is only convergent/useful for a and b approximately equal | |
469 | // (ideally |a-b| < 1). The series can also go divergent after a while | |
470 | // when b < 0, which limits precision to around that of double. In that | |
471 | // situation we return 0 to terminate the series as otherwise the divergent | |
472 | // terms will destroy all the bits in our result before they do eventually | |
473 | // converge again. One important use case for this series is for z < 0 | |
474 | // and |a| << |b| so that either b-a == b or at least most of the digits in a | |
475 | // are lost in the subtraction. Note that while you can easily convince yourself | |
476 | // that the result should be unity when b-a == b, in fact this is not (quite) | |
477 | // the case for large z. | |
478 | // | |
479 | hypergeometric_1F1_AS_13_3_6_series(const T& a, const T& b, const T& z, const T& b_minus_a, const Policy& pol_) | |
480 | : b_minus_a(b_minus_a), half_z(z / 2), poch_1(2 * b_minus_a - 1), poch_2(b_minus_a - a), b_poch(b), term(1), last_result(1), sign(1), n(0), cache_offset(-cache_size), scale(0), pol(pol_) | |
481 | { | |
482 | bessel_i_cache[cache_size - 1] = half_z > tools::log_max_value<T>() ? | |
483 | cyl_bessel_i_large_x_scaled(T(b_minus_a - 1.5f), half_z, scale, pol) : boost::math::cyl_bessel_i(b_minus_a - 1.5f, half_z, pol); | |
484 | refill_cache(); | |
485 | } | |
486 | T operator()() | |
487 | { | |
488 | BOOST_MATH_STD_USING | |
489 | if(n - cache_offset >= cache_size) | |
490 | refill_cache(); | |
491 | ||
492 | T result = term * (b_minus_a - 0.5f + n) * sign * bessel_i_cache[n - cache_offset]; | |
493 | ++n; | |
494 | term *= poch_1; | |
495 | poch_1 = (n == 1) ? T(2 * b_minus_a) : T(poch_1 + 1); | |
496 | term *= poch_2; | |
497 | poch_2 += 1; | |
498 | term /= n; | |
499 | term /= b_poch; | |
500 | b_poch += 1; | |
501 | sign = -sign; | |
502 | ||
503 | if ((fabs(result) > fabs(last_result)) && (n > 100)) | |
504 | return 0; // series has gone divergent! | |
505 | ||
506 | last_result = result; | |
507 | ||
508 | return result; | |
509 | } | |
510 | ||
511 | int scaling()const | |
512 | { | |
513 | return scale; | |
514 | } | |
515 | ||
516 | private: | |
517 | T b_minus_a, half_z, poch_1, poch_2, b_poch, term, last_result; | |
518 | int sign; | |
519 | int n, cache_offset, scale; | |
520 | const Policy& pol; | |
521 | std::array<T, cache_size> bessel_i_cache; | |
522 | ||
523 | void refill_cache() | |
524 | { | |
525 | BOOST_MATH_STD_USING | |
526 | // | |
527 | // We don't calculate a new bessel I value: instead start our iterator off | |
528 | // with an arbitrary small value, then when we get back to the last value in the previous cache | |
529 | // calculate the ratio and use it to renormalise all the values. This is more efficient, but | |
530 | // also avoids problems with I_v(x) underflowing to zero. | |
531 | // | |
532 | cache_offset += cache_size; | |
533 | T last_value = bessel_i_cache.back(); | |
534 | bessel_i_backwards_iterator<T> i(b_minus_a + cache_offset + (int)cache_size - 1.5f, half_z, tools::min_value<T>() * (fabs(last_value) > 1 ? last_value : 1) / tools::epsilon<T>()); | |
535 | ||
536 | for (int j = cache_size - 1; j >= 0; --j, ++i) | |
537 | { | |
538 | bessel_i_cache[j] = *i; | |
539 | // | |
540 | // Depending on the value of half_z, the values stored in the cache can grow so | |
541 | // large as to overflow, if that looks likely then we need to rescale all the | |
542 | // existing terms (most of which will then underflow to zero). In this situation | |
543 | // it's likely that our series will only need 1 or 2 terms of the series but we | |
544 | // can't be sure of that: | |
545 | // | |
546 | if((j < cache_size - 2) && (bessel_i_cache[j + 1] != 0) && (tools::max_value<T>() / fabs(64 * bessel_i_cache[j] / bessel_i_cache[j + 1]) < fabs(bessel_i_cache[j]))) | |
547 | { | |
548 | T rescale = pow(fabs(bessel_i_cache[j] / bessel_i_cache[j + 1]), j + 1) * 2; | |
549 | if (rescale > tools::max_value<T>()) | |
550 | rescale = tools::max_value<T>(); | |
551 | for (int k = j; k < cache_size; ++k) | |
552 | bessel_i_cache[k] /= rescale; | |
553 | i = bessel_i_backwards_iterator<T>(b_minus_a -0.5f + cache_offset + j, half_z, bessel_i_cache[j + 1], bessel_i_cache[j]); | |
554 | } | |
555 | } | |
556 | T ratio = last_value / *i; | |
557 | for (auto j = bessel_i_cache.begin(); j != bessel_i_cache.end(); ++j) | |
558 | *j *= ratio; | |
559 | } | |
560 | ||
561 | hypergeometric_1F1_AS_13_3_6_series(); | |
562 | hypergeometric_1F1_AS_13_3_6_series(const hypergeometric_1F1_AS_13_3_6_series&); | |
563 | hypergeometric_1F1_AS_13_3_6_series& operator=(const hypergeometric_1F1_AS_13_3_6_series&); | |
564 | }; | |
565 | ||
566 | template <class T, class Policy> | |
567 | T hypergeometric_1F1_AS_13_3_6(const T& a, const T& b, const T& z, const T& b_minus_a, const Policy& pol, int& log_scaling) | |
568 | { | |
569 | BOOST_MATH_STD_USING | |
570 | if(b_minus_a == 0) | |
571 | { | |
572 | // special case: M(a,a,z) == exp(z) | |
573 | int scale = itrunc(z, pol); | |
574 | log_scaling += scale; | |
575 | return exp(z - scale); | |
576 | } | |
577 | hypergeometric_1F1_AS_13_3_6_series<T, Policy> s(a, b, z, b_minus_a, pol); | |
578 | boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); | |
579 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); | |
580 | boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_AS_13_3_6<%1%>(%1%,%1%,%1%)", max_iter, pol); | |
581 | result *= boost::math::tgamma(b_minus_a - 0.5f) * pow(z / 4, -b_minus_a + 0.5f); | |
582 | int scale = itrunc(z / 2); | |
583 | log_scaling += scale; | |
584 | log_scaling += s.scaling(); | |
585 | result *= exp(z / 2 - scale); | |
586 | return result; | |
587 | } | |
588 | ||
589 | /****************************************************************************************************************/ | |
590 | // | |
591 | // The following are not used at present and are commented out for that reason: | |
592 | // | |
593 | /****************************************************************************************************************/ | |
594 | ||
595 | #if 0 | |
596 | ||
597 | template <class T, class Policy> | |
598 | struct hypergeometric_1F1_AS_13_3_8_series | |
599 | { | |
600 | // | |
601 | // TODO: store and cache Bessel function evaluations via backwards recurrence. | |
602 | // | |
603 | // The C term grows by at least an order of magnitude with each iteration, and | |
604 | // rate of growth is largely independent of the arguments. Free parameter h | |
605 | // seems to give accurate results for small values (almost zero) or h=1. | |
606 | // Convergence and accuracy, only when -a/z > 100, this appears to have no | |
607 | // or little benefit over 13.3.7 as it generally requires more iterations? | |
608 | // | |
609 | hypergeometric_1F1_AS_13_3_8_series(const T& a, const T& b, const T& z, const T& h, const Policy& pol_) | |
610 | : C_minus_2(1), C_minus_1(-b * h), C(b * (b + 1) * h * h / 2 - (2 * h - 1) * a / 2), | |
611 | bessel_arg(2 * sqrt(-a * z)), bessel_order(b - 1), power_term(std::pow(-a * z, (1 - b) / 2)), mult(z / std::sqrt(-a * z)), | |
612 | a_(a), b_(b), z_(z), h_(h), n(2), pol(pol_) | |
613 | { | |
614 | } | |
615 | T operator()() | |
616 | { | |
617 | // we actually return the n-2 term: | |
618 | T result = C_minus_2 * power_term * boost::math::cyl_bessel_j(bessel_order, bessel_arg, pol); | |
619 | bessel_order += 1; | |
620 | power_term *= mult; | |
621 | ++n; | |
622 | T C_next = ((1 - 2 * h_) * (n - 1) - b_ * h_) * C | |
623 | + ((1 - 2 * h_) * a_ - h_ * (h_ - 1) *(b_ + n - 2)) * C_minus_1 | |
624 | - h_ * (h_ - 1) * a_ * C_minus_2; | |
625 | C_next /= n; | |
626 | C_minus_2 = C_minus_1; | |
627 | C_minus_1 = C; | |
628 | C = C_next; | |
629 | return result; | |
630 | } | |
631 | T C, C_minus_1, C_minus_2, bessel_arg, bessel_order, power_term, mult, a_, b_, z_, h_; | |
632 | const Policy& pol; | |
633 | int n; | |
634 | ||
635 | typedef T result_type; | |
636 | }; | |
637 | ||
638 | template <class T, class Policy> | |
639 | T hypergeometric_1F1_AS_13_3_8(const T& a, const T& b, const T& z, const T& h, const Policy& pol) | |
640 | { | |
641 | BOOST_MATH_STD_USING | |
642 | T prefix = exp(h * z) * boost::math::tgamma(b); | |
643 | hypergeometric_1F1_AS_13_3_8_series<T, Policy> s(a, b, z, h, pol); | |
644 | boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); | |
645 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); | |
646 | boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_AS_13_3_8<%1%>(%1%,%1%,%1%)", max_iter, pol); | |
647 | result *= prefix; | |
648 | return result; | |
649 | } | |
650 | ||
651 | // | |
652 | // This is the series from https://dlmf.nist.gov/13.11 | |
653 | // It appears to be unusable for a,z < 0, and for | |
654 | // b < 0 appears to never be better than the defining series | |
655 | // for 1F1. | |
656 | // | |
657 | template <class T, class Policy> | |
658 | struct hypergeometric_1f1_13_11_1_series | |
659 | { | |
660 | typedef T result_type; | |
661 | ||
662 | hypergeometric_1f1_13_11_1_series(const T& a, const T& b, const T& z, const Policy& pol_) | |
663 | : term(1), two_a_minus_1_plus_s(2 * a - 1), two_a_minus_b_plus_s(2 * a - b), b_plus_s(b), a_minus_half_plus_s(a - 0.5f), half_z(z / 2), s(0), pol(pol_) | |
664 | { | |
665 | } | |
666 | T operator()() | |
667 | { | |
668 | T result = term * a_minus_half_plus_s * boost::math::cyl_bessel_i(a_minus_half_plus_s, half_z, pol); | |
669 | ||
670 | term *= two_a_minus_1_plus_s * two_a_minus_b_plus_s / (b_plus_s * ++s); | |
671 | two_a_minus_1_plus_s += 1; | |
672 | a_minus_half_plus_s += 1; | |
673 | two_a_minus_b_plus_s += 1; | |
674 | b_plus_s += 1; | |
675 | ||
676 | return result; | |
677 | } | |
678 | T term, two_a_minus_1_plus_s, two_a_minus_b_plus_s, b_plus_s, a_minus_half_plus_s, half_z; | |
679 | int s; | |
680 | const Policy& pol; | |
681 | }; | |
682 | ||
683 | template <class T, class Policy> | |
684 | T hypergeometric_1f1_13_11_1(const T& a, const T& b, const T& z, const Policy& pol, int& log_scaling) | |
685 | { | |
686 | BOOST_MATH_STD_USING | |
687 | bool use_logs = false; | |
688 | T prefix; | |
689 | int prefix_sgn = 1; | |
690 | if (true/*(a < boost::math::max_factorial<T>::value) && (a > 0)*/) | |
691 | prefix = boost::math::tgamma(a - 0.5f, pol); | |
692 | else | |
693 | { | |
694 | prefix = boost::math::lgamma(a - 0.5f, &prefix_sgn, pol); | |
695 | use_logs = true; | |
696 | } | |
697 | if (use_logs) | |
698 | { | |
699 | prefix += z / 2; | |
700 | prefix += log(z / 4) * (0.5f - a); | |
701 | } | |
702 | else if (z > 0) | |
703 | { | |
704 | prefix *= pow(z / 4, 0.5f - a); | |
705 | prefix *= exp(z / 2); | |
706 | } | |
707 | else | |
708 | { | |
709 | prefix *= exp(z / 2); | |
710 | prefix *= pow(z / 4, 0.5f - a); | |
711 | } | |
712 | ||
713 | hypergeometric_1f1_13_11_1_series<T, Policy> s(a, b, z, pol); | |
714 | boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); | |
715 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); | |
716 | boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_13_11_1<%1%>(%1%,%1%,%1%)", max_iter, pol); | |
717 | if (use_logs) | |
718 | { | |
719 | int scaling = itrunc(prefix); | |
720 | log_scaling += scaling; | |
721 | prefix -= scaling; | |
722 | result *= exp(prefix) * prefix_sgn; | |
723 | } | |
724 | else | |
725 | result *= prefix; | |
726 | ||
727 | return result; | |
728 | } | |
729 | ||
730 | #endif | |
731 | ||
732 | } } } // namespaces | |
733 | ||
734 | #endif // BOOST_MATH_HYPERGEOMETRIC_1F1_BESSEL_HPP |