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.