Daily bump.
[official-gcc.git] / libgcc / fp-bit.c
blobb3ab7bad52dc341096b7a1c2f3a3ac95ab79ef43
1 /* This is a software floating point library which can be used
2 for targets without hardware floating point.
3 Copyright (C) 1994-2021 Free Software Foundation, Inc.
5 This file is part of GCC.
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.
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.
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.
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/>. */
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.
30 The original code by Steve Chamberlain, hacked by Mark Eichin and Jim
31 Wilson, all of Cygnus Support. */
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. */
36 #include "tconfig.h"
37 #include "coretypes.h"
38 #include "tm.h"
39 #include "libgcc_tm.h"
40 #include "fp-bit.h"
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 */
58 /* We don't currently support extended floats (long doubles) on machines
59 without hardware to deal with them.
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. */
68 #ifdef DECLARE_LIBRARY_RENAMES
69 DECLARE_LIBRARY_RENAMES
70 #endif
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(); }
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 */
111 /* IEEE "special" number predicates */
113 #ifdef NO_NANS
115 #define nan() 0
116 #define isnan(x) 0
117 #define isinf(x) 0
118 #else
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
134 INLINE
135 static const fp_number_type *
136 makenan (void)
138 #ifdef TFLOAT
139 return & __thenan_tf;
140 #elif defined FLOAT
141 return & __thenan_sf;
142 #else
143 return & __thenan_df;
144 #endif
147 INLINE
148 static int
149 isnan (const fp_number_type *x)
151 return __builtin_expect (x->class == CLASS_SNAN || x->class == CLASS_QNAN,
155 INLINE
156 static int
157 isinf (const fp_number_type * x)
159 return __builtin_expect (x->class == CLASS_INFINITY, 0);
162 #endif /* NO_NANS */
164 INLINE
165 static int
166 iszero (const fp_number_type * x)
168 return x->class == CLASS_ZERO;
171 INLINE
172 static void
173 flip_sign ( fp_number_type * x)
175 x->sign = !x->sign;
178 /* Count leading zeroes in N. */
179 INLINE
180 static int
181 clzusi (USItype n)
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);
194 extern FLO_type pack_d (const fp_number_type * );
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)
200 FLO_union_type dst;
201 fractype fraction = src->fraction.ll; /* wasn't unsigned before? */
202 int sign = src->sign;
203 int exp = 0;
205 if (isnan (src))
207 exp = EXPMAX;
208 /* Restore the NaN's payload. */
209 fraction >>= NGARDS;
210 fraction &= QUIET_NAN - 1;
211 if (src->class == CLASS_QNAN || 1)
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
224 else if (isinf (src))
226 exp = EXPMAX;
227 fraction = 0;
229 else if (iszero (src))
231 exp = 0;
232 fraction = 0;
234 else if (fraction == 0)
236 exp = 0;
238 else
240 if (__builtin_expect (src->normal_exp < NORMAL_EXPMIN, 0))
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. */
253 int shift = NORMAL_EXPMIN - src->normal_exp;
255 exp = 0;
257 if (shift > FRAC_NBITS - NGARDS)
259 /* No point shifting, since it's more that 64 out. */
260 fraction = 0;
262 else
264 int lowbit = (fraction & (((fractype)1 << shift) - 1)) ? 1 : 0;
265 fraction = (fraction >> shift) | lowbit;
267 if ((fraction & GARDMASK) == GARDMSB)
269 if ((fraction & (1 << NGARDS)))
270 fraction += GARDROUND + 1;
272 else
274 /* Add to the guards to round up. */
275 fraction += GARDROUND;
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)
281 exp += 1;
283 fraction >>= NGARDS;
284 #endif /* NO_DENORMALS */
286 else if (__builtin_expect (src->normal_exp > EXPBIAS, 0))
288 exp = EXPMAX;
289 fraction = 0;
291 else
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)
299 if (fraction & (1 << NGARDS))
300 fraction += GARDROUND + 1;
302 else
304 /* Add a one to the guards to round up */
305 fraction += GARDROUND;
307 if (fraction >= IMPLICIT_2)
309 fraction >>= 1;
310 exp += 1;
312 fraction >>= NGARDS;
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 #if defined TFLOAT && defined HALFFRACBITS
321 halffractype high, low, unity;
322 int lowsign, lowexp;
324 unity = (halffractype) 1 << HALFFRACBITS;
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);
331 /* Get the initial sign and exponent of the low double. */
332 lowexp = exp - HALFFRACBITS - 1;
333 lowsign = sign;
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))
340 /* Round HIGH up and adjust LOW to match. */
341 high++;
342 if (high == unity)
344 /* May make it infinite, but that's OK. */
345 high = 0;
346 exp++;
348 low = unity * 2 - low;
349 lowsign ^= 1;
352 high |= (halffractype) exp << HALFFRACBITS;
353 high |= (halffractype) sign << (HALFFRACBITS + EXPBITS);
355 if (exp == EXPMAX || exp == 0 || low == 0)
356 low = 0;
357 else
359 while (lowexp > 0 && low < unity)
361 low <<= 1;
362 lowexp--;
365 if (lowexp <= 0)
367 halffractype roundmsb, round;
368 int shift;
370 shift = 1 - lowexp;
371 roundmsb = (1 << (shift - 1));
372 round = low & ((roundmsb << 1) - 1);
374 low >>= shift;
375 lowexp = 0;
377 if (round > roundmsb || (round == roundmsb && (low & 1) == 1))
379 low++;
380 if (low == unity)
381 /* LOW rounds up to the smallest normal number. */
382 lowexp++;
386 low &= unity - 1;
387 low |= (halffractype) lowexp << HALFFRACBITS;
388 low |= (halffractype) lowsign << (HALFFRACBITS + EXPBITS);
390 dst.value_raw = ((fractype) high << HALFSHIFT) | low;
392 #else
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
398 #if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT)
399 #ifdef TFLOAT
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;
408 #else
410 halffractype tmp = dst.words[0];
411 dst.words[0] = dst.words[1];
412 dst.words[1] = tmp;
414 #endif
415 #endif
417 return dst.value;
419 #endif
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)
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;
432 #if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT)
433 FLO_union_type swapped;
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
447 #if defined TFLOAT && defined HALFFRACBITS
449 halffractype high, low;
451 high = src->value_raw >> HALFSHIFT;
452 low = src->value_raw & (((fractype)1 << HALFSHIFT) - 1);
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;
459 if (exp != EXPMAX && exp != 0 && low != 0)
461 int lowexp = ((int)(low >> HALFFRACBITS)) & ((1 << EXPBITS) - 1);
462 int lowsign = ((int)(low >> (((HALFFRACBITS + EXPBITS))))) & 1;
463 int shift;
464 fractype xlow;
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
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--;
490 #else
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
496 dst->sign = sign;
497 if (exp == 0)
499 /* Hmm. Looks like 0 */
500 if (fraction == 0
501 #ifdef NO_DENORMALS
502 || 1
503 #endif
506 /* tastes like zero */
507 dst->class = CLASS_ZERO;
509 else
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;
517 dst->class = CLASS_NUMBER;
518 #if 1
519 while (fraction < IMPLICIT_1)
521 fraction <<= 1;
522 dst->normal_exp--;
524 #endif
525 dst->fraction.ll = fraction;
528 else if (__builtin_expect (exp == EXPMAX, 0))
530 /* Huge exponent*/
531 if (fraction == 0)
533 /* Attached to a zero fraction - means infinity */
534 dst->class = CLASS_INFINITY;
536 else
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
545 dst->class = CLASS_QNAN;
547 else
549 dst->class = CLASS_SNAN;
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;
557 else
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;
565 #endif /* L_unpack_df || L_unpack_sf */
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)
573 intfrac tfraction;
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;
581 if (isnan (a))
583 return a;
585 if (isnan (b))
587 return b;
589 if (isinf (a))
591 /* Adding infinities with opposite signs yields a NaN. */
592 if (isinf (b) && a->sign != b->sign)
593 return makenan ();
594 return a;
596 if (isinf (b))
598 return b;
600 if (iszero (b))
602 if (iszero (a))
604 *tmp = *a;
605 tmp->sign = a->sign & b->sign;
606 return tmp;
608 return a;
610 if (iszero (a))
612 return b;
615 /* Got two numbers. shift the smaller and increment the exponent till
616 they're the same */
618 int diff;
619 int sdiff;
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;
626 diff = a_normal_exp - b_normal_exp;
627 sdiff = diff;
629 if (diff < 0)
630 diff = -diff;
631 if (diff < FRAC_NBITS)
633 if (sdiff > 0)
635 b_normal_exp += diff;
636 LSHIFT (b_fraction, diff);
638 else if (sdiff < 0)
640 a_normal_exp += diff;
641 LSHIFT (a_fraction, diff);
644 else
646 /* Somethings's up.. choose the biggest */
647 if (a_normal_exp > b_normal_exp)
649 b_normal_exp = a_normal_exp;
650 b_fraction = 0;
652 else
654 a_normal_exp = b_normal_exp;
655 a_fraction = 0;
660 if (a->sign != b->sign)
662 if (a->sign)
664 tfraction = -a_fraction + b_fraction;
666 else
668 tfraction = a_fraction - b_fraction;
670 if (tfraction >= 0)
672 tmp->sign = 0;
673 tmp->normal_exp = a_normal_exp;
674 tmp->fraction.ll = tfraction;
676 else
678 tmp->sign = 1;
679 tmp->normal_exp = a_normal_exp;
680 tmp->fraction.ll = -tfraction;
682 /* and renormalize it */
684 while (tmp->fraction.ll < IMPLICIT_1 && tmp->fraction.ll)
686 tmp->fraction.ll <<= 1;
687 tmp->normal_exp--;
690 else
692 tmp->sign = a->sign;
693 tmp->normal_exp = a_normal_exp;
694 tmp->fraction.ll = a_fraction + b_fraction;
696 tmp->class = CLASS_NUMBER;
697 /* Now the fraction is added, we have to shift down to renormalize the
698 number */
700 if (tmp->fraction.ll >= IMPLICIT_2)
702 LSHIFT (tmp->fraction.ll, 1);
703 tmp->normal_exp++;
705 return tmp;
708 FLO_type
709 add (FLO_type arg_a, FLO_type arg_b)
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;
717 au.value = arg_a;
718 bu.value = arg_b;
720 unpack_d (&au, &a);
721 unpack_d (&bu, &b);
723 res = _fpadd_parts (&a, &b, &tmp);
725 return pack_d (res);
728 FLO_type
729 sub (FLO_type arg_a, FLO_type arg_b)
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;
737 au.value = arg_a;
738 bu.value = arg_b;
740 unpack_d (&au, &a);
741 unpack_d (&bu, &b);
743 b.sign ^= 1;
745 res = _fpadd_parts (&a, &b, &tmp);
747 return pack_d (res);
749 #endif /* L_addsub_sf || L_addsub_df */
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)
757 fractype low = 0;
758 fractype high = 0;
760 if (isnan (a))
762 a->sign = a->sign != b->sign;
763 return a;
765 if (isnan (b))
767 b->sign = a->sign != b->sign;
768 return b;
770 if (isinf (a))
772 if (iszero (b))
773 return makenan ();
774 a->sign = a->sign != b->sign;
775 return a;
777 if (isinf (b))
779 if (iszero (a))
781 return makenan ();
783 b->sign = a->sign != b->sign;
784 return b;
786 if (iszero (a))
788 a->sign = a->sign != b->sign;
789 return a;
791 if (iszero (b))
793 b->sign = a->sign != b->sign;
794 return b;
797 /* Calculate the mantissa by multiplying both numbers to get a
798 twice-as-wide number. */
800 #if defined(NO_DI_MODE) || defined(TFLOAT)
802 fractype x = a->fraction.ll;
803 fractype ylow = b->fraction.ll;
804 fractype yhigh = 0;
805 int bit;
807 /* ??? This does multiplies one bit at a time. Optimize. */
808 for (bit = 0; bit < FRAC_NBITS; bit++)
810 int carry;
812 if (x & 1)
814 carry = (low += ylow) < ylow;
815 high += yhigh + carry;
817 yhigh <<= 1;
818 if (ylow & FRACHIGH)
820 yhigh |= 1;
822 ylow <<= 1;
823 x >>= 1;
826 #elif defined(FLOAT)
827 /* Multiplying two USIs to get a UDI, we're safe. */
829 UDItype answer = (UDItype)a->fraction.ll * (UDItype)b->fraction.ll;
831 high = answer >> BITS_PER_SI;
832 low = answer;
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. */
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;
860 #endif
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)
868 tmp->normal_exp++;
869 if (high & 1)
871 low >>= 1;
872 low |= FRACHIGH;
874 high >>= 1;
876 while (high < IMPLICIT_1)
878 tmp->normal_exp--;
880 high <<= 1;
881 if (low & FRACHIGH)
882 high |= 1;
883 low <<= 1;
886 if ((high & GARDMASK) == GARDMSB)
888 if (high & (1 << NGARDS))
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). */
897 else if (low)
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;
905 /* Avoid further rounding in pack_d. */
906 high &= ~(fractype) GARDMASK;
909 tmp->fraction.ll = high;
910 tmp->class = CLASS_NUMBER;
911 return tmp;
914 FLO_type
915 multiply (FLO_type arg_a, FLO_type arg_b)
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;
923 au.value = arg_a;
924 bu.value = arg_b;
926 unpack_d (&au, &a);
927 unpack_d (&bu, &b);
929 res = _fpmul_parts (&a, &b, &tmp);
931 return pack_d (res);
933 #endif /* L_mul_sf || L_mul_df || L_mul_tf */
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)
940 fractype bit;
941 fractype numerator;
942 fractype denominator;
943 fractype quotient;
945 if (isnan (a))
947 return a;
949 if (isnan (b))
951 return b;
954 a->sign = a->sign ^ b->sign;
956 if (isinf (a) || iszero (a))
958 if (a->class == b->class)
959 return makenan ();
960 return a;
963 if (isinf (b))
965 a->fraction.ll = 0;
966 a->normal_exp = 0;
967 return a;
969 if (iszero (b))
971 a->class = CLASS_INFINITY;
972 return a;
975 /* Calculate the mantissa by multiplying both 64bit numbers to get a
976 128 bit number */
978 /* quotient =
979 ( numerator / denominator) * 2^(numerator exponent - denominator exponent)
982 a->normal_exp = a->normal_exp - b->normal_exp;
983 numerator = a->fraction.ll;
984 denominator = b->fraction.ll;
986 if (numerator < denominator)
988 /* Fraction will be less than 1.0 */
989 numerator *= 2;
990 a->normal_exp--;
992 bit = IMPLICIT_1;
993 quotient = 0;
994 /* ??? Does divide one bit at a time. Optimize. */
995 while (bit)
997 if (numerator >= denominator)
999 quotient |= bit;
1000 numerator -= denominator;
1002 bit >>= 1;
1003 numerator *= 2;
1006 if ((quotient & GARDMASK) == GARDMSB)
1008 if (quotient & (1 << NGARDS))
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. */
1015 else if (numerator)
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;
1023 /* Avoid further rounding in pack_d. */
1024 quotient &= ~(fractype) GARDMASK;
1028 a->fraction.ll = quotient;
1029 return (a);
1033 FLO_type
1034 divide (FLO_type arg_a, FLO_type arg_b)
1036 fp_number_type a;
1037 fp_number_type b;
1038 const fp_number_type *res;
1039 FLO_union_type au, bu;
1041 au.value = arg_a;
1042 bu.value = arg_b;
1044 unpack_d (&au, &a);
1045 unpack_d (&bu, &b);
1047 res = _fpdiv_parts (&a, &b);
1049 return pack_d (res);
1051 #endif /* L_div_sf || L_div_df */
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
1062 __fpcmp_parts (fp_number_type * a, fp_number_type * b)
1064 #if 0
1065 /* either nan -> unordered. Must be checked outside of this routine. */
1066 if (isnan (a) && isnan (b))
1068 return 1; /* still unordered! */
1070 #endif
1072 if (isnan (a) || isnan (b))
1074 return 1; /* how to indicate unordered compare? */
1076 if (isinf (a) && isinf (b))
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...
1087 return b->sign - a->sign;
1089 /* but not both... */
1090 if (isinf (a))
1092 return a->sign ? -1 : 1;
1094 if (isinf (b))
1096 return b->sign ? 1 : -1;
1098 if (iszero (a) && iszero (b))
1100 return 0;
1102 if (iszero (a))
1104 return b->sign ? 1 : -1;
1106 if (iszero (b))
1108 return a->sign ? -1 : 1;
1110 /* now both are "normal". */
1111 if (a->sign != b->sign)
1113 /* opposite signs */
1114 return a->sign ? -1 : 1;
1116 /* same sign; exponents? */
1117 if (a->normal_exp > b->normal_exp)
1119 return a->sign ? -1 : 1;
1121 if (a->normal_exp < b->normal_exp)
1123 return a->sign ? 1 : -1;
1125 /* same exponents; check size. */
1126 if (a->fraction.ll > b->fraction.ll)
1128 return a->sign ? -1 : 1;
1130 if (a->fraction.ll < b->fraction.ll)
1132 return a->sign ? 1 : -1;
1134 /* after all that, they're equal. */
1135 return 0;
1137 #endif
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)
1143 fp_number_type a;
1144 fp_number_type b;
1145 FLO_union_type au, bu;
1147 au.value = arg_a;
1148 bu.value = arg_b;
1150 unpack_d (&au, &a);
1151 unpack_d (&bu, &b);
1153 return __fpcmp_parts (&a, &b);
1155 #endif /* L_compare_sf || L_compare_df */
1157 /* These should be optimized for their specific tasks someday. */
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)
1163 fp_number_type a;
1164 fp_number_type b;
1165 FLO_union_type au, bu;
1167 au.value = arg_a;
1168 bu.value = arg_b;
1170 unpack_d (&au, &a);
1171 unpack_d (&bu, &b);
1173 if (isnan (&a) || isnan (&b))
1174 return 1; /* false, truth == 0 */
1176 return __fpcmp_parts (&a, &b) ;
1178 #endif /* L_eq_sf || L_eq_df */
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)
1184 fp_number_type a;
1185 fp_number_type b;
1186 FLO_union_type au, bu;
1188 au.value = arg_a;
1189 bu.value = arg_b;
1191 unpack_d (&au, &a);
1192 unpack_d (&bu, &b);
1194 if (isnan (&a) || isnan (&b))
1195 return 1; /* true, truth != 0 */
1197 return __fpcmp_parts (&a, &b) ;
1199 #endif /* L_ne_sf || L_ne_df */
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)
1205 fp_number_type a;
1206 fp_number_type b;
1207 FLO_union_type au, bu;
1209 au.value = arg_a;
1210 bu.value = arg_b;
1212 unpack_d (&au, &a);
1213 unpack_d (&bu, &b);
1215 if (isnan (&a) || isnan (&b))
1216 return -1; /* false, truth > 0 */
1218 return __fpcmp_parts (&a, &b);
1220 #endif /* L_gt_sf || L_gt_df */
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)
1226 fp_number_type a;
1227 fp_number_type b;
1228 FLO_union_type au, bu;
1230 au.value = arg_a;
1231 bu.value = arg_b;
1233 unpack_d (&au, &a);
1234 unpack_d (&bu, &b);
1236 if (isnan (&a) || isnan (&b))
1237 return -1; /* false, truth >= 0 */
1238 return __fpcmp_parts (&a, &b) ;
1240 #endif /* L_ge_sf || L_ge_df */
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)
1246 fp_number_type a;
1247 fp_number_type b;
1248 FLO_union_type au, bu;
1250 au.value = arg_a;
1251 bu.value = arg_b;
1253 unpack_d (&au, &a);
1254 unpack_d (&bu, &b);
1256 if (isnan (&a) || isnan (&b))
1257 return 1; /* false, truth < 0 */
1259 return __fpcmp_parts (&a, &b);
1261 #endif /* L_lt_sf || L_lt_df */
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)
1267 fp_number_type a;
1268 fp_number_type b;
1269 FLO_union_type au, bu;
1271 au.value = arg_a;
1272 bu.value = arg_b;
1274 unpack_d (&au, &a);
1275 unpack_d (&bu, &b);
1277 if (isnan (&a) || isnan (&b))
1278 return 1; /* false, truth <= 0 */
1280 return __fpcmp_parts (&a, &b) ;
1282 #endif /* L_le_sf || L_le_df */
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)
1288 fp_number_type a;
1289 fp_number_type b;
1290 FLO_union_type au, bu;
1292 au.value = arg_a;
1293 bu.value = arg_b;
1295 unpack_d (&au, &a);
1296 unpack_d (&bu, &b);
1298 return (isnan (&a) || isnan (&b));
1300 #endif /* L_unord_sf || L_unord_df */
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)
1306 fp_number_type in;
1308 in.class = CLASS_NUMBER;
1309 in.sign = arg_a < 0;
1310 if (!arg_a)
1312 in.class = CLASS_ZERO;
1314 else
1316 USItype uarg;
1317 int shift;
1318 in.normal_exp = FRACBITS + NGARDS;
1319 if (in.sign)
1321 /* Special case for minint, since there is no +ve integer
1322 representation for it */
1323 if (arg_a == (- MAX_SI_INT - 1))
1325 return (FLO_type)(- MAX_SI_INT - 1);
1327 uarg = (-arg_a);
1329 else
1330 uarg = arg_a;
1332 in.fraction.ll = uarg;
1333 shift = clzusi (uarg) - (BITS_PER_SI - 1 - FRACBITS - NGARDS);
1334 if (shift > 0)
1336 in.fraction.ll <<= shift;
1337 in.normal_exp -= shift;
1340 return pack_d (&in);
1342 #endif /* L_si_to_sf || L_si_to_df */
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)
1348 fp_number_type in;
1350 in.sign = 0;
1351 if (!arg_a)
1353 in.class = CLASS_ZERO;
1355 else
1357 int shift;
1358 in.class = CLASS_NUMBER;
1359 in.normal_exp = FRACBITS + NGARDS;
1360 in.fraction.ll = arg_a;
1362 shift = clzusi (arg_a) - (BITS_PER_SI - 1 - FRACBITS - NGARDS);
1363 if (shift < 0)
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;
1370 else if (shift > 0)
1372 in.fraction.ll <<= shift;
1373 in.normal_exp -= shift;
1376 return pack_d (&in);
1378 #endif
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)
1384 fp_number_type a;
1385 SItype tmp;
1386 FLO_union_type au;
1388 au.value = arg_a;
1389 unpack_d (&au, &a);
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);
1406 #endif /* L_sf_to_si || L_df_to_si */
1408 #if defined(L_tf_to_usi)
1409 USItype
1410 float_to_usi (FLO_type arg_a)
1412 fp_number_type a;
1413 FLO_union_type au;
1415 au.value = arg_a;
1416 unpack_d (&au, &a);
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);
1438 #endif /* L_tf_to_usi */
1440 #if defined(L_negate_sf) || defined(L_negate_df) || defined(L_negate_tf)
1441 FLO_type
1442 negate (FLO_type arg_a)
1444 fp_number_type a;
1445 FLO_union_type au;
1447 au.value = arg_a;
1448 unpack_d (&au, &a);
1450 flip_sign (&a);
1451 return pack_d (&a);
1453 #endif /* L_negate_sf || L_negate_df */
1455 #ifdef FLOAT
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)
1464 fp_number_type in;
1466 in.class = class;
1467 in.sign = sign;
1468 in.normal_exp = exp;
1469 in.fraction.ll = frac;
1470 return pack_d (&in);
1472 #endif /* L_make_sf */
1474 #ifndef FLOAT_ONLY
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. */
1481 #if defined(L_sf_to_df)
1482 DFtype
1483 sf_to_df (SFtype arg_a)
1485 fp_number_type in;
1486 FLO_union_type au;
1488 au.value = arg_a;
1489 unpack_d (&au, &in);
1491 return __make_dp (in.class, in.sign, in.normal_exp,
1492 ((UDItype) in.fraction.ll) << F_D_BITOFF);
1494 #endif /* L_sf_to_df */
1496 #if defined(L_sf_to_tf) && defined(TMODES)
1497 TFtype
1498 sf_to_tf (SFtype arg_a)
1500 fp_number_type in;
1501 FLO_union_type au;
1503 au.value = arg_a;
1504 unpack_d (&au, &in);
1506 return __make_tp (in.class, in.sign, in.normal_exp,
1507 ((UTItype) in.fraction.ll) << F_T_BITOFF);
1509 #endif /* L_sf_to_df */
1511 #endif /* ! FLOAT_ONLY */
1512 #endif /* FLOAT */
1514 #ifndef FLOAT
1516 extern SFtype __make_fp (fp_class_type, unsigned int, int, USItype);
1518 #if defined(L_make_df)
1519 DFtype
1520 __make_dp (fp_class_type class, unsigned int sign, int exp, UDItype frac)
1522 fp_number_type in;
1524 in.class = class;
1525 in.sign = sign;
1526 in.normal_exp = exp;
1527 in.fraction.ll = frac;
1528 return pack_d (&in);
1530 #endif /* L_make_df */
1532 #if defined(L_df_to_sf)
1533 SFtype
1534 df_to_sf (DFtype arg_a)
1536 fp_number_type in;
1537 USItype sffrac;
1538 FLO_union_type au;
1540 au.value = arg_a;
1541 unpack_d (&au, &in);
1543 sffrac = in.fraction.ll >> F_D_BITOFF;
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;
1550 return __make_fp (in.class, in.sign, in.normal_exp, sffrac);
1552 #endif /* L_df_to_sf */
1554 #if defined(L_df_to_tf) && defined(TMODES) \
1555 && !defined(FLOAT) && !defined(TFLOAT)
1556 TFtype
1557 df_to_tf (DFtype arg_a)
1559 fp_number_type in;
1560 FLO_union_type au;
1562 au.value = arg_a;
1563 unpack_d (&au, &in);
1565 return __make_tp (in.class, in.sign, in.normal_exp,
1566 ((UTItype) in.fraction.ll) << D_T_BITOFF);
1568 #endif /* L_sf_to_df */
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)
1578 fp_number_type in;
1580 in.class = class;
1581 in.sign = sign;
1582 in.normal_exp = exp;
1583 in.fraction.ll = frac;
1584 return pack_d (&in);
1586 #endif /* L_make_tf */
1588 #if defined(L_tf_to_df)
1589 DFtype
1590 tf_to_df (TFtype arg_a)
1592 fp_number_type in;
1593 UDItype sffrac;
1594 FLO_union_type au;
1596 au.value = arg_a;
1597 unpack_d (&au, &in);
1599 sffrac = in.fraction.ll >> D_T_BITOFF;
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;
1606 return __make_dp (in.class, in.sign, in.normal_exp, sffrac);
1608 #endif /* L_tf_to_df */
1610 #if defined(L_tf_to_sf)
1611 SFtype
1612 tf_to_sf (TFtype arg_a)
1614 fp_number_type in;
1615 USItype sffrac;
1616 FLO_union_type au;
1618 au.value = arg_a;
1619 unpack_d (&au, &in);
1621 sffrac = in.fraction.ll >> F_T_BITOFF;
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;
1628 return __make_fp (in.class, in.sign, in.normal_exp, sffrac);
1630 #endif /* L_tf_to_sf */
1631 #endif /* TFLOAT */
1633 #endif /* ! FLOAT */
1634 #endif /* !EXTENDED_FLOAT_STUBS */