diff --git a/include/boost/math/special_functions/detail/bessel_jn.hpp b/include/boost/math/special_functions/detail/bessel_jn.hpp index 4e683a7c3..c3a5def87 100644 --- a/include/boost/math/special_functions/detail/bessel_jn.hpp +++ b/include/boost/math/special_functions/detail/bessel_jn.hpp @@ -84,7 +84,11 @@ BOOST_MATH_GPU_ENABLED T bessel_jn(int n, T x, const Policy& pol) current = value; } } - else if((x < 1) || ((n > x * x / 4) && (x < 5))) + else if(x < 1) + { + return factor * bessel_j_small_z_series(T(n), x, pol); + } + else if((x < 5) && (n > x * x / 4)) { return factor * bessel_j_small_z_series(T(n), x, pol); } diff --git a/include/boost/math/special_functions/detail/bessel_jy.hpp b/include/boost/math/special_functions/detail/bessel_jy.hpp index 264b93e05..3b3e20805 100644 --- a/include/boost/math/special_functions/detail/bessel_jy.hpp +++ b/include/boost/math/special_functions/detail/bessel_jy.hpp @@ -327,7 +327,12 @@ namespace boost { namespace math { // x is positive until reflection W = T(2) / (x * pi()); // Wronskian T Yv_scale = 1; - if(((kind & need_y) == 0) && ((x < 1) || ((v > x * x / 4) && (x < 5)))) + + const bool kind_does_not_need_y { ((kind & need_y) == 0) }; + + const bool x_is_lt_one { (x < 1) }; + + if(kind_does_not_need_y && x_is_lt_one) { // // This series will actually converge rapidly for all small @@ -337,7 +342,17 @@ namespace boost { namespace math { Jv = bessel_j_small_z_series(v, x, pol); Yv = boost::math::numeric_limits::quiet_NaN(); } - else if((x < 1) && (u != 0) && (log(policies::get_epsilon() / 2) > v * log((x/2) * (x/2) / v))) + else if(kind_does_not_need_y && ((x < 5) && (v > x * x / 4))) + { + // + // This series will actually converge rapidly for all small + // x - say up to x < 20 - but the first few terms are large + // and divergent which leads to large errors :-( + // + Jv = bessel_j_small_z_series(v, x, pol); + Yv = boost::math::numeric_limits::quiet_NaN(); + } + else if(x_is_lt_one && (u != 0) && (log(policies::get_epsilon() / 2) > v * log((x/2) * (x/2) / v))) { // Evaluate using series representations. // This is particularly important for x << v as in this