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

Thu Jan 14 22:50:03 UTC 2010

The following reply was made to PR standards/142803; it has been noted by GNATS.

From: "Steven G. Kargl" <kargl at troutmask.apl.washington.edu>
To: Bruce Evans <brde at optusnet.com.au>
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 14:49:45 -0800 (PST)

 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/