From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 7752 invoked by alias); 15 Feb 2008 09:07:10 -0000 Received: (qmail 9148 invoked by uid 22791); 15 Feb 2008 00:31:28 -0000 X-Spam-Check-By: sourceware.org Date: Fri, 15 Feb 2008 09:07:00 -0000 From: Jason Stover To: Frank Reininghaus Cc: gsl-discuss@sourceware.org Subject: Re: Examining bug #21838 Message-ID: <20080215003043.GA7666@SDF.LONESTAR.ORG> Reply-To: jason@sakla.net References: <47B4C827.1040004@googlemail.com> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <47B4C827.1040004@googlemail.com> User-Agent: Mutt/1.4.2.1i Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org X-SW-Source: 2008-q1/txt/msg00032.txt.bz2 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