[Numeric-interest] fdlibm-comments from arundel at freebsd.org: IEEE 754 bug in fdlibm 5.3

[David Hough 754R work]

Subject: 	IEEE 754 bug in fdlibm 5.3
Date: 	Thu, 16 Sep 2010 20:00:27 +0000
From: 	Alexander Best <arundel at freebsd.org>
To: 	fdlibm-comments at sun.com



Hi there,

we have a bug report concering fdlibm in FreeBSD. The bug can be triggered
this way:

1) Compile the following code with `cc -o z testjn.c -lm`:

%

#include<stdio.h>
#include<math.h>

int
main(void)
{
double z;
int n;

z = 2.4048255576957729;
for (n = 2; n<  10; n++)
printf("%d %e\n", n, jn(n,z));
return (0);
}

%

2) Run it and witness the following output:

2 4.317548e-01
3 -inf
4 4.069027e-02
5 -inf
6 3.247664e-03
7 -inf
8 7.495602e-05
9 -inf

3) The problem: The value of a integer order Bessel function with argument
z = 2.4048255576957729 is never -inf.

4) The fix:

%

Index: e_jn.c
===================================================================
--- e_jn.c       (revision 204219)
+++ e_jn.c       (working copy)
@@ -200,7 +200,12 @@
                          }
                      }
                  }
-                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;

%

Please note that the report and the patch have been sent in by Steven G. Kargl.
The bug report itself is available under
http://www.freebsd.org/cgi/query-pr.cgi?pr=144306

cheers.
alex

-- 
a13x