Re: Examining bug #21838

[Jason Stover <jstover at sdf dot lonestar dot org>]

On Fri, Feb 15, 2008 at 12:00:55AM +0100, Frank Reininghaus wrote:
> 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).

The test cases have huge numbers of denominator degrees of freedom
(\nu_2 in the manual).  As \nu_2 --> infinity, the F density tends to
a gamma density. So it might be worthwhile to fix this by calling
gsl_cdf_gamma_[PQ] instead of beta_inc_AXPY() in such cases.

> 
> 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)



-- 
jstover@sdf.lonestar.org
SDF Public Access UNIX System - http://sdf.lonestar.org