diff --git a/libquadmath/ChangeLog b/libquadmath/ChangeLog index e169a3d8fd4..51c9ad2be9c 100644 --- a/libquadmath/ChangeLog +++ b/libquadmath/ChangeLog @@ -1,3 +1,17 @@ +2018-11-07 Joseph Myers + + * quadmath-imp.h (ieee854_float128): Use mantissa0, mantissa1, + mantissa2 and mantissa3 fields instead of mant_high and mant_low. + Change nan field to ieee_nan. + * update-quadmath.py (update_sources): Also update fmaq.c. + * math/nanq.c (nanq): Use ieee_nan field of union. + Zero-initialize f. Set quiet_nan field. + * printf/flt1282mpn.c, printf/printf_fphex.c, strtod/mpn2flt128.c, + strtod/strtoflt128.c: Use mantissa0, mantissa1, mantissa2 and + mantissa3 fields. Use ieee_nan and quiet_nan field. + * math/fmaq.c: Regenerate from glibc sources with + update-quadmath.py. + 2018-11-05 Joseph Myers PR libquadmath/68686 diff --git a/libquadmath/math/fmaq.c b/libquadmath/math/fmaq.c index 68a63cf8fd6..5fe4c39302c 100644 --- a/libquadmath/math/fmaq.c +++ b/libquadmath/math/fmaq.c @@ -1,5 +1,5 @@ /* Compute x * y + z as ternary operation. - Copyright (C) 2010-2017 Free Software Foundation, Inc. + Copyright (C) 2010-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek , 2010. @@ -18,16 +18,6 @@ . */ #include "quadmath-imp.h" -#include -#include -#ifdef HAVE_FENV_H -# include -# if defined HAVE_FEHOLDEXCEPT && defined HAVE_FESETROUND \ - && defined HAVE_FEUPDATEENV && defined HAVE_FETESTEXCEPT \ - && defined FE_TOWARDZERO && defined FE_INEXACT -# define USE_FENV_H -# endif -#endif /* This implementation uses rounding to odd to avoid problems with double rounding. See a paper by Boldo and Melquiond: @@ -73,7 +63,7 @@ fmaq (__float128 x, __float128 y, __float128 z) if (u.ieee.exponent + v.ieee.exponent > 0x7fff + IEEE854_FLOAT128_BIAS) return x * y; - /* If x * y is less than 1/4 of FLT128_DENORM_MIN, neither the + /* If x * y is less than 1/4 of FLT128_TRUE_MIN, neither the result nor whether there is underflow depends on its exact value, only on its sign. */ if (u.ieee.exponent + v.ieee.exponent @@ -94,8 +84,10 @@ fmaq (__float128 x, __float128 y, __float128 z) : (w.ieee.exponent == 0 || (w.ieee.exponent == 1 && w.ieee.negative != neg - && w.ieee.mant_low == 0 - && w.ieee.mant_high == 0))) + && w.ieee.mantissa3 == 0 + && w.ieee.mantissa2 == 0 + && w.ieee.mantissa1 == 0 + && w.ieee.mantissa0 == 0))) { __float128 force_underflow = x * y; math_force_eval (force_underflow); @@ -124,7 +116,7 @@ fmaq (__float128 x, __float128 y, __float128 z) very small, adjust them up to avoid spurious underflows, rather than down. */ if (u.ieee.exponent + v.ieee.exponent - <= IEEE854_FLOAT128_BIAS + FLT128_MANT_DIG) + <= IEEE854_FLOAT128_BIAS + 2 * FLT128_MANT_DIG) { if (u.ieee.exponent > v.ieee.exponent) u.ieee.exponent += 2 * FLT128_MANT_DIG + 2; @@ -181,17 +173,15 @@ fmaq (__float128 x, __float128 y, __float128 z) } /* Ensure correct sign of exact 0 + 0. */ - if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0)) + if (__glibc_unlikely ((x == 0 || y == 0) && z == 0)) { x = math_opt_barrier (x); return x * y + z; } -#ifdef USE_FENV_H fenv_t env; feholdexcept (&env); fesetround (FE_TONEAREST); -#endif /* Multiplication m1 + m2 = x * y using Dekker's algorithm. */ #define C ((1LL << (FLT128_MANT_DIG + 1) / 2) + 1) @@ -214,62 +204,46 @@ fmaq (__float128 x, __float128 y, __float128 z) /* Ensure the arithmetic is not scheduled after feclearexcept call. */ math_force_eval (m2); math_force_eval (a2); -#ifdef USE_FENV_H feclearexcept (FE_INEXACT); -#endif /* If the result is an exact zero, ensure it has the correct sign. */ if (a1 == 0 && m2 == 0) { -#ifdef USE_FENV_H feupdateenv (&env); -#endif /* Ensure that round-to-nearest value of z + m1 is not reused. */ z = math_opt_barrier (z); return z + m1; } -#ifdef USE_FENV_H fesetround (FE_TOWARDZERO); -#endif /* Perform m2 + a2 addition with round to odd. */ u.value = a2 + m2; - if (__builtin_expect (adjust == 0, 1)) + if (__glibc_likely (adjust == 0)) { -#ifdef USE_FENV_H - if ((u.ieee.mant_low & 1) == 0 && u.ieee.exponent != 0x7fff) - u.ieee.mant_low |= fetestexcept (FE_INEXACT) != 0; + if ((u.ieee.mantissa3 & 1) == 0 && u.ieee.exponent != 0x7fff) + u.ieee.mantissa3 |= fetestexcept (FE_INEXACT) != 0; feupdateenv (&env); -#endif /* Result is a1 + u.value. */ return a1 + u.value; } - else if (__builtin_expect (adjust > 0, 1)) + else if (__glibc_likely (adjust > 0)) { -#ifdef USE_FENV_H - if ((u.ieee.mant_low & 1) == 0 && u.ieee.exponent != 0x7fff) - u.ieee.mant_low |= fetestexcept (FE_INEXACT) != 0; + if ((u.ieee.mantissa3 & 1) == 0 && u.ieee.exponent != 0x7fff) + u.ieee.mantissa3 |= fetestexcept (FE_INEXACT) != 0; feupdateenv (&env); -#endif /* Result is a1 + u.value, scaled up. */ return (a1 + u.value) * 0x1p113Q; } else { -#ifdef USE_FENV_H - if ((u.ieee.mant_low & 1) == 0) - u.ieee.mant_low |= fetestexcept (FE_INEXACT) != 0; -#endif + if ((u.ieee.mantissa3 & 1) == 0) + u.ieee.mantissa3 |= fetestexcept (FE_INEXACT) != 0; v.value = a1 + u.value; /* Ensure the addition is not scheduled after fetestexcept call. */ - asm volatile ("" : : "m" (v.value)); -#ifdef USE_FENV_H + math_force_eval (v.value); int j = fetestexcept (FE_INEXACT) != 0; feupdateenv (&env); -#else - int j = 0; -#endif /* Ensure the following computations are performed in default rounding mode instead of just reusing the round to zero computation. */ asm volatile ("" : "=m" (u) : "m" (u)); @@ -281,11 +255,11 @@ fmaq (__float128 x, __float128 y, __float128 z) rounding will occur. */ if (v.ieee.exponent > 228) return (a1 + u.value) * 0x1p-228Q; - /* If v.value * 0x1p-228Q with round to zero is a subnormal above - or equal to FLT128_MIN / 2, then v.value * 0x1p-228Q shifts mantissa - down just by 1 bit, which means v.ieee.mant_low |= j would + /* If v.value * 0x1p-228L with round to zero is a subnormal above + or equal to FLT128_MIN / 2, then v.value * 0x1p-228L shifts mantissa + down just by 1 bit, which means v.ieee.mantissa3 |= j would change the round bit, not sticky or guard bit. - v.value * 0x1p-228Q never normalizes by shifting up, + v.value * 0x1p-228L never normalizes by shifting up, so round bit plus sticky bit should be already enough for proper rounding. */ if (v.ieee.exponent == 228) @@ -301,18 +275,18 @@ fmaq (__float128 x, __float128 y, __float128 z) if (w.ieee.exponent == 229) return w.value * 0x1p-228Q; } - /* v.ieee.mant_low & 2 is LSB bit of the result before rounding, - v.ieee.mant_low & 1 is the round bit and j is our sticky - bit. */ - w.value = 0.0Q; - w.ieee.mant_low = ((v.ieee.mant_low & 3) << 1) | j; + /* v.ieee.mantissa3 & 2 is LSB bit of the result before rounding, + v.ieee.mantissa3 & 1 is the round bit and j is our sticky + bit. */ + w.value = 0; + w.ieee.mantissa3 = ((v.ieee.mantissa3 & 3) << 1) | j; w.ieee.negative = v.ieee.negative; - v.ieee.mant_low &= ~3U; + v.ieee.mantissa3 &= ~3U; v.value *= 0x1p-228Q; w.value *= 0x1p-2Q; return v.value + w.value; } - v.ieee.mant_low |= j; + v.ieee.mantissa3 |= j; return v.value * 0x1p-228Q; } } diff --git a/libquadmath/math/nanq.c b/libquadmath/math/nanq.c index bace4706459..52786d906c5 100644 --- a/libquadmath/math/nanq.c +++ b/libquadmath/math/nanq.c @@ -4,8 +4,8 @@ __float128 nanq (const char *tagp __attribute__ ((unused))) { // FIXME -- we should use the argument - ieee854_float128 f; - f.ieee.exponent = 0x7fff; - f.ieee.mant_high = 0x1; + ieee854_float128 f = { 0 }; + f.ieee_nan.exponent = 0x7fff; + f.ieee_nan.quiet_nan = 0x1; return f.value; } diff --git a/libquadmath/printf/flt1282mpn.c b/libquadmath/printf/flt1282mpn.c index 0105314ef3a..a9a4c4fbf35 100644 --- a/libquadmath/printf/flt1282mpn.c +++ b/libquadmath/printf/flt1282mpn.c @@ -39,14 +39,14 @@ mpn_extract_flt128 (mp_ptr res_ptr, mp_size_t size, *expt = (int) u.ieee.exponent - IEEE854_FLOAT128_BIAS; #if BITS_PER_MP_LIMB == 32 - res_ptr[0] = u.ieee.mant_low; /* Low-order 32 bits of fraction. */ - res_ptr[1] = (u.ieee.mant_low >> 32); - res_ptr[2] = u.ieee.mant_high; - res_ptr[3] = (u.ieee.mant_high >> 32); /* High-order 32 bits. */ + res_ptr[0] = u.ieee.mantissa3; /* Low-order 32 bits of fraction. */ + res_ptr[1] = u.ieee.mantissa2; + res_ptr[2] = u.ieee.mantissa1; + res_ptr[3] = u.ieee.mantissa0; /* High-order 32 bits. */ #define N 4 #elif BITS_PER_MP_LIMB == 64 - res_ptr[0] = u.ieee.mant_low; - res_ptr[1] = u.ieee.mant_high; + res_ptr[0] = ((mp_limb_t) u.ieee.mantissa2 << 32) | u.ieee.mantissa3; + res_ptr[1] = ((mp_limb_t) u.ieee.mantissa0 << 32) | u.ieee.mantissa1; #define N 2 #else #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for" diff --git a/libquadmath/printf/printf_fphex.c b/libquadmath/printf/printf_fphex.c index fc960f38eb9..a40a6b00945 100644 --- a/libquadmath/printf/printf_fphex.c +++ b/libquadmath/printf/printf_fphex.c @@ -235,8 +235,10 @@ __quadmath_printf_fphex (struct __quadmath_printf_file *fp, assert (sizeof (long double) == 16); - num0 = fpnum.ieee.mant_high; - num1 = fpnum.ieee.mant_low; + num0 = (((unsigned long long int) fpnum.ieee.mantissa0) << 32 + | fpnum.ieee.mantissa1); + num1 = (((unsigned long long int) fpnum.ieee.mantissa2) << 32 + | fpnum.ieee.mantissa3); zero_mantissa = (num0|num1) == 0; diff --git a/libquadmath/quadmath-imp.h b/libquadmath/quadmath-imp.h index 8d22248504a..86b57878efc 100644 --- a/libquadmath/quadmath-imp.h +++ b/libquadmath/quadmath-imp.h @@ -96,11 +96,15 @@ typedef union #if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__ unsigned negative:1; unsigned exponent:15; - uint64_t mant_high:48; - uint64_t mant_low:64; + unsigned mantissa0:16; + unsigned mantissa1:32; + unsigned mantissa2:32; + unsigned mantissa3:32; #else - uint64_t mant_low:64; - uint64_t mant_high:48; + unsigned mantissa3:32; + unsigned mantissa2:32; + unsigned mantissa1:32; + unsigned mantissa0:16; unsigned exponent:15; unsigned negative:1; #endif @@ -142,16 +146,20 @@ typedef union unsigned negative:1; unsigned exponent:15; unsigned quiet_nan:1; - uint64_t mant_high:47; - uint64_t mant_low:64; + unsigned mantissa0:15; + unsigned mantissa1:32; + unsigned mantissa2:32; + unsigned mantissa3:32; #else - uint64_t mant_low:64; - uint64_t mant_high:47; + unsigned mantissa3:32; + unsigned mantissa2:32; + unsigned mantissa1:32; + unsigned mantissa0:15; unsigned quiet_nan:1; unsigned exponent:15; unsigned negative:1; #endif - } nan; + } ieee_nan; } ieee854_float128; diff --git a/libquadmath/strtod/mpn2flt128.c b/libquadmath/strtod/mpn2flt128.c index 844ae97d834..33cdb6232cd 100644 --- a/libquadmath/strtod/mpn2flt128.c +++ b/libquadmath/strtod/mpn2flt128.c @@ -34,15 +34,17 @@ mpn_construct_float128 (mp_srcptr frac_ptr, int expt, int sign) u.ieee.negative = sign; u.ieee.exponent = expt + IEEE854_FLOAT128_BIAS; #if BITS_PER_MP_LIMB == 32 - u.ieee.mant_low = (((uint64_t) frac_ptr[1]) << 32) - | (frac_ptr[0] & 0xffffffff); - u.ieee.mant_high = (((uint64_t) frac_ptr[3] - & (((mp_limb_t) 1 << (FLT128_MANT_DIG - 96)) - 1)) - << 32) | (frac_ptr[2] & 0xffffffff); + u.ieee.mantissa3 = frac_ptr[0]; + u.ieee.mantissa2 = frac_ptr[1]; + u.ieee.mantissa1 = frac_ptr[2]; + u.ieee.mantissa0 = frac_ptr[3] & (((mp_limb_t) 1 + << (FLT128_MANT_DIG - 96)) - 1); #elif BITS_PER_MP_LIMB == 64 - u.ieee.mant_low = frac_ptr[0]; - u.ieee.mant_high = frac_ptr[1] - & (((mp_limb_t) 1 << (FLT128_MANT_DIG - 64)) - 1); + u.ieee.mantissa3 = frac_ptr[0] & (((mp_limb_t) 1 << 32) - 1); + u.ieee.mantissa2 = frac_ptr[0] >> 32; + u.ieee.mantissa1 = frac_ptr[1] & (((mp_limb_t) 1 << 32) - 1); + u.ieee.mantissa0 = (frac_ptr[1] >> 32) & (((mp_limb_t) 1 + << (FLT128_MANT_DIG - 96)) - 1); #else #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for" #endif diff --git a/libquadmath/strtod/strtoflt128.c b/libquadmath/strtod/strtoflt128.c index acdf36e9e8a..cf2da4f1a5c 100644 --- a/libquadmath/strtod/strtoflt128.c +++ b/libquadmath/strtod/strtoflt128.c @@ -30,12 +30,15 @@ #endif #define MPN2FLOAT mpn_construct_float128 #define FLOAT_HUGE_VAL HUGE_VALQ -#define SET_MANTISSA(flt, mant) \ - do { ieee854_float128 u; \ - u.value = (flt); \ - u.ieee.mant_high = 0x800000000000ULL; \ - u.ieee.mant_low = mant; \ - (flt) = u.value; \ +#define SET_MANTISSA(flt, mant) \ + do { ieee854_float128 u; \ + u.value = (flt); \ + u.ieee_nan.mantissa0 = 0; \ + u.ieee_nan.mantissa1 = 0; \ + u.ieee_nan.mantissa2 = (mant) >> 32; \ + u.ieee_nan.mantissa3 = (mant); \ + u.ieee_nan.quiet_nan = 1; \ + (flt) = u.value; \ } while (0) static inline __attribute__((__always_inline__)) diff --git a/libquadmath/update-quadmath.py b/libquadmath/update-quadmath.py index ca6c9f0c7c7..d40b2724dd3 100755 --- a/libquadmath/update-quadmath.py +++ b/libquadmath/update-quadmath.py @@ -143,8 +143,7 @@ def update_sources(glibc_srcdir, quadmath_srcdir): # Replace all #includes with a single include of quadmath-imp.h. repl_map['(\n+#include[^\n]*)+\n+'] = '\n\n#include "quadmath-imp.h"\n\n' # Omitted from this list because code comes from more than one - # glibc source file: rem_pio2. Omitted because of a union not - # currently provided in libquadmath: fma. + # glibc source file: rem_pio2. ldbl_files = { 'e_acoshl.c': 'acoshq.c', 'e_acosl.c': 'acosq.c', 's_asinhl.c': 'asinhq.c', 'e_asinl.c': 'asinq.c', @@ -155,7 +154,7 @@ def update_sources(glibc_srcdir, quadmath_srcdir): 's_erfl.c': 'erfq.c', 's_expm1l.c': 'expm1q.c', 'e_expl.c': 'expq.c', 't_expl.h': 'expq_table.h', 's_fabsl.c': 'fabsq.c', 's_finitel.c': 'finiteq.c', 's_floorl.c': 'floorq.c', - 'e_fmodl.c': 'fmodq.c', 's_frexpl.c': 'frexpq.c', + 's_fmal.c': 'fmaq.c', 'e_fmodl.c': 'fmodq.c', 's_frexpl.c': 'frexpq.c', 'e_lgammal_r.c': 'lgammaq.c', 'lgamma_negl.c': 'lgammaq_neg.c', 'lgamma_productl.c': 'lgammaq_product.c', 'e_hypotl.c': 'hypotq.c', 'e_ilogbl.c': 'ilogbq.c', 's_isinfl.c': 'isinfq.c',