Merge changes between r17530 and r17700 from /fsf/trunk.

git-svn-id: svn://svn.eglibc.org/trunk@17701 7b3dc134-2b1b-0410-93df-9e9f96275f8d

16 files changed, 216 insertions, 253 deletions

diff --git a/libc/sysdeps/ieee754/dbl-64/e_atan2.c b/libc/sysdeps/ieee754/dbl-64/e_atan2.c
index dcef55f07..497afcab1 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_atan2.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_atan2.c
@@ -1,7 +1,7 @@
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
- * Copyright (C) 2001, 2011 Free Software Foundation
+ * Copyright (C) 2001-2012 Free Software Foundation
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published by
@@ -153,6 +153,13 @@ __ieee754_atan2(double y,double x) {
   /* if either x or y is extremely close to zero, scale abs(x), abs(y). */
   if (ax<twom500.d || ay<twom500.d) { ax*=two500.d;  ay*=two500.d; }
 
+  /* Likewise for large x and y.  */
+  if (ax > two500.d || ay > two500.d)
+    {
+      ax *= twom500.d;
+      ay *= twom500.d;
+    }
+
   /* x,y which are neither special nor extreme */
   if (ay<ax) {
     u=ay/ax;
diff --git a/libc/sysdeps/ieee754/dbl-64/e_exp.c b/libc/sysdeps/ieee754/dbl-64/e_exp.c
index cb8d9e8d9..5deba5e44 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_exp.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_exp.c
@@ -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
diff --git a/libc/sysdeps/ieee754/dbl-64/e_exp2.c b/libc/sysdeps/ieee754/dbl-64/e_exp2.c
index 4cf879b7f..e57ec9211 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_exp2.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_exp2.c
@@ -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);
-
-      /* 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);
-
-      /* 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];
-
-      /* 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.  */
-
-      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);
+      {
+	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.
+
+	   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.  */
+
+	/* '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);
+
+	/* 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);
+      }
 
+      /* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex).  */
       result = x22 * x + ex2_u.d;
 
       if (!unsafe)
diff --git a/libc/sysdeps/ieee754/dbl-64/e_pow.c b/libc/sysdeps/ieee754/dbl-64/e_pow.c
index 550633cf9..26ffaaddc 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_pow.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_pow.c
@@ -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;
   }
 
@@ -113,7 +111,7 @@ __ieee754_pow(double x, double y) {
     if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0)
 	|| (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000)
       return y;
-    if (ABS(y) > 1.0e20) return (y>0)?0:INF.x;
+    if (ABS(y) > 1.0e20) return (y>0)?0:1.0/ABS(x);
     k = checkint(y);
     if (k == -1)
       return y < 0 ? 1.0/x : x;
diff --git a/libc/sysdeps/ieee754/dbl-64/k_rem_pio2.c b/libc/sysdeps/ieee754/dbl-64/k_rem_pio2.c
index a1e0c6d72..fcf956afb 100644
--- a/libc/sysdeps/ieee754/dbl-64/k_rem_pio2.c
+++ b/libc/sysdeps/ieee754/dbl-64/k_rem_pio2.c
@@ -273,13 +273,16 @@ recompute:
 		y[0] = (ih==0)? fw: -fw;
 		break;
 	    case 1:
-	    case 2:
-		fw = 0.0;
-		for (i=jz;i>=0;i--) fw += fq[i];
-		y[0] = (ih==0)? fw: -fw;
-		fw = fq[0]-fw;
-		for (i=1;i<=jz;i++) fw += fq[i];
-		y[1] = (ih==0)? fw: -fw;
+	    case 2:;
+#if __FLT_EVAL_METHOD__ != 0
+		volatile
+#endif
+		double fv = 0.0;
+		for (i=jz;i>=0;i--) fv += fq[i];
+		y[0] = (ih==0)? fv: -fv;
+		fv = fq[0]-fv;
+		for (i=1;i<=jz;i++) fv += fq[i];
+		y[1] = (ih==0)? fv: -fv;
 		break;
 	    case 3:	/* painful */
 		for (i=jz;i>0;i--) {
diff --git a/libc/sysdeps/ieee754/dbl-64/s_fma.c b/libc/sysdeps/ieee754/dbl-64/s_fma.c
index a27e246a4..ab20a801a 100644
--- a/libc/sysdeps/ieee754/dbl-64/s_fma.c
+++ b/libc/sysdeps/ieee754/dbl-64/s_fma.c
@@ -149,35 +149,36 @@ __fma (double x, double y, double z)
 
   fenv_t env;
   libc_feholdexcept_setround (&env, FE_TOWARDZERO);
+
   /* Perform m2 + a2 addition with round to odd.  */
   u.d = a2 + m2;
 
+  if (__builtin_expect (adjust < 0, 0))
+    {
+      if ((u.ieee.mantissa1 & 1) == 0)
+	u.ieee.mantissa1 |= libc_fetestexcept (FE_INEXACT) != 0;
+      v.d = a1 + u.d;
+    }
+
+  /* Reset rounding mode and test for inexact simultaneously.  */
+  int j = libc_feupdateenv_test (&env, FE_INEXACT) != 0;
+
   if (__builtin_expect (adjust == 0, 1))
     {
       if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7ff)
-	u.ieee.mantissa1 |= libc_fetestexcept (FE_INEXACT) != 0;
-      libc_feupdateenv (&env);
+	u.ieee.mantissa1 |= j;
       /* 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 |= libc_fetestexcept (FE_INEXACT) != 0;
-      libc_feupdateenv (&env);
+	u.ieee.mantissa1 |= j;
       /* Result is a1 + u.d, scaled up.  */
       return (a1 + u.d) * 0x1p53;
     }
   else
     {
-      if ((u.ieee.mantissa1 & 1) == 0)
-	u.ieee.mantissa1 |= libc_fetestexcept (FE_INEXACT) != 0;
-      v.d = a1 + u.d;
-      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));
       /* If a1 + u.d is exact, the only rounding happens during
 	 scaling down.  */
       if (j == 0)
diff --git a/libc/sysdeps/ieee754/dbl-64/s_fmaf.c b/libc/sysdeps/ieee754/dbl-64/s_fmaf.c
index 00cd38270..7a939aaed 100644
--- a/libc/sysdeps/ieee754/dbl-64/s_fmaf.c
+++ b/libc/sysdeps/ieee754/dbl-64/s_fmaf.c
@@ -35,12 +35,18 @@ __fmaf (float x, float y, float z)
   /* Multiplication is always exact.  */
   double temp = (double) x * (double) y;
   union ieee754_double u;
-  libc_feholdexcept_setroundf (&env, FE_TOWARDZERO);
+
+  libc_feholdexcept_setround (&env, FE_TOWARDZERO);
+
   /* Perform addition with round to odd.  */
   u.d = temp + (double) z;
+
+  /* Reset rounding mode and test for inexact simultaneously.  */
+  int j = libc_feupdateenv_test (&env, FE_INEXACT) != 0;
+
   if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7ff)
-    u.ieee.mantissa1 |= libc_fetestexcept (FE_INEXACT) != 0;
-  libc_feupdateenv (&env);
+    u.ieee.mantissa1 |= j;
+
   /* And finally truncation with round to nearest.  */
   return (float) u.d;
 }
diff --git a/libc/sysdeps/ieee754/dbl-64/s_nearbyint.c b/libc/sysdeps/ieee754/dbl-64/s_nearbyint.c
index d780150a9..60afafded 100644
--- a/libc/sysdeps/ieee754/dbl-64/s_nearbyint.c
+++ b/libc/sysdeps/ieee754/dbl-64/s_nearbyint.c
@@ -38,18 +38,12 @@ double __nearbyint(double x)
 {
 	fenv_t env;
 	int32_t i0,j0,sx;
-	u_int32_t i,i1;
 	double w,t;
-	EXTRACT_WORDS(i0,i1,x);
+	GET_HIGH_WORD(i0,x);
 	sx = (i0>>31)&1;
 	j0 = ((i0>>20)&0x7ff)-0x3ff;
-	if(j0<20) {
+	if(j0<52) {
 	    if(j0<0) {
-		if(((i0&0x7fffffff)|i1)==0) return x;
-		i1 |= (i0&0x0fffff);
-		i0 &= 0xfffe0000;
-		i0 |= ((i1|-i1)>>12)&0x80000;
-		SET_HIGH_WORD(x,i0);
 		feholdexcept (&env);
 	        w = TWO52[sx]+x;
 	        t =  w-TWO52[sx];
@@ -57,32 +51,11 @@ double __nearbyint(double x)
 		GET_HIGH_WORD(i0,t);
 		SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
 	        return t;
-	    } else {
-		i = (0x000fffff)>>j0;
-		if(((i0&i)|i1)==0) return x; /* x is integral */
-		i>>=1;
-		if(((i0&i)|i1)!=0) {
-		  if (j0==19)
-		    i1 = 0x40000000;
-		  else if (j0<18)
-		    i0 = (i0&(~i))|((0x20000)>>j0);
-		  else
-		    {
-		      i0 &= ~i;
-		      i1 = 0x80000000;
-		    }
-		}
 	    }
-	} else if (j0>51) {
+	} else {
 	    if(j0==0x400) return x+x;	/* inf or NaN */
 	    else return x;		/* x is integral */
-	} else {
-	    i = ((u_int32_t)(0xffffffff))>>(j0-20);
-	    if((i1&i)==0) return x;	/* x is integral */
-	    i>>=1;
-	    if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
 	}
-	INSERT_WORDS(x,i0,i1);
 	feholdexcept (&env);
 	w = TWO52[sx]+x;
 	t = w-TWO52[sx];
diff --git a/libc/sysdeps/ieee754/dbl-64/s_rint.c b/libc/sysdeps/ieee754/dbl-64/s_rint.c
index c45d7848d..845890916 100644
--- a/libc/sysdeps/ieee754/dbl-64/s_rint.c
+++ b/libc/sysdeps/ieee754/dbl-64/s_rint.c
@@ -33,49 +33,22 @@ double
 __rint(double x)
 {
 	int32_t i0,j0,sx;
-	u_int32_t i,i1;
 	double w,t;
-	EXTRACT_WORDS(i0,i1,x);
+	GET_HIGH_WORD(i0,x);
 	sx = (i0>>31)&1;
 	j0 = ((i0>>20)&0x7ff)-0x3ff;
-	if(j0<20) {
+	if(j0<52) {
 	    if(j0<0) {
-		if(((i0&0x7fffffff)|i1)==0) return x;
-		i1 |= (i0&0x0fffff);
-		i0 &= 0xfffe0000;
-		i0 |= ((i1|-i1)>>12)&0x80000;
-		SET_HIGH_WORD(x,i0);
 		w = TWO52[sx]+x;
 		t =  w-TWO52[sx];
 		GET_HIGH_WORD(i0,t);
 		SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
 		return t;
-	    } else {
-		i = (0x000fffff)>>j0;
-		if(((i0&i)|i1)==0) return x; /* x is integral */
-		i>>=1;
-		if(((i0&i)|i1)!=0) {
-		    if (j0==19)
-		      i1 = 0x40000000;
-		    else if (j0<18)
-		      i0 = (i0&(~i))|((0x20000)>>j0);
-		    else
-		      {
-			i0 &= ~i;
-			i1 = 0x80000000;
-		      }
-		}
 	    }
-	} else if (j0>51) {
+	} else {
 	    if(j0==0x400) return x+x;	/* inf or NaN */
 	    else return x;		/* x is integral */
-	} else {
-	    i = ((u_int32_t)(0xffffffff))>>(j0-20);
-	    if((i1&i)==0) return x;	/* x is integral */
-	    i>>=1;
-	    if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
 	}
-	INSERT_WORDS(x,i0,i1);
 	w = TWO52[sx]+x;
 	return w-TWO52[sx];
 }
diff --git a/libc/sysdeps/ieee754/dbl-64/s_sin.c b/libc/sysdeps/ieee754/dbl-64/s_sin.c
index e3e3a2a96..7b9252f81 100644
--- a/libc/sysdeps/ieee754/dbl-64/s_sin.c
+++ b/libc/sysdeps/ieee754/dbl-64/s_sin.c
@@ -1,7 +1,7 @@
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
- * Copyright (C) 2001, 2009, 2011 Free Software Foundation
+ * Copyright (C) 2001-2012 Free Software Foundation
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published by
@@ -108,10 +108,9 @@ __sin(double x){
 #if 0
 	int4 nn;
 #endif
-	fenv_t env;
 	double retval = 0;
 
-	libc_feholdexcept_setround (&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 (&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 (&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 (&env);
   return retval;
 }
 
diff --git a/libc/sysdeps/ieee754/dbl-64/s_tan.c b/libc/sysdeps/ieee754/dbl-64/s_tan.c
index acff67c98..f8507eaa4 100644
--- a/libc/sysdeps/ieee754/dbl-64/s_tan.c
+++ b/libc/sysdeps/ieee754/dbl-64/s_tan.c
@@ -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 (&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 (&env);
   return retval;
 }
 
diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/math_private.h b/libc/sysdeps/ieee754/dbl-64/wordsize-64/math_private.h
new file mode 100644
index 000000000..b66085eb1
--- /dev/null
+++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/math_private.h
@@ -0,0 +1,35 @@
+#ifndef _MATH_PRIVATE_H_
+
+#include_next <math_private.h>
+
+#ifndef __isnan
+extern __always_inline int
+__isnan (double d)
+{
+  uint64_t di;
+  EXTRACT_WORDS64 (di, d);
+  return (di & 0x7fffffffffffffffull) > 0x7ff0000000000000ull;
+}
+#endif
+
+#ifndef __isinf_ns
+extern __always_inline int
+__isinf_ns (double d)
+{
+  uint64_t di;
+  EXTRACT_WORDS64 (di, d);
+  return (di & 0x7fffffffffffffffull) == 0x7ff0000000000000ull;
+}
+#endif
+
+#ifndef __finite
+extern __always_inline int
+__finite (double d)
+{
+  uint64_t di;
+  EXTRACT_WORDS64 (di, d);
+  return (di & 0x7fffffffffffffffull) < 0x7ff0000000000000ull;
+}
+#endif
+
+#endif /* _MATH_PRIVATE_H_ */
diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c
index f7c9ea568..4a72005d0 100644
--- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c
+++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c
@@ -1 +1,79 @@
-/* The code is the same as lround.  Use an alias, see l_round.c.  */
+/* Round double value to long long int.
+   Copyright (C) 1997-2012 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+#define lround __hidden_lround
+#define __lround __hidden___lround
+
+#include <math.h>
+
+#include <math_private.h>
+
+
+long long int
+__llround (double x)
+{
+  int32_t j0;
+  int64_t i0;
+  long long int result;
+  int sign;
+
+  EXTRACT_WORDS64 (i0, x);
+  j0 = ((i0 >> 52) & 0x7ff) - 0x3ff;
+  sign = i0 < 0 ? -1 : 1;
+  i0 &= UINT64_C(0xfffffffffffff);
+  i0 |= UINT64_C(0x10000000000000);
+
+  if (j0 < (int32_t) (8 * sizeof (long long int)) - 1)
+    {
+      if (j0 < 0)
+	return j0 < -1 ? 0 : sign;
+      else if (j0 >= 52)
+	result = i0 << (j0 - 52);
+      else
+	{
+	  i0 += UINT64_C(0x8000000000000) >> j0;
+
+	  result = i0 >> (52 - j0);
+	}
+    }
+  else
+    {
+      /* The number is too large.  It is left implementation defined
+	 what happens.  */
+      return (long long int) x;
+    }
+
+  return sign * result;
+}
+
+weak_alias (__llround, llround)
+#ifdef NO_LONG_DOUBLE
+strong_alias (__llround, __lroundl)
+weak_alias (__llround, lroundl)
+#endif
+
+/* long has the same width as long long on 64-bit machines.  */
+#undef lround
+#undef __lround
+strong_alias (__llround, __lround)
+weak_alias (__llround, lround)
+#ifdef NO_LONG_DOUBLE
+strong_alias (__llround, __llroundl)
+weak_alias (__llround, llroundl)
+#endif
diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c
index e3a1db0d5..1e06fcb16 100644
--- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c
+++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c
@@ -1,79 +1 @@
-/* Round double value to long int.
-   Copyright (C) 1997, 2004, 2009 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
-
-   The GNU C Library is free software; you can redistribute it and/or
-   modify it under the terms of the GNU Lesser General Public
-   License as published by the Free Software Foundation; either
-   version 2.1 of the License, or (at your option) any later version.
-
-   The GNU C Library is distributed in the hope that it will be useful,
-   but WITHOUT ANY WARRANTY; without even the implied warranty of
-   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-   Lesser General Public License for more details.
-
-   You should have received a copy of the GNU Lesser General Public
-   License along with the GNU C Library; if not, see
-   <http://www.gnu.org/licenses/>.  */
-
-#define llround __hidden_llround
-#define __llround __hidden___llround
-
-#include <math.h>
-
-#include <math_private.h>
-
-
-long int
-__lround (double x)
-{
-  int32_t j0;
-  int64_t i0;
-  long int result;
-  int sign;
-
-  EXTRACT_WORDS64 (i0, x);
-  j0 = ((i0 >> 52) & 0x7ff) - 0x3ff;
-  sign = i0 < 0 ? -1 : 1;
-  i0 &= UINT64_C(0xfffffffffffff);
-  i0 |= UINT64_C(0x10000000000000);
-
-  if (j0 < (int32_t) (8 * sizeof (long int)) - 1)
-    {
-      if (j0 < 0)
-	return j0 < -1 ? 0 : sign;
-      else if (j0 >= 52)
-	result = i0 << (j0 - 52);
-      else
-	{
-	  i0 += UINT64_C(0x8000000000000) >> j0;
-
-	  result = i0 >> (52 - j0);
-	}
-    }
-  else
-    {
-      /* The number is too large.  It is left implementation defined
-	 what happens.  */
-      return (long int) x;
-    }
-
-  return sign * result;
-}
-
-weak_alias (__lround, lround)
-#ifdef NO_LONG_DOUBLE
-strong_alias (__lround, __lroundl)
-weak_alias (__lround, lroundl)
-#endif
-
-/* long long has the same width as long on 64-bit machines.  */
-#undef llround
-#undef __llround
-strong_alias (__lround, __llround)
-weak_alias (__lround, llround)
-#ifdef NO_LONG_DOUBLE
-strong_alias (__lround, __llroundl)
-weak_alias (__lround, llroundl)
-#endif
+/* The code is the same as llround.  Use an alias, see s_llround.c.  */
diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c
index a3809984e..a58a6202e 100644
--- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c
+++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c
@@ -41,29 +41,17 @@ __nearbyint(double x)
 	j0 = ((i0>>52)&0x7ff)-0x3ff;
 	if(__builtin_expect(j0<52, 1)) {
 	    if(j0<0) {
-	      if((i0&UINT64_C(0x7fffffffffffffff))==0) return x;
-		uint64_t i = i0 & UINT64_C(0xfffffffffffff);
-		i0 &= UINT64_C(0xfffe000000000000);
-		i0 |= (((i|-i) >> 12) & UINT64_C(0x8000000000000));
-		INSERT_WORDS64(x,i0);
 		libc_feholdexcept (&env);
 		double w = TWO52[sx]+x;
 		double t =  w-TWO52[sx];
 		math_opt_barrier(t);
 		libc_fesetenv (&env);
 		return copysign(t, x);
-	    } else {
-		uint64_t i = UINT64_C(0x000fffffffffffff)>>j0;
-		if((i0&i)==0) return x; /* x is integral */
-		i>>=1;
-		if((i0&i)!=0)
-		    i0 = (i0&(~i))|(UINT64_C(0x4000000000000)>>j0);
 	    }
 	} else {
 	    if(j0==0x400) return x+x;	/* inf or NaN */
 	    else return x;		/* x is integral */
 	}
-	INSERT_WORDS64(x,i0);
 	libc_feholdexcept (&env);
 	double w = TWO52[sx]+x;
 	double t = w-TWO52[sx];
diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c
index 20eef3c5e..87b2339d4 100644
--- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c
+++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c
@@ -38,28 +38,16 @@ __rint(double x)
 	j0 = ((i0>>52)&0x7ff)-0x3ff;
 	if(j0<52) {
 	    if(j0<0) {
-		if((i0 & UINT64_C(0x7fffffffffffffff))==0) return x;
-		uint64_t i = i0 & UINT64_C(0xfffffffffffff);
-		i0 &= UINT64_C(0xfffe000000000000);
-		i0 |= (((i|-i) >> 12) & UINT64_C(0x8000000000000));
-		INSERT_WORDS64(x,i0);
 		double w = TWO52[sx]+x;
 		double t =  w-TWO52[sx];
 		EXTRACT_WORDS64(i0,t);
 		INSERT_WORDS64(t,(i0&UINT64_C(0x7fffffffffffffff))|(sx<<63));
 		return t;
-	    } else {
-		uint64_t i = UINT64_C(0x000fffffffffffff)>>j0;
-		if((i0&i)==0) return x; /* x is integral */
-		i>>=1;
-		if((i0&i)!=0)
-		    i0 = (i0&(~i))|(UINT64_C(0x4000000000000)>>j0);
 	    }
 	} else {
 	    if(j0==0x400) return x+x;	/* inf or NaN */
 	    else return x;		/* x is integral */
 	}
-	INSERT_WORDS64(x,i0);
 	double w = TWO52[sx]+x;
 	return w-TWO52[sx];
 }