From 964435a350e5ab1aec09ae6bc025be4efc04bfa0 Mon Sep 17 00:00:00 2001 From: John Maddock Date: Fri, 18 Oct 2013 08:56:42 +0000 Subject: [PATCH] Previous commit failed for types with an extended exponent range - use logarithms rather than assuming the result is zero. [SVN r86346] --- include/boost/math/special_functions/bessel.hpp | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/include/boost/math/special_functions/bessel.hpp b/include/boost/math/special_functions/bessel.hpp index e96d4e79d..ceca30916 100644 --- a/include/boost/math/special_functions/bessel.hpp +++ b/include/boost/math/special_functions/bessel.hpp @@ -50,11 +50,8 @@ struct sph_bessel_j_small_z_series_term mult = x / 2; if(v + 3 > max_factorial::value) { - // term = v * log(mult) - boost::math::lgamma(v+1+T(0.5f), Policy()); - // term = exp(term); - // Denominator increases faster than numerator each time v increases by one, - // so if tgamma would overflow then the result is necessarily zero. - term = 0; + 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());