12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298 |
- /* 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-2022 Free Software Foundation, Inc.
- This program 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 3 of the License, or
- (at your option) any later version.
- This program 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. If not, see <http://www.gnu.org/licenses/>. */
- /* 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
|