Bug reports for the GNU Scientific Library --from madongyu

[madongyu]

Dear GSL Developers：

   My name is Madongyu. I am a postgraduate student majoring in program 
analysis of computer science at Beijing University of Posts and 
Telecommunications. My research group is developing a symbolic execution tool 
called YUSE to detect floating-point exceptions in numerical code. We've found 
inputs that cause some of the GSL 2.7.1 functions to throw floating-point 
exceptions. I've written this mail to see if you consider any of them(8 
results) to be real bugs. If you could confirm whether any of these exceptions 
are useful warnings or actual bugs, we would very appreciate it. Because 
authoritative confirmation from GSL Developers to these exceptions is very 
important for us. We will release our tool as open source immediately if GSL 
Developers confirm for us that these are real bugs or useful warnings.

The purpose of our tool is to comprehensively detect floating-point exceptions 
in the program. For each category of floating-point exception we found, we 
listed 1 or 2. 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. Category 1：invalid operation exception 1.The function 
gsl_sf_conicalP_xlt1_large_neg_mu_e throws an invalid exception when its input 
tau=1.0, mu=2.0, x=3.0 at line 221 in file /gsl-2.7.1/specfunc/legendre_con.c: 
double p = x/sqrt(beta2*(1.0-x*x) + 1.0); Because the parameter of sqrt cannot 
be negative numbers, invalid exception occurs in this occasion.


2.The function gsl_cdf_laplace_Qinv throws an invalid exception when its input 
Q=-0.5 at line 67 in file /gsl-2.7.1/cdf/laplaceinv.c：x = -a * log(2*Q);Because 
the parameter of log cannot be negative numbers, invalid exception occurs in 
this occasion. Category 2：overflow 3.The function 
gsl_sf_bessel_Knu_scaled_asymp_unif_e throws a overflow exception when its 
input nu=1.45e-192, x=1.0 at line 398 in file /gsl-2.7.1/specfunc/bessel.c:  
double root_term = hypot(1.0,z);


4.The function gsl_sf_conicalP_xlt1_large_neg_mu_e throws a overflow exception 
when its input tau =1.23e189,mu=1e-2 at line 219 in file 
/gsl-2.7.1/specfunc/legendre_con.c: double beta2 = beta*beta; Category 
3：underflow 5.The function gsl_sf_bessel_Jnu_asympx_e throws an underflow 
exception when its input nu=1.23e-189, x=1.0 at line 217 in file 
/gsl-2.7.1/specfunc/bessel.c:  double mu = 4.0*nu*nu; Category 
4：division-by-zero 6.The function gsl_sf_bessel_Knu_scaled_asympx_e throws an 
division-by-zero floating-point exception when its input x=0.0 at line 318 in 
file /gsl-2.7.1/specfunc/bessel.c: double r = nu/x; Category 5：type conversion 
7.The function gsl_ran_gamma_knuth throws an overflow caused by type conversion 
exception when its input a=1.23e189 at line 44 in file 
/gsl-2.7.1/randist/gamma.c: unsigned int na = floor (a); Category 
6：inexact(rounding error) 8.The function hyperg_1F1_1 throws a rounding error 
exception which cause the while loop to not terminate when its input 
b=x=1.23e189 at line 378 in file /gsl-2.7.1/specfunc/hyperg_1F1.c: while(bp 
> b+0.1) {.The loop at 378 will not terminate when x reaches 1.23e189 (a 
very large value),because bp -= 1.0 at line 380 makes bp – 1.0 is rounded to bp 
itself, thus the while condition at bp > b+0.1 at 378 is satisfied all the 
time and can not terminate. I am very sorry about that I have disturbed you in 
your busy schedule. Authoritative confirmation from GSL Developers to these 
exceptions is very important for us. I am looking forward to your reply. Thank 
you very much.




Best wishes













From Madongyu

Beijing University of Posts and Telecommunications