From owner-freebsd-standards@FreeBSD.ORG  Thu Feb 25 11:10:02 2010
Return-Path: <owner-freebsd-standards@FreeBSD.ORG>
Delivered-To: freebsd-standards@hub.freebsd.org
Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34])
	by hub.freebsd.org (Postfix) with ESMTP id 9DC53106564A
	for <freebsd-standards@hub.freebsd.org>;
	Thu, 25 Feb 2010 11:10:02 +0000 (UTC)
	(envelope-from gnats@FreeBSD.org)
Received: from freefall.freebsd.org (freefall.freebsd.org
	[IPv6:2001:4f8:fff6::28])
	by mx1.freebsd.org (Postfix) with ESMTP id 7293A8FC12
	for <freebsd-standards@hub.freebsd.org>;
	Thu, 25 Feb 2010 11:10:02 +0000 (UTC)
Received: from freefall.freebsd.org (localhost [127.0.0.1])
	by freefall.freebsd.org (8.14.3/8.14.3) with ESMTP id o1PBA2BA055243
	for <freebsd-standards@freefall.freebsd.org>;
	Thu, 25 Feb 2010 11:10:02 GMT
	(envelope-from gnats@freefall.freebsd.org)
Received: (from gnats@localhost)
	by freefall.freebsd.org (8.14.3/8.14.3/Submit) id o1PBA247055242;
	Thu, 25 Feb 2010 11:10:02 GMT (envelope-from gnats)
Date: Thu, 25 Feb 2010 11:10:02 GMT
Message-Id: <201002251110.o1PBA247055242@freefall.freebsd.org>
To: freebsd-standards@FreeBSD.org
From: Bruce Evans <brde@optusnet.com.au>
Cc: 
Subject: Re: standards/142803: j0 Bessel function inaccurate near zeros of
 the function
X-BeenThere: freebsd-standards@freebsd.org
X-Mailman-Version: 2.1.5
Precedence: list
Reply-To: Bruce Evans <brde@optusnet.com.au>
List-Id: Standards compliance <freebsd-standards.freebsd.org>
List-Unsubscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=unsubscribe>
List-Archive: <http://lists.freebsd.org/pipermail/freebsd-standards>
List-Post: <mailto:freebsd-standards@freebsd.org>
List-Help: <mailto:freebsd-standards-request@freebsd.org?subject=help>
List-Subscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=subscribe>
X-List-Received-Date: Thu, 25 Feb 2010 11:10:02 -0000

The following reply was made to PR standards/142803; it has been noted by GNATS.

From: Bruce Evans <brde@optusnet.com.au>
To: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu>
Cc: Bruce Evans <brde@optusnet.com.au>, FreeBSD-gnats-submit@freebsd.org,
        freebsd-standards@freebsd.org
Subject: Re: standards/142803: j0 Bessel function inaccurate near zeros of
 the function
Date: Thu, 25 Feb 2010 22:09:09 +1100 (EST)

 On Mon, 22 Feb 2010, Steven G. Kargl wrote:
 
 > Since the my real-life work requires additional accuracy for
 > the Bessel function routines in libm, I spent sometime playing
 > with e_j0f.  I've come up with the patch that follows my signature.
 > It only deals with the first 32 zeros of j0(x), others can be
 > added to the table.  I've yet to come up with a general method
 > to improve the accuracy around all zeros.   To help you gauge
 > the improvement, one can look at
 
 Interesting.  I will reply more later...  I just browsed Abromowitz
 (sp?) and Stegun and saw that it gives methods for finding all the zeros
 and says that there are complex zeros.  It seems strong on formulae and
 weak on practical computational methods and algorithms (even for its
 time; of course the best methods now are different).
 
 > +float
 > __ieee754_j0f(float x)
 > {
 > 	float z, s,c,ss,cc,r,u,v;
 > -	int32_t hx,ix;
 > +	int32_t hx,ix,k;
 >
 > 	GET_FLOAT_WORD(hx,x);
 > 	ix = hx&0x7fffffff;
 > 	if(ix>=0x7f800000) return one/(x*x);
 > 	x = fabsf(x);
 > 	if(ix >= 0x40000000) {	/* |x| >= 2.0 */
 > +
 > +		if (x < 100.) {
 > +			k = (int)(x * invpi - 0.75);
 > +			if (fabsf(x - (k + 0.75) * pi - offset[k]) < 0.065)
 
 These and many other constants should by float constants.  Otherwise you
 get slowness and negatively useful extra accuracy (no useful extra accuracy,
 due to floats elsewhere, but more complicated error analysis).
 
 This reduction is similar to the one for the first few and/or first
 few million zeros of trig functions in e_rem_pio2*.c.
 
 > +				return (j0zerof(k, x));
 > +		}
 > +
 > 		s = sinf(x);
 > 		c = cosf(x);
 > 		ss = s-c;
 >
 
 Bruce