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