diff options
| -rw-r--r-- | ChangeLog | 10 | ||||
| -rw-r--r-- | math/math_private.h | 8 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/s_fma.c | 22 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/s_fmaf.c | 10 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/e_exp2f.c | 10 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/e_expf.c | 14 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/s_nearbyintf.c | 28 | ||||
| -rw-r--r-- | sysdeps/x86_64/fpu/math_private.h | 35 |
8 files changed, 81 insertions, 56 deletions
@@ -1,5 +1,15 @@ 2011-10-18 Ulrich Drepper <drepper@gmail.com> + * math/math_private.h: Define defaults for libc_fetestexcept and + libc_feupdateenv. + * sysdeps/ieee754/dbl-64/s_fma.c: Use libc_fe* interfaces. + * sysdeps/ieee754/dbl-64/s_fmaf.c: Likewise. + * sysdeps/ieee754/flt-32/e_exp2f.c: Likewise. + * sysdeps/ieee754/flt-32/e_expf.c: Likewise. + * sysdeps/ieee754/flt-32/s_nearbyintf.c: Likewise. + * sysdeps/x86_64/fpu/math_private.h: Define special versions of + libc_fetestexcept and libc_feupdateenv. + * math/math_private.h: Define defaults for libc_feholdexcept_setround, libc_feholdexcept_setroundf, libc_feholdexcept_setroundl. * sysdeps/ieee754/dbl-64/e_exp2.c: Use libc_feholdexcept_setround. diff --git a/math/math_private.h b/math/math_private.h index 38ff09e7e0..db733c8ac7 100644 --- a/math/math_private.h +++ b/math/math_private.h @@ -383,8 +383,16 @@ extern void __docos (double __x, double __dx, double __v[]); #define libc_feholdexcept_setroundl(e, r) \ do { feholdexcept (e); fesetround (r); } while (0) +#define libc_fetestexcept(e) fetestexcept (e) +#define libc_fetestexceptf(e) fetestexcept (e) +#define libc_fetestexceptl(e) fetestexcept (e) + #define libc_fesetenv(e) (void) fesetenv (e) #define libc_fesetenvf(e) (void) fesetenv (e) #define libc_fesetenvl(e) (void) fesetenv (e) +#define libc_feupdateenv(e) (void) feupdateenv (e) +#define libc_feupdateenvf(e) (void) feupdateenv (e) +#define libc_feupdateenvl(e) (void) feupdateenv (e) + #endif /* _MATH_PRIVATE_H_ */ diff --git a/sysdeps/ieee754/dbl-64/s_fma.c b/sysdeps/ieee754/dbl-64/s_fma.c index 3b0bfd5ce6..14c6503de2 100644 --- a/sysdeps/ieee754/dbl-64/s_fma.c +++ b/sysdeps/ieee754/dbl-64/s_fma.c @@ -1,5 +1,5 @@ /* Compute x * y + z as ternary operation. - Copyright (C) 2010 Free Software Foundation, Inc. + Copyright (C) 2010, 2011 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jakub@redhat.com>, 2010. @@ -22,6 +22,7 @@ #include <math.h> #include <fenv.h> #include <ieee754.h> +#include <math_private.h> /* This implementation uses rounding to odd to avoid problems with double rounding. See a paper by Boldo and Melquiond: @@ -47,7 +48,7 @@ __fma (double x, double y, double z) z rather than NaN. */ if (w.ieee.exponent == 0x7ff && u.ieee.exponent != 0x7ff - && v.ieee.exponent != 0x7ff) + && v.ieee.exponent != 0x7ff) return (z + x) + y; /* If x or y or z is Inf/NaN, or if fma will certainly overflow, or if x * y is less than half of DBL_DENORM_MIN, @@ -148,34 +149,33 @@ __fma (double x, double y, double z) double a2 = t1 + t2; fenv_t env; - feholdexcept (&env); - fesetround (FE_TOWARDZERO); + libc_feholdexcept_setround (&env, FE_TOWARDZERO); /* Perform m2 + a2 addition with round to odd. */ u.d = a2 + m2; if (__builtin_expect (adjust == 0, 1)) { if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7ff) - u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0; - feupdateenv (&env); + u.ieee.mantissa1 |= libc_fetestexcept (FE_INEXACT) != 0; + libc_feupdateenv (&env); /* Result is a1 + u.d. */ return a1 + u.d; } else if (__builtin_expect (adjust > 0, 1)) { if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7ff) - u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0; - feupdateenv (&env); + u.ieee.mantissa1 |= libc_fetestexcept (FE_INEXACT) != 0; + libc_feupdateenv (&env); /* Result is a1 + u.d, scaled up. */ return (a1 + u.d) * 0x1p53; } else { if ((u.ieee.mantissa1 & 1) == 0) - u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0; + u.ieee.mantissa1 |= libc_fetestexcept (FE_INEXACT) != 0; v.d = a1 + u.d; - int j = fetestexcept (FE_INEXACT) != 0; - feupdateenv (&env); + int j = libc_fetestexcept (FE_INEXACT) != 0; + libc_feupdateenv (&env); /* Ensure the following computations are performed in default rounding mode instead of just reusing the round to zero computation. */ asm volatile ("" : "=m" (u) : "m" (u)); diff --git a/sysdeps/ieee754/dbl-64/s_fmaf.c b/sysdeps/ieee754/dbl-64/s_fmaf.c index cd16cd1dce..dc748e5548 100644 --- a/sysdeps/ieee754/dbl-64/s_fmaf.c +++ b/sysdeps/ieee754/dbl-64/s_fmaf.c @@ -1,5 +1,5 @@ /* Compute x * y + z as ternary operation. - Copyright (C) 2010 Free Software Foundation, Inc. + Copyright (C) 2010, 2011 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jakub@redhat.com>, 2010. @@ -21,6 +21,7 @@ #include <math.h> #include <fenv.h> #include <ieee754.h> +#include <math_private.h> /* This implementation relies on double being more than twice as precise as float and uses rounding to odd in order to avoid problems @@ -35,13 +36,12 @@ __fmaf (float x, float y, float z) /* Multiplication is always exact. */ double temp = (double) x * (double) y; union ieee754_double u; - feholdexcept (&env); - fesetround (FE_TOWARDZERO); + libc_feholdexcept_setroundf (&env, FE_TOWARDZERO); /* Perform addition with round to odd. */ u.d = temp + (double) z; if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7ff) - u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0; - feupdateenv (&env); + u.ieee.mantissa1 |= libc_fetestexcept (FE_INEXACT) != 0; + libc_feupdateenv (&env); /* And finally truncation with round to nearest. */ return (float) u.d; } diff --git a/sysdeps/ieee754/flt-32/e_exp2f.c b/sysdeps/ieee754/flt-32/e_exp2f.c index 0703cea403..f07500d20d 100644 --- a/sysdeps/ieee754/flt-32/e_exp2f.c +++ b/sysdeps/ieee754/flt-32/e_exp2f.c @@ -57,11 +57,7 @@ __ieee754_exp2f (float x) union ieee754_float ex2_u, scale_u; fenv_t oldenv; - feholdexcept (&oldenv); -#ifdef FE_TONEAREST - /* If we don't have this, it's too bad. */ - fesetround (FE_TONEAREST); -#endif + libc_feholdexcept_setroundf (&oldenv, FE_TONEAREST); /* 1. Argument reduction. Choose integers ex, -128 <= t < 128, and some real @@ -104,7 +100,7 @@ __ieee754_exp2f (float x) x22 = (.24022656679f * x + .69314736128f) * ex2_u.f; /* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex). */ - fesetenv (&oldenv); + libc_fesetenv (&oldenv); result = x22 * x + ex2_u.f; @@ -116,7 +112,7 @@ __ieee754_exp2f (float x) /* Exceptional cases: */ else if (isless (x, himark)) { - if (__isinff (x)) + if (__isinf_nsf (x)) /* e^-inf == 0, with no error. */ return 0; else diff --git a/sysdeps/ieee754/flt-32/e_expf.c b/sysdeps/ieee754/flt-32/e_expf.c index 872d34bc5b..02105c4385 100644 --- a/sysdeps/ieee754/flt-32/e_expf.c +++ b/sysdeps/ieee754/flt-32/e_expf.c @@ -47,9 +47,6 @@ to perform an 'accurate table method' expf, because of the range reduction overhead (compare exp2f). */ -#ifndef _GNU_SOURCE -#define _GNU_SOURCE -#endif #include <float.h> #include <ieee754.h> #include <math.h> @@ -60,8 +57,8 @@ extern const float __exp_deltatable[178]; extern const double __exp_atable[355] /* __attribute__((mode(DF))) */; -static const volatile float TWOM100 = 7.88860905e-31; -static const volatile float TWO127 = 1.7014118346e+38; +static const float TWOM100 = 7.88860905e-31; +static const float TWO127 = 1.7014118346e+38; float __ieee754_expf (float x) @@ -86,10 +83,7 @@ __ieee754_expf (float x) union ieee754_double ex2_u; fenv_t oldenv; - feholdexcept (&oldenv); -#ifdef FE_TONEAREST - fesetround (FE_TONEAREST); -#endif + libc_feholdexcept_setroundf (&oldenv, FE_TONEAREST); /* Calculate n. */ n = x * M_1_LN2 + THREEp22; @@ -119,7 +113,7 @@ __ieee754_expf (float x) x22 = (0.5000000496709180453 * dx + 1.0000001192102037084) * dx + delta; /* Return result. */ - fesetenv (&oldenv); + libc_fesetenvf (&oldenv); result = x22 * ex2_u.d + ex2_u.d; return (float) result; diff --git a/sysdeps/ieee754/flt-32/s_nearbyintf.c b/sysdeps/ieee754/flt-32/s_nearbyintf.c index 7d6f262f51..04ef9ab20c 100644 --- a/sysdeps/ieee754/flt-32/s_nearbyintf.c +++ b/sysdeps/ieee754/flt-32/s_nearbyintf.c @@ -19,22 +19,14 @@ #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const float -#else -static float -#endif TWO23[2]={ 8.3886080000e+06, /* 0x4b000000 */ -8.3886080000e+06, /* 0xcb000000 */ }; -#ifdef __STDC__ - float __nearbyintf(float x) -#else - float __nearbyintf(x) - float x; -#endif +float +__nearbyintf(float x) { fenv_t env; int32_t i0,j0,sx; @@ -50,13 +42,13 @@ TWO23[2]={ i0 &= 0xfff00000; i0 |= ((i1|-i1)>>9)&0x400000; SET_FLOAT_WORD(x,i0); - feholdexcept (&env); - w = TWO23[sx]+x; - t = w-TWO23[sx]; - fesetenv (&env); + libc_feholdexceptf (&env); + w = TWO23[sx]+x; + t = w-TWO23[sx]; + libc_fesetenvf (&env); GET_FLOAT_WORD(i0,t); SET_FLOAT_WORD(t,(i0&0x7fffffff)|(sx<<31)); - return t; + return t; } else { i = (0x007fffff)>>j0; if((i0&i)==0) return x; /* x is integral */ @@ -64,14 +56,14 @@ TWO23[2]={ if((i0&i)!=0) i0 = (i0&(~i))|((0x100000)>>j0); } } else { - if(j0==0x80) return x+x; /* inf or NaN */ + if(__builtin_expect(j0==0x80, 0)) return x+x; /* inf or NaN */ else return x; /* x is integral */ } SET_FLOAT_WORD(x,i0); - feholdexcept (&env); + libc_feholdexceptf (&env); w = TWO23[sx]+x; t = w-TWO23[sx]; - fesetenv (&env); + libc_fesetenvf (&env); return t; } weak_alias (__nearbyintf, nearbyintf) diff --git a/sysdeps/x86_64/fpu/math_private.h b/sysdeps/x86_64/fpu/math_private.h index 28bd9cee2d..de4f121075 100644 --- a/sysdeps/x86_64/fpu/math_private.h +++ b/sysdeps/x86_64/fpu/math_private.h @@ -129,7 +129,8 @@ do { \ asm volatile ("stmxcsr %0" : "=m" (*&mxcsr)); \ (mxcsr & 0x6000) >> 3; \ }) -// #define libc_fegetroundf() fegetround () +#undef libc_fegetroundf +#define libc_fegetroundf() libc_fegetround () // #define libc_fegetroundl() fegetround () #undef libc_fesetround @@ -140,7 +141,8 @@ do { \ mxcsr = (mxcsr & ~0x6000) | ((r) << 3); \ asm volatile ("ldmxcsr %0" : : "m" (*&mxcsr)); \ } while (0) -// #define libc_fesetroundf(r) (void) fesetround (r) +#undef libc_fesetroundf +#define libc_fesetroundf(r) libc_fesetround (r) // #define libc_fesetroundl(r) (void) fesetround (r) #undef libc_feholdexcept @@ -152,7 +154,8 @@ do { \ mxcsr = (mxcsr | 0x1f80) & ~0x3f; \ asm volatile ("ldmxcsr %0" : : "m" (*&mxcsr)); \ } while (0) -// #define libc_feholdexceptf(e) (void) feholdexcept (e) +#undef libc_feholdexceptf +#define libc_feholdexceptf(e) libc_feholdexcept (e) // #define libc_feholdexceptl(e) (void) feholdexcept (e) #undef libc_feholdexcept_setround @@ -164,11 +167,33 @@ do { \ mxcsr = ((mxcsr | 0x1f80) & ~0x603f) | ((r) << 3); \ asm volatile ("ldmxcsr %0" : : "m" (*&mxcsr)); \ } while (0) -// #define libc_feholdexcept_setroundf(e, r) ... +#undef libc_feholdexcept_setroundf +#define libc_feholdexcept_setroundf(e, r) libc_feholdexcept_setround (e, r) // #define libc_feholdexcept_setroundl(e, r) ... +#undef libc_fetestexcept +#define libc_fetestexcept(e) \ + ({ unsigned int mxcsr; asm volatile ("stmxcsr %0" : "=m" (*&mxcsr)); \ + mxcsr & (e) & FE_ALL_EXCEPT; }) +#undef libc_fetestexceptf +#define libc_fetestexceptf(e) libc_fetestexcept (e) +// #define libc_fetestexceptl(e) fetestexcept (e) + #undef libc_fesetenv #define libc_fesetenv(e) \ asm volatile ("ldmxcsr %0" : : "m" ((e)->__mxcsr)) -// #define libc_fesetenvf(e) (void) fesetenv (e) +#undef libc_fesetenvf +#define libc_fesetenvf(e) libc_fesetenv (e) // #define libc_fesetenvl(e) (void) fesetenv (e) + +#undef libc_feupdateenv +#define libc_feupdateenv(e) \ + do { \ + unsigned int mxcsr; \ + asm volatile ("stmxcsr %0" : "=m" (*&mxcsr)); \ + asm volatile ("ldmxcsr %0" : : "m" ((e)->__mxcsr)); \ + feraiseexcept (mxcsr & FE_ALL_EXCEPT); \ + } while (0) +#undef libc_feupdateenvf +#define libc_feupdateenvf(e) libc_feupdateenv (e) +// #define libc_feupdateenvl(e) (void) feupdateenv (e) |