help-octave
[Top][All Lists]

## Re: Besseli Function

 From: David Bateman Subject: Re: Besseli Function Date: Tue, 9 Sep 2003 11:14:44 +0200 User-agent: Mutt/1.3.28i

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

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

[ ] Motorola Internal Use Only
[ ] Motorola Confidential Proprietary
```

patch1
Description: Text document

patch2
Description: Text document