From 05b227bdaea9a5f1faf08dad31221d8736f3659d Mon Sep 17 00:00:00 2001
From: "David S. Miller" <davem@davemloft.net>
Date: Sun, 18 Nov 2012 12:33:53 -0800
Subject: Correct tinyness handling in long-double and float y0/y1.

	With help from Joseph Myers.
	* sysdeps/ieee754/flt-32/e_j0f.c (__ieee754_y0f): Adjust tinyness
	cutoff to 2**-13.
	* sysdeps/ieee754/flt-32/e_j1f.c (__ieee754_y1f): Adjust tinyness
	cutoff to 2**-25.
	* sysdeps/ieee754/ldbl-128/e_j0l.c (U0): New constant.
	( __ieee754_y0l): Avoid arithmetic underflow when 'x' is very
	small.
	* sysdeps/ieee754/ldbl-128/e_j1l.c (__ieee754_y1l): Likewise.
	* math/libm-test.inc (y0_test): New tests.
	(y1_test): New tests.
	* sysdeps/i386/fpu/libm-test-ulps: Update.
	* sysdeps/x86_64/fpu/libm-test-ulps: Update.
	* sysdeps/sparc/fpu/libm-test-ulps: Update.
---
 sysdeps/ieee754/flt-32/e_j0f.c   | 2 +-
 sysdeps/ieee754/flt-32/e_j1f.c   | 2 +-
 sysdeps/ieee754/ldbl-128/e_j0l.c | 3 +++
 sysdeps/ieee754/ldbl-128/e_j1l.c | 2 ++
 4 files changed, 7 insertions(+), 2 deletions(-)

(limited to 'sysdeps/ieee754')

diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c
index 0729cd04e0..c4cabd584a 100644
--- a/sysdeps/ieee754/flt-32/e_j0f.c
+++ b/sysdeps/ieee754/flt-32/e_j0f.c
@@ -138,7 +138,7 @@ __ieee754_y0f(float x)
 		}
 		return z;
 	}
-	if(ix<=0x32000000) {	/* x < 2**-27 */
+	if(ix<=0x39800000) {	/* x < 2**-13 */
 	    return(u00 + tpi*__ieee754_logf(x));
 	}
 	z = x*x;
diff --git a/sysdeps/ieee754/flt-32/e_j1f.c b/sysdeps/ieee754/flt-32/e_j1f.c
index 30b7d8e250..cb9f97fa28 100644
--- a/sysdeps/ieee754/flt-32/e_j1f.c
+++ b/sysdeps/ieee754/flt-32/e_j1f.c
@@ -133,7 +133,7 @@ __ieee754_y1f(float x)
 		}
 		return z;
 	}
-	if(__builtin_expect(ix<=0x24800000, 0)) {    /* x < 2**-54 */
+	if(__builtin_expect(ix<=0x33000000, 0)) {    /* x < 2**-25 */
 	    return(-tpi/x);
 	}
 	z = x*x;
diff --git a/sysdeps/ieee754/ldbl-128/e_j0l.c b/sysdeps/ieee754/ldbl-128/e_j0l.c
index 112a8f3f9c..1b18289588 100644
--- a/sysdeps/ieee754/ldbl-128/e_j0l.c
+++ b/sysdeps/ieee754/ldbl-128/e_j0l.c
@@ -809,6 +809,7 @@ static long double Y0_2D[NY0_2D + 1] = {
  /* 1.000000000000000000000000000000000000000E0 */
 };
 
+static const long double U0 = -7.3804295108687225274343927948483016310862e-02L;
 
 /* Bessel function of the second kind, order zero.  */
 
@@ -831,6 +832,8 @@ long double
       return -HUGE_VALL + x;
     }
   xx = fabsl (x);
+  if (xx <= 0x1p-57)
+    return U0 + TWOOPI * __ieee754_logl (x);
   if (xx <= 2.0L)
     {
       /* 0 <= x <= 2 */
diff --git a/sysdeps/ieee754/ldbl-128/e_j1l.c b/sysdeps/ieee754/ldbl-128/e_j1l.c
index 1f62bd0920..f16343b26b 100644
--- a/sysdeps/ieee754/ldbl-128/e_j1l.c
+++ b/sysdeps/ieee754/ldbl-128/e_j1l.c
@@ -838,6 +838,8 @@ __ieee754_y1l (long double x)
       return -HUGE_VALL + x;
     }
   xx = fabsl (x);
+  if (xx <= 0x1p-114)
+    return -TWOOPI / x;
   if (xx <= 2.0L)
     {
       /* 0 <= x <= 2 */
-- 
cgit 1.4.1