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