From 1248c1c41508387ff282b208636737e8cdc9b5b0 Mon Sep 17 00:00:00 2001 From: Petr Baudis Date: Fri, 9 Sep 2011 22:16:10 -0400 Subject: Fix jn precision --- sysdeps/ieee754/ldbl-96/e_jnl.c | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) (limited to 'sysdeps/ieee754/ldbl-96/e_jnl.c') diff --git a/sysdeps/ieee754/ldbl-96/e_jnl.c b/sysdeps/ieee754/ldbl-96/e_jnl.c index 3d715d36aa..bedff7d566 100644 --- a/sysdeps/ieee754/ldbl-96/e_jnl.c +++ b/sysdeps/ieee754/ldbl-96/e_jnl.c @@ -281,7 +281,16 @@ __ieee754_jnl (n, x) } } } - b = (t * __ieee754_j0l (x) / b); + /* j0() and j1() suffer enormous loss of precision at and + * near zero; however, we know that their zero points never + * coincide, so just choose the one further away from zero. + */ + z = __ieee754_j0l (x); + w = __ieee754_j1l (x); + if (fabsl (z) >= fabsl (w)) + b = (t * z / b); + else + b = (t * w / a); } } if (sgn == 1) -- cgit 1.4.1