Create and use SET_RESTORE_ROUND{,_NOEX,_53BIT}{,F,L}.

This commit is contained in:
Richard Henderson 2012-03-10 08:55:53 -08:00
parent 7d2e8012cf
commit eb92c487b3
10 changed files with 221 additions and 118 deletions

View File

@ -1,5 +1,29 @@
2012-03-19 Richard Henderson <rth@twiddle.net>
* sysdeps/generic/math_private.h (libc_feholdsetround): New.
(libc_feholdsetroundf, libc_feholdsetroundl): New.
(libc_feresetround, libc_feresetroundf, libc_feresetroundl): New.
(libc_feresetround_noex): New.
(libc_feresetround_noexf): New.
(libc_feresetround_noexl): New.
(SET_RESTORE_ROUND, SET_RESTORE_ROUNDF, SET_RESTORE_ROUNDL): New.
(SET_RESTORE_ROUND_NOEX, SET_RESTORE_ROUND_NOEXF): New.
(SET_RESTORE_ROUND_NOEXL, SET_RESTORE_ROUND_53BIT): New.
* sysdeps/ieee754/dbl-64/e_exp.c (__ieee754_exp): Use
SET_RESTORE_ROUND.
* sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Likewise.
* sysdeps/ieee754/dbl-64/s_sin.c (__sin): Use SET_RESTORE_ROUND_53BIT.
(__cos): Likewise.
* sysdeps/ieee754/dbl-64/s_tan.c (__tan): Likewise.
* sysdeps/ieee754/dbl-64/e_exp2.c (__ieee754_exp2): Use
SET_RESTORE_ROUND_NOEX.
* sysdeps/ieee754/dbl-64/e_exp2f.c (__ieee754_exp2f): Use
SET_RESTORE_ROUND_NOEXF.
* sysdeps/ieee754/flt-32/e_expf.c (__ieee754_expf): Likewise.
* sysdeps/x86_64/fpu/math_private.h (libc_feholdsetround): New.
(libc_feholdsetroundf): New.
(libc_feresetround, libc_feresetroundf): New.
* sysdeps/i386/fpu/math_private.h: Include <fenv.h>, <fpu_control.h>.
(libc_feholdexcept_setround_53bit): Convert from macro to function.
(libc_feupdateenv_53bit): Likewise. Don't force _FPU_EXTENDED.

View File

@ -457,6 +457,75 @@ default_libc_feupdateenv (fenv_t *e)
# define libc_feupdateenv_53bit libc_feupdateenv
#endif
/* Save and set the rounding mode. The use of fenv_t to store the old mode
allows a target-specific version of this function to avoid converting the
rounding mode from the fpu format. By default we have no choice but to
manipulate the entire env. */
#ifndef libc_feholdsetround
# define libc_feholdsetround libc_feholdexcept_setround
#endif
#ifndef libc_feholdsetroundf
# define libc_feholdsetroundf libc_feholdexcept_setroundf
#endif
#ifndef libc_feholdsetroundl
# define libc_feholdsetroundl libc_feholdexcept_setroundl
#endif
/* ... and the reverse. */
#ifndef libc_feresetround
# define libc_feresetround libc_feupdateenv
#endif
#ifndef libc_feresetroundf
# define libc_feresetroundf libc_feupdateenvf
#endif
#ifndef libc_feresetroundl
# define libc_feresetroundl libc_feupdateenvl
#endif
/* ... and a version that may also discard exceptions. */
#ifndef libc_feresetround_noex
# define libc_feresetround_noex libc_fesetenv
#endif
#ifndef libc_feresetround_noexf
# define libc_feresetround_noexf libc_fesetenvf
#endif
#ifndef libc_feresetround_noexl
# define libc_feresetround_noexl libc_fesetenvl
#endif
/* Save and restore the rounding mode within a lexical block. */
#define SET_RESTORE_ROUND(RM) \
fenv_t __libc_save_rm __attribute__((cleanup(libc_feresetround))); \
libc_feholdsetround (&__libc_save_rm, (RM))
#define SET_RESTORE_ROUNDF(RM) \
fenv_t __libc_save_rm __attribute__((cleanup(libc_feresetroundf))); \
libc_feholdsetroundf (&__libc_save_rm, (RM))
#define SET_RESTORE_ROUNDL(RM) \
fenv_t __libc_save_rm __attribute__((cleanup(libc_feresetroundl))); \
libc_feholdsetroundl (&__libc_save_rm, (RM))
/* Save and restore the rounding mode within a lexical block, and also
the set of exceptions raised within the block may be discarded. */
#define SET_RESTORE_ROUND_NOEX(RM) \
fenv_t __libc_save_rm __attribute__((cleanup(libc_feresetround_noex))); \
libc_feholdsetround (&__libc_save_rm, (RM))
#define SET_RESTORE_ROUND_NOEXF(RM) \
fenv_t __libc_save_rm __attribute__((cleanup(libc_feresetround_noexf))); \
libc_feholdsetroundf (&__libc_save_rm, (RM))
#define SET_RESTORE_ROUND_NOEXL(RM) \
fenv_t __libc_save_rm __attribute__((cleanup(libc_feresetround_noexl))); \
libc_feholdsetroundl (&__libc_save_rm, (RM))
/* Like SET_RESTORE_ROUND, but also set rounding precision to 53 bits. */
#define SET_RESTORE_ROUND_53BIT(RM) \
fenv_t __libc_save_rm __attribute__((cleanup(libc_feupdateenv_53bit))); \
libc_feholdexcept_setround_53bit (&__libc_save_rm, (RM))
#define __nan(str) \
(__builtin_constant_p (str) && str[0] == '\0' ? NAN : __nan (str))
#define __nanf(str) \

View File

@ -59,10 +59,9 @@ __ieee754_exp(double x) {
int4 k;
#endif
int4 i,j,m,n,ex;
fenv_t env;
double retval;
libc_feholdexcept_setround (&env, FE_TONEAREST);
SET_RESTORE_ROUND (FE_TONEAREST);
junk1.x = x;
m = junk1.i[HIGH_HALF];
@ -157,7 +156,6 @@ __ieee754_exp(double x) {
else { retval = __slowexp(x); goto ret; }
}
ret:
libc_feupdateenv (&env);
return retval;
}
#ifndef __ieee754_exp

View File

@ -61,57 +61,56 @@ __ieee754_exp2 (double x)
int tval, unsafe;
double rx, x22, result;
union ieee754_double ex2_u, scale_u;
fenv_t oldenv;
libc_feholdexcept_setround (&oldenv, FE_TONEAREST);
{
SET_RESTORE_ROUND_NOEX (FE_TONEAREST);
/* 1. Argument reduction.
Choose integers ex, -256 <= t < 256, and some real
-1/1024 <= x1 <= 1024 so that
x = ex + t/512 + x1.
/* 1. Argument reduction.
Choose integers ex, -256 <= t < 256, and some real
-1/1024 <= x1 <= 1024 so that
x = ex + t/512 + x1.
First, calculate rx = ex + t/512. */
rx = x + THREEp42;
rx -= THREEp42;
x -= rx; /* Compute x=x1. */
/* Compute tval = (ex*512 + t)+256.
Now, t = (tval mod 512)-256 and ex=tval/512 [that's mod, NOT %; and
/-round-to-nearest not the usual c integer /]. */
tval = (int) (rx * 512.0 + 256.0);
First, calculate rx = ex + t/512. */
rx = x + THREEp42;
rx -= THREEp42;
x -= rx; /* Compute x=x1. */
/* Compute tval = (ex*512 + t)+256.
Now, t = (tval mod 512)-256 and ex=tval/512 [that's mod, NOT %;
and /-round-to-nearest not the usual c integer /]. */
tval = (int) (rx * 512.0 + 256.0);
/* 2. Adjust for accurate table entry.
Find e so that
x = ex + t/512 + e + x2
where -1e6 < e < 1e6, and
(double)(2^(t/512+e))
is accurate to one part in 2^-64. */
/* 2. Adjust for accurate table entry.
Find e so that
x = ex + t/512 + e + x2
where -1e6 < e < 1e6, and
(double)(2^(t/512+e))
is accurate to one part in 2^-64. */
/* 'tval & 511' is the same as 'tval%512' except that it's always
positive.
Compute x = x2. */
x -= exp2_deltatable[tval & 511];
/* 'tval & 511' is the same as 'tval%512' except that it's always
positive.
Compute x = x2. */
x -= exp2_deltatable[tval & 511];
/* 3. Compute ex2 = 2^(t/512+e+ex). */
ex2_u.d = exp2_accuratetable[tval & 511];
tval >>= 9;
unsafe = abs(tval) >= -DBL_MIN_EXP - 1;
ex2_u.ieee.exponent += tval >> unsafe;
scale_u.d = 1.0;
scale_u.ieee.exponent += tval - (tval >> unsafe);
/* 3. Compute ex2 = 2^(t/512+e+ex). */
ex2_u.d = exp2_accuratetable[tval & 511];
tval >>= 9;
unsafe = abs(tval) >= -DBL_MIN_EXP - 1;
ex2_u.ieee.exponent += tval >> unsafe;
scale_u.d = 1.0;
scale_u.ieee.exponent += tval - (tval >> unsafe);
/* 4. Approximate 2^x2 - 1, using a fourth-degree polynomial,
with maximum error in [-2^-10-2^-30,2^-10+2^-30]
less than 10^-19. */
/* 4. Approximate 2^x2 - 1, using a fourth-degree polynomial,
with maximum error in [-2^-10-2^-30,2^-10+2^-30]
less than 10^-19. */
x22 = (((.0096181293647031180
* x + .055504110254308625)
* x + .240226506959100583)
* x + .69314718055994495) * ex2_u.d;
math_opt_barrier (x22);
x22 = (((.0096181293647031180
* x + .055504110254308625)
* x + .240226506959100583)
* x + .69314718055994495) * ex2_u.d;
math_opt_barrier (x22);
}
/* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex). */
libc_fesetenv (&oldenv);
result = x22 * x + ex2_u.d;
if (!unsafe)

View File

@ -85,10 +85,9 @@ __ieee754_pow(double x, double y) {
(u.i[HIGH_HALF]==0 && u.i[LOW_HALF]!=0)) &&
/* 2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */
(v.i[HIGH_HALF]&0x7fffffff) < 0x4ff00000) { /* if y<-1 or y>1 */
fenv_t env;
double retval;
libc_feholdexcept_setround (&env, FE_TONEAREST);
SET_RESTORE_ROUND (FE_TONEAREST);
z = log1(x,&aa,&error); /* x^y =e^(y log (X)) */
t = y*134217729.0;
@ -105,7 +104,6 @@ __ieee754_pow(double x, double y) {
t = __exp1(a1,a2,1.9e16*error); /* return -10 or 0 if wasn't computed exactly */
retval = (t>0)?t:power1(x,y);
libc_feupdateenv (&env);
return retval;
}

View File

@ -108,10 +108,9 @@ __sin(double x){
#if 0
int4 nn;
#endif
fenv_t env;
double retval = 0;
libc_feholdexcept_setround_53bit (&env, FE_TONEAREST);
SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
u.x = x;
m = u.i[HIGH_HALF];
@ -365,7 +364,6 @@ __sin(double x){
}
ret:
libc_feupdateenv_53bit (&env);
return retval;
}
@ -383,10 +381,9 @@ __cos(double x)
mynumber u,v;
int4 k,m,n;
fenv_t env;
double retval = 0;
libc_feholdexcept_setround_53bit (&env, FE_TONEAREST);
SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
u.x = x;
m = u.i[HIGH_HALF];
@ -635,7 +632,6 @@ __cos(double x)
}
ret:
libc_feupdateenv_53bit (&env);
return retval;
}

View File

@ -68,13 +68,12 @@ tan(double x) {
mp_no mpy;
#endif
fenv_t env;
double retval;
int __branred(double, double *, double *);
int __mpranred(double, mp_no *, int);
libc_feholdexcept_setround_53bit (&env, FE_TONEAREST);
SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
/* x=+-INF, x=NaN */
num.d = x; ux = num.i[HIGH_HALF];
@ -503,7 +502,6 @@ tan(double x) {
goto ret;
ret:
libc_feupdateenv_53bit (&env);
return retval;
}

View File

@ -54,53 +54,52 @@ __ieee754_exp2f (float x)
int tval, unsafe;
float rx, x22, result;
union ieee754_float ex2_u, scale_u;
fenv_t oldenv;
libc_feholdexcept_setroundf (&oldenv, FE_TONEAREST);
{
SET_RESTORE_ROUND_NOEXF (FE_TONEAREST);
/* 1. Argument reduction.
Choose integers ex, -128 <= t < 128, and some real
-1/512 <= x1 <= 1/512 so that
x = ex + t/512 + x1.
/* 1. Argument reduction.
Choose integers ex, -128 <= t < 128, and some real
-1/512 <= x1 <= 1/512 so that
x = ex + t/512 + x1.
First, calculate rx = ex + t/256. */
rx = x + THREEp14;
rx -= THREEp14;
x -= rx; /* Compute x=x1. */
/* Compute tval = (ex*256 + t)+128.
Now, t = (tval mod 256)-128 and ex=tval/256 [that's mod, NOT %; and
/-round-to-nearest not the usual c integer /]. */
tval = (int) (rx * 256.0f + 128.0f);
First, calculate rx = ex + t/256. */
rx = x + THREEp14;
rx -= THREEp14;
x -= rx; /* Compute x=x1. */
/* Compute tval = (ex*256 + t)+128.
Now, t = (tval mod 256)-128 and ex=tval/256 [that's mod, NOT %;
and /-round-to-nearest not the usual c integer /]. */
tval = (int) (rx * 256.0f + 128.0f);
/* 2. Adjust for accurate table entry.
Find e so that
x = ex + t/256 + e + x2
where -7e-4 < e < 7e-4, and
(float)(2^(t/256+e))
is accurate to one part in 2^-64. */
/* 2. Adjust for accurate table entry.
Find e so that
x = ex + t/256 + e + x2
where -7e-4 < e < 7e-4, and
(float)(2^(t/256+e))
is accurate to one part in 2^-64. */
/* 'tval & 255' is the same as 'tval%256' except that it's always
positive.
Compute x = x2. */
x -= __exp2f_deltatable[tval & 255];
/* 'tval & 255' is the same as 'tval%256' except that it's always
positive.
Compute x = x2. */
x -= __exp2f_deltatable[tval & 255];
/* 3. Compute ex2 = 2^(t/255+e+ex). */
ex2_u.f = __exp2f_atable[tval & 255];
tval >>= 8;
unsafe = abs(tval) >= -FLT_MIN_EXP - 1;
ex2_u.ieee.exponent += tval >> unsafe;
scale_u.f = 1.0;
scale_u.ieee.exponent += tval - (tval >> unsafe);
/* 3. Compute ex2 = 2^(t/255+e+ex). */
ex2_u.f = __exp2f_atable[tval & 255];
tval >>= 8;
unsafe = abs(tval) >= -FLT_MIN_EXP - 1;
ex2_u.ieee.exponent += tval >> unsafe;
scale_u.f = 1.0;
scale_u.ieee.exponent += tval - (tval >> unsafe);
/* 4. Approximate 2^x2 - 1, using a second-degree polynomial,
with maximum error in [-2^-9 - 2^-14, 2^-9 + 2^-14]
less than 1.3e-10. */
/* 4. Approximate 2^x2 - 1, using a second-degree polynomial,
with maximum error in [-2^-9 - 2^-14, 2^-9 + 2^-14]
less than 1.3e-10. */
x22 = (.24022656679f * x + .69314736128f) * ex2_u.f;
x22 = (.24022656679f * x + .69314736128f) * ex2_u.f;
}
/* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex). */
libc_fesetenv (&oldenv);
result = x22 * x + ex2_u.f;
if (!unsafe)

View File

@ -80,40 +80,39 @@ __ieee754_expf (float x)
double x22, t, result, dx;
float n, delta;
union ieee754_double ex2_u;
fenv_t oldenv;
libc_feholdexcept_setroundf (&oldenv, FE_TONEAREST);
{
SET_RESTORE_ROUND_NOEXF (FE_TONEAREST);
/* Calculate n. */
n = x * M_1_LN2 + THREEp22;
n -= THREEp22;
dx = x - n*M_LN2;
/* Calculate n. */
n = x * M_1_LN2 + THREEp22;
n -= THREEp22;
dx = x - n*M_LN2;
/* Calculate t/512. */
t = dx + THREEp42;
t -= THREEp42;
dx -= t;
/* Calculate t/512. */
t = dx + THREEp42;
t -= THREEp42;
dx -= t;
/* Compute tval = t. */
tval = (int) (t * 512.0);
/* Compute tval = t. */
tval = (int) (t * 512.0);
if (t >= 0)
delta = - __exp_deltatable[tval];
else
delta = __exp_deltatable[-tval];
if (t >= 0)
delta = - __exp_deltatable[tval];
else
delta = __exp_deltatable[-tval];
/* Compute ex2 = 2^n e^(t/512+delta[t]). */
ex2_u.d = __exp_atable[tval+177];
ex2_u.ieee.exponent += (int) n;
/* Compute ex2 = 2^n e^(t/512+delta[t]). */
ex2_u.d = __exp_atable[tval+177];
ex2_u.ieee.exponent += (int) n;
/* Approximate e^(dx+delta) - 1, using a second-degree polynomial,
with maximum error in [-2^-10-2^-28,2^-10+2^-28]
less than 5e-11. */
x22 = (0.5000000496709180453 * dx + 1.0000001192102037084) * dx + delta;
/* Approximate e^(dx+delta) - 1, using a second-degree polynomial,
with maximum error in [-2^-10-2^-28,2^-10+2^-28]
less than 5e-11. */
x22 = (0.5000000496709180453 * dx + 1.0000001192102037084) * dx + delta;
}
/* Return result. */
libc_fesetenvf (&oldenv);
result = x22 * ex2_u.d + ex2_u.d;
return (float) result;
}

View File

@ -119,6 +119,29 @@ libc_feupdateenv (fenv_t *e)
#define libc_feupdateenv libc_feupdateenv
#define libc_feupdateenvf libc_feupdateenv
static __always_inline void
libc_feholdsetround (fenv_t *e, int r)
{
unsigned int mxcsr;
asm (STMXCSR " %0" : "=m" (*&mxcsr));
e->__mxcsr = mxcsr;
mxcsr = (mxcsr & ~0x6000) | (r << 3);
asm volatile (LDMXCSR " %0" : : "m" (*&mxcsr));
}
#define libc_feholdsetround libc_feholdsetround
#define libc_feholdsetroundf libc_feholdsetround
static __always_inline void
libc_feresetround (fenv_t *e)
{
unsigned int mxcsr;
asm (STMXCSR " %0" : "=m" (*&mxcsr));
mxcsr = (mxcsr & ~0x6000) | (e->__mxcsr & 0x6000);
asm volatile (LDMXCSR " %0" : : "m" (*&mxcsr));
}
#define libc_feresetround libc_feresetround
#define libc_feresetroundf libc_feresetround
#include_next <math_private.h>
extern __always_inline double