[src/trunk]: src/lib/libm/src fix accuracy problems in argument ranges where ...

[drochner <drochner%NetBSD.org@localhost>]

details:   https://anonhg.NetBSD.org/src/rev/1b6c76c71386
branches:  trunk
changeset: 759117:1b6c76c71386
user:      drochner <drochner%NetBSD.org@localhost>
date:      Mon Nov 29 15:10:06 2010 +0000

description:
fix accuracy problems in argument ranges where j0(x) is small, closes
PR lib/44170 by Henning Petersen
(originally from Steven G. Kargl per FreeBSD PR bin/144306)

diffstat:

 lib/libm/src/e_jn.c  |  9 +++++++--
 lib/libm/src/e_jnf.c |  9 +++++++--
 2 files changed, 14 insertions(+), 4 deletions(-)

diffs (52 lines):

diff -r 1dc7cfdcb23e -r 1b6c76c71386 lib/libm/src/e_jn.c
--- a/lib/libm/src/e_jn.c       Mon Nov 29 14:52:11 2010 +0000
+++ b/lib/libm/src/e_jn.c       Mon Nov 29 15:10:06 2010 +0000
@@ -12,7 +12,7 @@
 
 #include <sys/cdefs.h>
 #if defined(LIBM_SCCS) && !defined(lint)
-__RCSID("$NetBSD: e_jn.c,v 1.13 2007/08/20 16:01:39 drochner Exp $");
+__RCSID("$NetBSD: e_jn.c,v 1.14 2010/11/29 15:10:06 drochner Exp $");
 #endif
 
 /*
@@ -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 -r 1dc7cfdcb23e -r 1b6c76c71386 lib/libm/src/e_jnf.c
--- a/lib/libm/src/e_jnf.c      Mon Nov 29 14:52:11 2010 +0000
+++ b/lib/libm/src/e_jnf.c      Mon Nov 29 15:10:06 2010 +0000
@@ -15,7 +15,7 @@
 
 #include <sys/cdefs.h>
 #if defined(LIBM_SCCS) && !defined(lint)
-__RCSID("$NetBSD: e_jnf.c,v 1.10 2009/01/19 05:58:27 lukem Exp $");
+__RCSID("$NetBSD: e_jnf.c,v 1.11 2010/11/29 15:10:06 drochner Exp $");
 #endif
 
 #include "math.h"
@@ -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;