The following reply was made to PR bin/144306; it has been noted by GNATS.
From: "Steven G. Kargl" <ka...@troutmask.apl.washington.edu>
To: =?UTF-8?Q?Ulrich_Sp=C3=B6rlein?= <u...@spoerlein.net>
Cc: bug-follo...@freebsd.org
Subject: Re: bin/144306: [libm] [patch] Nasty bug in jn(3)
Date: Fri, 12 Nov 2010 10:01:54 -0800 (PST)
Ulrich Sp?rlein wrote:
> Hello Steve, I saw your updated patch for jnf(3) on the mailing list and
> was wondering what the correct values for your test case would look like?
> I'm asking because on Ubuntu I get slightly different values. Especially
> disturbing is, that for FreeBSD the jn/jnf values are identical in this
> case.
Apologies for a delayed response. It seems spamassasin decided
your email was junk mail and put your email into my junk mail
directory.
Bruce has already given you some answers towards
the correct values. You can also use MPFR to generate
a set of values for comparison.
Another way to look at the issue is to start near a
zero of j0f and compute say jnf(3,x) with x ranging
over a small set of values.
#include<stdio.h>
#include<math.h>
int
main(void)
{
float x;
x = 2.404820;
do {
x = nextafterf(x, 10.);
printf("%e %e %e\n", x, jnf(3,x), jn(3,(double)x));
} while(x < 2.404826);
return (0);
}
Without the patch, the 2nd column of the output is not
a smooth function, which is counter to the mathematical
properties of a bessel fcn. Note, the 3rd column is
jn(3,x) where I have applied the patch.
troutmask:kargl[202] cc -o z t.c -lm && ./z
2.404820e+00 1.989337e-01 1.989989e-01
2.404820e+00 1.977170e-01 1.989990e-01
2.404821e+00 1.997695e-01 1.989990e-01
2.404821e+00 1.977899e-01 1.989991e-01
2.404821e+00 2.007952e-01 1.989991e-01
2.404821e+00 1.986288e-01 1.989991e-01
2.404822e+00 1.970318e-01 1.989992e-01
2.404822e+00 1.999776e-01 1.989992e-01
2.404822e+00 1.976307e-01 1.989993e-01
2.404822e+00 2.017480e-01 1.989993e-01
2.404823e+00 1.987786e-01 1.989994e-01
2.404823e+00 1.961335e-01 1.989994e-01
2.404823e+00 1.999664e-01 1.989994e-01
2.404823e+00 2.053210e-01 1.989995e-01
(12 lines removed)
With the patch and with bde's correction for fabsf(), you get
troutmask:kargl[204] cc -o z t.c -lm && ./z
2.404820e+00 1.989989e-01 1.989989e-01
2.404820e+00 1.989990e-01 1.989990e-01
2.404821e+00 1.989990e-01 1.989990e-01
2.404821e+00 1.989991e-01 1.989991e-01
2.404821e+00 1.989991e-01 1.989991e-01
2.404821e+00 1.989991e-01 1.989991e-01
2.404822e+00 1.989992e-01 1.989992e-01
2.404822e+00 1.989992e-01 1.989992e-01
2.404822e+00 1.989993e-01 1.989993e-01
2.404822e+00 1.989994e-01 1.989993e-01
2.404823e+00 1.989994e-01 1.989994e-01
2.404823e+00 1.989994e-01 1.989994e-01
2.404823e+00 1.989995e-01 1.989994e-01
2.404823e+00 1.989995e-01 1.989995e-01
(12 lines removed)
>
> Ubuntu 10.04:
>
I could not find an on-line source code repository for Ubuntu.
I suspect that Ubuntu uses an algorithm for jn[f] based on the
code from fdlibm. This problem appears to date back to at
least 1995 or so.
Now to explain the problem and the solution. jn[f](n,x) for
certain ranges of x uses downward recursion to compute the
value of the function. The recursion sequence that is generated
is proportional to the actual desired value, so a normalization
step is taken. This normalization is j0[f](x) divided by
the the zeroth sequence member. As Bruce notes, near the
zeros of j0[f](x) the computed value can have giga-ULP inaccuracy.
I found for the 1st zero of j0f(x) only the leading decimal digit
is correct. The solution to the issue is fairly straight
forward. The zeros of j0(x) and j1(x) never coincide, so as
j0(x) approaches a zero, the normalization constant switches to
j1[f](x) divided by the 2nd sequence member. The expectation is
that j1[f](x) is an more accurately computed value.
--
Steve
http://troutmask.apl.washington.edu/~kargl/
_______________________________________________
freebsd-bugs@freebsd.org mailing list
http://lists.freebsd.org/mailman/listinfo/freebsd-bugs
To unsubscribe, send any mail to "freebsd-bugs-unsubscr...@freebsd.org"
