[Bug math/11589] New: jn() Bessel function is massively imprecise at some points

["pasky at suse dot cz" <sourceware-bugzilla at sourceware dot org>]

(Steve Kargl and Tobias Burnus outlined the problem and proposed the method of
fixing it. See also https://bugzilla.novell.com/show_bug.cgi?id=583475.)

There appears to be a really nasty bug in jn() from fdlibm, which
is the foundation for most libm implementations (including glibc libm).
The zeroth-order j0() and first-order j1() cylindrical Bessel functions are used
to recursively generate the jn() value, but only the zeroth-order Bessel
function is used to normalize it; however, each of the functions gets highly
imprecise (approaching "bogus") near its zero point, making the jn() value
itself bogus.

But in fact, the zero points of j0() and j1() never coincide, thus j1() should
be used in case it is more precise than j0(). (That is, simply when its value is
further from zero.)

As an example, 2.4048255576957729_8 is the first zero of j0():

#include <stdio.h>
#include <math.h>

int
main ()
{
  double r8 = 2.4048255576957729;
  printf("%f\n", jn (3, r8));
  return 0;
}

The proper value as calculated by Mathematica is 0.19899990535769... However,
jn() returns -inf on 64-bit arch, or 0.185007 on 32-bit arch. With the proposed
patch below, the returned value is 0.199000.

Steve Kargl went on to measure the systematic error of j0():

> http://troutmask.apl.washington.edu/~kargl/j0f.jpg
> http://troutmask.apl.washington.edu/~kargl/j0f_1.jpg
>
> These are the error in j0f(x) near x = 2.404... measured in
> ULP. I've marched through this zero using nextafterf(), and
> you get about 9 million ULP in error, ie., 1 significant 
> decimal digit.

-- 
           Summary: jn() Bessel function is massively imprecise at some
                    points
           Product: glibc
           Version: 2.12
            Status: NEW
          Severity: normal
          Priority: P2
         Component: math
        AssignedTo: aj at suse dot de
        ReportedBy: pasky at suse dot cz
                CC: glibc-bugs at sources dot redhat dot com


http://sourceware.org/bugzilla/show_bug.cgi?id=11589

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.