Re: [Boost-users] Boost_1_54_0: boost::math::sph_bessel(n, x) crash at x=0 n>0

Subject: Re: [Boost-users] Boost_1_54_0: boost::math::sph_bessel(n, x) crash at x=0 n>0
From: campagne (campagne_at_[hidden])
Date: 2013-10-17 14:52:48

Thanks
That's true that for x close to zero but positive there are also problems to the use of power series close to 0 is certainly the best thing to do. I am astonished that this not yet used.
Jean Eric

Envoyé depuis un mobile Samsung

-------- Message d'origine --------
De : John Maddock <john_at_[hidden]> 
Date :  
A : boost-users_at_[hidden] 
Objet : Re: [Boost-users] Boost_1_54_0: boost::math::sph_bessel(n,	x) crash at x=0 n>0 
 
> I'am running Boost 1_54_0 on Scientific Linux  slc-6.1,  gcc-4.4.5, 64 
> bits machine.
>
> I face a problem with computation of spherical bessel function, although 
> the problem seems well known here
> is a test program

Not well known to me - there are a number of fixes to the (cyclic) Bessel 
functions due out with 1.55, but I'm afraid this isn't one of them.  BTW the 
issue effects all x < 1, not just x == 0.

Here's the patch I'm about to commit to SVN:

Index: boost/math/special_functions/bessel.hpp
===================================================================
--- boost/math/special_functions/bessel.hpp (revision 86329)
+++ boost/math/special_functions/bessel.hpp (working copy)
@@ -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 neceesarily 
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:
    //
Index: libs/math/test/test_bessel_j.hpp
===================================================================
--- libs/math/test/test_bessel_j.hpp (revision 86329)
+++ libs/math/test/test_bessel_j.hpp (working copy)
@@ -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);

HTH, John. 

_______________________________________________
Boost-users mailing list
Boost-users_at_[hidden]
http://lists.boost.org/mailman/listinfo.cgi/boost-users