diff --git a/include/boost/math/special_functions/detail/bessel_ik.hpp b/include/boost/math/special_functions/detail/bessel_ik.hpp index 0c653b4753..2ffb0f060c 100644 --- a/include/boost/math/special_functions/detail/bessel_ik.hpp +++ b/include/boost/math/special_functions/detail/bessel_ik.hpp @@ -326,6 +326,11 @@ int bessel_ik(T v, T x, T* result_I, T* result_K, int kind, const Policy& pol) T scale = 1; T scale_sign = 1; + n = iround(v, pol); + u = v - n; // -1/2 <= u < 1/2 + BOOST_MATH_INSTRUMENT_VARIABLE(n); + BOOST_MATH_INSTRUMENT_VARIABLE(u); + if (((kind & need_i) == 0) && (fabs(4 * v * v - 25) / (8 * x) < tools::forth_root_epsilon())) { // A&S 9.7.2 @@ -337,11 +342,6 @@ int bessel_ik(T v, T x, T* result_I, T* result_K, int kind, const Policy& pol) } else { - n = iround(v, pol); - u = v - n; // -1/2 <= u < 1/2 - BOOST_MATH_INSTRUMENT_VARIABLE(n); - BOOST_MATH_INSTRUMENT_VARIABLE(u); - BOOST_MATH_ASSERT(x > 0); // Error handling for x <= 0 handled in cyl_bessel_i and cyl_bessel_k // x is positive until reflection @@ -414,6 +414,7 @@ int bessel_ik(T v, T x, T* result_I, T* result_K, int kind, const Policy& pol) } if (reflect) { + BOOST_MATH_ASSERT(fabs(v - n - u) < tools::forth_root_epsilon()); T z = (u + n % 2); T fact = (2 / pi()) * (boost::math::sin_pi(z, pol) * Kv); if(fact == 0) diff --git a/test/test_bessel_j.hpp b/test/test_bessel_j.hpp index 82106213ea..c54ad4fa8b 100644 --- a/test/test_bessel_j.hpp +++ b/test/test_bessel_j.hpp @@ -288,6 +288,7 @@ void test_bessel(T, const char* name) BOOST_CHECK_EQUAL(boost::math::cyl_bessel_j(T(0), std::numeric_limits::infinity()), T(0)); BOOST_CHECK_EQUAL(boost::math::cyl_bessel_j(T(2), std::numeric_limits::infinity()), T(0)); BOOST_CHECK_EQUAL(boost::math::cyl_bessel_j(T(1.25), std::numeric_limits::infinity()), T(0)); + BOOST_CHECK_EQUAL(boost::math::cyl_bessel_j(T(-1.25), std::numeric_limits::infinity()), T(0)); BOOST_CHECK_EQUAL(boost::math::sph_bessel(0, std::numeric_limits::infinity()), T(0)); BOOST_CHECK_EQUAL(boost::math::sph_bessel(1, std::numeric_limits::infinity()), T(0)); BOOST_CHECK_EQUAL(boost::math::sph_bessel(2, std::numeric_limits::infinity()), T(0));