Re: [Boost-users] spherical bessel functions

Subject: Re: [Boost-users] spherical bessel functions
From: John Maddock (boost.regex_at_[hidden])
Date: 2012-03-23 14:55:24

>It sounds like you've tracked it down. Thanks. David.

Here's the patch to the headers:

Index: bessel_jy.hpp
===================================================================
--- bessel_jy.hpp (revision 77307)
+++ bessel_jy.hpp (working copy)
@@ -491,6 +485,16 @@
            CF2_jy(u, x, &p, &q, pol);                  // continued 
fraction CF2_jy
            T t = u / x - fu;                   // t = J'/J
            gamma = (p - t) / q;
+           //
+           // We can't allow gamma to cancel out to zero competely as it 
messes up
+           // the subsequent logic.  So pretend that one bit didn't cancel 
out
+           // and set to a suitably small value.  The only test case we've 
been able to
+           // find for this, is when v = 8.5 and x = 4*PI.
+           //
+           if(gamma == 0)
+           {
+              gamma = u * tools::epsilon<T>() / x;
+           }
            Ju = sign(current) * sqrt(W / (q + gamma * (p - t)));

            Jv = Ju * ratio;                    // normalization

I'm re-running the full tests now before committing (there are changes to 
the tests not shown above), but the Bessel function tests all look good so 
far,

John.