[2/2] Remove dbl-64/wordsize-64
Commit Message
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
I'd expect this also to remove the ieee754/dbl-64/wordsize-64 entries from
Implies files.
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
>
@@ -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)
@@ -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)
@@ -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)
@@ -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;
}
@@ -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)
@@ -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;
@@ -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)
@@ -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
@@ -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
@@ -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)
@@ -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)
{
@@ -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)
@@ -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)
@@ -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)
@@ -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;
}
@@ -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
@@ -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
deleted file mode 100644
@@ -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)
deleted file mode 100644
@@ -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)
deleted file mode 100644
@@ -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)
deleted file mode 100644
@@ -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)
deleted file mode 100644
@@ -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)
deleted file mode 100644
@@ -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)
deleted file mode 100644
@@ -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)
deleted file mode 100644
@@ -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
deleted file mode 100644
@@ -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
deleted file mode 100644
@@ -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
deleted file mode 100644
@@ -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)
deleted file mode 100644
@@ -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)
deleted file mode 100644
@@ -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
deleted file mode 100644
@@ -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
deleted file mode 100644
@@ -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;
-}
deleted file mode 100644
@@ -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
deleted file mode 100644
@@ -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