Re: besseli: Am I missing something?

[Jordi Gutiérrez Hermoso]

On 11 January 2012 23:52, Robert T. Short <address@hidden> wrote:
> Before I file a bug report, maybe I should check that I am not all wet.
>
> If I recall my Bessel function theory,
>
> besseli(n,x) = besseli(-n,x)
>
> if n is an integer.

Yes, 9.6.27 on page 376 of Abramowitz & Stegun. Do you really have
this memorised? I can barely keep straight which letter goes with
which Bessel function.

> In octave 3.4.2 I get
>
> octave:3> besseli(1,1),besseli(-1,1)
> ans =  0.565159103992485
> ans =  0.565159103992485
>
> Looks good for n=1, but for n=10
>
> octave:4> besseli(10,1),besseli(-10,1)
> ans =  2.75294803983687e-10
> ans = -1.40610343427994e-07
>
> Not so good!
>
> And for n=100
>
> octave:5> besseli(100,1),besseli(-100,1)
> ans =  8.47367400813812e-189
> ans =  7.38028373423800e+170

This does seem like a considerable loss of precision... I compared
with the implementation in Scipy, which seems to get it right, and
with pari, which seems to error out due to a negative value of Gamma
(I wonder if it's trying to evaluate with the infinite series in terms
of Gamma, which is not what A&S recommends).

The code that actually does this is in libcruft, the AMOS library...
and there, well, FORTRANum est; non legitur... but before it gets to
AMOS, the Octave code attempts to use some identities to evaluate
besseli for negative alpha, right here:

    
http://hg.savannah.gnu.org/hgweb/octave/file/1eeb3a8d5708/liboctave/lo-specfun.cc#l831

Can you see a problem there? I can't follow the identities between I
and K that seem to be at play. We could simply specifically check for
integer alpha and just use besseli(-alpha, z) if alpha is negative,
but perhaps you can glance at A&S or Watson and see if you can find in
there an identity that will work better in general. ;-)

HTH,
- Jordi G. H.