Message ID | VE1PR08MB55992C16F841E86854BAC0A283C50@VE1PR08MB5599.eurprd08.prod.outlook.com |
---|---|
State | New |
Headers | show |
Series | [1/2] Remove dbl-64/wordsize-64 | expand |
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 >
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