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@johnmaddock.co.uk> Date : A : boost-users@lists.boost.org 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@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users