From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 20175 invoked by alias); 14 Feb 2008 23:01:31 -0000 Received: (qmail 20167 invoked by uid 22791); 14 Feb 2008 23:01:30 -0000 X-Spam-Check-By: sourceware.org Received: from nf-out-0910.google.com (HELO nf-out-0910.google.com) (64.233.182.186) by sourceware.org (qpsmtpd/0.31) with ESMTP; Thu, 14 Feb 2008 23:01:00 +0000 Received: by nf-out-0910.google.com with SMTP id b11so335408nfh.48 for ; Thu, 14 Feb 2008 15:00:57 -0800 (PST) Received: by 10.86.96.18 with SMTP id t18mr1807048fgb.13.1203030057635; Thu, 14 Feb 2008 15:00:57 -0800 (PST) Received: from ?192.168.178.21? ( [87.189.19.81]) by mx.google.com with ESMTPS id p38sm6359182fke.13.2008.02.14.15.00.55 (version=TLSv1/SSLv3 cipher=RC4-MD5); Thu, 14 Feb 2008 15:00:56 -0800 (PST) Message-ID: <47B4C827.1040004@googlemail.com> Date: Thu, 14 Feb 2008 23:01:00 -0000 From: Frank Reininghaus User-Agent: Thunderbird 2.0.0.6 (X11/20071022) MIME-Version: 1.0 To: gsl-discuss@sourceware.org Subject: Examining bug #21838 Content-Type: multipart/mixed; boundary="------------070807060801090909070801" 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/msg00031.txt.bz2 This is a multi-part message in MIME format. --------------070807060801090909070801 Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit Content-length: 2107 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) --------------070807060801090909070801 Content-Type: text/x-csrc; name="testcase-bug-21838.c" Content-Transfer-Encoding: 7bit Content-Disposition: inline; filename="testcase-bug-21838.c" Content-length: 843 #include #include #include /* 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; } --------------070807060801090909070801--