Use new internal libc_fe* interfaces in more functions

This commit is contained in:
Ulrich Drepper 2011-10-18 15:11:31 -04:00
parent 4855e3ddf5
commit d9a8d0abcc
8 changed files with 81 additions and 56 deletions

View File

@ -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.

View File

@ -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_ */

View File

@ -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));

View File

@ -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;
}

View File

@ -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

View File

@ -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;

View File

@ -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)

View File

@ -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)