Improve float range reduction accuracy near pi/2 (bug 21094) [committed]

Submitted by Joseph S. Myers on March 15, 2017, 10:01 p.m.

Details

Message ID alpine.DEB.2.20.1703152201150.25300@digraph.polyomino.org.uk
State New
Headers show

Commit Message

Joseph S. Myers March 15, 2017, 10:01 p.m.
Bug 21094 reports 3ulp errors of cosf and tanf for certain arguments
near pi/2 arising from the use of an insufficiently accurate range
reduction.  (To be clear, this is a quality-of-implementation issue
relating to the apparent intent of those particular cosf and tanf
implementations; 3ulp is within the general glibc accuracy goals, so
not inherently a bug.)

This patch fixes that error by making a wider range of cases use the
existing more accurate range reduction for arguments close to pi/2.
The wider range of values is still narrow enough for the "z -=
pio2_2;" in the more accurate case to be exact, as the code expects.

Tested for x86_64, x86 and mips64; no ulps updates needed (but at
least on mips64, the larger ulps were seen if the tests were added
without the substantive fix).  Committed.

(auto-libm-test-out-* diffs omitted below.)

2017-03-15  Joseph Myers  <joseph@codesourcery.com>

	[BZ #21094]
	* sysdeps/ieee754/flt-32/e_rem_pio2f.c (__ieee754_rem_pio2f): Use
	24+24+24-bit pi for wider range of values around pi/2.
	* math/auto-libm-test-in: Add more tests of cos and tan.
	* math/auto-libm-test-out-cos: Regenerated.
	* math/auto-libm-test-out-tan: Likewise.

Patch hide | download patch | download mbox

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index c8267aa..6a47044 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -4324,6 +4324,8 @@  cos -0x1.02e34cp+0
 cos 0xf.f0274p+4
 cos 0x3.042d88p+0
 cos 0x1.8475e5afd4481p+0
+cos 1.57079697
+cos -1.57079697
 
 cosh 0
 cosh -0
@@ -7236,6 +7238,8 @@  tan 0x2.091d68p+0
 tan -0x5.302ab9b18593264p+0
 tan 0x1.1ad374p+0
 tan -0x1.0d55b8p+0
+tan 1.57079697
+tan -1.57079697
 tan 0x1p-5
 tan 0x1p-10
 tan 0x1p-15
diff --git a/sysdeps/ieee754/flt-32/e_rem_pio2f.c b/sysdeps/ieee754/flt-32/e_rem_pio2f.c
index 0928373..c4d28c8 100644
--- a/sysdeps/ieee754/flt-32/e_rem_pio2f.c
+++ b/sysdeps/ieee754/flt-32/e_rem_pio2f.c
@@ -100,7 +100,7 @@  int32_t __ieee754_rem_pio2f(float x, float *y)
 	if(ix<0x4016cbe4) {  /* |x| < 3pi/4, special case with n=+-1 */
 	    if(hx>0) {
 		z = x - pio2_1;
-		if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
+		if((ix&0xffffffc0)!=0x3fc90fc0) { /* 24+24 bit pi OK */
 		    y[0] = z - pio2_1t;
 		    y[1] = (z-y[0])-pio2_1t;
 		} else {		/* near pi/2, use 24+24+24 bit pi */
@@ -111,7 +111,7 @@  int32_t __ieee754_rem_pio2f(float x, float *y)
 		return 1;
 	    } else {	/* negative x */
 		z = x + pio2_1;
-		if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
+		if((ix&0xffffffc0)!=0x3fc90fc0) { /* 24+24 bit pi OK */
 		    y[0] = z + pio2_1t;
 		    y[1] = (z-y[0])+pio2_1t;
 		} else {		/* near pi/2, use 24+24+24 bit pi */