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
