help-octave
[Top][All Lists]

## Re: Besseli Function

 From: David Bateman Subject: Re: Besseli Function Date: Mon, 8 Sep 2003 19:11:35 +0200 User-agent: Mutt/1.3.28i

```Why can't we just replace the code in lo-specfun.cc so that an overflow
without scaling, automatically forces a scaling. Something like

F77_FUNC (zbesi, ZBESI) (zr, zi, alpha, 0, 1, &yr, &yi, nz, ierr);

if (ierr == 2)
{
// Overflow condition. Try scaling the result
F77_FUNC (zbesi, ZBESI) (zr, zi, alpha, 1, 1, &yr, &yi, nz, ierr);

yr *= exp(z);
yi *= exp(z);

}

if (zi == 0.0 && zr > 0.0)
yi = 0.0;

retval = bessel_return_value (Complex (yr, yi), ierr);

Then the whole issue goes away. The patch the lo-specfun.cc is pretty
trivial. However, a lot of dead code will need to be removed from
a dozen different places, starting with besselj.cc and working downwards.

D.

According to John W. Eaton <address@hidden> (on 09/08/03):
> On  8-Sep-2003, Jose Otavio Bueno <address@hidden> wrote:
>
> | I'm having problems with the Besseli Functions in Octave. While with
> | Mathlab it is possible to solve the besseli functions with values up to
> | 700, with Octave the function returns a valid answer just for input
> | values below 81:
>
> Octave allows you to ask to have the result from besseli scaled by
> exp(-abs(x)).  So you can write
>
>   exp (abs (x)) * besseli (0, x, true);
>
> (the last argument to besseli can be anything, its presence just means
> "scale the result").  By doing this, Octave returns
>
>   octave:1> x = [1,10,81,82,700];
>   octave:2> exp(x).*besseli (0, x, true)
>   ans =
>
>       1.2661e+00    2.8157e+03    6.6864e+33    1.8064e+34   1.5296e+302
>
> It might be useful if Octave detected the overflow (the internal
> function zbesi provides that information) and then attempted to
> perform the scaled calculation instead.  Can someone please look at
> providing a patch for that?
>
> Also, the documentation for the scaling option could be better, so a
> patch for that would also be appreciated.
>
> 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
> -------------------------------------------------------------

--
Motorola CRM                                 +33 1 69 35 48 04 (Ph)
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax)
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as:

[x] General Business Information
[ ] Motorola Internal Use Only
[ ] Motorola Confidential Proprietary

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

```

reply via email to