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

[Steven G. Kargl]

>Number:         142803
>Category:       standards
>Synopsis:       j0 Bessel function inaccurate near zeros of the function
>Confidential:   no
>Severity:       serious
>Priority:       low
>Responsible:    freebsd-standards
>State:          open
>Quarter:        
>Keywords:       
>Date-Required:
>Class:          sw-bug
>Submitter-Id:   current-users
>Arrival-Date:   Thu Jan 14 01:50:01 UTC 2010
>Closed-Date:
>Last-Modified:
>Originator:     Steven G. Kargl
>Release:        FreeBSD 9.0-CURRENT amd64
>Organization:
APL-UW
>Environment:
System: FreeBSD troutmask.apl.washington.edu 9.0-CURRENT FreeBSD 9.0-CURRENT #3 r201269M: Mon Jan 4 13:52:03 PST 2010 kargl at troutmask.apl.washington.edu:/usr/obj/usr/src/sys/SPEW amd64


	
>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. 

    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
14.930918 -6.6019510E-09 -6.3911618E-09 -6.4815052E-09   8962.94   6722.88
18.071064  5.1734359E-09  5.3149818E-09  5.1532318E-09   1362.83  10910.50
21.211637 -1.5023797E-07 -1.5002509E-07 -1.5023348E-07   1213.07  56347.01
24.352472  1.2524691E-07  1.2516310E-07  1.2524570E-07     41.65   2834.53
27.493479  5.4330716E-08  5.4263626E-08  5.4331104E-08     18.77   3261.75
30.634607  1.2205560E-07  1.2203689E-07  1.2205546E-07      4.85    645.39
33.775822 -2.0213101E-07 -2.0206903E-07 -2.0213095E-07      5.48   6263.95
36.917099  8.4751605E-08  8.4749573E-08  8.4751581E-08      0.99     82.59
40.058426 -1.7484851E-08 -1.7475532E-08 -1.7484840E-08      0.91    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.0935554E-08 -2.0932896E-08 -2.0935556E-08      0.19    227.03
58.906982  1.5637660E-07  1.5634730E-07  1.5637661E-07      0.26    892.22
62.048470  3.5779891E-08  3.5787338E-08  3.5779888E-08      0.14    403.17

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.  

>How-To-Repeat:

>Fix:

I'll get to it eventually.


>Release-Note:
>Audit-Trail:
>Unformatted: