[Bug math/12292] New: Functions jn() and jnf() are not correctly computed.

["carlos at codesourcery dot com" <sourceware-bugzilla at sourceware dot org>]

http://sourceware.org/bugzilla/show_bug.cgi?id=12292

           Summary: Functions jn() and jnf() are not correctly computed.
           Product: glibc
           Version: 2.13
            Status: NEW
          Severity: normal
          Priority: P2
         Component: math
        AssignedTo: aj@suse.de
        ReportedBy: carlos@codesourcery.com
              Host: i686-pc-linux-gnu
            Target: i686-pc-linux-gnu
             Build: i686-pc-linux-gnu


Functions jn() and jnf() are not correctly computed near the zero for a bessel
function of the first kind of order 3.

The computation error is a result of a bug in the code for these two files:
sysdeps/ieee754/dbl-64/e_jn.c 
sysdeps/ieee754/flt-32/e_jnf.c.

Recently a bug was discovered in the FreeBSD libc:
http://www.freebsd.org/cgi/query-pr.cgi?pr=bin/144306
~~~
Log:
Fix bug in jn(3) and jnf(3) that led to -inf results

Explanation by Steve:
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 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 a more accurately computed value.

PR: bin/144306
Submitted by: Steven G. Kargl <kargl@troutmask.apl.washington.edu>
Reviewed by: bde
MFC after: 7 days

Modified:
head/lib/msun/src/e_jn.c
head/lib/msun/src/e_jnf.c
~~~

I have verified that the bug is present in glibc git head, and independently
verified using Octave 3.1 to compute the correct bessel function results.

The test case:
~~~
#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);
}
~~~
Should produce a smooth function near the zero of the bessel function of the
first kind with order 3.

It does not produce a smooth function, and in many steps it produces an
incorrect value (verified using Octave 3.1).

The FreeBSD libc patch:
~~~
--- head/lib/msun/src/e_jn.c     Sat Nov 13 10:38:06 2010        (r215236)
+++ head/lib/msun/src/e_jn.c     Sat Nov 13 10:54:10 2010        (r215237)
@@ -200,7 +200,12 @@ __ieee754_jn(int n, double x)
                         }
                     }
                 }
-                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;
Modified: head/lib/msun/src/e_jnf.c

==============================================================================
--- head/lib/msun/src/e_jnf.c    Sat Nov 13 10:38:06 2010        (r215236)
+++ head/lib/msun/src/e_jnf.c    Sat Nov 13 10:54:10 2010        (r215237)
@@ -152,7 +152,12 @@ __ieee754_jnf(int n, float x)
                         }
                     }
                 }
-                b = (t*__ieee754_j0f(x)/b);
+                z = __ieee754_j0f(x);
+                w = __ieee754_j1f(x);
+                if (fabsf(z) >= fabsf(w))
+                    b = (t*z/b);
+                else
+                    b = (t*w/a);
             }
         }
         if(sgn==1) return -b; else return b;
~~~
Should also be applied to:
sysdeps/ieee754/dbl-64/e_jn.c 
sysdeps/ieee754/flt-32/e_jnf.c.

Thank you.

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.