aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--ChangeLog10
-rw-r--r--math/math_private.h8
-rw-r--r--sysdeps/ieee754/dbl-64/s_fma.c22
-rw-r--r--sysdeps/ieee754/dbl-64/s_fmaf.c10
-rw-r--r--sysdeps/ieee754/flt-32/e_exp2f.c10
-rw-r--r--sysdeps/ieee754/flt-32/e_expf.c14
-rw-r--r--sysdeps/ieee754/flt-32/s_nearbyintf.c28
-rw-r--r--sysdeps/x86_64/fpu/math_private.h35
8 files changed, 81 insertions, 56 deletions
diff --git a/ChangeLog b/ChangeLog
index de8b2c3b98..852d45d800 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -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)