===================================================================
@@ -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 <jakub@redhat.com>, 2010.
@@ -18,16 +18,6 @@
<http://www.gnu.org/licenses/>. */
#include "quadmath-imp.h"
-#include <math.h>
-#include <float.h>
-#ifdef HAVE_FENV_H
-# include <fenv.h>
-# 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 @@
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 @@
: (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 @@
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 @@
}
/* 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 @@
/* 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 @@
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 @@
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;
}
}
===================================================================
@@ -4,8 +4,8 @@
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;
}
===================================================================
@@ -39,14 +39,14 @@
*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"
===================================================================
@@ -235,8 +235,10 @@
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;
===================================================================
@@ -96,11 +96,15 @@
#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 @@
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;
===================================================================
@@ -34,15 +34,17 @@
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
===================================================================
@@ -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__))
===================================================================
@@ -143,8 +143,7 @@
# 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 @@
'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',