From fb52d2de42ff7e1cbc61103a25fbc82c5e034a90 Mon Sep 17 00:00:00 2001 From: John Maddock Date: Thu, 24 Oct 2013 08:35:41 +0000 Subject: [PATCH] Fix initialization of power series so that we don't get a spurious overflow error from tgamma when the result is actually zero. [SVN r86415] --- include/boost/math/special_functions/bessel.hpp | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/include/boost/math/special_functions/bessel.hpp b/include/boost/math/special_functions/bessel.hpp index b1e03bfbe..ceca30916 100644 --- a/include/boost/math/special_functions/bessel.hpp +++ b/include/boost/math/special_functions/bessel.hpp @@ -48,7 +48,13 @@ struct sph_bessel_j_small_z_series_term { BOOST_MATH_STD_USING mult = x / 2; - term = pow(mult, T(v)) / boost::math::tgamma(v+1+T(0.5f), Policy()); + if(v + 3 > max_factorial::value) + { + term = v * log(mult) - boost::math::lgamma(v+1+T(0.5f), Policy()); + term = exp(term); + } + else + term = pow(mult, T(v)) / boost::math::tgamma(v+1+T(0.5f), Policy()); mult *= -mult; } T operator()() @@ -142,6 +148,11 @@ inline T sph_bessel_j_imp(unsigned n, T x, const Policy& pol) if(n == 0) return boost::math::sinc_pi(x, pol); // + // Special case for x == 0: + // + if(x == 0) + return 0; + // // When x is small we may end up with 0/0, use series evaluation // instead, especially as it converges rapidly: // @@ -752,3 +763,4 @@ inline OutputIterator cyl_neumann_zero(T v, #endif // BOOST_MATH_BESSEL_HPP +