]>
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_HYPERGEOMETRIC_1F1_LARGE_ABZ_HPP_ | |
9 | #define BOOST_HYPERGEOMETRIC_1F1_LARGE_ABZ_HPP_ | |
10 | ||
11 | #include <boost/math/special_functions/detail/hypergeometric_1F1_bessel.hpp> | |
12 | #include <boost/math/special_functions/detail/hypergeometric_series.hpp> | |
13 | #include <boost/math/special_functions/gamma.hpp> | |
1e59de90 | 14 | #include <boost/math/special_functions/trunc.hpp> |
92f5a8d4 TL |
15 | |
16 | namespace boost { namespace math { namespace detail { | |
17 | ||
18 | template <class T> | |
19 | inline bool is_negative_integer(const T& x) | |
20 | { | |
21 | using std::floor; | |
22 | return (x <= 0) && (floor(x) == x); | |
23 | } | |
24 | ||
25 | ||
26 | template <class T, class Policy> | |
27 | struct hypergeometric_1F1_igamma_series | |
28 | { | |
29 | enum{ cache_size = 64 }; | |
30 | ||
31 | typedef T result_type; | |
32 | hypergeometric_1F1_igamma_series(const T& alpha, const T& delta, const T& x, const Policy& pol) | |
33 | : delta_poch(-delta), alpha_poch(alpha), x(x), k(0), cache_offset(0), pol(pol) | |
34 | { | |
35 | BOOST_MATH_STD_USING | |
36 | T log_term = log(x) * -alpha; | |
1e59de90 | 37 | log_scaling = lltrunc(log_term - 3 - boost::math::tools::log_min_value<T>() / 50); |
92f5a8d4 TL |
38 | term = exp(log_term - log_scaling); |
39 | refill_cache(); | |
40 | } | |
41 | T operator()() | |
42 | { | |
43 | if (k - cache_offset >= cache_size) | |
44 | { | |
45 | cache_offset += cache_size; | |
46 | refill_cache(); | |
47 | } | |
48 | T result = term * gamma_cache[k - cache_offset]; | |
49 | term *= delta_poch * alpha_poch / (++k * x); | |
50 | delta_poch += 1; | |
51 | alpha_poch += 1; | |
52 | return result; | |
53 | } | |
54 | void refill_cache() | |
55 | { | |
56 | typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; | |
57 | ||
58 | gamma_cache[cache_size - 1] = boost::math::gamma_p(alpha_poch + ((int)cache_size - 1), x, pol); | |
59 | for (int i = cache_size - 1; i > 0; --i) | |
60 | { | |
61 | gamma_cache[i - 1] = gamma_cache[i] >= 1 ? T(1) : T(gamma_cache[i] + regularised_gamma_prefix(T(alpha_poch + (i - 1)), x, pol, lanczos_type()) / (alpha_poch + (i - 1))); | |
62 | } | |
63 | } | |
64 | T delta_poch, alpha_poch, x, term; | |
65 | T gamma_cache[cache_size]; | |
66 | int k; | |
1e59de90 | 67 | long long log_scaling; |
92f5a8d4 TL |
68 | int cache_offset; |
69 | Policy pol; | |
70 | }; | |
71 | ||
72 | template <class T, class Policy> | |
1e59de90 | 73 | T hypergeometric_1F1_igamma(const T& a, const T& b, const T& x, const T& b_minus_a, const Policy& pol, long long& log_scaling) |
92f5a8d4 TL |
74 | { |
75 | BOOST_MATH_STD_USING | |
76 | if (b_minus_a == 0) | |
77 | { | |
78 | // special case: M(a,a,z) == exp(z) | |
1e59de90 | 79 | long long scale = lltrunc(x, pol); |
92f5a8d4 TL |
80 | log_scaling += scale; |
81 | return exp(x - scale); | |
82 | } | |
83 | hypergeometric_1F1_igamma_series<T, Policy> s(b_minus_a, a - 1, x, pol); | |
84 | log_scaling += s.log_scaling; | |
1e59de90 | 85 | std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); |
92f5a8d4 TL |
86 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); |
87 | boost::math::policies::check_series_iterations<T>("boost::math::tgamma<%1%>(%1%,%1%)", max_iter, pol); | |
88 | T log_prefix = x + boost::math::lgamma(b, pol) - boost::math::lgamma(a, pol); | |
1e59de90 | 89 | long long scale = lltrunc(log_prefix); |
92f5a8d4 TL |
90 | log_scaling += scale; |
91 | return result * exp(log_prefix - scale); | |
92 | } | |
93 | ||
94 | template <class T, class Policy> | |
1e59de90 | 95 | T hypergeometric_1F1_shift_on_a(T h, const T& a_local, const T& b_local, const T& x, int a_shift, const Policy& pol, long long& log_scaling) |
92f5a8d4 TL |
96 | { |
97 | BOOST_MATH_STD_USING | |
98 | T a = a_local + a_shift; | |
99 | if (a_shift == 0) | |
100 | return h; | |
101 | else if (a_shift > 0) | |
102 | { | |
103 | // | |
104 | // Forward recursion on a is stable as long as 2a-b+z > 0. | |
105 | // If 2a-b+z < 0 then backwards recursion is stable even though | |
106 | // the function may be strictly increasing with a. Potentially | |
f67539c2 | 107 | // we may need to split the recurrence in 2 sections - one using |
92f5a8d4 TL |
108 | // forward recursion, and one backwards. |
109 | // | |
110 | // We will get the next seed value from the ratio | |
111 | // on b as that's stable and quick to compute. | |
112 | // | |
113 | ||
114 | T crossover_a = (b_local - x) / 2; | |
115 | int crossover_shift = itrunc(crossover_a - a_local); | |
116 | ||
117 | if (crossover_shift > 1) | |
118 | { | |
119 | // | |
120 | // Forwards recursion will start off unstable, but may switch to the stable direction later. | |
121 | // Start in the middle and go in both directions: | |
122 | // | |
123 | if (crossover_shift > a_shift) | |
124 | crossover_shift = a_shift; | |
125 | crossover_a = a_local + crossover_shift; | |
126 | boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(crossover_a, b_local, x); | |
1e59de90 | 127 | std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); |
92f5a8d4 TL |
128 | T b_ratio = boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter); |
129 | boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol); | |
130 | // | |
131 | // Convert to a ratio: | |
132 | // (1+a-b)M(a, b, z) - aM(a+1, b, z) + (b-1)M(a, b-1, z) = 0 | |
133 | // | |
134 | // hence: M(a+1,b,z) = ((1+a-b) / a) M(a,b,z) + ((b-1) / a) M(a,b,z)/b_ratio | |
135 | // | |
136 | T first = 1; | |
137 | T second = ((1 + crossover_a - b_local) / crossover_a) + ((b_local - 1) / crossover_a) / b_ratio; | |
138 | // | |
f67539c2 | 139 | // Recurse down to a_local, compare values and re-normalise first and second: |
92f5a8d4 TL |
140 | // |
141 | boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> a_coef(crossover_a, b_local, x); | |
1e59de90 | 142 | long long backwards_scale = 0; |
92f5a8d4 TL |
143 | T comparitor = boost::math::tools::apply_recurrence_relation_backward(a_coef, crossover_shift, second, first, &backwards_scale); |
144 | log_scaling -= backwards_scale; | |
145 | if ((h < 1) && (tools::max_value<T>() * h > comparitor)) | |
146 | { | |
147 | // Need to rescale! | |
1e59de90 | 148 | long long scale = lltrunc(log(h), pol) + 1; |
92f5a8d4 TL |
149 | h *= exp(T(-scale)); |
150 | log_scaling += scale; | |
151 | } | |
152 | comparitor /= h; | |
153 | first /= comparitor; | |
154 | second /= comparitor; | |
155 | // | |
156 | // Now we can recurse forwards for the rest of the range: | |
157 | // | |
158 | if (crossover_shift < a_shift) | |
159 | { | |
160 | boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> a_coef_2(crossover_a + 1, b_local, x); | |
161 | h = boost::math::tools::apply_recurrence_relation_forward(a_coef_2, a_shift - crossover_shift - 1, first, second, &log_scaling); | |
162 | } | |
163 | else | |
164 | h = first; | |
165 | } | |
166 | else | |
167 | { | |
168 | // | |
169 | // Regular case where forwards iteration is stable right from the start: | |
170 | // | |
171 | boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(a_local, b_local, x); | |
1e59de90 | 172 | std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); |
92f5a8d4 TL |
173 | T b_ratio = boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter); |
174 | boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol); | |
175 | // | |
176 | // Convert to a ratio: | |
177 | // (1+a-b)M(a, b, z) - aM(a+1, b, z) + (b-1)M(a, b-1, z) = 0 | |
178 | // | |
179 | // hence: M(a+1,b,z) = ((1+a-b) / a) M(a,b,z) + ((b-1) / a) M(a,b,z)/b_ratio | |
180 | // | |
181 | T second = ((1 + a_local - b_local) / a_local) * h + ((b_local - 1) / a_local) * h / b_ratio; | |
182 | boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> a_coef(a_local + 1, b_local, x); | |
183 | h = boost::math::tools::apply_recurrence_relation_forward(a_coef, --a_shift, h, second, &log_scaling); | |
184 | } | |
185 | } | |
186 | else | |
187 | { | |
188 | // | |
189 | // We've calculated h for a larger value of a than we want, and need to recurse down. | |
190 | // However, only forward iteration is stable, so calculate the ratio, compare values, | |
191 | // and normalise. Note that we calculate the ratio on b and convert to a since the | |
192 | // direction is the minimal solution for N->+INF. | |
193 | // | |
194 | // IMPORTANT: this is only currently called for a > b and therefore forwards iteration | |
195 | // is the only stable direction as we will only iterate down until a ~ b, but we | |
196 | // will check this with an assert: | |
197 | // | |
1e59de90 | 198 | BOOST_MATH_ASSERT(2 * a - b_local + x > 0); |
92f5a8d4 | 199 | boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(a, b_local, x); |
1e59de90 | 200 | std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); |
92f5a8d4 TL |
201 | T b_ratio = boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter); |
202 | boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol); | |
203 | // | |
204 | // Convert to a ratio: | |
205 | // (1+a-b)M(a, b, z) - aM(a+1, b, z) + (b-1)M(a, b-1, z) = 0 | |
206 | // | |
207 | // hence: M(a+1,b,z) = (1+a-b) / a M(a,b,z) + (b-1) / a M(a,b,z)/ (M(a,b,z)/M(a,b-1,z)) | |
208 | // | |
209 | T first = 1; // arbitrary value; | |
210 | T second = ((1 + a - b_local) / a) + ((b_local - 1) / a) * (1 / b_ratio); | |
211 | ||
212 | if (a_shift == -1) | |
213 | h = h / second; | |
214 | else | |
215 | { | |
216 | boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> a_coef(a + 1, b_local, x); | |
217 | T comparitor = boost::math::tools::apply_recurrence_relation_forward(a_coef, -(a_shift + 1), first, second); | |
218 | if (boost::math::tools::min_value<T>() * comparitor > h) | |
219 | { | |
220 | // Ooops, need to rescale h: | |
1e59de90 | 221 | long long rescale = lltrunc(log(fabs(h))); |
92f5a8d4 TL |
222 | T scale = exp(T(-rescale)); |
223 | h *= scale; | |
224 | log_scaling += rescale; | |
225 | } | |
226 | h = h / comparitor; | |
227 | } | |
228 | } | |
229 | return h; | |
230 | } | |
231 | ||
232 | template <class T, class Policy> | |
1e59de90 | 233 | T hypergeometric_1F1_shift_on_b(T h, const T& a, const T& b_local, const T& x, int b_shift, const Policy& pol, long long& log_scaling) |
92f5a8d4 TL |
234 | { |
235 | BOOST_MATH_STD_USING | |
236 | ||
237 | T b = b_local + b_shift; | |
238 | if (b_shift == 0) | |
239 | return h; | |
240 | else if (b_shift > 0) | |
241 | { | |
242 | // | |
243 | // We get here for b_shift > 0 when b > z. We can't use forward recursion on b - it's unstable, | |
244 | // so grab the ratio and work backwards to b - b_shift and normalise. | |
245 | // | |
246 | boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(a, b, x); | |
1e59de90 | 247 | std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); |
92f5a8d4 TL |
248 | |
249 | T first = 1; // arbitrary value; | |
250 | T second = 1 / boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter); | |
251 | boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol); | |
252 | if (b_shift == 1) | |
253 | h = h / second; | |
254 | else | |
255 | { | |
256 | // | |
257 | // Reset coefficients and recurse: | |
258 | // | |
259 | boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef_2(a, b - 1, x); | |
1e59de90 | 260 | long long local_scale = 0; |
92f5a8d4 TL |
261 | T comparitor = boost::math::tools::apply_recurrence_relation_backward(b_coef_2, --b_shift, first, second, &local_scale); |
262 | log_scaling -= local_scale; | |
263 | if (boost::math::tools::min_value<T>() * comparitor > h) | |
264 | { | |
265 | // Ooops, need to rescale h: | |
1e59de90 | 266 | long long rescale = lltrunc(log(fabs(h))); |
92f5a8d4 TL |
267 | T scale = exp(T(-rescale)); |
268 | h *= scale; | |
269 | log_scaling += rescale; | |
270 | } | |
271 | h = h / comparitor; | |
272 | } | |
273 | } | |
274 | else | |
275 | { | |
276 | T second; | |
277 | if (a == b_local) | |
278 | { | |
279 | // recurrence is trivial for a == b and method of ratios fails as the c-term goes to zero: | |
280 | second = -b_local * (1 - b_local - x) * h / (b_local * (b_local - 1)); | |
281 | } | |
282 | else | |
283 | { | |
1e59de90 | 284 | BOOST_MATH_ASSERT(!is_negative_integer(b - a)); |
92f5a8d4 | 285 | boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(a, b_local, x); |
1e59de90 | 286 | std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); |
92f5a8d4 TL |
287 | second = h / boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter); |
288 | boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol); | |
289 | } | |
290 | if (b_shift == -1) | |
291 | h = second; | |
292 | else | |
293 | { | |
294 | boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef_2(a, b_local - 1, x); | |
295 | h = boost::math::tools::apply_recurrence_relation_backward(b_coef_2, -(++b_shift), h, second, &log_scaling); | |
296 | } | |
297 | } | |
298 | return h; | |
299 | } | |
300 | ||
301 | ||
302 | template <class T, class Policy> | |
1e59de90 | 303 | T hypergeometric_1F1_large_igamma(const T& a, const T& b, const T& x, const T& b_minus_a, const Policy& pol, long long& log_scaling) |
92f5a8d4 TL |
304 | { |
305 | BOOST_MATH_STD_USING | |
306 | // | |
307 | // We need a < b < z in order to ensure there's at least a chance of convergence, | |
308 | // we can use recurrence relations to shift forwards on a+b or just a to achieve this, | |
309 | // for decent accuracy, try to keep 2b - 1 < a < 2b < z | |
310 | // | |
311 | int b_shift = b * 2 < x ? 0 : itrunc(b - x / 2); | |
312 | int a_shift = a > b - b_shift ? -itrunc(b - b_shift - a - 1) : -itrunc(b - b_shift - a); | |
313 | ||
314 | if (a_shift < 0) | |
315 | { | |
316 | // Might as well do all the shifting on b as scale a downwards: | |
317 | b_shift -= a_shift; | |
318 | a_shift = 0; | |
319 | } | |
320 | ||
321 | T a_local = a - a_shift; | |
322 | T b_local = b - b_shift; | |
323 | T b_minus_a_local = (a_shift == 0) && (b_shift == 0) ? b_minus_a : b_local - a_local; | |
1e59de90 | 324 | long long local_scaling = 0; |
92f5a8d4 TL |
325 | T h = hypergeometric_1F1_igamma(a_local, b_local, x, b_minus_a_local, pol, local_scaling); |
326 | log_scaling += local_scaling; | |
327 | ||
328 | // | |
329 | // Apply shifts on a and b as required: | |
330 | // | |
331 | h = hypergeometric_1F1_shift_on_a(h, a_local, b_local, x, a_shift, pol, log_scaling); | |
332 | h = hypergeometric_1F1_shift_on_b(h, a, b_local, x, b_shift, pol, log_scaling); | |
333 | ||
334 | return h; | |
335 | } | |
336 | ||
337 | template <class T, class Policy> | |
1e59de90 | 338 | T hypergeometric_1F1_large_series(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling) |
92f5a8d4 TL |
339 | { |
340 | BOOST_MATH_STD_USING | |
341 | // | |
342 | // We make a small, and b > z: | |
343 | // | |
344 | int a_shift(0), b_shift(0); | |
345 | if (a * z > b) | |
346 | { | |
347 | a_shift = itrunc(a) - 5; | |
348 | b_shift = b < z ? itrunc(b - z - 1) : 0; | |
349 | } | |
350 | // | |
351 | // If a_shift is trivially small, there's really not much point in losing | |
352 | // accuracy to save a couple of iterations: | |
353 | // | |
354 | if (a_shift < 5) | |
355 | a_shift = 0; | |
356 | T a_local = a - a_shift; | |
357 | T b_local = b - b_shift; | |
358 | T h = boost::math::detail::hypergeometric_1F1_generic_series(a_local, b_local, z, pol, log_scaling, "hypergeometric_1F1_large_series<%1%>(a,b,z)"); | |
359 | // | |
360 | // Apply shifts on a and b as required: | |
361 | // | |
362 | if (a_shift && (a_local == 0)) | |
363 | { | |
364 | // | |
365 | // Shifting on a via method of ratios in hypergeometric_1F1_shift_on_a fails when | |
366 | // a_local == 0. However, the value of h calculated was trivial (unity), so | |
367 | // calculate a second 1F1 for a == 1 and recurse as normal: | |
368 | // | |
1e59de90 | 369 | long long scale = 0; |
92f5a8d4 TL |
370 | T h2 = boost::math::detail::hypergeometric_1F1_generic_series(T(a_local + 1), b_local, z, pol, scale, "hypergeometric_1F1_large_series<%1%>(a,b,z)"); |
371 | if (scale != log_scaling) | |
372 | { | |
373 | h2 *= exp(T(scale - log_scaling)); | |
374 | } | |
375 | boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> coef(a_local + 1, b_local, z); | |
376 | h = boost::math::tools::apply_recurrence_relation_forward(coef, a_shift - 1, h, h2, &log_scaling); | |
377 | h = hypergeometric_1F1_shift_on_b(h, a, b_local, z, b_shift, pol, log_scaling); | |
378 | } | |
379 | else | |
380 | { | |
381 | h = hypergeometric_1F1_shift_on_a(h, a_local, b_local, z, a_shift, pol, log_scaling); | |
382 | h = hypergeometric_1F1_shift_on_b(h, a, b_local, z, b_shift, pol, log_scaling); | |
383 | } | |
384 | return h; | |
385 | } | |
386 | ||
387 | template <class T, class Policy> | |
1e59de90 | 388 | T hypergeometric_1F1_large_13_3_6_series(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling) |
92f5a8d4 TL |
389 | { |
390 | BOOST_MATH_STD_USING | |
391 | // | |
392 | // A&S 13.3.6 is good only when a ~ b, but isn't too fussy on the size of z. | |
393 | // So shift b to match a (b shifting seems to be more stable via method of ratios). | |
394 | // | |
395 | int b_shift = itrunc(b - a); | |
396 | T b_local = b - b_shift; | |
397 | T h = boost::math::detail::hypergeometric_1F1_AS_13_3_6(a, b_local, z, T(b_local - a), pol, log_scaling); | |
398 | return hypergeometric_1F1_shift_on_b(h, a, b_local, z, b_shift, pol, log_scaling); | |
399 | } | |
400 | ||
401 | template <class T, class Policy> | |
1e59de90 | 402 | T hypergeometric_1F1_large_abz(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling) |
92f5a8d4 TL |
403 | { |
404 | BOOST_MATH_STD_USING | |
405 | // | |
406 | // This is the selection logic to pick the "best" method. | |
f67539c2 | 407 | // We have a,b,z >> 0 and need to compute the approximate cost of each method |
92f5a8d4 TL |
408 | // and then select whichever wins out. |
409 | // | |
410 | enum method | |
411 | { | |
412 | method_series = 0, | |
413 | method_shifted_series, | |
414 | method_gamma, | |
415 | method_bessel | |
416 | }; | |
417 | // | |
418 | // Cost of direct series, is the approx number of terms required until we hit convergence: | |
419 | // | |
420 | T current_cost = (sqrt(16 * z * (3 * a + z) + 9 * b * b - 24 * b * z) - 3 * b + 4 * z) / 6; | |
421 | method current_method = method_series; | |
422 | // | |
423 | // Cost of shifted series, is the number of recurrences required to move to a zone where | |
424 | // the series is convergent right from the start. | |
425 | // Note that recurrence relations fail for very small b, and too many recurrences on a | |
426 | // will completely destroy all our digits. | |
427 | // Also note that the method fails when b-a is a negative integer unless b is already | |
428 | // larger than z and thus does not need shifting. | |
429 | // | |
430 | T cost = a + ((b < z) ? T(z - b) : T(0)); | |
431 | if((b > 1) && (cost < current_cost) && ((b > z) || !is_negative_integer(b-a))) | |
432 | { | |
433 | current_method = method_shifted_series; | |
434 | current_cost = cost; | |
435 | } | |
436 | // | |
437 | // Cost for gamma function method is the number of recurrences required to move it | |
438 | // into a convergent zone, note that recurrence relations fail for very small b. | |
439 | // Also add on a fudge factor to account for the fact that this method is both | |
440 | // more expensive to compute (requires gamma functions), and less accurate than the | |
441 | // methods above: | |
442 | // | |
443 | T b_shift = fabs(b * 2 < z ? T(0) : T(b - z / 2)); | |
444 | T a_shift = fabs(a > b - b_shift ? T(-(b - b_shift - a - 1)) : T(-(b - b_shift - a))); | |
445 | cost = 1000 + b_shift + a_shift; | |
446 | if((b > 1) && (cost <= current_cost)) | |
447 | { | |
448 | current_method = method_gamma; | |
449 | current_cost = cost; | |
450 | } | |
451 | // | |
452 | // Cost for bessel approximation is the number of recurrences required to make a ~ b, | |
453 | // Note that recurrence relations fail for very small b. We also have issue with large | |
454 | // z: either overflow/numeric instability or else the series goes divergent. We seem to be | |
455 | // OK for z smaller than log_max_value<Quad> though, maybe we can stretch a little further | |
456 | // but that's not clear... | |
457 | // Also need to add on a fudge factor to the cost to account for the fact that we need | |
458 | // to calculate the Bessel functions, this is not quite as high as the gamma function | |
f67539c2 | 459 | // method above as this is generally more accurate and so preferred if the methods are close: |
92f5a8d4 TL |
460 | // |
461 | cost = 50 + fabs(b - a); | |
462 | if((b > 1) && (cost <= current_cost) && (z < tools::log_max_value<T>()) && (z < 11356) && (b - a != 0.5f)) | |
463 | { | |
464 | current_method = method_bessel; | |
465 | current_cost = cost; | |
466 | } | |
467 | ||
468 | switch (current_method) | |
469 | { | |
470 | case method_series: | |
471 | return detail::hypergeometric_1F1_generic_series(a, b, z, pol, log_scaling, "hypergeometric_1f1_large_abz<%1%>(a,b,z)"); | |
472 | case method_shifted_series: | |
473 | return detail::hypergeometric_1F1_large_series(a, b, z, pol, log_scaling); | |
474 | case method_gamma: | |
475 | return detail::hypergeometric_1F1_large_igamma(a, b, z, T(b - a), pol, log_scaling); | |
476 | case method_bessel: | |
477 | return detail::hypergeometric_1F1_large_13_3_6_series(a, b, z, pol, log_scaling); | |
478 | } | |
479 | return 0; // We don't get here! | |
480 | } | |
481 | ||
482 | } } } // namespaces | |
483 | ||
484 | #endif // BOOST_HYPERGEOMETRIC_1F1_LARGE_ABZ_HPP_ |