1*c87b03e5Sespie /* This is a software floating point library which can be used
2*c87b03e5Sespie for targets without hardware floating point.
3*c87b03e5Sespie Copyright (C) 1994, 1995, 1996, 1997, 1998, 2000, 2001, 2002, 2003
4*c87b03e5Sespie Free Software Foundation, Inc.
5*c87b03e5Sespie
6*c87b03e5Sespie This file is free software; you can redistribute it and/or modify it
7*c87b03e5Sespie under the terms of the GNU General Public License as published by the
8*c87b03e5Sespie Free Software Foundation; either version 2, or (at your option) any
9*c87b03e5Sespie later version.
10*c87b03e5Sespie
11*c87b03e5Sespie In addition to the permissions in the GNU General Public License, the
12*c87b03e5Sespie Free Software Foundation gives you unlimited permission to link the
13*c87b03e5Sespie compiled version of this file with other programs, and to distribute
14*c87b03e5Sespie those programs without any restriction coming from the use of this
15*c87b03e5Sespie file. (The General Public License restrictions do apply in other
16*c87b03e5Sespie respects; for example, they cover modification of the file, and
17*c87b03e5Sespie distribution when not linked into another program.)
18*c87b03e5Sespie
19*c87b03e5Sespie This file is distributed in the hope that it will be useful, but
20*c87b03e5Sespie WITHOUT ANY WARRANTY; without even the implied warranty of
21*c87b03e5Sespie MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22*c87b03e5Sespie General Public License for more details.
23*c87b03e5Sespie
24*c87b03e5Sespie You should have received a copy of the GNU General Public License
25*c87b03e5Sespie along with this program; see the file COPYING. If not, write to
26*c87b03e5Sespie the Free Software Foundation, 59 Temple Place - Suite 330,
27*c87b03e5Sespie Boston, MA 02111-1307, USA. */
28*c87b03e5Sespie
29*c87b03e5Sespie /* As a special exception, if you link this library with other files,
30*c87b03e5Sespie some of which are compiled with GCC, to produce an executable,
31*c87b03e5Sespie this library does not by itself cause the resulting executable
32*c87b03e5Sespie to be covered by the GNU General Public License.
33*c87b03e5Sespie This exception does not however invalidate any other reasons why
34*c87b03e5Sespie the executable file might be covered by the GNU General Public License. */
35*c87b03e5Sespie
36*c87b03e5Sespie /* This implements IEEE 754 format arithmetic, but does not provide a
37*c87b03e5Sespie mechanism for setting the rounding mode, or for generating or handling
38*c87b03e5Sespie exceptions.
39*c87b03e5Sespie
40*c87b03e5Sespie The original code by Steve Chamberlain, hacked by Mark Eichin and Jim
41*c87b03e5Sespie Wilson, all of Cygnus Support. */
42*c87b03e5Sespie
43*c87b03e5Sespie /* The intended way to use this file is to make two copies, add `#define FLOAT'
44*c87b03e5Sespie to one copy, then compile both copies and add them to libgcc.a. */
45*c87b03e5Sespie
46*c87b03e5Sespie #include "tconfig.h"
47*c87b03e5Sespie #include "fp-bit.h"
48*c87b03e5Sespie
49*c87b03e5Sespie /* The following macros can be defined to change the behavior of this file:
50*c87b03e5Sespie FLOAT: Implement a `float', aka SFmode, fp library. If this is not
51*c87b03e5Sespie defined, then this file implements a `double', aka DFmode, fp library.
52*c87b03e5Sespie FLOAT_ONLY: Used with FLOAT, to implement a `float' only library, i.e.
53*c87b03e5Sespie don't include float->double conversion which requires the double library.
54*c87b03e5Sespie This is useful only for machines which can't support doubles, e.g. some
55*c87b03e5Sespie 8-bit processors.
56*c87b03e5Sespie CMPtype: Specify the type that floating point compares should return.
57*c87b03e5Sespie This defaults to SItype, aka int.
58*c87b03e5Sespie US_SOFTWARE_GOFAST: This makes all entry points use the same names as the
59*c87b03e5Sespie US Software goFast library.
60*c87b03e5Sespie _DEBUG_BITFLOAT: This makes debugging the code a little easier, by adding
61*c87b03e5Sespie two integers to the FLO_union_type.
62*c87b03e5Sespie NO_DENORMALS: Disable handling of denormals.
63*c87b03e5Sespie NO_NANS: Disable nan and infinity handling
64*c87b03e5Sespie SMALL_MACHINE: Useful when operations on QIs and HIs are faster
65*c87b03e5Sespie than on an SI */
66*c87b03e5Sespie
67*c87b03e5Sespie /* We don't currently support extended floats (long doubles) on machines
68*c87b03e5Sespie without hardware to deal with them.
69*c87b03e5Sespie
70*c87b03e5Sespie These stubs are just to keep the linker from complaining about unresolved
71*c87b03e5Sespie references which can be pulled in from libio & libstdc++, even if the
72*c87b03e5Sespie user isn't using long doubles. However, they may generate an unresolved
73*c87b03e5Sespie external to abort if abort is not used by the function, and the stubs
74*c87b03e5Sespie are referenced from within libc, since libgcc goes before and after the
75*c87b03e5Sespie system library. */
76*c87b03e5Sespie
77*c87b03e5Sespie #ifdef DECLARE_LIBRARY_RENAMES
78*c87b03e5Sespie DECLARE_LIBRARY_RENAMES
79*c87b03e5Sespie #endif
80*c87b03e5Sespie
81*c87b03e5Sespie #ifdef EXTENDED_FLOAT_STUBS
82*c87b03e5Sespie extern void abort (void);
__extendsfxf2(void)83*c87b03e5Sespie void __extendsfxf2 (void) { abort(); }
__extenddfxf2(void)84*c87b03e5Sespie void __extenddfxf2 (void) { abort(); }
__truncxfdf2(void)85*c87b03e5Sespie void __truncxfdf2 (void) { abort(); }
__truncxfsf2(void)86*c87b03e5Sespie void __truncxfsf2 (void) { abort(); }
__fixxfsi(void)87*c87b03e5Sespie void __fixxfsi (void) { abort(); }
__floatsixf(void)88*c87b03e5Sespie void __floatsixf (void) { abort(); }
__addxf3(void)89*c87b03e5Sespie void __addxf3 (void) { abort(); }
__subxf3(void)90*c87b03e5Sespie void __subxf3 (void) { abort(); }
__mulxf3(void)91*c87b03e5Sespie void __mulxf3 (void) { abort(); }
__divxf3(void)92*c87b03e5Sespie void __divxf3 (void) { abort(); }
__negxf2(void)93*c87b03e5Sespie void __negxf2 (void) { abort(); }
__eqxf2(void)94*c87b03e5Sespie void __eqxf2 (void) { abort(); }
__nexf2(void)95*c87b03e5Sespie void __nexf2 (void) { abort(); }
__gtxf2(void)96*c87b03e5Sespie void __gtxf2 (void) { abort(); }
__gexf2(void)97*c87b03e5Sespie void __gexf2 (void) { abort(); }
__lexf2(void)98*c87b03e5Sespie void __lexf2 (void) { abort(); }
__ltxf2(void)99*c87b03e5Sespie void __ltxf2 (void) { abort(); }
100*c87b03e5Sespie
__extendsftf2(void)101*c87b03e5Sespie void __extendsftf2 (void) { abort(); }
__extenddftf2(void)102*c87b03e5Sespie void __extenddftf2 (void) { abort(); }
__trunctfdf2(void)103*c87b03e5Sespie void __trunctfdf2 (void) { abort(); }
__trunctfsf2(void)104*c87b03e5Sespie void __trunctfsf2 (void) { abort(); }
__fixtfsi(void)105*c87b03e5Sespie void __fixtfsi (void) { abort(); }
__floatsitf(void)106*c87b03e5Sespie void __floatsitf (void) { abort(); }
__addtf3(void)107*c87b03e5Sespie void __addtf3 (void) { abort(); }
__subtf3(void)108*c87b03e5Sespie void __subtf3 (void) { abort(); }
__multf3(void)109*c87b03e5Sespie void __multf3 (void) { abort(); }
__divtf3(void)110*c87b03e5Sespie void __divtf3 (void) { abort(); }
__negtf2(void)111*c87b03e5Sespie void __negtf2 (void) { abort(); }
__eqtf2(void)112*c87b03e5Sespie void __eqtf2 (void) { abort(); }
__netf2(void)113*c87b03e5Sespie void __netf2 (void) { abort(); }
__gttf2(void)114*c87b03e5Sespie void __gttf2 (void) { abort(); }
__getf2(void)115*c87b03e5Sespie void __getf2 (void) { abort(); }
__letf2(void)116*c87b03e5Sespie void __letf2 (void) { abort(); }
__lttf2(void)117*c87b03e5Sespie void __lttf2 (void) { abort(); }
118*c87b03e5Sespie #else /* !EXTENDED_FLOAT_STUBS, rest of file */
119*c87b03e5Sespie
120*c87b03e5Sespie /* IEEE "special" number predicates */
121*c87b03e5Sespie
122*c87b03e5Sespie #ifdef NO_NANS
123*c87b03e5Sespie
124*c87b03e5Sespie #define nan() 0
125*c87b03e5Sespie #define isnan(x) 0
126*c87b03e5Sespie #define isinf(x) 0
127*c87b03e5Sespie #else
128*c87b03e5Sespie
129*c87b03e5Sespie #if defined L_thenan_sf
130*c87b03e5Sespie const fp_number_type __thenan_sf = { CLASS_SNAN, 0, 0, {(fractype) 0} };
131*c87b03e5Sespie #elif defined L_thenan_df
132*c87b03e5Sespie const fp_number_type __thenan_df = { CLASS_SNAN, 0, 0, {(fractype) 0} };
133*c87b03e5Sespie #elif defined L_thenan_tf
134*c87b03e5Sespie const fp_number_type __thenan_tf = { CLASS_SNAN, 0, 0, {(fractype) 0} };
135*c87b03e5Sespie #elif defined TFLOAT
136*c87b03e5Sespie extern const fp_number_type __thenan_tf;
137*c87b03e5Sespie #elif defined FLOAT
138*c87b03e5Sespie extern const fp_number_type __thenan_sf;
139*c87b03e5Sespie #else
140*c87b03e5Sespie extern const fp_number_type __thenan_df;
141*c87b03e5Sespie #endif
142*c87b03e5Sespie
143*c87b03e5Sespie INLINE
144*c87b03e5Sespie static fp_number_type *
145*c87b03e5Sespie nan (void)
146*c87b03e5Sespie {
147*c87b03e5Sespie /* Discard the const qualifier... */
148*c87b03e5Sespie #ifdef TFLOAT
149*c87b03e5Sespie return (fp_number_type *) (& __thenan_tf);
150*c87b03e5Sespie #elif defined FLOAT
151*c87b03e5Sespie return (fp_number_type *) (& __thenan_sf);
152*c87b03e5Sespie #else
153*c87b03e5Sespie return (fp_number_type *) (& __thenan_df);
154*c87b03e5Sespie #endif
155*c87b03e5Sespie }
156*c87b03e5Sespie
157*c87b03e5Sespie INLINE
158*c87b03e5Sespie static int
159*c87b03e5Sespie isnan ( fp_number_type * x)
160*c87b03e5Sespie {
161*c87b03e5Sespie return x->class == CLASS_SNAN || x->class == CLASS_QNAN;
162*c87b03e5Sespie }
163*c87b03e5Sespie
164*c87b03e5Sespie INLINE
165*c87b03e5Sespie static int
166*c87b03e5Sespie isinf ( fp_number_type * x)
167*c87b03e5Sespie {
168*c87b03e5Sespie return x->class == CLASS_INFINITY;
169*c87b03e5Sespie }
170*c87b03e5Sespie
171*c87b03e5Sespie #endif /* NO_NANS */
172*c87b03e5Sespie
173*c87b03e5Sespie INLINE
174*c87b03e5Sespie static int
175*c87b03e5Sespie iszero ( fp_number_type * x)
176*c87b03e5Sespie {
177*c87b03e5Sespie return x->class == CLASS_ZERO;
178*c87b03e5Sespie }
179*c87b03e5Sespie
180*c87b03e5Sespie INLINE
181*c87b03e5Sespie static void
182*c87b03e5Sespie flip_sign ( fp_number_type * x)
183*c87b03e5Sespie {
184*c87b03e5Sespie x->sign = !x->sign;
185*c87b03e5Sespie }
186*c87b03e5Sespie
187*c87b03e5Sespie extern FLO_type pack_d ( fp_number_type * );
188*c87b03e5Sespie
189*c87b03e5Sespie #if defined(L_pack_df) || defined(L_pack_sf) || defined(L_pack_tf)
190*c87b03e5Sespie FLO_type
191*c87b03e5Sespie pack_d ( fp_number_type * src)
192*c87b03e5Sespie {
193*c87b03e5Sespie FLO_union_type dst;
194*c87b03e5Sespie fractype fraction = src->fraction.ll; /* wasn't unsigned before? */
195*c87b03e5Sespie int sign = src->sign;
196*c87b03e5Sespie int exp = 0;
197*c87b03e5Sespie
198*c87b03e5Sespie if (LARGEST_EXPONENT_IS_NORMAL (FRAC_NBITS) && (isnan (src) || isinf (src)))
199*c87b03e5Sespie {
200*c87b03e5Sespie /* We can't represent these values accurately. By using the
201*c87b03e5Sespie largest possible magnitude, we guarantee that the conversion
202*c87b03e5Sespie of infinity is at least as big as any finite number. */
203*c87b03e5Sespie exp = EXPMAX;
204*c87b03e5Sespie fraction = ((fractype) 1 << FRACBITS) - 1;
205*c87b03e5Sespie }
206*c87b03e5Sespie else if (isnan (src))
207*c87b03e5Sespie {
208*c87b03e5Sespie exp = EXPMAX;
209*c87b03e5Sespie if (src->class == CLASS_QNAN || 1)
210*c87b03e5Sespie {
211*c87b03e5Sespie fraction |= QUIET_NAN;
212*c87b03e5Sespie }
213*c87b03e5Sespie }
214*c87b03e5Sespie else if (isinf (src))
215*c87b03e5Sespie {
216*c87b03e5Sespie exp = EXPMAX;
217*c87b03e5Sespie fraction = 0;
218*c87b03e5Sespie }
219*c87b03e5Sespie else if (iszero (src))
220*c87b03e5Sespie {
221*c87b03e5Sespie exp = 0;
222*c87b03e5Sespie fraction = 0;
223*c87b03e5Sespie }
224*c87b03e5Sespie else if (fraction == 0)
225*c87b03e5Sespie {
226*c87b03e5Sespie exp = 0;
227*c87b03e5Sespie }
228*c87b03e5Sespie else
229*c87b03e5Sespie {
230*c87b03e5Sespie if (src->normal_exp < NORMAL_EXPMIN)
231*c87b03e5Sespie {
232*c87b03e5Sespie #ifdef NO_DENORMALS
233*c87b03e5Sespie /* Go straight to a zero representation if denormals are not
234*c87b03e5Sespie supported. The denormal handling would be harmless but
235*c87b03e5Sespie isn't unnecessary. */
236*c87b03e5Sespie exp = 0;
237*c87b03e5Sespie fraction = 0;
238*c87b03e5Sespie #else /* NO_DENORMALS */
239*c87b03e5Sespie /* This number's exponent is too low to fit into the bits
240*c87b03e5Sespie available in the number, so we'll store 0 in the exponent and
241*c87b03e5Sespie shift the fraction to the right to make up for it. */
242*c87b03e5Sespie
243*c87b03e5Sespie int shift = NORMAL_EXPMIN - src->normal_exp;
244*c87b03e5Sespie
245*c87b03e5Sespie exp = 0;
246*c87b03e5Sespie
247*c87b03e5Sespie if (shift > FRAC_NBITS - NGARDS)
248*c87b03e5Sespie {
249*c87b03e5Sespie /* No point shifting, since it's more that 64 out. */
250*c87b03e5Sespie fraction = 0;
251*c87b03e5Sespie }
252*c87b03e5Sespie else
253*c87b03e5Sespie {
254*c87b03e5Sespie int lowbit = (fraction & (((fractype)1 << shift) - 1)) ? 1 : 0;
255*c87b03e5Sespie fraction = (fraction >> shift) | lowbit;
256*c87b03e5Sespie }
257*c87b03e5Sespie if ((fraction & GARDMASK) == GARDMSB)
258*c87b03e5Sespie {
259*c87b03e5Sespie if ((fraction & (1 << NGARDS)))
260*c87b03e5Sespie fraction += GARDROUND + 1;
261*c87b03e5Sespie }
262*c87b03e5Sespie else
263*c87b03e5Sespie {
264*c87b03e5Sespie /* Add to the guards to round up. */
265*c87b03e5Sespie fraction += GARDROUND;
266*c87b03e5Sespie }
267*c87b03e5Sespie /* Perhaps the rounding means we now need to change the
268*c87b03e5Sespie exponent, because the fraction is no longer denormal. */
269*c87b03e5Sespie if (fraction >= IMPLICIT_1)
270*c87b03e5Sespie {
271*c87b03e5Sespie exp += 1;
272*c87b03e5Sespie }
273*c87b03e5Sespie fraction >>= NGARDS;
274*c87b03e5Sespie #endif /* NO_DENORMALS */
275*c87b03e5Sespie }
276*c87b03e5Sespie else if (!LARGEST_EXPONENT_IS_NORMAL (FRAC_NBITS)
277*c87b03e5Sespie && src->normal_exp > EXPBIAS)
278*c87b03e5Sespie {
279*c87b03e5Sespie exp = EXPMAX;
280*c87b03e5Sespie fraction = 0;
281*c87b03e5Sespie }
282*c87b03e5Sespie else
283*c87b03e5Sespie {
284*c87b03e5Sespie exp = src->normal_exp + EXPBIAS;
285*c87b03e5Sespie if (!ROUND_TOWARDS_ZERO)
286*c87b03e5Sespie {
287*c87b03e5Sespie /* IF the gard bits are the all zero, but the first, then we're
288*c87b03e5Sespie half way between two numbers, choose the one which makes the
289*c87b03e5Sespie lsb of the answer 0. */
290*c87b03e5Sespie if ((fraction & GARDMASK) == GARDMSB)
291*c87b03e5Sespie {
292*c87b03e5Sespie if (fraction & (1 << NGARDS))
293*c87b03e5Sespie fraction += GARDROUND + 1;
294*c87b03e5Sespie }
295*c87b03e5Sespie else
296*c87b03e5Sespie {
297*c87b03e5Sespie /* Add a one to the guards to round up */
298*c87b03e5Sespie fraction += GARDROUND;
299*c87b03e5Sespie }
300*c87b03e5Sespie if (fraction >= IMPLICIT_2)
301*c87b03e5Sespie {
302*c87b03e5Sespie fraction >>= 1;
303*c87b03e5Sespie exp += 1;
304*c87b03e5Sespie }
305*c87b03e5Sespie }
306*c87b03e5Sespie fraction >>= NGARDS;
307*c87b03e5Sespie
308*c87b03e5Sespie if (LARGEST_EXPONENT_IS_NORMAL (FRAC_NBITS) && exp > EXPMAX)
309*c87b03e5Sespie {
310*c87b03e5Sespie /* Saturate on overflow. */
311*c87b03e5Sespie exp = EXPMAX;
312*c87b03e5Sespie fraction = ((fractype) 1 << FRACBITS) - 1;
313*c87b03e5Sespie }
314*c87b03e5Sespie }
315*c87b03e5Sespie }
316*c87b03e5Sespie
317*c87b03e5Sespie /* We previously used bitfields to store the number, but this doesn't
318*c87b03e5Sespie handle little/big endian systems conveniently, so use shifts and
319*c87b03e5Sespie masks */
320*c87b03e5Sespie #ifdef FLOAT_BIT_ORDER_MISMATCH
321*c87b03e5Sespie dst.bits.fraction = fraction;
322*c87b03e5Sespie dst.bits.exp = exp;
323*c87b03e5Sespie dst.bits.sign = sign;
324*c87b03e5Sespie #else
325*c87b03e5Sespie # if defined TFLOAT && defined HALFFRACBITS
326*c87b03e5Sespie {
327*c87b03e5Sespie halffractype high, low;
328*c87b03e5Sespie
329*c87b03e5Sespie high = (fraction >> (FRACBITS - HALFFRACBITS));
330*c87b03e5Sespie high &= (((fractype)1) << HALFFRACBITS) - 1;
331*c87b03e5Sespie high |= ((fractype) (exp & ((1 << EXPBITS) - 1))) << HALFFRACBITS;
332*c87b03e5Sespie high |= ((fractype) (sign & 1)) << (HALFFRACBITS | EXPBITS);
333*c87b03e5Sespie
334*c87b03e5Sespie low = (halffractype)fraction &
335*c87b03e5Sespie ((((halffractype)1) << (FRACBITS - HALFFRACBITS)) - 1);
336*c87b03e5Sespie
337*c87b03e5Sespie if (exp == EXPMAX || exp == 0 || low == 0)
338*c87b03e5Sespie low = 0;
339*c87b03e5Sespie else
340*c87b03e5Sespie {
341*c87b03e5Sespie exp -= HALFFRACBITS + 1;
342*c87b03e5Sespie
343*c87b03e5Sespie while (exp > 0
344*c87b03e5Sespie && low < ((halffractype)1 << HALFFRACBITS))
345*c87b03e5Sespie {
346*c87b03e5Sespie low <<= 1;
347*c87b03e5Sespie exp--;
348*c87b03e5Sespie }
349*c87b03e5Sespie
350*c87b03e5Sespie if (exp <= 0)
351*c87b03e5Sespie {
352*c87b03e5Sespie halffractype roundmsb, round;
353*c87b03e5Sespie
354*c87b03e5Sespie exp = -exp + 1;
355*c87b03e5Sespie
356*c87b03e5Sespie roundmsb = (1 << (exp - 1));
357*c87b03e5Sespie round = low & ((roundmsb << 1) - 1);
358*c87b03e5Sespie
359*c87b03e5Sespie low >>= exp;
360*c87b03e5Sespie exp = 0;
361*c87b03e5Sespie
362*c87b03e5Sespie if (round > roundmsb || (round == roundmsb && (low & 1)))
363*c87b03e5Sespie {
364*c87b03e5Sespie low++;
365*c87b03e5Sespie if (low >= ((halffractype)1 << HALFFRACBITS))
366*c87b03e5Sespie /* We don't shift left, since it has just become the
367*c87b03e5Sespie smallest normal number, whose implicit 1 bit is
368*c87b03e5Sespie now indicated by the non-zero exponent. */
369*c87b03e5Sespie exp++;
370*c87b03e5Sespie }
371*c87b03e5Sespie }
372*c87b03e5Sespie
373*c87b03e5Sespie low &= ((halffractype)1 << HALFFRACBITS) - 1;
374*c87b03e5Sespie low |= ((fractype) (exp & ((1 << EXPBITS) - 1))) << HALFFRACBITS;
375*c87b03e5Sespie low |= ((fractype) (sign & 1)) << (HALFFRACBITS | EXPBITS);
376*c87b03e5Sespie }
377*c87b03e5Sespie
378*c87b03e5Sespie dst.value_raw = (((fractype) high) << HALFSHIFT) | low;
379*c87b03e5Sespie }
380*c87b03e5Sespie # else
381*c87b03e5Sespie dst.value_raw = fraction & ((((fractype)1) << FRACBITS) - (fractype)1);
382*c87b03e5Sespie dst.value_raw |= ((fractype) (exp & ((1 << EXPBITS) - 1))) << FRACBITS;
383*c87b03e5Sespie dst.value_raw |= ((fractype) (sign & 1)) << (FRACBITS | EXPBITS);
384*c87b03e5Sespie # endif
385*c87b03e5Sespie #endif
386*c87b03e5Sespie
387*c87b03e5Sespie #if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT)
388*c87b03e5Sespie #ifdef TFLOAT
389*c87b03e5Sespie {
390*c87b03e5Sespie qrtrfractype tmp1 = dst.words[0];
391*c87b03e5Sespie qrtrfractype tmp2 = dst.words[1];
392*c87b03e5Sespie dst.words[0] = dst.words[3];
393*c87b03e5Sespie dst.words[1] = dst.words[2];
394*c87b03e5Sespie dst.words[2] = tmp2;
395*c87b03e5Sespie dst.words[3] = tmp1;
396*c87b03e5Sespie }
397*c87b03e5Sespie #else
398*c87b03e5Sespie {
399*c87b03e5Sespie halffractype tmp = dst.words[0];
400*c87b03e5Sespie dst.words[0] = dst.words[1];
401*c87b03e5Sespie dst.words[1] = tmp;
402*c87b03e5Sespie }
403*c87b03e5Sespie #endif
404*c87b03e5Sespie #endif
405*c87b03e5Sespie
406*c87b03e5Sespie return dst.value;
407*c87b03e5Sespie }
408*c87b03e5Sespie #endif
409*c87b03e5Sespie
410*c87b03e5Sespie #if defined(L_unpack_df) || defined(L_unpack_sf) || defined(L_unpack_tf)
411*c87b03e5Sespie void
412*c87b03e5Sespie unpack_d (FLO_union_type * src, fp_number_type * dst)
413*c87b03e5Sespie {
414*c87b03e5Sespie /* We previously used bitfields to store the number, but this doesn't
415*c87b03e5Sespie handle little/big endian systems conveniently, so use shifts and
416*c87b03e5Sespie masks */
417*c87b03e5Sespie fractype fraction;
418*c87b03e5Sespie int exp;
419*c87b03e5Sespie int sign;
420*c87b03e5Sespie
421*c87b03e5Sespie #if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT)
422*c87b03e5Sespie FLO_union_type swapped;
423*c87b03e5Sespie
424*c87b03e5Sespie #ifdef TFLOAT
425*c87b03e5Sespie swapped.words[0] = src->words[3];
426*c87b03e5Sespie swapped.words[1] = src->words[2];
427*c87b03e5Sespie swapped.words[2] = src->words[1];
428*c87b03e5Sespie swapped.words[3] = src->words[0];
429*c87b03e5Sespie #else
430*c87b03e5Sespie swapped.words[0] = src->words[1];
431*c87b03e5Sespie swapped.words[1] = src->words[0];
432*c87b03e5Sespie #endif
433*c87b03e5Sespie src = &swapped;
434*c87b03e5Sespie #endif
435*c87b03e5Sespie
436*c87b03e5Sespie #ifdef FLOAT_BIT_ORDER_MISMATCH
437*c87b03e5Sespie fraction = src->bits.fraction;
438*c87b03e5Sespie exp = src->bits.exp;
439*c87b03e5Sespie sign = src->bits.sign;
440*c87b03e5Sespie #else
441*c87b03e5Sespie # if defined TFLOAT && defined HALFFRACBITS
442*c87b03e5Sespie {
443*c87b03e5Sespie halffractype high, low;
444*c87b03e5Sespie
445*c87b03e5Sespie high = src->value_raw >> HALFSHIFT;
446*c87b03e5Sespie low = src->value_raw & (((fractype)1 << HALFSHIFT) - 1);
447*c87b03e5Sespie
448*c87b03e5Sespie fraction = high & ((((fractype)1) << HALFFRACBITS) - 1);
449*c87b03e5Sespie fraction <<= FRACBITS - HALFFRACBITS;
450*c87b03e5Sespie exp = ((int)(high >> HALFFRACBITS)) & ((1 << EXPBITS) - 1);
451*c87b03e5Sespie sign = ((int)(high >> (((HALFFRACBITS + EXPBITS))))) & 1;
452*c87b03e5Sespie
453*c87b03e5Sespie if (exp != EXPMAX && exp != 0 && low != 0)
454*c87b03e5Sespie {
455*c87b03e5Sespie int lowexp = ((int)(low >> HALFFRACBITS)) & ((1 << EXPBITS) - 1);
456*c87b03e5Sespie int lowsign = ((int)(low >> (((HALFFRACBITS + EXPBITS))))) & 1;
457*c87b03e5Sespie int shift;
458*c87b03e5Sespie fractype xlow;
459*c87b03e5Sespie
460*c87b03e5Sespie xlow = low & ((((fractype)1) << HALFFRACBITS) - 1);
461*c87b03e5Sespie if (lowexp)
462*c87b03e5Sespie xlow |= (((halffractype)1) << HALFFRACBITS);
463*c87b03e5Sespie else
464*c87b03e5Sespie lowexp = 1;
465*c87b03e5Sespie shift = (FRACBITS - HALFFRACBITS) - (exp - lowexp);
466*c87b03e5Sespie if (shift > 0)
467*c87b03e5Sespie xlow <<= shift;
468*c87b03e5Sespie else if (shift < 0)
469*c87b03e5Sespie xlow >>= -shift;
470*c87b03e5Sespie if (sign == lowsign)
471*c87b03e5Sespie fraction += xlow;
472*c87b03e5Sespie else
473*c87b03e5Sespie fraction -= xlow;
474*c87b03e5Sespie }
475*c87b03e5Sespie }
476*c87b03e5Sespie # else
477*c87b03e5Sespie fraction = src->value_raw & ((((fractype)1) << FRACBITS) - 1);
478*c87b03e5Sespie exp = ((int)(src->value_raw >> FRACBITS)) & ((1 << EXPBITS) - 1);
479*c87b03e5Sespie sign = ((int)(src->value_raw >> (FRACBITS + EXPBITS))) & 1;
480*c87b03e5Sespie # endif
481*c87b03e5Sespie #endif
482*c87b03e5Sespie
483*c87b03e5Sespie dst->sign = sign;
484*c87b03e5Sespie if (exp == 0)
485*c87b03e5Sespie {
486*c87b03e5Sespie /* Hmm. Looks like 0 */
487*c87b03e5Sespie if (fraction == 0
488*c87b03e5Sespie #ifdef NO_DENORMALS
489*c87b03e5Sespie || 1
490*c87b03e5Sespie #endif
491*c87b03e5Sespie )
492*c87b03e5Sespie {
493*c87b03e5Sespie /* tastes like zero */
494*c87b03e5Sespie dst->class = CLASS_ZERO;
495*c87b03e5Sespie }
496*c87b03e5Sespie else
497*c87b03e5Sespie {
498*c87b03e5Sespie /* Zero exponent with nonzero fraction - it's denormalized,
499*c87b03e5Sespie so there isn't a leading implicit one - we'll shift it so
500*c87b03e5Sespie it gets one. */
501*c87b03e5Sespie dst->normal_exp = exp - EXPBIAS + 1;
502*c87b03e5Sespie fraction <<= NGARDS;
503*c87b03e5Sespie
504*c87b03e5Sespie dst->class = CLASS_NUMBER;
505*c87b03e5Sespie #if 1
506*c87b03e5Sespie while (fraction < IMPLICIT_1)
507*c87b03e5Sespie {
508*c87b03e5Sespie fraction <<= 1;
509*c87b03e5Sespie dst->normal_exp--;
510*c87b03e5Sespie }
511*c87b03e5Sespie #endif
512*c87b03e5Sespie dst->fraction.ll = fraction;
513*c87b03e5Sespie }
514*c87b03e5Sespie }
515*c87b03e5Sespie else if (!LARGEST_EXPONENT_IS_NORMAL (FRAC_NBITS) && exp == EXPMAX)
516*c87b03e5Sespie {
517*c87b03e5Sespie /* Huge exponent*/
518*c87b03e5Sespie if (fraction == 0)
519*c87b03e5Sespie {
520*c87b03e5Sespie /* Attached to a zero fraction - means infinity */
521*c87b03e5Sespie dst->class = CLASS_INFINITY;
522*c87b03e5Sespie }
523*c87b03e5Sespie else
524*c87b03e5Sespie {
525*c87b03e5Sespie /* Nonzero fraction, means nan */
526*c87b03e5Sespie if (fraction & QUIET_NAN)
527*c87b03e5Sespie {
528*c87b03e5Sespie dst->class = CLASS_QNAN;
529*c87b03e5Sespie }
530*c87b03e5Sespie else
531*c87b03e5Sespie {
532*c87b03e5Sespie dst->class = CLASS_SNAN;
533*c87b03e5Sespie }
534*c87b03e5Sespie /* Keep the fraction part as the nan number */
535*c87b03e5Sespie dst->fraction.ll = fraction;
536*c87b03e5Sespie }
537*c87b03e5Sespie }
538*c87b03e5Sespie else
539*c87b03e5Sespie {
540*c87b03e5Sespie /* Nothing strange about this number */
541*c87b03e5Sespie dst->normal_exp = exp - EXPBIAS;
542*c87b03e5Sespie dst->class = CLASS_NUMBER;
543*c87b03e5Sespie dst->fraction.ll = (fraction << NGARDS) | IMPLICIT_1;
544*c87b03e5Sespie }
545*c87b03e5Sespie }
546*c87b03e5Sespie #endif /* L_unpack_df || L_unpack_sf */
547*c87b03e5Sespie
548*c87b03e5Sespie #if defined(L_addsub_sf) || defined(L_addsub_df) || defined(L_addsub_tf)
549*c87b03e5Sespie static fp_number_type *
550*c87b03e5Sespie _fpadd_parts (fp_number_type * a,
551*c87b03e5Sespie fp_number_type * b,
552*c87b03e5Sespie fp_number_type * tmp)
553*c87b03e5Sespie {
554*c87b03e5Sespie intfrac tfraction;
555*c87b03e5Sespie
556*c87b03e5Sespie /* Put commonly used fields in local variables. */
557*c87b03e5Sespie int a_normal_exp;
558*c87b03e5Sespie int b_normal_exp;
559*c87b03e5Sespie fractype a_fraction;
560*c87b03e5Sespie fractype b_fraction;
561*c87b03e5Sespie
562*c87b03e5Sespie if (isnan (a))
563*c87b03e5Sespie {
564*c87b03e5Sespie return a;
565*c87b03e5Sespie }
566*c87b03e5Sespie if (isnan (b))
567*c87b03e5Sespie {
568*c87b03e5Sespie return b;
569*c87b03e5Sespie }
570*c87b03e5Sespie if (isinf (a))
571*c87b03e5Sespie {
572*c87b03e5Sespie /* Adding infinities with opposite signs yields a NaN. */
573*c87b03e5Sespie if (isinf (b) && a->sign != b->sign)
574*c87b03e5Sespie return nan ();
575*c87b03e5Sespie return a;
576*c87b03e5Sespie }
577*c87b03e5Sespie if (isinf (b))
578*c87b03e5Sespie {
579*c87b03e5Sespie return b;
580*c87b03e5Sespie }
581*c87b03e5Sespie if (iszero (b))
582*c87b03e5Sespie {
583*c87b03e5Sespie if (iszero (a))
584*c87b03e5Sespie {
585*c87b03e5Sespie *tmp = *a;
586*c87b03e5Sespie tmp->sign = a->sign & b->sign;
587*c87b03e5Sespie return tmp;
588*c87b03e5Sespie }
589*c87b03e5Sespie return a;
590*c87b03e5Sespie }
591*c87b03e5Sespie if (iszero (a))
592*c87b03e5Sespie {
593*c87b03e5Sespie return b;
594*c87b03e5Sespie }
595*c87b03e5Sespie
596*c87b03e5Sespie /* Got two numbers. shift the smaller and increment the exponent till
597*c87b03e5Sespie they're the same */
598*c87b03e5Sespie {
599*c87b03e5Sespie int diff;
600*c87b03e5Sespie
601*c87b03e5Sespie a_normal_exp = a->normal_exp;
602*c87b03e5Sespie b_normal_exp = b->normal_exp;
603*c87b03e5Sespie a_fraction = a->fraction.ll;
604*c87b03e5Sespie b_fraction = b->fraction.ll;
605*c87b03e5Sespie
606*c87b03e5Sespie diff = a_normal_exp - b_normal_exp;
607*c87b03e5Sespie
608*c87b03e5Sespie if (diff < 0)
609*c87b03e5Sespie diff = -diff;
610*c87b03e5Sespie if (diff < FRAC_NBITS)
611*c87b03e5Sespie {
612*c87b03e5Sespie /* ??? This does shifts one bit at a time. Optimize. */
613*c87b03e5Sespie while (a_normal_exp > b_normal_exp)
614*c87b03e5Sespie {
615*c87b03e5Sespie b_normal_exp++;
616*c87b03e5Sespie LSHIFT (b_fraction);
617*c87b03e5Sespie }
618*c87b03e5Sespie while (b_normal_exp > a_normal_exp)
619*c87b03e5Sespie {
620*c87b03e5Sespie a_normal_exp++;
621*c87b03e5Sespie LSHIFT (a_fraction);
622*c87b03e5Sespie }
623*c87b03e5Sespie }
624*c87b03e5Sespie else
625*c87b03e5Sespie {
626*c87b03e5Sespie /* Somethings's up.. choose the biggest */
627*c87b03e5Sespie if (a_normal_exp > b_normal_exp)
628*c87b03e5Sespie {
629*c87b03e5Sespie b_normal_exp = a_normal_exp;
630*c87b03e5Sespie b_fraction = 0;
631*c87b03e5Sespie }
632*c87b03e5Sespie else
633*c87b03e5Sespie {
634*c87b03e5Sespie a_normal_exp = b_normal_exp;
635*c87b03e5Sespie a_fraction = 0;
636*c87b03e5Sespie }
637*c87b03e5Sespie }
638*c87b03e5Sespie }
639*c87b03e5Sespie
640*c87b03e5Sespie if (a->sign != b->sign)
641*c87b03e5Sespie {
642*c87b03e5Sespie if (a->sign)
643*c87b03e5Sespie {
644*c87b03e5Sespie tfraction = -a_fraction + b_fraction;
645*c87b03e5Sespie }
646*c87b03e5Sespie else
647*c87b03e5Sespie {
648*c87b03e5Sespie tfraction = a_fraction - b_fraction;
649*c87b03e5Sespie }
650*c87b03e5Sespie if (tfraction >= 0)
651*c87b03e5Sespie {
652*c87b03e5Sespie tmp->sign = 0;
653*c87b03e5Sespie tmp->normal_exp = a_normal_exp;
654*c87b03e5Sespie tmp->fraction.ll = tfraction;
655*c87b03e5Sespie }
656*c87b03e5Sespie else
657*c87b03e5Sespie {
658*c87b03e5Sespie tmp->sign = 1;
659*c87b03e5Sespie tmp->normal_exp = a_normal_exp;
660*c87b03e5Sespie tmp->fraction.ll = -tfraction;
661*c87b03e5Sespie }
662*c87b03e5Sespie /* and renormalize it */
663*c87b03e5Sespie
664*c87b03e5Sespie while (tmp->fraction.ll < IMPLICIT_1 && tmp->fraction.ll)
665*c87b03e5Sespie {
666*c87b03e5Sespie tmp->fraction.ll <<= 1;
667*c87b03e5Sespie tmp->normal_exp--;
668*c87b03e5Sespie }
669*c87b03e5Sespie }
670*c87b03e5Sespie else
671*c87b03e5Sespie {
672*c87b03e5Sespie tmp->sign = a->sign;
673*c87b03e5Sespie tmp->normal_exp = a_normal_exp;
674*c87b03e5Sespie tmp->fraction.ll = a_fraction + b_fraction;
675*c87b03e5Sespie }
676*c87b03e5Sespie tmp->class = CLASS_NUMBER;
677*c87b03e5Sespie /* Now the fraction is added, we have to shift down to renormalize the
678*c87b03e5Sespie number */
679*c87b03e5Sespie
680*c87b03e5Sespie if (tmp->fraction.ll >= IMPLICIT_2)
681*c87b03e5Sespie {
682*c87b03e5Sespie LSHIFT (tmp->fraction.ll);
683*c87b03e5Sespie tmp->normal_exp++;
684*c87b03e5Sespie }
685*c87b03e5Sespie return tmp;
686*c87b03e5Sespie
687*c87b03e5Sespie }
688*c87b03e5Sespie
689*c87b03e5Sespie FLO_type
690*c87b03e5Sespie add (FLO_type arg_a, FLO_type arg_b)
691*c87b03e5Sespie {
692*c87b03e5Sespie fp_number_type a;
693*c87b03e5Sespie fp_number_type b;
694*c87b03e5Sespie fp_number_type tmp;
695*c87b03e5Sespie fp_number_type *res;
696*c87b03e5Sespie FLO_union_type au, bu;
697*c87b03e5Sespie
698*c87b03e5Sespie au.value = arg_a;
699*c87b03e5Sespie bu.value = arg_b;
700*c87b03e5Sespie
701*c87b03e5Sespie unpack_d (&au, &a);
702*c87b03e5Sespie unpack_d (&bu, &b);
703*c87b03e5Sespie
704*c87b03e5Sespie res = _fpadd_parts (&a, &b, &tmp);
705*c87b03e5Sespie
706*c87b03e5Sespie return pack_d (res);
707*c87b03e5Sespie }
708*c87b03e5Sespie
709*c87b03e5Sespie FLO_type
710*c87b03e5Sespie sub (FLO_type arg_a, FLO_type arg_b)
711*c87b03e5Sespie {
712*c87b03e5Sespie fp_number_type a;
713*c87b03e5Sespie fp_number_type b;
714*c87b03e5Sespie fp_number_type tmp;
715*c87b03e5Sespie fp_number_type *res;
716*c87b03e5Sespie FLO_union_type au, bu;
717*c87b03e5Sespie
718*c87b03e5Sespie au.value = arg_a;
719*c87b03e5Sespie bu.value = arg_b;
720*c87b03e5Sespie
721*c87b03e5Sespie unpack_d (&au, &a);
722*c87b03e5Sespie unpack_d (&bu, &b);
723*c87b03e5Sespie
724*c87b03e5Sespie b.sign ^= 1;
725*c87b03e5Sespie
726*c87b03e5Sespie res = _fpadd_parts (&a, &b, &tmp);
727*c87b03e5Sespie
728*c87b03e5Sespie return pack_d (res);
729*c87b03e5Sespie }
730*c87b03e5Sespie #endif /* L_addsub_sf || L_addsub_df */
731*c87b03e5Sespie
732*c87b03e5Sespie #if defined(L_mul_sf) || defined(L_mul_df) || defined(L_mul_tf)
733*c87b03e5Sespie static inline __attribute__ ((__always_inline__)) fp_number_type *
734*c87b03e5Sespie _fpmul_parts ( fp_number_type * a,
735*c87b03e5Sespie fp_number_type * b,
736*c87b03e5Sespie fp_number_type * tmp)
737*c87b03e5Sespie {
738*c87b03e5Sespie fractype low = 0;
739*c87b03e5Sespie fractype high = 0;
740*c87b03e5Sespie
741*c87b03e5Sespie if (isnan (a))
742*c87b03e5Sespie {
743*c87b03e5Sespie a->sign = a->sign != b->sign;
744*c87b03e5Sespie return a;
745*c87b03e5Sespie }
746*c87b03e5Sespie if (isnan (b))
747*c87b03e5Sespie {
748*c87b03e5Sespie b->sign = a->sign != b->sign;
749*c87b03e5Sespie return b;
750*c87b03e5Sespie }
751*c87b03e5Sespie if (isinf (a))
752*c87b03e5Sespie {
753*c87b03e5Sespie if (iszero (b))
754*c87b03e5Sespie return nan ();
755*c87b03e5Sespie a->sign = a->sign != b->sign;
756*c87b03e5Sespie return a;
757*c87b03e5Sespie }
758*c87b03e5Sespie if (isinf (b))
759*c87b03e5Sespie {
760*c87b03e5Sespie if (iszero (a))
761*c87b03e5Sespie {
762*c87b03e5Sespie return nan ();
763*c87b03e5Sespie }
764*c87b03e5Sespie b->sign = a->sign != b->sign;
765*c87b03e5Sespie return b;
766*c87b03e5Sespie }
767*c87b03e5Sespie if (iszero (a))
768*c87b03e5Sespie {
769*c87b03e5Sespie a->sign = a->sign != b->sign;
770*c87b03e5Sespie return a;
771*c87b03e5Sespie }
772*c87b03e5Sespie if (iszero (b))
773*c87b03e5Sespie {
774*c87b03e5Sespie b->sign = a->sign != b->sign;
775*c87b03e5Sespie return b;
776*c87b03e5Sespie }
777*c87b03e5Sespie
778*c87b03e5Sespie /* Calculate the mantissa by multiplying both numbers to get a
779*c87b03e5Sespie twice-as-wide number. */
780*c87b03e5Sespie {
781*c87b03e5Sespie #if defined(NO_DI_MODE) || defined(TFLOAT)
782*c87b03e5Sespie {
783*c87b03e5Sespie fractype x = a->fraction.ll;
784*c87b03e5Sespie fractype ylow = b->fraction.ll;
785*c87b03e5Sespie fractype yhigh = 0;
786*c87b03e5Sespie int bit;
787*c87b03e5Sespie
788*c87b03e5Sespie /* ??? This does multiplies one bit at a time. Optimize. */
789*c87b03e5Sespie for (bit = 0; bit < FRAC_NBITS; bit++)
790*c87b03e5Sespie {
791*c87b03e5Sespie int carry;
792*c87b03e5Sespie
793*c87b03e5Sespie if (x & 1)
794*c87b03e5Sespie {
795*c87b03e5Sespie carry = (low += ylow) < ylow;
796*c87b03e5Sespie high += yhigh + carry;
797*c87b03e5Sespie }
798*c87b03e5Sespie yhigh <<= 1;
799*c87b03e5Sespie if (ylow & FRACHIGH)
800*c87b03e5Sespie {
801*c87b03e5Sespie yhigh |= 1;
802*c87b03e5Sespie }
803*c87b03e5Sespie ylow <<= 1;
804*c87b03e5Sespie x >>= 1;
805*c87b03e5Sespie }
806*c87b03e5Sespie }
807*c87b03e5Sespie #elif defined(FLOAT)
808*c87b03e5Sespie /* Multiplying two USIs to get a UDI, we're safe. */
809*c87b03e5Sespie {
810*c87b03e5Sespie UDItype answer = (UDItype)a->fraction.ll * (UDItype)b->fraction.ll;
811*c87b03e5Sespie
812*c87b03e5Sespie high = answer >> BITS_PER_SI;
813*c87b03e5Sespie low = answer;
814*c87b03e5Sespie }
815*c87b03e5Sespie #else
816*c87b03e5Sespie /* fractype is DImode, but we need the result to be twice as wide.
817*c87b03e5Sespie Assuming a widening multiply from DImode to TImode is not
818*c87b03e5Sespie available, build one by hand. */
819*c87b03e5Sespie {
820*c87b03e5Sespie USItype nl = a->fraction.ll;
821*c87b03e5Sespie USItype nh = a->fraction.ll >> BITS_PER_SI;
822*c87b03e5Sespie USItype ml = b->fraction.ll;
823*c87b03e5Sespie USItype mh = b->fraction.ll >> BITS_PER_SI;
824*c87b03e5Sespie UDItype pp_ll = (UDItype) ml * nl;
825*c87b03e5Sespie UDItype pp_hl = (UDItype) mh * nl;
826*c87b03e5Sespie UDItype pp_lh = (UDItype) ml * nh;
827*c87b03e5Sespie UDItype pp_hh = (UDItype) mh * nh;
828*c87b03e5Sespie UDItype res2 = 0;
829*c87b03e5Sespie UDItype res0 = 0;
830*c87b03e5Sespie UDItype ps_hh__ = pp_hl + pp_lh;
831*c87b03e5Sespie if (ps_hh__ < pp_hl)
832*c87b03e5Sespie res2 += (UDItype)1 << BITS_PER_SI;
833*c87b03e5Sespie pp_hl = (UDItype)(USItype)ps_hh__ << BITS_PER_SI;
834*c87b03e5Sespie res0 = pp_ll + pp_hl;
835*c87b03e5Sespie if (res0 < pp_ll)
836*c87b03e5Sespie res2++;
837*c87b03e5Sespie res2 += (ps_hh__ >> BITS_PER_SI) + pp_hh;
838*c87b03e5Sespie high = res2;
839*c87b03e5Sespie low = res0;
840*c87b03e5Sespie }
841*c87b03e5Sespie #endif
842*c87b03e5Sespie }
843*c87b03e5Sespie
844*c87b03e5Sespie tmp->normal_exp = a->normal_exp + b->normal_exp
845*c87b03e5Sespie + FRAC_NBITS - (FRACBITS + NGARDS);
846*c87b03e5Sespie tmp->sign = a->sign != b->sign;
847*c87b03e5Sespie while (high >= IMPLICIT_2)
848*c87b03e5Sespie {
849*c87b03e5Sespie tmp->normal_exp++;
850*c87b03e5Sespie if (high & 1)
851*c87b03e5Sespie {
852*c87b03e5Sespie low >>= 1;
853*c87b03e5Sespie low |= FRACHIGH;
854*c87b03e5Sespie }
855*c87b03e5Sespie high >>= 1;
856*c87b03e5Sespie }
857*c87b03e5Sespie while (high < IMPLICIT_1)
858*c87b03e5Sespie {
859*c87b03e5Sespie tmp->normal_exp--;
860*c87b03e5Sespie
861*c87b03e5Sespie high <<= 1;
862*c87b03e5Sespie if (low & FRACHIGH)
863*c87b03e5Sespie high |= 1;
864*c87b03e5Sespie low <<= 1;
865*c87b03e5Sespie }
866*c87b03e5Sespie /* rounding is tricky. if we only round if it won't make us round later. */
867*c87b03e5Sespie #if 0
868*c87b03e5Sespie if (low & FRACHIGH2)
869*c87b03e5Sespie {
870*c87b03e5Sespie if (((high & GARDMASK) != GARDMSB)
871*c87b03e5Sespie && (((high + 1) & GARDMASK) == GARDMSB))
872*c87b03e5Sespie {
873*c87b03e5Sespie /* don't round, it gets done again later. */
874*c87b03e5Sespie }
875*c87b03e5Sespie else
876*c87b03e5Sespie {
877*c87b03e5Sespie high++;
878*c87b03e5Sespie }
879*c87b03e5Sespie }
880*c87b03e5Sespie #endif
881*c87b03e5Sespie if (!ROUND_TOWARDS_ZERO && (high & GARDMASK) == GARDMSB)
882*c87b03e5Sespie {
883*c87b03e5Sespie if (high & (1 << NGARDS))
884*c87b03e5Sespie {
885*c87b03e5Sespie /* half way, so round to even */
886*c87b03e5Sespie high += GARDROUND + 1;
887*c87b03e5Sespie }
888*c87b03e5Sespie else if (low)
889*c87b03e5Sespie {
890*c87b03e5Sespie /* but we really weren't half way */
891*c87b03e5Sespie high += GARDROUND + 1;
892*c87b03e5Sespie }
893*c87b03e5Sespie }
894*c87b03e5Sespie tmp->fraction.ll = high;
895*c87b03e5Sespie tmp->class = CLASS_NUMBER;
896*c87b03e5Sespie return tmp;
897*c87b03e5Sespie }
898*c87b03e5Sespie
899*c87b03e5Sespie FLO_type
900*c87b03e5Sespie multiply (FLO_type arg_a, FLO_type arg_b)
901*c87b03e5Sespie {
902*c87b03e5Sespie fp_number_type a;
903*c87b03e5Sespie fp_number_type b;
904*c87b03e5Sespie fp_number_type tmp;
905*c87b03e5Sespie fp_number_type *res;
906*c87b03e5Sespie FLO_union_type au, bu;
907*c87b03e5Sespie
908*c87b03e5Sespie au.value = arg_a;
909*c87b03e5Sespie bu.value = arg_b;
910*c87b03e5Sespie
911*c87b03e5Sespie unpack_d (&au, &a);
912*c87b03e5Sespie unpack_d (&bu, &b);
913*c87b03e5Sespie
914*c87b03e5Sespie res = _fpmul_parts (&a, &b, &tmp);
915*c87b03e5Sespie
916*c87b03e5Sespie return pack_d (res);
917*c87b03e5Sespie }
918*c87b03e5Sespie #endif /* L_mul_sf || L_mul_df */
919*c87b03e5Sespie
920*c87b03e5Sespie #if defined(L_div_sf) || defined(L_div_df) || defined(L_div_tf)
921*c87b03e5Sespie static inline __attribute__ ((__always_inline__)) fp_number_type *
922*c87b03e5Sespie _fpdiv_parts (fp_number_type * a,
923*c87b03e5Sespie fp_number_type * b)
924*c87b03e5Sespie {
925*c87b03e5Sespie fractype bit;
926*c87b03e5Sespie fractype numerator;
927*c87b03e5Sespie fractype denominator;
928*c87b03e5Sespie fractype quotient;
929*c87b03e5Sespie
930*c87b03e5Sespie if (isnan (a))
931*c87b03e5Sespie {
932*c87b03e5Sespie return a;
933*c87b03e5Sespie }
934*c87b03e5Sespie if (isnan (b))
935*c87b03e5Sespie {
936*c87b03e5Sespie return b;
937*c87b03e5Sespie }
938*c87b03e5Sespie
939*c87b03e5Sespie a->sign = a->sign ^ b->sign;
940*c87b03e5Sespie
941*c87b03e5Sespie if (isinf (a) || iszero (a))
942*c87b03e5Sespie {
943*c87b03e5Sespie if (a->class == b->class)
944*c87b03e5Sespie return nan ();
945*c87b03e5Sespie return a;
946*c87b03e5Sespie }
947*c87b03e5Sespie
948*c87b03e5Sespie if (isinf (b))
949*c87b03e5Sespie {
950*c87b03e5Sespie a->fraction.ll = 0;
951*c87b03e5Sespie a->normal_exp = 0;
952*c87b03e5Sespie return a;
953*c87b03e5Sespie }
954*c87b03e5Sespie if (iszero (b))
955*c87b03e5Sespie {
956*c87b03e5Sespie a->class = CLASS_INFINITY;
957*c87b03e5Sespie return a;
958*c87b03e5Sespie }
959*c87b03e5Sespie
960*c87b03e5Sespie /* Calculate the mantissa by multiplying both 64bit numbers to get a
961*c87b03e5Sespie 128 bit number */
962*c87b03e5Sespie {
963*c87b03e5Sespie /* quotient =
964*c87b03e5Sespie ( numerator / denominator) * 2^(numerator exponent - denominator exponent)
965*c87b03e5Sespie */
966*c87b03e5Sespie
967*c87b03e5Sespie a->normal_exp = a->normal_exp - b->normal_exp;
968*c87b03e5Sespie numerator = a->fraction.ll;
969*c87b03e5Sespie denominator = b->fraction.ll;
970*c87b03e5Sespie
971*c87b03e5Sespie if (numerator < denominator)
972*c87b03e5Sespie {
973*c87b03e5Sespie /* Fraction will be less than 1.0 */
974*c87b03e5Sespie numerator *= 2;
975*c87b03e5Sespie a->normal_exp--;
976*c87b03e5Sespie }
977*c87b03e5Sespie bit = IMPLICIT_1;
978*c87b03e5Sespie quotient = 0;
979*c87b03e5Sespie /* ??? Does divide one bit at a time. Optimize. */
980*c87b03e5Sespie while (bit)
981*c87b03e5Sespie {
982*c87b03e5Sespie if (numerator >= denominator)
983*c87b03e5Sespie {
984*c87b03e5Sespie quotient |= bit;
985*c87b03e5Sespie numerator -= denominator;
986*c87b03e5Sespie }
987*c87b03e5Sespie bit >>= 1;
988*c87b03e5Sespie numerator *= 2;
989*c87b03e5Sespie }
990*c87b03e5Sespie
991*c87b03e5Sespie if (!ROUND_TOWARDS_ZERO && (quotient & GARDMASK) == GARDMSB)
992*c87b03e5Sespie {
993*c87b03e5Sespie if (quotient & (1 << NGARDS))
994*c87b03e5Sespie {
995*c87b03e5Sespie /* half way, so round to even */
996*c87b03e5Sespie quotient += GARDROUND + 1;
997*c87b03e5Sespie }
998*c87b03e5Sespie else if (numerator)
999*c87b03e5Sespie {
1000*c87b03e5Sespie /* but we really weren't half way, more bits exist */
1001*c87b03e5Sespie quotient += GARDROUND + 1;
1002*c87b03e5Sespie }
1003*c87b03e5Sespie }
1004*c87b03e5Sespie
1005*c87b03e5Sespie a->fraction.ll = quotient;
1006*c87b03e5Sespie return (a);
1007*c87b03e5Sespie }
1008*c87b03e5Sespie }
1009*c87b03e5Sespie
1010*c87b03e5Sespie FLO_type
1011*c87b03e5Sespie divide (FLO_type arg_a, FLO_type arg_b)
1012*c87b03e5Sespie {
1013*c87b03e5Sespie fp_number_type a;
1014*c87b03e5Sespie fp_number_type b;
1015*c87b03e5Sespie fp_number_type *res;
1016*c87b03e5Sespie FLO_union_type au, bu;
1017*c87b03e5Sespie
1018*c87b03e5Sespie au.value = arg_a;
1019*c87b03e5Sespie bu.value = arg_b;
1020*c87b03e5Sespie
1021*c87b03e5Sespie unpack_d (&au, &a);
1022*c87b03e5Sespie unpack_d (&bu, &b);
1023*c87b03e5Sespie
1024*c87b03e5Sespie res = _fpdiv_parts (&a, &b);
1025*c87b03e5Sespie
1026*c87b03e5Sespie return pack_d (res);
1027*c87b03e5Sespie }
1028*c87b03e5Sespie #endif /* L_div_sf || L_div_df */
1029*c87b03e5Sespie
1030*c87b03e5Sespie #if defined(L_fpcmp_parts_sf) || defined(L_fpcmp_parts_df) \
1031*c87b03e5Sespie || defined(L_fpcmp_parts_tf)
1032*c87b03e5Sespie /* according to the demo, fpcmp returns a comparison with 0... thus
1033*c87b03e5Sespie a<b -> -1
1034*c87b03e5Sespie a==b -> 0
1035*c87b03e5Sespie a>b -> +1
1036*c87b03e5Sespie */
1037*c87b03e5Sespie
1038*c87b03e5Sespie int
1039*c87b03e5Sespie __fpcmp_parts (fp_number_type * a, fp_number_type * b)
1040*c87b03e5Sespie {
1041*c87b03e5Sespie #if 0
1042*c87b03e5Sespie /* either nan -> unordered. Must be checked outside of this routine. */
1043*c87b03e5Sespie if (isnan (a) && isnan (b))
1044*c87b03e5Sespie {
1045*c87b03e5Sespie return 1; /* still unordered! */
1046*c87b03e5Sespie }
1047*c87b03e5Sespie #endif
1048*c87b03e5Sespie
1049*c87b03e5Sespie if (isnan (a) || isnan (b))
1050*c87b03e5Sespie {
1051*c87b03e5Sespie return 1; /* how to indicate unordered compare? */
1052*c87b03e5Sespie }
1053*c87b03e5Sespie if (isinf (a) && isinf (b))
1054*c87b03e5Sespie {
1055*c87b03e5Sespie /* +inf > -inf, but +inf != +inf */
1056*c87b03e5Sespie /* b \a| +inf(0)| -inf(1)
1057*c87b03e5Sespie ______\+--------+--------
1058*c87b03e5Sespie +inf(0)| a==b(0)| a<b(-1)
1059*c87b03e5Sespie -------+--------+--------
1060*c87b03e5Sespie -inf(1)| a>b(1) | a==b(0)
1061*c87b03e5Sespie -------+--------+--------
1062*c87b03e5Sespie So since unordered must be nonzero, just line up the columns...
1063*c87b03e5Sespie */
1064*c87b03e5Sespie return b->sign - a->sign;
1065*c87b03e5Sespie }
1066*c87b03e5Sespie /* but not both... */
1067*c87b03e5Sespie if (isinf (a))
1068*c87b03e5Sespie {
1069*c87b03e5Sespie return a->sign ? -1 : 1;
1070*c87b03e5Sespie }
1071*c87b03e5Sespie if (isinf (b))
1072*c87b03e5Sespie {
1073*c87b03e5Sespie return b->sign ? 1 : -1;
1074*c87b03e5Sespie }
1075*c87b03e5Sespie if (iszero (a) && iszero (b))
1076*c87b03e5Sespie {
1077*c87b03e5Sespie return 0;
1078*c87b03e5Sespie }
1079*c87b03e5Sespie if (iszero (a))
1080*c87b03e5Sespie {
1081*c87b03e5Sespie return b->sign ? 1 : -1;
1082*c87b03e5Sespie }
1083*c87b03e5Sespie if (iszero (b))
1084*c87b03e5Sespie {
1085*c87b03e5Sespie return a->sign ? -1 : 1;
1086*c87b03e5Sespie }
1087*c87b03e5Sespie /* now both are "normal". */
1088*c87b03e5Sespie if (a->sign != b->sign)
1089*c87b03e5Sespie {
1090*c87b03e5Sespie /* opposite signs */
1091*c87b03e5Sespie return a->sign ? -1 : 1;
1092*c87b03e5Sespie }
1093*c87b03e5Sespie /* same sign; exponents? */
1094*c87b03e5Sespie if (a->normal_exp > b->normal_exp)
1095*c87b03e5Sespie {
1096*c87b03e5Sespie return a->sign ? -1 : 1;
1097*c87b03e5Sespie }
1098*c87b03e5Sespie if (a->normal_exp < b->normal_exp)
1099*c87b03e5Sespie {
1100*c87b03e5Sespie return a->sign ? 1 : -1;
1101*c87b03e5Sespie }
1102*c87b03e5Sespie /* same exponents; check size. */
1103*c87b03e5Sespie if (a->fraction.ll > b->fraction.ll)
1104*c87b03e5Sespie {
1105*c87b03e5Sespie return a->sign ? -1 : 1;
1106*c87b03e5Sespie }
1107*c87b03e5Sespie if (a->fraction.ll < b->fraction.ll)
1108*c87b03e5Sespie {
1109*c87b03e5Sespie return a->sign ? 1 : -1;
1110*c87b03e5Sespie }
1111*c87b03e5Sespie /* after all that, they're equal. */
1112*c87b03e5Sespie return 0;
1113*c87b03e5Sespie }
1114*c87b03e5Sespie #endif
1115*c87b03e5Sespie
1116*c87b03e5Sespie #if defined(L_compare_sf) || defined(L_compare_df) || defined(L_compoare_tf)
1117*c87b03e5Sespie CMPtype
1118*c87b03e5Sespie compare (FLO_type arg_a, FLO_type arg_b)
1119*c87b03e5Sespie {
1120*c87b03e5Sespie fp_number_type a;
1121*c87b03e5Sespie fp_number_type b;
1122*c87b03e5Sespie FLO_union_type au, bu;
1123*c87b03e5Sespie
1124*c87b03e5Sespie au.value = arg_a;
1125*c87b03e5Sespie bu.value = arg_b;
1126*c87b03e5Sespie
1127*c87b03e5Sespie unpack_d (&au, &a);
1128*c87b03e5Sespie unpack_d (&bu, &b);
1129*c87b03e5Sespie
1130*c87b03e5Sespie return __fpcmp_parts (&a, &b);
1131*c87b03e5Sespie }
1132*c87b03e5Sespie #endif /* L_compare_sf || L_compare_df */
1133*c87b03e5Sespie
1134*c87b03e5Sespie #ifndef US_SOFTWARE_GOFAST
1135*c87b03e5Sespie
1136*c87b03e5Sespie /* These should be optimized for their specific tasks someday. */
1137*c87b03e5Sespie
1138*c87b03e5Sespie #if defined(L_eq_sf) || defined(L_eq_df) || defined(L_eq_tf)
1139*c87b03e5Sespie CMPtype
1140*c87b03e5Sespie _eq_f2 (FLO_type arg_a, FLO_type arg_b)
1141*c87b03e5Sespie {
1142*c87b03e5Sespie fp_number_type a;
1143*c87b03e5Sespie fp_number_type b;
1144*c87b03e5Sespie FLO_union_type au, bu;
1145*c87b03e5Sespie
1146*c87b03e5Sespie au.value = arg_a;
1147*c87b03e5Sespie bu.value = arg_b;
1148*c87b03e5Sespie
1149*c87b03e5Sespie unpack_d (&au, &a);
1150*c87b03e5Sespie unpack_d (&bu, &b);
1151*c87b03e5Sespie
1152*c87b03e5Sespie if (isnan (&a) || isnan (&b))
1153*c87b03e5Sespie return 1; /* false, truth == 0 */
1154*c87b03e5Sespie
1155*c87b03e5Sespie return __fpcmp_parts (&a, &b) ;
1156*c87b03e5Sespie }
1157*c87b03e5Sespie #endif /* L_eq_sf || L_eq_df */
1158*c87b03e5Sespie
1159*c87b03e5Sespie #if defined(L_ne_sf) || defined(L_ne_df) || defined(L_ne_tf)
1160*c87b03e5Sespie CMPtype
1161*c87b03e5Sespie _ne_f2 (FLO_type arg_a, FLO_type arg_b)
1162*c87b03e5Sespie {
1163*c87b03e5Sespie fp_number_type a;
1164*c87b03e5Sespie fp_number_type b;
1165*c87b03e5Sespie FLO_union_type au, bu;
1166*c87b03e5Sespie
1167*c87b03e5Sespie au.value = arg_a;
1168*c87b03e5Sespie bu.value = arg_b;
1169*c87b03e5Sespie
1170*c87b03e5Sespie unpack_d (&au, &a);
1171*c87b03e5Sespie unpack_d (&bu, &b);
1172*c87b03e5Sespie
1173*c87b03e5Sespie if (isnan (&a) || isnan (&b))
1174*c87b03e5Sespie return 1; /* true, truth != 0 */
1175*c87b03e5Sespie
1176*c87b03e5Sespie return __fpcmp_parts (&a, &b) ;
1177*c87b03e5Sespie }
1178*c87b03e5Sespie #endif /* L_ne_sf || L_ne_df */
1179*c87b03e5Sespie
1180*c87b03e5Sespie #if defined(L_gt_sf) || defined(L_gt_df) || defined(L_gt_tf)
1181*c87b03e5Sespie CMPtype
1182*c87b03e5Sespie _gt_f2 (FLO_type arg_a, FLO_type arg_b)
1183*c87b03e5Sespie {
1184*c87b03e5Sespie fp_number_type a;
1185*c87b03e5Sespie fp_number_type b;
1186*c87b03e5Sespie FLO_union_type au, bu;
1187*c87b03e5Sespie
1188*c87b03e5Sespie au.value = arg_a;
1189*c87b03e5Sespie bu.value = arg_b;
1190*c87b03e5Sespie
1191*c87b03e5Sespie unpack_d (&au, &a);
1192*c87b03e5Sespie unpack_d (&bu, &b);
1193*c87b03e5Sespie
1194*c87b03e5Sespie if (isnan (&a) || isnan (&b))
1195*c87b03e5Sespie return -1; /* false, truth > 0 */
1196*c87b03e5Sespie
1197*c87b03e5Sespie return __fpcmp_parts (&a, &b);
1198*c87b03e5Sespie }
1199*c87b03e5Sespie #endif /* L_gt_sf || L_gt_df */
1200*c87b03e5Sespie
1201*c87b03e5Sespie #if defined(L_ge_sf) || defined(L_ge_df) || defined(L_ge_tf)
1202*c87b03e5Sespie CMPtype
1203*c87b03e5Sespie _ge_f2 (FLO_type arg_a, FLO_type arg_b)
1204*c87b03e5Sespie {
1205*c87b03e5Sespie fp_number_type a;
1206*c87b03e5Sespie fp_number_type b;
1207*c87b03e5Sespie FLO_union_type au, bu;
1208*c87b03e5Sespie
1209*c87b03e5Sespie au.value = arg_a;
1210*c87b03e5Sespie bu.value = arg_b;
1211*c87b03e5Sespie
1212*c87b03e5Sespie unpack_d (&au, &a);
1213*c87b03e5Sespie unpack_d (&bu, &b);
1214*c87b03e5Sespie
1215*c87b03e5Sespie if (isnan (&a) || isnan (&b))
1216*c87b03e5Sespie return -1; /* false, truth >= 0 */
1217*c87b03e5Sespie return __fpcmp_parts (&a, &b) ;
1218*c87b03e5Sespie }
1219*c87b03e5Sespie #endif /* L_ge_sf || L_ge_df */
1220*c87b03e5Sespie
1221*c87b03e5Sespie #if defined(L_lt_sf) || defined(L_lt_df) || defined(L_lt_tf)
1222*c87b03e5Sespie CMPtype
1223*c87b03e5Sespie _lt_f2 (FLO_type arg_a, FLO_type arg_b)
1224*c87b03e5Sespie {
1225*c87b03e5Sespie fp_number_type a;
1226*c87b03e5Sespie fp_number_type b;
1227*c87b03e5Sespie FLO_union_type au, bu;
1228*c87b03e5Sespie
1229*c87b03e5Sespie au.value = arg_a;
1230*c87b03e5Sespie bu.value = arg_b;
1231*c87b03e5Sespie
1232*c87b03e5Sespie unpack_d (&au, &a);
1233*c87b03e5Sespie unpack_d (&bu, &b);
1234*c87b03e5Sespie
1235*c87b03e5Sespie if (isnan (&a) || isnan (&b))
1236*c87b03e5Sespie return 1; /* false, truth < 0 */
1237*c87b03e5Sespie
1238*c87b03e5Sespie return __fpcmp_parts (&a, &b);
1239*c87b03e5Sespie }
1240*c87b03e5Sespie #endif /* L_lt_sf || L_lt_df */
1241*c87b03e5Sespie
1242*c87b03e5Sespie #if defined(L_le_sf) || defined(L_le_df) || defined(L_le_tf)
1243*c87b03e5Sespie CMPtype
1244*c87b03e5Sespie _le_f2 (FLO_type arg_a, FLO_type arg_b)
1245*c87b03e5Sespie {
1246*c87b03e5Sespie fp_number_type a;
1247*c87b03e5Sespie fp_number_type b;
1248*c87b03e5Sespie FLO_union_type au, bu;
1249*c87b03e5Sespie
1250*c87b03e5Sespie au.value = arg_a;
1251*c87b03e5Sespie bu.value = arg_b;
1252*c87b03e5Sespie
1253*c87b03e5Sespie unpack_d (&au, &a);
1254*c87b03e5Sespie unpack_d (&bu, &b);
1255*c87b03e5Sespie
1256*c87b03e5Sespie if (isnan (&a) || isnan (&b))
1257*c87b03e5Sespie return 1; /* false, truth <= 0 */
1258*c87b03e5Sespie
1259*c87b03e5Sespie return __fpcmp_parts (&a, &b) ;
1260*c87b03e5Sespie }
1261*c87b03e5Sespie #endif /* L_le_sf || L_le_df */
1262*c87b03e5Sespie
1263*c87b03e5Sespie #endif /* ! US_SOFTWARE_GOFAST */
1264*c87b03e5Sespie
1265*c87b03e5Sespie #if defined(L_unord_sf) || defined(L_unord_df) || defined(L_unord_tf)
1266*c87b03e5Sespie CMPtype
1267*c87b03e5Sespie _unord_f2 (FLO_type arg_a, FLO_type arg_b)
1268*c87b03e5Sespie {
1269*c87b03e5Sespie fp_number_type a;
1270*c87b03e5Sespie fp_number_type b;
1271*c87b03e5Sespie FLO_union_type au, bu;
1272*c87b03e5Sespie
1273*c87b03e5Sespie au.value = arg_a;
1274*c87b03e5Sespie bu.value = arg_b;
1275*c87b03e5Sespie
1276*c87b03e5Sespie unpack_d (&au, &a);
1277*c87b03e5Sespie unpack_d (&bu, &b);
1278*c87b03e5Sespie
1279*c87b03e5Sespie return (isnan (&a) || isnan (&b));
1280*c87b03e5Sespie }
1281*c87b03e5Sespie #endif /* L_unord_sf || L_unord_df */
1282*c87b03e5Sespie
1283*c87b03e5Sespie #if defined(L_si_to_sf) || defined(L_si_to_df) || defined(L_si_to_tf)
1284*c87b03e5Sespie FLO_type
1285*c87b03e5Sespie si_to_float (SItype arg_a)
1286*c87b03e5Sespie {
1287*c87b03e5Sespie fp_number_type in;
1288*c87b03e5Sespie
1289*c87b03e5Sespie in.class = CLASS_NUMBER;
1290*c87b03e5Sespie in.sign = arg_a < 0;
1291*c87b03e5Sespie if (!arg_a)
1292*c87b03e5Sespie {
1293*c87b03e5Sespie in.class = CLASS_ZERO;
1294*c87b03e5Sespie }
1295*c87b03e5Sespie else
1296*c87b03e5Sespie {
1297*c87b03e5Sespie in.normal_exp = FRACBITS + NGARDS;
1298*c87b03e5Sespie if (in.sign)
1299*c87b03e5Sespie {
1300*c87b03e5Sespie /* Special case for minint, since there is no +ve integer
1301*c87b03e5Sespie representation for it */
1302*c87b03e5Sespie if (arg_a == (- MAX_SI_INT - 1))
1303*c87b03e5Sespie {
1304*c87b03e5Sespie return (FLO_type)(- MAX_SI_INT - 1);
1305*c87b03e5Sespie }
1306*c87b03e5Sespie in.fraction.ll = (-arg_a);
1307*c87b03e5Sespie }
1308*c87b03e5Sespie else
1309*c87b03e5Sespie in.fraction.ll = arg_a;
1310*c87b03e5Sespie
1311*c87b03e5Sespie while (in.fraction.ll < ((fractype)1 << (FRACBITS + NGARDS)))
1312*c87b03e5Sespie {
1313*c87b03e5Sespie in.fraction.ll <<= 1;
1314*c87b03e5Sespie in.normal_exp -= 1;
1315*c87b03e5Sespie }
1316*c87b03e5Sespie }
1317*c87b03e5Sespie return pack_d (&in);
1318*c87b03e5Sespie }
1319*c87b03e5Sespie #endif /* L_si_to_sf || L_si_to_df */
1320*c87b03e5Sespie
1321*c87b03e5Sespie #if defined(L_usi_to_sf) || defined(L_usi_to_df) || defined(L_usi_to_tf)
1322*c87b03e5Sespie FLO_type
1323*c87b03e5Sespie usi_to_float (USItype arg_a)
1324*c87b03e5Sespie {
1325*c87b03e5Sespie fp_number_type in;
1326*c87b03e5Sespie
1327*c87b03e5Sespie in.sign = 0;
1328*c87b03e5Sespie if (!arg_a)
1329*c87b03e5Sespie {
1330*c87b03e5Sespie in.class = CLASS_ZERO;
1331*c87b03e5Sespie }
1332*c87b03e5Sespie else
1333*c87b03e5Sespie {
1334*c87b03e5Sespie in.class = CLASS_NUMBER;
1335*c87b03e5Sespie in.normal_exp = FRACBITS + NGARDS;
1336*c87b03e5Sespie in.fraction.ll = arg_a;
1337*c87b03e5Sespie
1338*c87b03e5Sespie while (in.fraction.ll > ((fractype)1 << (FRACBITS + NGARDS)))
1339*c87b03e5Sespie {
1340*c87b03e5Sespie in.fraction.ll >>= 1;
1341*c87b03e5Sespie in.normal_exp += 1;
1342*c87b03e5Sespie }
1343*c87b03e5Sespie while (in.fraction.ll < ((fractype)1 << (FRACBITS + NGARDS)))
1344*c87b03e5Sespie {
1345*c87b03e5Sespie in.fraction.ll <<= 1;
1346*c87b03e5Sespie in.normal_exp -= 1;
1347*c87b03e5Sespie }
1348*c87b03e5Sespie }
1349*c87b03e5Sespie return pack_d (&in);
1350*c87b03e5Sespie }
1351*c87b03e5Sespie #endif
1352*c87b03e5Sespie
1353*c87b03e5Sespie #if defined(L_sf_to_si) || defined(L_df_to_si) || defined(L_tf_to_si)
1354*c87b03e5Sespie SItype
1355*c87b03e5Sespie float_to_si (FLO_type arg_a)
1356*c87b03e5Sespie {
1357*c87b03e5Sespie fp_number_type a;
1358*c87b03e5Sespie SItype tmp;
1359*c87b03e5Sespie FLO_union_type au;
1360*c87b03e5Sespie
1361*c87b03e5Sespie au.value = arg_a;
1362*c87b03e5Sespie unpack_d (&au, &a);
1363*c87b03e5Sespie
1364*c87b03e5Sespie if (iszero (&a))
1365*c87b03e5Sespie return 0;
1366*c87b03e5Sespie if (isnan (&a))
1367*c87b03e5Sespie return 0;
1368*c87b03e5Sespie /* get reasonable MAX_SI_INT... */
1369*c87b03e5Sespie if (isinf (&a))
1370*c87b03e5Sespie return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT;
1371*c87b03e5Sespie /* it is a number, but a small one */
1372*c87b03e5Sespie if (a.normal_exp < 0)
1373*c87b03e5Sespie return 0;
1374*c87b03e5Sespie if (a.normal_exp > BITS_PER_SI - 2)
1375*c87b03e5Sespie return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT;
1376*c87b03e5Sespie tmp = a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1377*c87b03e5Sespie return a.sign ? (-tmp) : (tmp);
1378*c87b03e5Sespie }
1379*c87b03e5Sespie #endif /* L_sf_to_si || L_df_to_si */
1380*c87b03e5Sespie
1381*c87b03e5Sespie #if defined(L_sf_to_usi) || defined(L_df_to_usi) || defined(L_tf_to_usi)
1382*c87b03e5Sespie #if defined US_SOFTWARE_GOFAST || defined(L_tf_to_usi)
1383*c87b03e5Sespie /* While libgcc2.c defines its own __fixunssfsi and __fixunsdfsi routines,
1384*c87b03e5Sespie we also define them for GOFAST because the ones in libgcc2.c have the
1385*c87b03e5Sespie wrong names and I'd rather define these here and keep GOFAST CYG-LOC's
1386*c87b03e5Sespie out of libgcc2.c. We can't define these here if not GOFAST because then
1387*c87b03e5Sespie there'd be duplicate copies. */
1388*c87b03e5Sespie
1389*c87b03e5Sespie USItype
1390*c87b03e5Sespie float_to_usi (FLO_type arg_a)
1391*c87b03e5Sespie {
1392*c87b03e5Sespie fp_number_type a;
1393*c87b03e5Sespie FLO_union_type au;
1394*c87b03e5Sespie
1395*c87b03e5Sespie au.value = arg_a;
1396*c87b03e5Sespie unpack_d (&au, &a);
1397*c87b03e5Sespie
1398*c87b03e5Sespie if (iszero (&a))
1399*c87b03e5Sespie return 0;
1400*c87b03e5Sespie if (isnan (&a))
1401*c87b03e5Sespie return 0;
1402*c87b03e5Sespie /* it is a negative number */
1403*c87b03e5Sespie if (a.sign)
1404*c87b03e5Sespie return 0;
1405*c87b03e5Sespie /* get reasonable MAX_USI_INT... */
1406*c87b03e5Sespie if (isinf (&a))
1407*c87b03e5Sespie return MAX_USI_INT;
1408*c87b03e5Sespie /* it is a number, but a small one */
1409*c87b03e5Sespie if (a.normal_exp < 0)
1410*c87b03e5Sespie return 0;
1411*c87b03e5Sespie if (a.normal_exp > BITS_PER_SI - 1)
1412*c87b03e5Sespie return MAX_USI_INT;
1413*c87b03e5Sespie else if (a.normal_exp > (FRACBITS + NGARDS))
1414*c87b03e5Sespie return a.fraction.ll << (a.normal_exp - (FRACBITS + NGARDS));
1415*c87b03e5Sespie else
1416*c87b03e5Sespie return a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1417*c87b03e5Sespie }
1418*c87b03e5Sespie #endif /* US_SOFTWARE_GOFAST */
1419*c87b03e5Sespie #endif /* L_sf_to_usi || L_df_to_usi */
1420*c87b03e5Sespie
1421*c87b03e5Sespie #if defined(L_negate_sf) || defined(L_negate_df) || defined(L_negate_tf)
1422*c87b03e5Sespie FLO_type
1423*c87b03e5Sespie negate (FLO_type arg_a)
1424*c87b03e5Sespie {
1425*c87b03e5Sespie fp_number_type a;
1426*c87b03e5Sespie FLO_union_type au;
1427*c87b03e5Sespie
1428*c87b03e5Sespie au.value = arg_a;
1429*c87b03e5Sespie unpack_d (&au, &a);
1430*c87b03e5Sespie
1431*c87b03e5Sespie flip_sign (&a);
1432*c87b03e5Sespie return pack_d (&a);
1433*c87b03e5Sespie }
1434*c87b03e5Sespie #endif /* L_negate_sf || L_negate_df */
1435*c87b03e5Sespie
1436*c87b03e5Sespie #ifdef FLOAT
1437*c87b03e5Sespie
1438*c87b03e5Sespie #if defined(L_make_sf)
1439*c87b03e5Sespie SFtype
1440*c87b03e5Sespie __make_fp(fp_class_type class,
1441*c87b03e5Sespie unsigned int sign,
1442*c87b03e5Sespie int exp,
1443*c87b03e5Sespie USItype frac)
1444*c87b03e5Sespie {
1445*c87b03e5Sespie fp_number_type in;
1446*c87b03e5Sespie
1447*c87b03e5Sespie in.class = class;
1448*c87b03e5Sespie in.sign = sign;
1449*c87b03e5Sespie in.normal_exp = exp;
1450*c87b03e5Sespie in.fraction.ll = frac;
1451*c87b03e5Sespie return pack_d (&in);
1452*c87b03e5Sespie }
1453*c87b03e5Sespie #endif /* L_make_sf */
1454*c87b03e5Sespie
1455*c87b03e5Sespie #ifndef FLOAT_ONLY
1456*c87b03e5Sespie
1457*c87b03e5Sespie /* This enables one to build an fp library that supports float but not double.
1458*c87b03e5Sespie Otherwise, we would get an undefined reference to __make_dp.
1459*c87b03e5Sespie This is needed for some 8-bit ports that can't handle well values that
1460*c87b03e5Sespie are 8-bytes in size, so we just don't support double for them at all. */
1461*c87b03e5Sespie
1462*c87b03e5Sespie #if defined(L_sf_to_df)
1463*c87b03e5Sespie DFtype
1464*c87b03e5Sespie sf_to_df (SFtype arg_a)
1465*c87b03e5Sespie {
1466*c87b03e5Sespie fp_number_type in;
1467*c87b03e5Sespie FLO_union_type au;
1468*c87b03e5Sespie
1469*c87b03e5Sespie au.value = arg_a;
1470*c87b03e5Sespie unpack_d (&au, &in);
1471*c87b03e5Sespie
1472*c87b03e5Sespie return __make_dp (in.class, in.sign, in.normal_exp,
1473*c87b03e5Sespie ((UDItype) in.fraction.ll) << F_D_BITOFF);
1474*c87b03e5Sespie }
1475*c87b03e5Sespie #endif /* L_sf_to_df */
1476*c87b03e5Sespie
1477*c87b03e5Sespie #if defined(L_sf_to_tf) && defined(TMODES)
1478*c87b03e5Sespie TFtype
1479*c87b03e5Sespie sf_to_tf (SFtype arg_a)
1480*c87b03e5Sespie {
1481*c87b03e5Sespie fp_number_type in;
1482*c87b03e5Sespie FLO_union_type au;
1483*c87b03e5Sespie
1484*c87b03e5Sespie au.value = arg_a;
1485*c87b03e5Sespie unpack_d (&au, &in);
1486*c87b03e5Sespie
1487*c87b03e5Sespie return __make_tp (in.class, in.sign, in.normal_exp,
1488*c87b03e5Sespie ((UTItype) in.fraction.ll) << F_T_BITOFF);
1489*c87b03e5Sespie }
1490*c87b03e5Sespie #endif /* L_sf_to_df */
1491*c87b03e5Sespie
1492*c87b03e5Sespie #endif /* ! FLOAT_ONLY */
1493*c87b03e5Sespie #endif /* FLOAT */
1494*c87b03e5Sespie
1495*c87b03e5Sespie #ifndef FLOAT
1496*c87b03e5Sespie
1497*c87b03e5Sespie extern SFtype __make_fp (fp_class_type, unsigned int, int, USItype);
1498*c87b03e5Sespie
1499*c87b03e5Sespie #if defined(L_make_df)
1500*c87b03e5Sespie DFtype
1501*c87b03e5Sespie __make_dp (fp_class_type class, unsigned int sign, int exp, UDItype frac)
1502*c87b03e5Sespie {
1503*c87b03e5Sespie fp_number_type in;
1504*c87b03e5Sespie
1505*c87b03e5Sespie in.class = class;
1506*c87b03e5Sespie in.sign = sign;
1507*c87b03e5Sespie in.normal_exp = exp;
1508*c87b03e5Sespie in.fraction.ll = frac;
1509*c87b03e5Sespie return pack_d (&in);
1510*c87b03e5Sespie }
1511*c87b03e5Sespie #endif /* L_make_df */
1512*c87b03e5Sespie
1513*c87b03e5Sespie #if defined(L_df_to_sf)
1514*c87b03e5Sespie SFtype
1515*c87b03e5Sespie df_to_sf (DFtype arg_a)
1516*c87b03e5Sespie {
1517*c87b03e5Sespie fp_number_type in;
1518*c87b03e5Sespie USItype sffrac;
1519*c87b03e5Sespie FLO_union_type au;
1520*c87b03e5Sespie
1521*c87b03e5Sespie au.value = arg_a;
1522*c87b03e5Sespie unpack_d (&au, &in);
1523*c87b03e5Sespie
1524*c87b03e5Sespie sffrac = in.fraction.ll >> F_D_BITOFF;
1525*c87b03e5Sespie
1526*c87b03e5Sespie /* We set the lowest guard bit in SFFRAC if we discarded any non
1527*c87b03e5Sespie zero bits. */
1528*c87b03e5Sespie if ((in.fraction.ll & (((USItype) 1 << F_D_BITOFF) - 1)) != 0)
1529*c87b03e5Sespie sffrac |= 1;
1530*c87b03e5Sespie
1531*c87b03e5Sespie return __make_fp (in.class, in.sign, in.normal_exp, sffrac);
1532*c87b03e5Sespie }
1533*c87b03e5Sespie #endif /* L_df_to_sf */
1534*c87b03e5Sespie
1535*c87b03e5Sespie #if defined(L_df_to_tf) && defined(TMODES) \
1536*c87b03e5Sespie && !defined(FLOAT) && !defined(TFLOAT)
1537*c87b03e5Sespie TFtype
1538*c87b03e5Sespie df_to_tf (DFtype arg_a)
1539*c87b03e5Sespie {
1540*c87b03e5Sespie fp_number_type in;
1541*c87b03e5Sespie FLO_union_type au;
1542*c87b03e5Sespie
1543*c87b03e5Sespie au.value = arg_a;
1544*c87b03e5Sespie unpack_d (&au, &in);
1545*c87b03e5Sespie
1546*c87b03e5Sespie return __make_tp (in.class, in.sign, in.normal_exp,
1547*c87b03e5Sespie ((UTItype) in.fraction.ll) << D_T_BITOFF);
1548*c87b03e5Sespie }
1549*c87b03e5Sespie #endif /* L_sf_to_df */
1550*c87b03e5Sespie
1551*c87b03e5Sespie #ifdef TFLOAT
1552*c87b03e5Sespie #if defined(L_make_tf)
1553*c87b03e5Sespie TFtype
1554*c87b03e5Sespie __make_tp(fp_class_type class,
1555*c87b03e5Sespie unsigned int sign,
1556*c87b03e5Sespie int exp,
1557*c87b03e5Sespie UTItype frac)
1558*c87b03e5Sespie {
1559*c87b03e5Sespie fp_number_type in;
1560*c87b03e5Sespie
1561*c87b03e5Sespie in.class = class;
1562*c87b03e5Sespie in.sign = sign;
1563*c87b03e5Sespie in.normal_exp = exp;
1564*c87b03e5Sespie in.fraction.ll = frac;
1565*c87b03e5Sespie return pack_d (&in);
1566*c87b03e5Sespie }
1567*c87b03e5Sespie #endif /* L_make_tf */
1568*c87b03e5Sespie
1569*c87b03e5Sespie #if defined(L_tf_to_df)
1570*c87b03e5Sespie DFtype
1571*c87b03e5Sespie tf_to_df (TFtype arg_a)
1572*c87b03e5Sespie {
1573*c87b03e5Sespie fp_number_type in;
1574*c87b03e5Sespie UDItype sffrac;
1575*c87b03e5Sespie FLO_union_type au;
1576*c87b03e5Sespie
1577*c87b03e5Sespie au.value = arg_a;
1578*c87b03e5Sespie unpack_d (&au, &in);
1579*c87b03e5Sespie
1580*c87b03e5Sespie sffrac = in.fraction.ll >> D_T_BITOFF;
1581*c87b03e5Sespie
1582*c87b03e5Sespie /* We set the lowest guard bit in SFFRAC if we discarded any non
1583*c87b03e5Sespie zero bits. */
1584*c87b03e5Sespie if ((in.fraction.ll & (((UTItype) 1 << D_T_BITOFF) - 1)) != 0)
1585*c87b03e5Sespie sffrac |= 1;
1586*c87b03e5Sespie
1587*c87b03e5Sespie return __make_dp (in.class, in.sign, in.normal_exp, sffrac);
1588*c87b03e5Sespie }
1589*c87b03e5Sespie #endif /* L_tf_to_df */
1590*c87b03e5Sespie
1591*c87b03e5Sespie #if defined(L_tf_to_sf)
1592*c87b03e5Sespie SFtype
1593*c87b03e5Sespie tf_to_sf (TFtype arg_a)
1594*c87b03e5Sespie {
1595*c87b03e5Sespie fp_number_type in;
1596*c87b03e5Sespie USItype sffrac;
1597*c87b03e5Sespie FLO_union_type au;
1598*c87b03e5Sespie
1599*c87b03e5Sespie au.value = arg_a;
1600*c87b03e5Sespie unpack_d (&au, &in);
1601*c87b03e5Sespie
1602*c87b03e5Sespie sffrac = in.fraction.ll >> F_T_BITOFF;
1603*c87b03e5Sespie
1604*c87b03e5Sespie /* We set the lowest guard bit in SFFRAC if we discarded any non
1605*c87b03e5Sespie zero bits. */
1606*c87b03e5Sespie if ((in.fraction.ll & (((UTItype) 1 << F_T_BITOFF) - 1)) != 0)
1607*c87b03e5Sespie sffrac |= 1;
1608*c87b03e5Sespie
1609*c87b03e5Sespie return __make_fp (in.class, in.sign, in.normal_exp, sffrac);
1610*c87b03e5Sespie }
1611*c87b03e5Sespie #endif /* L_tf_to_sf */
1612*c87b03e5Sespie #endif /* TFLOAT */
1613*c87b03e5Sespie
1614*c87b03e5Sespie #endif /* ! FLOAT */
1615*c87b03e5Sespie #endif /* !EXTENDED_FLOAT_STUBS */
1616