[Bug-gsl] Bug in sys/ldfrexp.c

[Chris Mihelich]

Good evening,

I wish to report and correct the following bugs in sys/ldfrexp.c:

(1) gsl_ldexp computes a power of 2 (p2) that may overflow or
underflow even though the correct answer is representable;
consequently the return value x * p2 may be spuriously zero or
infinity.

(2) gsl_frexp can hang (this is what bit me in the rear at work
yesterday) if its call to gsl_ldexp spuriously returns zero or
infinity as in (1).  The first while loop divides infinity by 2.0 over
and over and wonders why it doesn't get any smaller; the second while
loop has an analogous problem with zero.

I have attached a tgz ball containing revised copies of the following
files hoping that you will correct these faults in the next release:

(a) sys/ldfrexp.c, rewritten from scratch.  You will want to put
suitable identifying comments at the top.
(b) sys/test.c, to which I have added very thorough tests for
gsl_ldexp and gsl_frexp.
(c) test/results.c, to which I have added a function
gsl_test_rel_subnormals_ok that is just like gsl_test_rel except for
not automatically failing if an argument is subnormal.  The additions
in (b) need this.
(d) err/gsl_test.h, to which I have added the declaration of
gsl_test_rel_subnormals_ok.

Thank you for continuing to maintain GSL.  Aside from this cute little
glitch, it has been a very helpful tool for me in my day job.

Chris Mihelich (address@hidden, 609-658-5827)
"Never made it as a wise man"

gsl_ldfrexp_bug.tar.gz
Description: GNU Zip compressed data