diff options
Diffstat (limited to 'sysdeps/ieee754')
-rw-r--r-- | sysdeps/ieee754/dbl-64/s_roundeven.c | 106 | ||||
-rw-r--r-- | sysdeps/ieee754/dbl-64/wordsize-64/s_roundeven.c | 72 | ||||
-rw-r--r-- | sysdeps/ieee754/flt-32/s_roundevenf.c | 68 | ||||
-rw-r--r-- | sysdeps/ieee754/ldbl-128/s_roundevenl.c | 102 | ||||
-rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_roundevenl.c | 69 | ||||
-rw-r--r-- | sysdeps/ieee754/ldbl-96/s_roundevenl.c | 124 | ||||
-rw-r--r-- | sysdeps/ieee754/ldbl-opt/Makefile | 4 | ||||
-rw-r--r-- | sysdeps/ieee754/ldbl-opt/nldbl-roundeven.c | 26 |
8 files changed, 570 insertions, 1 deletions
diff --git a/sysdeps/ieee754/dbl-64/s_roundeven.c b/sysdeps/ieee754/dbl-64/s_roundeven.c new file mode 100644 index 0000000000..27e15e35ad --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_roundeven.c @@ -0,0 +1,106 @@ +/* Round to nearest integer value, rounding halfway cases to even. + dbl-64 version. + Copyright (C) 2016 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 + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <stdint.h> + +#define BIAS 0x3ff +#define MANT_DIG 53 +#define MAX_EXP (2 * BIAS + 1) + +double +roundeven (double x) +{ + uint32_t hx, lx, uhx; + EXTRACT_WORDS (hx, lx, x); + uhx = hx & 0x7fffffff; + int exponent = uhx >> (MANT_DIG - 1 - 32); + 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 + 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; + 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; + } + else + { + /* Rounds to 0. */ + hx &= 0x80000000; + lx = 0; + } + INSERT_WORDS (x, hx, lx); + return x; +} +hidden_def (roundeven) +#ifdef NO_LONG_DOUBLE +weak_alias (roundeven, roundevenl) +#endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_roundeven.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_roundeven.c new file mode 100644 index 0000000000..1fa3ef319d --- /dev/null +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_roundeven.c @@ -0,0 +1,72 @@ +/* Round to nearest integer value, rounding halfway cases to even. + dbl-64/wordsize-64 version. + Copyright (C) 2016 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 + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.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) +#ifdef NO_LONG_DOUBLE +weak_alias (roundeven, roundevenl) +#endif diff --git a/sysdeps/ieee754/flt-32/s_roundevenf.c b/sysdeps/ieee754/flt-32/s_roundevenf.c new file mode 100644 index 0000000000..db2e0487d6 --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_roundevenf.c @@ -0,0 +1,68 @@ +/* Round to nearest integer value, rounding halfway cases to even. + flt-32 version. + Copyright (C) 2016 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 + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <stdint.h> + +#define BIAS 0x7f +#define MANT_DIG 24 +#define MAX_EXP (2 * BIAS + 1) + +float +roundevenf (float x) +{ + uint32_t ix, ux; + GET_FLOAT_WORD (ix, x); + ux = ix & 0x7fffffff; + 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; + uint32_t half_bit = 1U << half_pos; + uint32_t int_bit = 1U << 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 > 0x3f000000) + /* Interval (0.5, 1). */ + ix = (ix & 0x80000000) | 0x3f800000; + else + /* Rounds to 0. */ + ix &= 0x80000000; + SET_FLOAT_WORD (x, ix); + return x; +} diff --git a/sysdeps/ieee754/ldbl-128/s_roundevenl.c b/sysdeps/ieee754/ldbl-128/s_roundevenl.c new file mode 100644 index 0000000000..de042acf2d --- /dev/null +++ b/sysdeps/ieee754/ldbl-128/s_roundevenl.c @@ -0,0 +1,102 @@ +/* Round to nearest integer value, rounding halfway cases to even. + ldbl-128 version. + Copyright (C) 2016 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 + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <stdint.h> + +#define BIAS 0x3fff +#define MANT_DIG 113 +#define MAX_EXP (2 * BIAS + 1) + +_Float128 +roundevenl (_Float128 x) +{ + uint64_t hx, lx, uhx; + GET_LDOUBLE_WORDS64 (hx, lx, x); + uhx = hx & 0x7fffffffffffffffULL; + int exponent = uhx >> (MANT_DIG - 1 - 64); + 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 + MANT_DIG - 64) + { + /* 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; + uint64_t half_bit = 1ULL << half_pos; + uint64_t int_bit = 1ULL << 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 - 65) + { + /* Not necessarily an integer; integer bit is bottom of high + word, half bit is top of low word. */ + if (((hx & 1) | (lx & 0x7fffffffffffffffULL)) != 0) + { + lx += 0x8000000000000000ULL; + hx += lx < 0x8000000000000000ULL; + } + 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 - 65) - exponent; + int half_pos = int_pos - 1; + uint64_t half_bit = 1ULL << half_pos; + uint64_t int_bit = 1ULL << 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 > 0x3ffe000000000000ULL || lx != 0)) + { + /* Interval (0.5, 1). */ + hx = (hx & 0x8000000000000000ULL) | 0x3fff000000000000ULL; + lx = 0; + } + else + { + /* Rounds to 0. */ + hx &= 0x8000000000000000ULL; + lx = 0; + } + SET_LDOUBLE_WORDS64 (x, hx, lx); + return x; +} diff --git a/sysdeps/ieee754/ldbl-128ibm/s_roundevenl.c b/sysdeps/ieee754/ldbl-128ibm/s_roundevenl.c new file mode 100644 index 0000000000..cd91d3905c --- /dev/null +++ b/sysdeps/ieee754/ldbl-128ibm/s_roundevenl.c @@ -0,0 +1,69 @@ +/* Round to nearest integer value, rounding halfway cases to even. + ldbl-128ibm version. + Copyright (C) 2016 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 + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> + +long double +roundevenl (long double x) +{ + double xh, xl, hi; + + ldbl_unpack (x, &xh, &xl); + + if (xh != 0 && isfinite (xh)) + { + hi = roundeven (xh); + if (hi != xh) + { + /* The high part is not an integer; the low part only + affects the result if the high part is exactly half way + between two integers and the low part is nonzero in the + opposite direction to the rounding of the high part. */ + double diff = hi - xh; + if (fabs (diff) == 0.5) + { + if (xl < 0 && diff > 0) + xh = hi - 1; + else if (xl > 0 && diff < 0) + xh = hi + 1; + else + xh = hi; + } + else + xh = hi; + xl = 0; + } + else + { + /* The high part is a nonzero integer. Rounding the low + part to nearest, ties round to even, is always correct, + as a high part that is an odd integer together with a low + part with magnitude 0.5 is not a valid long double. */ + xl = roundeven (xl); + xh = hi; + ldbl_canonicalize_int (&xh, &xl); + } + } + else + /* Quiet signaling NaN arguments. */ + xh += xh; + + return ldbl_pack (xh, xl); +} diff --git a/sysdeps/ieee754/ldbl-96/s_roundevenl.c b/sysdeps/ieee754/ldbl-96/s_roundevenl.c new file mode 100644 index 0000000000..f635b1f683 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_roundevenl.c @@ -0,0 +1,124 @@ +/* Round to nearest integer value, rounding halfway cases to even. + ldbl-96 version. + Copyright (C) 2016 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 + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <stdint.h> + +#define BIAS 0x3fff +#define MANT_DIG 64 +#define MAX_EXP (2 * BIAS + 1) + +long double +roundevenl (long double x) +{ + uint16_t se; + uint32_t hx, lx; + GET_LDOUBLE_WORDS (se, hx, lx, x); + int exponent = se & 0x7fff; + 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 + 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) + { + /* No need to test whether HALF_BIT is set. */ + lx += half_bit; + if (lx < half_bit) + { + hx++; + if (hx == 0) + { + hx = 0x80000000; + se++; + } + } + } + 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; + if (lx < 0x80000000) + { + hx++; + if (hx == 0) + { + hx = 0x80000000; + se++; + } + } + } + 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. */ + int int_pos = (BIAS + MANT_DIG - 33) - 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; + if (hx < half_bit) + { + hx = 0x80000000; + se++; + } + } + hx &= ~(int_bit - 1); + lx = 0; + } + else if (exponent == BIAS - 1 && (hx > 0x80000000 || lx != 0)) + { + /* Interval (0.5, 1). */ + se = (se & 0x8000) | 0x3fff; + hx = 0x80000000; + lx = 0; + } + else + { + /* Rounds to 0. */ + se &= 0x8000; + hx = 0; + lx = 0; + } + SET_LDOUBLE_WORDS (x, se, hx, lx); + return x; +} diff --git a/sysdeps/ieee754/ldbl-opt/Makefile b/sysdeps/ieee754/ldbl-opt/Makefile index 99aa6c6e0b..006d0d59ee 100644 --- a/sysdeps/ieee754/ldbl-opt/Makefile +++ b/sysdeps/ieee754/ldbl-opt/Makefile @@ -43,7 +43,8 @@ libnldbl-calls = asprintf dprintf fprintf fscanf fwprintf fwscanf iovfscanf \ isoc99_wscanf isoc99_fwscanf isoc99_swscanf \ isoc99_vwscanf isoc99_vfwscanf isoc99_vswscanf \ nextup nextdown totalorder totalordermag getpayload \ - canonicalize setpayload setpayloadsig llogb fmaxmag fminmag + canonicalize setpayload setpayloadsig llogb fmaxmag fminmag \ + roundeven libnldbl-routines = $(libnldbl-calls:%=nldbl-%) libnldbl-inhibit-o = $(object-suffixes) libnldbl-static-only-routines = $(libnldbl-routines) @@ -138,6 +139,7 @@ CFLAGS-nldbl-remainder.c = -fno-builtin-remainderl -fno-builtin-dreml CFLAGS-nldbl-remquo.c = -fno-builtin-remquol CFLAGS-nldbl-rint.c = -fno-builtin-rintl CFLAGS-nldbl-round.c = -fno-builtin-roundl +CFLAGS-nldbl-roundeven.c = -fno-builtin-roundevenl CFLAGS-nldbl-scalb.c = -fno-builtin-scalbl CFLAGS-nldbl-scalbln.c = -fno-builtin-scalblnl CFLAGS-nldbl-scalbn.c = -fno-builtin-scalbnl diff --git a/sysdeps/ieee754/ldbl-opt/nldbl-roundeven.c b/sysdeps/ieee754/ldbl-opt/nldbl-roundeven.c new file mode 100644 index 0000000000..cd52e7376f --- /dev/null +++ b/sysdeps/ieee754/ldbl-opt/nldbl-roundeven.c @@ -0,0 +1,26 @@ +/* Compatibility routine for IEEE double as long double for roundeven. + Copyright (C) 2016 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 + <http://www.gnu.org/licenses/>. */ + +#include "nldbl-compat.h" + +double +attribute_hidden +roundevenl (double x) +{ + return roundeven (x); +} |