1 /* This is a software floating point library which can be used instead of
2 the floating point routines in libgcc1.c for targets without hardware
3 floating point. */
4
5 /* Copyright (C) 1994-2024 Free Software Foundation, Inc.
6
7 This program is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 3 of the License, or
10 (at your option) any later version.
11
12 This program is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program. If not, see <http://www.gnu.org/licenses/>. */
19
20 /* As a special exception, if you link this library with other files,
21 some of which are compiled with GCC, to produce an executable,
22 this library does not by itself cause the resulting executable
23 to be covered by the GNU General Public License.
24 This exception does not however invalidate any other reasons why
25 the executable file might be covered by the GNU General Public License. */
26
27 /* This implements IEEE 754 format arithmetic, but does not provide a
28 mechanism for setting the rounding mode, or for generating or handling
29 exceptions.
30
31 The original code by Steve Chamberlain, hacked by Mark Eichin and Jim
32 Wilson, all of Cygnus Support. */
33
34 /* The intended way to use this file is to make two copies, add `#define FLOAT'
35 to one copy, then compile both copies and add them to libgcc.a. */
36
37 /* The following macros can be defined to change the behaviour of this file:
38 FLOAT: Implement a `float', aka SFmode, fp library. If this is not
39 defined, then this file implements a `double', aka DFmode, fp library.
40 FLOAT_ONLY: Used with FLOAT, to implement a `float' only library, i.e.
41 don't include float->double conversion which requires the double library.
42 This is useful only for machines which can't support doubles, e.g. some
43 8-bit processors.
44 CMPtype: Specify the type that floating point compares should return.
45 This defaults to SItype, aka int.
46 US_SOFTWARE_GOFAST: This makes all entry points use the same names as the
47 US Software goFast library. If this is not defined, the entry points use
48 the same names as libgcc1.c.
49 _DEBUG_BITFLOAT: This makes debugging the code a little easier, by adding
50 two integers to the FLO_union_type.
51 NO_NANS: Disable nan and infinity handling
52 SMALL_MACHINE: Useful when operations on QIs and HIs are faster
53 than on an SI */
54
55 #ifndef SFtype
56 typedef SFtype __attribute__ ((mode (SF)));
57 #endif
58 #ifndef DFtype
59 typedef DFtype __attribute__ ((mode (DF)));
60 #endif
61
62 #ifndef HItype
63 typedef int HItype __attribute__ ((mode (HI)));
64 #endif
65 #ifndef SItype
66 typedef int SItype __attribute__ ((mode (SI)));
67 #endif
68 #ifndef DItype
69 typedef int DItype __attribute__ ((mode (DI)));
70 #endif
71
72 /* The type of the result of a fp compare */
73 #ifndef CMPtype
74 #define CMPtype SItype
75 #endif
76
77 #ifndef UHItype
78 typedef unsigned int UHItype __attribute__ ((mode (HI)));
79 #endif
80 #ifndef USItype
81 typedef unsigned int USItype __attribute__ ((mode (SI)));
82 #endif
83 #ifndef UDItype
84 typedef unsigned int UDItype __attribute__ ((mode (DI)));
85 #endif
86
87 #define MAX_SI_INT ((SItype) ((unsigned) (~0)>>1))
88 #define MAX_USI_INT ((USItype) ~0)
89
90
91 #ifdef FLOAT_ONLY
92 #define NO_DI_MODE
93 #endif
94
95 #ifdef FLOAT
96 # define NGARDS 7L
97 # define GARDROUND 0x3f
98 # define GARDMASK 0x7f
99 # define GARDMSB 0x40
100 # define EXPBITS 8
101 # define EXPBIAS 127
102 # define FRACBITS 23
103 # define EXPMAX (0xff)
104 # define QUIET_NAN 0x100000L
105 # define FRAC_NBITS 32
106 # define FRACHIGH 0x80000000L
107 # define FRACHIGH2 0xc0000000L
108 typedef USItype fractype;
109 typedef UHItype halffractype;
110 typedef SFtype FLO_type;
111 typedef SItype intfrac;
112
113 #else
114 # define PREFIXFPDP dp
115 # define PREFIXSFDF df
116 # define NGARDS 8L
117 # define GARDROUND 0x7f
118 # define GARDMASK 0xff
119 # define GARDMSB 0x80
120 # define EXPBITS 11
121 # define EXPBIAS 1023
122 # define FRACBITS 52
123 # define EXPMAX (0x7ff)
124 # define QUIET_NAN 0x8000000000000LL
125 # define FRAC_NBITS 64
126 # define FRACHIGH 0x8000000000000000LL
127 # define FRACHIGH2 0xc000000000000000LL
128 typedef UDItype fractype;
129 typedef USItype halffractype;
130 typedef DFtype FLO_type;
131 typedef DItype intfrac;
132 #endif
133
134 #ifdef US_SOFTWARE_GOFAST
135 # ifdef FLOAT
136 # define add fpadd
137 # define sub fpsub
138 # define multiply fpmul
139 # define divide fpdiv
140 # define compare fpcmp
141 # define si_to_float sitofp
142 # define float_to_si fptosi
143 # define float_to_usi fptoui
144 # define negate __negsf2
145 # define sf_to_df fptodp
146 # define dptofp dptofp
147 #else
148 # define add dpadd
149 # define sub dpsub
150 # define multiply dpmul
151 # define divide dpdiv
152 # define compare dpcmp
153 # define si_to_float litodp
154 # define float_to_si dptoli
155 # define float_to_usi dptoul
156 # define negate __negdf2
157 # define df_to_sf dptofp
158 #endif
159 #else
160 # ifdef FLOAT
161 # define add __addsf3
162 # define sub __subsf3
163 # define multiply __mulsf3
164 # define divide __divsf3
165 # define compare __cmpsf2
166 # define _eq_f2 __eqsf2
167 # define _ne_f2 __nesf2
168 # define _gt_f2 __gtsf2
169 # define _ge_f2 __gesf2
170 # define _lt_f2 __ltsf2
171 # define _le_f2 __lesf2
172 # define si_to_float __floatsisf
173 # define float_to_si __fixsfsi
174 # define float_to_usi __fixunssfsi
175 # define negate __negsf2
176 # define sf_to_df __extendsfdf2
177 #else
178 # define add __adddf3
179 # define sub __subdf3
180 # define multiply __muldf3
181 # define divide __divdf3
182 # define compare __cmpdf2
183 # define _eq_f2 __eqdf2
184 # define _ne_f2 __nedf2
185 # define _gt_f2 __gtdf2
186 # define _ge_f2 __gedf2
187 # define _lt_f2 __ltdf2
188 # define _le_f2 __ledf2
189 # define si_to_float __floatsidf
190 # define float_to_si __fixdfsi
191 # define float_to_usi __fixunsdfsi
192 # define negate __negdf2
193 # define df_to_sf __truncdfsf2
194 # endif
195 #endif
196
197
198 #ifndef INLINE
199 #define INLINE __inline__
200 #endif
201
202 /* Preserve the sticky-bit when shifting fractions to the right. */
203 #define LSHIFT(a) { a = (a & 1) | (a >> 1); }
204
205 /* numeric parameters */
206 /* F_D_BITOFF is the number of bits offset between the MSB of the mantissa
207 of a float and of a double. Assumes there are only two float types.
208 (double::FRAC_BITS+double::NGARGS-(float::FRAC_BITS-float::NGARDS))
209 */
210 #define F_D_BITOFF (52+8-(23+7))
211
212
213 #define NORMAL_EXPMIN (-(EXPBIAS)+1)
214 #define IMPLICIT_1 (1LL<<(FRACBITS+NGARDS))
215 #define IMPLICIT_2 (1LL<<(FRACBITS+1+NGARDS))
216
217 /* common types */
218
219 typedef enum
220 {
221 CLASS_SNAN,
222 CLASS_QNAN,
223 CLASS_ZERO,
224 CLASS_NUMBER,
225 CLASS_INFINITY
226 } fp_class_type;
227
228 typedef struct
229 {
230 #ifdef SMALL_MACHINE
231 char class;
232 unsigned char sign;
233 short normal_exp;
234 #else
235 fp_class_type class;
236 unsigned int sign;
237 int normal_exp;
238 #endif
239
240 union
241 {
242 fractype ll;
243 halffractype l[2];
244 } fraction;
245 } fp_number_type;
246
247 typedef union
248 {
249 FLO_type value;
250 #ifdef _DEBUG_BITFLOAT
251 int l[2];
252 #endif
253 struct
254 {
255 #ifndef FLOAT_BIT_ORDER_MISMATCH
256 unsigned int sign:1 ATTRIBUTE_PACKED;
257 unsigned int exp:EXPBITS ATTRIBUTE_PACKED;
258 fractype fraction:FRACBITS ATTRIBUTE_PACKED;
259 #else
260 fractype fraction:FRACBITS ATTRIBUTE_PACKED;
261 unsigned int exp:EXPBITS ATTRIBUTE_PACKED;
262 unsigned int sign:1 ATTRIBUTE_PACKED;
263 #endif
264 }
265 bits;
266 }
267 FLO_union_type;
268
269
270 /* end of header */
271
272 /* IEEE "special" number predicates */
273
274 #ifdef NO_NANS
275
276 #define nan() 0
277 #define isnan(x) 0
278 #define isinf(x) 0
279 #else
280
281 INLINE
282 static fp_number_type *
nan()283 nan ()
284 {
285 static fp_number_type thenan;
286
287 return &thenan;
288 }
289
290 INLINE
291 static int
isnan(fp_number_type * x)292 isnan ( fp_number_type * x)
293 {
294 return x->class == CLASS_SNAN || x->class == CLASS_QNAN;
295 }
296
297 INLINE
298 static int
isinf(fp_number_type * x)299 isinf ( fp_number_type * x)
300 {
301 return x->class == CLASS_INFINITY;
302 }
303
304 #endif
305
306 INLINE
307 static int
iszero(fp_number_type * x)308 iszero ( fp_number_type * x)
309 {
310 return x->class == CLASS_ZERO;
311 }
312
313 INLINE
314 static void
flip_sign(fp_number_type * x)315 flip_sign ( fp_number_type * x)
316 {
317 x->sign = !x->sign;
318 }
319
320 static FLO_type
pack_d(fp_number_type * src)321 pack_d ( fp_number_type * src)
322 {
323 FLO_union_type dst;
324 fractype fraction = src->fraction.ll; /* wasn't unsigned before? */
325
326 dst.bits.sign = src->sign;
327
328 if (isnan (src))
329 {
330 dst.bits.exp = EXPMAX;
331 dst.bits.fraction = src->fraction.ll;
332 if (src->class == CLASS_QNAN || 1)
333 {
334 dst.bits.fraction |= QUIET_NAN;
335 }
336 }
337 else if (isinf (src))
338 {
339 dst.bits.exp = EXPMAX;
340 dst.bits.fraction = 0;
341 }
342 else if (iszero (src))
343 {
344 dst.bits.exp = 0;
345 dst.bits.fraction = 0;
346 }
347 else if (fraction == 0)
348 {
349 dst.value = 0;
350 }
351 else
352 {
353 if (src->normal_exp < NORMAL_EXPMIN)
354 {
355 /* This number's exponent is too low to fit into the bits
356 available in the number, so we'll store 0 in the exponent and
357 shift the fraction to the right to make up for it. */
358
359 int shift = NORMAL_EXPMIN - src->normal_exp;
360
361 dst.bits.exp = 0;
362
363 if (shift > FRAC_NBITS - NGARDS)
364 {
365 /* No point shifting, since it's more that 64 out. */
366 fraction = 0;
367 }
368 else
369 {
370 /* Shift by the value */
371 fraction >>= shift;
372 }
373 fraction >>= NGARDS;
374 dst.bits.fraction = fraction;
375 }
376 else if (src->normal_exp > EXPBIAS)
377 {
378 dst.bits.exp = EXPMAX;
379 dst.bits.fraction = 0;
380 }
381 else
382 {
383 dst.bits.exp = src->normal_exp + EXPBIAS;
384 /* IF the gard bits are the all zero, but the first, then we're
385 half way between two numbers, choose the one which makes the
386 lsb of the answer 0. */
387 if ((fraction & GARDMASK) == GARDMSB)
388 {
389 if (fraction & (1 << NGARDS))
390 fraction += GARDROUND + 1;
391 }
392 else
393 {
394 /* Add a one to the guards to round up */
395 fraction += GARDROUND;
396 }
397 if (fraction >= IMPLICIT_2)
398 {
399 fraction >>= 1;
400 dst.bits.exp += 1;
401 }
402 fraction >>= NGARDS;
403 dst.bits.fraction = fraction;
404 }
405 }
406 return dst.value;
407 }
408
409 static void
unpack_d(FLO_union_type * src,fp_number_type * dst)410 unpack_d (FLO_union_type * src, fp_number_type * dst)
411 {
412 fractype fraction = src->bits.fraction;
413
414 dst->sign = src->bits.sign;
415 if (src->bits.exp == 0)
416 {
417 /* Hmm. Looks like 0 */
418 if (fraction == 0)
419 {
420 /* tastes like zero */
421 dst->class = CLASS_ZERO;
422 }
423 else
424 {
425 /* Zero exponent with non zero fraction - it's denormalized,
426 so there isn't a leading implicit one - we'll shift it so
427 it gets one. */
428 dst->normal_exp = src->bits.exp - EXPBIAS + 1;
429 fraction <<= NGARDS;
430
431 dst->class = CLASS_NUMBER;
432 #if 1
433 while (fraction < IMPLICIT_1)
434 {
435 fraction <<= 1;
436 dst->normal_exp--;
437 }
438 #endif
439 dst->fraction.ll = fraction;
440 }
441 }
442 else if (src->bits.exp == EXPMAX)
443 {
444 /* Huge exponent*/
445 if (fraction == 0)
446 {
447 /* Attached to a zero fraction - means infinity */
448 dst->class = CLASS_INFINITY;
449 }
450 else
451 {
452 /* Non zero fraction, means nan */
453 if (dst->sign)
454 {
455 dst->class = CLASS_SNAN;
456 }
457 else
458 {
459 dst->class = CLASS_QNAN;
460 }
461 /* Keep the fraction part as the nan number */
462 dst->fraction.ll = fraction;
463 }
464 }
465 else
466 {
467 /* Nothing strange about this number */
468 dst->normal_exp = src->bits.exp - EXPBIAS;
469 dst->class = CLASS_NUMBER;
470 dst->fraction.ll = (fraction << NGARDS) | IMPLICIT_1;
471 }
472 }
473
474 static fp_number_type *
_fpadd_parts(fp_number_type * a,fp_number_type * b,fp_number_type * tmp)475 _fpadd_parts (fp_number_type * a,
476 fp_number_type * b,
477 fp_number_type * tmp)
478 {
479 intfrac tfraction;
480
481 /* Put commonly used fields in local variables. */
482 int a_normal_exp;
483 int b_normal_exp;
484 fractype a_fraction;
485 fractype b_fraction;
486
487 if (isnan (a))
488 {
489 return a;
490 }
491 if (isnan (b))
492 {
493 return b;
494 }
495 if (isinf (a))
496 {
497 /* Adding infinities with opposite signs yields a NaN. */
498 if (isinf (b) && a->sign != b->sign)
499 return nan ();
500 return a;
501 }
502 if (isinf (b))
503 {
504 return b;
505 }
506 if (iszero (b))
507 {
508 return a;
509 }
510 if (iszero (a))
511 {
512 return b;
513 }
514
515 /* Got two numbers. shift the smaller and increment the exponent till
516 they're the same */
517 {
518 int diff;
519
520 a_normal_exp = a->normal_exp;
521 b_normal_exp = b->normal_exp;
522 a_fraction = a->fraction.ll;
523 b_fraction = b->fraction.ll;
524
525 diff = a_normal_exp - b_normal_exp;
526
527 if (diff < 0)
528 diff = -diff;
529 if (diff < FRAC_NBITS)
530 {
531 /* ??? This does shifts one bit at a time. Optimize. */
532 while (a_normal_exp > b_normal_exp)
533 {
534 b_normal_exp++;
535 LSHIFT (b_fraction);
536 }
537 while (b_normal_exp > a_normal_exp)
538 {
539 a_normal_exp++;
540 LSHIFT (a_fraction);
541 }
542 }
543 else
544 {
545 /* Somethings's up.. choose the biggest */
546 if (a_normal_exp > b_normal_exp)
547 {
548 b_normal_exp = a_normal_exp;
549 b_fraction = 0;
550 }
551 else
552 {
553 a_normal_exp = b_normal_exp;
554 a_fraction = 0;
555 }
556 }
557 }
558
559 if (a->sign != b->sign)
560 {
561 if (a->sign)
562 {
563 tfraction = -a_fraction + b_fraction;
564 }
565 else
566 {
567 tfraction = a_fraction - b_fraction;
568 }
569 if (tfraction > 0)
570 {
571 tmp->sign = 0;
572 tmp->normal_exp = a_normal_exp;
573 tmp->fraction.ll = tfraction;
574 }
575 else
576 {
577 tmp->sign = 1;
578 tmp->normal_exp = a_normal_exp;
579 tmp->fraction.ll = -tfraction;
580 }
581 /* and renormalize it */
582
583 while (tmp->fraction.ll < IMPLICIT_1 && tmp->fraction.ll)
584 {
585 tmp->fraction.ll <<= 1;
586 tmp->normal_exp--;
587 }
588 }
589 else
590 {
591 tmp->sign = a->sign;
592 tmp->normal_exp = a_normal_exp;
593 tmp->fraction.ll = a_fraction + b_fraction;
594 }
595 tmp->class = CLASS_NUMBER;
596 /* Now the fraction is added, we have to shift down to renormalize the
597 number */
598
599 if (tmp->fraction.ll >= IMPLICIT_2)
600 {
601 LSHIFT (tmp->fraction.ll);
602 tmp->normal_exp++;
603 }
604 return tmp;
605
606 }
607
608 FLO_type
add(FLO_type arg_a,FLO_type arg_b)609 add (FLO_type arg_a, FLO_type arg_b)
610 {
611 fp_number_type a;
612 fp_number_type b;
613 fp_number_type tmp;
614 fp_number_type *res;
615
616 unpack_d ((FLO_union_type *) & arg_a, &a);
617 unpack_d ((FLO_union_type *) & arg_b, &b);
618
619 res = _fpadd_parts (&a, &b, &tmp);
620
621 return pack_d (res);
622 }
623
624 FLO_type
sub(FLO_type arg_a,FLO_type arg_b)625 sub (FLO_type arg_a, FLO_type arg_b)
626 {
627 fp_number_type a;
628 fp_number_type b;
629 fp_number_type tmp;
630 fp_number_type *res;
631
632 unpack_d ((FLO_union_type *) & arg_a, &a);
633 unpack_d ((FLO_union_type *) & arg_b, &b);
634
635 b.sign ^= 1;
636
637 res = _fpadd_parts (&a, &b, &tmp);
638
639 return pack_d (res);
640 }
641
642 static fp_number_type *
_fpmul_parts(fp_number_type * a,fp_number_type * b,fp_number_type * tmp)643 _fpmul_parts ( fp_number_type * a,
644 fp_number_type * b,
645 fp_number_type * tmp)
646 {
647 fractype low = 0;
648 fractype high = 0;
649
650 if (isnan (a))
651 {
652 a->sign = a->sign != b->sign;
653 return a;
654 }
655 if (isnan (b))
656 {
657 b->sign = a->sign != b->sign;
658 return b;
659 }
660 if (isinf (a))
661 {
662 if (iszero (b))
663 return nan ();
664 a->sign = a->sign != b->sign;
665 return a;
666 }
667 if (isinf (b))
668 {
669 if (iszero (a))
670 {
671 return nan ();
672 }
673 b->sign = a->sign != b->sign;
674 return b;
675 }
676 if (iszero (a))
677 {
678 a->sign = a->sign != b->sign;
679 return a;
680 }
681 if (iszero (b))
682 {
683 b->sign = a->sign != b->sign;
684 return b;
685 }
686
687 /* Calculate the mantissa by multiplying both 64bit numbers to get a
688 128 bit number */
689 {
690 fractype x = a->fraction.ll;
691 fractype ylow = b->fraction.ll;
692 fractype yhigh = 0;
693 int bit;
694
695 #if defined(NO_DI_MODE)
696 {
697 /* ??? This does multiplies one bit at a time. Optimize. */
698 for (bit = 0; bit < FRAC_NBITS; bit++)
699 {
700 int carry;
701
702 if (x & 1)
703 {
704 carry = (low += ylow) < ylow;
705 high += yhigh + carry;
706 }
707 yhigh <<= 1;
708 if (ylow & FRACHIGH)
709 {
710 yhigh |= 1;
711 }
712 ylow <<= 1;
713 x >>= 1;
714 }
715 }
716 #elif defined(FLOAT)
717 {
718 /* Multiplying two 32 bit numbers to get a 64 bit number on
719 a machine with DI, so we're safe */
720
721 DItype answer = (DItype)(a->fraction.ll) * (DItype)(b->fraction.ll);
722
723 high = answer >> 32;
724 low = answer;
725 }
726 #else
727 /* Doing a 64*64 to 128 */
728 {
729 UDItype nl = a->fraction.ll & 0xffffffff;
730 UDItype nh = a->fraction.ll >> 32;
731 UDItype ml = b->fraction.ll & 0xffffffff;
732 UDItype mh = b->fraction.ll >>32;
733 UDItype pp_ll = ml * nl;
734 UDItype pp_hl = mh * nl;
735 UDItype pp_lh = ml * nh;
736 UDItype pp_hh = mh * nh;
737 UDItype res2 = 0;
738 UDItype res0 = 0;
739 UDItype ps_hh__ = pp_hl + pp_lh;
740 if (ps_hh__ < pp_hl)
741 res2 += 0x100000000LL;
742 pp_hl = (ps_hh__ << 32) & 0xffffffff00000000LL;
743 res0 = pp_ll + pp_hl;
744 if (res0 < pp_ll)
745 res2++;
746 res2 += ((ps_hh__ >> 32) & 0xffffffffL) + pp_hh;
747 high = res2;
748 low = res0;
749 }
750 #endif
751 }
752
753 tmp->normal_exp = a->normal_exp + b->normal_exp;
754 tmp->sign = a->sign != b->sign;
755 #ifdef FLOAT
756 tmp->normal_exp += 2; /* ??????????????? */
757 #else
758 tmp->normal_exp += 4; /* ??????????????? */
759 #endif
760 while (high >= IMPLICIT_2)
761 {
762 tmp->normal_exp++;
763 if (high & 1)
764 {
765 low >>= 1;
766 low |= FRACHIGH;
767 }
768 high >>= 1;
769 }
770 while (high < IMPLICIT_1)
771 {
772 tmp->normal_exp--;
773
774 high <<= 1;
775 if (low & FRACHIGH)
776 high |= 1;
777 low <<= 1;
778 }
779 /* rounding is tricky. if we only round if it won't make us round later. */
780 #if 0
781 if (low & FRACHIGH2)
782 {
783 if (((high & GARDMASK) != GARDMSB)
784 && (((high + 1) & GARDMASK) == GARDMSB))
785 {
786 /* don't round, it gets done again later. */
787 }
788 else
789 {
790 high++;
791 }
792 }
793 #endif
794 if ((high & GARDMASK) == GARDMSB)
795 {
796 if (high & (1 << NGARDS))
797 {
798 /* half way, so round to even */
799 high += GARDROUND + 1;
800 }
801 else if (low)
802 {
803 /* but we really weren't half way */
804 high += GARDROUND + 1;
805 }
806 }
807 tmp->fraction.ll = high;
808 tmp->class = CLASS_NUMBER;
809 return tmp;
810 }
811
812 FLO_type
multiply(FLO_type arg_a,FLO_type arg_b)813 multiply (FLO_type arg_a, FLO_type arg_b)
814 {
815 fp_number_type a;
816 fp_number_type b;
817 fp_number_type tmp;
818 fp_number_type *res;
819
820 unpack_d ((FLO_union_type *) & arg_a, &a);
821 unpack_d ((FLO_union_type *) & arg_b, &b);
822
823 res = _fpmul_parts (&a, &b, &tmp);
824
825 return pack_d (res);
826 }
827
828 static fp_number_type *
_fpdiv_parts(fp_number_type * a,fp_number_type * b,fp_number_type * tmp)829 _fpdiv_parts (fp_number_type * a,
830 fp_number_type * b,
831 fp_number_type * tmp)
832 {
833 fractype low = 0;
834 fractype high = 0;
835 fractype r0, r1, y0, y1, bit;
836 fractype q;
837 fractype numerator;
838 fractype denominator;
839 fractype quotient;
840 fractype remainder;
841
842 if (isnan (a))
843 {
844 return a;
845 }
846 if (isnan (b))
847 {
848 return b;
849 }
850 if (isinf (a) || iszero (a))
851 {
852 if (a->class == b->class)
853 return nan ();
854 return a;
855 }
856 a->sign = a->sign ^ b->sign;
857
858 if (isinf (b))
859 {
860 a->fraction.ll = 0;
861 a->normal_exp = 0;
862 return a;
863 }
864 if (iszero (b))
865 {
866 a->class = CLASS_INFINITY;
867 return b;
868 }
869
870 /* Calculate the mantissa by multiplying both 64bit numbers to get a
871 128 bit number */
872 {
873 int carry;
874 intfrac d0, d1; /* weren't unsigned before ??? */
875
876 /* quotient =
877 ( numerator / denominator) * 2^(numerator exponent - denominator exponent)
878 */
879
880 a->normal_exp = a->normal_exp - b->normal_exp;
881 numerator = a->fraction.ll;
882 denominator = b->fraction.ll;
883
884 if (numerator < denominator)
885 {
886 /* Fraction will be less than 1.0 */
887 numerator *= 2;
888 a->normal_exp--;
889 }
890 bit = IMPLICIT_1;
891 quotient = 0;
892 /* ??? Does divide one bit at a time. Optimize. */
893 while (bit)
894 {
895 if (numerator >= denominator)
896 {
897 quotient |= bit;
898 numerator -= denominator;
899 }
900 bit >>= 1;
901 numerator *= 2;
902 }
903
904 if ((quotient & GARDMASK) == GARDMSB)
905 {
906 if (quotient & (1 << NGARDS))
907 {
908 /* half way, so round to even */
909 quotient += GARDROUND + 1;
910 }
911 else if (numerator)
912 {
913 /* but we really weren't half way, more bits exist */
914 quotient += GARDROUND + 1;
915 }
916 }
917
918 a->fraction.ll = quotient;
919 return (a);
920 }
921 }
922
923 FLO_type
divide(FLO_type arg_a,FLO_type arg_b)924 divide (FLO_type arg_a, FLO_type arg_b)
925 {
926 fp_number_type a;
927 fp_number_type b;
928 fp_number_type tmp;
929 fp_number_type *res;
930
931 unpack_d ((FLO_union_type *) & arg_a, &a);
932 unpack_d ((FLO_union_type *) & arg_b, &b);
933
934 res = _fpdiv_parts (&a, &b, &tmp);
935
936 return pack_d (res);
937 }
938
939 /* according to the demo, fpcmp returns a comparison with 0... thus
940 a<b -> -1
941 a==b -> 0
942 a>b -> +1
943 */
944
945 static int
_fpcmp_parts(fp_number_type * a,fp_number_type * b)946 _fpcmp_parts (fp_number_type * a, fp_number_type * b)
947 {
948 #if 0
949 /* either nan -> unordered. Must be checked outside of this routine. */
950 if (isnan (a) && isnan (b))
951 {
952 return 1; /* still unordered! */
953 }
954 #endif
955
956 if (isnan (a) || isnan (b))
957 {
958 return 1; /* how to indicate unordered compare? */
959 }
960 if (isinf (a) && isinf (b))
961 {
962 /* +inf > -inf, but +inf != +inf */
963 /* b \a| +inf(0)| -inf(1)
964 ______\+--------+--------
965 +inf(0)| a==b(0)| a<b(-1)
966 -------+--------+--------
967 -inf(1)| a>b(1) | a==b(0)
968 -------+--------+--------
969 So since unordered must be non zero, just line up the columns...
970 */
971 return b->sign - a->sign;
972 }
973 /* but not both... */
974 if (isinf (a))
975 {
976 return a->sign ? -1 : 1;
977 }
978 if (isinf (b))
979 {
980 return b->sign ? 1 : -1;
981 }
982 if (iszero (a) && iszero (b))
983 {
984 return 0;
985 }
986 if (iszero (a))
987 {
988 return b->sign ? 1 : -1;
989 }
990 if (iszero (b))
991 {
992 return a->sign ? -1 : 1;
993 }
994 /* now both are "normal". */
995 if (a->sign != b->sign)
996 {
997 /* opposite signs */
998 return a->sign ? -1 : 1;
999 }
1000 /* same sign; exponents? */
1001 if (a->normal_exp > b->normal_exp)
1002 {
1003 return a->sign ? -1 : 1;
1004 }
1005 if (a->normal_exp < b->normal_exp)
1006 {
1007 return a->sign ? 1 : -1;
1008 }
1009 /* same exponents; check size. */
1010 if (a->fraction.ll > b->fraction.ll)
1011 {
1012 return a->sign ? -1 : 1;
1013 }
1014 if (a->fraction.ll < b->fraction.ll)
1015 {
1016 return a->sign ? 1 : -1;
1017 }
1018 /* after all that, they're equal. */
1019 return 0;
1020 }
1021
1022 CMPtype
compare(FLO_type arg_a,FLO_type arg_b)1023 compare (FLO_type arg_a, FLO_type arg_b)
1024 {
1025 fp_number_type a;
1026 fp_number_type b;
1027
1028 unpack_d ((FLO_union_type *) & arg_a, &a);
1029 unpack_d ((FLO_union_type *) & arg_b, &b);
1030
1031 return _fpcmp_parts (&a, &b);
1032 }
1033
1034 #ifndef US_SOFTWARE_GOFAST
1035
1036 /* These should be optimized for their specific tasks someday. */
1037
1038 CMPtype
_eq_f2(FLO_type arg_a,FLO_type arg_b)1039 _eq_f2 (FLO_type arg_a, FLO_type arg_b)
1040 {
1041 fp_number_type a;
1042 fp_number_type b;
1043
1044 unpack_d ((FLO_union_type *) & arg_a, &a);
1045 unpack_d ((FLO_union_type *) & arg_b, &b);
1046
1047 if (isnan (&a) || isnan (&b))
1048 return 1; /* false, truth == 0 */
1049
1050 return _fpcmp_parts (&a, &b) ;
1051 }
1052
1053 CMPtype
_ne_f2(FLO_type arg_a,FLO_type arg_b)1054 _ne_f2 (FLO_type arg_a, FLO_type arg_b)
1055 {
1056 fp_number_type a;
1057 fp_number_type b;
1058
1059 unpack_d ((FLO_union_type *) & arg_a, &a);
1060 unpack_d ((FLO_union_type *) & arg_b, &b);
1061
1062 if (isnan (&a) || isnan (&b))
1063 return 1; /* true, truth != 0 */
1064
1065 return _fpcmp_parts (&a, &b) ;
1066 }
1067
1068 CMPtype
_gt_f2(FLO_type arg_a,FLO_type arg_b)1069 _gt_f2 (FLO_type arg_a, FLO_type arg_b)
1070 {
1071 fp_number_type a;
1072 fp_number_type b;
1073
1074 unpack_d ((FLO_union_type *) & arg_a, &a);
1075 unpack_d ((FLO_union_type *) & arg_b, &b);
1076
1077 if (isnan (&a) || isnan (&b))
1078 return -1; /* false, truth > 0 */
1079
1080 return _fpcmp_parts (&a, &b);
1081 }
1082
1083 CMPtype
_ge_f2(FLO_type arg_a,FLO_type arg_b)1084 _ge_f2 (FLO_type arg_a, FLO_type arg_b)
1085 {
1086 fp_number_type a;
1087 fp_number_type b;
1088
1089 unpack_d ((FLO_union_type *) & arg_a, &a);
1090 unpack_d ((FLO_union_type *) & arg_b, &b);
1091
1092 if (isnan (&a) || isnan (&b))
1093 return -1; /* false, truth >= 0 */
1094 return _fpcmp_parts (&a, &b) ;
1095 }
1096
1097 CMPtype
_lt_f2(FLO_type arg_a,FLO_type arg_b)1098 _lt_f2 (FLO_type arg_a, FLO_type arg_b)
1099 {
1100 fp_number_type a;
1101 fp_number_type b;
1102
1103 unpack_d ((FLO_union_type *) & arg_a, &a);
1104 unpack_d ((FLO_union_type *) & arg_b, &b);
1105
1106 if (isnan (&a) || isnan (&b))
1107 return 1; /* false, truth < 0 */
1108
1109 return _fpcmp_parts (&a, &b);
1110 }
1111
1112 CMPtype
_le_f2(FLO_type arg_a,FLO_type arg_b)1113 _le_f2 (FLO_type arg_a, FLO_type arg_b)
1114 {
1115 fp_number_type a;
1116 fp_number_type b;
1117
1118 unpack_d ((FLO_union_type *) & arg_a, &a);
1119 unpack_d ((FLO_union_type *) & arg_b, &b);
1120
1121 if (isnan (&a) || isnan (&b))
1122 return 1; /* false, truth <= 0 */
1123
1124 return _fpcmp_parts (&a, &b) ;
1125 }
1126
1127 #endif /* ! US_SOFTWARE_GOFAST */
1128
1129 FLO_type
si_to_float(SItype arg_a)1130 si_to_float (SItype arg_a)
1131 {
1132 fp_number_type in;
1133
1134 in.class = CLASS_NUMBER;
1135 in.sign = arg_a < 0;
1136 if (!arg_a)
1137 {
1138 in.class = CLASS_ZERO;
1139 }
1140 else
1141 {
1142 in.normal_exp = FRACBITS + NGARDS;
1143 if (in.sign)
1144 {
1145 /* Special case for minint, since there is no +ve integer
1146 representation for it */
1147 if (arg_a == 0x80000000)
1148 {
1149 return -2147483648.0;
1150 }
1151 in.fraction.ll = (-arg_a);
1152 }
1153 else
1154 in.fraction.ll = arg_a;
1155
1156 while (in.fraction.ll < (1LL << (FRACBITS + NGARDS)))
1157 {
1158 in.fraction.ll <<= 1;
1159 in.normal_exp -= 1;
1160 }
1161 }
1162 return pack_d (&in);
1163 }
1164
1165 SItype
float_to_si(FLO_type arg_a)1166 float_to_si (FLO_type arg_a)
1167 {
1168 fp_number_type a;
1169 SItype tmp;
1170
1171 unpack_d ((FLO_union_type *) & arg_a, &a);
1172 if (iszero (&a))
1173 return 0;
1174 if (isnan (&a))
1175 return 0;
1176 /* get reasonable MAX_SI_INT... */
1177 if (isinf (&a))
1178 return a.sign ? MAX_SI_INT : (-MAX_SI_INT)-1;
1179 /* it is a number, but a small one */
1180 if (a.normal_exp < 0)
1181 return 0;
1182 if (a.normal_exp > 30)
1183 return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT;
1184 tmp = a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1185 return a.sign ? (-tmp) : (tmp);
1186 }
1187
1188 #ifdef US_SOFTWARE_GOFAST
1189 /* While libgcc2.c defines its own __fixunssfsi and __fixunsdfsi routines,
1190 we also define them for GOFAST because the ones in libgcc2.c have the
1191 wrong names and I'd rather define these here and keep GOFAST CYG-LOC's
1192 out of libgcc2.c. We can't define these here if not GOFAST because then
1193 there'd be duplicate copies. */
1194
1195 USItype
float_to_usi(FLO_type arg_a)1196 float_to_usi (FLO_type arg_a)
1197 {
1198 fp_number_type a;
1199
1200 unpack_d ((FLO_union_type *) & arg_a, &a);
1201 if (iszero (&a))
1202 return 0;
1203 if (isnan (&a))
1204 return 0;
1205 /* get reasonable MAX_USI_INT... */
1206 if (isinf (&a))
1207 return a.sign ? MAX_USI_INT : 0;
1208 /* it is a negative number */
1209 if (a.sign)
1210 return 0;
1211 /* it is a number, but a small one */
1212 if (a.normal_exp < 0)
1213 return 0;
1214 if (a.normal_exp > 31)
1215 return MAX_USI_INT;
1216 else if (a.normal_exp > (FRACBITS + NGARDS))
1217 return a.fraction.ll << ((FRACBITS + NGARDS) - a.normal_exp);
1218 else
1219 return a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
1220 }
1221 #endif
1222
1223 FLO_type
negate(FLO_type arg_a)1224 negate (FLO_type arg_a)
1225 {
1226 fp_number_type a;
1227
1228 unpack_d ((FLO_union_type *) & arg_a, &a);
1229 flip_sign (&a);
1230 return pack_d (&a);
1231 }
1232
1233 #ifdef FLOAT
1234
1235 SFtype
__make_fp(fp_class_type class,unsigned int sign,int exp,USItype frac)1236 __make_fp(fp_class_type class,
1237 unsigned int sign,
1238 int exp,
1239 USItype frac)
1240 {
1241 fp_number_type in;
1242
1243 in.class = class;
1244 in.sign = sign;
1245 in.normal_exp = exp;
1246 in.fraction.ll = frac;
1247 return pack_d (&in);
1248 }
1249
1250 #ifndef FLOAT_ONLY
1251
1252 /* This enables one to build an fp library that supports float but not double.
1253 Otherwise, we would get an undefined reference to __make_dp.
1254 This is needed for some 8-bit ports that can't handle well values that
1255 are 8-bytes in size, so we just don't support double for them at all. */
1256
1257 extern DFtype __make_dp (fp_class_type, unsigned int, int, UDItype frac);
1258
1259 DFtype
sf_to_df(SFtype arg_a)1260 sf_to_df (SFtype arg_a)
1261 {
1262 fp_number_type in;
1263
1264 unpack_d ((FLO_union_type *) & arg_a, &in);
1265 return __make_dp (in.class, in.sign, in.normal_exp,
1266 ((UDItype) in.fraction.ll) << F_D_BITOFF);
1267 }
1268
1269 #endif
1270 #endif
1271
1272 #ifndef FLOAT
1273
1274 extern SFtype __make_fp (fp_class_type, unsigned int, int, USItype);
1275
1276 DFtype
__make_dp(fp_class_type class,unsigned int sign,int exp,UDItype frac)1277 __make_dp (fp_class_type class, unsigned int sign, int exp, UDItype frac)
1278 {
1279 fp_number_type in;
1280
1281 in.class = class;
1282 in.sign = sign;
1283 in.normal_exp = exp;
1284 in.fraction.ll = frac;
1285 return pack_d (&in);
1286 }
1287
1288 SFtype
df_to_sf(DFtype arg_a)1289 df_to_sf (DFtype arg_a)
1290 {
1291 fp_number_type in;
1292
1293 unpack_d ((FLO_union_type *) & arg_a, &in);
1294 return __make_fp (in.class, in.sign, in.normal_exp,
1295 in.fraction.ll >> F_D_BITOFF);
1296 }
1297
1298 #endif
1299