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

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/