Date: Thu, 25 Feb 2010 22:09:09 +1100 (EST) From: Bruce Evans <brde@optusnet.com.au> To: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu> Cc: FreeBSD-gnats-submit@freebsd.org, freebsd-standards@freebsd.org Subject: Re: standards/142803: j0 Bessel function inaccurate near zeros of the function Message-ID: <20100225214825.E1337@besplex.bde.org> In-Reply-To: <201002230021.o1N0LLNb002353@troutmask.apl.washington.edu> References: <201002230021.o1N0LLNb002353@troutmask.apl.washington.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20100225214825.E1337>