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