/* This is a software floating point library which can be used
for targets without hardware floating point.
- Copyright (C) 1994, 1995, 1996, 1997, 1998, 2000, 2001, 2002, 2003
- Free Software Foundation, Inc.
+ Copyright (C) 1994, 1995, 1996, 1997, 1998, 2000, 2001, 2002, 2003,
+ 2004, 2005 Free Software Foundation, Inc.
-This file is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 2, or (at your option) any
-later version.
+This file is part of GCC.
+
+GCC is free software; you can redistribute it and/or modify it under
+the terms of the GNU General Public License as published by the Free
+Software Foundation; either version 2, or (at your option) any later
+version.
In addition to the permissions in the GNU General Public License, the
Free Software Foundation gives you unlimited permission to link the
-compiled version of this file with other programs, and to distribute
-those programs without any restriction coming from the use of this
-file. (The General Public License restrictions do apply in other
-respects; for example, they cover modification of the file, and
-distribution when not linked into another program.)
-
-This file 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
-General Public License for more details.
+compiled version of this file into combinations with other programs,
+and to distribute those combinations without any restriction coming
+from the use of this file. (The General Public License restrictions
+do apply in other respects; for example, they cover modification of
+the file, and distribution when not linked into a combine
+executable.)
+
+GCC 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 General Public License
+for more details.
You should have received a copy of the GNU General Public License
-along with this program; see the file COPYING. If not, write to
-the Free Software Foundation, 59 Temple Place - Suite 330,
-Boston, MA 02111-1307, USA. */
-
-/* As a special exception, if you link this library with other files,
- some of which are compiled with GCC, to produce an executable,
- this library does not by itself cause the resulting executable
- to be covered by the GNU General Public License.
- This exception does not however invalidate any other reasons why
- the executable file might be covered by the GNU General Public License. */
+along with GCC; see the file COPYING. If not, write to the Free
+Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
+02110-1301, USA. */
/* This implements IEEE 754 format arithmetic, but does not provide a
mechanism for setting the rounding mode, or for generating or handling
#include "tconfig.h"
#include "coretypes.h"
#include "tm.h"
-#include "fp-bit.h"
+#include "config/fp-bit.h"
/* The following macros can be defined to change the behavior of this file:
FLOAT: Implement a `float', aka SFmode, fp library. If this is not
static int
isnan ( fp_number_type * x)
{
- return x->class == CLASS_SNAN || x->class == CLASS_QNAN;
+ return __builtin_expect (x->class == CLASS_SNAN || x->class == CLASS_QNAN,
+ 0);
}
INLINE
static int
isinf ( fp_number_type * x)
{
- return x->class == CLASS_INFINITY;
+ return __builtin_expect (x->class == CLASS_INFINITY, 0);
}
#endif /* NO_NANS */
x->sign = !x->sign;
}
+/* Count leading zeroes in N. */
+INLINE
+static int
+clzusi (USItype n)
+{
+ extern int __clzsi2 (USItype);
+ if (sizeof (USItype) == sizeof (unsigned int))
+ return __builtin_clz (n);
+ else if (sizeof (USItype) == sizeof (unsigned long))
+ return __builtin_clzl (n);
+ else if (sizeof (USItype) == sizeof (unsigned long long))
+ return __builtin_clzll (n);
+ else
+ return __clzsi2 (n);
+}
+
extern FLO_type pack_d ( fp_number_type * );
#if defined(L_pack_df) || defined(L_pack_sf) || defined(L_pack_tf)
exp = EXPMAX;
if (src->class == CLASS_QNAN || 1)
{
+#ifdef QUIET_NAN_NEGATED
+ fraction |= QUIET_NAN - 1;
+#else
fraction |= QUIET_NAN;
+#endif
}
}
else if (isinf (src))
}
else
{
- if (src->normal_exp < NORMAL_EXPMIN)
+ if (__builtin_expect (src->normal_exp < NORMAL_EXPMIN, 0))
{
#ifdef NO_DENORMALS
/* Go straight to a zero representation if denormals are not
#endif /* NO_DENORMALS */
}
else if (!LARGEST_EXPONENT_IS_NORMAL (FRAC_NBITS)
- && src->normal_exp > EXPBIAS)
+ && __builtin_expect (src->normal_exp > EXPBIAS, 0))
{
exp = EXPMAX;
fraction = 0;
dst.bits.exp = exp;
dst.bits.sign = sign;
#else
+# if defined TFLOAT && defined HALFFRACBITS
+ {
+ halffractype high, low, unity;
+ int lowsign, lowexp;
+
+ unity = (halffractype) 1 << HALFFRACBITS;
+
+ /* Set HIGH to the high double's significand, masking out the implicit 1.
+ Set LOW to the low double's full significand. */
+ high = (fraction >> (FRACBITS - HALFFRACBITS)) & (unity - 1);
+ low = fraction & (unity * 2 - 1);
+
+ /* Get the initial sign and exponent of the low double. */
+ lowexp = exp - HALFFRACBITS - 1;
+ lowsign = sign;
+
+ /* HIGH should be rounded like a normal double, making |LOW| <=
+ 0.5 ULP of HIGH. Assume round-to-nearest. */
+ if (exp < EXPMAX)
+ if (low > unity || (low == unity && (high & 1) == 1))
+ {
+ /* Round HIGH up and adjust LOW to match. */
+ high++;
+ if (high == unity)
+ {
+ /* May make it infinite, but that's OK. */
+ high = 0;
+ exp++;
+ }
+ low = unity * 2 - low;
+ lowsign ^= 1;
+ }
+
+ high |= (halffractype) exp << HALFFRACBITS;
+ high |= (halffractype) sign << (HALFFRACBITS + EXPBITS);
+
+ if (exp == EXPMAX || exp == 0 || low == 0)
+ low = 0;
+ else
+ {
+ while (lowexp > 0 && low < unity)
+ {
+ low <<= 1;
+ lowexp--;
+ }
+
+ if (lowexp <= 0)
+ {
+ halffractype roundmsb, round;
+ int shift;
+
+ shift = 1 - lowexp;
+ roundmsb = (1 << (shift - 1));
+ round = low & ((roundmsb << 1) - 1);
+
+ low >>= shift;
+ lowexp = 0;
+
+ if (round > roundmsb || (round == roundmsb && (low & 1) == 1))
+ {
+ low++;
+ if (low == unity)
+ /* LOW rounds up to the smallest normal number. */
+ lowexp++;
+ }
+ }
+
+ low &= unity - 1;
+ low |= (halffractype) lowexp << HALFFRACBITS;
+ low |= (halffractype) lowsign << (HALFFRACBITS + EXPBITS);
+ }
+ dst.value_raw = ((fractype) high << HALFSHIFT) | low;
+ }
+# else
dst.value_raw = fraction & ((((fractype)1) << FRACBITS) - (fractype)1);
dst.value_raw |= ((fractype) (exp & ((1 << EXPBITS) - 1))) << FRACBITS;
dst.value_raw |= ((fractype) (sign & 1)) << (FRACBITS | EXPBITS);
+# endif
#endif
#if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT)
exp = src->bits.exp;
sign = src->bits.sign;
#else
+# if defined TFLOAT && defined HALFFRACBITS
+ {
+ halffractype high, low;
+
+ high = src->value_raw >> HALFSHIFT;
+ low = src->value_raw & (((fractype)1 << HALFSHIFT) - 1);
+
+ fraction = high & ((((fractype)1) << HALFFRACBITS) - 1);
+ fraction <<= FRACBITS - HALFFRACBITS;
+ exp = ((int)(high >> HALFFRACBITS)) & ((1 << EXPBITS) - 1);
+ sign = ((int)(high >> (((HALFFRACBITS + EXPBITS))))) & 1;
+
+ if (exp != EXPMAX && exp != 0 && low != 0)
+ {
+ int lowexp = ((int)(low >> HALFFRACBITS)) & ((1 << EXPBITS) - 1);
+ int lowsign = ((int)(low >> (((HALFFRACBITS + EXPBITS))))) & 1;
+ int shift;
+ fractype xlow;
+
+ xlow = low & ((((fractype)1) << HALFFRACBITS) - 1);
+ if (lowexp)
+ xlow |= (((halffractype)1) << HALFFRACBITS);
+ else
+ lowexp = 1;
+ shift = (FRACBITS - HALFFRACBITS) - (exp - lowexp);
+ if (shift > 0)
+ xlow <<= shift;
+ else if (shift < 0)
+ xlow >>= -shift;
+ if (sign == lowsign)
+ fraction += xlow;
+ else if (fraction >= xlow)
+ fraction -= xlow;
+ else
+ {
+ /* The high part is a power of two but the full number is lower.
+ This code will leave the implicit 1 in FRACTION, but we'd
+ have added that below anyway. */
+ fraction = (((fractype) 1 << FRACBITS) - xlow) << 1;
+ exp--;
+ }
+ }
+ }
+# else
fraction = src->value_raw & ((((fractype)1) << FRACBITS) - 1);
exp = ((int)(src->value_raw >> FRACBITS)) & ((1 << EXPBITS) - 1);
sign = ((int)(src->value_raw >> (FRACBITS + EXPBITS))) & 1;
+# endif
#endif
dst->sign = sign;
dst->fraction.ll = fraction;
}
}
- else if (!LARGEST_EXPONENT_IS_NORMAL (FRAC_NBITS) && exp == EXPMAX)
+ else if (!LARGEST_EXPONENT_IS_NORMAL (FRAC_NBITS)
+ && __builtin_expect (exp == EXPMAX, 0))
{
/* Huge exponent*/
if (fraction == 0)
else
{
/* Nonzero fraction, means nan */
+#ifdef QUIET_NAN_NEGATED
+ if ((fraction & QUIET_NAN) == 0)
+#else
if (fraction & QUIET_NAN)
+#endif
{
dst->class = CLASS_QNAN;
}
they're the same */
{
int diff;
+ int sdiff;
a_normal_exp = a->normal_exp;
b_normal_exp = b->normal_exp;
b_fraction = b->fraction.ll;
diff = a_normal_exp - b_normal_exp;
+ sdiff = diff;
if (diff < 0)
diff = -diff;
if (diff < FRAC_NBITS)
{
- /* ??? This does shifts one bit at a time. Optimize. */
- while (a_normal_exp > b_normal_exp)
+ if (sdiff > 0)
{
- b_normal_exp++;
- LSHIFT (b_fraction);
+ b_normal_exp += diff;
+ LSHIFT (b_fraction, diff);
}
- while (b_normal_exp > a_normal_exp)
+ else if (sdiff < 0)
{
- a_normal_exp++;
- LSHIFT (a_fraction);
+ a_normal_exp += diff;
+ LSHIFT (a_fraction, diff);
}
}
else
if (tmp->fraction.ll >= IMPLICIT_2)
{
- LSHIFT (tmp->fraction.ll);
+ LSHIFT (tmp->fraction.ll, 1);
tmp->normal_exp++;
}
return tmp;
high |= 1;
low <<= 1;
}
- /* rounding is tricky. if we only round if it won't make us round later. */
-#if 0
- if (low & FRACHIGH2)
- {
- if (((high & GARDMASK) != GARDMSB)
- && (((high + 1) & GARDMASK) == GARDMSB))
- {
- /* don't round, it gets done again later. */
- }
- else
- {
- high++;
- }
- }
-#endif
+
if (!ROUND_TOWARDS_ZERO && (high & GARDMASK) == GARDMSB)
{
if (high & (1 << NGARDS))
{
- /* half way, so round to even */
- high += GARDROUND + 1;
+ /* Because we're half way, we would round to even by adding
+ GARDROUND + 1, except that's also done in the packing
+ function, and rounding twice will lose precision and cause
+ the result to be too far off. Example: 32-bit floats with
+ bit patterns 0xfff * 0x3f800400 ~= 0xfff (less than 0.5ulp
+ off), not 0x1000 (more than 0.5ulp off). */
}
else if (low)
{
- /* but we really weren't half way */
+ /* We're a further than half way by a small amount corresponding
+ to the bits set in "low". Knowing that, we round here and
+ not in pack_d, because there we don't have "low" available
+ anymore. */
high += GARDROUND + 1;
+
+ /* Avoid further rounding in pack_d. */
+ high &= ~(fractype) GARDMASK;
}
}
tmp->fraction.ll = high;
return pack_d (res);
}
-#endif /* L_mul_sf || L_mul_df */
+#endif /* L_mul_sf || L_mul_df || L_mul_tf */
#if defined(L_div_sf) || defined(L_div_df) || defined(L_div_tf)
static inline __attribute__ ((__always_inline__)) fp_number_type *
{
if (quotient & (1 << NGARDS))
{
- /* half way, so round to even */
- quotient += GARDROUND + 1;
+ /* Because we're half way, we would round to even by adding
+ GARDROUND + 1, except that's also done in the packing
+ function, and rounding twice will lose precision and cause
+ the result to be too far off. */
}
else if (numerator)
{
- /* but we really weren't half way, more bits exist */
+ /* We're a further than half way by the small amount
+ corresponding to the bits set in "numerator". Knowing
+ that, we round here and not in pack_d, because there we
+ don't have "numerator" available anymore. */
quotient += GARDROUND + 1;
+
+ /* Avoid further rounding in pack_d. */
+ quotient &= ~(fractype) GARDMASK;
}
}
}
else
{
+ USItype uarg;
+ int shift;
in.normal_exp = FRACBITS + NGARDS;
if (in.sign)
{
{
return (FLO_type)(- MAX_SI_INT - 1);
}
- in.fraction.ll = (-arg_a);
+ uarg = (-arg_a);
}
else
- in.fraction.ll = arg_a;
+ uarg = arg_a;
- while (in.fraction.ll < ((fractype)1 << (FRACBITS + NGARDS)))
+ in.fraction.ll = uarg;
+ shift = clzusi (uarg) - (BITS_PER_SI - 1 - FRACBITS - NGARDS);
+ if (shift > 0)
{
- in.fraction.ll <<= 1;
- in.normal_exp -= 1;
+ in.fraction.ll <<= shift;
+ in.normal_exp -= shift;
}
}
return pack_d (&in);
}
else
{
+ int shift;
in.class = CLASS_NUMBER;
in.normal_exp = FRACBITS + NGARDS;
in.fraction.ll = arg_a;
- while (in.fraction.ll > ((fractype)1 << (FRACBITS + NGARDS)))
- {
- in.fraction.ll >>= 1;
- in.normal_exp += 1;
- }
- while (in.fraction.ll < ((fractype)1 << (FRACBITS + NGARDS)))
+ shift = clzusi (arg_a) - (BITS_PER_SI - 1 - FRACBITS - NGARDS);
+ if (shift < 0)
+ {
+ fractype guard = in.fraction.ll & (((fractype)1 << -shift) - 1);
+ in.fraction.ll >>= -shift;
+ in.fraction.ll |= (guard != 0);
+ in.normal_exp -= shift;
+ }
+ else if (shift > 0)
{
- in.fraction.ll <<= 1;
- in.normal_exp -= 1;
+ in.fraction.ll <<= shift;
+ in.normal_exp -= shift;
}
}
return pack_d (&in);
OSDN Git Service
