Re: Besseli Function

[David Bateman]

Since I put my foot in my mouth saying this was trivial, I tried it
out last night. I have two versions of a patch, but a question on
which to choose. The problem is the AMOS package doesn't signal
underflows correctly. For example consider

octave:1> x = [-90,90];
octave:2> y0 = besselk(0,x);
octave:3> y1 = exp(-x).*besselk(0,x,true);

Without my original patch this gives

y0 =

  Inf + Infi    0 +   0i

y1 =

  0.0000e+00 - 1.6145e+38i  1.0810e-40 + 0.0000e+00i

With the patch it gives 

y0 =

  0.0000e+00 - 1.6145e+38i  0.0000e+00 + 0.0000e+00i

y1 =

  0.0000e+00 - 1.6145e+38i  1.0810e-40 + 0.0000e+00i


The problem is that the AMOS function zbesk returns ierr=2 for an overflow,
but does not signal an underflow condition. So the first patch is half
a solution.

My question is then, why not just used the scaled calculation always? A quick
test with the unpatched code shows that there seems to be little speed 
difference. 

octave:1> x = 20*rand(100,100)-10; tic; y0 = airy(0,x); toc
ans = 0.14151
octave:2> x = 20*rand(100,100)-10; tic; y1 = exp( - 2. / 3. .* x .* 
sqrt(x)).*airy(0,x,true); toc
ans = 0.13462

Is there possibly a problem with precision? If not I propose to just used
the scaled calculations always internally, and return the scaled or unscaled 
values depending on the the optional argument..

Frankly, without changing the underlying Fortran code, I don't see any
other way to treat this problem than the two solutions proposed. Barring 
any precision issues I propose that the second solution is better (faster,
treats underflows).

Regards
David

According to John W. Eaton <address@hidden> (on 09/08/03):
> On  8-Sep-2003, David Bateman <address@hidden> wrote:
> 
> | 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.
> 
> Why not leave the option for scaling, but try scaling if it was not
> requested and overflow has occurred?
> 
> But in any case, I'm asking that someone who wants this feature
> provide a complete patch rather than waiting for me to get around to
> doing it.
> 
> Thanks,
> 
> jwe

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

patch1
Description: Text document

patch2
Description: Text document