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: JE Campagne (campagne_at_[hidden])
Date: 2013-10-18 04:08:39

Hi John

1) in your patch for T operator()() part, what's the "n" variable, 
looking at 1_54_0 version there is no mention to this variable.
2) is there a way to get your current version of the bessel.hpp file ?

Thanks
Jean-Eric

Le 17/10/2013 20:20, John Maddock a écrit :
>> 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