Boost_1_54_0: boost::math::sph_bessel(n, x) crash at x=0 n>0
Dear developpers,
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
{
int L;
try {
for(L=100;L<2010;L++) {
boost::math::sph_bessel(L, 0.);
}
} catch (std::exception& sex) {
cout << "Error j_"<
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.
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@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
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.
There are no changes to operator()(), the n variable appears in the next part of the diff in inline T sph_bessel_j_imp(unsigned n, T x, const Policy& pol)
2) is there a way to get your current version of the bessel.hpp file ?
Attached. Note I had to change the previous fix somewhat as it failed for types with extended exponent ranges (I realised that just as I was dozing off to sleep last night!)
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.
Sorry I don't understand: the problem was occurring inside the power series for small x. John.
Your fix pass all my personnal tests. Thanks. Sorry for the trouble, according to the limit close to zero. Bye the way do you know when 1_55 version will be publicly available? Jean-Eric Le 18/10/2013 11:01, John Maddock a écrit :
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.
There are no changes to operator()(), the n variable appears in the next part of the diff in
inline T sph_bessel_j_imp(unsigned n, T x, const Policy& pol)
2) is there a way to get your current version of the bessel.hpp file ?
Attached.
Note I had to change the previous fix somewhat as it failed for types with extended exponent ranges (I realised that just as I was dozing off to sleep last night!)
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.
Sorry I don't understand: the problem was occurring inside the power series for small x.
John.
_______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
From: Boost-users [mailto:boost-users-bounces@lists.boost.org] On Behalf Of JE Campagne Sent: Friday, October 18, 2013 11:22 AM To: boost-users@lists.boost.org Subject: Re: [Boost-users] Boost_1_54_0: boost::math::sph_bessel(n, x) crash at x=0 n>0 Your fix pass all my personnal tests. Thanks. J Sorry for the trouble, according to the limit close to zero. Ensuring a 'good failure' is nearly as troublesome as calculating the value! Bye the way do you know when 1_55 version will be publicly available? There are delays trying to ensure that Boost 1.55 compiles with the latest MS compiler. But probably only a week or two? Paul
-----Original Message----- From: Boost-users [mailto:boost-users-bounces@lists.boost.org] On Behalf Of JE Campagne Sent: Friday, October 18, 2013 9:09 AM To: boost-users@lists.boost.org Subject: Re: [Boost-users] Boost_1_54_0: boost::math::sph_bessel(n, x) crash at x=0 n>0
Hi
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.
Ask the organ grinder ;-)
2) is there a way to get your current version of the bessel.hpp file ?
http://svn.boost.org/svn/boost/trunk/boost/math/special_functions/bessel.hpp in case you want to try this now. (There might be other dependencies but I guess not). Revision: 86346 Author: johnmaddock Date: 18 October 2013 09:56:42 Message: Previous commit failed for types with an extended exponent range - use logarithms rather than assuming the result is zero. ---- Modified : /trunk/boost/math/special_functions/bessel.hpp Revision: 86343 Author: johnmaddock Date: 17 October 2013 19:21:56 Message: Fix for sph_bessel when v is large and x is small. ---- Modified : /trunk/boost/math/special_functions/bessel.hpp Modified : /trunk/libs/math/test/test_bessel_j.hpp Paul --- Paul A. Bristow, Prizet Farmhouse, Kendal LA8 8AB UK +44 1539 561830 07714330204 pbristow@hetp.u-net.com
participants (3)
-
JE Campagne -
John Maddock -
Paul A. Bristow