This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix float range reduction problems (bug 14283)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: <libc-alpha at sourceware dot org>
- Date: Fri, 22 Jun 2012 20:53:23 +0000
- Subject: Fix float range reduction problems (bug 14283)
This patch (for 2.17) fixes bug 14283 (inaccurate range reduction for
sinf and cosf). There were two problems I found from the sample tests
I extracted from the provided log: a mistake in adapting the range
reduction code from double to float meant a shift was off-by-one, so
causing loss of precision for values close to a multiple of pi/2, and
the logic producing the final array of float failed to allow for
FLT_EVAL_METHOD the way the double implementation of range reduction
does, causing bad results for 32-bit x86 only.
Tested x86_64 and x86; ulps updates were only needed for x86 in that
testing (the more accurate range reduction ends up adding 1ulp to the
errors for some j1 and y0 tests with large inputs).
2012-06-22 Joseph Myers <joseph@codesourcery.com>
[BZ #14283]
* sysdeps/ieee754/flt-32/k_rem_pio2f.c (__kernel_rem_pio2f): Shift
by 7 not 8 to examine high bit of fractional part. Use volatile
variables when splitting into final array of floats if
__FLT_EVAL_METHOD__ != 0.
* math/libm-test.inc (cos_test): Add another test.
(sin_test): Likewise.
* sysdeps/i386/fpu/libm-test-ulps: Update.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 8e4d02e..e5afa06 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -2589,6 +2589,8 @@ cos_test (void)
TEST_f_f (cos, 0.80190127184058835, 0.69534156199418473);
#endif
+ TEST_f_f (cos, 0x1.442f74p+15, 2.4407839902314016628485779006274989801517e-06L);
+
#ifndef TEST_FLOAT
TEST_f_f (cos, 1e22, 0.5232147853951389454975944733847094921409L);
TEST_f_f (cos, 0x1p1023, -0.826369834614147994500785680811743734805L);
@@ -7687,6 +7689,8 @@ sin_test (void)
TEST_f_f (sin, 0x1p65, -0.047183876212354673805106149805700013943218L);
TEST_f_f (sin, -0x1p65, 0.047183876212354673805106149805700013943218L);
+ TEST_f_f (sin, 0x1.7f4134p+103, -6.6703229329788657073304190650534846045235e-08L);
+
#ifdef TEST_DOUBLE
TEST_f_f (sin, 0.80190127184058835, 0.71867942238767868);
TEST_f_f (sin, 2.522464e-1, 2.4957989804940911e-1);
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 38a69e6..4661aea 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -1656,8 +1656,8 @@ Test "j1 (0.75) == 0.349243602174862192523281016426251335":
double: 1
idouble: 1
Test "j1 (0x1.3ffp+74) == 1.818984347516051243459364437186082741567e-12":
-float: 1
-ifloat: 1
+float: 2
+ifloat: 2
ildouble: 1
ldouble: 1
Test "j1 (0x1.ff00000000002p+840) == 1.846591691699331493194965158699937660696e-127":
@@ -2359,8 +2359,8 @@ float: 1
idouble: 1
ifloat: 1
Test "y0 (0x1.3ffp+74) == 1.818984347516051243459467456433028748678e-12":
-float: 1
-ifloat: 1
+float: 2
+ifloat: 2
ildouble: 1
ldouble: 1
Test "y0 (0x1.ff00000000002p+840) == 1.846591691699331493194965158699937660696e-127":
@@ -2956,9 +2956,9 @@ ldouble: 2
Function: "j1":
double: 2
-float: 1
+float: 2
idouble: 2
-ifloat: 1
+ifloat: 2
ildouble: 1
ldouble: 1
@@ -3128,9 +3128,9 @@ ldouble: 1
Function: "y0":
double: 2
-float: 1
+float: 2
idouble: 2
-ifloat: 1
+ifloat: 2
ildouble: 1
ldouble: 1
diff --git a/sysdeps/ieee754/flt-32/k_rem_pio2f.c b/sysdeps/ieee754/flt-32/k_rem_pio2f.c
index 06c2f37..7989c72 100644
--- a/sysdeps/ieee754/flt-32/k_rem_pio2f.c
+++ b/sysdeps/ieee754/flt-32/k_rem_pio2f.c
@@ -88,7 +88,7 @@ recompute:
iq[jz-1] -= i<<(8-q0);
ih = iq[jz-1]>>(7-q0);
}
- else if(q0==0) ih = iq[jz-1]>>8;
+ else if(q0==0) ih = iq[jz-1]>>7;
else if(z>=(float)0.5) ih=2;
if(ih>0) { /* q > 0.5 */
@@ -166,24 +166,33 @@ recompute:
y[0] = (ih==0)? fw: -fw;
break;
case 1:
- case 2:
- fw = 0.0;
- for (i=jz;i>=0;i--) fw += fq[i];
- y[0] = (ih==0)? fw: -fw;
- fw = fq[0]-fw;
- for (i=1;i<=jz;i++) fw += fq[i];
- y[1] = (ih==0)? fw: -fw;
+ case 2:;
+#if __FLT_EVAL_METHOD__ != 0
+ volatile
+#endif
+ float fv = 0.0;
+ for (i=jz;i>=0;i--) fv += fq[i];
+ y[0] = (ih==0)? fv: -fv;
+ fv = fq[0]-fv;
+ for (i=1;i<=jz;i++) fv += fq[i];
+ y[1] = (ih==0)? fv: -fv;
break;
case 3: /* painful */
for (i=jz;i>0;i--) {
- fw = fq[i-1]+fq[i];
- fq[i] += fq[i-1]-fw;
- fq[i-1] = fw;
+#if __FLT_EVAL_METHOD__ != 0
+ volatile
+#endif
+ float fv = fq[i-1]+fq[i];
+ fq[i] += fq[i-1]-fv;
+ fq[i-1] = fv;
}
for (i=jz;i>1;i--) {
- fw = fq[i-1]+fq[i];
- fq[i] += fq[i-1]-fw;
- fq[i-1] = fw;
+#if __FLT_EVAL_METHOD__ != 0
+ volatile
+#endif
+ float fv = fq[i-1]+fq[i];
+ fq[i] += fq[i-1]-fv;
+ fq[i-1] = fv;
}
for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
if(ih==0) {
--
Joseph S. Myers
joseph@codesourcery.com