[Boost-commit] svn:boost r86343 - in trunk: boost/math/special_functions libs/math/test

Subject: [Boost-commit] svn:boost r86343 - in trunk: boost/math/special_functions libs/math/test
From: john_at_[hidden]
Date: 2013-10-17 14:21:56

Author: johnmaddock
Date: 2013-10-17 14:21:56 EDT (Thu, 17 Oct 2013)
New Revision: 86343
URL: http://svn.boost.org/trac/boost/changeset/86343

Log:
Fix for sph_bessel when v is large and x is small.

Text files modified: 
   trunk/boost/math/special_functions/bessel.hpp |    17 ++++++++++++++++-                       
   trunk/libs/math/test/test_bessel_j.hpp        |     7 +++++++                                 
   2 files changed, 23 insertions(+), 1 deletions(-)

Modified: trunk/boost/math/special_functions/bessel.hpp
==============================================================================
--- trunk/boost/math/special_functions/bessel.hpp	Thu Oct 17 13:34:31 2013	(r86342)
+++ trunk/boost/math/special_functions/bessel.hpp	2013-10-17 14:21:56 EDT (Thu, 17 Oct 2013)	(r86343)
@@ -48,7 +48,16 @@
    {
       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<T>::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;
+      }
+      else
+         term = pow(mult, T(v)) / boost::math::tgamma(v+1+T(0.5f), Policy());
       mult *= -mult;
    }
    T operator()()
@@ -142,6 +151,11 @@
    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 +766,4 @@
 
 #endif // BOOST_MATH_BESSEL_HPP
 
+

Modified: trunk/libs/math/test/test_bessel_j.hpp
==============================================================================
--- trunk/libs/math/test/test_bessel_j.hpp	Thu Oct 17 13:34:31 2013	(r86342)
+++ trunk/libs/math/test/test_bessel_j.hpp	2013-10-17 14:21:56 EDT (Thu, 17 Oct 2013)	(r86343)
@@ -262,6 +262,13 @@
     do_test_sph_bessel_j<T>(sph_bessel_data, name, "Bessel j: Random Data");
 
     //
+    // Some special cases:
+    //
+    BOOST_CHECK_EQUAL(boost::math::sph_bessel(0, T(0)), T(1));
+    BOOST_CHECK_EQUAL(boost::math::sph_bessel(1, T(0)), T(0));
+    BOOST_CHECK_EQUAL(boost::math::sph_bessel(100000, T(0)), T(0));
+
+    //
     // Special cases that are errors:
     //
     BOOST_CHECK_THROW(boost::math::cyl_bessel_j(T(-2.5), T(0)), std::domain_error);