octave-maintainers
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: besseli: Am I missing something?


From: Jussi Lehtola
Subject: Re: besseli: Am I missing something?
Date: Thu, 12 Jan 2012 11:00:03 +0200

On Thu, 12 Jan 2012 00:39:14 -0500
Jordi Gutiérrez Hermoso <address@hidden> wrote:
> > 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:

I'm not sure why this is done in Octave at all, as there is another GNU
project which specializes in the evaluation of scientific functions.
Furthermore, an Octave package already exists with interfaces these
two: octave-gsl.

octave:6> bessel_In(10,-1)
ans =  2.75294803983687e-10
octave:7> bessel_In(10,1)
ans =  2.75294803983687e-10
octave:8> bessel_In(100,1)
ans =  8.47367400813807e-189
octave:9> bessel_In(100,-1)
ans =  8.47367400813807e-189

I really think that you should switch to using GSL. IMHO there's no
sense in maintaining two versions of the wheel. The code in GSL is
probably much better tested for weird use cases.
-- 
Jussi Lehtola
Fedora Project Contributor
address@hidden


reply via email to

[Prev in Thread] Current Thread [Next in Thread]