From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 12208 invoked by alias); 17 Aug 2010 20:38:42 -0000 Received: (qmail 12200 invoked by uid 22791); 17 Aug 2010 20:38:41 -0000 X-SWARE-Spam-Status: No, hits=0.2 required=5.0 tests=AWL,BAYES_20,TW_SQ,T_RP_MATCHES_RCVD X-Spam-Check-By: sourceware.org Received: from fencepost.gnu.org (HELO fencepost.gnu.org) (140.186.70.10) by sourceware.org (qpsmtpd/0.43rc1) with ESMTP; Tue, 17 Aug 2010 20:38:37 +0000 Received: from localhost ([127.0.0.1]:52810 helo=fencepost.gnu.org) by fencepost.gnu.org with esmtp (Exim 4.69) (envelope-from ) id 1OlSva-00069S-Tb; Tue, 17 Aug 2010 16:38:34 -0400 Date: Tue, 17 Aug 2010 20:38:00 -0000 Message-ID: <43mxsl6klg.wl%bjg@gnu.org> From: Brian Gough To: Thanh Vo Cc: gsl-discuss@sourceware.org, bug-gsl@gnu.org, "Earl T. Barr" Subject: Re: [Bug-gsl] Potential Bugs? In-Reply-To: References: User-Agent: Wanderlust/2.15.6 (Almost Unreal) Emacs/23.2 Mule/6.0 (HANACHIRUSATO) MIME-Version: 1.0 (generated by SEMI 1.14.6 - "Maruoka") Content-Type: text/plain; charset=US-ASCII 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: 2010-q3/txt/msg00003.txt.bz2 At Wed, 11 Aug 2010 17:14:34 -0700, Thanh Vo wrote: > My name's Thanh Vo. My research group is developing a tool to detect > floating point exceptions in numerical code. We've found inputs that > cause some of > the GSL 1.14 functions to throw floating-point exceptions. I've written > this post to see if you consider any of them to be real bugs. Hi, Thanks for your email. I have made some comments below. This sounds like a useful tool, I hope you will be able to release it as free software. > > 1. The function airy_aie throws a divide by zero exception when its > input x = 0 at line 618 in file gsl/specfunc/airy.c : double z = > 2.0/(x*sqx) - 1.0;. It throws an invalid exception when x = -1 at > line 617: double sqx = sqrt(x);. > > 2. The function airy_bie throws a divide by zero exception when its > input x = 0 at line 635 in file gsl/specfunc/airy.c: double z = > ATR/(x*sqx) + BTR;. It throws an invalid exception when x=-1 at line > 634: double sqx = sqrt(x);. These would be bugs but the functions are static and only called with x > 1 or x > 2 so the exceptions cannot occur in practice. > 3. The function gsl_sf_asymp_thetanu_corr_e throws a divide by zero > exception when its input x = 0 and nu != 0 at the line 181 in file > gsl/specfunc/bessel_amp_phase.c: const double r = 2.0*nu/x;. It > throws an invalid exception when x=0 and nu=0 at line 181: const > double r = 2.0*nu/x;. > > 4. The function gsl_sf_bessel_asymp_Mnu_e throws a divide by zero > exception when its input x=0 and nu !=0 at the line 167 in file > gsl/specfunc/bessel_amp_phase.c: const double r = 2.0*nu/x;. It > throws an invalid exception when x=0 and nu=0 at line 167: const > double r = 2.0*nu/x;. It's a valid warning. The functions are publicly accessible, but not documented. Within the library they are only called with large values of x >= 1000 so the exception does not occur in practice. Ideally they should be static functions. > 5. The function gsl_sf_bessel_J0 throws a overflow exception when its > input x = 5.85177e+145 at the line 88 in the file > gsl/specfunc/bessel_J0.c: const double z = 32.0/(y*y) - 1.0;. It > throws an underflow exception when x = 6.8481e-155 at the line 81: > result->err = y*y;. The overflow exception is a valid warning. It's also a bug in the sense that it is avoidable (although the end result is correct, z=-1 since 1/inf=0). It would be better written as (32/y)/y. For the underflow it's a valid warning. There isn't really any alternative to returning zero in that case though. -- Brian Gough