This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

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


Bug 14155 is spurious underflow exceptions from Bessel functions for
large arguments.  (The correct results for large x are roughly
constant * sin or cos (x + constant) / sqrt (x), so no underflow
exceptions should occur based on the final result.)

There are various places underflows may occur in the intermediate
calculations that cause the failures listed in that bug.  This patch
fixes problems for the double version where underflows occur in
calculating the intermediate functions P and Q (in particular, x**-12
gets computed while calculating Q).  Appropriate approximations are
used for P and Q for arguments at least 0x1p28 and above to avoid the
underflows.

For sufficiently large x - 0x1p129 and above - the code already has a
cut-off to avoid calculating P and Q at all, which means the
approximations -0.125 / x and 0.375 / x can't themselves cause
underflows calculating Q.  This cut-off is heuristically reasonable
for the point beyond which Q can be neglected (based on expecting
around 0x1p-64 to be the least absolute value of sin or cos for large
arguments representable in double).

The float versions use a cut-off 0x1p17, which is less heuristically
justifiable but should still only affect values near zeroes of the
Bessel functions where these implementations are intrinsically
inaccurate anyway (bugs 14469-14472), and should serve to avoid
underflows (the float underflow for jn in bug 14155 probably comes
from the recurrence to compute jn).  ldbl-96 uses 0x1p129, which may
not really be enough heuristically (0x1p143 or so might be safer - 143
= 64 + 79, number of mantissa bits plus total number of significant
bits in representation) but again should avoid underflows and only
affect values where the code is substantially inaccurate anyway.
ldbl-128 and ldbl-128ibm share a completely different implementation
with no such cut-off, which I propose to fix separately.

Tested x86_64 and x86.

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.

diff --git a/math/libm-test.inc b/math/libm-test.inc
index 6ac3cd2..132a540 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 f393a76..d641a09 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 cba4d46..cca5f20 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;}

-- 
Joseph S. Myers
joseph@codesourcery.com


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]