diff options
Diffstat (limited to 'boost/math/special_functions/detail/bessel_ik.hpp')
| -rw-r--r-- | boost/math/special_functions/detail/bessel_ik.hpp | 23 |
1 files changed, 22 insertions, 1 deletions
diff --git a/boost/math/special_functions/detail/bessel_ik.hpp b/boost/math/special_functions/detail/bessel_ik.hpp index a589673ffb..10118d9715 100644 --- a/boost/math/special_functions/detail/bessel_ik.hpp +++ b/boost/math/special_functions/detail/bessel_ik.hpp @@ -234,6 +234,7 @@ int CF2_ik(T v, T x, T* Kv, T* Kv1, const Policy& pol) BOOST_MATH_INSTRUMENT_VARIABLE(b); BOOST_MATH_INSTRUMENT_VARIABLE(D); BOOST_MATH_INSTRUMENT_VARIABLE(f); + for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++) // starting from 2 { // continued fraction f = z1 / z0 @@ -250,10 +251,27 @@ int CF2_ik(T v, T x, T* Kv, T* Kv1, const Policy& pol) C *= -a / k; Q += C * q; S += Q * delta; + // + // Under some circumstances q can grow very small and C very + // large, leading to under/overflow. This is particularly an + // issue for types which have many digits precision but a narrow + // exponent range. A typical example being a "double double" type. + // To avoid this situation we can normalise q (and related prev/current) + // and C. All other variables remain unchanged in value. A typical + // test case occurs when x is close to 2, for example cyl_bessel_k(9.125, 2.125). + // + if(q < tools::epsilon<T>()) + { + C *= q; + prev /= q; + current /= q; + q = 1; + } // S converges slower than f BOOST_MATH_INSTRUMENT_VARIABLE(Q * delta); BOOST_MATH_INSTRUMENT_VARIABLE(abs(S) * tolerance); + BOOST_MATH_INSTRUMENT_VARIABLE(S); if (abs(Q * delta) < abs(S) * tolerance) { break; @@ -261,7 +279,10 @@ int CF2_ik(T v, T x, T* Kv, T* Kv1, const Policy& pol) } policies::check_series_iterations<T>("boost::math::bessel_ik<%1%>(%1%,%1%) in CF2_ik", k, pol); - *Kv = sqrt(pi<T>() / (2 * x)) * exp(-x) / S; + if(x >= tools::log_max_value<T>()) + *Kv = exp(0.5f * log(pi<T>() / (2 * x)) - x - log(S)); + else + *Kv = sqrt(pi<T>() / (2 * x)) * exp(-x) / S; *Kv1 = *Kv * (0.5f + v + x + (v * v - 0.25f) * f) / x; BOOST_MATH_INSTRUMENT_VARIABLE(*Kv); BOOST_MATH_INSTRUMENT_VARIABLE(*Kv1); |