Re: [Bug-gsl] Potential Bugs?

[Brian Gough <bjg@gnu.org>]

At Wed, 11 Aug 2010 17:14:34 -0700,
Thanh Vo wrote:
> My name's Thanh Vo. My research group is developing a tool to detect
> floating point exceptions in numerical code. We've found inputs that
> cause some of
> the GSL 1.14 functions to throw floating-point exceptions.  I've written
> this post to see if you consider any of them to be real bugs.

Hi,

Thanks for your email.  I have made some comments below.  This sounds
like a useful tool, I hope you will be able to release it as free
software.

> 
> 1.  The function airy_aie throws a divide by zero exception when its
> input x = 0 at line 618 in file gsl/specfunc/airy.c : double z   =
> 2.0/(x*sqx) - 1.0;.  It throws an invalid exception when x = -1 at
> line 617:  double sqx = sqrt(x);.
> 
> 2.  The function airy_bie throws a divide by zero exception when its
> input x = 0 at line 635 in file gsl/specfunc/airy.c: double z   =
> ATR/(x*sqx) + BTR;. It throws an invalid exception when x=-1 at line
> 634: double sqx = sqrt(x);.

These would be bugs but the functions are static and only called with
x > 1 or x > 2 so the exceptions cannot occur in practice.

> 3.  The function gsl_sf_asymp_thetanu_corr_e throws a divide by zero
> exception when its input x = 0 and nu != 0 at the line 181 in file
> gsl/specfunc/bessel_amp_phase.c: const double r  = 2.0*nu/x;. It
> throws an invalid exception when x=0 and nu=0 at line 181: const
> double r  = 2.0*nu/x;.
>
> 4.  The function gsl_sf_bessel_asymp_Mnu_e throws a divide by zero
> exception when its input x=0 and nu !=0 at the line 167 in file
> gsl/specfunc/bessel_amp_phase.c: const double r  = 2.0*nu/x;. It
> throws an invalid exception when x=0 and nu=0 at line 167: const
> double r  = 2.0*nu/x;.

It's a valid warning. The functions are publicly accessible, but not
documented.  Within the library they are only called with large values
of x >= 1000 so the exception does not occur in practice.  Ideally
they should be static functions.
 
> 5.  The function gsl_sf_bessel_J0 throws a overflow exception when its
> input x = 5.85177e+145 at the line 88 in the file
> gsl/specfunc/bessel_J0.c: const double z = 32.0/(y*y) - 1.0;. It
> throws an underflow exception when x = 6.8481e-155 at the line 81:
> result->err = y*y;.

The overflow exception is a valid warning.  It's also a bug in the
sense that it is avoidable (although the end result is correct, z=-1
since 1/inf=0).  It would be better written as (32/y)/y.

For the underflow it's a valid warning.  There isn't really any
alternative to returning zero in that case though.

-- 
Brian Gough