diff mbox series

[2/2] Remove dbl-64/wordsize-64

Message ID VE1PR08MB55992C16F841E86854BAC0A283C50@VE1PR08MB5599.eurprd08.prod.outlook.com
State Superseded
Headers show
Series [1/2] Remove dbl-64/wordsize-64 | expand

Commit Message

Wilco Dijkstra Dec. 16, 2020, 5 p.m. UTC
Remove the wordsize-64 implementations by merging them into the main dbl-64
directory. The second patch just moves all wordsize-64 files without any changes
(diff looks huge but is literally git mv -f wordsize-64/* .)

GLIBC regress passes on AArch64 and Arm. Buildmanyglibc passes. OK for commit?

---

Comments

Joseph Myers Dec. 16, 2020, 6:34 p.m. UTC | #1
I'd expect this also to remove the ieee754/dbl-64/wordsize-64 entries from 
Implies files.
Adhemerval Zanella Dec. 21, 2020, 5:52 p.m. UTC | #2
On 16/12/2020 14:00, Wilco Dijkstra via Libc-alpha wrote:
> Remove the wordsize-64 implementations by merging them into the main dbl-64
> directory. The second patch just moves all wordsize-64 files without any changes
> (diff looks huge but is literally git mv -f wordsize-64/* .)
> 
> GLIBC regress passes on AArch64 and Arm. Buildmanyglibc passes. OK for commit?

LGTM with the fix noted by Joseph (remove ieee754/dbl-64/wordsize-64 entries).

I did a sniff test on powerpc32 (which without any --with-cpu should set
FIX_INT_FP_CONVERT_ZERO to 1) and I saw no regressions.

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

> 
> ---
> 
> diff --git a/sysdeps/ieee754/dbl-64/e_acosh.c b/sysdeps/ieee754/dbl-64/e_acosh.c
> index 75df0ab5ef15a858c469370142ca119485337f33..a241366f308abb6e11da80f68d86998708d79847 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 6c78a3a4e9b5037f9f3976a5a98b46a1494ffe6c..4f41ca2c92b37263e5684f3e485db6675f2ba61f 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 f6a095ba82905f94bc834776ba0877497328e9ee..52a86874484011f567a6759324ce941a89e77625 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 44a4bd2faa9792c68ac883c19da2dbfb8070616f..b89064fb7c41dd857d216b911aeb3935605c6d38 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 c96a86966563043d184c7df9096aa41d41d4d830..b1d14354e05037c029cae989d044997273196ac7 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-2020 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 5a055be35a4e2f94c2691655e338981f8cefeb27..d541f0f1b6b00ed78d2ec6f05086f5c053841f2b 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-2020 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 85578a27c5ddab2bb41e0d0095fbb28a0b525e6e..c849d11ab1c8256a4190ad38a33ce39cb83b09c6 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 8e8f94e66ac49c428136039f3b54cf6862b376ce..1d9c6c5b1676920c951fdb56cf133bfc64195405 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 17bcb69d3110a9999076e0ae8d40b943e513d565..c0cebe57b798c910f2f3cdc02e86cbecb14285a3 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-2020 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 722511c64ac180a08c35d3f777b45dfc2335935e..8d14e78ef006e59dea0316221692dac572e0e139 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 4a55c6a3558ec381283acbd68f3c0d0794b43785..279285f40fff1cda31357d266131d752628f3558 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 ac8c64e229ca6590b9b4301b79c05c0f0dc08d5e..f6b592a368199679fb078d545771219ac794f46c 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-2020 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 0e3d732e48842bb69941f98a39afa457aff6d3c6..071c9d7794ad398006f3e4f050e2509555721b18 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 cf4d6846ee451d682a43a6849a36f50f4456d4b5..4491227f3e3b5cf431564146b4aadc9cc206339e 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 d43491c3d74a028dc34a26059094397e43f5e156..dda177c5ee7a7a53878c7db575e288d9a021870b 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-2020 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 c603c135103586d380a8f1ddf015ad878cc325fb..acc629a00ffbcb8ebcadc38caadfe46d3d8189b8 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-2020 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 82ea71af64d99c4cc2ff8b0bd7054ee16eee36a0..a60cf57129d80c50e6e8608dc252db68106cc47d 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-2020 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 a241366f308abb6e11da80f68d86998708d79847..0000000000000000000000000000000000000000
> --- 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 4f41ca2c92b37263e5684f3e485db6675f2ba61f..0000000000000000000000000000000000000000
> --- 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 52a86874484011f567a6759324ce941a89e77625..0000000000000000000000000000000000000000
> --- 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 b89064fb7c41dd857d216b911aeb3935605c6d38..0000000000000000000000000000000000000000
> --- 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 b1d14354e05037c029cae989d044997273196ac7..0000000000000000000000000000000000000000
> --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c
> +++ /dev/null
> @@ -1,66 +0,0 @@
> -/* Copyright (C) 2011-2020 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 d541f0f1b6b00ed78d2ec6f05086f5c053841f2b..0000000000000000000000000000000000000000
> --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_getpayload.c
> +++ /dev/null
> @@ -1,38 +0,0 @@
> -/* Get NaN payload.
> -   Copyright (C) 2016-2020 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 c849d11ab1c8256a4190ad38a33ce39cb83b09c6..0000000000000000000000000000000000000000
> --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_issignaling.c
> +++ /dev/null
> @@ -1,43 +0,0 @@
> -/* Test for signaling NaN.
> -   Copyright (C) 2013-2020 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 1d9c6c5b1676920c951fdb56cf133bfc64195405..0000000000000000000000000000000000000000
> --- 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-2020 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 c0cebe57b798c910f2f3cdc02e86cbecb14285a3..0000000000000000000000000000000000000000
> --- 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-2020 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 8d14e78ef006e59dea0316221692dac572e0e139..0000000000000000000000000000000000000000
> --- 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 279285f40fff1cda31357d266131d752628f3558..0000000000000000000000000000000000000000
> --- 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-2020 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 f6b592a368199679fb078d545771219ac794f46c..0000000000000000000000000000000000000000
> --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_roundeven.c
> +++ /dev/null
> @@ -1,70 +0,0 @@
> -/* Round to nearest integer value, rounding halfway cases to even.
> -   Copyright (C) 2016-2020 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 071c9d7794ad398006f3e4f050e2509555721b18..0000000000000000000000000000000000000000
> --- 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 4491227f3e3b5cf431564146b4aadc9cc206339e..0000000000000000000000000000000000000000
> --- 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 dda177c5ee7a7a53878c7db575e288d9a021870b..0000000000000000000000000000000000000000
> --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_setpayload_main.c
> +++ /dev/null
> @@ -1,54 +0,0 @@
> -/* Set NaN payload.
> -   Copyright (C) 2016-2020 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 acc629a00ffbcb8ebcadc38caadfe46d3d8189b8..0000000000000000000000000000000000000000
> --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_totalorder.c
> +++ /dev/null
> @@ -1,76 +0,0 @@
> -/* Total order operation.
> -   Copyright (C) 2016-2020 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 a60cf57129d80c50e6e8608dc252db68106cc47d..0000000000000000000000000000000000000000
> --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_totalordermag.c
> +++ /dev/null
> @@ -1,73 +0,0 @@
> -/* Total order operation on absolute values.
> -   Copyright (C) 2016-2020 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
>
diff mbox series

Patch

diff --git a/sysdeps/ieee754/dbl-64/e_acosh.c b/sysdeps/ieee754/dbl-64/e_acosh.c
index 75df0ab5ef15a858c469370142ca119485337f33..a241366f308abb6e11da80f68d86998708d79847 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 6c78a3a4e9b5037f9f3976a5a98b46a1494ffe6c..4f41ca2c92b37263e5684f3e485db6675f2ba61f 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 f6a095ba82905f94bc834776ba0877497328e9ee..52a86874484011f567a6759324ce941a89e77625 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 44a4bd2faa9792c68ac883c19da2dbfb8070616f..b89064fb7c41dd857d216b911aeb3935605c6d38 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 c96a86966563043d184c7df9096aa41d41d4d830..b1d14354e05037c029cae989d044997273196ac7 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-2020 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 5a055be35a4e2f94c2691655e338981f8cefeb27..d541f0f1b6b00ed78d2ec6f05086f5c053841f2b 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-2020 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 85578a27c5ddab2bb41e0d0095fbb28a0b525e6e..c849d11ab1c8256a4190ad38a33ce39cb83b09c6 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 8e8f94e66ac49c428136039f3b54cf6862b376ce..1d9c6c5b1676920c951fdb56cf133bfc64195405 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 17bcb69d3110a9999076e0ae8d40b943e513d565..c0cebe57b798c910f2f3cdc02e86cbecb14285a3 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-2020 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 722511c64ac180a08c35d3f777b45dfc2335935e..8d14e78ef006e59dea0316221692dac572e0e139 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 4a55c6a3558ec381283acbd68f3c0d0794b43785..279285f40fff1cda31357d266131d752628f3558 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 ac8c64e229ca6590b9b4301b79c05c0f0dc08d5e..f6b592a368199679fb078d545771219ac794f46c 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-2020 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 0e3d732e48842bb69941f98a39afa457aff6d3c6..071c9d7794ad398006f3e4f050e2509555721b18 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 cf4d6846ee451d682a43a6849a36f50f4456d4b5..4491227f3e3b5cf431564146b4aadc9cc206339e 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 d43491c3d74a028dc34a26059094397e43f5e156..dda177c5ee7a7a53878c7db575e288d9a021870b 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-2020 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 c603c135103586d380a8f1ddf015ad878cc325fb..acc629a00ffbcb8ebcadc38caadfe46d3d8189b8 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-2020 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 82ea71af64d99c4cc2ff8b0bd7054ee16eee36a0..a60cf57129d80c50e6e8608dc252db68106cc47d 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-2020 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 a241366f308abb6e11da80f68d86998708d79847..0000000000000000000000000000000000000000
--- 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 4f41ca2c92b37263e5684f3e485db6675f2ba61f..0000000000000000000000000000000000000000
--- 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 52a86874484011f567a6759324ce941a89e77625..0000000000000000000000000000000000000000
--- 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 b89064fb7c41dd857d216b911aeb3935605c6d38..0000000000000000000000000000000000000000
--- 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 b1d14354e05037c029cae989d044997273196ac7..0000000000000000000000000000000000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c
+++ /dev/null
@@ -1,66 +0,0 @@ 
-/* Copyright (C) 2011-2020 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 d541f0f1b6b00ed78d2ec6f05086f5c053841f2b..0000000000000000000000000000000000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_getpayload.c
+++ /dev/null
@@ -1,38 +0,0 @@ 
-/* Get NaN payload.
-   Copyright (C) 2016-2020 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 c849d11ab1c8256a4190ad38a33ce39cb83b09c6..0000000000000000000000000000000000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_issignaling.c
+++ /dev/null
@@ -1,43 +0,0 @@ 
-/* Test for signaling NaN.
-   Copyright (C) 2013-2020 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 1d9c6c5b1676920c951fdb56cf133bfc64195405..0000000000000000000000000000000000000000
--- 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-2020 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 c0cebe57b798c910f2f3cdc02e86cbecb14285a3..0000000000000000000000000000000000000000
--- 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-2020 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 8d14e78ef006e59dea0316221692dac572e0e139..0000000000000000000000000000000000000000
--- 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 279285f40fff1cda31357d266131d752628f3558..0000000000000000000000000000000000000000
--- 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-2020 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 f6b592a368199679fb078d545771219ac794f46c..0000000000000000000000000000000000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_roundeven.c
+++ /dev/null
@@ -1,70 +0,0 @@ 
-/* Round to nearest integer value, rounding halfway cases to even.
-   Copyright (C) 2016-2020 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 071c9d7794ad398006f3e4f050e2509555721b18..0000000000000000000000000000000000000000
--- 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 4491227f3e3b5cf431564146b4aadc9cc206339e..0000000000000000000000000000000000000000
--- 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 dda177c5ee7a7a53878c7db575e288d9a021870b..0000000000000000000000000000000000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_setpayload_main.c
+++ /dev/null
@@ -1,54 +0,0 @@ 
-/* Set NaN payload.
-   Copyright (C) 2016-2020 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 acc629a00ffbcb8ebcadc38caadfe46d3d8189b8..0000000000000000000000000000000000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_totalorder.c
+++ /dev/null
@@ -1,76 +0,0 @@ 
-/* Total order operation.
-   Copyright (C) 2016-2020 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 a60cf57129d80c50e6e8608dc252db68106cc47d..0000000000000000000000000000000000000000
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_totalordermag.c
+++ /dev/null
@@ -1,73 +0,0 @@ 
-/* Total order operation on absolute values.
-   Copyright (C) 2016-2020 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