standards/142803: j0 Bessel function inaccurate near zeros of the function

[Bruce Evans]

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

From: Bruce Evans <brde at optusnet.com.au>
To: "Steven G. Kargl" <kargl at troutmask.apl.washington.edu>
Cc: Bruce Evans <brde at optusnet.com.au>, FreeBSD-gnats-submit at freebsd.org,
        freebsd-standards at 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