Fix spurious underflow exceptions for Bessel functions for double (bug 14155).

4 files changed, 23 insertions, 8 deletions

diff --git a/ChangeLog b/ChangeLog
index 8fb2ba76a4..2a73b8b61b 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,16 @@
+2013-03-14  Joseph Myers  <joseph@codesourcery.com>
+
+	[BZ #14155]
+	* sysdeps/ieee754/dbl-64/e_j0.c (pzero): Return 1.0 for arguments
+	0x1p28 and above.
+	(qzero): Return -0.125 / x for arguments 0x1p28 and above.
+	* sysdeps/ieee754/dbl-64/e_j1.c (pzero): Return 1.0 for arguments
+	0x1p28 and above.
+	(qzero): Return 0.375 / x for arguments 0x1p28 and above.
+	* math/libm-test.inc (j0_test): Do not allow one spurious
+	underflow exception.
+	(y1_test): Likewise.
+
 2013-03-14  Siddhesh Poyarekar  <siddhesh@redhat.com>
 
 	* manual/Makefile (chapters): Add nptl.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 6ac3cd2d0b..132a540e28 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -6239,8 +6239,7 @@ j0_test (void)
   TEST_f_f (j0, 4.0, -3.9714980986384737228659076845169804197562E-1L);
   TEST_f_f (j0, -4.0, -3.9714980986384737228659076845169804197562E-1L);
 
-  /* Bug 14155: spurious exception may occur.  */
-  TEST_f_f (j0, 0x1.d7ce3ap+107L, 2.775523647291230802651040996274861694514e-17L, UNDERFLOW_EXCEPTION_OK);
+  TEST_f_f (j0, 0x1.d7ce3ap+107L, 2.775523647291230802651040996274861694514e-17L);
 
 #ifndef TEST_FLOAT
   /* Bug 14155: spurious exception may occur.  */
@@ -10494,8 +10493,7 @@ y1_test (void)
   TEST_f_f (y1, 8.0, -0.158060461731247494255555266187483550L);
   TEST_f_f (y1, 10.0, 0.249015424206953883923283474663222803L);
 
-  /* Bug 14155: spurious exception may occur.  */
-  TEST_f_f (y1, 0x1.27e204p+99L, -8.881610148467797208469612080785210013461e-16L, UNDERFLOW_EXCEPTION_OK);
+  TEST_f_f (y1, 0x1.27e204p+99L, -8.881610148467797208469612080785210013461e-16L);
 
 #ifndef TEST_FLOAT
   /* Bug 14155: spurious exception may occur.  */
diff --git a/sysdeps/ieee754/dbl-64/e_j0.c b/sysdeps/ieee754/dbl-64/e_j0.c
index f393a762b2..d641a09149 100644
--- a/sysdeps/ieee754/dbl-64/e_j0.c
+++ b/sysdeps/ieee754/dbl-64/e_j0.c
@@ -293,7 +293,8 @@ pzero(double x)
 	int32_t ix;
 	GET_HIGH_WORD(ix,x);
 	ix &= 0x7fffffff;
-	if(ix>=0x40200000)     {p = pR8; q= pS8;}
+	if (ix>=0x41b00000)    {return one;}
+	else if(ix>=0x40200000){p = pR8; q= pS8;}
 	else if(ix>=0x40122E8B){p = pR5; q= pS5;}
 	else if(ix>=0x4006DB6D){p = pR3; q= pS3;}
 	else if(ix>=0x40000000){p = pR2; q= pS2;}
@@ -400,7 +401,8 @@ qzero(double x)
 	int32_t ix;
 	GET_HIGH_WORD(ix,x);
 	ix &= 0x7fffffff;
-	if(ix>=0x40200000)     {p = qR8; q= qS8;}
+	if (ix>=0x41b00000)    {return -.125/x;}
+	else if(ix>=0x40200000){p = qR8; q= qS8;}
 	else if(ix>=0x40122E8B){p = qR5; q= qS5;}
 	else if(ix>=0x4006DB6D){p = qR3; q= qS3;}
 	else if(ix>=0x40000000){p = qR2; q= qS2;}
diff --git a/sysdeps/ieee754/dbl-64/e_j1.c b/sysdeps/ieee754/dbl-64/e_j1.c
index cba4d46b18..cca5f20b4f 100644
--- a/sysdeps/ieee754/dbl-64/e_j1.c
+++ b/sysdeps/ieee754/dbl-64/e_j1.c
@@ -291,7 +291,8 @@ pone(double x)
 	int32_t ix;
 	GET_HIGH_WORD(ix,x);
 	ix &= 0x7fffffff;
-	if(ix>=0x40200000)     {p = pr8; q= ps8;}
+	if (ix>=0x41b00000)    {return one;}
+	else if(ix>=0x40200000){p = pr8; q= ps8;}
 	else if(ix>=0x40122E8B){p = pr5; q= ps5;}
 	else if(ix>=0x4006DB6D){p = pr3; q= ps3;}
 	else if(ix>=0x40000000){p = pr2; q= ps2;}
@@ -399,7 +400,8 @@ qone(double x)
 	int32_t ix;
 	GET_HIGH_WORD(ix,x);
 	ix &= 0x7fffffff;
-	if(ix>=0x40200000)     {p = qr8; q= qs8;}
+	if (ix>=0x41b00000)    {return .375/x;}
+	else if(ix>=0x40200000){p = qr8; q= qs8;}
 	else if(ix>=0x40122E8B){p = qr5; q= qs5;}
 	else if(ix>=0x4006DB6D){p = qr3; q= qs3;}
 	else if(ix>=0x40000000){p = qr2; q= qs2;}