[Bug-gsl] Potential Bugs?

[Thanh Vo]

Hi,

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.

For each exception, we describe the name of function, the type of
exception thrown, the inputs that trigger the exception, and specific
location of the code that causes the exception.

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);.

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

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

If anyone could confirm whether any of these exceptions is an actual
bug, I would appreciate it.

Looking forward to hearing from you,

Best regards,

Thanh Vo.