This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Re: Fix for bz#13658 - using fsincos on x86-64
On Wed, 7 Mar 2012, Joseph S. Myers wrote:
> On Wed, 7 Mar 2012, Andreas Jaeger wrote:
>
> > test-double fails with a lot of tests, starting with:
> > testing double (without inline functions)
> > Failure: Test: cos (M_PI_6l * 2.0) == 0.5
> > Result:
> > is: 4.99861241344128237607e-01 0x1.ffdba01071a860000000p-2
> > should be: 5.00000000000000000000e-01 0x1.00000000000000000000p-1
> > difference: 1.38758655871762393019e-04 0x1.22ff7c72bd0000000000p-13
> > ulp : 1249826861757.0000
> > max.ulp : 1.0000
> >
> > I need some more time to investigate this - and would appreciate any help,
>
> One possibility is that the C implementations of sin and cos only work
> with rounding to 53 bits, not with x87 excess precision. If so, then
> setting the FPU to 53-bit mode should help (and the libc_feupdateenv calls
> should handle restoring 64-bit mode), though some new macro would be
> needed to handle the environment setup for the affected functions (or
> you could have an i386 version of s_sin.c that includes headers, redefines
> libc_feholdexcept_setround then includes the main s_sin.c).
Something like this. Together with removing the files you listed in
<http://sourceware.org/ml/libc-alpha/2012-03/msg00164.html>, this appears
to work on x86 (some tests of j1/jn/yn increase their ulps from 1ulp to
2ulp to me, but I think that's fine as a side-effect of slightly more
accurate sincos; there are no other math/ failures) - though I haven't
tested with any new cases. Defining new macros for setting 53-bit
precision and restoring extended precision is a bit x86-specific, but then
I think it's OK for the functions to declare their requirements regarding
excess precision like that, just as they declare that they need
round-to-nearest mode. Changing the control word twice on both entry and
exit from the functions isn't optimal, but that's probably best optimized
once Richard's changes to how this rounding mode save / restore works have
gone in.
I hope we can combine this with your changes (plus ulps updates) to get
something working for x86 and x86_64, sin, cos and sincos, that can go in
fairly soon as a basis for further incremental work.
diff --git a/sysdeps/generic/math_private.h b/sysdeps/generic/math_private.h
index 777762d..be1e4d2 100644
--- a/sysdeps/generic/math_private.h
+++ b/sysdeps/generic/math_private.h
@@ -370,6 +370,9 @@ extern void __docos (double __x, double __dx, double __v[]);
#define libc_feholdexcept_setroundl(e, r) \
do { feholdexcept (e); fesetround (r); } while (0)
+#define libc_feholdexcept_setround_53bit(e, r) \
+ libc_feholdexcept_setround (e, r)
+
#define libc_fetestexcept(e) fetestexcept (e)
#define libc_fetestexceptf(e) fetestexcept (e)
#define libc_fetestexceptl(e) fetestexcept (e)
@@ -382,6 +385,8 @@ extern void __docos (double __x, double __dx, double __v[]);
#define libc_feupdateenvf(e) (void) feupdateenv (e)
#define libc_feupdateenvl(e) (void) feupdateenv (e)
+#define libc_feupdateenv_53bit(e) libc_feupdateenv (e)
+
#define __nan(str) \
(__builtin_constant_p (str) && str[0] == '\0' ? NAN : __nan (str))
#define __nanf(str) \
diff --git a/sysdeps/i386/fpu/math_private.h b/sysdeps/i386/fpu/math_private.h
index 5253998..a0aaab8 100644
--- a/sysdeps/i386/fpu/math_private.h
+++ b/sysdeps/i386/fpu/math_private.h
@@ -16,4 +16,37 @@ do \
while (0)
#include_next <math_private.h>
+
+#ifdef NEED_53BIT_PRECISION
+
+# include <fpu_control.h>
+
+# undef libc_feholdexcept_setround_53bit
+# define libc_feholdexcept_setround_53bit(e, r) \
+ do \
+ { \
+ fpu_control_t cw; \
+ libc_feholdexcept_setround (e, r); \
+ _FPU_GETCW (cw); \
+ cw &= ~(fpu_control_t) _FPU_EXTENDED; \
+ cw |= _FPU_DOUBLE; \
+ _FPU_SETCW (cw); \
+ } \
+ while (0)
+
+# undef libc_feupdateenv_53bit
+# define libc_feupdateenv_53bit(e) \
+ do \
+ { \
+ fpu_control_t cw; \
+ libc_feupdateenv (e); \
+ _FPU_GETCW (cw); \
+ cw &= ~(fpu_control_t) _FPU_EXTENDED; \
+ cw |= _FPU_EXTENDED; \
+ _FPU_SETCW (cw); \
+ } \
+ while (0)
+
+#endif
+
#endif
diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c
index e3e3a2a..4b4b675 100644
--- a/sysdeps/ieee754/dbl-64/s_sin.c
+++ b/sysdeps/ieee754/dbl-64/s_sin.c
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001, 2009, 2011 Free Software Foundation
+ * Copyright (C) 2001-2012 Free Software Foundation
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
@@ -111,7 +111,7 @@ __sin(double x){
fenv_t env;
double retval = 0;
- libc_feholdexcept_setround (&env, FE_TONEAREST);
+ libc_feholdexcept_setround_53bit (&env, FE_TONEAREST);
u.x = x;
m = u.i[HIGH_HALF];
@@ -365,7 +365,7 @@ __sin(double x){
}
ret:
- libc_feupdateenv (&env);
+ libc_feupdateenv_53bit (&env);
return retval;
}
@@ -386,7 +386,7 @@ __cos(double x)
fenv_t env;
double retval = 0;
- libc_feholdexcept_setround (&env, FE_TONEAREST);
+ libc_feholdexcept_setround_53bit (&env, FE_TONEAREST);
u.x = x;
m = u.i[HIGH_HALF];
@@ -635,7 +635,7 @@ __cos(double x)
}
ret:
- libc_feupdateenv (&env);
+ libc_feupdateenv_53bit (&env);
return retval;
}
--
Joseph S. Myers
joseph@codesourcery.com