[Boost-commit] svn:boost r84951 - trunk/boost/math/special_functions/detail

Subject: [Boost-commit] svn:boost r84951 - trunk/boost/math/special_functions/detail
From: john_at_[hidden]
Date: 2013-07-04 11:31:32

Author: johnmaddock
Date: 2013-07-04 11:31:32 EDT (Thu, 04 Jul 2013)
New Revision: 84951
URL: http://svn.boost.org/trac/boost/changeset/84951

Log:
Fix overflow/underflow errors when x is very close to 2.

Text files modified: 
   trunk/boost/math/special_functions/detail/bessel_ik.hpp |    17 +++++++++++++++++                       
   1 files changed, 17 insertions(+), 0 deletions(-)

Modified: trunk/boost/math/special_functions/detail/bessel_ik.hpp
==============================================================================
--- trunk/boost/math/special_functions/detail/bessel_ik.hpp	Thu Jul  4 05:13:23 2013	(r84950)
+++ trunk/boost/math/special_functions/detail/bessel_ik.hpp	2013-07-04 11:31:32 EDT (Thu, 04 Jul 2013)	(r84951)
@@ -234,6 +234,7 @@
     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,6 +251,22 @@
         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);