Re: [Boost-users] spherical bessel functions

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

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