Re: Besseli Function

[David Bateman]

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

-- 
David Bateman                                address@hidden
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
-------------------------------------------------------------