lib/44170: Bug in jn and jnf.

[henning.petersen%t-online.de@localhost]

>Number:         44170
>Category:       lib
>Synopsis:       Bug in jn and jnf.
>Confidential:   no
>Severity:       serious
>Priority:       medium
>Responsible:    lib-bug-people
>State:          open
>Class:          sw-bug
>Submitter-Id:   net
>Arrival-Date:   Mon Nov 29 10:20:01 +0000 2010
>Originator:     Henning Petersen
>Release:        NetBSD-current
>Organization:
>Environment:
>Description:
With program

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

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

z = 2.4048255576957729;
for(n = 2; n < 10; n++)
{
printf("%d %e\n", n, jn(n,z));
}
return 0;
}
i have following errors.

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

>How-To-Repeat:

>Fix:
diff -u -r1.13 e_jn.c
--- lib/libm/src/e_jn.c 20 Aug 2007 16:01:39 -0000      1.13
+++ lib/libm/src/e_jn.c 29 Nov 2010 07:27:52 -0000
@@ -203,7 +203,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;
diff -u -r1.10 e_jnf.c
--- lib/libm/src/e_jnf.c        19 Jan 2009 05:58:27 -0000      1.10
+++ lib/libm/src/e_jnf.c        29 Nov 2010 07:27:52 -0000
@@ -157,7 +157,12 @@
                        }
                    }
                }
-               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;