standards/142803: j0 Bessel function inaccurate near
zeros of the function
Steven G. Kargl
kargl at troutmask.apl.washington.edu
Thu Jan 14 22:49:39 UTC 2010
Bruce Evans wrote:
> 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.
Yes, it is very hard, but apparently you do not use bessel
functions in your everyday life. :)
I only discover this issue because I need bessel functions
of complex arguments and I found my routines have issues
in the vicinity of zeros. So, I decided to look at the
libm routines.
> > 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.
I forgot to mention that 'my err' and 'libm err' are
in units of epsilon (ie, FLT_EPSILON for j0f).
> > 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.
>
(snip)
> 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?).
The 4000-bit setting for mpfr was a hold over from testing mpfr_j0
against my ascending series implementation of j0 with mpfr
primitives. As few as 128-bits is sufficient to achieve the
following:
2.404825 5.6434398E-08 5.9634296E-08 5.6434400E-08 0.05 152824.59
5.520078 2.4476657E-08 2.4153294E-08 2.4476659E-08 0.10 18878.52
8.653728 1.0355303E-07 1.0359805E-07 1.0355306E-07 0.86 1694.47
11.791534 -3.5291243E-09 -3.5193941E-09 -3.5301714E-09 75.93 781.53
14.930918 -6.4815082E-09 -6.3911618E-09 -6.4815052E-09 0.23 6722.88
18.071064 5.1532352E-09 5.3149818E-09 5.1532318E-09 0.23 10910.50
21.211637 -1.5023349E-07 -1.5002509E-07 -1.5023348E-07 2.70 56347.01
24.352472 1.2524569E-07 1.2516310E-07 1.2524570E-07 0.28 2834.53
27.493479 5.4331110E-08 5.4263626E-08 5.4331104E-08 0.29 3261.75
30.634607 1.2205545E-07 1.2203689E-07 1.2205546E-07 0.09 645.39
33.775822 -2.0213095E-07 -2.0206903E-07 -2.0213095E-07 0.27 6263.95
36.917099 8.4751576E-08 8.4749573E-08 8.4751581E-08 0.18 82.59
40.058426 -1.7484838E-08 -1.7475532E-08 -1.7484840E-08 0.12 767.56
43.199791 -9.2091398E-08 -9.2135146E-08 -9.2091406E-08 2.47 13530.51
46.341187 2.1663259E-07 2.1664336E-07 2.1663259E-07 0.16 268.90
49.482609 -1.2502527E-07 -1.2504157E-07 -1.2502526E-07 2.69 23512.60
52.624050 1.8706569E-07 1.8707487E-07 1.8706569E-07 0.01 251.43
55.765511 -2.0935557E-08 -2.0932896E-08 -2.0935556E-08 0.10 227.04
58.906982 1.5637660E-07 1.5634730E-07 1.5637661E-07 0.28 892.23
62.048470 3.5779891E-08 3.5787338E-08 3.5779899E-08 0.42 402.61
As I suspected by adding additional terms to the asymptotic
approximation and performing all computations with double
precision, reduces 'my err' (5th column). The value at
x=11.7... is the best I can get. The asymptotic approximations
contain divergent series and additional terms do not help.
--
Steve
http://troutmask.apl.washington.edu/~kargl/
More information about the freebsd-standards
mailing list