[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Bug-gsl] Potential Bugs?
From: |
Brian Gough |
Subject: |
Re: [Bug-gsl] Potential Bugs? |
Date: |
Tue, 17 Aug 2010 21:38:19 +0100 |
User-agent: |
Wanderlust/2.15.6 (Almost Unreal) Emacs/23.2 Mule/6.0 (HANACHIRUSATO) |
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