boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Fix for sph_bessel when v is large and x is small.

Browse Source 
[SVN r86343]
This commit is contained in:
John Maddock
2013-10-17 18:21:56 +00:00
parent 7660f2dd3c
commit 57b40ac443
2 changed files with 23 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

17
include/boost/math/special_functions/bessel.hpp
View File

@@ -48,7 +48,16 @@ struct sph_bessel_j_small_z_series_term
{
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 necessarily 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 @@ inline T sph_bessel_j_imp(unsigned n, T x, const Policy& pol)
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:
//
@@ -752,3 +766,4 @@ inline OutputIterator cyl_neumann_zero(T v,
#endif // BOOST_MATH_BESSEL_HPP

7
test/test_bessel_j.hpp
View File

@@ -261,6 +261,13 @@ void test_bessel(T, const char* name)
#include "sph_bessel_data.ipp"
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:
//