| /* 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, 1995, 1996, 1997 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, 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. */ |
| |
| /* 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 */ |
| |
| /* We don't currently support extended floats (long doubles) on machines |
| without hardware to deal with them. |
| |
| These stubs are just to keep the linker from complaining about unresolved |
| references which can be pulled in from libio & libstdc++, even if the |
| user isn't using long doubles. However, they may generate an unresolved |
| external to abort if abort is not used by the function, and the stubs |
| are referenced from within libc, since libgcc goes before and after the |
| system library. */ |
| |
| #ifdef EXTENDED_FLOAT_STUBS |
| __truncxfsf2 (){ abort(); } |
| __extendsfxf2 (){ abort(); } |
| __addxf3 (){ abort(); } |
| __divxf3 (){ abort(); } |
| __eqxf2 (){ abort(); } |
| __extenddfxf2 (){ abort(); } |
| __gtxf2 (){ abort(); } |
| __lexf2 (){ abort(); } |
| __ltxf2 (){ abort(); } |
| __mulxf3 (){ abort(); } |
| __negxf2 (){ abort(); } |
| __nexf2 (){ abort(); } |
| __subxf3 (){ abort(); } |
| __truncxfdf2 (){ abort(); } |
| |
| __trunctfsf2 (){ abort(); } |
| __extendsftf2 (){ abort(); } |
| __addtf3 (){ abort(); } |
| __divtf3 (){ abort(); } |
| __eqtf2 (){ abort(); } |
| __extenddftf2 (){ abort(); } |
| __gttf2 (){ abort(); } |
| __letf2 (){ abort(); } |
| __lttf2 (){ abort(); } |
| __multf3 (){ abort(); } |
| __negtf2 (){ abort(); } |
| __netf2 (){ abort(); } |
| __subtf3 (){ abort(); } |
| __trunctfdf2 (){ abort(); } |
| #else /* !EXTENDED_FLOAT_STUBS, rest of file */ |
| |
| |
| typedef SFtype __attribute__ ((mode (SF))); |
| typedef DFtype __attribute__ ((mode (DF))); |
| |
| typedef int HItype __attribute__ ((mode (HI))); |
| typedef int SItype __attribute__ ((mode (SI))); |
| typedef int DItype __attribute__ ((mode (DI))); |
| |
| /* The type of the result of a fp compare */ |
| #ifndef CMPtype |
| #define CMPtype SItype |
| #endif |
| |
| typedef unsigned int UHItype __attribute__ ((mode (HI))); |
| typedef unsigned int USItype __attribute__ ((mode (SI))); |
| typedef unsigned int UDItype __attribute__ ((mode (DI))); |
| |
| #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 |
| # define pack_d pack_f |
| # define unpack_d unpack_f |
| 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 |
| |
| |
| #define INLINE __inline__ |
| |
| /* 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; |
| fractype value_raw; |
| |
| #ifndef FLOAT |
| halffractype words[2]; |
| #endif |
| |
| #ifdef FLOAT_BIT_ORDER_MISMATCH |
| struct |
| { |
| fractype fraction:FRACBITS __attribute__ ((packed)); |
| unsigned int exp:EXPBITS __attribute__ ((packed)); |
| unsigned int sign:1 __attribute__ ((packed)); |
| } |
| bits; |
| #endif |
| |
| #ifdef _DEBUG_BITFLOAT |
| struct |
| { |
| unsigned int sign:1 __attribute__ ((packed)); |
| unsigned int exp:EXPBITS __attribute__ ((packed)); |
| fractype fraction:FRACBITS __attribute__ ((packed)); |
| } |
| bits_big_endian; |
| |
| struct |
| { |
| fractype fraction:FRACBITS __attribute__ ((packed)); |
| unsigned int exp:EXPBITS __attribute__ ((packed)); |
| unsigned int sign:1 __attribute__ ((packed)); |
| } |
| bits_little_endian; |
| #endif |
| } |
| 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? */ |
| int sign = src->sign; |
| int exp = 0; |
| |
| if (isnan (src)) |
| { |
| exp = EXPMAX; |
| if (src->class == CLASS_QNAN || 1) |
| { |
| fraction |= QUIET_NAN; |
| } |
| } |
| else if (isinf (src)) |
| { |
| exp = EXPMAX; |
| fraction = 0; |
| } |
| else if (iszero (src)) |
| { |
| exp = 0; |
| fraction = 0; |
| } |
| else if (fraction == 0) |
| { |
| exp = 0; |
| sign = 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; |
| |
| 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; |
| } |
| else if (src->normal_exp > EXPBIAS) |
| { |
| exp = EXPMAX; |
| fraction = 0; |
| } |
| else |
| { |
| 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; |
| exp += 1; |
| } |
| fraction >>= NGARDS; |
| } |
| } |
| |
| /* We previously used bitfields to store the number, but this doesn't |
| handle little/big endian systems conviently, so use shifts and |
| masks */ |
| #ifdef FLOAT_BIT_ORDER_MISMATCH |
| dst.bits.fraction = fraction; |
| dst.bits.exp = exp; |
| dst.bits.sign = sign; |
| #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 |
| |
| #if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT) |
| { |
| halffractype tmp = dst.words[0]; |
| dst.words[0] = dst.words[1]; |
| dst.words[1] = tmp; |
| } |
| #endif |
| |
| return dst.value; |
| } |
| |
| static void |
| unpack_d (FLO_union_type * src, fp_number_type * dst) |
| { |
| /* We previously used bitfields to store the number, but this doesn't |
| handle little/big endian systems conviently, so use shifts and |
| masks */ |
| fractype fraction; |
| int exp; |
| int sign; |
| |
| #if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT) |
| FLO_union_type swapped; |
| |
| swapped.words[0] = src->words[1]; |
| swapped.words[1] = src->words[0]; |
| src = &swapped; |
| #endif |
| |
| #ifdef FLOAT_BIT_ORDER_MISMATCH |
| fraction = src->bits.fraction; |
| exp = src->bits.exp; |
| sign = src->bits.sign; |
| #else |
| fraction = src->value_raw & ((((fractype)1) << FRACBITS) - (fractype)1); |
| exp = ((int)(src->value_raw >> FRACBITS)) & ((1 << EXPBITS) - 1); |
| sign = ((int)(src->value_raw >> (FRACBITS + EXPBITS))) & 1; |
| #endif |
| |
| dst->sign = sign; |
| if (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 = 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 (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 (fraction & QUIET_NAN) |
| { |
| dst->class = CLASS_QNAN; |
| } |
| else |
| { |
| dst->class = CLASS_SNAN; |
| } |
| /* Keep the fraction part as the nan number */ |
| dst->fraction.ll = fraction; |
| } |
| } |
| else |
| { |
| /* Nothing strange about this number */ |
| dst->normal_exp = 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; |
| } |
| |
| a->sign = a->sign ^ b->sign; |
| |
| if (isinf (a) || iszero (a)) |
| { |
| if (a->class == b->class) |
| return nan (); |
| return a; |
| } |
| |
| 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)-1 : MAX_SI_INT; |
| /* 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; |
| /* it is a negative number */ |
| if (a.sign) |
| return 0; |
| /* get reasonable MAX_USI_INT... */ |
| if (isinf (&a)) |
| return MAX_USI_INT; |
| /* 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 << (a.normal_exp - (FRACBITS + NGARDS)); |
| 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; |
| USItype sffrac; |
| |
| unpack_d ((FLO_union_type *) & arg_a, &in); |
| |
| sffrac = in.fraction.ll >> F_D_BITOFF; |
| |
| /* We set the lowest guard bit in SFFRAC if we discarded any non |
| zero bits. */ |
| if ((in.fraction.ll & (((USItype) 1 << F_D_BITOFF) - 1)) != 0) |
| sffrac |= 1; |
| |
| return __make_fp (in.class, in.sign, in.normal_exp, sffrac); |
| } |
| |
| #endif |
| #endif /* !EXTENDED_FLOAT_STUBS */ |