From bug-octave-request at bevo dot che dot wisc dot edu  Thu Jun  7 15:09:16 2001
Subject: besseli output range apparently limited to about 10^33
From: "John W. Eaton" <jwe at bevo dot che dot wisc dot edu>
To: Daniel Alt <daniel dot alt at ttpcom dot com>
Cc: bug-octave at bevo dot che dot wisc dot edu
Date: Thu, 7 Jun 2001 15:08:29 -0500

On  4-Jun-2001, Daniel Alt <daniel dot alt at ttpcom dot com> wrote:

| "Bug" report for Octave 2.1.31 configured for MS Windows
| (octave-windows-2000oct25d.exe downloaded 1 June 2001 from
| sourceforge.net)
| 
| (I attempted, unsuccessfully, to use Octave's bug_report command to
| generate this email.  My PC doesn't have Unix shell command emulators
| installed, so bug_report fails with "not found" errors for grep, cat,
| chmod, etc.)
| 
| It might be stretching the definition to call this a bug, but besseli's
| output goes to infinity at around 10^33 (instead of at the maximum
| floating point value of about 10^308).  Presumably this is because some
| intermediate result exceeds the floating point range, but both Matlab
| and Microsoft Excel can calculate besseli values of at least 10^300.
| 
| Here is a bit of an Octave transcript showing the transition:
| 
| octave:87> besseli(19,81.3565)
| ans =  1.0327e+33
| octave:88> besseli(19,81.357)
| ans = Inf + Infi
| 
| As one might expect, besselj has similar behaviour for equivalent
| inputs:
| 
| octave:92> besselj(19,81.357i)
| ans = Inf + Infi
| octave:93> besselj(19,81.3565i)
| ans = -6.3238e+16 - 1.0327e+33i
| 
| (One question on this besselj result:  Shouldn't ans be purely imaginary
| in the second case?)

Octave uses the subroutines provided by the "Amos" library from Netlib
to compute Bessel functions.

You can ask for the results to be scaled to avoid overflow.  For
example, I noticed that

  besseli (19, 700)

fails, but

  exp (700) * besseli (19, 700, 1)

produces the same result as Matlab does for besseli (19, 700).

I'm not sure why the Fortran code doesn't get this right.

I'm not sure that it would be correct for Octave to simply call the
Fortran code and request the scaled result, then multiply that by the
inverse scale factor before returning.  It seems that it might not
always be as accurate.

I'd appreciate it if anyone could shed some light on this and maybe
provide a fix.

Thanks,

jwe



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------