aboutsummaryrefslogtreecommitdiff
path: root/sim/ppc/dp-bit.c
diff options
context:
space:
mode:
Diffstat (limited to 'sim/ppc/dp-bit.c')
-rw-r--r--sim/ppc/dp-bit.c1307
1 files changed, 0 insertions, 1307 deletions
diff --git a/sim/ppc/dp-bit.c b/sim/ppc/dp-bit.c
deleted file mode 100644
index 3313132..0000000
--- a/sim/ppc/dp-bit.c
+++ /dev/null
@@ -1,1307 +0,0 @@
-/* This is a software floating point library which can be used instead of
- the floating point routines in libgcc1.c for targets without hardware
- floating point. */
-
-/* Copyright (C) 1994 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.
-
-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.
-
-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, 675 Mass Ave, Cambridge, MA 02139, 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. */
-
-/* This implements IEEE 754 format arithmetic, but does not provide a
- mechanism for setting the rounding mode, or for generating or handling
- exceptions.
-
- The original code by Steve Chamberlain, hacked by Mark Eichin and Jim
- Wilson, all of Cygnus Support. */
-
-/* The intended way to use this file is to make two copies, add `#define FLOAT'
- to one copy, then compile both copies and add them to libgcc.a. */
-
-/* The following macros can be defined to change the behaviour of this file:
- FLOAT: Implement a `float', aka SFmode, fp library. If this is not
- defined, then this file implements a `double', aka DFmode, fp library.
- FLOAT_ONLY: Used with FLOAT, to implement a `float' only library, i.e.
- don't include float->double conversion which requires the double library.
- This is useful only for machines which can't support doubles, e.g. some
- 8-bit processors.
- CMPtype: Specify the type that floating point compares should return.
- This defaults to SItype, aka int.
- US_SOFTWARE_GOFAST: This makes all entry points use the same names as the
- US Software goFast library. If this is not defined, the entry points use
- the same names as libgcc1.c.
- _DEBUG_BITFLOAT: This makes debugging the code a little easier, by adding
- two integers to the FLO_union_type.
- NO_NANS: Disable nan and infinity handling
- SMALL_MACHINE: Useful when operations on QIs and HIs are faster
- than on an SI */
-
-#ifndef SFtype
-typedef SFtype __attribute__ ((mode (SF)));
-#endif
-#ifndef DFtype
-typedef DFtype __attribute__ ((mode (DF)));
-#endif
-
-#ifndef HItype
-typedef int HItype __attribute__ ((mode (HI)));
-#endif
-#ifndef SItype
-typedef int SItype __attribute__ ((mode (SI)));
-#endif
-#ifndef DItype
-typedef int DItype __attribute__ ((mode (DI)));
-#endif
-
-/* The type of the result of a fp compare */
-#ifndef CMPtype
-#define CMPtype SItype
-#endif
-
-#ifndef UHItype
-typedef unsigned int UHItype __attribute__ ((mode (HI)));
-#endif
-#ifndef USItype
-typedef unsigned int USItype __attribute__ ((mode (SI)));
-#endif
-#ifndef UDItype
-typedef unsigned int UDItype __attribute__ ((mode (DI)));
-#endif
-
-#define MAX_SI_INT ((SItype) ((unsigned) (~0)>>1))
-#define MAX_USI_INT ((USItype) ~0)
-
-
-#ifdef FLOAT_ONLY
-#define NO_DI_MODE
-#endif
-
-#ifdef FLOAT
-# define NGARDS 7L
-# define GARDROUND 0x3f
-# define GARDMASK 0x7f
-# define GARDMSB 0x40
-# define EXPBITS 8
-# define EXPBIAS 127
-# define FRACBITS 23
-# define EXPMAX (0xff)
-# define QUIET_NAN 0x100000L
-# define FRAC_NBITS 32
-# define FRACHIGH 0x80000000L
-# define FRACHIGH2 0xc0000000L
- typedef USItype fractype;
- typedef UHItype halffractype;
- typedef SFtype FLO_type;
- typedef SItype intfrac;
-
-#else
-# define PREFIXFPDP dp
-# define PREFIXSFDF df
-# define NGARDS 8L
-# define GARDROUND 0x7f
-# define GARDMASK 0xff
-# define GARDMSB 0x80
-# define EXPBITS 11
-# define EXPBIAS 1023
-# define FRACBITS 52
-# define EXPMAX (0x7ff)
-# define QUIET_NAN 0x8000000000000LL
-# define FRAC_NBITS 64
-# define FRACHIGH 0x8000000000000000LL
-# define FRACHIGH2 0xc000000000000000LL
- typedef UDItype fractype;
- typedef USItype halffractype;
- typedef DFtype FLO_type;
- typedef DItype intfrac;
-#endif
-
-#ifdef US_SOFTWARE_GOFAST
-# ifdef FLOAT
-# define add fpadd
-# define sub fpsub
-# define multiply fpmul
-# define divide fpdiv
-# define compare fpcmp
-# define si_to_float sitofp
-# define float_to_si fptosi
-# define float_to_usi fptoui
-# define negate __negsf2
-# define sf_to_df fptodp
-# define dptofp dptofp
-#else
-# define add dpadd
-# define sub dpsub
-# define multiply dpmul
-# define divide dpdiv
-# define compare dpcmp
-# define si_to_float litodp
-# define float_to_si dptoli
-# define float_to_usi dptoul
-# define negate __negdf2
-# define df_to_sf dptofp
-#endif
-#else
-# ifdef FLOAT
-# define add __addsf3
-# define sub __subsf3
-# define multiply __mulsf3
-# define divide __divsf3
-# define compare __cmpsf2
-# define _eq_f2 __eqsf2
-# define _ne_f2 __nesf2
-# define _gt_f2 __gtsf2
-# define _ge_f2 __gesf2
-# define _lt_f2 __ltsf2
-# define _le_f2 __lesf2
-# define si_to_float __floatsisf
-# define float_to_si __fixsfsi
-# define float_to_usi __fixunssfsi
-# define negate __negsf2
-# define sf_to_df __extendsfdf2
-#else
-# define add __adddf3
-# define sub __subdf3
-# define multiply __muldf3
-# define divide __divdf3
-# define compare __cmpdf2
-# define _eq_f2 __eqdf2
-# define _ne_f2 __nedf2
-# define _gt_f2 __gtdf2
-# define _ge_f2 __gedf2
-# define _lt_f2 __ltdf2
-# define _le_f2 __ledf2
-# define si_to_float __floatsidf
-# define float_to_si __fixdfsi
-# define float_to_usi __fixunsdfsi
-# define negate __negdf2
-# define df_to_sf __truncdfsf2
-# endif
-#endif
-
-
-#ifndef INLINE
-#define INLINE __inline__
-#endif
-
-/* Preserve the sticky-bit when shifting fractions to the right. */
-#define LSHIFT(a) { a = (a & 1) | (a >> 1); }
-
-/* numeric parameters */
-/* F_D_BITOFF is the number of bits offset between the MSB of the mantissa
- of a float and of a double. Assumes there are only two float types.
- (double::FRAC_BITS+double::NGARGS-(float::FRAC_BITS-float::NGARDS))
- */
-#define F_D_BITOFF (52+8-(23+7))
-
-
-#define NORMAL_EXPMIN (-(EXPBIAS)+1)
-#define IMPLICIT_1 (1LL<<(FRACBITS+NGARDS))
-#define IMPLICIT_2 (1LL<<(FRACBITS+1+NGARDS))
-
-/* common types */
-
-typedef enum
-{
- CLASS_SNAN,
- CLASS_QNAN,
- CLASS_ZERO,
- CLASS_NUMBER,
- CLASS_INFINITY
-} fp_class_type;
-
-typedef struct
-{
-#ifdef SMALL_MACHINE
- char class;
- unsigned char sign;
- short normal_exp;
-#else
- fp_class_type class;
- unsigned int sign;
- int normal_exp;
-#endif
-
- union
- {
- fractype ll;
- halffractype l[2];
- } fraction;
-} fp_number_type;
-
-typedef union
-{
- FLO_type value;
-#ifdef _DEBUG_BITFLOAT
- int l[2];
-#endif
- struct
- {
-#ifndef FLOAT_BIT_ORDER_MISMATCH
- unsigned int sign:1 __attribute__ ((packed));
- unsigned int exp:EXPBITS __attribute__ ((packed));
- fractype fraction:FRACBITS __attribute__ ((packed));
-#else
- fractype fraction:FRACBITS __attribute__ ((packed));
- unsigned int exp:EXPBITS __attribute__ ((packed));
- unsigned int sign:1 __attribute__ ((packed));
-#endif
- }
- bits;
-}
-FLO_union_type;
-
-
-/* end of header */
-
-/* IEEE "special" number predicates */
-
-#ifdef NO_NANS
-
-#define nan() 0
-#define isnan(x) 0
-#define isinf(x) 0
-#else
-
-INLINE
-static fp_number_type *
-nan ()
-{
- static fp_number_type thenan;
-
- return &thenan;
-}
-
-INLINE
-static int
-isnan ( fp_number_type * x)
-{
- return x->class == CLASS_SNAN || x->class == CLASS_QNAN;
-}
-
-INLINE
-static int
-isinf ( fp_number_type * x)
-{
- return x->class == CLASS_INFINITY;
-}
-
-#endif
-
-INLINE
-static int
-iszero ( fp_number_type * x)
-{
- return x->class == CLASS_ZERO;
-}
-
-INLINE
-static void
-flip_sign ( fp_number_type * x)
-{
- x->sign = !x->sign;
-}
-
-static FLO_type
-pack_d ( fp_number_type * src)
-{
- FLO_union_type dst;
- fractype fraction = src->fraction.ll; /* wasn't unsigned before? */
-
- dst.bits.sign = src->sign;
-
- if (isnan (src))
- {
- dst.bits.exp = EXPMAX;
- dst.bits.fraction = src->fraction.ll;
- if (src->class == CLASS_QNAN || 1)
- {
- dst.bits.fraction |= QUIET_NAN;
- }
- }
- else if (isinf (src))
- {
- dst.bits.exp = EXPMAX;
- dst.bits.fraction = 0;
- }
- else if (iszero (src))
- {
- dst.bits.exp = 0;
- dst.bits.fraction = 0;
- }
- else if (fraction == 0)
- {
- dst.value = 0;
- }
- else
- {
- if (src->normal_exp < NORMAL_EXPMIN)
- {
- /* This number's exponent is too low to fit into the bits
- available in the number, so we'll store 0 in the exponent and
- shift the fraction to the right to make up for it. */
-
- int shift = NORMAL_EXPMIN - src->normal_exp;
-
- dst.bits.exp = 0;
-
- if (shift > FRAC_NBITS - NGARDS)
- {
- /* No point shifting, since it's more that 64 out. */
- fraction = 0;
- }
- else
- {
- /* Shift by the value */
- fraction >>= shift;
- }
- fraction >>= NGARDS;
- dst.bits.fraction = fraction;
- }
- else if (src->normal_exp > EXPBIAS)
- {
- dst.bits.exp = EXPMAX;
- dst.bits.fraction = 0;
- }
- else
- {
- dst.bits.exp = src->normal_exp + EXPBIAS;
- /* IF the gard bits are the all zero, but the first, then we're
- half way between two numbers, choose the one which makes the
- lsb of the answer 0. */
- if ((fraction & GARDMASK) == GARDMSB)
- {
- if (fraction & (1 << NGARDS))
- fraction += GARDROUND + 1;
- }
- else
- {
- /* Add a one to the guards to round up */
- fraction += GARDROUND;
- }
- if (fraction >= IMPLICIT_2)
- {
- fraction >>= 1;
- dst.bits.exp += 1;
- }
- fraction >>= NGARDS;
- dst.bits.fraction = fraction;
- }
- }
- return dst.value;
-}
-
-static void
-unpack_d (FLO_union_type * src, fp_number_type * dst)
-{
- fractype fraction = src->bits.fraction;
-
- dst->sign = src->bits.sign;
- if (src->bits.exp == 0)
- {
- /* Hmm. Looks like 0 */
- if (fraction == 0)
- {
- /* tastes like zero */
- dst->class = CLASS_ZERO;
- }
- else
- {
- /* Zero exponent with non zero fraction - it's denormalized,
- so there isn't a leading implicit one - we'll shift it so
- it gets one. */
- dst->normal_exp = src->bits.exp - EXPBIAS + 1;
- fraction <<= NGARDS;
-
- dst->class = CLASS_NUMBER;
-#if 1
- while (fraction < IMPLICIT_1)
- {
- fraction <<= 1;
- dst->normal_exp--;
- }
-#endif
- dst->fraction.ll = fraction;
- }
- }
- else if (src->bits.exp == EXPMAX)
- {
- /* Huge exponent*/
- if (fraction == 0)
- {
- /* Attached to a zero fraction - means infinity */
- dst->class = CLASS_INFINITY;
- }
- else
- {
- /* Non zero fraction, means nan */
- if (dst->sign)
- {
- dst->class = CLASS_SNAN;
- }
- else
- {
- dst->class = CLASS_QNAN;
- }
- /* Keep the fraction part as the nan number */
- dst->fraction.ll = fraction;
- }
- }
- else
- {
- /* Nothing strange about this number */
- dst->normal_exp = src->bits.exp - EXPBIAS;
- dst->class = CLASS_NUMBER;
- dst->fraction.ll = (fraction << NGARDS) | IMPLICIT_1;
- }
-}
-
-static fp_number_type *
-_fpadd_parts (fp_number_type * a,
- fp_number_type * b,
- fp_number_type * tmp)
-{
- intfrac tfraction;
-
- /* Put commonly used fields in local variables. */
- int a_normal_exp;
- int b_normal_exp;
- fractype a_fraction;
- fractype b_fraction;
-
- if (isnan (a))
- {
- return a;
- }
- if (isnan (b))
- {
- return b;
- }
- if (isinf (a))
- {
- /* Adding infinities with opposite signs yields a NaN. */
- if (isinf (b) && a->sign != b->sign)
- return nan ();
- return a;
- }
- if (isinf (b))
- {
- return b;
- }
- if (iszero (b))
- {
- return a;
- }
- if (iszero (a))
- {
- return b;
- }
-
- /* Got two numbers. shift the smaller and increment the exponent till
- they're the same */
- {
- int diff;
-
- a_normal_exp = a->normal_exp;
- b_normal_exp = b->normal_exp;
- a_fraction = a->fraction.ll;
- b_fraction = b->fraction.ll;
-
- diff = a_normal_exp - b_normal_exp;
-
- 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)
- {
- b_normal_exp++;
- LSHIFT (b_fraction);
- }
- while (b_normal_exp > a_normal_exp)
- {
- a_normal_exp++;
- LSHIFT (a_fraction);
- }
- }
- else
- {
- /* Somethings's up.. choose the biggest */
- if (a_normal_exp > b_normal_exp)
- {
- b_normal_exp = a_normal_exp;
- b_fraction = 0;
- }
- else
- {
- a_normal_exp = b_normal_exp;
- a_fraction = 0;
- }
- }
- }
-
- if (a->sign != b->sign)
- {
- if (a->sign)
- {
- tfraction = -a_fraction + b_fraction;
- }
- else
- {
- tfraction = a_fraction - b_fraction;
- }
- if (tfraction > 0)
- {
- tmp->sign = 0;
- tmp->normal_exp = a_normal_exp;
- tmp->fraction.ll = tfraction;
- }
- else
- {
- tmp->sign = 1;
- tmp->normal_exp = a_normal_exp;
- tmp->fraction.ll = -tfraction;
- }
- /* and renormalize it */
-
- while (tmp->fraction.ll < IMPLICIT_1 && tmp->fraction.ll)
- {
- tmp->fraction.ll <<= 1;
- tmp->normal_exp--;
- }
- }
- else
- {
- tmp->sign = a->sign;
- tmp->normal_exp = a_normal_exp;
- tmp->fraction.ll = a_fraction + b_fraction;
- }
- tmp->class = CLASS_NUMBER;
- /* Now the fraction is added, we have to shift down to renormalize the
- number */
-
- if (tmp->fraction.ll >= IMPLICIT_2)
- {
- LSHIFT (tmp->fraction.ll);
- tmp->normal_exp++;
- }
- return tmp;
-
-}
-
-FLO_type
-add (FLO_type arg_a, FLO_type arg_b)
-{
- fp_number_type a;
- fp_number_type b;
- fp_number_type tmp;
- fp_number_type *res;
-
- unpack_d ((FLO_union_type *) & arg_a, &a);
- unpack_d ((FLO_union_type *) & arg_b, &b);
-
- res = _fpadd_parts (&a, &b, &tmp);
-
- return pack_d (res);
-}
-
-FLO_type
-sub (FLO_type arg_a, FLO_type arg_b)
-{
- fp_number_type a;
- fp_number_type b;
- fp_number_type tmp;
- fp_number_type *res;
-
- unpack_d ((FLO_union_type *) & arg_a, &a);
- unpack_d ((FLO_union_type *) & arg_b, &b);
-
- b.sign ^= 1;
-
- res = _fpadd_parts (&a, &b, &tmp);
-
- return pack_d (res);
-}
-
-static fp_number_type *
-_fpmul_parts ( fp_number_type * a,
- fp_number_type * b,
- fp_number_type * tmp)
-{
- fractype low = 0;
- fractype high = 0;
-
- if (isnan (a))
- {
- a->sign = a->sign != b->sign;
- return a;
- }
- if (isnan (b))
- {
- b->sign = a->sign != b->sign;
- return b;
- }
- if (isinf (a))
- {
- if (iszero (b))
- return nan ();
- a->sign = a->sign != b->sign;
- return a;
- }
- if (isinf (b))
- {
- if (iszero (a))
- {
- return nan ();
- }
- b->sign = a->sign != b->sign;
- return b;
- }
- if (iszero (a))
- {
- a->sign = a->sign != b->sign;
- return a;
- }
- if (iszero (b))
- {
- b->sign = a->sign != b->sign;
- return b;
- }
-
- /* Calculate the mantissa by multiplying both 64bit numbers to get a
- 128 bit number */
- {
- fractype x = a->fraction.ll;
- fractype ylow = b->fraction.ll;
- fractype yhigh = 0;
- int bit;
-
-#if defined(NO_DI_MODE)
- {
- /* ??? This does multiplies one bit at a time. Optimize. */
- for (bit = 0; bit < FRAC_NBITS; bit++)
- {
- int carry;
-
- if (x & 1)
- {
- carry = (low += ylow) < ylow;
- high += yhigh + carry;
- }
- yhigh <<= 1;
- if (ylow & FRACHIGH)
- {
- yhigh |= 1;
- }
- ylow <<= 1;
- x >>= 1;
- }
- }
-#elif defined(FLOAT)
- {
- /* Multiplying two 32 bit numbers to get a 64 bit number on
- a machine with DI, so we're safe */
-
- DItype answer = (DItype)(a->fraction.ll) * (DItype)(b->fraction.ll);
-
- high = answer >> 32;
- low = answer;
- }
-#else
- /* Doing a 64*64 to 128 */
- {
- UDItype nl = a->fraction.ll & 0xffffffff;
- UDItype nh = a->fraction.ll >> 32;
- UDItype ml = b->fraction.ll & 0xffffffff;
- UDItype mh = b->fraction.ll >>32;
- UDItype pp_ll = ml * nl;
- UDItype pp_hl = mh * nl;
- UDItype pp_lh = ml * nh;
- UDItype pp_hh = mh * nh;
- UDItype res2 = 0;
- UDItype res0 = 0;
- UDItype ps_hh__ = pp_hl + pp_lh;
- if (ps_hh__ < pp_hl)
- res2 += 0x100000000LL;
- pp_hl = (ps_hh__ << 32) & 0xffffffff00000000LL;
- res0 = pp_ll + pp_hl;
- if (res0 < pp_ll)
- res2++;
- res2 += ((ps_hh__ >> 32) & 0xffffffffL) + pp_hh;
- high = res2;
- low = res0;
- }
-#endif
- }
-
- tmp->normal_exp = a->normal_exp + b->normal_exp;
- tmp->sign = a->sign != b->sign;
-#ifdef FLOAT
- tmp->normal_exp += 2; /* ??????????????? */
-#else
- tmp->normal_exp += 4; /* ??????????????? */
-#endif
- while (high >= IMPLICIT_2)
- {
- tmp->normal_exp++;
- if (high & 1)
- {
- low >>= 1;
- low |= FRACHIGH;
- }
- high >>= 1;
- }
- while (high < IMPLICIT_1)
- {
- tmp->normal_exp--;
-
- high <<= 1;
- if (low & FRACHIGH)
- 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 ((high & GARDMASK) == GARDMSB)
- {
- if (high & (1 << NGARDS))
- {
- /* half way, so round to even */
- high += GARDROUND + 1;
- }
- else if (low)
- {
- /* but we really weren't half way */
- high += GARDROUND + 1;
- }
- }
- tmp->fraction.ll = high;
- tmp->class = CLASS_NUMBER;
- return tmp;
-}
-
-FLO_type
-multiply (FLO_type arg_a, FLO_type arg_b)
-{
- fp_number_type a;
- fp_number_type b;
- fp_number_type tmp;
- fp_number_type *res;
-
- unpack_d ((FLO_union_type *) & arg_a, &a);
- unpack_d ((FLO_union_type *) & arg_b, &b);
-
- res = _fpmul_parts (&a, &b, &tmp);
-
- return pack_d (res);
-}
-
-static fp_number_type *
-_fpdiv_parts (fp_number_type * a,
- fp_number_type * b,
- fp_number_type * tmp)
-{
- fractype low = 0;
- fractype high = 0;
- fractype r0, r1, y0, y1, bit;
- fractype q;
- fractype numerator;
- fractype denominator;
- fractype quotient;
- fractype remainder;
-
- if (isnan (a))
- {
- return a;
- }
- if (isnan (b))
- {
- return b;
- }
- if (isinf (a) || iszero (a))
- {
- if (a->class == b->class)
- return nan ();
- return a;
- }
- a->sign = a->sign ^ b->sign;
-
- if (isinf (b))
- {
- a->fraction.ll = 0;
- a->normal_exp = 0;
- return a;
- }
- if (iszero (b))
- {
- a->class = CLASS_INFINITY;
- return b;
- }
-
- /* Calculate the mantissa by multiplying both 64bit numbers to get a
- 128 bit number */
- {
- int carry;
- intfrac d0, d1; /* weren't unsigned before ??? */
-
- /* quotient =
- ( numerator / denominator) * 2^(numerator exponent - denominator exponent)
- */
-
- a->normal_exp = a->normal_exp - b->normal_exp;
- numerator = a->fraction.ll;
- denominator = b->fraction.ll;
-
- if (numerator < denominator)
- {
- /* Fraction will be less than 1.0 */
- numerator *= 2;
- a->normal_exp--;
- }
- bit = IMPLICIT_1;
- quotient = 0;
- /* ??? Does divide one bit at a time. Optimize. */
- while (bit)
- {
- if (numerator >= denominator)
- {
- quotient |= bit;
- numerator -= denominator;
- }
- bit >>= 1;
- numerator *= 2;
- }
-
- if ((quotient & GARDMASK) == GARDMSB)
- {
- if (quotient & (1 << NGARDS))
- {
- /* half way, so round to even */
- quotient += GARDROUND + 1;
- }
- else if (numerator)
- {
- /* but we really weren't half way, more bits exist */
- quotient += GARDROUND + 1;
- }
- }
-
- a->fraction.ll = quotient;
- return (a);
- }
-}
-
-FLO_type
-divide (FLO_type arg_a, FLO_type arg_b)
-{
- fp_number_type a;
- fp_number_type b;
- fp_number_type tmp;
- fp_number_type *res;
-
- unpack_d ((FLO_union_type *) & arg_a, &a);
- unpack_d ((FLO_union_type *) & arg_b, &b);
-
- res = _fpdiv_parts (&a, &b, &tmp);
-
- return pack_d (res);
-}
-
-/* according to the demo, fpcmp returns a comparison with 0... thus
- a<b -> -1
- a==b -> 0
- a>b -> +1
- */
-
-static int
-_fpcmp_parts (fp_number_type * a, fp_number_type * b)
-{
-#if 0
- /* either nan -> unordered. Must be checked outside of this routine. */
- if (isnan (a) && isnan (b))
- {
- return 1; /* still unordered! */
- }
-#endif
-
- if (isnan (a) || isnan (b))
- {
- return 1; /* how to indicate unordered compare? */
- }
- if (isinf (a) && isinf (b))
- {
- /* +inf > -inf, but +inf != +inf */
- /* b \a| +inf(0)| -inf(1)
- ______\+--------+--------
- +inf(0)| a==b(0)| a<b(-1)
- -------+--------+--------
- -inf(1)| a>b(1) | a==b(0)
- -------+--------+--------
- So since unordered must be non zero, just line up the columns...
- */
- return b->sign - a->sign;
- }
- /* but not both... */
- if (isinf (a))
- {
- return a->sign ? -1 : 1;
- }
- if (isinf (b))
- {
- return b->sign ? 1 : -1;
- }
- if (iszero (a) && iszero (b))
- {
- return 0;
- }
- if (iszero (a))
- {
- return b->sign ? 1 : -1;
- }
- if (iszero (b))
- {
- return a->sign ? -1 : 1;
- }
- /* now both are "normal". */
- if (a->sign != b->sign)
- {
- /* opposite signs */
- return a->sign ? -1 : 1;
- }
- /* same sign; exponents? */
- if (a->normal_exp > b->normal_exp)
- {
- return a->sign ? -1 : 1;
- }
- if (a->normal_exp < b->normal_exp)
- {
- return a->sign ? 1 : -1;
- }
- /* same exponents; check size. */
- if (a->fraction.ll > b->fraction.ll)
- {
- return a->sign ? -1 : 1;
- }
- if (a->fraction.ll < b->fraction.ll)
- {
- return a->sign ? 1 : -1;
- }
- /* after all that, they're equal. */
- return 0;
-}
-
-CMPtype
-compare (FLO_type arg_a, FLO_type arg_b)
-{
- fp_number_type a;
- fp_number_type b;
-
- unpack_d ((FLO_union_type *) & arg_a, &a);
- unpack_d ((FLO_union_type *) & arg_b, &b);
-
- return _fpcmp_parts (&a, &b);
-}
-
-#ifndef US_SOFTWARE_GOFAST
-
-/* These should be optimized for their specific tasks someday. */
-
-CMPtype
-_eq_f2 (FLO_type arg_a, FLO_type arg_b)
-{
- fp_number_type a;
- fp_number_type b;
-
- unpack_d ((FLO_union_type *) & arg_a, &a);
- unpack_d ((FLO_union_type *) & arg_b, &b);
-
- if (isnan (&a) || isnan (&b))
- return 1; /* false, truth == 0 */
-
- return _fpcmp_parts (&a, &b) ;
-}
-
-CMPtype
-_ne_f2 (FLO_type arg_a, FLO_type arg_b)
-{
- fp_number_type a;
- fp_number_type b;
-
- unpack_d ((FLO_union_type *) & arg_a, &a);
- unpack_d ((FLO_union_type *) & arg_b, &b);
-
- if (isnan (&a) || isnan (&b))
- return 1; /* true, truth != 0 */
-
- return _fpcmp_parts (&a, &b) ;
-}
-
-CMPtype
-_gt_f2 (FLO_type arg_a, FLO_type arg_b)
-{
- fp_number_type a;
- fp_number_type b;
-
- unpack_d ((FLO_union_type *) & arg_a, &a);
- unpack_d ((FLO_union_type *) & arg_b, &b);
-
- if (isnan (&a) || isnan (&b))
- return -1; /* false, truth > 0 */
-
- return _fpcmp_parts (&a, &b);
-}
-
-CMPtype
-_ge_f2 (FLO_type arg_a, FLO_type arg_b)
-{
- fp_number_type a;
- fp_number_type b;
-
- unpack_d ((FLO_union_type *) & arg_a, &a);
- unpack_d ((FLO_union_type *) & arg_b, &b);
-
- if (isnan (&a) || isnan (&b))
- return -1; /* false, truth >= 0 */
- return _fpcmp_parts (&a, &b) ;
-}
-
-CMPtype
-_lt_f2 (FLO_type arg_a, FLO_type arg_b)
-{
- fp_number_type a;
- fp_number_type b;
-
- unpack_d ((FLO_union_type *) & arg_a, &a);
- unpack_d ((FLO_union_type *) & arg_b, &b);
-
- if (isnan (&a) || isnan (&b))
- return 1; /* false, truth < 0 */
-
- return _fpcmp_parts (&a, &b);
-}
-
-CMPtype
-_le_f2 (FLO_type arg_a, FLO_type arg_b)
-{
- fp_number_type a;
- fp_number_type b;
-
- unpack_d ((FLO_union_type *) & arg_a, &a);
- unpack_d ((FLO_union_type *) & arg_b, &b);
-
- if (isnan (&a) || isnan (&b))
- return 1; /* false, truth <= 0 */
-
- return _fpcmp_parts (&a, &b) ;
-}
-
-#endif /* ! US_SOFTWARE_GOFAST */
-
-FLO_type
-si_to_float (SItype arg_a)
-{
- fp_number_type in;
-
- in.class = CLASS_NUMBER;
- in.sign = arg_a < 0;
- if (!arg_a)
- {
- in.class = CLASS_ZERO;
- }
- else
- {
- in.normal_exp = FRACBITS + NGARDS;
- if (in.sign)
- {
- /* Special case for minint, since there is no +ve integer
- representation for it */
- if (arg_a == 0x80000000)
- {
- return -2147483648.0;
- }
- in.fraction.ll = (-arg_a);
- }
- else
- in.fraction.ll = arg_a;
-
- while (in.fraction.ll < (1LL << (FRACBITS + NGARDS)))
- {
- in.fraction.ll <<= 1;
- in.normal_exp -= 1;
- }
- }
- return pack_d (&in);
-}
-
-SItype
-float_to_si (FLO_type arg_a)
-{
- fp_number_type a;
- SItype tmp;
-
- unpack_d ((FLO_union_type *) & arg_a, &a);
- if (iszero (&a))
- return 0;
- if (isnan (&a))
- return 0;
- /* get reasonable MAX_SI_INT... */
- if (isinf (&a))
- return a.sign ? MAX_SI_INT : (-MAX_SI_INT)-1;
- /* it is a number, but a small one */
- if (a.normal_exp < 0)
- return 0;
- if (a.normal_exp > 30)
- return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT;
- tmp = a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
- return a.sign ? (-tmp) : (tmp);
-}
-
-#ifdef US_SOFTWARE_GOFAST
-/* While libgcc2.c defines its own __fixunssfsi and __fixunsdfsi routines,
- we also define them for GOFAST because the ones in libgcc2.c have the
- wrong names and I'd rather define these here and keep GOFAST CYG-LOC's
- out of libgcc2.c. We can't define these here if not GOFAST because then
- there'd be duplicate copies. */
-
-USItype
-float_to_usi (FLO_type arg_a)
-{
- fp_number_type a;
-
- unpack_d ((FLO_union_type *) & arg_a, &a);
- if (iszero (&a))
- return 0;
- if (isnan (&a))
- return 0;
- /* get reasonable MAX_USI_INT... */
- if (isinf (&a))
- return a.sign ? MAX_USI_INT : 0;
- /* it is a negative number */
- if (a.sign)
- return 0;
- /* it is a number, but a small one */
- if (a.normal_exp < 0)
- return 0;
- if (a.normal_exp > 31)
- return MAX_USI_INT;
- else if (a.normal_exp > (FRACBITS + NGARDS))
- return a.fraction.ll << ((FRACBITS + NGARDS) - a.normal_exp);
- else
- return a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
-}
-#endif
-
-FLO_type
-negate (FLO_type arg_a)
-{
- fp_number_type a;
-
- unpack_d ((FLO_union_type *) & arg_a, &a);
- flip_sign (&a);
- return pack_d (&a);
-}
-
-#ifdef FLOAT
-
-SFtype
-__make_fp(fp_class_type class,
- unsigned int sign,
- int exp,
- USItype frac)
-{
- fp_number_type in;
-
- in.class = class;
- in.sign = sign;
- in.normal_exp = exp;
- in.fraction.ll = frac;
- return pack_d (&in);
-}
-
-#ifndef FLOAT_ONLY
-
-/* This enables one to build an fp library that supports float but not double.
- Otherwise, we would get an undefined reference to __make_dp.
- This is needed for some 8-bit ports that can't handle well values that
- are 8-bytes in size, so we just don't support double for them at all. */
-
-extern DFtype __make_dp (fp_class_type, unsigned int, int, UDItype frac);
-
-DFtype
-sf_to_df (SFtype arg_a)
-{
- fp_number_type in;
-
- unpack_d ((FLO_union_type *) & arg_a, &in);
- return __make_dp (in.class, in.sign, in.normal_exp,
- ((UDItype) in.fraction.ll) << F_D_BITOFF);
-}
-
-#endif
-#endif
-
-#ifndef FLOAT
-
-extern SFtype __make_fp (fp_class_type, unsigned int, int, USItype);
-
-DFtype
-__make_dp (fp_class_type class, unsigned int sign, int exp, UDItype frac)
-{
- fp_number_type in;
-
- in.class = class;
- in.sign = sign;
- in.normal_exp = exp;
- in.fraction.ll = frac;
- return pack_d (&in);
-}
-
-SFtype
-df_to_sf (DFtype arg_a)
-{
- fp_number_type in;
-
- unpack_d ((FLO_union_type *) & arg_a, &in);
- return __make_fp (in.class, in.sign, in.normal_exp,
- in.fraction.ll >> F_D_BITOFF);
-}
-
-#endif