Examining bug #21838

[Frank Reininghaus <frank78ac at googlemail dot com>]

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)

#include <gsl/gsl_cdf.h>
#include <stdio.h>

#include <gsl/gsl_math.h>

/* tolerance taken from cdf/test.c */
#define TEST_TOL6  (1048576.0*GSL_DBL_EPSILON)

int main () {
  double P [4];  /* P values calculated using gsl_cdf_fdist_P () */
  double eP [4]; /* expected values taken from cdf/test.c */
  
  P [0] = gsl_cdf_fdist_P (3.479082213465832574, 1, 4040712);
  eP [0] = 0.93785072763723411967;

  P [1] = gsl_cdf_fdist_P (3.002774644786533109, 1, 4040712);
  eP [1] = 0.91687787379476055771;

  P [2] = gsl_cdf_fdist_P (3.000854441173130827, 1, 4040712);
  eP [2] = 0.91677930719813578619;

  P [3] = gsl_cdf_fdist_P (3.000064021622133037, 1, 4040712);
  eP [3] = 0.91678021757407962678;

  int i;
  
  for (i = 0; i < 4; i++)
    printf ("(P - eP) [%i] = %.16e, tolerance = %.16e\n",
	    i, eP [i] - P [i], TEST_TOL6);

  return 0;
}