From 1248c1c41508387ff282b208636737e8cdc9b5b0 Mon Sep 17 00:00:00 2001
From: Petr Baudis <pasky@suse.cz>
Date: Fri, 9 Sep 2011 22:16:10 -0400
Subject: Fix jn precision

---
 sysdeps/ieee754/flt-32/e_jnf.c | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

(limited to 'sysdeps/ieee754/flt-32')

diff --git a/sysdeps/ieee754/flt-32/e_jnf.c b/sysdeps/ieee754/flt-32/e_jnf.c
index de2e53d..dd3d551 100644
--- a/sysdeps/ieee754/flt-32/e_jnf.c
+++ b/sysdeps/ieee754/flt-32/e_jnf.c
@@ -165,7 +165,16 @@ static float zero  =  0.0000000000e+00;
 			}
 		    }
 		}
-		b = (t*__ieee754_j0f(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_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;
-- 
cgit v1.1