bin/144306: [libm] [patch] Nasty bug in jn(3)

[Steven G. Kargl]

The following reply was made to PR bin/144306; it has been noted by GNATS.

From: "Steven G. Kargl" <kargl at troutmask.apl.washington.edu>
To: =?UTF-8?Q?Ulrich_Sp=C3=B6rlein?= <uqs at spoerlein.net>
Cc: bug-followup at 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/