From owner-freebsd-bugs@FreeBSD.ORG Fri Nov 12 18:10:12 2010 Return-Path: Delivered-To: freebsd-bugs@hub.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 4E40B106564A for ; Fri, 12 Nov 2010 18:10:12 +0000 (UTC) (envelope-from gnats@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [IPv6:2001:4f8:fff6::28]) by mx1.freebsd.org (Postfix) with ESMTP id 20D1D8FC24 for ; Fri, 12 Nov 2010 18:10:12 +0000 (UTC) Received: from freefall.freebsd.org (localhost [127.0.0.1]) by freefall.freebsd.org (8.14.4/8.14.4) with ESMTP id oACIACPD016097 for ; Fri, 12 Nov 2010 18:10:12 GMT (envelope-from gnats@freefall.freebsd.org) Received: (from gnats@localhost) by freefall.freebsd.org (8.14.4/8.14.4/Submit) id oACIABor016096; Fri, 12 Nov 2010 18:10:11 GMT (envelope-from gnats) Date: Fri, 12 Nov 2010 18:10:11 GMT Message-Id: <201011121810.oACIABor016096@freefall.freebsd.org> To: freebsd-bugs@FreeBSD.org From: "Steven G. Kargl" Cc: Subject: Re: bin/144306: [libm] [patch] Nasty bug in jn(3) X-BeenThere: freebsd-bugs@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list Reply-To: "Steven G. Kargl" List-Id: Bug reports List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Fri, 12 Nov 2010 18:10:12 -0000 The following reply was made to PR bin/144306; it has been noted by GNATS. From: "Steven G. Kargl" To: =?UTF-8?Q?Ulrich_Sp=C3=B6rlein?= Cc: bug-followup@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 #include 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/