Fix spurious jnf underflows (bug 14155).

6 files changed, 68 insertions, 10 deletions

diff --git a/ChangeLog b/ChangeLog
index 608e133..34f7c4a 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,15 @@
+2013-09-02  Joseph Myers  <joseph@codesourcery.com>
+
+	[BZ #14155]
+	* sysdeps/ieee754/flt-32/e_jnf.c (__ieee754_jnf): Use double for
+	intermediate calculations in recurrence.
+	(__ieee754_ynf): Likewise.
+	* math/libm-test.inc (jn_test_data): Do not allow spurious
+	underflow exception.  Add more tests.
+	(yn_test_data): Add more tests.
+	* sysdeps/i386/fpu/libm-test-ulps: Update.
+	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
+
 2013-09-02   Ondřej Bílka  <neleai@seznam.cz>
 
 	* sysdeps/powerpc/powerpc32/power6/memcpy.S: Fix typo.
diff --git a/NEWS b/NEWS
index 4fe4b63..17e8eb6 100644
--- a/NEWS
+++ b/NEWS
@@ -9,8 +9,8 @@ Version 2.19
 
 * The following bugs are resolved with this release:
 
-  14699, 15522, 15531, 15532, 15736, 15749, 15797, 15867, 15890, 15897,
-  15905, 15909.
+  14155, 14699, 15522, 15531, 15532, 15736, 15749, 15797, 15867, 15890,
+  15897, 15905, 15909.
 
 * CVE-2013-4237 The readdir_r function could write more than NAME_MAX bytes
   to the d_name member of struct dirent, or omit the terminating NUL
diff --git a/math/libm-test.inc b/math/libm-test.inc
index e534fc0..fb5e977 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -9913,8 +9913,14 @@ static const struct test_if_f_data jn_test_data[] =
     TEST_if_f (jn, 8, 2.4048255576957729L, 0.92165786705344923232879022467054148E-4L),
     TEST_if_f (jn, 9, 2.4048255576957729L, 0.12517270977961513005428966643852564E-4L),
 
-    /* Bug 14155: spurious exception may occur.  */
-    TEST_if_f (jn, 2, 0x1.ffff62p+99L, -4.43860668048170034334926693188979974489e-16L, UNDERFLOW_EXCEPTION_OK),
+    TEST_if_f (jn, 2, 0x1.ffff62p+99L, -4.43860668048170034334926693188979974489e-16L),
+    TEST_if_f (jn, 2, 0x1p127L, -6.0784021821505059176832624052765568656702e-20L),
+#ifndef TEST_FLOAT
+    TEST_if_f (jn, 2, 0x1p1023L, 1.5665258060609012834424478437196679802783e-155L),
+#endif
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+    TEST_if_f (jn, 2, 0x1p16383L, -9.5859502826270374691362975419147645151233e-2467L),
+#endif
   };
 
 static void
@@ -14531,6 +14537,15 @@ static const struct test_if_f_data yn_test_data[] =
     /* Check whether yn returns correct value for LDBL_MIN, DBL_MIN,
        and FLT_MIN.  See Bug 14173.  */
     TEST_if_f (yn, 10, min_value, minus_infty, OVERFLOW_EXCEPTION|ERRNO_ERANGE),
+
+    TEST_if_f (yn, 2, 0x1.ffff62p+99L, -5.5244413477397111790415387179517953221757e-16L),
+    TEST_if_f (yn, 2, 0x1p127L, 6.8569250690166637098111268958532649249771e-21L),
+#ifndef TEST_FLOAT
+    TEST_if_f (yn, 2, 0x1p1023L, -8.2687542933709649327986678723012001545638e-155L),
+#endif
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+    TEST_if_f (yn, 2, 0x1p16383L, 3.8895531955766020648617743624167352352217e-2467L),
+#endif
   };
 
 static void
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 530dbd7..8244863 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -5635,9 +5635,9 @@ ildouble: 1
 ldouble: 1
 Test "jn (10, 10.0)":
 double: 1
-float: 1
+float: 2
 idouble: 1
-ifloat: 1
+ifloat: 2
 ildouble: 2
 ldouble: 2
 Test "jn (10, 2.0)":
@@ -5648,6 +5648,14 @@ float: 1
 ifloat: 1
 ildouble: 1
 ldouble: 1
+Test "jn (2, 0x1p1023)":
+double: 1
+idouble: 1
+Test "jn (2, 0x1p127)":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
 Test "jn (2, 2.4048255576957729)":
 double: 1
 float: 1
@@ -6844,6 +6852,14 @@ ifloat: 1
 Test "yn (10, 2.0)":
 float: 3
 ifloat: 3
+Test "yn (2, 0x1.ffff62p+99)":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "yn (2, 0x1p127)":
+float: 2
+ifloat: 2
 Test "yn (3, 0.125)":
 ildouble: 1
 ldouble: 1
diff --git a/sysdeps/ieee754/flt-32/e_jnf.c b/sysdeps/ieee754/flt-32/e_jnf.c
index ad26d7e..5984d94 100644
--- a/sysdeps/ieee754/flt-32/e_jnf.c
+++ b/sysdeps/ieee754/flt-32/e_jnf.c
@@ -54,7 +54,7 @@ __ieee754_jnf(int n, float x)
 	    b = __ieee754_j1f(x);
 	    for(i=1;i<n;i++){
 		temp = b;
-		b = b*((float)(i+i)/x) - a; /* avoid underflow */
+		b = b*((double)(i+i)/x) - a; /* avoid underflow */
 		a = temp;
 	    }
 	} else {
@@ -196,7 +196,7 @@ __ieee754_ynf(int n, float x)
 	GET_FLOAT_WORD(ib,b);
 	for(i=1;i<n&&ib!=0xff800000;i++){
 	    temp = b;
-	    b = ((float)(i+i)/x)*b - a;
+	    b = ((double)(i+i)/x)*b - a;
 	    GET_FLOAT_WORD(ib,b);
 	    a = temp;
 	}
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index d02618a..477eedc 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -6403,6 +6403,11 @@ idouble: 2
 ifloat: 2
 ildouble: 1
 ldouble: 1
+Test "jn (2, 0x1p127)":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
 Test "jn (2, 2.4048255576957729)":
 double: 2
 float: 1
@@ -7728,6 +7733,16 @@ double: 3
 float: 1
 idouble: 3
 ifloat: 1
+Test "yn (2, 0x1.ffff62p+99)":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "yn (2, 0x1p127)":
+double: 1
+float: 3
+idouble: 1
+ifloat: 3
 Test "yn (3, 0.125)":
 double: 1
 idouble: 1
@@ -8428,9 +8443,9 @@ ldouble: 2
 
 Function: "yn":
 double: 3
-float: 2
+float: 3
 idouble: 3
-ifloat: 2
+ifloat: 3
 ildouble: 4
 ldouble: 4