Date: Thu, 14 Jan 2010 21:26:49 +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: <20100114200057.C62558@delplex.bde.org> In-Reply-To: <201001140143.o0E1hOm8040083@troutmask.apl.washington.edu> References: <201001140143.o0E1hOm8040083@troutmask.apl.washington.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
On Wed, 13 Jan 2010, Steven G. Kargl wrote: >> Description: > > The j0 bessel function supplied by libm is fairly inaccurate at > arguments at and near a zero of the function. Here's the first > 30 zeros computed by j0f, my implementation of j0f, a 4000-bit > significand computed via MPFR (the assumed exact value), and the > relative absolute error. This is a very hard and relatively unimportant problem. The corresponding problem for trigonometric functions is fairly hard and relatively very important. It is solved in fdlibm (and this in FreeBSD libm) essentially by using a table of all the relevant zeros, with the necessary several thousands of binary digits of accuracy for many of the zeros. Since trigonometric functions are periodic, it is not necessary to have a table of approximations (e.g., polynomials) for each of the zeros of interest. There are O(2**100) zeros of interest for 128-bit long doubles, so simple tables wouldn't work even for the zeros. The corresponding problem for lgamma() is relatively simple and relatively important. It is solved in places like the Intel ia64 libm by using a table of all the relevant zeros and a table of approximations for each of the zeros. There are only about 70 relevant zeros (since zeroes less than about -70 are so close to the (pseudo-)poles that they cannot be represented in double precision (unlike the poles which, since they are at the negative integers, can be represented down to about -10**53 in double precision)). For the corresponding problem for at least the simplest bessel function j0(), the zeros are distributed like the zeros of sin(), but they aren't exactly at multiples of Pi like those of sin(), so just finding all the relevant ones to the necessary thousands of binary digits of accuracy is a hard problem. At least j0() is bounded like sin() (and unlike lgamma()); thus we know that there are O(2**100) relevant zeros, so it is impossible to find them all in advance. Since bessel functions aren't periodic, it is also unclear how they could be calculated accurately near the zeros without using a different special approximation near every zero. In general, such calculations need to be done in extra precision even for the linear term, giving relatively minor extra complications. I once tried to caclulate lgamma(z) near just the first zero of lgamma(), using only extra precision applied to standard formalas (Stirling approximation and not reflection since extra precision was only easy to obtain than for the former), to see how far I could get using only extra precision. It took approximately doubled double precision (double precision with most intermediate calculations done in extra precision) to get near float precision for lgamma(). It might be worth making bessel functions accurate near just the first few zeros. > x my j0f(x) libm j0f(x) MPFR j0 my err libm err > 2.404825 5.6434434E-08 5.9634296E-08 5.6434400E-08 1.64 152824.59 > 5.520078 2.4476664E-08 2.4153294E-08 2.4476659E-08 0.31 18878.52 > 8.653728 1.0355323E-07 1.0359805E-07 1.0355306E-07 6.36 1694.47 > 11.791534 -2.4511966E-09 -3.5193941E-09 -3.5301714E-09 78243.14 781.53 Hmm. > ... > 62.048470 3.5779891E-08 3.5787338E-08 3.5779888E-08 0.14 403.17 The largest errors across all 2**32 float values for the 1-parameter float precision bessel functions in libm is about 3.7 Mulps for j0() and about 17 Gulps for y1(): %%% j0: max_er = 0x6f8686c5c9f26 3654467.3863, avg_er = 0.418, #>=1:0.5 = 671562746:1332944948 j1: max_er = 0x449026e9c286f 2246675.4566, avg_er = 0.425, #>=1:0.5 = 674510414:1347770568 lgamma:max_er = 0xe4cf242400000 7497618.0703, avg_er = 0.274, #>=1:0.5 = 70033452:508702222 y0: max_er = 0x93a2340c00000 4837658.0234, avg_er = 0.380, #>=1:0.5 = 594207097:1303826516 y1: max_er = 0x7ffa2068256f72ab 17176789825.1699, avg_er = 5.188, #>=1:0.5 = 459137173:1213136103 %%% > Note, my j0f(x) currently uses double precision to accumulate intermediate > values. Below x = 10, I use the ascending series to compute the value. > Above x = 10, I'm using an asymptotic approximation. I haven't investigated > whether additional terms in an asymptotic approximation would pull 'my err' > for x = 11, 14, 18, and 21 closer to the exact value. Stirling asymptotic is what barely worked (with significant extra precision) for the first zero of lgamma(). lgamma() is especially bad for higher zeros of lgamma(), since there is a pole very near the zero. I guess the behaviour for j0() is not so bad since there are no poles, but the asymptotic method also seems to work not so well (compared with a power series method) on lgamma(x) for x > 1 where there are no poles. Anyway, if you can get anywhere near < 10 ulp error near all zeros using only an asymptotic method, then that would be good. Then the asymptotic method would also be capable of locating the zeros very accurately. But I would be very surprised if this worked. I know of nothing similar for reducing mod Pi for trigonometric functions, which seems a simpler problem. I would expect it to at best involve thousands of binary digits in the tables for the asymptotic method, and corresponding thousands of digits of precision in the calculation (4000 as for mfpr enough for the 2**100th zero?). Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20100114200057.C62558>