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

Bruce Evans brde at optusnet.com.au
Thu Jan 14 10:30:05 UTC 2010


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: 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, 14 Jan 2010 21:26:49 +1100 (EST)

 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


More information about the freebsd-standards mailing list