diff options
Diffstat (limited to 'sysdeps/ieee754/ldbl-128ibm/math_ldbl.h')
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/math_ldbl.h | 290 |
1 files changed, 0 insertions, 290 deletions
diff --git a/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h b/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h deleted file mode 100644 index 8f2984e924..0000000000 --- a/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h +++ /dev/null @@ -1,290 +0,0 @@ -/* Manipulation of the bit representation of 'long double' quantities. - Copyright (C) 2006-2017 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/>. */ - -#ifndef _MATH_LDBL_H_ -#define _MATH_LDBL_H_ 1 - -#include <ieee754.h> -#include <stdint.h> - -/* To suit our callers we return *hi64 and *lo64 as if they came from - an ieee854 112 bit mantissa, that is, 48 bits in *hi64 (plus one - implicit bit) and 64 bits in *lo64. */ - -static inline void -ldbl_extract_mantissa (int64_t *hi64, uint64_t *lo64, int *exp, long double x) -{ - /* We have 105 bits of mantissa plus one implicit digit. Since - 106 bits are representable we use the first implicit digit for - the number before the decimal point and the second implicit bit - as bit 53 of the mantissa. */ - uint64_t hi, lo; - union ibm_extended_long_double u; - - u.ld = x; - *exp = u.d[0].ieee.exponent - IEEE754_DOUBLE_BIAS; - - lo = ((uint64_t) u.d[1].ieee.mantissa0 << 32) | u.d[1].ieee.mantissa1; - hi = ((uint64_t) u.d[0].ieee.mantissa0 << 32) | u.d[0].ieee.mantissa1; - - if (u.d[0].ieee.exponent != 0) - { - int ediff; - - /* If not a denormal or zero then we have an implicit 53rd bit. */ - hi |= (uint64_t) 1 << 52; - - if (u.d[1].ieee.exponent != 0) - lo |= (uint64_t) 1 << 52; - else - /* A denormal is to be interpreted as having a biased exponent - of 1. */ - lo = lo << 1; - - /* We are going to shift 4 bits out of hi later, because we only - want 48 bits in *hi64. That means we want 60 bits in lo, but - we currently only have 53. Shift the value up. */ - lo = lo << 7; - - /* The lower double is normalized separately from the upper. - We may need to adjust the lower mantissa to reflect this. - The difference between the exponents can be larger than 53 - when the low double is much less than 1ULP of the upper - (in which case there are significant bits, all 0's or all - 1's, between the two significands). The difference between - the exponents can be less than 53 when the upper double - exponent is nearing its minimum value (in which case the low - double is denormal ie. has an exponent of zero). */ - ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent - 53; - if (ediff > 0) - { - if (ediff < 64) - lo = lo >> ediff; - else - lo = 0; - } - else if (ediff < 0) - lo = lo << -ediff; - - if (u.d[0].ieee.negative != u.d[1].ieee.negative - && lo != 0) - { - hi--; - lo = ((uint64_t) 1 << 60) - lo; - if (hi < (uint64_t) 1 << 52) - { - /* We have a borrow from the hidden bit, so shift left 1. */ - hi = (hi << 1) | (lo >> 59); - lo = (((uint64_t) 1 << 60) - 1) & (lo << 1); - *exp = *exp - 1; - } - } - } - else - /* If the larger magnitude double is denormal then the smaller - one must be zero. */ - hi = hi << 1; - - *lo64 = (hi << 60) | lo; - *hi64 = hi >> 4; -} - -static inline long double -ldbl_insert_mantissa (int sign, int exp, int64_t hi64, uint64_t lo64) -{ - union ibm_extended_long_double u; - int expnt2; - uint64_t hi, lo; - - u.d[0].ieee.negative = sign; - u.d[1].ieee.negative = sign; - u.d[0].ieee.exponent = exp + IEEE754_DOUBLE_BIAS; - u.d[1].ieee.exponent = 0; - expnt2 = exp - 53 + IEEE754_DOUBLE_BIAS; - - /* Expect 113 bits (112 bits + hidden) right justified in two longs. - The low order 53 bits (52 + hidden) go into the lower double */ - lo = (lo64 >> 7) & (((uint64_t) 1 << 53) - 1); - /* The high order 53 bits (52 + hidden) go into the upper double */ - hi = lo64 >> 60; - hi |= hi64 << 4; - - if (lo != 0) - { - int lzcount; - - /* hidden bit of low double controls rounding of the high double. - If hidden is '1' and either the explicit mantissa is non-zero - or hi is odd, then round up hi and adjust lo (2nd mantissa) - plus change the sign of the low double to compensate. */ - if ((lo & ((uint64_t) 1 << 52)) != 0 - && ((hi & 1) != 0 || (lo & (((uint64_t) 1 << 52) - 1)) != 0)) - { - hi++; - if ((hi & ((uint64_t) 1 << 53)) != 0) - { - hi = hi >> 1; - u.d[0].ieee.exponent++; - } - u.d[1].ieee.negative = !sign; - lo = ((uint64_t) 1 << 53) - lo; - } - - /* Normalize the low double. Shift the mantissa left until - the hidden bit is '1' and adjust the exponent accordingly. */ - - if (sizeof (lo) == sizeof (long)) - lzcount = __builtin_clzl (lo); - else if ((lo >> 32) != 0) - lzcount = __builtin_clzl ((long) (lo >> 32)); - else - lzcount = __builtin_clzl ((long) lo) + 32; - lzcount = lzcount - (64 - 53); - lo <<= lzcount; - expnt2 -= lzcount; - - if (expnt2 >= 1) - /* Not denormal. */ - u.d[1].ieee.exponent = expnt2; - else - { - /* Is denormal. Note that biased exponent of 0 is treated - as if it was 1, hence the extra shift. */ - if (expnt2 > -53) - lo >>= 1 - expnt2; - else - lo = 0; - } - } - else - u.d[1].ieee.negative = 0; - - u.d[1].ieee.mantissa1 = lo; - u.d[1].ieee.mantissa0 = lo >> 32; - u.d[0].ieee.mantissa1 = hi; - u.d[0].ieee.mantissa0 = hi >> 32; - return u.ld; -} - -/* Handy utility functions to pack/unpack/cononicalize and find the nearbyint - of long double implemented as double double. */ -static inline long double -default_ldbl_pack (double a, double aa) -{ - union ibm_extended_long_double u; - u.d[0].d = a; - u.d[1].d = aa; - return u.ld; -} - -static inline void -default_ldbl_unpack (long double l, double *a, double *aa) -{ - union ibm_extended_long_double u; - u.ld = l; - *a = u.d[0].d; - *aa = u.d[1].d; -} - -#ifndef ldbl_pack -# define ldbl_pack default_ldbl_pack -#endif -#ifndef ldbl_unpack -# define ldbl_unpack default_ldbl_unpack -#endif - -/* Extract high double. */ -#define ldbl_high(x) ((double) x) - -/* Convert a finite long double to canonical form. - Does not handle +/-Inf properly. */ -static inline void -ldbl_canonicalize (double *a, double *aa) -{ - double xh, xl; - - xh = *a + *aa; - xl = (*a - xh) + *aa; - *a = xh; - *aa = xl; -} - -/* Simple inline nearbyint (double) function. - Only works in the default rounding mode - but is useful in long double rounding functions. */ -static inline double -ldbl_nearbyint (double a) -{ - double two52 = 0x1p52; - - if (__glibc_likely ((__builtin_fabs (a) < two52))) - { - if (__glibc_likely ((a > 0.0))) - { - a += two52; - a -= two52; - } - else if (__glibc_likely ((a < 0.0))) - { - a = two52 - a; - a = -(a - two52); - } - } - return a; -} - -/* Canonicalize a result from an integer rounding function, in any - rounding mode. *A and *AA are finite and integers, with *A being - nonzero; if the result is not already canonical, *AA is plus or - minus a power of 2 that does not exceed the least set bit in - *A. */ -static inline void -ldbl_canonicalize_int (double *a, double *aa) -{ - /* Previously we used EXTRACT_WORDS64 from math_private.h, but in order - to avoid including internal headers we duplicate that code here. */ - uint64_t ax, aax; - union { double value; uint64_t word; } extractor; - extractor.value = *a; - ax = extractor.word; - extractor.value = *aa; - aax = extractor.word; - - int expdiff = ((ax >> 52) & 0x7ff) - ((aax >> 52) & 0x7ff); - if (expdiff <= 53) - { - if (expdiff == 53) - { - /* Half way between two double values; noncanonical iff the - low bit of A's mantissa is 1. */ - if ((ax & 1) != 0) - { - *a += 2 * *aa; - *aa = -*aa; - } - } - else - { - /* The sum can be represented in a single double. */ - *a += *aa; - *aa = 0; - } - } -} - -#endif /* math_ldbl.h */ |