mirrors/glibc
1
0
Fork 0
mirror of https://sourceware.org/git/glibc.git synced 2024-12-04 15:13:53 +08:00
Code Issues Actions Packages Projects Releases Wiki Activity

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

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

Reviewed-by: Adhemerval Zanella  <adhemerval.zanella@linaro.org>
This commit is contained in:
Wilco Dijkstra 2021-01-07 15:26:26 +00:00
parent caa884dda7
commit 9e97f239ea
40 changed files with 427 additions and 1848 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
sysdeps/aarch64/Implies
View File

@ -1,5 +1,4 @@
wordsize-64
ieee754/ldbl-128
ieee754/dbl-64/wordsize-64
ieee754/dbl-64
ieee754/flt-32

1
sysdeps/alpha/Implies
View File

@ -1,6 +1,5 @@
wordsize-64
# Alpha uses IEEE 754 single, double and quad precision floating point.
ieee754/ldbl-128
ieee754/dbl-64/wordsize-64
ieee754/dbl-64
ieee754/flt-32

50
sysdeps/ieee754/dbl-64/e_acosh.c
View File

@ -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 */
int64_t hx;
EXTRACT_WORDS64 (hx, x);
if (hx > INT64_C (0x4000000000000000))
{
return (x - x) / (x - x);
}
else if (hx >= 0x41b00000) /* x > 2**28 */
{
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)

75
sysdeps/ieee754/dbl-64/e_cosh.c
View File

@ -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] */
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| 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 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)

212
sysdeps/ieee754/dbl-64/e_fmod.c
View File

@ -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;
/* 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
{
for (ix = -1022, i = (hx << 11); i > 0; i <<= 1)
ix -= 1;
}
}
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)
/* 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;
}
else
{
for (iy = -1022, i = (hy << 11); i > 0; i <<= 1)
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 */
}
}
else
iy = (hy >> 20) - 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;
}
}
/* 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;
}
else
{
if ((hz | lz) == 0) /* return sign(x)*0 */
return Zero[(uint32_t) sx >> 31];
hx = hz + hz + (lz >> 31); lx = lz + lz;
}
}
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;
}
else if (n <= 31)
{
lx = (hx << (32 - n)) | (lx >> n); hx = sx;
}
else
{
lx = hx >> (n - 32); hx = sx;
}
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)

42
sysdeps/ieee754/dbl-64/e_log10.c
View File

@ -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;
}

81
sysdeps/ieee754/dbl-64/s_frexp.c
View File

@ -1,21 +1,28 @@
/* @(#)s_frexp.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* Copyright (C) 2011-2021 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
#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
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 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)

15
sysdeps/ieee754/dbl-64/s_getpayload.c
View File

@ -1,4 +1,4 @@
/* Get NaN payload. dbl-64 version.
/* Get NaN payload.
Copyright (C) 2016-2021 Free Software Foundation, Inc.
This file is part of the GNU C Library.
@ -16,22 +16,21 @@
License along with the GNU C Library; if not, see
<https://www.gnu.org/licenses/>. */
#include <fix-int-fp-convert-zero.h>
#include <math.h>
#include <math_private.h>
#include <libm-alias-double.h>
#include <stdint.h>
#include <fix-int-fp-convert-zero.h>
double
__getpayload (const double *x)
{
uint32_t hx, lx;
EXTRACT_WORDS (hx, lx, *x);
if ((hx & 0x7ff00000) != 0x7ff00000
|| ((hx & 0xfffff) | lx) == 0)
uint64_t ix;
EXTRACT_WORDS64 (ix, *x);
if ((ix & 0x7ff0000000000000ULL) != 0x7ff0000000000000ULL
|| (ix & 0xfffffffffffffULL) == 0)
return -1;
hx &= 0x7ffff;
uint64_t ix = ((uint64_t) hx << 32) | lx;
ix &= 0x7ffffffffffffULL;
if (FIX_INT_FP_CONVERT_ZERO && ix == 0)
return 0.0f;
return (double) ix;

14
sysdeps/ieee754/dbl-64/s_issignaling.c
View File

@ -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)

49
sysdeps/ieee754/dbl-64/s_llround.c
View File

@ -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;
i0 += UINT64_C(0x8000000000000) >> 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;
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

57
sysdeps/ieee754/dbl-64/s_lround.c
View File

@ -1,7 +1,6 @@
/* Round double value to long int.
Copyright (C) 1997-2021 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
@ -25,55 +24,41 @@
#include <libm-alias-double.h>
#include <fix-fp-int-convert-overflow.h>
/* For LP64, lround is an alias for llround. */
#ifndef _LP64
long int
__lround (double x)
{
int32_t j0;
uint32_t i1, i0;
int64_t i0;
long int result;
int sign;
EXTRACT_WORDS (i0, i1, x);
j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
sign = (i0 & 0x80000000) != 0 ? -1 : 1;
i0 &= 0xfffff;
i0 |= 0x100000;
EXTRACT_WORDS64 (i0, x);
j0 = ((i0 >> 52) & 0x7ff) - 0x3ff;
sign = i0 < 0 ? -1 : 1;
i0 &= UINT64_C(0xfffffffffffff);
i0 |= UINT64_C(0x10000000000000);
if (j0 < 20)
if (j0 < (int32_t) (8 * sizeof (long int)) - 1)
{
if (j0 < 0)
return j0 < -1 ? 0 : sign;
else if (j0 >= 52)
result = i0 << (j0 - 52);
else
{
i0 += 0x80000 >> j0;
i0 += UINT64_C(0x8000000000000) >> j0;
result = i0 >> (20 - j0);
}
}
else if (j0 < (int32_t) (8 * sizeof (long int)) - 1)
{
if (j0 >= 52)
result = ((long int) i0 << (j0 - 20)) | ((long int) i1 << (j0 - 52));
else
{
uint32_t j = i1 + (0x80000000 >> (j0 - 20));
if (j < i1)
++i0;
if (j0 == 20)
result = (long int) i0;
else
{
result = ((long int) i0 << (j0 - 20)) | (j >> (52 - j0));
result = i0 >> (52 - j0);
#ifdef FE_INVALID
if (sizeof (long int) == 4
&& sign == 1
&& result == LONG_MIN)
/* Rounding brought the value out of range. */
feraiseexcept (FE_INVALID);
if (sizeof (long int) == 4
&& sign == 1
&& result == LONG_MIN)
/* Rounding brought the value out of range. */
feraiseexcept (FE_INVALID);
#endif
}
}
}
else
@ -92,8 +77,8 @@ __lround (double x)
return sign == 1 ? LONG_MAX : LONG_MIN;
}
else if (!FIX_DBL_LONG_CONVERT_OVERFLOW
&& sizeof (long int) == 4
&& x <= (double) LONG_MIN - 0.5)
&& sizeof (long int) == 4
&& x <= (double) LONG_MIN - 0.5)
{
/* If truncation produces LONG_MIN, the cast will not raise
the exception, but may raise "inexact". */
@ -108,3 +93,5 @@ __lround (double x)
}
libm_alias_double (__lround, lround)
#endif

80
sysdeps/ieee754/dbl-64/s_modf.c
View File

@ -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)

43
sysdeps/ieee754/dbl-64/s_remquo.c
View File

@ -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)
{

79
sysdeps/ieee754/dbl-64/s_roundeven.c
View File

@ -1,5 +1,4 @@
/* Round to nearest integer value, rounding halfway cases to even.
dbl-64 version.
Copyright (C) 2016-2021 Free Software Foundation, Inc.
This file is part of the GNU C Library.
@ -29,10 +28,10 @@
double
__roundeven (double x)
{
uint32_t hx, lx, uhx;
EXTRACT_WORDS (hx, lx, x);
uhx = hx & 0x7fffffff;
int exponent = uhx >> (MANT_DIG - 1 - 32);
uint64_t ix, ux;
EXTRACT_WORDS64 (ix, x);
ux = ix & 0x7fffffffffffffffULL;
int exponent = ux >> (MANT_DIG - 1);
if (exponent >= BIAS + MANT_DIG - 1)
{
/* Integer, infinity or NaN. */
@ -42,63 +41,29 @@ __roundeven (double x)
else
return x;
}
else if (exponent >= BIAS + MANT_DIG - 32)
{
/* Not necessarily an integer; integer bit is in low word.
Locate the bits with exponents 0 and -1. */
int int_pos = (BIAS + MANT_DIG - 1) - exponent;
int half_pos = int_pos - 1;
uint32_t half_bit = 1U << half_pos;
uint32_t int_bit = 1U << int_pos;
if ((lx & (int_bit | (half_bit - 1))) != 0)
{
/* Carry into the exponent works correctly. No need to test
whether HALF_BIT is set. */
lx += half_bit;
hx += lx < half_bit;
}
lx &= ~(int_bit - 1);
}
else if (exponent == BIAS + MANT_DIG - 33)
{
/* Not necessarily an integer; integer bit is bottom of high
word, half bit is top of low word. */
if (((hx & 1) | (lx & 0x7fffffff)) != 0)
{
lx += 0x80000000;
hx += lx < 0x80000000;
}
lx = 0;
}
else if (exponent >= BIAS)
{
/* At least 1; not necessarily an integer, integer bit and half
bit are in the high word. Locate the bits with exponents 0
and -1 (when the unbiased exponent is 0, the bit with
exponent 0 is implicit, but as the bias is odd it is OK to
take it from the low bit of the exponent). */
int int_pos = (BIAS + MANT_DIG - 33) - exponent;
/* At least 1; not necessarily an integer. Locate the bits with
exponents 0 and -1 (when the unbiased exponent is 0, the bit
with exponent 0 is implicit, but as the bias is odd it is OK
to take it from the low bit of the exponent). */
int int_pos = (BIAS + MANT_DIG - 1) - exponent;
int half_pos = int_pos - 1;
uint32_t half_bit = 1U << half_pos;
uint32_t int_bit = 1U << int_pos;
if (((hx & (int_bit | (half_bit - 1))) | lx) != 0)
hx += half_bit;
hx &= ~(int_bit - 1);
lx = 0;
}
else if (exponent == BIAS - 1 && (uhx > 0x3fe00000 || lx != 0))
{
/* Interval (0.5, 1). */
hx = (hx & 0x80000000) | 0x3ff00000;
lx = 0;
uint64_t half_bit = 1ULL << half_pos;
uint64_t int_bit = 1ULL << int_pos;
if ((ix & (int_bit | (half_bit - 1))) != 0)
/* Carry into the exponent works correctly. No need to test
whether HALF_BIT is set. */
ix += half_bit;
ix &= ~(int_bit - 1);
}
else if (exponent == BIAS - 1 && ux > 0x3fe0000000000000ULL)
/* Interval (0.5, 1). */
ix = (ix & 0x8000000000000000ULL) | 0x3ff0000000000000ULL;
else
{
/* Rounds to 0. */
hx &= 0x80000000;
lx = 0;
}
INSERT_WORDS (x, hx, lx);
/* Rounds to 0. */
ix &= 0x8000000000000000ULL;
INSERT_WORDS64 (x, ix);
return x;
}
hidden_def (__roundeven)

63
sysdeps/ieee754/dbl-64/s_scalbln.c
View File

@ -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)

63
sysdeps/ieee754/dbl-64/s_scalbn.c
View File

@ -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)

42
sysdeps/ieee754/dbl-64/s_setpayload_main.c
View File

@ -1,4 +1,4 @@
/* Set NaN payload. dbl-64 version.
/* Set NaN payload.
Copyright (C) 2016-2021 Free Software Foundation, Inc.
This file is part of the GNU C Library.
@ -30,41 +30,25 @@
int
FUNC (double *x, double payload)
{
uint32_t hx, lx;
EXTRACT_WORDS (hx, lx, payload);
int exponent = hx >> (EXPLICIT_MANT_DIG - 32);
uint64_t ix;
EXTRACT_WORDS64 (ix, payload);
int exponent = ix >> EXPLICIT_MANT_DIG;
/* Test if argument is (a) negative or too large; (b) too small,
except for 0 when allowed; (c) not an integer. */
if (exponent >= BIAS + PAYLOAD_DIG
|| (exponent < BIAS && !(SET_HIGH_BIT && hx == 0 && lx == 0)))
|| (exponent < BIAS && !(SET_HIGH_BIT && ix == 0))
|| (ix & ((1ULL << (BIAS + EXPLICIT_MANT_DIG - exponent)) - 1)) != 0)
{
INSERT_WORDS (*x, 0, 0);
INSERT_WORDS64 (*x, 0);
return 1;
}
int shift = BIAS + EXPLICIT_MANT_DIG - exponent;
if (shift < 32
? (lx & ((1U << shift) - 1)) != 0
: (lx != 0 || (hx & ((1U << (shift - 32)) - 1)) != 0))
if (ix != 0)
{
INSERT_WORDS (*x, 0, 0);
return 1;
ix &= (1ULL << EXPLICIT_MANT_DIG) - 1;
ix |= 1ULL << EXPLICIT_MANT_DIG;
ix >>= BIAS + EXPLICIT_MANT_DIG - exponent;
}
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;
}
}
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;
}

32
sysdeps/ieee754/dbl-64/s_totalorder.c
View File

@ -1,4 +1,4 @@
/* Total order operation. dbl-64 version.
/* Total order operation.
Copyright (C) 2016-2021 Free Software Foundation, Inc.
This file is part of the GNU C Library.
@ -18,8 +18,8 @@
#include <math.h>
#include <math_private.h>
#include <libm-alias-double.h>
#include <nan-high-order-bit.h>
#include <libm-alias-double.h>
#include <stdint.h>
#include <shlib-compat.h>
#include <first-versions.h>
@ -27,30 +27,26 @@
int
__totalorder (const double *x, const double *y)
{
int32_t hx, hy;
uint32_t lx, ly;
EXTRACT_WORDS (hx, lx, *x);
EXTRACT_WORDS (hy, ly, *y);
int64_t ix, iy;
EXTRACT_WORDS64 (ix, *x);
EXTRACT_WORDS64 (iy, *y);
#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN
uint32_t uhx = hx & 0x7fffffff, uhy = hy & 0x7fffffff;
/* For the preferred quiet NaN convention, this operation is a
comparison of the representations of the arguments interpreted as
sign-magnitude integers. If both arguments are NaNs, invert the
quiet/signaling bit so comparing that way works. */
if ((uhx > 0x7ff00000 || (uhx == 0x7ff00000 && lx != 0))
&& (uhy > 0x7ff00000 || (uhy == 0x7ff00000 && ly != 0)))
if ((ix & 0x7fffffffffffffffULL) > 0x7ff0000000000000ULL
&& (iy & 0x7fffffffffffffffULL) > 0x7ff0000000000000ULL)
{
hx ^= 0x00080000;
hy ^= 0x00080000;
ix ^= 0x0008000000000000ULL;
iy ^= 0x0008000000000000ULL;
}
#endif
uint32_t hx_sign = hx >> 31;
uint32_t hy_sign = hy >> 31;
hx ^= hx_sign >> 1;
lx ^= hx_sign;
hy ^= hy_sign >> 1;
ly ^= hy_sign;
return hx < hy || (hx == hy && lx <= ly);
uint64_t ix_sign = ix >> 63;
uint64_t iy_sign = iy >> 63;
ix ^= ix_sign >> 1;
iy ^= iy_sign >> 1;
return ix <= iy;
}
#ifdef SHARED
# define CONCATX(x, y) x ## y

24
sysdeps/ieee754/dbl-64/s_totalordermag.c
View File

@ -1,4 +1,4 @@
/* Total order operation on absolute values. dbl-64 version.
/* Total order operation on absolute values.
Copyright (C) 2016-2021 Free Software Foundation, Inc.
This file is part of the GNU C Library.
@ -18,8 +18,8 @@
#include <math.h>
#include <math_private.h>
#include <libm-alias-double.h>
#include <nan-high-order-bit.h>
#include <libm-alias-double.h>
#include <stdint.h>
#include <shlib-compat.h>
#include <first-versions.h>
@ -27,25 +27,23 @@
int
__totalordermag (const double *x, const double *y)
{
uint32_t hx, hy;
uint32_t lx, ly;
EXTRACT_WORDS (hx, lx, *x);
EXTRACT_WORDS (hy, ly, *y);
hx &= 0x7fffffff;
hy &= 0x7fffffff;
uint64_t ix, iy;
EXTRACT_WORDS64 (ix, *x);
EXTRACT_WORDS64 (iy, *y);
ix &= 0x7fffffffffffffffULL;
iy &= 0x7fffffffffffffffULL;
#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN
/* For the preferred quiet NaN convention, this operation is a
comparison of the representations of the absolute values of the
arguments. If both arguments are NaNs, invert the
quiet/signaling bit so comparing that way works. */
if ((hx > 0x7ff00000 || (hx == 0x7ff00000 && lx != 0))
&& (hy > 0x7ff00000 || (hy == 0x7ff00000 && ly != 0)))
if (ix > 0x7ff0000000000000ULL && iy > 0x7ff0000000000000ULL)
{
hx ^= 0x00080000;
hy ^= 0x00080000;
ix ^= 0x0008000000000000ULL;
iy ^= 0x0008000000000000ULL;
}
#endif
return hx < hy || (hx == hy && lx <= ly);
return ix <= iy;
}
#ifdef SHARED
# define CONCATX(x, y) x ## y

68
sysdeps/ieee754/dbl-64/wordsize-64/e_acosh.c
View File

@ -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)

85
sysdeps/ieee754/dbl-64/wordsize-64/e_cosh.c
View File

@ -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)

106
sysdeps/ieee754/dbl-64/wordsize-64/e_fmod.c
View File

@ -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)

90
sysdeps/ieee754/dbl-64/wordsize-64/e_log10.c
View File

@ -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)

66
sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c
View File

@ -1,66 +0,0 @@
/* Copyright (C) 2011-2021 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, see
<https://www.gnu.org/licenses/>. */
#include <inttypes.h>
#include <math.h>
#include <math_private.h>
#include <libm-alias-double.h>
/*
* for non-zero, finite x
* x = frexp(arg,&exp);
* return a double fp quantity x such that 0.5 <= |x| <1.0
* and the corresponding binary exponent "exp". That is
* arg = x*2^exp.
* If arg is inf, 0.0, or NaN, then frexp(arg,&exp) returns arg
* with *exp=0.
*/
double
__frexp (double x, int *eptr)
{
int64_t ix;
EXTRACT_WORDS64 (ix, x);
int32_t ex = 0x7ff & (ix >> 52);
int e = 0;
if (__glibc_likely (ex != 0x7ff && x != 0.0))
{
/* Not zero and finite. */
e = ex - 1022;
if (__glibc_unlikely (ex == 0))
{
/* Subnormal. */
x *= 0x1p54;
EXTRACT_WORDS64 (ix, x);
ex = 0x7ff & (ix >> 52);
e = ex - 1022 - 54;
}
ix = (ix & INT64_C (0x800fffffffffffff)) | INT64_C (0x3fe0000000000000);
INSERT_WORDS64 (x, ix);
}
else
/* Quiet signaling NaNs. */
x += x;
*eptr = e;
return x;
}
libm_alias_double (__frexp, frexp)

38
sysdeps/ieee754/dbl-64/wordsize-64/s_getpayload.c
View File

@ -1,38 +0,0 @@
/* Get NaN payload. dbl-64/wordsize-64 version.
Copyright (C) 2016-2021 Free Software Foundation, Inc.
This file is part of the GNU C Library.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, see
<https://www.gnu.org/licenses/>. */
#include <math.h>
#include <math_private.h>
#include <libm-alias-double.h>
#include <stdint.h>
#include <fix-int-fp-convert-zero.h>
double
__getpayload (const double *x)
{
uint64_t ix;
EXTRACT_WORDS64 (ix, *x);
if ((ix & 0x7ff0000000000000ULL) != 0x7ff0000000000000ULL
|| (ix & 0xfffffffffffffULL) == 0)
return -1;
ix &= 0x7ffffffffffffULL;
if (FIX_INT_FP_CONVERT_ZERO && ix == 0)
return 0.0f;
return (double) ix;
}
libm_alias_double (__getpayload, getpayload)

43
sysdeps/ieee754/dbl-64/wordsize-64/s_issignaling.c
View File

@ -1,43 +0,0 @@
/* Test for signaling NaN.
Copyright (C) 2013-2021 Free Software Foundation, Inc.
This file is part of the GNU C Library.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, see
<https://www.gnu.org/licenses/>. */
#include <math.h>
#include <math_private.h>
#include <nan-high-order-bit.h>
int
__issignaling (double x)
{
uint64_t xi;
EXTRACT_WORDS64 (xi, x);
#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN
/* We only have to care about the high-order bit of x's significand, because
having it set (sNaN) already makes the significand different from that
used to designate infinity. */
return (xi & UINT64_C (0x7ff8000000000000)) == UINT64_C (0x7ff8000000000000);
#else
/* To keep the following comparison simple, toggle the quiet/signaling bit,
so that it is set for sNaNs. This is inverse to IEEE 754-2008 (as well as
common practice for IEEE 754-1985). */
xi ^= UINT64_C (0x0008000000000000);
/* We have to compare for greater (instead of greater or equal), because x's
significand being all-zero designates infinity not NaN. */
return (xi & UINT64_C (0x7fffffffffffffff)) > UINT64_C (0x7ff8000000000000);
#endif
}
libm_hidden_def (__issignaling)

85
sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c
View File

@ -1,85 +0,0 @@
/* Round double value to long long int.
Copyright (C) 1997-2021 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, see
<https://www.gnu.org/licenses/>. */
#define lround __hidden_lround
#define __lround __hidden___lround
#include <fenv.h>
#include <limits.h>
#include <math.h>
#include <sysdep.h>
#include <math_private.h>
#include <libm-alias-double.h>
#include <fix-fp-int-convert-overflow.h>
long long int
__llround (double x)
{
int32_t j0;
int64_t i0;
long long int result;
int sign;
EXTRACT_WORDS64 (i0, x);
j0 = ((i0 >> 52) & 0x7ff) - 0x3ff;
sign = i0 < 0 ? -1 : 1;
i0 &= UINT64_C(0xfffffffffffff);
i0 |= UINT64_C(0x10000000000000);
if (j0 < (int32_t) (8 * sizeof (long long int)) - 1)
{
if (j0 < 0)
return j0 < -1 ? 0 : sign;
else if (j0 >= 52)
result = i0 << (j0 - 52);
else
{
i0 += UINT64_C(0x8000000000000) >> j0;
result = i0 >> (52 - j0);
}
}
else
{
#ifdef FE_INVALID
/* The number is too large. Unless it rounds to LLONG_MIN,
FE_INVALID must be raised and the return value is
unspecified. */
if (FIX_DBL_LLONG_CONVERT_OVERFLOW && x != (double) LLONG_MIN)
{
feraiseexcept (FE_INVALID);
return sign == 1 ? LLONG_MAX : LLONG_MIN;
}
#endif
return (long long int) x;
}
return sign * result;
}
libm_alias_double (__llround, llround)
/* long has the same width as long long on LP64 machines, so use an alias. */
#undef lround
#undef __lround
#ifdef _LP64
strong_alias (__llround, __lround)
libm_alias_double (__lround, lround)
#endif

97
sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c
View File

@ -1,97 +0,0 @@
/* Round double value to long int.
Copyright (C) 1997-2021 Free Software Foundation, Inc.
This file is part of the GNU C Library.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, see
<https://www.gnu.org/licenses/>. */
#include <fenv.h>
#include <limits.h>
#include <math.h>
#include <math_private.h>
#include <libm-alias-double.h>
#include <fix-fp-int-convert-overflow.h>
/* For LP64, lround is an alias for llround. */
#ifndef _LP64
long int
__lround (double x)
{
int32_t j0;
int64_t i0;
long int result;
int sign;
EXTRACT_WORDS64 (i0, x);
j0 = ((i0 >> 52) & 0x7ff) - 0x3ff;
sign = i0 < 0 ? -1 : 1;
i0 &= UINT64_C(0xfffffffffffff);
i0 |= UINT64_C(0x10000000000000);
if (j0 < (int32_t) (8 * sizeof (long int)) - 1)
{
if (j0 < 0)
return j0 < -1 ? 0 : sign;
else if (j0 >= 52)
result = i0 << (j0 - 52);
else
{
i0 += UINT64_C(0x8000000000000) >> j0;
result = i0 >> (52 - j0);
#ifdef FE_INVALID
if (sizeof (long int) == 4
&& sign == 1
&& result == LONG_MIN)
/* Rounding brought the value out of range. */
feraiseexcept (FE_INVALID);
#endif
}
}
else
{
/* The number is too large. Unless it rounds to LONG_MIN,
FE_INVALID must be raised and the return value is
unspecified. */
#ifdef FE_INVALID
if (FIX_DBL_LONG_CONVERT_OVERFLOW
&& !(sign == -1
&& (sizeof (long int) == 4
? x > (double) LONG_MIN - 0.5
: x >= (double) LONG_MIN)))
{
feraiseexcept (FE_INVALID);
return sign == 1 ? LONG_MAX : LONG_MIN;
}
else if (!FIX_DBL_LONG_CONVERT_OVERFLOW
&& sizeof (long int) == 4
&& x <= (double) LONG_MIN - 0.5)
{
/* If truncation produces LONG_MIN, the cast will not raise
the exception, but may raise "inexact". */
feraiseexcept (FE_INVALID);
return LONG_MIN;
}
#endif
return (long int) x;
}
return sign * result;
}
libm_alias_double (__lround, lround)
#endif

65
sysdeps/ieee754/dbl-64/wordsize-64/s_modf.c
View File

@ -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

111
sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c
View File

@ -1,111 +0,0 @@
/* Compute remainder and a congruent to the quotient.
Copyright (C) 1997-2021 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, see
<https://www.gnu.org/licenses/>. */
#include <math.h>
#include <math_private.h>
#include <libm-alias-double.h>
#include <stdint.h>
static const double zero = 0.0;
double
__remquo (double x, double y, int *quo)
{
int64_t hx, hy;
uint64_t sx, qs;
int cquo;
EXTRACT_WORDS64 (hx, x);
EXTRACT_WORDS64 (hy, y);
sx = hx & UINT64_C(0x8000000000000000);
qs = sx ^ (hy & UINT64_C(0x8000000000000000));
hy &= UINT64_C(0x7fffffffffffffff);
hx &= UINT64_C(0x7fffffffffffffff);
/* Purge off exception values. */
if (__glibc_unlikely (hy == 0))
return (x * y) / (x * y); /* y = 0 */
if (__builtin_expect (hx >= UINT64_C(0x7ff0000000000000) /* x not finite */
|| hy > UINT64_C(0x7ff0000000000000), 0))/* y is NaN */
return (x * y) / (x * y);
if (hy <= UINT64_C(0x7fbfffffffffffff))
x = __ieee754_fmod (x, 8 * y); /* now x < 8y */
if (__glibc_unlikely (hx == hy))
{
*quo = qs ? -1 : 1;
return zero * x;
}
x = fabs (x);
INSERT_WORDS64 (y, hy);
cquo = 0;
if (hy <= UINT64_C(0x7fcfffffffffffff) && x >= 4 * y)
{
x -= 4 * y;
cquo += 4;
}
if (hy <= UINT64_C(0x7fdfffffffffffff) && x >= 2 * y)
{
x -= 2 * y;
cquo += 2;
}
if (hy < UINT64_C(0x0020000000000000))
{
if (x + x > y)
{
x -= y;
++cquo;
if (x + x >= y)
{
x -= y;
++cquo;
}
}
}
else
{
double y_half = 0.5 * y;
if (x > y_half)
{
x -= y;
++cquo;
if (x >= y_half)
{
x -= y;
++cquo;
}
}
}
*quo = qs ? -cquo : cquo;
/* Ensure correct sign of zero result in round-downward mode. */
if (x == 0.0)
x = 0.0;
if (sx)
x = -x;
return x;
}
libm_alias_double (__remquo, remquo)

71
sysdeps/ieee754/dbl-64/wordsize-64/s_roundeven.c
View File

@ -1,71 +0,0 @@
/* Round to nearest integer value, rounding halfway cases to even.
dbl-64/wordsize-64 version.
Copyright (C) 2016-2021 Free Software Foundation, Inc.
This file is part of the GNU C Library.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, see
<https://www.gnu.org/licenses/>. */
#include <math.h>
#include <math_private.h>
#include <libm-alias-double.h>
#include <stdint.h>
#define BIAS 0x3ff
#define MANT_DIG 53
#define MAX_EXP (2 * BIAS + 1)
double
__roundeven (double x)
{
uint64_t ix, ux;
EXTRACT_WORDS64 (ix, x);
ux = ix & 0x7fffffffffffffffULL;
int exponent = ux >> (MANT_DIG - 1);
if (exponent >= BIAS + MANT_DIG - 1)
{
/* Integer, infinity or NaN. */
if (exponent == MAX_EXP)
/* Infinity or NaN; quiet signaling NaNs. */
return x + x;
else
return x;
}
else if (exponent >= BIAS)
{
/* At least 1; not necessarily an integer. Locate the bits with
exponents 0 and -1 (when the unbiased exponent is 0, the bit
with exponent 0 is implicit, but as the bias is odd it is OK
to take it from the low bit of the exponent). */
int int_pos = (BIAS + MANT_DIG - 1) - exponent;
int half_pos = int_pos - 1;
uint64_t half_bit = 1ULL << half_pos;
uint64_t int_bit = 1ULL << int_pos;
if ((ix & (int_bit | (half_bit - 1))) != 0)
/* Carry into the exponent works correctly. No need to test
whether HALF_BIT is set. */
ix += half_bit;
ix &= ~(int_bit - 1);
}
else if (exponent == BIAS - 1 && ux > 0x3fe0000000000000ULL)
/* Interval (0.5, 1). */
ix = (ix & 0x8000000000000000ULL) | 0x3ff0000000000000ULL;
else
/* Rounds to 0. */
ix &= 0x8000000000000000ULL;
INSERT_WORDS64 (x, ix);
return x;
}
hidden_def (__roundeven)
libm_alias_double (__roundeven, roundeven)

60
sysdeps/ieee754/dbl-64/wordsize-64/s_scalbln.c
View File

@ -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

60
sysdeps/ieee754/dbl-64/wordsize-64/s_scalbn.c
View File

@ -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

54
sysdeps/ieee754/dbl-64/wordsize-64/s_setpayload_main.c
View File

@ -1,54 +0,0 @@
/* Set NaN payload. dbl-64/wordsize-64 version.
Copyright (C) 2016-2021 Free Software Foundation, Inc.
This file is part of the GNU C Library.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, see
<https://www.gnu.org/licenses/>. */
#include <math.h>
#include <math_private.h>
#include <libm-alias-double.h>
#include <nan-high-order-bit.h>
#include <stdint.h>
#define SET_HIGH_BIT (HIGH_ORDER_BIT_IS_SET_FOR_SNAN ? SIG : !SIG)
#define BIAS 0x3ff
#define PAYLOAD_DIG 51
#define EXPLICIT_MANT_DIG 52
int
FUNC (double *x, double payload)
{
uint64_t ix;
EXTRACT_WORDS64 (ix, payload);
int exponent = ix >> EXPLICIT_MANT_DIG;
/* Test if argument is (a) negative or too large; (b) too small,
except for 0 when allowed; (c) not an integer. */
if (exponent >= BIAS + PAYLOAD_DIG
|| (exponent < BIAS && !(SET_HIGH_BIT && ix == 0))
|| (ix & ((1ULL << (BIAS + EXPLICIT_MANT_DIG - exponent)) - 1)) != 0)
{
INSERT_WORDS64 (*x, 0);
return 1;
}
if (ix != 0)
{
ix &= (1ULL << EXPLICIT_MANT_DIG) - 1;
ix |= 1ULL << EXPLICIT_MANT_DIG;
ix >>= BIAS + EXPLICIT_MANT_DIG - exponent;
}
ix |= 0x7ff0000000000000ULL | (SET_HIGH_BIT ? 0x8000000000000ULL : 0);
INSERT_WORDS64 (*x, ix);
return 0;
}

76
sysdeps/ieee754/dbl-64/wordsize-64/s_totalorder.c
View File

@ -1,76 +0,0 @@
/* Total order operation. dbl-64/wordsize-64 version.
Copyright (C) 2016-2021 Free Software Foundation, Inc.
This file is part of the GNU C Library.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, see
<https://www.gnu.org/licenses/>. */
#include <math.h>
#include <math_private.h>
#include <nan-high-order-bit.h>
#include <libm-alias-double.h>
#include <stdint.h>
#include <shlib-compat.h>
#include <first-versions.h>
int
__totalorder (const double *x, const double *y)
{
int64_t ix, iy;
EXTRACT_WORDS64 (ix, *x);
EXTRACT_WORDS64 (iy, *y);
#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN
/* For the preferred quiet NaN convention, this operation is a
comparison of the representations of the arguments interpreted as
sign-magnitude integers. If both arguments are NaNs, invert the
quiet/signaling bit so comparing that way works. */
if ((ix & 0x7fffffffffffffffULL) > 0x7ff0000000000000ULL
&& (iy & 0x7fffffffffffffffULL) > 0x7ff0000000000000ULL)
{
ix ^= 0x0008000000000000ULL;
iy ^= 0x0008000000000000ULL;
}
#endif
uint64_t ix_sign = ix >> 63;
uint64_t iy_sign = iy >> 63;
ix ^= ix_sign >> 1;
iy ^= iy_sign >> 1;
return ix <= iy;
}
#ifdef SHARED
# define CONCATX(x, y) x ## y
# define CONCAT(x, y) CONCATX (x, y)
# define UNIQUE_ALIAS(name) CONCAT (name, __COUNTER__)
# define do_symbol(orig_name, name, aliasname) \
strong_alias (orig_name, name) \
versioned_symbol (libm, name, aliasname, GLIBC_2_31)
# undef weak_alias
# define weak_alias(name, aliasname) \
do_symbol (name, UNIQUE_ALIAS (name), aliasname);
#endif
libm_alias_double (__totalorder, totalorder)
#if SHLIB_COMPAT (libm, GLIBC_2_25, GLIBC_2_31)
int
attribute_compat_text_section
__totalorder_compat (double x, double y)
{
return __totalorder (&x, &y);
}
#undef do_symbol
#define do_symbol(orig_name, name, aliasname) \
strong_alias (orig_name, name) \
compat_symbol (libm, name, aliasname, \
CONCAT (FIRST_VERSION_libm_, aliasname))
libm_alias_double (__totalorder_compat, totalorder)
#endif

73
sysdeps/ieee754/dbl-64/wordsize-64/s_totalordermag.c
View File

@ -1,73 +0,0 @@
/* Total order operation on absolute values. dbl-64/wordsize-64 version.
Copyright (C) 2016-2021 Free Software Foundation, Inc.
This file is part of the GNU C Library.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, see
<https://www.gnu.org/licenses/>. */
#include <math.h>
#include <math_private.h>
#include <nan-high-order-bit.h>
#include <libm-alias-double.h>
#include <stdint.h>
#include <shlib-compat.h>
#include <first-versions.h>
int
__totalordermag (const double *x, const double *y)
{
uint64_t ix, iy;
EXTRACT_WORDS64 (ix, *x);
EXTRACT_WORDS64 (iy, *y);
ix &= 0x7fffffffffffffffULL;
iy &= 0x7fffffffffffffffULL;
#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN
/* For the preferred quiet NaN convention, this operation is a
comparison of the representations of the absolute values of the
arguments. If both arguments are NaNs, invert the
quiet/signaling bit so comparing that way works. */
if (ix > 0x7ff0000000000000ULL && iy > 0x7ff0000000000000ULL)
{
ix ^= 0x0008000000000000ULL;
iy ^= 0x0008000000000000ULL;
}
#endif
return ix <= iy;
}
#ifdef SHARED
# define CONCATX(x, y) x ## y
# define CONCAT(x, y) CONCATX (x, y)
# define UNIQUE_ALIAS(name) CONCAT (name, __COUNTER__)
# define do_symbol(orig_name, name, aliasname) \
strong_alias (orig_name, name) \
versioned_symbol (libm, name, aliasname, GLIBC_2_31)
# undef weak_alias
# define weak_alias(name, aliasname) \
do_symbol (name, UNIQUE_ALIAS (name), aliasname);
#endif
libm_alias_double (__totalordermag, totalordermag)
#if SHLIB_COMPAT (libm, GLIBC_2_25, GLIBC_2_31)
int
attribute_compat_text_section
__totalordermag_compat (double x, double y)
{
return __totalordermag (&x, &y);
}
#undef do_symbol
#define do_symbol(orig_name, name, aliasname) \
strong_alias (orig_name, name) \
compat_symbol (libm, name, aliasname, \
CONCAT (FIRST_VERSION_libm_, aliasname))
libm_alias_double (__totalordermag_compat, totalordermag)
#endif

1
sysdeps/mips/mips64/Implies
View File

@ -1,5 +1,4 @@
# MIPS uses IEEE 754 floating point.
mips/ieee754
ieee754/flt-32
ieee754/dbl-64/wordsize-64
ieee754/dbl-64

1
sysdeps/s390/s390-64/Implies
View File

@ -1,2 +1 @@
wordsize-64
ieee754/dbl-64/wordsize-64

1
sysdeps/sparc/sparc64/Implies
View File

@ -1,6 +1,5 @@
wordsize-64
# SPARC uses IEEE 754 floating point.
ieee754/ldbl-128
ieee754/dbl-64/wordsize-64
ieee754/dbl-64
ieee754/flt-32

1
sysdeps/x86_64/Implies
View File

@ -1,6 +1,5 @@
x86
ieee754/float128
ieee754/ldbl-96
ieee754/dbl-64/wordsize-64
ieee754/dbl-64
ieee754/flt-32