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

2 files changed, 68 insertions, 70 deletions

diff --git a/sysdeps/ieee754/flt-32/e_exp2f.c b/sysdeps/ieee754/flt-32/e_exp2f.c
index e728e6ec74..267d81b23f 100644
--- a/sysdeps/ieee754/flt-32/e_exp2f.c
+++ b/sysdeps/ieee754/flt-32/e_exp2f.c
@@ -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);
-
-      /* 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);
-
-      /* 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];
-
-      /* 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.  */
-
-      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);
+      {
+	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.
+
+	   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.  */
+
+	/* '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);
+
+	/* 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;
+      }
 
+      /* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex).  */
       result = x22 * x + ex2_u.f;
 
       if (!unsafe)
diff --git a/sysdeps/ieee754/flt-32/e_expf.c b/sysdeps/ieee754/flt-32/e_expf.c
index e69e7f6ae0..57aff16ab7 100644
--- a/sysdeps/ieee754/flt-32/e_expf.c
+++ b/sysdeps/ieee754/flt-32/e_expf.c
@@ -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;
     }