diff options
Diffstat (limited to 'sim/ppc/dp-bit.c')
| -rw-r--r-- | sim/ppc/dp-bit.c | 1307 |
1 files changed, 1307 insertions, 0 deletions
diff --git a/sim/ppc/dp-bit.c b/sim/ppc/dp-bit.c new file mode 100644 index 0000000..3313132 --- /dev/null +++ b/sim/ppc/dp-bit.c @@ -0,0 +1,1307 @@ +/* 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 |