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_1F1_HPP
11 #define BOOST_MATH_HYPERGEOMETRIC_1F1_HPP
12
13 #include <boost/math/tools/config.hpp>
14 #include <boost/math/policies/policy.hpp>
15 #include <boost/math/policies/error_handling.hpp>
16 #include <boost/math/special_functions/detail/hypergeometric_series.hpp>
17 #include <boost/math/special_functions/detail/hypergeometric_asym.hpp>
18 #include <boost/math/special_functions/detail/hypergeometric_rational.hpp>
19 #include <boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp>
20 #include <boost/math/special_functions/detail/hypergeometric_1F1_by_ratios.hpp>
21 #include <boost/math/special_functions/detail/hypergeometric_pade.hpp>
22 #include <boost/math/special_functions/detail/hypergeometric_1F1_bessel.hpp>
23 #include <boost/math/special_functions/detail/hypergeometric_1F1_scaled_series.hpp>
24 #include <boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp>
25 #include <boost/math/special_functions/detail/hypergeometric_1F1_addition_theorems_on_z.hpp>
26 #include <boost/math/special_functions/detail/hypergeometric_1F1_large_abz.hpp>
27 #include <boost/math/special_functions/detail/hypergeometric_1F1_small_a_negative_b_by_r.hpp>
28 #include <boost/math/special_functions/detail/hypergeometric_1F1_negative_b_regions.hpp>
29
30 namespace boost { namespace math { namespace detail {
31
32 // check when 1F1 series can't decay to polynom
33 template <class T>
check_hypergeometric_1F1_parameters(const T & a,const T & b)34 inline bool check_hypergeometric_1F1_parameters(const T& a, const T& b)
35 {
36 BOOST_MATH_STD_USING
37
38 if ((b <= 0) && (b == floor(b)))
39 {
40 if ((a >= 0) || (a < b) || (a != floor(a)))
41 return false;
42 }
43
44 return true;
45 }
46
47 template <class T, class Policy>
48 T hypergeometric_1F1_divergent_fallback(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling)
49 {
50 BOOST_MATH_STD_USING
51 const char* function = "hypergeometric_1F1_divergent_fallback<%1%>(%1%,%1%,%1%)";
52 //
53 // We get here if either:
54 // 1) We decide up front that Tricomi's method won't work, or:
55 // 2) We've called Tricomi's method and it's failed.
56 //
57 if (b > 0)
58 {
59 // Commented out since recurrence seems to always be better?
60 #if 0
61 if ((z < b) && (a > -50))
62 // Might as well use a recurrence in preference to z-recurrence:
63 return hypergeometric_1F1_backward_recurrence_for_negative_a(a, b, z, pol, function, log_scaling);
64 T z_limit = fabs((2 * a - b) / (sqrt(fabs(a))));
65 int k = 1 + itrunc(z - z_limit);
66 // If k is too large we destroy all the digits in the result:
67 T convergence_at_50 = (b - a + 50) * k / (z * 50);
68 if ((k > 0) && (k < 50) && (fabs(convergence_at_50) < 1) && (z > z_limit))
69 {
70 return boost::math::detail::hypergeometric_1f1_recurrence_on_z_minus_zero(a, b, T(z - k), k, pol, log_scaling);
71 }
72 #endif
73 if (z < b)
74 return hypergeometric_1F1_backward_recurrence_for_negative_a(a, b, z, pol, function, log_scaling);
75 else
76 return hypergeometric_1F1_backwards_recursion_on_b_for_negative_a(a, b, z, pol, function, log_scaling);
77 }
78 else // b < 0
79 {
80 if (a < 0)
81 {
82 if ((b < a) && (z < -b / 4))
83 return hypergeometric_1F1_from_function_ratio_negative_ab(a, b, z, pol, log_scaling);
84 else
85 {
86 //
87 // Solve (a+n)z/((b+n)n) == 1 for n, the number of iterations till the series starts to converge.
88 // If this is well away from the origin then it's probably better to use the series to evaluate this.
89 // Note that if sqr is negative then we have no solution, so assign an arbitrarily large value to the
90 // number of iterations.
91 //
92 bool can_use_recursion = (z - b + 100 < boost::math::policies::get_max_series_iterations<Policy>()) && (100 - a < boost::math::policies::get_max_series_iterations<Policy>());
93 T sqr = 4 * a * z + b * b - 2 * b * z + z * z;
94 T iterations_to_convergence = sqr > 0 ? T(0.5f * (-sqrt(sqr) - b + z)) : T(-a - b);
95 if(can_use_recursion && ((std::max)(a, b) + iterations_to_convergence > -300))
96 return hypergeometric_1F1_backwards_recursion_on_b_for_negative_a(a, b, z, pol, function, log_scaling);
97 //
98 // When a < b and if we fall through to the series, then we get divergent behaviour when b crosses the origin
99 // so ideally we would pick another method. Otherwise the terms immediately after b crosses the origin may
100 // suffer catastrophic cancellation....
101 //
102 if((a < b) && can_use_recursion)
103 return hypergeometric_1F1_backwards_recursion_on_b_for_negative_a(a, b, z, pol, function, log_scaling);
104 }
105 }
106 else
107 {
108 //
109 // Start by getting the domain of the recurrence relations, we get either:
110 // -1 Backwards recursion is stable and the CF will converge to double precision.
111 // +1 Forwards recursion is stable and the CF will converge to double precision.
112 // 0 No man's land, we're not far enough away from the crossover point to get double precision from either CF.
113 //
114 // At higher than double precision we need to be further away from the crossover location to
115 // get full converge, but it's not clear how much further - indeed at quad precision it's
116 // basically impossible to ever get forwards iteration to work. Backwards seems to work
117 // OK as long as a > 1 whatever the precision tbough.
118 //
119 int domain = hypergeometric_1F1_negative_b_recurrence_region(a, b, z);
120 if ((domain < 0) && ((a > 1) || (boost::math::policies::digits<T, Policy>() <= 64)))
121 return hypergeometric_1F1_from_function_ratio_negative_b(a, b, z, pol, log_scaling);
122 else if (domain > 0)
123 {
124 if (boost::math::policies::digits<T, Policy>() <= 64)
125 return hypergeometric_1F1_from_function_ratio_negative_b_forwards(a, b, z, pol, log_scaling);
126 try
127 {
128 return hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling);
129 }
130 catch (const evaluation_error&)
131 {
132 //
133 // The series failed, try the recursions instead and hope we get at least double precision:
134 //
135 return hypergeometric_1F1_from_function_ratio_negative_b_forwards(a, b, z, pol, log_scaling);
136 }
137 }
138 //
139 // We could fall back to Tricomi's approximation if we're in the transition zone
140 // between the above two regions. However, I've been unable to find any examples
141 // where this is better than the series, and there are many cases where it leads to
142 // quite grievous errors.
143 /*
144 else if (allow_tricomi)
145 {
146 T aa = a < 1 ? T(1) : a;
147 if (z < fabs((2 * aa - b) / (sqrt(fabs(aa * b)))))
148 return hypergeometric_1F1_AS_13_3_7_tricomi(a, b, z, pol, log_scaling);
149 }
150 */
151 }
152 }
153
154 // If we get here, then we've run out of methods to try, use the checked series which will
155 // raise an error if the result is garbage:
156 return hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling);
157 }
158
159 template <class T>
is_convergent_negative_z_series(const T & a,const T & b,const T & z,const T & b_minus_a)160 bool is_convergent_negative_z_series(const T& a, const T& b, const T& z, const T& b_minus_a)
161 {
162 BOOST_MATH_STD_USING
163 //
164 // Filter out some cases we don't want first:
165 //
166 if((b_minus_a > 0) && (b > 0))
167 {
168 if (a < 0)
169 return false;
170 }
171 //
172 // Generic check: we have small initial divergence and are convergent after 10 terms:
173 //
174 if ((fabs(z * a / b) < 2) && (fabs(z * (a + 10) / ((b + 10) * 10)) < 1))
175 {
176 // Double check for divergence when we cross the origin on a and b:
177 if (a < 0)
178 {
179 T n = 300 - floor(a);
180 if (fabs((a + n) * z / ((b + n) * n)) < 1)
181 {
182 if (b < 0)
183 {
184 T m = 3 - floor(b);
185 if (fabs((a + m) * z / ((b + m) * m)) < 1)
186 return true;
187 }
188 else
189 return true;
190 }
191 }
192 else if (b < 0)
193 {
194 T n = 3 - floor(b);
195 if (fabs((a + n) * z / ((b + n) * n)) < 1)
196 return true;
197 }
198 }
199 if ((b > 0) && (a < 0))
200 {
201 //
202 // For a and z both negative, we're OK with some initial divergence as long as
203 // it occurs before we hit the origin, as to start with all the terms have the
204 // same sign.
205 //
206 // https://www.wolframalpha.com/input/?i=solve+(a%2Bn)z+%2F+((b%2Bn)n)+%3D%3D+1+for+n
207 //
208 T sqr = 4 * a * z + b * b - 2 * b * z + z * z;
209 T iterations_to_convergence = sqr > 0 ? T(0.5f * (-sqrt(sqr) - b + z)) : T(-a + b);
210 if (iterations_to_convergence < 0)
211 iterations_to_convergence = 0.5f * (sqrt(sqr) - b + z);
212 if (a + iterations_to_convergence < -50)
213 {
214 // Need to check for divergence when we cross the origin on a:
215 if (a > -1)
216 return true;
217 T n = 300 - floor(a);
218 if(fabs((a + n) * z / ((b + n) * n)) < 1)
219 return true;
220 }
221 }
222 return false;
223 }
224
225 template <class T>
cyl_bessel_i_shrinkage_rate(const T & z)226 inline T cyl_bessel_i_shrinkage_rate(const T& z)
227 {
228 // Approximately the ratio I_10.5(z/2) / I_9.5(z/2), this gives us an idea of how quickly
229 // the Bessel terms in A&S 13.6.4 are converging:
230 if (z < -160)
231 return 1;
232 if (z < -40)
233 return 0.75f;
234 if (z < -20)
235 return 0.5f;
236 if (z < -7)
237 return 0.25f;
238 if (z < -2)
239 return 0.1f;
240 return 0.05f;
241 }
242
243 template <class T>
hypergeometric_1F1_is_13_3_6_region(const T & a,const T & b,const T & z)244 inline bool hypergeometric_1F1_is_13_3_6_region(const T& a, const T& b, const T& z)
245 {
246 BOOST_MATH_STD_USING
247 if(fabs(a) == 0.5)
248 return false;
249 if ((z < 0) && (fabs(10 * a / b) < 1) && (fabs(a) < 50))
250 {
251 T shrinkage = cyl_bessel_i_shrinkage_rate(z);
252 // We want the first term not too divergent, and convergence by term 10:
253 if ((fabs((2 * a - 1) * (2 * a - b) / b) < 2) && (fabs(shrinkage * (2 * a + 9) * (2 * a - b + 10) / (10 * (b + 10))) < 0.75))
254 return true;
255 }
256 return false;
257 }
258
259 template <class T>
hypergeometric_1F1_need_kummer_reflection(const T & a,const T & b,const T & z)260 inline bool hypergeometric_1F1_need_kummer_reflection(const T& a, const T& b, const T& z)
261 {
262 BOOST_MATH_STD_USING
263 //
264 // Check to see if we should apply Kummer's relation or not:
265 //
266 if (z > 0)
267 return false;
268 if (z < -1)
269 return true;
270 //
271 // When z is small and negative, things get more complex.
272 // More often than not we do not need apply Kummer's relation and the
273 // series is convergent as is, but we do need to check:
274 //
275 if (a > 0)
276 {
277 if (b > 0)
278 {
279 return fabs((a + 10) * z / (10 * (b + 10))) < 1; // Is the 10'th term convergent?
280 }
281 else
282 {
283 return true; // Likely to be divergent as b crosses the origin
284 }
285 }
286 else // a < 0
287 {
288 if (b > 0)
289 {
290 return false; // Terms start off all positive and then by the time a crosses the origin we *must* be convergent.
291 }
292 else
293 {
294 return true; // Likely to be divergent as b crosses the origin, but hard to rationalise about!
295 }
296 }
297 }
298
299
300 template <class T, class Policy>
301 T hypergeometric_1F1_imp(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling)
302 {
303 BOOST_MATH_STD_USING // exp, fabs, sqrt
304
305 static const char* const function = "boost::math::hypergeometric_1F1<%1%,%1%,%1%>(%1%,%1%,%1%)";
306
307 if ((z == 0) || (a == 0))
308 return T(1);
309
310 // undefined result:
311 if (!detail::check_hypergeometric_1F1_parameters(a, b))
312 return policies::raise_domain_error<T>(
313 function,
314 "Function is indeterminate for negative integer b = %1%.",
315 b,
316 pol);
317
318 // other checks:
319 if (a == -1)
320 return 1 - (z / b);
321
322 const T b_minus_a = b - a;
323
324 // 0f0 a == b case;
325 if (b_minus_a == 0)
326 {
327 long long scale = lltrunc(z, pol);
328 log_scaling += scale;
329 return exp(z - scale);
330 }
331 // Special case for b-a = -1, we don't use for small a as it throws the digits of a away and leads to large errors:
332 if ((b_minus_a == -1) && (fabs(a) > 0.5))
333 {
334 // for negative small integer a it is reasonable to use truncated series - polynomial
335 if ((a < 0) && (a == ceil(a)) && (a > -50))
336 return detail::hypergeometric_1F1_generic_series(a, b, z, pol, log_scaling, function);
337
338 return (b + z) * exp(z) / b;
339 }
340
341 if ((a == 1) && (b == 2))
342 return boost::math::expm1(z, pol) / z;
343
344 if ((b - a == b) && (fabs(z / b) < policies::get_epsilon<T, Policy>()))
345 return 1;
346 //
347 // Special case for A&S 13.3.6:
348 //
349 if (z < 0)
350 {
351 if (hypergeometric_1F1_is_13_3_6_region(a, b, z))
352 {
353 // a is tiny compared to b, and z < 0
354 // 13.3.6 appears to be the most efficient and often the most accurate method.
355 T r = boost::math::detail::hypergeometric_1F1_AS_13_3_6(b_minus_a, b, T(-z), a, pol, log_scaling);
356 long long scale = lltrunc(z, pol);
357 log_scaling += scale;
358 return r * exp(z - scale);
359 }
360 if ((b < 0) && (fabs(a) < 1e-2))
361 {
362 //
363 // This is a tricky area, potentially we have no good method at all:
364 //
365 if (b - ceil(b) == a)
366 {
367 // Fractional parts of a and b are genuinely equal, we might as well
368 // apply Kummer's relation and get a truncated series:
369 long long scaling = lltrunc(z);
370 T r = exp(z - scaling) * detail::hypergeometric_1F1_imp<T>(b_minus_a, b, -z, pol, log_scaling);
371 log_scaling += scaling;
372 return r;
373 }
374 if ((b < -1) && (max_b_for_1F1_small_a_negative_b_by_ratio(z) < b))
375 return hypergeometric_1F1_small_a_negative_b_by_ratio(a, b, z, pol, log_scaling);
376 if ((b > -1) && (b < -0.5f))
377 {
378 // Recursion is meta-stable:
379 T first = hypergeometric_1F1_imp(a, T(b + 2), z, pol);
380 T second = hypergeometric_1F1_imp(a, T(b + 1), z, pol);
381 return tools::apply_recurrence_relation_backward(hypergeometric_1F1_recurrence_small_b_coefficients<T>(a, b, z, 1), 1, first, second);
382 }
383 //
384 // We've got nothing left but 13.3.6, even though it may be initially divergent:
385 //
386 T r = boost::math::detail::hypergeometric_1F1_AS_13_3_6(b_minus_a, b, T(-z), a, pol, log_scaling);
387 long long scale = lltrunc(z, pol);
388 log_scaling += scale;
389 return r * exp(z - scale);
390 }
391 }
392 //
393 // Asymptotic expansion for large z
394 // TODO: check region for higher precision types.
395 // Use recurrence relations to move to this region when a and b are also large.
396 //
397 if (detail::hypergeometric_1F1_asym_region(a, b, z, pol))
398 {
399 long long saved_scale = log_scaling;
400 try
401 {
402 return hypergeometric_1F1_asym_large_z_series(a, b, z, pol, log_scaling);
403 }
404 catch (const evaluation_error&)
405 {
406 }
407 //
408 // Very occasionally our convergence criteria don't quite go to full precision
409 // and we have to try another method:
410 //
411 log_scaling = saved_scale;
412 }
413
414 if ((fabs(a * z / b) < 3.5) && (fabs(z * 100) < fabs(b)) && ((fabs(a) > 1e-2) || (b < -5)))
415 return detail::hypergeometric_1F1_rational(a, b, z, pol);
416
417 if (hypergeometric_1F1_need_kummer_reflection(a, b, z))
418 {
419 if (a == 1)
420 return detail::hypergeometric_1F1_pade(b, z, pol);
421 if (is_convergent_negative_z_series(a, b, z, b_minus_a))
422 {
423 if ((boost::math::sign(b_minus_a) == boost::math::sign(b)) && ((b > 0) || (b < -200)))
424 {
425 // Series is close enough to convergent that we should be OK,
426 // In this domain b - a ~ b and since 1F1[a, a, z] = e^z 1F1[b-a, b, -z]
427 // and 1F1[a, a, -z] = e^-z the result must necessarily be somewhere near unity.
428 // We have to rule out b small and negative because if b crosses the origin early
429 // in the series (before we're pretty much converged) then all bets are off.
430 // Note that this can go badly wrong when b and z are both large and negative,
431 // in that situation the series goes in waves of large and small values which
432 // may or may not cancel out. Likewise the initial part of the series may or may
433 // not converge, and even if it does may or may not give a correct answer!
434 // For example 1F1[-small, -1252.5, -1043.7] can loose up to ~800 digits due to
435 // cancellation and is basically incalculable via this method.
436 return hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling);
437 }
438 }
439 // Let's otherwise make z positive (almost always)
440 // by Kummer's transformation
441 // (we also don't transform if z belongs to [-1,0])
442 long long scaling = lltrunc(z);
443 T r = exp(z - scaling) * detail::hypergeometric_1F1_imp<T>(b_minus_a, b, -z, pol, log_scaling);
444 log_scaling += scaling;
445 return r;
446 }
447 //
448 // Check for initial divergence:
449 //
450 bool series_is_divergent = (a + 1) * z / (b + 1) < -1;
451 if (series_is_divergent && (a < 0) && (b < 0) && (a > -1))
452 series_is_divergent = false; // Best off taking the series in this situation
453 //
454 // If series starts off non-divergent, and becomes divergent later
455 // then it's because both a and b are negative, so check for later
456 // divergence as well:
457 //
458 if (!series_is_divergent && (a < 0) && (b < 0) && (b > a))
459 {
460 //
461 // We need to exclude situations where we're over the initial "hump"
462 // in the series terms (ie series has already converged by the time
463 // b crosses the origin:
464 //
465 //T fa = fabs(a);
466 //T fb = fabs(b);
467 T convergence_point = sqrt((a - 1) * (a - b)) - a;
468 if (-b < convergence_point)
469 {
470 T n = -floor(b);
471 series_is_divergent = (a + n) * z / ((b + n) * n) < -1;
472 }
473 }
474 else if (!series_is_divergent && (b < 0) && (a > 0))
475 {
476 // Series almost always become divergent as b crosses the origin:
477 series_is_divergent = true;
478 }
479 if (series_is_divergent && (b < -1) && (b > -5) && (a > b))
480 series_is_divergent = false; // don't bother with divergence, series will be OK
481
482 //
483 // Test for alternating series due to negative a,
484 // in particular, see if the series is initially divergent
485 // If so use the recurrence relation on a:
486 //
487 if (series_is_divergent)
488 {
489 if((a < 0) && (floor(a) == a) && (-a < policies::get_max_series_iterations<Policy>()))
490 // This works amazingly well for negative integer a:
491 return hypergeometric_1F1_backward_recurrence_for_negative_a(a, b, z, pol, function, log_scaling);
492 //
493 // In what follows we have to set limits on how large z can be otherwise
494 // the Bessel series become large and divergent and all the digits cancel out.
495 // The criteria are distinctly empiracle rather than based on a firm analysis
496 // of the terms in the series.
497 //
498 if (b > 0)
499 {
500 T z_limit = fabs((2 * a - b) / (sqrt(fabs(a))));
501 if ((z < z_limit) && hypergeometric_1F1_is_tricomi_viable_positive_b(a, b, z))
502 return detail::hypergeometric_1F1_AS_13_3_7_tricomi(a, b, z, pol, log_scaling);
503 }
504 else // b < 0
505 {
506 if (a < 0)
507 {
508 T z_limit = fabs((2 * a - b) / (sqrt(fabs(a))));
509 //
510 // I hate these hard limits, but they're about the best we can do to try and avoid
511 // Bessel function internal failures: these will be caught and handled
512 // but up the expense of this function call:
513 //
514 if (((z < z_limit) || (a > -500)) && ((b > -500) || (b - 2 * a > 0)) && (z < -a))
515 {
516 //
517 // Outside this domain we will probably get better accuracy from the recursive methods.
518 //
519 if(!(((a < b) && (z > -b)) || (z > z_limit)))
520 return detail::hypergeometric_1F1_AS_13_3_7_tricomi(a, b, z, pol, log_scaling);
521 //
522 // When b and z are both very small, we get large errors from the recurrence methods
523 // in the fallbacks. Tricomi seems to work well here, as does direct series evaluation
524 // at least some of the time. Picking the right method is not easy, and sometimes this
525 // is much worse than the fallback. Overall though, it's a reasonable choice that keeps
526 // the very worst errors under control.
527 //
528 if(b > -1)
529 return detail::hypergeometric_1F1_AS_13_3_7_tricomi(a, b, z, pol, log_scaling);
530 }
531 }
532 //
533 // We previously used Tricomi here, but it appears to be worse than
534 // the recurrence-based algorithms in hypergeometric_1F1_divergent_fallback.
535 /*
536 else
537 {
538 T aa = a < 1 ? T(1) : a;
539 if (z < fabs((2 * aa - b) / (sqrt(fabs(aa * b)))))
540 return detail::hypergeometric_1F1_AS_13_3_7_tricomi(a, b, z, pol, log_scaling);
541 }*/
542 }
543
544 return hypergeometric_1F1_divergent_fallback(a, b, z, pol, log_scaling);
545 }
546
547 if (hypergeometric_1F1_is_13_3_6_region(b_minus_a, b, T(-z)))
548 {
549 // b_minus_a is tiny compared to b, and -z < 0
550 // 13.3.6 appears to be the most efficient and often the most accurate method.
551 return boost::math::detail::hypergeometric_1F1_AS_13_3_6(a, b, z, b_minus_a, pol, log_scaling);
552 }
553 #if 0
554 if ((a > 0) && (b > 0) && (a * z / b > 2))
555 {
556 //
557 // Series is initially divergent and slow to converge, see if applying
558 // Kummer's relation can improve things:
559 //
560 if (is_convergent_negative_z_series(b_minus_a, b, T(-z), b_minus_a))
561 {
562 long long scaling = lltrunc(z);
563 T r = exp(z - scaling) * detail::hypergeometric_1F1_checked_series_impl(b_minus_a, b, T(-z), pol, log_scaling);
564 log_scaling += scaling;
565 return r;
566 }
567
568 }
569 #endif
570 if ((a > 0) && (b > 0) && (a * z > 50))
571 return detail::hypergeometric_1F1_large_abz(a, b, z, pol, log_scaling);
572
573 if (b < 0)
574 return detail::hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling);
575
576 return detail::hypergeometric_1F1_generic_series(a, b, z, pol, log_scaling, function);
577 }
578
579 template <class T, class Policy>
hypergeometric_1F1_imp(const T & a,const T & b,const T & z,const Policy & pol)580 inline T hypergeometric_1F1_imp(const T& a, const T& b, const T& z, const Policy& pol)
581 {
582 BOOST_MATH_STD_USING // exp, fabs, sqrt
583 long long log_scaling = 0;
584 T result = hypergeometric_1F1_imp(a, b, z, pol, log_scaling);
585 //
586 // Actual result will be result * e^log_scaling.
587 //
588 static const thread_local long long max_scaling = lltrunc(boost::math::tools::log_max_value<T>()) - 2;
589 static const thread_local T max_scale_factor = exp(T(max_scaling));
590
591 while (log_scaling > max_scaling)
592 {
593 result *= max_scale_factor;
594 log_scaling -= max_scaling;
595 }
596 while (log_scaling < -max_scaling)
597 {
598 result /= max_scale_factor;
599 log_scaling += max_scaling;
600 }
601 if (log_scaling)
602 result *= exp(T(log_scaling));
603 return result;
604 }
605
606 template <class T, class Policy>
log_hypergeometric_1F1_imp(const T & a,const T & b,const T & z,int * sign,const Policy & pol)607 inline T log_hypergeometric_1F1_imp(const T& a, const T& b, const T& z, int* sign, const Policy& pol)
608 {
609 BOOST_MATH_STD_USING // exp, fabs, sqrt
610 long long log_scaling = 0;
611 T result = hypergeometric_1F1_imp(a, b, z, pol, log_scaling);
612 if (sign)
613 *sign = result < 0 ? -1 : 1;
614 result = log(fabs(result)) + log_scaling;
615 return result;
616 }
617
618 template <class T, class Policy>
hypergeometric_1F1_regularized_imp(const T & a,const T & b,const T & z,const Policy & pol)619 inline T hypergeometric_1F1_regularized_imp(const T& a, const T& b, const T& z, const Policy& pol)
620 {
621 BOOST_MATH_STD_USING // exp, fabs, sqrt
622 long long log_scaling = 0;
623 T result = hypergeometric_1F1_imp(a, b, z, pol, log_scaling);
624 //
625 // Actual result will be result * e^log_scaling / tgamma(b).
626 //
627 int result_sign = 1;
628 T scale = log_scaling - boost::math::lgamma(b, &result_sign, pol);
629
630 static const thread_local T max_scaling = boost::math::tools::log_max_value<T>() - 2;
631 static const thread_local T max_scale_factor = exp(max_scaling);
632
633 while (scale > max_scaling)
634 {
635 result *= max_scale_factor;
636 scale -= max_scaling;
637 }
638 while (scale < -max_scaling)
639 {
640 result /= max_scale_factor;
641 scale += max_scaling;
642 }
643 if (scale != 0)
644 result *= exp(scale);
645 return result * result_sign;
646 }
647
648 } // namespace detail
649
650 template <class T1, class T2, class T3, class Policy>
hypergeometric_1F1(T1 a,T2 b,T3 z,const Policy &)651 inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_1F1(T1 a, T2 b, T3 z, const Policy& /* pol */)
652 {
653 BOOST_FPU_EXCEPTION_GUARD
654 typedef typename tools::promote_args<T1, T2, T3>::type result_type;
655 typedef typename policies::evaluation<result_type, Policy>::type value_type;
656 typedef typename policies::normalise<
657 Policy,
658 policies::promote_float<false>,
659 policies::promote_double<false>,
660 policies::discrete_quantile<>,
661 policies::assert_undefined<> >::type forwarding_policy;
662 return policies::checked_narrowing_cast<result_type, Policy>(
663 detail::hypergeometric_1F1_imp<value_type>(
664 static_cast<value_type>(a),
665 static_cast<value_type>(b),
666 static_cast<value_type>(z),
667 forwarding_policy()),
668 "boost::math::hypergeometric_1F1<%1%>(%1%,%1%,%1%)");
669 }
670
671 template <class T1, class T2, class T3>
hypergeometric_1F1(T1 a,T2 b,T3 z)672 inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_1F1(T1 a, T2 b, T3 z)
673 {
674 return hypergeometric_1F1(a, b, z, policies::policy<>());
675 }
676
677 template <class T1, class T2, class T3, class Policy>
hypergeometric_1F1_regularized(T1 a,T2 b,T3 z,const Policy &)678 inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_1F1_regularized(T1 a, T2 b, T3 z, const Policy& /* pol */)
679 {
680 BOOST_FPU_EXCEPTION_GUARD
681 typedef typename tools::promote_args<T1, T2, T3>::type result_type;
682 typedef typename policies::evaluation<result_type, Policy>::type value_type;
683 typedef typename policies::normalise<
684 Policy,
685 policies::promote_float<false>,
686 policies::promote_double<false>,
687 policies::discrete_quantile<>,
688 policies::assert_undefined<> >::type forwarding_policy;
689 return policies::checked_narrowing_cast<result_type, Policy>(
690 detail::hypergeometric_1F1_regularized_imp<value_type>(
691 static_cast<value_type>(a),
692 static_cast<value_type>(b),
693 static_cast<value_type>(z),
694 forwarding_policy()),
695 "boost::math::hypergeometric_1F1<%1%>(%1%,%1%,%1%)");
696 }
697
698 template <class T1, class T2, class T3>
hypergeometric_1F1_regularized(T1 a,T2 b,T3 z)699 inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_1F1_regularized(T1 a, T2 b, T3 z)
700 {
701 return hypergeometric_1F1_regularized(a, b, z, policies::policy<>());
702 }
703
704 template <class T1, class T2, class T3, class Policy>
log_hypergeometric_1F1(T1 a,T2 b,T3 z,const Policy &)705 inline typename tools::promote_args<T1, T2, T3>::type log_hypergeometric_1F1(T1 a, T2 b, T3 z, const Policy& /* pol */)
706 {
707 BOOST_FPU_EXCEPTION_GUARD
708 typedef typename tools::promote_args<T1, T2, T3>::type result_type;
709 typedef typename policies::evaluation<result_type, Policy>::type value_type;
710 typedef typename policies::normalise<
711 Policy,
712 policies::promote_float<false>,
713 policies::promote_double<false>,
714 policies::discrete_quantile<>,
715 policies::assert_undefined<> >::type forwarding_policy;
716 return policies::checked_narrowing_cast<result_type, Policy>(
717 detail::log_hypergeometric_1F1_imp<value_type>(
718 static_cast<value_type>(a),
719 static_cast<value_type>(b),
720 static_cast<value_type>(z),
721 0,
722 forwarding_policy()),
723 "boost::math::hypergeometric_1F1<%1%>(%1%,%1%,%1%)");
724 }
725
726 template <class T1, class T2, class T3>
log_hypergeometric_1F1(T1 a,T2 b,T3 z)727 inline typename tools::promote_args<T1, T2, T3>::type log_hypergeometric_1F1(T1 a, T2 b, T3 z)
728 {
729 return log_hypergeometric_1F1(a, b, z, policies::policy<>());
730 }
731
732 template <class T1, class T2, class T3, class Policy>
log_hypergeometric_1F1(T1 a,T2 b,T3 z,int * sign,const Policy &)733 inline typename tools::promote_args<T1, T2, T3>::type log_hypergeometric_1F1(T1 a, T2 b, T3 z, int* sign, const Policy& /* pol */)
734 {
735 BOOST_FPU_EXCEPTION_GUARD
736 typedef typename tools::promote_args<T1, T2, T3>::type result_type;
737 typedef typename policies::evaluation<result_type, Policy>::type value_type;
738 typedef typename policies::normalise<
739 Policy,
740 policies::promote_float<false>,
741 policies::promote_double<false>,
742 policies::discrete_quantile<>,
743 policies::assert_undefined<> >::type forwarding_policy;
744 return policies::checked_narrowing_cast<result_type, Policy>(
745 detail::log_hypergeometric_1F1_imp<value_type>(
746 static_cast<value_type>(a),
747 static_cast<value_type>(b),
748 static_cast<value_type>(z),
749 sign,
750 forwarding_policy()),
751 "boost::math::hypergeometric_1F1<%1%>(%1%,%1%,%1%)");
752 }
753
754 template <class T1, class T2, class T3>
log_hypergeometric_1F1(T1 a,T2 b,T3 z,int * sign)755 inline typename tools::promote_args<T1, T2, T3>::type log_hypergeometric_1F1(T1 a, T2 b, T3 z, int* sign)
756 {
757 return log_hypergeometric_1F1(a, b, z, sign, policies::policy<>());
758 }
759
760
761 } } // namespace boost::math
762
763 #endif // BOOST_MATH_HYPERGEOMETRIC_HPP
764