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.