/* 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.
<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:
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
: (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);
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;
}
/* 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)
/* 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));
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)
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;
}
}
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
# 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',
'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',
projects / gcc.git / commitdiff