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

[Bruce Evans]

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