Hi everyone,

I just tried to examine bug #21838 in gsl_cdf_fdist_P () (source file 
cdf/fdist.c), reported in

http://lists.gnu.org/archive/html/bug-gsl/2007-09/msg00023.html

with the attached testcase program (testcases, expected values and 
tolerance taken from cdf/test.c).

I examined the first testcase: It seems that the origin of the NaN is in 
beta_cont_frac () in the source file cdf/beta_inc.c. The loop aborts 
after the maximum number of iterations because fabs (delta_frac - 1.0) 
does not get small enough, and NaN is returned.

The funny thing is that I could influence this behaviour in a way I 
didn't expect:

1. Linking to the current git version of GSL which was compiled with the 
default setting -O2 yielded:

(P - eP) [0] = nan, tolerance = 2.3283064365386963e-10
(P - eP) [1] = nan, tolerance = 2.3283064365386963e-10
(P - eP) [2] = nan, tolerance = 2.3283064365386963e-10
(P - eP) [3] = 4.1520314247533996e-05, tolerance = 2.3283064365386963e-10

2. Linking to GSL which was compiled without optimisation (-O0):

(P - eP) [0] = nan, tolerance = 2.3283064365386963e-10
(P - eP) [1] = nan, tolerance = 2.3283064365386963e-10
(P - eP) [2] = 1.8060553053089734e-11, tolerance = 2.3283064365386963e-10
(P - eP) [3] = nan, tolerance = 2.3283064365386963e-10

3. Running it with Valgrind (independent of the optimisation settings):

(P - eP) [0] = 6.9538819147396680e-12, tolerance = 2.3283064365386963e-10
(P - eP) [1] = -7.3852035598065413e-12, tolerance = 2.3283064365386963e-10
(P - eP) [2] = -1.1277645484142340e-12, tolerance = 2.3283064365386963e-10
(P - eP) [3] = 4.1520320904209207e-05, tolerance = 2.3283064365386963e-10

Maybe these (at least for me) surprising observations could be used to 
isolate the cause of the bug. Does anyone have a good idea?

Regards,
Frank

P.S.: I compiled and ran the program on two systems I have access to and 
got exactly the same results:

1.
* Athlon 64 X2 in 32 bit mode
* Kubuntu 7.10
* gcc (GCC) 4.1.3 20070929 (prerelease) (Ubuntu 4.1.2-16ubuntu2))

2.
* Pentium III
* Opensuse 10.2
* gcc (GCC) 4.1.2 20061115 (prerelease) (SUSE Linux)