> 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" :-(

>>Ok, I could swear it worked for me but like you said, if the input is nudged a bit the issue goes away (or returns).

> 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.

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

_____________________________________________
Boost-users mailing list
Boost-users@lists.boost.org
http://lists.boost.org/mailman/listinfo.cgi/boost-users