1. I tried your pseudo code below Y[8.5](4*Pi)
using my value of pi (I believe that it was in an earlier post) and did not get a bad answer. so if this is really a cancellation error could it be related to the precision of pi (how many decimals did you use)?
I used: boost::math::cyl_neumann(8.5, 4 * pi); using your value of pi to reproduce the issue - I did say "exquisite" cancellation though - the values have to be just right - down to the last bit - to cause the issue. In fact the cancellation problem occurs after quite a bit of calculation so you have to be very unlucky indeed for things to work out "just so" :-(
2. Since we are calling bessel of a half integer order I would suspect that the Gamma( ) function is called. Could a bad value from that function caused the error? If so I might expect the debug statement to trace back that far but they didn't.
No it's an error in the logic inside bessel_jy in bessel_jy.hpp, I'm testing a fix now.
I have GSL on a linux machine at work and plan to try that too.
GSL doesn't support bessel J or Y with non-integer orders (unless they're changed since I last looked), they do support the spherical bessel functions but use a different calculation method so it probably works there. Still testing yours, John.