111
|
1 /* This is a software floating point library which can be used
|
|
2 for targets without hardware floating point.
|
145
|
3 Copyright (C) 1994-2020 Free Software Foundation, Inc.
|
111
|
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 */
|
145
|
319 #if defined TFLOAT && defined HALFFRACBITS
|
111
|
320 {
|
|
321 halffractype high, low, unity;
|
|
322 int lowsign, lowexp;
|
|
323
|
|
324 unity = (halffractype) 1 << HALFFRACBITS;
|
|
325
|
|
326 /* Set HIGH to the high double's significand, masking out the implicit 1.
|
|
327 Set LOW to the low double's full significand. */
|
|
328 high = (fraction >> (FRACBITS - HALFFRACBITS)) & (unity - 1);
|
|
329 low = fraction & (unity * 2 - 1);
|
|
330
|
|
331 /* Get the initial sign and exponent of the low double. */
|
|
332 lowexp = exp - HALFFRACBITS - 1;
|
|
333 lowsign = sign;
|
|
334
|
|
335 /* HIGH should be rounded like a normal double, making |LOW| <=
|
|
336 0.5 ULP of HIGH. Assume round-to-nearest. */
|
|
337 if (exp < EXPMAX)
|
|
338 if (low > unity || (low == unity && (high & 1) == 1))
|
|
339 {
|
|
340 /* Round HIGH up and adjust LOW to match. */
|
|
341 high++;
|
|
342 if (high == unity)
|
|
343 {
|
|
344 /* May make it infinite, but that's OK. */
|
|
345 high = 0;
|
|
346 exp++;
|
|
347 }
|
|
348 low = unity * 2 - low;
|
|
349 lowsign ^= 1;
|
|
350 }
|
|
351
|
|
352 high |= (halffractype) exp << HALFFRACBITS;
|
|
353 high |= (halffractype) sign << (HALFFRACBITS + EXPBITS);
|
|
354
|
|
355 if (exp == EXPMAX || exp == 0 || low == 0)
|
|
356 low = 0;
|
|
357 else
|
|
358 {
|
|
359 while (lowexp > 0 && low < unity)
|
|
360 {
|
|
361 low <<= 1;
|
|
362 lowexp--;
|
|
363 }
|
|
364
|
|
365 if (lowexp <= 0)
|
|
366 {
|
|
367 halffractype roundmsb, round;
|
|
368 int shift;
|
|
369
|
|
370 shift = 1 - lowexp;
|
|
371 roundmsb = (1 << (shift - 1));
|
|
372 round = low & ((roundmsb << 1) - 1);
|
|
373
|
|
374 low >>= shift;
|
|
375 lowexp = 0;
|
|
376
|
|
377 if (round > roundmsb || (round == roundmsb && (low & 1) == 1))
|
|
378 {
|
|
379 low++;
|
|
380 if (low == unity)
|
|
381 /* LOW rounds up to the smallest normal number. */
|
|
382 lowexp++;
|
|
383 }
|
|
384 }
|
|
385
|
|
386 low &= unity - 1;
|
|
387 low |= (halffractype) lowexp << HALFFRACBITS;
|
|
388 low |= (halffractype) lowsign << (HALFFRACBITS + EXPBITS);
|
|
389 }
|
|
390 dst.value_raw = ((fractype) high << HALFSHIFT) | low;
|
|
391 }
|
145
|
392 #else
|
111
|
393 dst.value_raw = fraction & ((((fractype)1) << FRACBITS) - (fractype)1);
|
|
394 dst.value_raw |= ((fractype) (exp & ((1 << EXPBITS) - 1))) << FRACBITS;
|
|
395 dst.value_raw |= ((fractype) (sign & 1)) << (FRACBITS | EXPBITS);
|
|
396 #endif
|
|
397
|
|
398 #if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT)
|
|
399 #ifdef TFLOAT
|
|
400 {
|
|
401 qrtrfractype tmp1 = dst.words[0];
|
|
402 qrtrfractype tmp2 = dst.words[1];
|
|
403 dst.words[0] = dst.words[3];
|
|
404 dst.words[1] = dst.words[2];
|
|
405 dst.words[2] = tmp2;
|
|
406 dst.words[3] = tmp1;
|
|
407 }
|
|
408 #else
|
|
409 {
|
|
410 halffractype tmp = dst.words[0];
|
|
411 dst.words[0] = dst.words[1];
|
|
412 dst.words[1] = tmp;
|
|
413 }
|
|
414 #endif
|
|
415 #endif
|
|
416
|
|
417 return dst.value;
|
|
418 }
|
|
419 #endif
|
|
420
|
|
421 #if defined(L_unpack_df) || defined(L_unpack_sf) || defined(L_unpack_tf)
|
|
422 void
|
|
423 unpack_d (FLO_union_type * src, fp_number_type * dst)
|
|
424 {
|
|
425 /* We previously used bitfields to store the number, but this doesn't
|
|
426 handle little/big endian systems conveniently, so use shifts and
|
|
427 masks */
|
|
428 fractype fraction;
|
|
429 int exp;
|
|
430 int sign;
|
|
431
|
|
432 #if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT)
|
|
433 FLO_union_type swapped;
|
|
434
|
|
435 #ifdef TFLOAT
|
|
436 swapped.words[0] = src->words[3];
|
|
437 swapped.words[1] = src->words[2];
|
|
438 swapped.words[2] = src->words[1];
|
|
439 swapped.words[3] = src->words[0];
|
|
440 #else
|
|
441 swapped.words[0] = src->words[1];
|
|
442 swapped.words[1] = src->words[0];
|
|
443 #endif
|
|
444 src = &swapped;
|
|
445 #endif
|
|
446
|
145
|
447 #if defined TFLOAT && defined HALFFRACBITS
|
111
|
448 {
|
|
449 halffractype high, low;
|
|
450
|
|
451 high = src->value_raw >> HALFSHIFT;
|
|
452 low = src->value_raw & (((fractype)1 << HALFSHIFT) - 1);
|
|
453
|
|
454 fraction = high & ((((fractype)1) << HALFFRACBITS) - 1);
|
|
455 fraction <<= FRACBITS - HALFFRACBITS;
|
|
456 exp = ((int)(high >> HALFFRACBITS)) & ((1 << EXPBITS) - 1);
|
|
457 sign = ((int)(high >> (((HALFFRACBITS + EXPBITS))))) & 1;
|
|
458
|
|
459 if (exp != EXPMAX && exp != 0 && low != 0)
|
|
460 {
|
|
461 int lowexp = ((int)(low >> HALFFRACBITS)) & ((1 << EXPBITS) - 1);
|
|
462 int lowsign = ((int)(low >> (((HALFFRACBITS + EXPBITS))))) & 1;
|
|
463 int shift;
|
|
464 fractype xlow;
|
|
465
|
|
466 xlow = low & ((((fractype)1) << HALFFRACBITS) - 1);
|
|
467 if (lowexp)
|
|
468 xlow |= (((halffractype)1) << HALFFRACBITS);
|
|
469 else
|
|
470 lowexp = 1;
|
|
471 shift = (FRACBITS - HALFFRACBITS) - (exp - lowexp);
|
|
472 if (shift > 0)
|
|
473 xlow <<= shift;
|
|
474 else if (shift < 0)
|
|
475 xlow >>= -shift;
|
|
476 if (sign == lowsign)
|
|
477 fraction += xlow;
|
|
478 else if (fraction >= xlow)
|
|
479 fraction -= xlow;
|
|
480 else
|
|
481 {
|
|
482 /* The high part is a power of two but the full number is lower.
|
|
483 This code will leave the implicit 1 in FRACTION, but we'd
|
|
484 have added that below anyway. */
|
|
485 fraction = (((fractype) 1 << FRACBITS) - xlow) << 1;
|
|
486 exp--;
|
|
487 }
|
|
488 }
|
|
489 }
|
145
|
490 #else
|
111
|
491 fraction = src->value_raw & ((((fractype)1) << FRACBITS) - 1);
|
|
492 exp = ((int)(src->value_raw >> FRACBITS)) & ((1 << EXPBITS) - 1);
|
|
493 sign = ((int)(src->value_raw >> (FRACBITS + EXPBITS))) & 1;
|
|
494 #endif
|
|
495
|
|
496 dst->sign = sign;
|
|
497 if (exp == 0)
|
|
498 {
|
|
499 /* Hmm. Looks like 0 */
|
|
500 if (fraction == 0
|
|
501 #ifdef NO_DENORMALS
|
|
502 || 1
|
|
503 #endif
|
|
504 )
|
|
505 {
|
|
506 /* tastes like zero */
|
|
507 dst->class = CLASS_ZERO;
|
|
508 }
|
|
509 else
|
|
510 {
|
|
511 /* Zero exponent with nonzero fraction - it's denormalized,
|
|
512 so there isn't a leading implicit one - we'll shift it so
|
|
513 it gets one. */
|
|
514 dst->normal_exp = exp - EXPBIAS + 1;
|
|
515 fraction <<= NGARDS;
|
|
516
|
|
517 dst->class = CLASS_NUMBER;
|
|
518 #if 1
|
|
519 while (fraction < IMPLICIT_1)
|
|
520 {
|
|
521 fraction <<= 1;
|
|
522 dst->normal_exp--;
|
|
523 }
|
|
524 #endif
|
|
525 dst->fraction.ll = fraction;
|
|
526 }
|
|
527 }
|
|
528 else if (__builtin_expect (exp == EXPMAX, 0))
|
|
529 {
|
|
530 /* Huge exponent*/
|
|
531 if (fraction == 0)
|
|
532 {
|
|
533 /* Attached to a zero fraction - means infinity */
|
|
534 dst->class = CLASS_INFINITY;
|
|
535 }
|
|
536 else
|
|
537 {
|
|
538 /* Nonzero fraction, means nan */
|
|
539 #ifdef QUIET_NAN_NEGATED
|
|
540 if ((fraction & QUIET_NAN) == 0)
|
|
541 #else
|
|
542 if (fraction & QUIET_NAN)
|
|
543 #endif
|
|
544 {
|
|
545 dst->class = CLASS_QNAN;
|
|
546 }
|
|
547 else
|
|
548 {
|
|
549 dst->class = CLASS_SNAN;
|
|
550 }
|
|
551 /* Now that we know which kind of NaN we got, discard the
|
|
552 quiet/signaling bit, but do preserve the NaN payload. */
|
|
553 fraction &= ~QUIET_NAN;
|
|
554 dst->fraction.ll = fraction << NGARDS;
|
|
555 }
|
|
556 }
|
|
557 else
|
|
558 {
|
|
559 /* Nothing strange about this number */
|
|
560 dst->normal_exp = exp - EXPBIAS;
|
|
561 dst->class = CLASS_NUMBER;
|
|
562 dst->fraction.ll = (fraction << NGARDS) | IMPLICIT_1;
|
|
563 }
|
|
564 }
|
|
565 #endif /* L_unpack_df || L_unpack_sf */
|
|
566
|
|
567 #if defined(L_addsub_sf) || defined(L_addsub_df) || defined(L_addsub_tf)
|
|
568 static const fp_number_type *
|
|
569 _fpadd_parts (fp_number_type * a,
|
|
570 fp_number_type * b,
|
|
571 fp_number_type * tmp)
|
|
572 {
|
|
573 intfrac tfraction;
|
|
574
|
|
575 /* Put commonly used fields in local variables. */
|
|
576 int a_normal_exp;
|
|
577 int b_normal_exp;
|
|
578 fractype a_fraction;
|
|
579 fractype b_fraction;
|
|
580
|
|
581 if (isnan (a))
|
|
582 {
|
|
583 return a;
|
|
584 }
|
|
585 if (isnan (b))
|
|
586 {
|
|
587 return b;
|
|
588 }
|
|
589 if (isinf (a))
|
|
590 {
|
|
591 /* Adding infinities with opposite signs yields a NaN. */
|
|
592 if (isinf (b) && a->sign != b->sign)
|
|
593 return makenan ();
|
|
594 return a;
|
|
595 }
|
|
596 if (isinf (b))
|
|
597 {
|
|
598 return b;
|
|
599 }
|
|
600 if (iszero (b))
|
|
601 {
|
|
602 if (iszero (a))
|
|
603 {
|
|
604 *tmp = *a;
|
|
605 tmp->sign = a->sign & b->sign;
|
|
606 return tmp;
|
|
607 }
|
|
608 return a;
|
|
609 }
|
|
610 if (iszero (a))
|
|
611 {
|
|
612 return b;
|
|
613 }
|
|
614
|
|
615 /* Got two numbers. shift the smaller and increment the exponent till
|
|
616 they're the same */
|
|
617 {
|
|
618 int diff;
|
|
619 int sdiff;
|
|
620
|
|
621 a_normal_exp = a->normal_exp;
|
|
622 b_normal_exp = b->normal_exp;
|
|
623 a_fraction = a->fraction.ll;
|
|
624 b_fraction = b->fraction.ll;
|
|
625
|
|
626 diff = a_normal_exp - b_normal_exp;
|
|
627 sdiff = diff;
|
|
628
|
|
629 if (diff < 0)
|
|
630 diff = -diff;
|
|
631 if (diff < FRAC_NBITS)
|
|
632 {
|
|
633 if (sdiff > 0)
|
|
634 {
|
|
635 b_normal_exp += diff;
|
|
636 LSHIFT (b_fraction, diff);
|
|
637 }
|
|
638 else if (sdiff < 0)
|
|
639 {
|
|
640 a_normal_exp += diff;
|
|
641 LSHIFT (a_fraction, diff);
|
|
642 }
|
|
643 }
|
|
644 else
|
|
645 {
|
|
646 /* Somethings's up.. choose the biggest */
|
|
647 if (a_normal_exp > b_normal_exp)
|
|
648 {
|
|
649 b_normal_exp = a_normal_exp;
|
|
650 b_fraction = 0;
|
|
651 }
|
|
652 else
|
|
653 {
|
|
654 a_normal_exp = b_normal_exp;
|
|
655 a_fraction = 0;
|
|
656 }
|
|
657 }
|
|
658 }
|
|
659
|
|
660 if (a->sign != b->sign)
|
|
661 {
|
|
662 if (a->sign)
|
|
663 {
|
|
664 tfraction = -a_fraction + b_fraction;
|
|
665 }
|
|
666 else
|
|
667 {
|
|
668 tfraction = a_fraction - b_fraction;
|
|
669 }
|
|
670 if (tfraction >= 0)
|
|
671 {
|
|
672 tmp->sign = 0;
|
|
673 tmp->normal_exp = a_normal_exp;
|
|
674 tmp->fraction.ll = tfraction;
|
|
675 }
|
|
676 else
|
|
677 {
|
|
678 tmp->sign = 1;
|
|
679 tmp->normal_exp = a_normal_exp;
|
|
680 tmp->fraction.ll = -tfraction;
|
|
681 }
|
|
682 /* and renormalize it */
|
|
683
|
|
684 while (tmp->fraction.ll < IMPLICIT_1 && tmp->fraction.ll)
|
|
685 {
|
|
686 tmp->fraction.ll <<= 1;
|
|
687 tmp->normal_exp--;
|
|
688 }
|
|
689 }
|
|
690 else
|
|
691 {
|
|
692 tmp->sign = a->sign;
|
|
693 tmp->normal_exp = a_normal_exp;
|
|
694 tmp->fraction.ll = a_fraction + b_fraction;
|
|
695 }
|
|
696 tmp->class = CLASS_NUMBER;
|
|
697 /* Now the fraction is added, we have to shift down to renormalize the
|
|
698 number */
|
|
699
|
|
700 if (tmp->fraction.ll >= IMPLICIT_2)
|
|
701 {
|
|
702 LSHIFT (tmp->fraction.ll, 1);
|
|
703 tmp->normal_exp++;
|
|
704 }
|
|
705 return tmp;
|
|
706 }
|
|
707
|
|
708 FLO_type
|
|
709 add (FLO_type arg_a, FLO_type arg_b)
|
|
710 {
|
|
711 fp_number_type a;
|
|
712 fp_number_type b;
|
|
713 fp_number_type tmp;
|
|
714 const fp_number_type *res;
|
|
715 FLO_union_type au, bu;
|
|
716
|
|
717 au.value = arg_a;
|
|
718 bu.value = arg_b;
|
|
719
|
|
720 unpack_d (&au, &a);
|
|
721 unpack_d (&bu, &b);
|
|
722
|
|
723 res = _fpadd_parts (&a, &b, &tmp);
|
|
724
|
|
725 return pack_d (res);
|
|
726 }
|
|
727
|
|
728 FLO_type
|
|
729 sub (FLO_type arg_a, FLO_type arg_b)
|
|
730 {
|
|
731 fp_number_type a;
|
|
732 fp_number_type b;
|
|
733 fp_number_type tmp;
|
|
734 const fp_number_type *res;
|
|
735 FLO_union_type au, bu;
|
|
736
|
|
737 au.value = arg_a;
|
|
738 bu.value = arg_b;
|
|
739
|
|
740 unpack_d (&au, &a);
|
|
741 unpack_d (&bu, &b);
|
|
742
|
|
743 b.sign ^= 1;
|
|
744
|
|
745 res = _fpadd_parts (&a, &b, &tmp);
|
|
746
|
|
747 return pack_d (res);
|
|
748 }
|
|
749 #endif /* L_addsub_sf || L_addsub_df */
|
|
750
|
|
751 #if defined(L_mul_sf) || defined(L_mul_df) || defined(L_mul_tf)
|
|
752 static inline __attribute__ ((__always_inline__)) const fp_number_type *
|
|
753 _fpmul_parts ( fp_number_type * a,
|
|
754 fp_number_type * b,
|
|
755 fp_number_type * tmp)
|
|
756 {
|
|
757 fractype low = 0;
|
|
758 fractype high = 0;
|
|
759
|
|
760 if (isnan (a))
|
|
761 {
|
|
762 a->sign = a->sign != b->sign;
|
|
763 return a;
|
|
764 }
|
|
765 if (isnan (b))
|
|
766 {
|
|
767 b->sign = a->sign != b->sign;
|
|
768 return b;
|
|
769 }
|
|
770 if (isinf (a))
|
|
771 {
|
|
772 if (iszero (b))
|
|
773 return makenan ();
|
|
774 a->sign = a->sign != b->sign;
|
|
775 return a;
|
|
776 }
|
|
777 if (isinf (b))
|
|
778 {
|
|
779 if (iszero (a))
|
|
780 {
|
|
781 return makenan ();
|
|
782 }
|
|
783 b->sign = a->sign != b->sign;
|
|
784 return b;
|
|
785 }
|
|
786 if (iszero (a))
|
|
787 {
|
|
788 a->sign = a->sign != b->sign;
|
|
789 return a;
|
|
790 }
|
|
791 if (iszero (b))
|
|
792 {
|
|
793 b->sign = a->sign != b->sign;
|
|
794 return b;
|
|
795 }
|
|
796
|
|
797 /* Calculate the mantissa by multiplying both numbers to get a
|
|
798 twice-as-wide number. */
|
|
799 {
|
|
800 #if defined(NO_DI_MODE) || defined(TFLOAT)
|
|
801 {
|
|
802 fractype x = a->fraction.ll;
|
|
803 fractype ylow = b->fraction.ll;
|
|
804 fractype yhigh = 0;
|
|
805 int bit;
|
|
806
|
|
807 /* ??? This does multiplies one bit at a time. Optimize. */
|
|
808 for (bit = 0; bit < FRAC_NBITS; bit++)
|
|
809 {
|
|
810 int carry;
|
|
811
|
|
812 if (x & 1)
|
|
813 {
|
|
814 carry = (low += ylow) < ylow;
|
|
815 high += yhigh + carry;
|
|
816 }
|
|
817 yhigh <<= 1;
|
|
818 if (ylow & FRACHIGH)
|
|
819 {
|
|
820 yhigh |= 1;
|
|
821 }
|
|
822 ylow <<= 1;
|
|
823 x >>= 1;
|
|
824 }
|
|
825 }
|
|
826 #elif defined(FLOAT)
|
|
827 /* Multiplying two USIs to get a UDI, we're safe. */
|
|
828 {
|
|
829 UDItype answer = (UDItype)a->fraction.ll * (UDItype)b->fraction.ll;
|
|
830
|
|
831 high = answer >> BITS_PER_SI;
|
|
832 low = answer;
|
|
833 }
|
|
834 #else
|
|
835 /* fractype is DImode, but we need the result to be twice as wide.
|
|
836 Assuming a widening multiply from DImode to TImode is not
|
|
837 available, build one by hand. */
|
|
838 {
|
|
839 USItype nl = a->fraction.ll;
|
|
840 USItype nh = a->fraction.ll >> BITS_PER_SI;
|
|
841 USItype ml = b->fraction.ll;
|
|
842 USItype mh = b->fraction.ll >> BITS_PER_SI;
|
|
843 UDItype pp_ll = (UDItype) ml * nl;
|
|
844 UDItype pp_hl = (UDItype) mh * nl;
|
|
845 UDItype pp_lh = (UDItype) ml * nh;
|
|
846 UDItype pp_hh = (UDItype) mh * nh;
|
|
847 UDItype res2 = 0;
|
|
848 UDItype res0 = 0;
|
|
849 UDItype ps_hh__ = pp_hl + pp_lh;
|
|
850 if (ps_hh__ < pp_hl)
|
|
851 res2 += (UDItype)1 << BITS_PER_SI;
|
|
852 pp_hl = (UDItype)(USItype)ps_hh__ << BITS_PER_SI;
|
|
853 res0 = pp_ll + pp_hl;
|
|
854 if (res0 < pp_ll)
|
|
855 res2++;
|
|
856 res2 += (ps_hh__ >> BITS_PER_SI) + pp_hh;
|
|
857 high = res2;
|
|
858 low = res0;
|
|
859 }
|
|
860 #endif
|
|
861 }
|
|
862
|
|
863 tmp->normal_exp = a->normal_exp + b->normal_exp
|
|
864 + FRAC_NBITS - (FRACBITS + NGARDS);
|
|
865 tmp->sign = a->sign != b->sign;
|
|
866 while (high >= IMPLICIT_2)
|
|
867 {
|
|
868 tmp->normal_exp++;
|
|
869 if (high & 1)
|
|
870 {
|
|
871 low >>= 1;
|
|
872 low |= FRACHIGH;
|
|
873 }
|
|
874 high >>= 1;
|
|
875 }
|
|
876 while (high < IMPLICIT_1)
|
|
877 {
|
|
878 tmp->normal_exp--;
|
|
879
|
|
880 high <<= 1;
|
|
881 if (low & FRACHIGH)
|
|
882 high |= 1;
|
|
883 low <<= 1;
|
|
884 }
|
|
885
|
|
886 if ((high & GARDMASK) == GARDMSB)
|
|
887 {
|
|
888 if (high & (1 << NGARDS))
|
|
889 {
|
|
890 /* Because we're half way, we would round to even by adding
|
|
891 GARDROUND + 1, except that's also done in the packing
|
|
892 function, and rounding twice will lose precision and cause
|
|
893 the result to be too far off. Example: 32-bit floats with
|
|
894 bit patterns 0xfff * 0x3f800400 ~= 0xfff (less than 0.5ulp
|
|
895 off), not 0x1000 (more than 0.5ulp off). */
|
|
896 }
|
|
897 else if (low)
|
|
898 {
|
|
899 /* We're a further than half way by a small amount corresponding
|
|
900 to the bits set in "low". Knowing that, we round here and
|
|
901 not in pack_d, because there we don't have "low" available
|
|
902 anymore. */
|
|
903 high += GARDROUND + 1;
|
|
904
|
|
905 /* Avoid further rounding in pack_d. */
|
|
906 high &= ~(fractype) GARDMASK;
|
|
907 }
|
|
908 }
|
|
909 tmp->fraction.ll = high;
|
|
910 tmp->class = CLASS_NUMBER;
|
|
911 return tmp;
|
|
912 }
|
|
913
|
|
914 FLO_type
|
|
915 multiply (FLO_type arg_a, FLO_type arg_b)
|
|
916 {
|
|
917 fp_number_type a;
|
|
918 fp_number_type b;
|
|
919 fp_number_type tmp;
|
|
920 const fp_number_type *res;
|
|
921 FLO_union_type au, bu;
|
|
922
|
|
923 au.value = arg_a;
|
|
924 bu.value = arg_b;
|
|
925
|
|
926 unpack_d (&au, &a);
|
|
927 unpack_d (&bu, &b);
|
|
928
|
|
929 res = _fpmul_parts (&a, &b, &tmp);
|
|
930
|
|
931 return pack_d (res);
|
|
932 }
|
|
933 #endif /* L_mul_sf || L_mul_df || L_mul_tf */
|
|
934
|
|
935 #if defined(L_div_sf) || defined(L_div_df) || defined(L_div_tf)
|
|
936 static inline __attribute__ ((__always_inline__)) const fp_number_type *
|
|
937 _fpdiv_parts (fp_number_type * a,
|
|
938 fp_number_type * b)
|
|
939 {
|
|
940 fractype bit;
|
|
941 fractype numerator;
|
|
942 fractype denominator;
|
|
943 fractype quotient;
|
|
944
|
|
945 if (isnan (a))
|
|
946 {
|
|
947 return a;
|
|
948 }
|
|
949 if (isnan (b))
|
|
950 {
|
|
951 return b;
|
|
952 }
|
|
953
|
|
954 a->sign = a->sign ^ b->sign;
|
|
955
|
|
956 if (isinf (a) || iszero (a))
|
|
957 {
|
|
958 if (a->class == b->class)
|
|
959 return makenan ();
|
|
960 return a;
|
|
961 }
|
|
962
|
|
963 if (isinf (b))
|
|
964 {
|
|
965 a->fraction.ll = 0;
|
|
966 a->normal_exp = 0;
|
|
967 return a;
|
|
968 }
|
|
969 if (iszero (b))
|
|
970 {
|
|
971 a->class = CLASS_INFINITY;
|
|
972 return a;
|
|
973 }
|
|
974
|
|
975 /* Calculate the mantissa by multiplying both 64bit numbers to get a
|
|
976 128 bit number */
|
|
977 {
|
|
978 /* quotient =
|
|
979 ( numerator / denominator) * 2^(numerator exponent - denominator exponent)
|
|
980 */
|
|
981
|
|
982 a->normal_exp = a->normal_exp - b->normal_exp;
|
|
983 numerator = a->fraction.ll;
|
|
984 denominator = b->fraction.ll;
|
|
985
|
|
986 if (numerator < denominator)
|
|
987 {
|
|
988 /* Fraction will be less than 1.0 */
|
|
989 numerator *= 2;
|
|
990 a->normal_exp--;
|
|
991 }
|
|
992 bit = IMPLICIT_1;
|
|
993 quotient = 0;
|
|
994 /* ??? Does divide one bit at a time. Optimize. */
|
|
995 while (bit)
|
|
996 {
|
|
997 if (numerator >= denominator)
|
|
998 {
|
|
999 quotient |= bit;
|
|
1000 numerator -= denominator;
|
|
1001 }
|
|
1002 bit >>= 1;
|
|
1003 numerator *= 2;
|
|
1004 }
|
|
1005
|
|
1006 if ((quotient & GARDMASK) == GARDMSB)
|
|
1007 {
|
|
1008 if (quotient & (1 << NGARDS))
|
|
1009 {
|
|
1010 /* Because we're half way, we would round to even by adding
|
|
1011 GARDROUND + 1, except that's also done in the packing
|
|
1012 function, and rounding twice will lose precision and cause
|
|
1013 the result to be too far off. */
|
|
1014 }
|
|
1015 else if (numerator)
|
|
1016 {
|
|
1017 /* We're a further than half way by the small amount
|
|
1018 corresponding to the bits set in "numerator". Knowing
|
|
1019 that, we round here and not in pack_d, because there we
|
|
1020 don't have "numerator" available anymore. */
|
|
1021 quotient += GARDROUND + 1;
|
|
1022
|
|
1023 /* Avoid further rounding in pack_d. */
|
|
1024 quotient &= ~(fractype) GARDMASK;
|
|
1025 }
|
|
1026 }
|
|
1027
|
|
1028 a->fraction.ll = quotient;
|
|
1029 return (a);
|
|
1030 }
|
|
1031 }
|
|
1032
|
|
1033 FLO_type
|
|
1034 divide (FLO_type arg_a, FLO_type arg_b)
|
|
1035 {
|
|
1036 fp_number_type a;
|
|
1037 fp_number_type b;
|
|
1038 const fp_number_type *res;
|
|
1039 FLO_union_type au, bu;
|
|
1040
|
|
1041 au.value = arg_a;
|
|
1042 bu.value = arg_b;
|
|
1043
|
|
1044 unpack_d (&au, &a);
|
|
1045 unpack_d (&bu, &b);
|
|
1046
|
|
1047 res = _fpdiv_parts (&a, &b);
|
|
1048
|
|
1049 return pack_d (res);
|
|
1050 }
|
|
1051 #endif /* L_div_sf || L_div_df */
|
|
1052
|
|
1053 #if defined(L_fpcmp_parts_sf) || defined(L_fpcmp_parts_df) \
|
|
1054 || defined(L_fpcmp_parts_tf)
|
|
1055 /* according to the demo, fpcmp returns a comparison with 0... thus
|
|
1056 a<b -> -1
|
|
1057 a==b -> 0
|
|
1058 a>b -> +1
|
|
1059 */
|
|
1060
|
|
1061 int
|
|
1062 __fpcmp_parts (fp_number_type * a, fp_number_type * b)
|
|
1063 {
|
|
1064 #if 0
|
|
1065 /* either nan -> unordered. Must be checked outside of this routine. */
|
|
1066 if (isnan (a) && isnan (b))
|
|
1067 {
|
|
1068 return 1; /* still unordered! */
|
|
1069 }
|
|
1070 #endif
|
|
1071
|
|
1072 if (isnan (a) || isnan (b))
|
|
1073 {
|
|
1074 return 1; /* how to indicate unordered compare? */
|
|
1075 }
|
|
1076 if (isinf (a) && isinf (b))
|
|
1077 {
|
|
1078 /* +inf > -inf, but +inf != +inf */
|
|
1079 /* b \a| +inf(0)| -inf(1)
|
|
1080 ______\+--------+--------
|
|
1081 +inf(0)| a==b(0)| a<b(-1)
|
|
1082 -------+--------+--------
|
|
1083 -inf(1)| a>b(1) | a==b(0)
|
|
1084 -------+--------+--------
|
|
1085 So since unordered must be nonzero, just line up the columns...
|
|
1086 */
|
|
1087 return b->sign - a->sign;
|
|
1088 }
|
|
1089 /* but not both... */
|
|
1090 if (isinf (a))
|
|
1091 {
|
|
1092 return a->sign ? -1 : 1;
|
|
1093 }
|
|
1094 if (isinf (b))
|
|
1095 {
|
|
1096 return b->sign ? 1 : -1;
|
|
1097 }
|
|
1098 if (iszero (a) && iszero (b))
|
|
1099 {
|
|
1100 return 0;
|
|
1101 }
|
|
1102 if (iszero (a))
|
|
1103 {
|
|
1104 return b->sign ? 1 : -1;
|
|
1105 }
|
|
1106 if (iszero (b))
|
|
1107 {
|
|
1108 return a->sign ? -1 : 1;
|
|
1109 }
|
|
1110 /* now both are "normal". */
|
|
1111 if (a->sign != b->sign)
|
|
1112 {
|
|
1113 /* opposite signs */
|
|
1114 return a->sign ? -1 : 1;
|
|
1115 }
|
|
1116 /* same sign; exponents? */
|
|
1117 if (a->normal_exp > b->normal_exp)
|
|
1118 {
|
|
1119 return a->sign ? -1 : 1;
|
|
1120 }
|
|
1121 if (a->normal_exp < b->normal_exp)
|
|
1122 {
|
|
1123 return a->sign ? 1 : -1;
|
|
1124 }
|
|
1125 /* same exponents; check size. */
|
|
1126 if (a->fraction.ll > b->fraction.ll)
|
|
1127 {
|
|
1128 return a->sign ? -1 : 1;
|
|
1129 }
|
|
1130 if (a->fraction.ll < b->fraction.ll)
|
|
1131 {
|
|
1132 return a->sign ? 1 : -1;
|
|
1133 }
|
|
1134 /* after all that, they're equal. */
|
|
1135 return 0;
|
|
1136 }
|
|
1137 #endif
|
|
1138
|
|
1139 #if defined(L_compare_sf) || defined(L_compare_df) || defined(L_compoare_tf)
|
|
1140 CMPtype
|
|
1141 compare (FLO_type arg_a, FLO_type arg_b)
|
|
1142 {
|
|
1143 fp_number_type a;
|
|
1144 fp_number_type b;
|
|
1145 FLO_union_type au, bu;
|
|
1146
|
|
1147 au.value = arg_a;
|
|
1148 bu.value = arg_b;
|
|
1149
|
|
1150 unpack_d (&au, &a);
|
|
1151 unpack_d (&bu, &b);
|
|
1152
|
|
1153 return __fpcmp_parts (&a, &b);
|
|
1154 }
|
|
1155 #endif /* L_compare_sf || L_compare_df */
|
|
1156
|
|
1157 /* These should be optimized for their specific tasks someday. */
|
|
1158
|
|
1159 #if defined(L_eq_sf) || defined(L_eq_df) || defined(L_eq_tf)
|
|
1160 CMPtype
|
|
1161 _eq_f2 (FLO_type arg_a, FLO_type arg_b)
|
|
1162 {
|
|
1163 fp_number_type a;
|
|
1164 fp_number_type b;
|
|
1165 FLO_union_type au, bu;
|
|
1166
|
|
1167 au.value = arg_a;
|
|
1168 bu.value = arg_b;
|
|
1169
|
|
1170 unpack_d (&au, &a);
|
|
1171 unpack_d (&bu, &b);
|
|
1172
|
|
1173 if (isnan (&a) || isnan (&b))
|
|
1174 return 1; /* false, truth == 0 */
|
|
1175
|
|
1176 return __fpcmp_parts (&a, &b) ;
|
|
1177 }
|
|
1178 #endif /* L_eq_sf || L_eq_df */
|
|
1179
|
|
1180 #if defined(L_ne_sf) || defined(L_ne_df) || defined(L_ne_tf)
|
|
1181 CMPtype
|
|
1182 _ne_f2 (FLO_type arg_a, FLO_type arg_b)
|
|
1183 {
|
|
1184 fp_number_type a;
|
|
1185 fp_number_type b;
|
|
1186 FLO_union_type au, bu;
|
|
1187
|
|
1188 au.value = arg_a;
|
|
1189 bu.value = arg_b;
|
|
1190
|
|
1191 unpack_d (&au, &a);
|
|
1192 unpack_d (&bu, &b);
|
|
1193
|
|
1194 if (isnan (&a) || isnan (&b))
|
|
1195 return 1; /* true, truth != 0 */
|
|
1196
|
|
1197 return __fpcmp_parts (&a, &b) ;
|
|
1198 }
|
|
1199 #endif /* L_ne_sf || L_ne_df */
|
|
1200
|
|
1201 #if defined(L_gt_sf) || defined(L_gt_df) || defined(L_gt_tf)
|
|
1202 CMPtype
|
|
1203 _gt_f2 (FLO_type arg_a, FLO_type arg_b)
|
|
1204 {
|
|
1205 fp_number_type a;
|
|
1206 fp_number_type b;
|
|
1207 FLO_union_type au, bu;
|
|
1208
|
|
1209 au.value = arg_a;
|
|
1210 bu.value = arg_b;
|
|
1211
|
|
1212 unpack_d (&au, &a);
|
|
1213 unpack_d (&bu, &b);
|
|
1214
|
|
1215 if (isnan (&a) || isnan (&b))
|
|
1216 return -1; /* false, truth > 0 */
|
|
1217
|
|
1218 return __fpcmp_parts (&a, &b);
|
|
1219 }
|
|
1220 #endif /* L_gt_sf || L_gt_df */
|
|
1221
|
|
1222 #if defined(L_ge_sf) || defined(L_ge_df) || defined(L_ge_tf)
|
|
1223 CMPtype
|
|
1224 _ge_f2 (FLO_type arg_a, FLO_type arg_b)
|
|
1225 {
|
|
1226 fp_number_type a;
|
|
1227 fp_number_type b;
|
|
1228 FLO_union_type au, bu;
|
|
1229
|
|
1230 au.value = arg_a;
|
|
1231 bu.value = arg_b;
|
|
1232
|
|
1233 unpack_d (&au, &a);
|
|
1234 unpack_d (&bu, &b);
|
|
1235
|
|
1236 if (isnan (&a) || isnan (&b))
|
|
1237 return -1; /* false, truth >= 0 */
|
|
1238 return __fpcmp_parts (&a, &b) ;
|
|
1239 }
|
|
1240 #endif /* L_ge_sf || L_ge_df */
|
|
1241
|
|
1242 #if defined(L_lt_sf) || defined(L_lt_df) || defined(L_lt_tf)
|
|
1243 CMPtype
|
|
1244 _lt_f2 (FLO_type arg_a, FLO_type arg_b)
|
|
1245 {
|
|
1246 fp_number_type a;
|
|
1247 fp_number_type b;
|
|
1248 FLO_union_type au, bu;
|
|
1249
|
|
1250 au.value = arg_a;
|
|
1251 bu.value = arg_b;
|
|
1252
|
|
1253 unpack_d (&au, &a);
|
|
1254 unpack_d (&bu, &b);
|
|
1255
|
|
1256 if (isnan (&a) || isnan (&b))
|
|
1257 return 1; /* false, truth < 0 */
|
|
1258
|
|
1259 return __fpcmp_parts (&a, &b);
|
|
1260 }
|
|
1261 #endif /* L_lt_sf || L_lt_df */
|
|
1262
|
|
1263 #if defined(L_le_sf) || defined(L_le_df) || defined(L_le_tf)
|
|
1264 CMPtype
|
|
1265 _le_f2 (FLO_type arg_a, FLO_type arg_b)
|
|
1266 {
|
|
1267 fp_number_type a;
|
|
1268 fp_number_type b;
|
|
1269 FLO_union_type au, bu;
|
|
1270
|
|
1271 au.value = arg_a;
|
|
1272 bu.value = arg_b;
|
|
1273
|
|
1274 unpack_d (&au, &a);
|
|
1275 unpack_d (&bu, &b);
|
|
1276
|
|
1277 if (isnan (&a) || isnan (&b))
|
|
1278 return 1; /* false, truth <= 0 */
|
|
1279
|
|
1280 return __fpcmp_parts (&a, &b) ;
|
|
1281 }
|
|
1282 #endif /* L_le_sf || L_le_df */
|
|
1283
|
|
1284 #if defined(L_unord_sf) || defined(L_unord_df) || defined(L_unord_tf)
|
|
1285 CMPtype
|
|
1286 _unord_f2 (FLO_type arg_a, FLO_type arg_b)
|
|
1287 {
|
|
1288 fp_number_type a;
|
|
1289 fp_number_type b;
|
|
1290 FLO_union_type au, bu;
|
|
1291
|
|
1292 au.value = arg_a;
|
|
1293 bu.value = arg_b;
|
|
1294
|
|
1295 unpack_d (&au, &a);
|
|
1296 unpack_d (&bu, &b);
|
|
1297
|
|
1298 return (isnan (&a) || isnan (&b));
|
|
1299 }
|
|
1300 #endif /* L_unord_sf || L_unord_df */
|
|
1301
|
|
1302 #if defined(L_si_to_sf) || defined(L_si_to_df) || defined(L_si_to_tf)
|
|
1303 FLO_type
|
|
1304 si_to_float (SItype arg_a)
|
|
1305 {
|
|
1306 fp_number_type in;
|
|
1307
|
|
1308 in.class = CLASS_NUMBER;
|
|
1309 in.sign = arg_a < 0;
|
|
1310 if (!arg_a)
|
|
1311 {
|
|
1312 in.class = CLASS_ZERO;
|
|
1313 }
|
|
1314 else
|
|
1315 {
|
|
1316 USItype uarg;
|
|
1317 int shift;
|
|
1318 in.normal_exp = FRACBITS + NGARDS;
|
|
1319 if (in.sign)
|
|
1320 {
|
|
1321 /* Special case for minint, since there is no +ve integer
|
|
1322 representation for it */
|
|
1323 if (arg_a == (- MAX_SI_INT - 1))
|
|
1324 {
|
|
1325 return (FLO_type)(- MAX_SI_INT - 1);
|
|
1326 }
|
|
1327 uarg = (-arg_a);
|
|
1328 }
|
|
1329 else
|
|
1330 uarg = arg_a;
|
|
1331
|
|
1332 in.fraction.ll = uarg;
|
|
1333 shift = clzusi (uarg) - (BITS_PER_SI - 1 - FRACBITS - NGARDS);
|
|
1334 if (shift > 0)
|
|
1335 {
|
|
1336 in.fraction.ll <<= shift;
|
|
1337 in.normal_exp -= shift;
|
|
1338 }
|
|
1339 }
|
|
1340 return pack_d (&in);
|
|
1341 }
|
|
1342 #endif /* L_si_to_sf || L_si_to_df */
|
|
1343
|
|
1344 #if defined(L_usi_to_sf) || defined(L_usi_to_df) || defined(L_usi_to_tf)
|
|
1345 FLO_type
|
|
1346 usi_to_float (USItype arg_a)
|
|
1347 {
|
|
1348 fp_number_type in;
|
|
1349
|
|
1350 in.sign = 0;
|
|
1351 if (!arg_a)
|
|
1352 {
|
|
1353 in.class = CLASS_ZERO;
|
|
1354 }
|
|
1355 else
|
|
1356 {
|
|
1357 int shift;
|
|
1358 in.class = CLASS_NUMBER;
|
|
1359 in.normal_exp = FRACBITS + NGARDS;
|
|
1360 in.fraction.ll = arg_a;
|
|
1361
|
|
1362 shift = clzusi (arg_a) - (BITS_PER_SI - 1 - FRACBITS - NGARDS);
|
|
1363 if (shift < 0)
|
|
1364 {
|
|
1365 fractype guard = in.fraction.ll & (((fractype)1 << -shift) - 1);
|
|
1366 in.fraction.ll >>= -shift;
|
|
1367 in.fraction.ll |= (guard != 0);
|
|
1368 in.normal_exp -= shift;
|
|
1369 }
|
|
1370 else if (shift > 0)
|
|
1371 {
|
|
1372 in.fraction.ll <<= shift;
|
|
1373 in.normal_exp -= shift;
|
|
1374 }
|
|
1375 }
|
|
1376 return pack_d (&in);
|
|
1377 }
|
|
1378 #endif
|
|
1379
|
|
1380 #if defined(L_sf_to_si) || defined(L_df_to_si) || defined(L_tf_to_si)
|
|
1381 SItype
|
|
1382 float_to_si (FLO_type arg_a)
|
|
1383 {
|
|
1384 fp_number_type a;
|
|
1385 SItype tmp;
|
|
1386 FLO_union_type au;
|
|
1387
|
|
1388 au.value = arg_a;
|
|
1389 unpack_d (&au, &a);
|
|
1390
|
|
1391 if (iszero (&a))
|
|
1392 return 0;
|
|
1393 if (isnan (&a))
|
|
1394 return 0;
|
|
1395 /* get reasonable MAX_SI_INT... */
|
|
1396 if (isinf (&a))
|
|
1397 return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT;
|
|
1398 /* it is a number, but a small one */
|
|
1399 if (a.normal_exp < 0)
|
|
1400 return 0;
|
|
1401 if (a.normal_exp > BITS_PER_SI - 2)
|
|
1402 return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT;
|
|
1403 tmp = a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
|
|
1404 return a.sign ? (-tmp) : (tmp);
|
|
1405 }
|
|
1406 #endif /* L_sf_to_si || L_df_to_si */
|
|
1407
|
|
1408 #if defined(L_tf_to_usi)
|
|
1409 USItype
|
|
1410 float_to_usi (FLO_type arg_a)
|
|
1411 {
|
|
1412 fp_number_type a;
|
|
1413 FLO_union_type au;
|
|
1414
|
|
1415 au.value = arg_a;
|
|
1416 unpack_d (&au, &a);
|
|
1417
|
|
1418 if (iszero (&a))
|
|
1419 return 0;
|
|
1420 if (isnan (&a))
|
|
1421 return 0;
|
|
1422 /* it is a negative number */
|
|
1423 if (a.sign)
|
|
1424 return 0;
|
|
1425 /* get reasonable MAX_USI_INT... */
|
|
1426 if (isinf (&a))
|
|
1427 return MAX_USI_INT;
|
|
1428 /* it is a number, but a small one */
|
|
1429 if (a.normal_exp < 0)
|
|
1430 return 0;
|
|
1431 if (a.normal_exp > BITS_PER_SI - 1)
|
|
1432 return MAX_USI_INT;
|
|
1433 else if (a.normal_exp > (FRACBITS + NGARDS))
|
|
1434 return a.fraction.ll << (a.normal_exp - (FRACBITS + NGARDS));
|
|
1435 else
|
|
1436 return a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
|
|
1437 }
|
|
1438 #endif /* L_tf_to_usi */
|
|
1439
|
|
1440 #if defined(L_negate_sf) || defined(L_negate_df) || defined(L_negate_tf)
|
|
1441 FLO_type
|
|
1442 negate (FLO_type arg_a)
|
|
1443 {
|
|
1444 fp_number_type a;
|
|
1445 FLO_union_type au;
|
|
1446
|
|
1447 au.value = arg_a;
|
|
1448 unpack_d (&au, &a);
|
|
1449
|
|
1450 flip_sign (&a);
|
|
1451 return pack_d (&a);
|
|
1452 }
|
|
1453 #endif /* L_negate_sf || L_negate_df */
|
|
1454
|
|
1455 #ifdef FLOAT
|
|
1456
|
|
1457 #if defined(L_make_sf)
|
|
1458 SFtype
|
|
1459 __make_fp(fp_class_type class,
|
|
1460 unsigned int sign,
|
|
1461 int exp,
|
|
1462 USItype frac)
|
|
1463 {
|
|
1464 fp_number_type in;
|
|
1465
|
|
1466 in.class = class;
|
|
1467 in.sign = sign;
|
|
1468 in.normal_exp = exp;
|
|
1469 in.fraction.ll = frac;
|
|
1470 return pack_d (&in);
|
|
1471 }
|
|
1472 #endif /* L_make_sf */
|
|
1473
|
|
1474 #ifndef FLOAT_ONLY
|
|
1475
|
|
1476 /* This enables one to build an fp library that supports float but not double.
|
|
1477 Otherwise, we would get an undefined reference to __make_dp.
|
|
1478 This is needed for some 8-bit ports that can't handle well values that
|
|
1479 are 8-bytes in size, so we just don't support double for them at all. */
|
|
1480
|
|
1481 #if defined(L_sf_to_df)
|
|
1482 DFtype
|
|
1483 sf_to_df (SFtype arg_a)
|
|
1484 {
|
|
1485 fp_number_type in;
|
|
1486 FLO_union_type au;
|
|
1487
|
|
1488 au.value = arg_a;
|
|
1489 unpack_d (&au, &in);
|
|
1490
|
|
1491 return __make_dp (in.class, in.sign, in.normal_exp,
|
|
1492 ((UDItype) in.fraction.ll) << F_D_BITOFF);
|
|
1493 }
|
|
1494 #endif /* L_sf_to_df */
|
|
1495
|
|
1496 #if defined(L_sf_to_tf) && defined(TMODES)
|
|
1497 TFtype
|
|
1498 sf_to_tf (SFtype arg_a)
|
|
1499 {
|
|
1500 fp_number_type in;
|
|
1501 FLO_union_type au;
|
|
1502
|
|
1503 au.value = arg_a;
|
|
1504 unpack_d (&au, &in);
|
|
1505
|
|
1506 return __make_tp (in.class, in.sign, in.normal_exp,
|
|
1507 ((UTItype) in.fraction.ll) << F_T_BITOFF);
|
|
1508 }
|
|
1509 #endif /* L_sf_to_df */
|
|
1510
|
|
1511 #endif /* ! FLOAT_ONLY */
|
|
1512 #endif /* FLOAT */
|
|
1513
|
|
1514 #ifndef FLOAT
|
|
1515
|
|
1516 extern SFtype __make_fp (fp_class_type, unsigned int, int, USItype);
|
|
1517
|
|
1518 #if defined(L_make_df)
|
|
1519 DFtype
|
|
1520 __make_dp (fp_class_type class, unsigned int sign, int exp, UDItype frac)
|
|
1521 {
|
|
1522 fp_number_type in;
|
|
1523
|
|
1524 in.class = class;
|
|
1525 in.sign = sign;
|
|
1526 in.normal_exp = exp;
|
|
1527 in.fraction.ll = frac;
|
|
1528 return pack_d (&in);
|
|
1529 }
|
|
1530 #endif /* L_make_df */
|
|
1531
|
|
1532 #if defined(L_df_to_sf)
|
|
1533 SFtype
|
|
1534 df_to_sf (DFtype arg_a)
|
|
1535 {
|
|
1536 fp_number_type in;
|
|
1537 USItype sffrac;
|
|
1538 FLO_union_type au;
|
|
1539
|
|
1540 au.value = arg_a;
|
|
1541 unpack_d (&au, &in);
|
|
1542
|
|
1543 sffrac = in.fraction.ll >> F_D_BITOFF;
|
|
1544
|
|
1545 /* We set the lowest guard bit in SFFRAC if we discarded any non
|
|
1546 zero bits. */
|
|
1547 if ((in.fraction.ll & (((USItype) 1 << F_D_BITOFF) - 1)) != 0)
|
|
1548 sffrac |= 1;
|
|
1549
|
|
1550 return __make_fp (in.class, in.sign, in.normal_exp, sffrac);
|
|
1551 }
|
|
1552 #endif /* L_df_to_sf */
|
|
1553
|
|
1554 #if defined(L_df_to_tf) && defined(TMODES) \
|
|
1555 && !defined(FLOAT) && !defined(TFLOAT)
|
|
1556 TFtype
|
|
1557 df_to_tf (DFtype arg_a)
|
|
1558 {
|
|
1559 fp_number_type in;
|
|
1560 FLO_union_type au;
|
|
1561
|
|
1562 au.value = arg_a;
|
|
1563 unpack_d (&au, &in);
|
|
1564
|
|
1565 return __make_tp (in.class, in.sign, in.normal_exp,
|
|
1566 ((UTItype) in.fraction.ll) << D_T_BITOFF);
|
|
1567 }
|
|
1568 #endif /* L_sf_to_df */
|
|
1569
|
|
1570 #ifdef TFLOAT
|
|
1571 #if defined(L_make_tf)
|
|
1572 TFtype
|
|
1573 __make_tp(fp_class_type class,
|
|
1574 unsigned int sign,
|
|
1575 int exp,
|
|
1576 UTItype frac)
|
|
1577 {
|
|
1578 fp_number_type in;
|
|
1579
|
|
1580 in.class = class;
|
|
1581 in.sign = sign;
|
|
1582 in.normal_exp = exp;
|
|
1583 in.fraction.ll = frac;
|
|
1584 return pack_d (&in);
|
|
1585 }
|
|
1586 #endif /* L_make_tf */
|
|
1587
|
|
1588 #if defined(L_tf_to_df)
|
|
1589 DFtype
|
|
1590 tf_to_df (TFtype arg_a)
|
|
1591 {
|
|
1592 fp_number_type in;
|
|
1593 UDItype sffrac;
|
|
1594 FLO_union_type au;
|
|
1595
|
|
1596 au.value = arg_a;
|
|
1597 unpack_d (&au, &in);
|
|
1598
|
|
1599 sffrac = in.fraction.ll >> D_T_BITOFF;
|
|
1600
|
|
1601 /* We set the lowest guard bit in SFFRAC if we discarded any non
|
|
1602 zero bits. */
|
|
1603 if ((in.fraction.ll & (((UTItype) 1 << D_T_BITOFF) - 1)) != 0)
|
|
1604 sffrac |= 1;
|
|
1605
|
|
1606 return __make_dp (in.class, in.sign, in.normal_exp, sffrac);
|
|
1607 }
|
|
1608 #endif /* L_tf_to_df */
|
|
1609
|
|
1610 #if defined(L_tf_to_sf)
|
|
1611 SFtype
|
|
1612 tf_to_sf (TFtype arg_a)
|
|
1613 {
|
|
1614 fp_number_type in;
|
|
1615 USItype sffrac;
|
|
1616 FLO_union_type au;
|
|
1617
|
|
1618 au.value = arg_a;
|
|
1619 unpack_d (&au, &in);
|
|
1620
|
|
1621 sffrac = in.fraction.ll >> F_T_BITOFF;
|
|
1622
|
|
1623 /* We set the lowest guard bit in SFFRAC if we discarded any non
|
|
1624 zero bits. */
|
|
1625 if ((in.fraction.ll & (((UTItype) 1 << F_T_BITOFF) - 1)) != 0)
|
|
1626 sffrac |= 1;
|
|
1627
|
|
1628 return __make_fp (in.class, in.sign, in.normal_exp, sffrac);
|
|
1629 }
|
|
1630 #endif /* L_tf_to_sf */
|
|
1631 #endif /* TFLOAT */
|
|
1632
|
|
1633 #endif /* ! FLOAT */
|
|
1634 #endif /* !EXTENDED_FLOAT_STUBS */
|