Remove dbl-64/wordsize-64 (part 2)

Remove the wordsize-64 implementations by merging them into the main dbl-64
directory.  The second patch just moves all wordsize-64 files and removes a
few wordsize-64 uses in comments and Implies files.

Reviewed-by: Adhemerval Zanella  <adhemerval.zanella@linaro.org>

34 files changed, 422 insertions, 1837 deletions

diff --git a/sysdeps/ieee754/dbl-64/e_acosh.c b/sysdeps/ieee754/dbl-64/e_acosh.c
index 75df0ab5ef..a241366f30 100644
--- a/sysdeps/ieee754/dbl-64/e_acosh.c
+++ b/sysdeps/ieee754/dbl-64/e_acosh.c
@@ -1,4 +1,4 @@
-/* @(#)e_acosh.c 5.1 93/09/24 */
+/* Optimized for 64-bit by Ulrich Drepper <drepper@gmail.com>, 2012 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -29,42 +29,40 @@
 #include <libm-alias-finite.h>
 
 static const double
-  one = 1.0,
-  ln2 = 6.93147180559945286227e-01;    /* 0x3FE62E42, 0xFEFA39EF */
+one	= 1.0,
+ln2	= 6.93147180559945286227e-01;  /* 0x3FE62E42, 0xFEFA39EF */
 
 double
 __ieee754_acosh (double x)
 {
-  double t;
-  int32_t hx;
-  uint32_t lx;
-  EXTRACT_WORDS (hx, lx, x);
-  if (hx < 0x3ff00000)                  /* x < 1 */
-    {
-      return (x - x) / (x - x);
-    }
-  else if (hx >= 0x41b00000)            /* x > 2**28 */
+  int64_t hx;
+  EXTRACT_WORDS64 (hx, x);
+
+  if (hx > INT64_C (0x4000000000000000))
     {
-      if (hx >= 0x7ff00000)             /* x is inf of NaN */
+      if (__glibc_unlikely (hx >= INT64_C (0x41b0000000000000)))
 	{
-	  return x + x;
+	  /* x > 2**28 */
+	  if (hx >= INT64_C (0x7ff0000000000000))
+	    /* x is inf of NaN */
+	    return x + x;
+	  else
+	    return __ieee754_log (x) + ln2;/* acosh(huge)=log(2x) */
 	}
-      else
-	return __ieee754_log (x) + ln2;         /* acosh(huge)=log(2x) */
-    }
-  else if (((hx - 0x3ff00000) | lx) == 0)
-    {
-      return 0.0;                       /* acosh(1) = 0 */
-    }
-  else if (hx > 0x40000000)             /* 2**28 > x > 2 */
-    {
-      t = x * x;
+
+      /* 2**28 > x > 2 */
+      double t = x * x;
       return __ieee754_log (2.0 * x - one / (x + sqrt (t - one)));
     }
-  else                                  /* 1<x<2 */
+  else if (__glibc_likely (hx > INT64_C (0x3ff0000000000000)))
     {
-      t = x - one;
+      /* 1<x<2 */
+      double t = x - one;
       return __log1p (t + sqrt (2.0 * t + t * t));
     }
+  else if (__glibc_likely (hx == INT64_C (0x3ff0000000000000)))
+    return 0.0;				/* acosh(1) = 0 */
+  else					/* x < 1 */
+    return (x - x) / (x - x);
 }
 libm_alias_finite (__ieee754_acosh, __acosh)
diff --git a/sysdeps/ieee754/dbl-64/e_cosh.c b/sysdeps/ieee754/dbl-64/e_cosh.c
index 6c78a3a4e9..4f41ca2c92 100644
--- a/sysdeps/ieee754/dbl-64/e_cosh.c
+++ b/sysdeps/ieee754/dbl-64/e_cosh.c
@@ -32,59 +32,54 @@
  */
 
 #include <math.h>
-#include <math-narrow-eval.h>
 #include <math_private.h>
 #include <libm-alias-finite.h>
 
-static const double one = 1.0, half = 0.5, huge = 1.0e300;
+static const double one = 1.0, half=0.5, huge = 1.0e300;
 
 double
 __ieee754_cosh (double x)
 {
-  double t, w;
-  int32_t ix;
-  uint32_t lx;
+	double t,w;
+	int32_t ix;
 
-  /* High word of |x|. */
-  GET_HIGH_WORD (ix, x);
-  ix &= 0x7fffffff;
+    /* High word of |x|. */
+	GET_HIGH_WORD(ix,x);
+	ix &= 0x7fffffff;
 
-  /* |x| in [0,22] */
-  if (ix < 0x40360000)
-    {
-      /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
-      if (ix < 0x3fd62e43)
-	{
-	  if (ix < 0x3c800000)
-	    return one;                                   /* cosh(tiny) = 1 */
-	  t = __expm1 (fabs (x));
-	  w = one + t;
-	  return one + (t * t) / (w + w);
-	}
+    /* |x| in [0,22] */
+	if (ix < 0x40360000) {
+	    /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
+		if(ix<0x3fd62e43) {
+		    if (ix<0x3c800000)			/* cosh(tiny) = 1 */
+		      return one;
+		    t = __expm1(fabs(x));
+		    w = one+t;
+		    return one+(t*t)/(w+w);
+		}
 
-      /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
-      t = __ieee754_exp (fabs (x));
-      return half * t + half / t;
-    }
+	    /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
+		t = __ieee754_exp(fabs(x));
+		return half*t+half/t;
+	}
 
-  /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
-  if (ix < 0x40862e42)
-    return half * __ieee754_exp (fabs (x));
+    /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
+	if (ix < 0x40862e42)  return half*__ieee754_exp(fabs(x));
 
-  /* |x| in [log(maxdouble), overflowthresold] */
-  GET_LOW_WORD (lx, x);
-  if (ix < 0x408633ce || ((ix == 0x408633ce) && (lx <= (uint32_t) 0x8fb9f87d)))
-    {
-      w = __ieee754_exp (half * fabs (x));
-      t = half * w;
-      return t * w;
-    }
+    /* |x| in [log(maxdouble), overflowthresold] */
+	int64_t fix;
+	EXTRACT_WORDS64(fix, x);
+	fix &= UINT64_C(0x7fffffffffffffff);
+	if (fix <= UINT64_C(0x408633ce8fb9f87d)) {
+	    w = __ieee754_exp(half*fabs(x));
+	    t = half*w;
+	    return t*w;
+	}
 
-  /* x is INF or NaN */
-  if (ix >= 0x7ff00000)
-    return x * x;
+    /* x is INF or NaN */
+	if(ix>=0x7ff00000) return x*x;
 
-  /* |x| > overflowthresold, cosh(x) overflow */
-  return math_narrow_eval (huge * huge);
+    /* |x| > overflowthresold, cosh(x) overflow */
+	return huge*huge;
 }
 libm_alias_finite (__ieee754_cosh, __cosh)
diff --git a/sysdeps/ieee754/dbl-64/e_fmod.c b/sysdeps/ieee754/dbl-64/e_fmod.c
index f6a095ba82..52a8687448 100644
--- a/sysdeps/ieee754/dbl-64/e_fmod.c
+++ b/sysdeps/ieee754/dbl-64/e_fmod.c
@@ -1,3 +1,4 @@
+/* Rewritten for 64-bit machines by Ulrich Drepper <drepper@gmail.com>.  */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -17,158 +18,89 @@
 
 #include <math.h>
 #include <math_private.h>
+#include <stdint.h>
 #include <libm-alias-finite.h>
 
-static const double one = 1.0, Zero[] = { 0.0, -0.0, };
+static const double one = 1.0, Zero[] = {0.0, -0.0,};
 
 double
 __ieee754_fmod (double x, double y)
 {
-  int32_t n, hx, hy, hz, ix, iy, sx, i;
-  uint32_t lx, ly, lz;
+	int32_t n,ix,iy;
+	int64_t hx,hy,hz,sx,i;
 
-  EXTRACT_WORDS (hx, lx, x);
-  EXTRACT_WORDS (hy, ly, y);
-  sx = hx & 0x80000000;                 /* sign of x */
-  hx ^= sx;                     /* |x| */
-  hy &= 0x7fffffff;             /* |y| */
+	EXTRACT_WORDS64(hx,x);
+	EXTRACT_WORDS64(hy,y);
+	sx = hx&UINT64_C(0x8000000000000000);	/* sign of x */
+	hx ^=sx;				/* |x| */
+	hy &= UINT64_C(0x7fffffffffffffff);	/* |y| */
 
-  /* purge off exception values */
-  if ((hy | ly) == 0 || (hx >= 0x7ff00000) ||   /* y=0,or x not finite */
-      ((hy | ((ly | -ly) >> 31)) > 0x7ff00000)) /* or y is NaN */
-    return (x * y) / (x * y);
-  if (hx <= hy)
-    {
-      if ((hx < hy) || (lx < ly))
-	return x;                               /* |x|<|y| return x */
-      if (lx == ly)
-	return Zero[(uint32_t) sx >> 31];      /* |x|=|y| return x*0*/
-    }
-
-  /* determine ix = ilogb(x) */
-  if (__glibc_unlikely (hx < 0x00100000))                  /* subnormal x */
-    {
-      if (hx == 0)
-	{
-	  for (ix = -1043, i = lx; i > 0; i <<= 1)
-	    ix -= 1;
-	}
-      else
-	{
-	  for (ix = -1022, i = (hx << 11); i > 0; i <<= 1)
-	    ix -= 1;
+    /* purge off exception values */
+	if(__builtin_expect(hy==0
+			    || hx >= UINT64_C(0x7ff0000000000000)
+			    || hy > UINT64_C(0x7ff0000000000000), 0))
+	  /* y=0,or x not finite or y is NaN */
+	    return (x*y)/(x*y);
+	if(__builtin_expect(hx<=hy, 0)) {
+	    if(hx<hy) return x;	/* |x|<|y| return x */
+	    return Zero[(uint64_t)sx>>63];	/* |x|=|y| return x*0*/
 	}
-    }
-  else
-    ix = (hx >> 20) - 1023;
 
-  /* determine iy = ilogb(y) */
-  if (__glibc_unlikely (hy < 0x00100000))                  /* subnormal y */
-    {
-      if (hy == 0)
-	{
-	  for (iy = -1043, i = ly; i > 0; i <<= 1)
-	    iy -= 1;
-	}
-      else
-	{
-	  for (iy = -1022, i = (hy << 11); i > 0; i <<= 1)
-	    iy -= 1;
-	}
-    }
-  else
-    iy = (hy >> 20) - 1023;
+    /* determine ix = ilogb(x) */
+	if(__builtin_expect(hx<UINT64_C(0x0010000000000000), 0)) {
+	  /* subnormal x */
+	  for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
+	} else ix = (hx>>52)-1023;
 
-  /* set up {hx,lx}, {hy,ly} and align y to x */
-  if (__glibc_likely (ix >= -1022))
-    hx = 0x00100000 | (0x000fffff & hx);
-  else                  /* subnormal x, shift x to normal */
-    {
-      n = -1022 - ix;
-      if (n <= 31)
-	{
-	  hx = (hx << n) | (lx >> (32 - n));
-	  lx <<= n;
-	}
-      else
-	{
-	  hx = lx << (n - 32);
-	  lx = 0;
-	}
-    }
-  if (__glibc_likely (iy >= -1022))
-    hy = 0x00100000 | (0x000fffff & hy);
-  else                  /* subnormal y, shift y to normal */
-    {
-      n = -1022 - iy;
-      if (n <= 31)
-	{
-	  hy = (hy << n) | (ly >> (32 - n));
-	  ly <<= n;
-	}
-      else
-	{
-	  hy = ly << (n - 32);
-	  ly = 0;
-	}
-    }
+    /* determine iy = ilogb(y) */
+	if(__builtin_expect(hy<UINT64_C(0x0010000000000000), 0)) {	/* subnormal y */
+	  for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
+	} else iy = (hy>>52)-1023;
 
-  /* fix point fmod */
-  n = ix - iy;
-  while (n--)
-    {
-      hz = hx - hy; lz = lx - ly; if (lx < ly)
-	hz -= 1;
-      if (hz < 0)
-	{
-	  hx = hx + hx + (lx >> 31); lx = lx + lx;
+    /* set up hx, hy and align y to x */
+	if(__builtin_expect(ix >= -1022, 1))
+	    hx = UINT64_C(0x0010000000000000)|(UINT64_C(0x000fffffffffffff)&hx);
+	else {		/* subnormal x, shift x to normal */
+	    n = -1022-ix;
+	    hx<<=n;
 	}
-      else
-	{
-	  if ((hz | lz) == 0)           /* return sign(x)*0 */
-	    return Zero[(uint32_t) sx >> 31];
-	  hx = hz + hz + (lz >> 31); lx = lz + lz;
+	if(__builtin_expect(iy >= -1022, 1))
+	    hy = UINT64_C(0x0010000000000000)|(UINT64_C(0x000fffffffffffff)&hy);
+	else {		/* subnormal y, shift y to normal */
+	    n = -1022-iy;
+	    hy<<=n;
 	}
-    }
-  hz = hx - hy; lz = lx - ly; if (lx < ly)
-    hz -= 1;
-  if (hz >= 0)
-    {
-      hx = hz; lx = lz;
-    }
 
-  /* convert back to floating value and restore the sign */
-  if ((hx | lx) == 0)                   /* return sign(x)*0 */
-    return Zero[(uint32_t) sx >> 31];
-  while (hx < 0x00100000)               /* normalize x */
-    {
-      hx = hx + hx + (lx >> 31); lx = lx + lx;
-      iy -= 1;
-    }
-  if (__glibc_likely (iy >= -1022))              /* normalize output */
-    {
-      hx = ((hx - 0x00100000) | ((iy + 1023) << 20));
-      INSERT_WORDS (x, hx | sx, lx);
-    }
-  else                          /* subnormal output */
-    {
-      n = -1022 - iy;
-      if (n <= 20)
-	{
-	  lx = (lx >> n) | ((uint32_t) hx << (32 - n));
-	  hx >>= n;
+    /* fix point fmod */
+	n = ix - iy;
+	while(n--) {
+	    hz=hx-hy;
+	    if(hz<0){hx = hx+hx;}
+	    else {
+		if(hz==0)		/* return sign(x)*0 */
+		    return Zero[(uint64_t)sx>>63];
+		hx = hz+hz;
+	    }
 	}
-      else if (n <= 31)
-	{
-	  lx = (hx << (32 - n)) | (lx >> n); hx = sx;
+	hz=hx-hy;
+	if(hz>=0) {hx=hz;}
+
+    /* convert back to floating value and restore the sign */
+	if(hx==0)			/* return sign(x)*0 */
+	    return Zero[(uint64_t)sx>>63];
+	while(hx<UINT64_C(0x0010000000000000)) {	/* normalize x */
+	    hx = hx+hx;
+	    iy -= 1;
 	}
-      else
-	{
-	  lx = hx >> (n - 32); hx = sx;
+	if(__builtin_expect(iy>= -1022, 1)) {	/* normalize output */
+	  hx = ((hx-UINT64_C(0x0010000000000000))|((uint64_t)(iy+1023)<<52));
+	    INSERT_WORDS64(x,hx|sx);
+	} else {		/* subnormal output */
+	    n = -1022 - iy;
+	    hx>>=n;
+	    INSERT_WORDS64(x,hx|sx);
+	    x *= one;		/* create necessary signal */
 	}
-      INSERT_WORDS (x, hx | sx, lx);
-      x *= one;                 /* create necessary signal */
-    }
-  return x;                     /* exact output */
+	return x;		/* exact output */
 }
 libm_alias_finite (__ieee754_fmod, __fmod)
diff --git a/sysdeps/ieee754/dbl-64/e_log10.c b/sysdeps/ieee754/dbl-64/e_log10.c
index 44a4bd2faa..b89064fb7c 100644
--- a/sysdeps/ieee754/dbl-64/e_log10.c
+++ b/sysdeps/ieee754/dbl-64/e_log10.c
@@ -44,44 +44,46 @@
  */
 
 #include <math.h>
-#include <math_private.h>
 #include <fix-int-fp-convert-zero.h>
+#include <math_private.h>
+#include <stdint.h>
 #include <libm-alias-finite.h>
 
-static const double two54 = 1.80143985094819840000e+16;         /* 0x43500000, 0x00000000 */
-static const double ivln10 = 4.34294481903251816668e-01;        /* 0x3FDBCB7B, 0x1526E50E */
-static const double log10_2hi = 3.01029995663611771306e-01;     /* 0x3FD34413, 0x509F6000 */
-static const double log10_2lo = 3.69423907715893078616e-13;     /* 0x3D59FEF3, 0x11F12B36 */
+static const double two54 = 1.80143985094819840000e+16;		/* 0x4350000000000000 */
+static const double ivln10 = 4.34294481903251816668e-01;	/* 0x3FDBCB7B1526E50E */
+static const double log10_2hi = 3.01029995663611771306e-01;	/* 0x3FD34413509F6000 */
+static const double log10_2lo = 3.69423907715893078616e-13;	/* 0x3D59FEF311F12B36 */
 
 double
 __ieee754_log10 (double x)
 {
   double y, z;
-  int32_t i, k, hx;
-  uint32_t lx;
+  int64_t i, hx;
+  int32_t k;
 
-  EXTRACT_WORDS (hx, lx, x);
+  EXTRACT_WORDS64 (hx, x);
 
   k = 0;
-  if (hx < 0x00100000)
-    {                           /* x < 2**-1022  */
-      if (__glibc_unlikely (((hx & 0x7fffffff) | lx) == 0))
-	return -two54 / fabs (x);	/* log(+-0)=-inf  */
+  if (hx < INT64_C(0x0010000000000000))
+    {				/* x < 2**-1022  */
+      if (__glibc_unlikely ((hx & UINT64_C(0x7fffffffffffffff)) == 0))
+	return -two54 / fabs (x);	/* log(+-0)=-inf */
       if (__glibc_unlikely (hx < 0))
-	return (x - x) / (x - x);       /* log(-#) = NaN */
+	return (x - x) / (x - x);	/* log(-#) = NaN */
       k -= 54;
-      x *= two54;               /* subnormal number, scale up x */
-      GET_HIGH_WORD (hx, x);
+      x *= two54;		/* subnormal number, scale up x */
+      EXTRACT_WORDS64 (hx, x);
     }
-  if (__glibc_unlikely (hx >= 0x7ff00000))
+  /* scale up resulted in a NaN number  */
+  if (__glibc_unlikely (hx >= UINT64_C(0x7ff0000000000000)))
     return x + x;
-  k += (hx >> 20) - 1023;
-  i = ((uint32_t) k & 0x80000000) >> 31;
-  hx = (hx & 0x000fffff) | ((0x3ff - i) << 20);
+  k += (hx >> 52) - 1023;
+  i = ((uint64_t) k & UINT64_C(0x8000000000000000)) >> 63;
+  hx = (hx & UINT64_C(0x000fffffffffffff)) | ((0x3ff - i) << 52);
   y = (double) (k + i);
   if (FIX_INT_FP_CONVERT_ZERO && y == 0.0)
     y = 0.0;
-  SET_HIGH_WORD (x, hx);
+  INSERT_WORDS64 (x, hx);
   z = y * log10_2lo + ivln10 * __ieee754_log (x);
   return z + y * log10_2hi;
 }
diff --git a/sysdeps/ieee754/dbl-64/s_frexp.c b/sysdeps/ieee754/dbl-64/s_frexp.c
index c96a869665..f6ddf4aaee 100644
--- a/sysdeps/ieee754/dbl-64/s_frexp.c
+++ b/sysdeps/ieee754/dbl-64/s_frexp.c
@@ -1,21 +1,28 @@
-/* @(#)s_frexp.c 5.1 93/09/24 */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
+/* Copyright (C) 2011-2021 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   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.
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: s_frexp.c,v 1.9 1995/05/10 20:47:24 jtc Exp $";
-#endif
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <https://www.gnu.org/licenses/>.  */
+
+#include <inttypes.h>
+#include <math.h>
+#include <math_private.h>
+#include <libm-alias-double.h>
 
 /*
- * for non-zero x
+ * for non-zero, finite x
  *	x = frexp(arg,&exp);
  * return a double fp quantity x such that 0.5 <= |x| <1.0
  * and the corresponding binary exponent "exp". That is
@@ -24,32 +31,36 @@ static char rcsid[] = "$NetBSD: s_frexp.c,v 1.9 1995/05/10 20:47:24 jtc Exp $";
  * with *exp=0.
  */
 
-#include <math.h>
-#include <math_private.h>
-#include <libm-alias-double.h>
-
-static const double
-  two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */
 
 double
 __frexp (double x, int *eptr)
 {
-  int32_t hx, ix, lx;
-  EXTRACT_WORDS (hx, lx, x);
-  ix = 0x7fffffff & hx;
-  *eptr = 0;
-  if (ix >= 0x7ff00000 || ((ix | lx) == 0))
-    return x + x;                                           /* 0,inf,nan */
-  if (ix < 0x00100000)                  /* subnormal */
+  int64_t ix;
+  EXTRACT_WORDS64 (ix, x);
+  int32_t ex = 0x7ff & (ix >> 52);
+  int e = 0;
+
+  if (__glibc_likely (ex != 0x7ff && x != 0.0))
     {
-      x *= two54;
-      GET_HIGH_WORD (hx, x);
-      ix = hx & 0x7fffffff;
-      *eptr = -54;
+      /* Not zero and finite.  */
+      e = ex - 1022;
+      if (__glibc_unlikely (ex == 0))
+	{
+	  /* Subnormal.  */
+	  x *= 0x1p54;
+	  EXTRACT_WORDS64 (ix, x);
+	  ex = 0x7ff & (ix >> 52);
+	  e = ex - 1022 - 54;
+	}
+
+      ix = (ix & INT64_C (0x800fffffffffffff)) | INT64_C (0x3fe0000000000000);
+      INSERT_WORDS64 (x, ix);
     }
-  *eptr += (ix >> 20) - 1022;
-  hx = (hx & 0x800fffff) | 0x3fe00000;
-  SET_HIGH_WORD (x, hx);
+  else
+    /* Quiet signaling NaNs.  */
+    x += x;
+
+  *eptr = e;
   return x;
 }
 libm_alias_double (__frexp, frexp)
diff --git a/sysdeps/ieee754/dbl-64/s_getpayload.c b/sysdeps/ieee754/dbl-64/s_getpayload.c
index 30eae8b2e0..d095d2077d 100644
--- a/sysdeps/ieee754/dbl-64/s_getpayload.c
+++ b/sysdeps/ieee754/dbl-64/s_getpayload.c
@@ -1,4 +1,4 @@
-/* Get NaN payload.  dbl-64 version.
+/* Get NaN payload.
    Copyright (C) 2016-2021 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
@@ -16,22 +16,21 @@
    License along with the GNU C Library; if not, see
    <https://www.gnu.org/licenses/>.  */
 
-#include <fix-int-fp-convert-zero.h>
 #include <math.h>
 #include <math_private.h>
 #include <libm-alias-double.h>
 #include <stdint.h>
+#include <fix-int-fp-convert-zero.h>
 
 double
 __getpayload (const double *x)
 {
-  uint32_t hx, lx;
-  EXTRACT_WORDS (hx, lx, *x);
-  if ((hx & 0x7ff00000) != 0x7ff00000
-      || ((hx & 0xfffff) | lx) == 0)
+  uint64_t ix;
+  EXTRACT_WORDS64 (ix, *x);
+  if ((ix & 0x7ff0000000000000ULL) != 0x7ff0000000000000ULL
+      || (ix & 0xfffffffffffffULL) == 0)
     return -1;
-  hx &= 0x7ffff;
-  uint64_t ix = ((uint64_t) hx << 32) | lx;
+  ix &= 0x7ffffffffffffULL;
   if (FIX_INT_FP_CONVERT_ZERO && ix == 0)
     return 0.0f;
   return (double) ix;
diff --git a/sysdeps/ieee754/dbl-64/s_issignaling.c b/sysdeps/ieee754/dbl-64/s_issignaling.c
index d3344e5116..5fb6fbc7d4 100644
--- a/sysdeps/ieee754/dbl-64/s_issignaling.c
+++ b/sysdeps/ieee754/dbl-64/s_issignaling.c
@@ -23,25 +23,21 @@
 int
 __issignaling (double x)
 {
+  uint64_t xi;
+  EXTRACT_WORDS64 (xi, x);
 #if HIGH_ORDER_BIT_IS_SET_FOR_SNAN
-  uint32_t hxi;
-  GET_HIGH_WORD (hxi, x);
   /* We only have to care about the high-order bit of x's significand, because
      having it set (sNaN) already makes the significand different from that
      used to designate infinity.  */
-  return (hxi & 0x7ff80000) == 0x7ff80000;
+  return (xi & UINT64_C (0x7ff8000000000000)) == UINT64_C (0x7ff8000000000000);
 #else
-  uint32_t hxi, lxi;
-  EXTRACT_WORDS (hxi, lxi, x);
   /* To keep the following comparison simple, toggle the quiet/signaling bit,
      so that it is set for sNaNs.  This is inverse to IEEE 754-2008 (as well as
      common practice for IEEE 754-1985).  */
-  hxi ^= 0x00080000;
-  /* If lxi != 0, then set any suitable bit of the significand in hxi.  */
-  hxi |= (lxi | -lxi) >> 31;
+  xi ^= UINT64_C (0x0008000000000000);
   /* We have to compare for greater (instead of greater or equal), because x's
      significand being all-zero designates infinity not NaN.  */
-  return (hxi & 0x7fffffff) > 0x7ff80000;
+  return (xi & UINT64_C (0x7fffffffffffffff)) > UINT64_C (0x7ff8000000000000);
 #endif
 }
 libm_hidden_def (__issignaling)
diff --git a/sysdeps/ieee754/dbl-64/s_llround.c b/sysdeps/ieee754/dbl-64/s_llround.c
index 69a55862b6..7020fd0156 100644
--- a/sysdeps/ieee754/dbl-64/s_llround.c
+++ b/sysdeps/ieee754/dbl-64/s_llround.c
@@ -17,54 +17,43 @@
    License along with the GNU C Library; if not, see
    <https://www.gnu.org/licenses/>.  */
 
+#define lround __hidden_lround
+#define __lround __hidden___lround
+
 #include <fenv.h>
 #include <limits.h>
 #include <math.h>
+#include <sysdep.h>
 
 #include <math_private.h>
 #include <libm-alias-double.h>
 #include <fix-fp-int-convert-overflow.h>
 
-
 long long int
 __llround (double x)
 {
   int32_t j0;
-  uint32_t i1, i0;
+  int64_t i0;
   long long int result;
   int sign;
 
-  EXTRACT_WORDS (i0, i1, x);
-  j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
-  sign = (i0 & 0x80000000) != 0 ? -1 : 1;
-  i0 &= 0xfffff;
-  i0 |= 0x100000;
+  EXTRACT_WORDS64 (i0, x);
+  j0 = ((i0 >> 52) & 0x7ff) - 0x3ff;
+  sign = i0 < 0 ? -1 : 1;
+  i0 &= UINT64_C(0xfffffffffffff);
+  i0 |= UINT64_C(0x10000000000000);
 
-  if (j0 < 20)
+  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 += 0x80000 >> j0;
-
-	  result = i0 >> (20 - j0);
-	}
-    }
-  else if (j0 < (int32_t) (8 * sizeof (long long int)) - 1)
-    {
-      if (j0 >= 52)
-	result = (((long long int) i0 << 32) | i1) << (j0 - 52);
-      else
-	{
-	  uint32_t j = i1 + (0x80000000 >> (j0 - 20));
-	  if (j < i1)
-	    ++i0;
+	  i0 += UINT64_C(0x8000000000000) >> j0;
 
-	  if (j0 == 20)
-	    result = (long long int) i0;
-	  else
-	    result = ((long long int) i0 << (j0 - 20)) | (j >> (52 - j0));
+	  result = i0 >> (52 - j0);
 	}
     }
   else
@@ -86,3 +75,11 @@ __llround (double x)
 }
 
 libm_alias_double (__llround, llround)
+
+/* long has the same width as long long on LP64 machines, so use an alias.  */
+#undef lround
+#undef __lround
+#ifdef _LP64
+strong_alias (__llround, __lround)
+libm_alias_double (__lround, lround)
+#endif
diff --git a/sysdeps/ieee754/dbl-64/s_lround.c b/sysdeps/ieee754/dbl-64/s_lround.c
index c7d097ee5d..5284c4da09 100644
--- a/sysdeps/ieee754/dbl-64/s_lround.c
+++ b/sysdeps/ieee754/dbl-64/s_lround.c
@@ -1,7 +1,6 @@
 /* Round double value to long int.
    Copyright (C) 1997-2021 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
@@ -25,55 +24,41 @@
 #include <libm-alias-double.h>
 #include <fix-fp-int-convert-overflow.h>
 
+/* For LP64, lround is an alias for llround.  */
+#ifndef _LP64
 
 long int
 __lround (double x)
 {
   int32_t j0;
-  uint32_t i1, i0;
+  int64_t i0;
   long int result;
   int sign;
 
-  EXTRACT_WORDS (i0, i1, x);
-  j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
-  sign = (i0 & 0x80000000) != 0 ? -1 : 1;
-  i0 &= 0xfffff;
-  i0 |= 0x100000;
+  EXTRACT_WORDS64 (i0, x);
+  j0 = ((i0 >> 52) & 0x7ff) - 0x3ff;
+  sign = i0 < 0 ? -1 : 1;
+  i0 &= UINT64_C(0xfffffffffffff);
+  i0 |= UINT64_C(0x10000000000000);
 
-  if (j0 < 20)
+  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 += 0x80000 >> j0;
+	  i0 += UINT64_C(0x8000000000000) >> j0;
 
-	  result = i0 >> (20 - j0);
-	}
-    }
-  else if (j0 < (int32_t) (8 * sizeof (long int)) - 1)
-    {
-      if (j0 >= 52)
-	result = ((long int) i0 << (j0 - 20)) | ((long int) i1 << (j0 - 52));
-      else
-	{
-	  uint32_t j = i1 + (0x80000000 >> (j0 - 20));
-	  if (j < i1)
-	    ++i0;
-
-	  if (j0 == 20)
-	    result = (long int) i0;
-	  else
-	    {
-	      result = ((long int) i0 << (j0 - 20)) | (j >> (52 - j0));
+	  result = i0 >> (52 - j0);
 #ifdef FE_INVALID
-	      if (sizeof (long int) == 4
-		  && sign == 1
-		  && result == LONG_MIN)
-		/* Rounding brought the value out of range.  */
-		feraiseexcept (FE_INVALID);
+	  if (sizeof (long int) == 4
+	      && sign == 1
+	      && result == LONG_MIN)
+	    /* Rounding brought the value out of range.  */
+	    feraiseexcept (FE_INVALID);
 #endif
-	    }
 	}
     }
   else
@@ -92,8 +77,8 @@ __lround (double x)
 	  return sign == 1 ? LONG_MAX : LONG_MIN;
 	}
       else if (!FIX_DBL_LONG_CONVERT_OVERFLOW
-	       && sizeof (long int) == 4
-	       && x <= (double) LONG_MIN - 0.5)
+	  && sizeof (long int) == 4
+	  && x <= (double) LONG_MIN - 0.5)
 	{
 	  /* If truncation produces LONG_MIN, the cast will not raise
 	     the exception, but may raise "inexact".  */
@@ -108,3 +93,5 @@ __lround (double x)
 }
 
 libm_alias_double (__lround, lround)
+
+#endif
diff --git a/sysdeps/ieee754/dbl-64/s_modf.c b/sysdeps/ieee754/dbl-64/s_modf.c
index 722511c64a..8d14e78ef0 100644
--- a/sysdeps/ieee754/dbl-64/s_modf.c
+++ b/sysdeps/ieee754/dbl-64/s_modf.c
@@ -1,3 +1,4 @@
+/* Rewritten for 64-bit machines by Ulrich Drepper <drepper@gmail.com>.  */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -22,63 +23,42 @@
 #include <math.h>
 #include <math_private.h>
 #include <libm-alias-double.h>
+#include <stdint.h>
 
 static const double one = 1.0;
 
 double
-__modf (double x, double *iptr)
+__modf(double x, double *iptr)
 {
-  int32_t i0, i1, j0;
-  uint32_t i;
-  EXTRACT_WORDS (i0, i1, x);
-  j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;    /* exponent of x */
-  if (j0 < 20)                          /* integer part in high x */
-    {
-      if (j0 < 0)                       /* |x|<1 */
-	{
-	  INSERT_WORDS (*iptr, i0 & 0x80000000, 0);     /* *iptr = +-0 */
-	  return x;
-	}
-      else
-	{
-	  i = (0x000fffff) >> j0;
-	  if (((i0 & i) | i1) == 0)             /* x is integral */
-	    {
-	      *iptr = x;
-	      INSERT_WORDS (x, i0 & 0x80000000, 0);     /* return +-0 */
-	      return x;
-	    }
-	  else
-	    {
-	      INSERT_WORDS (*iptr, i0 & (~i), 0);
-	      return x - *iptr;
+	int64_t i0;
+	int32_t j0;
+	EXTRACT_WORDS64(i0,x);
+	j0 = ((i0>>52)&0x7ff)-0x3ff;	/* exponent of x */
+	if(j0<52) {			/* integer part in x */
+	    if(j0<0) {			/* |x|<1 */
+		/* *iptr = +-0 */
+		INSERT_WORDS64(*iptr,i0&UINT64_C(0x8000000000000000));
+		return x;
+	    } else {
+		uint64_t i = UINT64_C(0x000fffffffffffff)>>j0;
+		if((i0&i)==0) {		/* x is integral */
+		    *iptr = x;
+		    /* return +-0 */
+		    INSERT_WORDS64(x,i0&UINT64_C(0x8000000000000000));
+		    return x;
+		} else {
+		    INSERT_WORDS64(*iptr,i0&(~i));
+		    return x - *iptr;
+		}
 	    }
+	} else { /* no fraction part */
+	    *iptr = x*one;
+	    /* We must handle NaNs separately.  */
+	    if (j0 == 0x400 && (i0 & UINT64_C(0xfffffffffffff)))
+	      return x*one;
+	    INSERT_WORDS64(x,i0&UINT64_C(0x8000000000000000));	/* return +-0 */
+	    return x;
 	}
-    }
-  else if (__glibc_unlikely (j0 > 51))              /* no fraction part */
-    {
-      *iptr = x * one;
-      /* We must handle NaNs separately.  */
-      if (j0 == 0x400 && ((i0 & 0xfffff) | i1))
-	return x * one;
-      INSERT_WORDS (x, i0 & 0x80000000, 0);     /* return +-0 */
-      return x;
-    }
-  else                                  /* fraction part in low x */
-    {
-      i = ((uint32_t) (0xffffffff)) >> (j0 - 20);
-      if ((i1 & i) == 0)                /* x is integral */
-	{
-	  *iptr = x;
-	  INSERT_WORDS (x, i0 & 0x80000000, 0);         /* return +-0 */
-	  return x;
-	}
-      else
-	{
-	  INSERT_WORDS (*iptr, i0, i1 & (~i));
-	  return x - *iptr;
-	}
-    }
 }
 #ifndef __modf
 libm_alias_double (__modf, modf)
diff --git a/sysdeps/ieee754/dbl-64/s_remquo.c b/sysdeps/ieee754/dbl-64/s_remquo.c
index 928c379e39..cbaa7f79a2 100644
--- a/sysdeps/ieee754/dbl-64/s_remquo.c
+++ b/sysdeps/ieee754/dbl-64/s_remquo.c
@@ -21,7 +21,7 @@
 
 #include <math_private.h>
 #include <libm-alias-double.h>
-
+#include <stdint.h>
 
 static const double zero = 0.0;
 
@@ -29,50 +29,49 @@ static const double zero = 0.0;
 double
 __remquo (double x, double y, int *quo)
 {
-  int32_t hx, hy;
-  uint32_t sx, lx, ly;
-  int cquo, qs;
+  int64_t hx, hy;
+  uint64_t sx, qs;
+  int cquo;
 
-  EXTRACT_WORDS (hx, lx, x);
-  EXTRACT_WORDS (hy, ly, y);
-  sx = hx & 0x80000000;
-  qs = sx ^ (hy & 0x80000000);
-  hy &= 0x7fffffff;
-  hx &= 0x7fffffff;
+  EXTRACT_WORDS64 (hx, x);
+  EXTRACT_WORDS64 (hy, y);
+  sx = hx & UINT64_C(0x8000000000000000);
+  qs = sx ^ (hy & UINT64_C(0x8000000000000000));
+  hy &= UINT64_C(0x7fffffffffffffff);
+  hx &= UINT64_C(0x7fffffffffffffff);
 
   /* Purge off exception values.  */
-  if ((hy | ly) == 0)
-    return (x * y) / (x * y);                   /* y = 0 */
-  if ((hx >= 0x7ff00000)                        /* x not finite */
-      || ((hy >= 0x7ff00000)                    /* p is NaN */
-	  && (((hy - 0x7ff00000) | ly) != 0)))
+  if (__glibc_unlikely (hy == 0))
+    return (x * y) / (x * y);			/* y = 0 */
+  if (__builtin_expect (hx >= UINT64_C(0x7ff0000000000000) /* x not finite */
+			|| hy > UINT64_C(0x7ff0000000000000), 0))/* y is NaN */
     return (x * y) / (x * y);
 
-  if (hy <= 0x7fbfffff)
-    x = __ieee754_fmod (x, 8 * y);              /* now x < 8y */
+  if (hy <= UINT64_C(0x7fbfffffffffffff))
+    x = __ieee754_fmod (x, 8 * y);		/* now x < 8y */
 
-  if (((hx - hy) | (lx - ly)) == 0)
+  if (__glibc_unlikely (hx == hy))
     {
       *quo = qs ? -1 : 1;
       return zero * x;
     }
 
   x = fabs (x);
-  y = fabs (y);
+  INSERT_WORDS64 (y, hy);
   cquo = 0;
 
-  if (hy <= 0x7fcfffff && x >= 4 * y)
+  if (hy <= UINT64_C(0x7fcfffffffffffff) && x >= 4 * y)
     {
       x -= 4 * y;
       cquo += 4;
     }
-  if (hy <= 0x7fdfffff && x >= 2 * y)
+  if (hy <= UINT64_C(0x7fdfffffffffffff) && x >= 2 * y)
     {
       x -= 2 * y;
       cquo += 2;
     }
 
-  if (hy < 0x00200000)
+  if (hy < UINT64_C(0x0020000000000000))
     {
       if (x + x > y)
 	{
diff --git a/sysdeps/ieee754/dbl-64/s_roundeven.c b/sysdeps/ieee754/dbl-64/s_roundeven.c
index b7f4bd63ee..943b2c634c 100644
--- a/sysdeps/ieee754/dbl-64/s_roundeven.c
+++ b/sysdeps/ieee754/dbl-64/s_roundeven.c
@@ -1,5 +1,4 @@
 /* Round to nearest integer value, rounding halfway cases to even.
-   dbl-64 version.
    Copyright (C) 2016-2021 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
@@ -29,10 +28,10 @@
 double
 __roundeven (double x)
 {
-  uint32_t hx, lx, uhx;
-  EXTRACT_WORDS (hx, lx, x);
-  uhx = hx & 0x7fffffff;
-  int exponent = uhx >> (MANT_DIG - 1 - 32);
+  uint64_t ix, ux;
+  EXTRACT_WORDS64 (ix, x);
+  ux = ix & 0x7fffffffffffffffULL;
+  int exponent = ux >> (MANT_DIG - 1);
   if (exponent >= BIAS + MANT_DIG - 1)
     {
       /* Integer, infinity or NaN.  */
@@ -42,63 +41,29 @@ __roundeven (double x)
       else
 	return x;
     }
-  else if (exponent >= BIAS + MANT_DIG - 32)
-    {
-      /* Not necessarily an integer; integer bit is in low word.
-	 Locate the bits with exponents 0 and -1.  */
-      int int_pos = (BIAS + MANT_DIG - 1) - exponent;
-      int half_pos = int_pos - 1;
-      uint32_t half_bit = 1U << half_pos;
-      uint32_t int_bit = 1U << int_pos;
-      if ((lx & (int_bit | (half_bit - 1))) != 0)
-	{
-	  /* Carry into the exponent works correctly.  No need to test
-	     whether HALF_BIT is set.  */
-	  lx += half_bit;
-	  hx += lx < half_bit;
-	}
-      lx &= ~(int_bit - 1);
-    }
-  else if (exponent == BIAS + MANT_DIG - 33)
-    {
-      /* Not necessarily an integer; integer bit is bottom of high
-	 word, half bit is top of low word.  */
-      if (((hx & 1) | (lx & 0x7fffffff)) != 0)
-	{
-	  lx += 0x80000000;
-	  hx += lx < 0x80000000;
-	}
-      lx = 0;
-    }
   else if (exponent >= BIAS)
     {
-      /* At least 1; not necessarily an integer, integer bit and half
-	 bit are in the high word.  Locate the bits with exponents 0
-	 and -1 (when the unbiased exponent is 0, the bit with
-	 exponent 0 is implicit, but as the bias is odd it is OK to
-	 take it from the low bit of the exponent).  */
-      int int_pos = (BIAS + MANT_DIG - 33) - exponent;
+      /* At least 1; not necessarily an integer.  Locate the bits with
+	 exponents 0 and -1 (when the unbiased exponent is 0, the bit
+	 with exponent 0 is implicit, but as the bias is odd it is OK
+	 to take it from the low bit of the exponent).  */
+      int int_pos = (BIAS + MANT_DIG - 1) - exponent;
       int half_pos = int_pos - 1;
-      uint32_t half_bit = 1U << half_pos;
-      uint32_t int_bit = 1U << int_pos;
-      if (((hx & (int_bit | (half_bit - 1))) | lx) != 0)
-	hx += half_bit;
-      hx &= ~(int_bit - 1);
-      lx = 0;
-    }
-  else if (exponent == BIAS - 1 && (uhx > 0x3fe00000 || lx != 0))
-    {
-      /* Interval (0.5, 1).  */
-      hx = (hx & 0x80000000) | 0x3ff00000;
-      lx = 0;
+      uint64_t half_bit = 1ULL << half_pos;
+      uint64_t int_bit = 1ULL << int_pos;
+      if ((ix & (int_bit | (half_bit - 1))) != 0)
+	/* Carry into the exponent works correctly.  No need to test
+	   whether HALF_BIT is set.  */
+	ix += half_bit;
+      ix &= ~(int_bit - 1);
     }
+  else if (exponent == BIAS - 1 && ux > 0x3fe0000000000000ULL)
+    /* Interval (0.5, 1).  */
+    ix = (ix & 0x8000000000000000ULL) | 0x3ff0000000000000ULL;
   else
-    {
-      /* Rounds to 0.  */
-      hx &= 0x80000000;
-      lx = 0;
-    }
-  INSERT_WORDS (x, hx, lx);
+    /* Rounds to 0.  */
+    ix &= 0x8000000000000000ULL;
+  INSERT_WORDS64 (x, ix);
   return x;
 }
 hidden_def (__roundeven)
diff --git a/sysdeps/ieee754/dbl-64/s_scalbln.c b/sysdeps/ieee754/dbl-64/s_scalbln.c
index 0e3d732e48..071c9d7794 100644
--- a/sysdeps/ieee754/dbl-64/s_scalbln.c
+++ b/sysdeps/ieee754/dbl-64/s_scalbln.c
@@ -20,43 +20,40 @@
 #include <math_private.h>
 
 static const double
-  two54 = 1.80143985094819840000e+16,  /* 0x43500000, 0x00000000 */
-  twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
-  huge = 1.0e+300,
-  tiny = 1.0e-300;
+two54   =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
+twom54  =  5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
+huge   = 1.0e+300,
+tiny   = 1.0e-300;
 
 double
 __scalbln (double x, long int n)
 {
-  int32_t k, hx, lx;
-  EXTRACT_WORDS (hx, lx, x);
-  k = (hx & 0x7ff00000) >> 20;                  /* extract exponent */
-  if (__glibc_unlikely (k == 0))                   /* 0 or subnormal x */
-    {
-      if ((lx | (hx & 0x7fffffff)) == 0)
-	return x;                                  /* +-0 */
-      x *= two54;
-      GET_HIGH_WORD (hx, x);
-      k = ((hx & 0x7ff00000) >> 20) - 54;
-    }
-  if (__glibc_unlikely (k == 0x7ff))
-    return x + x;                                       /* NaN or Inf */
-  if (__glibc_unlikely (n < -50000))
-    return tiny * copysign (tiny, x);   /*underflow*/
-  if (__glibc_unlikely (n > 50000 || k + n > 0x7fe))
-    return huge * copysign (huge, x);   /* overflow  */
-  /* Now k and n are bounded we know that k = k+n does not
-     overflow.  */
-  k = k + n;
-  if (__glibc_likely (k > 0))                    /* normal result */
-    {
-      SET_HIGH_WORD (x, (hx & 0x800fffff) | (k << 20)); return x;
-    }
-  if (k <= -54)
-    return tiny * copysign (tiny, x);         /*underflow*/
-  k += 54;                                      /* subnormal result */
-  SET_HIGH_WORD (x, (hx & 0x800fffff) | (k << 20));
-  return x * twom54;
+	int64_t ix;
+	int64_t k;
+	EXTRACT_WORDS64(ix,x);
+	k = (ix >> 52) & 0x7ff;			/* extract exponent */
+	if (__builtin_expect(k==0, 0)) {	/* 0 or subnormal x */
+	    if ((ix & UINT64_C(0xfffffffffffff))==0) return x; /* +-0 */
+	    x *= two54;
+	    EXTRACT_WORDS64(ix,x);
+	    k = ((ix >> 52) & 0x7ff) - 54;
+	    }
+	if (__builtin_expect(k==0x7ff, 0)) return x+x;	/* NaN or Inf */
+	if (__builtin_expect(n< -50000, 0))
+	  return tiny*copysign(tiny,x); /*underflow*/
+	if (__builtin_expect(n> 50000 || k+n > 0x7fe, 0))
+	  return huge*copysign(huge,x); /* overflow  */
+	/* Now k and n are bounded we know that k = k+n does not
+	   overflow.  */
+	k = k+n;
+	if (__builtin_expect(k > 0, 1))		/* normal result */
+	    {INSERT_WORDS64(x,(ix&UINT64_C(0x800fffffffffffff))|(k<<52));
+	      return x;}
+	if (k <= -54)
+	  return tiny*copysign(tiny,x);	/*underflow*/
+	k += 54;				/* subnormal result */
+	INSERT_WORDS64(x,(ix&INT64_C(0x800fffffffffffff))|(k<<52));
+	return x*twom54;
 }
 #ifdef NO_LONG_DOUBLE
 strong_alias (__scalbln, __scalblnl)
diff --git a/sysdeps/ieee754/dbl-64/s_scalbn.c b/sysdeps/ieee754/dbl-64/s_scalbn.c
index cf4d6846ee..4491227f3e 100644
--- a/sysdeps/ieee754/dbl-64/s_scalbn.c
+++ b/sysdeps/ieee754/dbl-64/s_scalbn.c
@@ -20,43 +20,40 @@
 #include <math_private.h>
 
 static const double
-  two54 = 1.80143985094819840000e+16,  /* 0x43500000, 0x00000000 */
-  twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
-  huge = 1.0e+300,
-  tiny = 1.0e-300;
+two54   =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
+twom54  =  5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
+huge   = 1.0e+300,
+tiny   = 1.0e-300;
 
 double
 __scalbn (double x, int n)
 {
-  int32_t k, hx, lx;
-  EXTRACT_WORDS (hx, lx, x);
-  k = (hx & 0x7ff00000) >> 20;                  /* extract exponent */
-  if (__glibc_unlikely (k == 0))                   /* 0 or subnormal x */
-    {
-      if ((lx | (hx & 0x7fffffff)) == 0)
-	return x;                                  /* +-0 */
-      x *= two54;
-      GET_HIGH_WORD (hx, x);
-      k = ((hx & 0x7ff00000) >> 20) - 54;
-    }
-  if (__glibc_unlikely (k == 0x7ff))
-    return x + x;                                       /* NaN or Inf */
-  if (__glibc_unlikely (n < -50000))
-    return tiny * copysign (tiny, x);   /*underflow*/
-  if (__glibc_unlikely (n > 50000 || k + n > 0x7fe))
-    return huge * copysign (huge, x);   /* overflow  */
-  /* Now k and n are bounded we know that k = k+n does not
-     overflow.  */
-  k = k + n;
-  if (__glibc_likely (k > 0))                    /* normal result */
-    {
-      SET_HIGH_WORD (x, (hx & 0x800fffff) | (k << 20)); return x;
-    }
-  if (k <= -54)
-    return tiny * copysign (tiny, x);         /*underflow*/
-  k += 54;                                      /* subnormal result */
-  SET_HIGH_WORD (x, (hx & 0x800fffff) | (k << 20));
-  return x * twom54;
+	int64_t ix;
+	int64_t k;
+	EXTRACT_WORDS64(ix,x);
+	k = (ix >> 52) & 0x7ff;			/* extract exponent */
+	if (__builtin_expect(k==0, 0)) {	/* 0 or subnormal x */
+	    if ((ix & UINT64_C(0xfffffffffffff))==0) return x; /* +-0 */
+	    x *= two54;
+	    EXTRACT_WORDS64(ix,x);
+	    k = ((ix >> 52) & 0x7ff) - 54;
+	    }
+	if (__builtin_expect(k==0x7ff, 0)) return x+x;	/* NaN or Inf */
+	if (__builtin_expect(n< -50000, 0))
+	  return tiny*copysign(tiny,x); /*underflow*/
+	if (__builtin_expect(n> 50000 || k+n > 0x7fe, 0))
+	  return huge*copysign(huge,x); /* overflow  */
+	/* Now k and n are bounded we know that k = k+n does not
+	   overflow.  */
+	k = k+n;
+	if (__builtin_expect(k > 0, 1))		/* normal result */
+	    {INSERT_WORDS64(x,(ix&UINT64_C(0x800fffffffffffff))|(k<<52));
+	      return x;}
+	if (k <= -54)
+	  return tiny*copysign(tiny,x);	/*underflow*/
+	k += 54;				/* subnormal result */
+	INSERT_WORDS64(x,(ix&INT64_C(0x800fffffffffffff))|(k<<52));
+	return x*twom54;
 }
 #ifdef NO_LONG_DOUBLE
 strong_alias (__scalbn, __scalbnl)
diff --git a/sysdeps/ieee754/dbl-64/s_setpayload_main.c b/sysdeps/ieee754/dbl-64/s_setpayload_main.c
index 0b0a295d6c..e0014a3b09 100644
--- a/sysdeps/ieee754/dbl-64/s_setpayload_main.c
+++ b/sysdeps/ieee754/dbl-64/s_setpayload_main.c
@@ -1,4 +1,4 @@
-/* Set NaN payload.  dbl-64 version.
+/* Set NaN payload.
    Copyright (C) 2016-2021 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
@@ -30,41 +30,25 @@
 int
 FUNC (double *x, double payload)
 {
-  uint32_t hx, lx;
-  EXTRACT_WORDS (hx, lx, payload);
-  int exponent = hx >> (EXPLICIT_MANT_DIG - 32);
+  uint64_t ix;
+  EXTRACT_WORDS64 (ix, payload);
+  int exponent = ix >> EXPLICIT_MANT_DIG;
   /* Test if argument is (a) negative or too large; (b) too small,
      except for 0 when allowed; (c) not an integer.  */
   if (exponent >= BIAS + PAYLOAD_DIG
-      || (exponent < BIAS && !(SET_HIGH_BIT && hx == 0 && lx == 0)))
+      || (exponent < BIAS && !(SET_HIGH_BIT && ix == 0))
+      || (ix & ((1ULL << (BIAS + EXPLICIT_MANT_DIG - exponent)) - 1)) != 0)
     {
-      INSERT_WORDS (*x, 0, 0);
+      INSERT_WORDS64 (*x, 0);
       return 1;
     }
-  int shift = BIAS + EXPLICIT_MANT_DIG - exponent;
-  if (shift < 32
-      ? (lx & ((1U << shift) - 1)) != 0
-      : (lx != 0 || (hx & ((1U << (shift - 32)) - 1)) != 0))
+  if (ix != 0)
     {
-      INSERT_WORDS (*x, 0, 0);
-      return 1;
-    }
-  if (exponent != 0)
-    {
-      hx &= (1U << (EXPLICIT_MANT_DIG - 32)) - 1;
-      hx |= 1U << (EXPLICIT_MANT_DIG - 32);
-      if (shift >= 32)
-	{
-	  lx = hx >> (shift - 32);
-	  hx = 0;
-	}
-      else if (shift != 0)
-	{
-	  lx = (lx >> shift) | (hx << (32 - shift));
-	  hx >>= shift;
-	}
+      ix &= (1ULL << EXPLICIT_MANT_DIG) - 1;
+      ix |= 1ULL << EXPLICIT_MANT_DIG;
+      ix >>= BIAS + EXPLICIT_MANT_DIG - exponent;
     }
-  hx |= 0x7ff00000 | (SET_HIGH_BIT ? 0x80000 : 0);
-  INSERT_WORDS (*x, hx, lx);
+  ix |= 0x7ff0000000000000ULL | (SET_HIGH_BIT ? 0x8000000000000ULL : 0);
+  INSERT_WORDS64 (*x, ix);
   return 0;
 }
diff --git a/sysdeps/ieee754/dbl-64/s_totalorder.c b/sysdeps/ieee754/dbl-64/s_totalorder.c
index ace32e0827..13bde9e538 100644
--- a/sysdeps/ieee754/dbl-64/s_totalorder.c
+++ b/sysdeps/ieee754/dbl-64/s_totalorder.c
@@ -1,4 +1,4 @@
-/* Total order operation.  dbl-64 version.
+/* Total order operation.
    Copyright (C) 2016-2021 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
@@ -18,8 +18,8 @@
 
 #include <math.h>
 #include <math_private.h>
-#include <libm-alias-double.h>
 #include <nan-high-order-bit.h>
+#include <libm-alias-double.h>
 #include <stdint.h>
 #include <shlib-compat.h>
 #include <first-versions.h>
@@ -27,30 +27,26 @@
 int
 __totalorder (const double *x, const double *y)
 {
-  int32_t hx, hy;
-  uint32_t lx, ly;
-  EXTRACT_WORDS (hx, lx, *x);
-  EXTRACT_WORDS (hy, ly, *y);
+  int64_t ix, iy;
+  EXTRACT_WORDS64 (ix, *x);
+  EXTRACT_WORDS64 (iy, *y);
 #if HIGH_ORDER_BIT_IS_SET_FOR_SNAN
-  uint32_t uhx = hx & 0x7fffffff, uhy = hy & 0x7fffffff;
   /* For the preferred quiet NaN convention, this operation is a
      comparison of the representations of the arguments interpreted as
      sign-magnitude integers.  If both arguments are NaNs, invert the
      quiet/signaling bit so comparing that way works.  */
-  if ((uhx > 0x7ff00000 || (uhx == 0x7ff00000 && lx != 0))
-      && (uhy > 0x7ff00000 || (uhy == 0x7ff00000 && ly != 0)))
+  if ((ix & 0x7fffffffffffffffULL) > 0x7ff0000000000000ULL
+      && (iy & 0x7fffffffffffffffULL) > 0x7ff0000000000000ULL)
     {
-      hx ^= 0x00080000;
-      hy ^= 0x00080000;
+      ix ^= 0x0008000000000000ULL;
+      iy ^= 0x0008000000000000ULL;
     }
 #endif
-  uint32_t hx_sign = hx >> 31;
-  uint32_t hy_sign = hy >> 31;
-  hx ^= hx_sign >> 1;
-  lx ^= hx_sign;
-  hy ^= hy_sign >> 1;
-  ly ^= hy_sign;
-  return hx < hy || (hx == hy && lx <= ly);
+  uint64_t ix_sign = ix >> 63;
+  uint64_t iy_sign = iy >> 63;
+  ix ^= ix_sign >> 1;
+  iy ^= iy_sign >> 1;
+  return ix <= iy;
 }
 #ifdef SHARED
 # define CONCATX(x, y) x ## y
diff --git a/sysdeps/ieee754/dbl-64/s_totalordermag.c b/sysdeps/ieee754/dbl-64/s_totalordermag.c
index e6efc387a2..fd8aade28c 100644
--- a/sysdeps/ieee754/dbl-64/s_totalordermag.c
+++ b/sysdeps/ieee754/dbl-64/s_totalordermag.c
@@ -1,4 +1,4 @@
-/* Total order operation on absolute values.  dbl-64 version.
+/* Total order operation on absolute values.
    Copyright (C) 2016-2021 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
@@ -18,8 +18,8 @@
 
 #include <math.h>
 #include <math_private.h>
-#include <libm-alias-double.h>
 #include <nan-high-order-bit.h>
+#include <libm-alias-double.h>
 #include <stdint.h>
 #include <shlib-compat.h>
 #include <first-versions.h>
@@ -27,25 +27,23 @@
 int
 __totalordermag (const double *x, const double *y)
 {
-  uint32_t hx, hy;
-  uint32_t lx, ly;
-  EXTRACT_WORDS (hx, lx, *x);
-  EXTRACT_WORDS (hy, ly, *y);
-  hx &= 0x7fffffff;
-  hy &= 0x7fffffff;
+  uint64_t ix, iy;
+  EXTRACT_WORDS64 (ix, *x);
+  EXTRACT_WORDS64 (iy, *y);
+  ix &= 0x7fffffffffffffffULL;
+  iy &= 0x7fffffffffffffffULL;
 #if HIGH_ORDER_BIT_IS_SET_FOR_SNAN
   /* For the preferred quiet NaN convention, this operation is a
      comparison of the representations of the absolute values of the
      arguments.  If both arguments are NaNs, invert the
      quiet/signaling bit so comparing that way works.  */
-  if ((hx > 0x7ff00000 || (hx == 0x7ff00000 && lx != 0))
-      && (hy > 0x7ff00000 || (hy == 0x7ff00000 && ly != 0)))
+  if (ix > 0x7ff0000000000000ULL && iy > 0x7ff0000000000000ULL)
     {
-      hx ^= 0x00080000;
-      hy ^= 0x00080000;
+      ix ^= 0x0008000000000000ULL;
+      iy ^= 0x0008000000000000ULL;
     }
 #endif
-  return hx < hy || (hx == hy && lx <= ly);
+  return ix <= iy;
 }
 #ifdef SHARED
 # define CONCATX(x, y) x ## y
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/e_acosh.c b/sysdeps/ieee754/dbl-64/wordsize-64/e_acosh.c
deleted file mode 100644
index a241366f30..0000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/e_acosh.c
+++ /dev/null
@@ -1,68 +0,0 @@
-/* Optimized for 64-bit by Ulrich Drepper <drepper@gmail.com>, 2012 */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-/* __ieee754_acosh(x)
- * Method :
- *	Based on
- *		acosh(x) = log [ x + sqrt(x*x-1) ]
- *	we have
- *		acosh(x) := log(x)+ln2,	if x is large; else
- *		acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
- *		acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
- *
- * Special cases:
- *	acosh(x) is NaN with signal if x<1.
- *	acosh(NaN) is NaN without signal.
- */
-
-#include <math.h>
-#include <math_private.h>
-#include <libm-alias-finite.h>
-
-static const double
-one	= 1.0,
-ln2	= 6.93147180559945286227e-01;  /* 0x3FE62E42, 0xFEFA39EF */
-
-double
-__ieee754_acosh (double x)
-{
-  int64_t hx;
-  EXTRACT_WORDS64 (hx, x);
-
-  if (hx > INT64_C (0x4000000000000000))
-    {
-      if (__glibc_unlikely (hx >= INT64_C (0x41b0000000000000)))
-	{
-	  /* x > 2**28 */
-	  if (hx >= INT64_C (0x7ff0000000000000))
-	    /* x is inf of NaN */
-	    return x + x;
-	  else
-	    return __ieee754_log (x) + ln2;/* acosh(huge)=log(2x) */
-	}
-
-      /* 2**28 > x > 2 */
-      double t = x * x;
-      return __ieee754_log (2.0 * x - one / (x + sqrt (t - one)));
-    }
-  else if (__glibc_likely (hx > INT64_C (0x3ff0000000000000)))
-    {
-      /* 1<x<2 */
-      double t = x - one;
-      return __log1p (t + sqrt (2.0 * t + t * t));
-    }
-  else if (__glibc_likely (hx == INT64_C (0x3ff0000000000000)))
-    return 0.0;				/* acosh(1) = 0 */
-  else					/* x < 1 */
-    return (x - x) / (x - x);
-}
-libm_alias_finite (__ieee754_acosh, __acosh)
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/e_cosh.c b/sysdeps/ieee754/dbl-64/wordsize-64/e_cosh.c
deleted file mode 100644
index 4f41ca2c92..0000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/e_cosh.c
+++ /dev/null
@@ -1,85 +0,0 @@
-/* Optimized by Ulrich Drepper <drepper@gmail.com>, 2011 */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-/* __ieee754_cosh(x)
- * Method :
- * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
- *	1. Replace x by |x| (cosh(x) = cosh(-x)).
- *	2.
- *							[ exp(x) - 1 ]^2
- *	    0        <= x <= ln2/2  :  cosh(x) := 1 + -------------------
- *							   2*exp(x)
- *
- *						  exp(x) +  1/exp(x)
- *	    ln2/2    <= x <= 22     :  cosh(x) := -------------------
- *							  2
- *	    22       <= x <= lnovft :  cosh(x) := exp(x)/2
- *	    lnovft   <= x <= ln2ovft:  cosh(x) := exp(x/2)/2 * exp(x/2)
- *	    ln2ovft  <  x	    :  cosh(x) := huge*huge (overflow)
- *
- * Special cases:
- *	cosh(x) is |x| if x is +INF, -INF, or NaN.
- *	only cosh(0)=1 is exact for finite x.
- */
-
-#include <math.h>
-#include <math_private.h>
-#include <libm-alias-finite.h>
-
-static const double one = 1.0, half=0.5, huge = 1.0e300;
-
-double
-__ieee754_cosh (double x)
-{
-	double t,w;
-	int32_t ix;
-
-    /* High word of |x|. */
-	GET_HIGH_WORD(ix,x);
-	ix &= 0x7fffffff;
-
-    /* |x| in [0,22] */
-	if (ix < 0x40360000) {
-	    /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
-		if(ix<0x3fd62e43) {
-		    if (ix<0x3c800000)			/* cosh(tiny) = 1 */
-		      return one;
-		    t = __expm1(fabs(x));
-		    w = one+t;
-		    return one+(t*t)/(w+w);
-		}
-
-	    /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
-		t = __ieee754_exp(fabs(x));
-		return half*t+half/t;
-	}
-
-    /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
-	if (ix < 0x40862e42)  return half*__ieee754_exp(fabs(x));
-
-    /* |x| in [log(maxdouble), overflowthresold] */
-	int64_t fix;
-	EXTRACT_WORDS64(fix, x);
-	fix &= UINT64_C(0x7fffffffffffffff);
-	if (fix <= UINT64_C(0x408633ce8fb9f87d)) {
-	    w = __ieee754_exp(half*fabs(x));
-	    t = half*w;
-	    return t*w;
-	}
-
-    /* x is INF or NaN */
-	if(ix>=0x7ff00000) return x*x;
-
-    /* |x| > overflowthresold, cosh(x) overflow */
-	return huge*huge;
-}
-libm_alias_finite (__ieee754_cosh, __cosh)
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/e_fmod.c b/sysdeps/ieee754/dbl-64/wordsize-64/e_fmod.c
deleted file mode 100644
index 52a8687448..0000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/e_fmod.c
+++ /dev/null
@@ -1,106 +0,0 @@
-/* Rewritten for 64-bit machines by Ulrich Drepper <drepper@gmail.com>.  */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-/*
- * __ieee754_fmod(x,y)
- * Return x mod y in exact arithmetic
- * Method: shift and subtract
- */
-
-#include <math.h>
-#include <math_private.h>
-#include <stdint.h>
-#include <libm-alias-finite.h>
-
-static const double one = 1.0, Zero[] = {0.0, -0.0,};
-
-double
-__ieee754_fmod (double x, double y)
-{
-	int32_t n,ix,iy;
-	int64_t hx,hy,hz,sx,i;
-
-	EXTRACT_WORDS64(hx,x);
-	EXTRACT_WORDS64(hy,y);
-	sx = hx&UINT64_C(0x8000000000000000);	/* sign of x */
-	hx ^=sx;				/* |x| */
-	hy &= UINT64_C(0x7fffffffffffffff);	/* |y| */
-
-    /* purge off exception values */
-	if(__builtin_expect(hy==0
-			    || hx >= UINT64_C(0x7ff0000000000000)
-			    || hy > UINT64_C(0x7ff0000000000000), 0))
-	  /* y=0,or x not finite or y is NaN */
-	    return (x*y)/(x*y);
-	if(__builtin_expect(hx<=hy, 0)) {
-	    if(hx<hy) return x;	/* |x|<|y| return x */
-	    return Zero[(uint64_t)sx>>63];	/* |x|=|y| return x*0*/
-	}
-
-    /* determine ix = ilogb(x) */
-	if(__builtin_expect(hx<UINT64_C(0x0010000000000000), 0)) {
-	  /* subnormal x */
-	  for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
-	} else ix = (hx>>52)-1023;
-
-    /* determine iy = ilogb(y) */
-	if(__builtin_expect(hy<UINT64_C(0x0010000000000000), 0)) {	/* subnormal y */
-	  for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
-	} else iy = (hy>>52)-1023;
-
-    /* set up hx, hy and align y to x */
-	if(__builtin_expect(ix >= -1022, 1))
-	    hx = UINT64_C(0x0010000000000000)|(UINT64_C(0x000fffffffffffff)&hx);
-	else {		/* subnormal x, shift x to normal */
-	    n = -1022-ix;
-	    hx<<=n;
-	}
-	if(__builtin_expect(iy >= -1022, 1))
-	    hy = UINT64_C(0x0010000000000000)|(UINT64_C(0x000fffffffffffff)&hy);
-	else {		/* subnormal y, shift y to normal */
-	    n = -1022-iy;
-	    hy<<=n;
-	}
-
-    /* fix point fmod */
-	n = ix - iy;
-	while(n--) {
-	    hz=hx-hy;
-	    if(hz<0){hx = hx+hx;}
-	    else {
-		if(hz==0)		/* return sign(x)*0 */
-		    return Zero[(uint64_t)sx>>63];
-		hx = hz+hz;
-	    }
-	}
-	hz=hx-hy;
-	if(hz>=0) {hx=hz;}
-
-    /* convert back to floating value and restore the sign */
-	if(hx==0)			/* return sign(x)*0 */
-	    return Zero[(uint64_t)sx>>63];
-	while(hx<UINT64_C(0x0010000000000000)) {	/* normalize x */
-	    hx = hx+hx;
-	    iy -= 1;
-	}
-	if(__builtin_expect(iy>= -1022, 1)) {	/* normalize output */
-	  hx = ((hx-UINT64_C(0x0010000000000000))|((uint64_t)(iy+1023)<<52));
-	    INSERT_WORDS64(x,hx|sx);
-	} else {		/* subnormal output */
-	    n = -1022 - iy;
-	    hx>>=n;
-	    INSERT_WORDS64(x,hx|sx);
-	    x *= one;		/* create necessary signal */
-	}
-	return x;		/* exact output */
-}
-libm_alias_finite (__ieee754_fmod, __fmod)
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/e_log10.c b/sysdeps/ieee754/dbl-64/wordsize-64/e_log10.c
deleted file mode 100644
index b89064fb7c..0000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/e_log10.c
+++ /dev/null
@@ -1,90 +0,0 @@
-/* @(#)e_log10.c 5.1 93/09/24 */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-/* __ieee754_log10(x)
- * Return the base 10 logarithm of x
- *
- * Method :
- *	Let log10_2hi = leading 40 bits of log10(2) and
- *	    log10_2lo = log10(2) - log10_2hi,
- *	    ivln10   = 1/log(10) rounded.
- *	Then
- *		n = ilogb(x),
- *		if(n<0)  n = n+1;
- *		x = scalbn(x,-n);
- *		log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x))
- *
- * Note 1:
- *	To guarantee log10(10**n)=n, where 10**n is normal, the rounding
- *	mode must set to Round-to-Nearest.
- * Note 2:
- *	[1/log(10)] rounded to 53 bits has error  .198   ulps;
- *	log10 is monotonic at all binary break points.
- *
- * Special cases:
- *	log10(x) is NaN with signal if x < 0;
- *	log10(+INF) is +INF with no signal; log10(0) is -INF with signal;
- *	log10(NaN) is that NaN with no signal;
- *	log10(10**N) = N  for N=0,1,...,22.
- *
- * Constants:
- * The hexadecimal values are the intended ones for the following constants.
- * The decimal values may be used, provided that the compiler will convert
- * from decimal to binary accurately enough to produce the hexadecimal values
- * shown.
- */
-
-#include <math.h>
-#include <fix-int-fp-convert-zero.h>
-#include <math_private.h>
-#include <stdint.h>
-#include <libm-alias-finite.h>
-
-static const double two54 = 1.80143985094819840000e+16;		/* 0x4350000000000000 */
-static const double ivln10 = 4.34294481903251816668e-01;	/* 0x3FDBCB7B1526E50E */
-static const double log10_2hi = 3.01029995663611771306e-01;	/* 0x3FD34413509F6000 */
-static const double log10_2lo = 3.69423907715893078616e-13;	/* 0x3D59FEF311F12B36 */
-
-double
-__ieee754_log10 (double x)
-{
-  double y, z;
-  int64_t i, hx;
-  int32_t k;
-
-  EXTRACT_WORDS64 (hx, x);
-
-  k = 0;
-  if (hx < INT64_C(0x0010000000000000))
-    {				/* x < 2**-1022  */
-      if (__glibc_unlikely ((hx & UINT64_C(0x7fffffffffffffff)) == 0))
-	return -two54 / fabs (x);	/* log(+-0)=-inf */
-      if (__glibc_unlikely (hx < 0))
-	return (x - x) / (x - x);	/* log(-#) = NaN */
-      k -= 54;
-      x *= two54;		/* subnormal number, scale up x */
-      EXTRACT_WORDS64 (hx, x);
-    }
-  /* scale up resulted in a NaN number  */
-  if (__glibc_unlikely (hx >= UINT64_C(0x7ff0000000000000)))
-    return x + x;
-  k += (hx >> 52) - 1023;
-  i = ((uint64_t) k & UINT64_C(0x8000000000000000)) >> 63;
-  hx = (hx & UINT64_C(0x000fffffffffffff)) | ((0x3ff - i) << 52);
-  y = (double) (k + i);
-  if (FIX_INT_FP_CONVERT_ZERO && y == 0.0)
-    y = 0.0;
-  INSERT_WORDS64 (x, hx);
-  z = y * log10_2lo + ivln10 * __ieee754_log (x);
-  return z + y * log10_2hi;
-}
-libm_alias_finite (__ieee754_log10, __log10)
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c
deleted file mode 100644
index f6ddf4aaee..0000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c
+++ /dev/null
@@ -1,66 +0,0 @@
-/* Copyright (C) 2011-2021 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
-
-   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
-   <https://www.gnu.org/licenses/>.  */
-
-#include <inttypes.h>
-#include <math.h>
-#include <math_private.h>
-#include <libm-alias-double.h>
-
-/*
- * for non-zero, finite x
- *	x = frexp(arg,&exp);
- * return a double fp quantity x such that 0.5 <= |x| <1.0
- * and the corresponding binary exponent "exp". That is
- *	arg = x*2^exp.
- * If arg is inf, 0.0, or NaN, then frexp(arg,&exp) returns arg
- * with *exp=0.
- */
-
-
-double
-__frexp (double x, int *eptr)
-{
-  int64_t ix;
-  EXTRACT_WORDS64 (ix, x);
-  int32_t ex = 0x7ff & (ix >> 52);
-  int e = 0;
-
-  if (__glibc_likely (ex != 0x7ff && x != 0.0))
-    {
-      /* Not zero and finite.  */
-      e = ex - 1022;
-      if (__glibc_unlikely (ex == 0))
-	{
-	  /* Subnormal.  */
-	  x *= 0x1p54;
-	  EXTRACT_WORDS64 (ix, x);
-	  ex = 0x7ff & (ix >> 52);
-	  e = ex - 1022 - 54;
-	}
-
-      ix = (ix & INT64_C (0x800fffffffffffff)) | INT64_C (0x3fe0000000000000);
-      INSERT_WORDS64 (x, ix);
-    }
-  else
-    /* Quiet signaling NaNs.  */
-    x += x;
-
-  *eptr = e;
-  return x;
-}
-libm_alias_double (__frexp, frexp)
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_getpayload.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_getpayload.c
deleted file mode 100644
index 5e4ccd9ad1..0000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_getpayload.c
+++ /dev/null
@@ -1,38 +0,0 @@
-/* Get NaN payload.  dbl-64/wordsize-64 version.
-   Copyright (C) 2016-2021 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-
-   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
-   <https://www.gnu.org/licenses/>.  */
-
-#include <math.h>
-#include <math_private.h>
-#include <libm-alias-double.h>
-#include <stdint.h>
-#include <fix-int-fp-convert-zero.h>
-
-double
-__getpayload (const double *x)
-{
-  uint64_t ix;
-  EXTRACT_WORDS64 (ix, *x);
-  if ((ix & 0x7ff0000000000000ULL) != 0x7ff0000000000000ULL
-      || (ix & 0xfffffffffffffULL) == 0)
-    return -1;
-  ix &= 0x7ffffffffffffULL;
-  if (FIX_INT_FP_CONVERT_ZERO && ix == 0)
-    return 0.0f;
-  return (double) ix;
-}
-libm_alias_double (__getpayload, getpayload)
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_issignaling.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_issignaling.c
deleted file mode 100644
index 5fb6fbc7d4..0000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_issignaling.c
+++ /dev/null
@@ -1,43 +0,0 @@
-/* Test for signaling NaN.
-   Copyright (C) 2013-2021 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-
-   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
-   <https://www.gnu.org/licenses/>.  */
-
-#include <math.h>
-#include <math_private.h>
-#include <nan-high-order-bit.h>
-
-int
-__issignaling (double x)
-{
-  uint64_t xi;
-  EXTRACT_WORDS64 (xi, x);
-#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN
-  /* We only have to care about the high-order bit of x's significand, because
-     having it set (sNaN) already makes the significand different from that
-     used to designate infinity.  */
-  return (xi & UINT64_C (0x7ff8000000000000)) == UINT64_C (0x7ff8000000000000);
-#else
-  /* To keep the following comparison simple, toggle the quiet/signaling bit,
-     so that it is set for sNaNs.  This is inverse to IEEE 754-2008 (as well as
-     common practice for IEEE 754-1985).  */
-  xi ^= UINT64_C (0x0008000000000000);
-  /* We have to compare for greater (instead of greater or equal), because x's
-     significand being all-zero designates infinity not NaN.  */
-  return (xi & UINT64_C (0x7fffffffffffffff)) > UINT64_C (0x7ff8000000000000);
-#endif
-}
-libm_hidden_def (__issignaling)
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c
deleted file mode 100644
index 7020fd0156..0000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c
+++ /dev/null
@@ -1,85 +0,0 @@
-/* Round double value to long long int.
-   Copyright (C) 1997-2021 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
-   <https://www.gnu.org/licenses/>.  */
-
-#define lround __hidden_lround
-#define __lround __hidden___lround
-
-#include <fenv.h>
-#include <limits.h>
-#include <math.h>
-#include <sysdep.h>
-
-#include <math_private.h>
-#include <libm-alias-double.h>
-#include <fix-fp-int-convert-overflow.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
-    {
-#ifdef FE_INVALID
-      /* The number is too large.  Unless it rounds to LLONG_MIN,
-	 FE_INVALID must be raised and the return value is
-	 unspecified.  */
-      if (FIX_DBL_LLONG_CONVERT_OVERFLOW && x != (double) LLONG_MIN)
-	{
-	  feraiseexcept (FE_INVALID);
-	  return sign == 1 ? LLONG_MAX : LLONG_MIN;
-	}
-#endif
-      return (long long int) x;
-    }
-
-  return sign * result;
-}
-
-libm_alias_double (__llround, llround)
-
-/* long has the same width as long long on LP64 machines, so use an alias.  */
-#undef lround
-#undef __lround
-#ifdef _LP64
-strong_alias (__llround, __lround)
-libm_alias_double (__lround, lround)
-#endif
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c
deleted file mode 100644
index 5284c4da09..0000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c
+++ /dev/null
@@ -1,97 +0,0 @@
-/* Round double value to long int.
-   Copyright (C) 1997-2021 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-
-   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
-   <https://www.gnu.org/licenses/>.  */
-
-#include <fenv.h>
-#include <limits.h>
-#include <math.h>
-
-#include <math_private.h>
-#include <libm-alias-double.h>
-#include <fix-fp-int-convert-overflow.h>
-
-/* For LP64, lround is an alias for llround.  */
-#ifndef _LP64
-
-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);
-#ifdef FE_INVALID
-	  if (sizeof (long int) == 4
-	      && sign == 1
-	      && result == LONG_MIN)
-	    /* Rounding brought the value out of range.  */
-	    feraiseexcept (FE_INVALID);
-#endif
-	}
-    }
-  else
-    {
-      /* The number is too large.  Unless it rounds to LONG_MIN,
-	 FE_INVALID must be raised and the return value is
-	 unspecified.  */
-#ifdef FE_INVALID
-      if (FIX_DBL_LONG_CONVERT_OVERFLOW
-	  && !(sign == -1
-	       && (sizeof (long int) == 4
-		   ? x > (double) LONG_MIN - 0.5
-		   : x >= (double) LONG_MIN)))
-	{
-	  feraiseexcept (FE_INVALID);
-	  return sign == 1 ? LONG_MAX : LONG_MIN;
-	}
-      else if (!FIX_DBL_LONG_CONVERT_OVERFLOW
-	  && sizeof (long int) == 4
-	  && x <= (double) LONG_MIN - 0.5)
-	{
-	  /* If truncation produces LONG_MIN, the cast will not raise
-	     the exception, but may raise "inexact".  */
-	  feraiseexcept (FE_INVALID);
-	  return LONG_MIN;
-	}
-#endif
-      return (long int) x;
-    }
-
-  return sign * result;
-}
-
-libm_alias_double (__lround, lround)
-
-#endif
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_modf.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_modf.c
deleted file mode 100644
index 8d14e78ef0..0000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_modf.c
+++ /dev/null
@@ -1,65 +0,0 @@
-/* Rewritten for 64-bit machines by Ulrich Drepper <drepper@gmail.com>.  */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-/*
- * modf(double x, double *iptr)
- * return fraction part of x, and return x's integral part in *iptr.
- * Method:
- *	Bit twiddling.
- *
- * Exception:
- *	No exception.
- */
-
-#include <math.h>
-#include <math_private.h>
-#include <libm-alias-double.h>
-#include <stdint.h>
-
-static const double one = 1.0;
-
-double
-__modf(double x, double *iptr)
-{
-	int64_t i0;
-	int32_t j0;
-	EXTRACT_WORDS64(i0,x);
-	j0 = ((i0>>52)&0x7ff)-0x3ff;	/* exponent of x */
-	if(j0<52) {			/* integer part in x */
-	    if(j0<0) {			/* |x|<1 */
-		/* *iptr = +-0 */
-		INSERT_WORDS64(*iptr,i0&UINT64_C(0x8000000000000000));
-		return x;
-	    } else {
-		uint64_t i = UINT64_C(0x000fffffffffffff)>>j0;
-		if((i0&i)==0) {		/* x is integral */
-		    *iptr = x;
-		    /* return +-0 */
-		    INSERT_WORDS64(x,i0&UINT64_C(0x8000000000000000));
-		    return x;
-		} else {
-		    INSERT_WORDS64(*iptr,i0&(~i));
-		    return x - *iptr;
-		}
-	    }
-	} else { /* no fraction part */
-	    *iptr = x*one;
-	    /* We must handle NaNs separately.  */
-	    if (j0 == 0x400 && (i0 & UINT64_C(0xfffffffffffff)))
-	      return x*one;
-	    INSERT_WORDS64(x,i0&UINT64_C(0x8000000000000000));	/* return +-0 */
-	    return x;
-	}
-}
-#ifndef __modf
-libm_alias_double (__modf, modf)
-#endif
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c
deleted file mode 100644
index cbaa7f79a2..0000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c
+++ /dev/null
@@ -1,111 +0,0 @@
-/* Compute remainder and a congruent to the quotient.
-   Copyright (C) 1997-2021 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
-   <https://www.gnu.org/licenses/>.  */
-
-#include <math.h>
-
-#include <math_private.h>
-#include <libm-alias-double.h>
-#include <stdint.h>
-
-static const double zero = 0.0;
-
-
-double
-__remquo (double x, double y, int *quo)
-{
-  int64_t hx, hy;
-  uint64_t sx, qs;
-  int cquo;
-
-  EXTRACT_WORDS64 (hx, x);
-  EXTRACT_WORDS64 (hy, y);
-  sx = hx & UINT64_C(0x8000000000000000);
-  qs = sx ^ (hy & UINT64_C(0x8000000000000000));
-  hy &= UINT64_C(0x7fffffffffffffff);
-  hx &= UINT64_C(0x7fffffffffffffff);
-
-  /* Purge off exception values.  */
-  if (__glibc_unlikely (hy == 0))
-    return (x * y) / (x * y);			/* y = 0 */
-  if (__builtin_expect (hx >= UINT64_C(0x7ff0000000000000) /* x not finite */
-			|| hy > UINT64_C(0x7ff0000000000000), 0))/* y is NaN */
-    return (x * y) / (x * y);
-
-  if (hy <= UINT64_C(0x7fbfffffffffffff))
-    x = __ieee754_fmod (x, 8 * y);		/* now x < 8y */
-
-  if (__glibc_unlikely (hx == hy))
-    {
-      *quo = qs ? -1 : 1;
-      return zero * x;
-    }
-
-  x = fabs (x);
-  INSERT_WORDS64 (y, hy);
-  cquo = 0;
-
-  if (hy <= UINT64_C(0x7fcfffffffffffff) && x >= 4 * y)
-    {
-      x -= 4 * y;
-      cquo += 4;
-    }
-  if (hy <= UINT64_C(0x7fdfffffffffffff) && x >= 2 * y)
-    {
-      x -= 2 * y;
-      cquo += 2;
-    }
-
-  if (hy < UINT64_C(0x0020000000000000))
-    {
-      if (x + x > y)
-	{
-	  x -= y;
-	  ++cquo;
-	  if (x + x >= y)
-	    {
-	      x -= y;
-	      ++cquo;
-	    }
-	}
-    }
-  else
-    {
-      double y_half = 0.5 * y;
-      if (x > y_half)
-	{
-	  x -= y;
-	  ++cquo;
-	  if (x >= y_half)
-	    {
-	      x -= y;
-	      ++cquo;
-	    }
-	}
-    }
-
-  *quo = qs ? -cquo : cquo;
-
-  /* Ensure correct sign of zero result in round-downward mode.  */
-  if (x == 0.0)
-    x = 0.0;
-  if (sx)
-    x = -x;
-  return x;
-}
-libm_alias_double (__remquo, remquo)
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_roundeven.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_roundeven.c
deleted file mode 100644
index 289804c840..0000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_roundeven.c
+++ /dev/null
@@ -1,71 +0,0 @@
-/* Round to nearest integer value, rounding halfway cases to even.
-   dbl-64/wordsize-64 version.
-   Copyright (C) 2016-2021 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-
-   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
-   <https://www.gnu.org/licenses/>.  */
-
-#include <math.h>
-#include <math_private.h>
-#include <libm-alias-double.h>
-#include <stdint.h>
-
-#define BIAS 0x3ff
-#define MANT_DIG 53
-#define MAX_EXP (2 * BIAS + 1)
-
-double
-__roundeven (double x)
-{
-  uint64_t ix, ux;
-  EXTRACT_WORDS64 (ix, x);
-  ux = ix & 0x7fffffffffffffffULL;
-  int exponent = ux >> (MANT_DIG - 1);
-  if (exponent >= BIAS + MANT_DIG - 1)
-    {
-      /* Integer, infinity or NaN.  */
-      if (exponent == MAX_EXP)
-	/* Infinity or NaN; quiet signaling NaNs.  */
-	return x + x;
-      else
-	return x;
-    }
-  else if (exponent >= BIAS)
-    {
-      /* At least 1; not necessarily an integer.  Locate the bits with
-	 exponents 0 and -1 (when the unbiased exponent is 0, the bit
-	 with exponent 0 is implicit, but as the bias is odd it is OK
-	 to take it from the low bit of the exponent).  */
-      int int_pos = (BIAS + MANT_DIG - 1) - exponent;
-      int half_pos = int_pos - 1;
-      uint64_t half_bit = 1ULL << half_pos;
-      uint64_t int_bit = 1ULL << int_pos;
-      if ((ix & (int_bit | (half_bit - 1))) != 0)
-	/* Carry into the exponent works correctly.  No need to test
-	   whether HALF_BIT is set.  */
-	ix += half_bit;
-      ix &= ~(int_bit - 1);
-    }
-  else if (exponent == BIAS - 1 && ux > 0x3fe0000000000000ULL)
-    /* Interval (0.5, 1).  */
-    ix = (ix & 0x8000000000000000ULL) | 0x3ff0000000000000ULL;
-  else
-    /* Rounds to 0.  */
-    ix &= 0x8000000000000000ULL;
-  INSERT_WORDS64 (x, ix);
-  return x;
-}
-hidden_def (__roundeven)
-libm_alias_double (__roundeven, roundeven)
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_scalbln.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_scalbln.c
deleted file mode 100644
index 071c9d7794..0000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_scalbln.c
+++ /dev/null
@@ -1,60 +0,0 @@
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-/*
- * scalbn (double x, int n)
- * scalbn(x,n) returns x* 2**n  computed by  exponent
- * manipulation rather than by actually performing an
- * exponentiation or a multiplication.
- */
-
-#include <math.h>
-#include <math_private.h>
-
-static const double
-two54   =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
-twom54  =  5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
-huge   = 1.0e+300,
-tiny   = 1.0e-300;
-
-double
-__scalbln (double x, long int n)
-{
-	int64_t ix;
-	int64_t k;
-	EXTRACT_WORDS64(ix,x);
-	k = (ix >> 52) & 0x7ff;			/* extract exponent */
-	if (__builtin_expect(k==0, 0)) {	/* 0 or subnormal x */
-	    if ((ix & UINT64_C(0xfffffffffffff))==0) return x; /* +-0 */
-	    x *= two54;
-	    EXTRACT_WORDS64(ix,x);
-	    k = ((ix >> 52) & 0x7ff) - 54;
-	    }
-	if (__builtin_expect(k==0x7ff, 0)) return x+x;	/* NaN or Inf */
-	if (__builtin_expect(n< -50000, 0))
-	  return tiny*copysign(tiny,x); /*underflow*/
-	if (__builtin_expect(n> 50000 || k+n > 0x7fe, 0))
-	  return huge*copysign(huge,x); /* overflow  */
-	/* Now k and n are bounded we know that k = k+n does not
-	   overflow.  */
-	k = k+n;
-	if (__builtin_expect(k > 0, 1))		/* normal result */
-	    {INSERT_WORDS64(x,(ix&UINT64_C(0x800fffffffffffff))|(k<<52));
-	      return x;}
-	if (k <= -54)
-	  return tiny*copysign(tiny,x);	/*underflow*/
-	k += 54;				/* subnormal result */
-	INSERT_WORDS64(x,(ix&INT64_C(0x800fffffffffffff))|(k<<52));
-	return x*twom54;
-}
-#ifdef NO_LONG_DOUBLE
-strong_alias (__scalbln, __scalblnl)
-#endif
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_scalbn.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_scalbn.c
deleted file mode 100644
index 4491227f3e..0000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_scalbn.c
+++ /dev/null
@@ -1,60 +0,0 @@
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-/*
- * scalbn (double x, int n)
- * scalbn(x,n) returns x* 2**n  computed by  exponent
- * manipulation rather than by actually performing an
- * exponentiation or a multiplication.
- */
-
-#include <math.h>
-#include <math_private.h>
-
-static const double
-two54   =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
-twom54  =  5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
-huge   = 1.0e+300,
-tiny   = 1.0e-300;
-
-double
-__scalbn (double x, int n)
-{
-	int64_t ix;
-	int64_t k;
-	EXTRACT_WORDS64(ix,x);
-	k = (ix >> 52) & 0x7ff;			/* extract exponent */
-	if (__builtin_expect(k==0, 0)) {	/* 0 or subnormal x */
-	    if ((ix & UINT64_C(0xfffffffffffff))==0) return x; /* +-0 */
-	    x *= two54;
-	    EXTRACT_WORDS64(ix,x);
-	    k = ((ix >> 52) & 0x7ff) - 54;
-	    }
-	if (__builtin_expect(k==0x7ff, 0)) return x+x;	/* NaN or Inf */
-	if (__builtin_expect(n< -50000, 0))
-	  return tiny*copysign(tiny,x); /*underflow*/
-	if (__builtin_expect(n> 50000 || k+n > 0x7fe, 0))
-	  return huge*copysign(huge,x); /* overflow  */
-	/* Now k and n are bounded we know that k = k+n does not
-	   overflow.  */
-	k = k+n;
-	if (__builtin_expect(k > 0, 1))		/* normal result */
-	    {INSERT_WORDS64(x,(ix&UINT64_C(0x800fffffffffffff))|(k<<52));
-	      return x;}
-	if (k <= -54)
-	  return tiny*copysign(tiny,x);	/*underflow*/
-	k += 54;				/* subnormal result */
-	INSERT_WORDS64(x,(ix&INT64_C(0x800fffffffffffff))|(k<<52));
-	return x*twom54;
-}
-#ifdef NO_LONG_DOUBLE
-strong_alias (__scalbn, __scalbnl)
-#endif
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_setpayload_main.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_setpayload_main.c
deleted file mode 100644
index b622a50985..0000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_setpayload_main.c
+++ /dev/null
@@ -1,54 +0,0 @@
-/* Set NaN payload.  dbl-64/wordsize-64 version.
-   Copyright (C) 2016-2021 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-
-   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
-   <https://www.gnu.org/licenses/>.  */
-
-#include <math.h>
-#include <math_private.h>
-#include <libm-alias-double.h>
-#include <nan-high-order-bit.h>
-#include <stdint.h>
-
-#define SET_HIGH_BIT (HIGH_ORDER_BIT_IS_SET_FOR_SNAN ? SIG : !SIG)
-#define BIAS 0x3ff
-#define PAYLOAD_DIG 51
-#define EXPLICIT_MANT_DIG 52
-
-int
-FUNC (double *x, double payload)
-{
-  uint64_t ix;
-  EXTRACT_WORDS64 (ix, payload);
-  int exponent = ix >> EXPLICIT_MANT_DIG;
-  /* Test if argument is (a) negative or too large; (b) too small,
-     except for 0 when allowed; (c) not an integer.  */
-  if (exponent >= BIAS + PAYLOAD_DIG
-      || (exponent < BIAS && !(SET_HIGH_BIT && ix == 0))
-      || (ix & ((1ULL << (BIAS + EXPLICIT_MANT_DIG - exponent)) - 1)) != 0)
-    {
-      INSERT_WORDS64 (*x, 0);
-      return 1;
-    }
-  if (ix != 0)
-    {
-      ix &= (1ULL << EXPLICIT_MANT_DIG) - 1;
-      ix |= 1ULL << EXPLICIT_MANT_DIG;
-      ix >>= BIAS + EXPLICIT_MANT_DIG - exponent;
-    }
-  ix |= 0x7ff0000000000000ULL | (SET_HIGH_BIT ? 0x8000000000000ULL : 0);
-  INSERT_WORDS64 (*x, ix);
-  return 0;
-}
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_totalorder.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_totalorder.c
deleted file mode 100644
index 1e83b45008..0000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_totalorder.c
+++ /dev/null
@@ -1,76 +0,0 @@
-/* Total order operation.  dbl-64/wordsize-64 version.
-   Copyright (C) 2016-2021 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-
-   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
-   <https://www.gnu.org/licenses/>.  */
-
-#include <math.h>
-#include <math_private.h>
-#include <nan-high-order-bit.h>
-#include <libm-alias-double.h>
-#include <stdint.h>
-#include <shlib-compat.h>
-#include <first-versions.h>
-
-int
-__totalorder (const double *x, const double *y)
-{
-  int64_t ix, iy;
-  EXTRACT_WORDS64 (ix, *x);
-  EXTRACT_WORDS64 (iy, *y);
-#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN
-  /* For the preferred quiet NaN convention, this operation is a
-     comparison of the representations of the arguments interpreted as
-     sign-magnitude integers.  If both arguments are NaNs, invert the
-     quiet/signaling bit so comparing that way works.  */
-  if ((ix & 0x7fffffffffffffffULL) > 0x7ff0000000000000ULL
-      && (iy & 0x7fffffffffffffffULL) > 0x7ff0000000000000ULL)
-    {
-      ix ^= 0x0008000000000000ULL;
-      iy ^= 0x0008000000000000ULL;
-    }
-#endif
-  uint64_t ix_sign = ix >> 63;
-  uint64_t iy_sign = iy >> 63;
-  ix ^= ix_sign >> 1;
-  iy ^= iy_sign >> 1;
-  return ix <= iy;
-}
-#ifdef SHARED
-# define CONCATX(x, y) x ## y
-# define CONCAT(x, y) CONCATX (x, y)
-# define UNIQUE_ALIAS(name) CONCAT (name, __COUNTER__)
-# define do_symbol(orig_name, name, aliasname)		\
-  strong_alias (orig_name, name)			\
-  versioned_symbol (libm, name, aliasname, GLIBC_2_31)
-# undef weak_alias
-# define weak_alias(name, aliasname)			\
-  do_symbol (name, UNIQUE_ALIAS (name), aliasname);
-#endif
-libm_alias_double (__totalorder, totalorder)
-#if SHLIB_COMPAT (libm, GLIBC_2_25, GLIBC_2_31)
-int
-attribute_compat_text_section
-__totalorder_compat (double x, double y)
-{
-  return __totalorder (&x, &y);
-}
-#undef do_symbol
-#define do_symbol(orig_name, name, aliasname)			\
-  strong_alias (orig_name, name)				\
-  compat_symbol (libm, name, aliasname,				\
-		 CONCAT (FIRST_VERSION_libm_, aliasname))
-libm_alias_double (__totalorder_compat, totalorder)
-#endif
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_totalordermag.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_totalordermag.c
deleted file mode 100644
index 2da739782a..0000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_totalordermag.c
+++ /dev/null
@@ -1,73 +0,0 @@
-/* Total order operation on absolute values.  dbl-64/wordsize-64 version.
-   Copyright (C) 2016-2021 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-
-   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
-   <https://www.gnu.org/licenses/>.  */
-
-#include <math.h>
-#include <math_private.h>
-#include <nan-high-order-bit.h>
-#include <libm-alias-double.h>
-#include <stdint.h>
-#include <shlib-compat.h>
-#include <first-versions.h>
-
-int
-__totalordermag (const double *x, const double *y)
-{
-  uint64_t ix, iy;
-  EXTRACT_WORDS64 (ix, *x);
-  EXTRACT_WORDS64 (iy, *y);
-  ix &= 0x7fffffffffffffffULL;
-  iy &= 0x7fffffffffffffffULL;
-#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN
-  /* For the preferred quiet NaN convention, this operation is a
-     comparison of the representations of the absolute values of the
-     arguments.  If both arguments are NaNs, invert the
-     quiet/signaling bit so comparing that way works.  */
-  if (ix > 0x7ff0000000000000ULL && iy > 0x7ff0000000000000ULL)
-    {
-      ix ^= 0x0008000000000000ULL;
-      iy ^= 0x0008000000000000ULL;
-    }
-#endif
-  return ix <= iy;
-}
-#ifdef SHARED
-# define CONCATX(x, y) x ## y
-# define CONCAT(x, y) CONCATX (x, y)
-# define UNIQUE_ALIAS(name) CONCAT (name, __COUNTER__)
-# define do_symbol(orig_name, name, aliasname)		\
-  strong_alias (orig_name, name)			\
-  versioned_symbol (libm, name, aliasname, GLIBC_2_31)
-# undef weak_alias
-# define weak_alias(name, aliasname)			\
-  do_symbol (name, UNIQUE_ALIAS (name), aliasname);
-#endif
-libm_alias_double (__totalordermag, totalordermag)
-#if SHLIB_COMPAT (libm, GLIBC_2_25, GLIBC_2_31)
-int
-attribute_compat_text_section
-__totalordermag_compat (double x, double y)
-{
-  return __totalordermag (&x, &y);
-}
-#undef do_symbol
-#define do_symbol(orig_name, name, aliasname)			\
-  strong_alias (orig_name, name)				\
-  compat_symbol (libm, name, aliasname,				\
-		 CONCAT (FIRST_VERSION_libm_, aliasname))
-libm_alias_double (__totalordermag_compat, totalordermag)
-#endif