[Bug 583475] New: glibc/libm: Wrong result for Bessel function jn on x86-64
http://bugzilla.novell.com/show_bug.cgi?id=583475
http://bugzilla.novell.com/show_bug.cgi?id=583475#c0
Summary: glibc/libm: Wrong result for Bessel function jn on
x86-64
Classification: openSUSE
Product: openSUSE 11.3
Version: Factory
Platform: Other
OS/Version: Other
Status: NEW
Severity: Normal
Priority: P5 - None
Component: Basesystem
AssignedTo: bnc-team-screening@forge.provo.novell.com
ReportedBy: burnus@gmx.de
QAContact: qa@suse.de
Found By: ---
Blocker: ---
User-Agent: Mozilla/5.0 (X11; U; Linux x86_64; de; rv:1.9.2.0)
Gecko/20100115 SUSE/3.6.0-1.6 Firefox/3.6
The Bessel function "jn" gives -INF on x86-64; the 32bit version and the
float/long double version work OK, as the example shows below. As I am more
familiar with Fortran, I only included the float/long double version as Fortran
version.
Note: The compile-time result using MPFR is OK; the long double version is also
fine. However, the float version is a bit imprecise - especially on x86-64.
$ cat bessel.c
#include
034212501569217985245804444265895
Reproducible: Always -- Configure bugmail: http://bugzilla.novell.com/userprefs.cgi?tab=email ------- You are receiving this mail because: ------- You are on the CC list for the bug.
http://bugzilla.novell.com/show_bug.cgi?id=583475
http://bugzilla.novell.com/show_bug.cgi?id=583475#c
yang xiaoyu
http://bugzilla.novell.com/show_bug.cgi?id=583475
http://bugzilla.novell.com/show_bug.cgi?id=583475#c1
--- Comment #1 from Tobias Burnus
There appears to be a really nasty bug in jn(3) from fdlibm, which is the foundation for most libm implementations (e.g., FreeBSD's libm).
2.4048255576957729_8 is the first zero of the zeroth order cylindrical Bessel function. jn(3) uses j0(3) to normalize its recursively generated result.
There are other problems with the Bessel function routines in fdlibm; namely, an extreme loss of precision at values near a zero.
If you libm is based on fdlibm, this patch fixes the problem with jn(3).
[Note: I have not checked whether this applies to GLIBC but it seems to be plausible that GLIBC does something similar. The rather bad result for "jnf" is another issue, which should also be investigated.] Index: e_jn.c =================================================================== --- e_jn.c (revision 204219) +++ e_jn.c (working copy) @@ -200,7 +200,12 @@ } } } - b = (t*__ieee754_j0(x)/b); + z = __ieee754_j0(x); + w = __ieee754_j1(x); + if (fabs(z) >= fabs(w)) + b = (t*z/b); + else + b = (t*w/a); } } if(sgn==1) return -b; else return b; -- Configure bugmail: http://bugzilla.novell.com/userprefs.cgi?tab=email ------- You are receiving this mail because: ------- You are on the CC list for the bug.
http://bugzilla.novell.com/show_bug.cgi?id=583475
http://bugzilla.novell.com/show_bug.cgi?id=583475#c2
--- Comment #2 from Tobias Burnus
The simple explanation is that jn(n,x) uses a downward recursion to compute a sequence of values f(0,x), f(1,x), f(3,x), ... f(n,x). These values are proportional to the actual Bessel function values, so one normalizes the the value f(3,x) by f(3,x) * j0(x) / f(0,x). Well, for 'x = a zero of j0(x)', you end up with 'f(3,x) * nonzero value / 0' because j0(x) has an enormous loss of precision at and near a zero.
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.
My patch works because it chooses the normalization based on whether |j0(x)| >= |j1(x)| and we know that the zeros of j0(x) and j1(x) never coincide.
If you fixed it, can you push the patch upstream? Thanks! -- Configure bugmail: http://bugzilla.novell.com/userprefs.cgi?tab=email ------- You are receiving this mail because: ------- You are on the CC list for the bug.
http://bugzilla.novell.com/show_bug.cgi?id=583475
http://bugzilla.novell.com/show_bug.cgi?id=583475#c3
Petr Baudis
http://bugzilla.novell.com/show_bug.cgi?id=583475
http://bugzilla.novell.com/show_bug.cgi?id=583475#c4
Tobias Burnus
is Steve planning to push this upstream?
To my knowledge he is only working on FreeBSD's libm and not on GLIBC or the other various libm versions, which are at the end seemingly based on SUN's fdlibm version (available at http://www.netlib.org/fdlibm/e_jn.c). Thus,I think the answer is: No for glibc upstream. -- Configure bugmail: http://bugzilla.novell.com/userprefs.cgi?tab=email ------- You are receiving this mail because: ------- You are on the CC list for the bug.
http://bugzilla.novell.com/show_bug.cgi?id=583475
http://bugzilla.novell.com/show_bug.cgi?id=583475#c5
Petr Baudis
participants (1)
-
bugzilla_noreply＠novell.com