Daily bump.
[official-gcc.git] / gcc / floatlib.c
blob929aef7353c1eb2bf5c76a76b1ab3d95c083b401
1 /*
2 ** libgcc support for software floating point.
3 ** Copyright (C) 1991 by Pipeline Associates, Inc. All rights reserved.
4 ** Permission is granted to do *anything* you want with this file,
5 ** commercial or otherwise, provided this message remains intact. So there!
6 ** I would appreciate receiving any updates/patches/changes that anyone
7 ** makes, and am willing to be the repository for said changes (am I
8 ** making a big mistake?).
10 Warning! Only single-precision is actually implemented. This file
11 won't really be much use until double-precision is supported.
13 However, once that is done, this file might make possible
14 cross-compilation for an IEEE target machine from a non-IEEE
15 host such as a VAX.
17 If you'd like to work on completing this, please talk to rms@gnu.ai.mit.edu.
19 --> Double precision floating support added by James Carlson on 20 April 1998.
22 ** Pat Wood
23 ** Pipeline Associates, Inc.
24 ** pipeline!phw@motown.com or
25 ** sun!pipeline!phw or
26 ** uunet!motown!pipeline!phw
28 ** 05/01/91 -- V1.0 -- first release to gcc mailing lists
29 ** 05/04/91 -- V1.1 -- added float and double prototypes and return values
30 ** -- fixed problems with adding and subtracting zero
31 ** -- fixed rounding in truncdfsf2
32 ** -- fixed SWAP define and tested on 386
36 ** The following are routines that replace the libgcc soft floating point
37 ** routines that are called automatically when -msoft-float is selected.
38 ** The support single and double precision IEEE format, with provisions
39 ** for byte-swapped machines (tested on 386). Some of the double-precision
40 ** routines work at full precision, but most of the hard ones simply punt
41 ** and call the single precision routines, producing a loss of accuracy.
42 ** long long support is not assumed or included.
43 ** Overall accuracy is close to IEEE (actually 68882) for single-precision
44 ** arithmetic. I think there may still be a 1 in 1000 chance of a bit
45 ** being rounded the wrong way during a multiply. I'm not fussy enough to
46 ** bother with it, but if anyone is, knock yourself out.
48 ** Efficiency has only been addressed where it was obvious that something
49 ** would make a big difference. Anyone who wants to do this right for
50 ** best speed should go in and rewrite in assembler.
52 ** I have tested this only on a 68030 workstation and 386/ix integrated
53 ** in with -msoft-float.
56 /* the following deal with IEEE single-precision numbers */
57 #define EXCESS 126
58 #define SIGNBIT 0x80000000
59 #define HIDDEN (1 << 23)
60 #define SIGN(fp) ((fp) & SIGNBIT)
61 #define EXP(fp) (((fp) >> 23) & 0xFF)
62 #define MANT(fp) (((fp) & 0x7FFFFF) | HIDDEN)
63 #define PACK(s,e,m) ((s) | ((e) << 23) | (m))
65 /* the following deal with IEEE double-precision numbers */
66 #define EXCESSD 1022
67 #define HIDDEND (1 << 20)
68 #define EXPD(fp) (((fp.l.upper) >> 20) & 0x7FF)
69 #define SIGND(fp) ((fp.l.upper) & SIGNBIT)
70 #define MANTD(fp) (((((fp.l.upper) & 0xFFFFF) | HIDDEND) << 10) | \
71 (fp.l.lower >> 22))
72 #define HIDDEND_LL ((long long)1 << 52)
73 #define MANTD_LL(fp) ((fp.ll & (HIDDEND_LL-1)) | HIDDEND_LL)
74 #define PACKD_LL(s,e,m) (((long long)((s)+((e)<<20))<<32)|(m))
76 /* define SWAP for 386/960 reverse-byte-order brain-damaged CPUs */
77 union double_long {
78 double d;
79 #ifdef SWAP
80 struct {
81 unsigned long lower;
82 long upper;
83 } l;
84 #else
85 struct {
86 long upper;
87 unsigned long lower;
88 } l;
89 #endif
90 long long ll;
93 union float_long
95 float f;
96 long l;
99 /* add two floats */
100 float
101 __addsf3 (float a1, float a2)
103 long mant1, mant2;
104 union float_long fl1, fl2;
105 int exp1, exp2;
106 int sign = 0;
108 fl1.f = a1;
109 fl2.f = a2;
111 /* check for zero args */
112 if (!fl1.l) {
113 fl1.f = fl2.f;
114 goto test_done;
116 if (!fl2.l)
117 goto test_done;
119 exp1 = EXP (fl1.l);
120 exp2 = EXP (fl2.l);
122 if (exp1 > exp2 + 25)
123 goto test_done;
124 if (exp2 > exp1 + 25) {
125 fl1.f = fl2.f;
126 goto test_done;
129 /* do everything in excess precision so's we can round later */
130 mant1 = MANT (fl1.l) << 6;
131 mant2 = MANT (fl2.l) << 6;
133 if (SIGN (fl1.l))
134 mant1 = -mant1;
135 if (SIGN (fl2.l))
136 mant2 = -mant2;
138 if (exp1 > exp2)
140 mant2 >>= exp1 - exp2;
142 else
144 mant1 >>= exp2 - exp1;
145 exp1 = exp2;
147 mant1 += mant2;
149 if (mant1 < 0)
151 mant1 = -mant1;
152 sign = SIGNBIT;
154 else if (!mant1) {
155 fl1.f = 0;
156 goto test_done;
159 /* normalize up */
160 while (!(mant1 & 0xE0000000))
162 mant1 <<= 1;
163 exp1--;
166 /* normalize down? */
167 if (mant1 & (1 << 30))
169 mant1 >>= 1;
170 exp1++;
173 /* round to even */
174 mant1 += (mant1 & 0x40) ? 0x20 : 0x1F;
176 /* normalize down? */
177 if (mant1 & (1 << 30))
179 mant1 >>= 1;
180 exp1++;
183 /* lose extra precision */
184 mant1 >>= 6;
186 /* turn off hidden bit */
187 mant1 &= ~HIDDEN;
189 /* pack up and go home */
190 fl1.l = PACK (sign, exp1, mant1);
191 test_done:
192 return (fl1.f);
195 /* subtract two floats */
196 float
197 __subsf3 (float a1, float a2)
199 union float_long fl1, fl2;
201 fl1.f = a1;
202 fl2.f = a2;
204 /* check for zero args */
205 if (!fl2.l)
206 return (fl1.f);
207 if (!fl1.l)
208 return (-fl2.f);
210 /* twiddle sign bit and add */
211 fl2.l ^= SIGNBIT;
212 return __addsf3 (a1, fl2.f);
215 /* compare two floats */
216 long
217 __cmpsf2 (float a1, float a2)
219 union float_long fl1, fl2;
221 fl1.f = a1;
222 fl2.f = a2;
224 if (SIGN (fl1.l) && SIGN (fl2.l))
226 fl1.l ^= SIGNBIT;
227 fl2.l ^= SIGNBIT;
229 if (fl1.l < fl2.l)
230 return (-1);
231 if (fl1.l > fl2.l)
232 return (1);
233 return (0);
236 /* multiply two floats */
237 float
238 __mulsf3 (float a1, float a2)
240 union float_long fl1, fl2;
241 unsigned long result;
242 int exp;
243 int sign;
245 fl1.f = a1;
246 fl2.f = a2;
248 if (!fl1.l || !fl2.l) {
249 fl1.f = 0;
250 goto test_done;
253 /* compute sign and exponent */
254 sign = SIGN (fl1.l) ^ SIGN (fl2.l);
255 exp = EXP (fl1.l) - EXCESS;
256 exp += EXP (fl2.l);
258 fl1.l = MANT (fl1.l);
259 fl2.l = MANT (fl2.l);
261 /* the multiply is done as one 16x16 multiply and two 16x8 multiples */
262 result = (fl1.l >> 8) * (fl2.l >> 8);
263 result += ((fl1.l & 0xFF) * (fl2.l >> 8)) >> 8;
264 result += ((fl2.l & 0xFF) * (fl1.l >> 8)) >> 8;
266 result >>= 2;
267 if (result & 0x20000000)
269 /* round */
270 result += 0x20;
271 result >>= 6;
273 else
275 /* round */
276 result += 0x10;
277 result >>= 5;
278 exp--;
280 if (result & (HIDDEN<<1)) {
281 result >>= 1;
282 exp++;
285 result &= ~HIDDEN;
287 /* pack up and go home */
288 fl1.l = PACK (sign, exp, result);
289 test_done:
290 return (fl1.f);
293 /* divide two floats */
294 float
295 __divsf3 (float a1, float a2)
297 union float_long fl1, fl2;
298 int result;
299 int mask;
300 int exp, sign;
302 fl1.f = a1;
303 fl2.f = a2;
305 /* subtract exponents */
306 exp = EXP (fl1.l) - EXP (fl2.l) + EXCESS;
308 /* compute sign */
309 sign = SIGN (fl1.l) ^ SIGN (fl2.l);
311 /* divide by zero??? */
312 if (!fl2.l)
313 /* return NaN or -NaN */
314 return (sign ? 0xFFFFFFFF : 0x7FFFFFFF);
316 /* numerator zero??? */
317 if (!fl1.l)
318 return (0);
320 /* now get mantissas */
321 fl1.l = MANT (fl1.l);
322 fl2.l = MANT (fl2.l);
324 /* this assures we have 25 bits of precision in the end */
325 if (fl1.l < fl2.l)
327 fl1.l <<= 1;
328 exp--;
331 /* now we perform repeated subtraction of fl2.l from fl1.l */
332 mask = 0x1000000;
333 result = 0;
334 while (mask)
336 if (fl1.l >= fl2.l)
338 result |= mask;
339 fl1.l -= fl2.l;
341 fl1.l <<= 1;
342 mask >>= 1;
345 /* round */
346 result += 1;
348 /* normalize down */
349 exp++;
350 result >>= 1;
352 result &= ~HIDDEN;
354 /* pack up and go home */
355 fl1.l = PACK (sign, exp, result);
356 return (fl1.f);
359 /* convert int to double */
360 double
361 __floatsidf (long a1)
363 int sign = 0, exp = 31 + EXCESSD;
364 union double_long dl;
366 if (!a1)
368 dl.l.upper = dl.l.lower = 0;
369 return (dl.d);
372 if (a1 < 0)
374 sign = SIGNBIT;
375 a1 = -a1;
378 while (a1 < 0x1000000)
380 a1 <<= 4;
381 exp -= 4;
384 while (a1 < 0x40000000)
386 a1 <<= 1;
387 exp--;
390 /* pack up and go home */
391 dl.l.upper = sign;
392 dl.l.upper |= exp << 20;
393 dl.l.upper |= (a1 >> 10) & ~HIDDEND;
394 dl.l.lower = a1 << 22;
396 return (dl.d);
399 double
400 __floatdidf (long long a1)
402 int exp = 63 + EXCESSD;
403 union double_long dl;
405 dl.l.upper = dl.l.lower = 0;
406 if (a1 == 0)
407 return (dl.d);
409 if (a1 < 0) {
410 dl.l.upper = SIGNBIT;
411 a1 = -a1;
414 while (a1 < (long long)1<<54) {
415 a1 <<= 8;
416 exp -= 8;
418 while (a1 < (long long)1<<62) {
419 a1 <<= 1;
420 exp -= 1;
423 /* pack up and go home */
424 dl.ll |= (a1 >> 10) & ~HIDDEND_LL;
425 dl.l.upper |= exp << 20;
427 return (dl.d);
430 float
431 __floatsisf (long a1)
433 (float)__floatsidf(a1);
436 float
437 __floatdisf (long long a1)
439 (float)__floatdidf(a1);
442 /* negate a float */
443 float
444 __negsf2 (float a1)
446 union float_long fl1;
448 fl1.f = a1;
449 if (!fl1.l)
450 return (0);
452 fl1.l ^= SIGNBIT;
453 return (fl1.f);
456 /* negate a double */
457 double
458 __negdf2 (double a1)
460 union double_long dl1;
462 dl1.d = a1;
464 if (!dl1.l.upper && !dl1.l.lower)
465 return (dl1.d);
467 dl1.l.upper ^= SIGNBIT;
468 return (dl1.d);
471 /* convert float to double */
472 double
473 __extendsfdf2 (float a1)
475 union float_long fl1;
476 union double_long dl;
477 int exp;
479 fl1.f = a1;
481 if (!fl1.l)
483 dl.l.upper = dl.l.lower = 0;
484 return (dl.d);
487 dl.l.upper = SIGN (fl1.l);
488 exp = EXP (fl1.l) - EXCESS + EXCESSD;
489 dl.l.upper |= exp << 20;
490 dl.l.upper |= (MANT (fl1.l) & ~HIDDEN) >> 3;
491 dl.l.lower = MANT (fl1.l) << 29;
493 return (dl.d);
496 /* convert double to float */
497 float
498 __truncdfsf2 (double a1)
500 int exp;
501 long mant;
502 union float_long fl;
503 union double_long dl1;
505 dl1.d = a1;
507 if (!dl1.l.upper && !dl1.l.lower)
508 return (float)(0);
510 exp = EXPD (dl1) - EXCESSD + EXCESS;
512 /* shift double mantissa 6 bits so we can round */
513 mant = MANTD (dl1) >> 6;
515 /* now round and shift down */
516 mant += 1;
517 mant >>= 1;
519 /* did the round overflow? */
520 if (mant & 0xFE000000)
522 mant >>= 1;
523 exp++;
526 mant &= ~HIDDEN;
528 /* pack up and go home */
529 fl.l = PACK (SIGND (dl1), exp, mant);
530 return (fl.f);
533 /* compare two doubles */
534 long
535 __cmpdf2 (double a1, double a2)
537 union double_long dl1, dl2;
539 dl1.d = a1;
540 dl2.d = a2;
542 if (SIGND (dl1) && SIGND (dl2))
544 dl1.l.upper ^= SIGNBIT;
545 dl2.l.upper ^= SIGNBIT;
547 if (dl1.l.upper < dl2.l.upper)
548 return (-1);
549 if (dl1.l.upper > dl2.l.upper)
550 return (1);
551 if (dl1.l.lower < dl2.l.lower)
552 return (-1);
553 if (dl1.l.lower > dl2.l.lower)
554 return (1);
555 return (0);
558 /* convert double to int */
559 long
560 __fixdfsi (double a1)
562 union double_long dl1;
563 int exp;
564 long l;
566 dl1.d = a1;
568 if (!dl1.l.upper && !dl1.l.lower)
569 return (0);
571 exp = EXPD (dl1) - EXCESSD - 31;
572 l = MANTD (dl1);
574 if (exp > 0)
575 return SIGND(dl1) ? (1<<31) : ((1ul<<31)-1);
577 /* shift down until exp = 0 or l = 0 */
578 if (exp < 0 && exp > -32 && l)
579 l >>= -exp;
580 else
581 return (0);
583 return (SIGND (dl1) ? -l : l);
586 /* convert double to int */
587 long long
588 __fixdfdi (double a1)
590 union double_long dl1;
591 int exp;
592 long long l;
594 dl1.d = a1;
596 if (!dl1.l.upper && !dl1.l.lower)
597 return (0);
599 exp = EXPD (dl1) - EXCESSD - 64;
600 l = MANTD_LL(dl1);
602 if (exp > 0) {
603 l = (long long)1<<63;
604 if (!SIGND(dl1))
605 l--;
606 return l;
609 /* shift down until exp = 0 or l = 0 */
610 if (exp < 0 && exp > -64 && l)
611 l >>= -exp;
612 else
613 return (0);
615 return (SIGND (dl1) ? -l : l);
618 /* convert double to unsigned int */
619 unsigned long
620 __fixunsdfsi (double a1)
622 union double_long dl1;
623 int exp;
624 unsigned long l;
626 dl1.d = a1;
628 if (!dl1.l.upper && !dl1.l.lower)
629 return (0);
631 exp = EXPD (dl1) - EXCESSD - 32;
632 l = (((((dl1.l.upper) & 0xFFFFF) | HIDDEND) << 11) | (dl1.l.lower >> 21));
634 if (exp > 0)
635 return (0xFFFFFFFFul); /* largest integer */
637 /* shift down until exp = 0 or l = 0 */
638 if (exp < 0 && exp > -32 && l)
639 l >>= -exp;
640 else
641 return (0);
643 return (l);
646 /* convert double to unsigned int */
647 unsigned long long
648 __fixunsdfdi (double a1)
650 union double_long dl1;
651 int exp;
652 unsigned long long l;
654 dl1.d = a1;
656 if (dl1.ll == 0)
657 return (0);
659 exp = EXPD (dl1) - EXCESSD - 64;
661 l = dl1.ll;
663 if (exp > 0)
664 return (unsigned long long)-1;
666 /* shift down until exp = 0 or l = 0 */
667 if (exp < 0 && exp > -64 && l)
668 l >>= -exp;
669 else
670 return (0);
672 return (l);
675 /* addtwo doubles */
676 double
677 __adddf3 (double a1, double a2)
679 long long mant1, mant2;
680 union double_long fl1, fl2;
681 int exp1, exp2;
682 int sign = 0;
684 fl1.d = a1;
685 fl2.d = a2;
687 /* check for zero args */
688 if (!fl2.ll)
689 goto test_done;
690 if (!fl1.ll) {
691 fl1.d = fl2.d;
692 goto test_done;
695 exp1 = EXPD(fl1);
696 exp2 = EXPD(fl2);
698 if (exp1 > exp2 + 54)
699 goto test_done;
700 if (exp2 > exp1 + 54) {
701 fl1.d = fl2.d;
702 goto test_done;
705 /* do everything in excess precision so's we can round later */
706 mant1 = MANTD_LL(fl1) << 9;
707 mant2 = MANTD_LL(fl2) << 9;
709 if (SIGND(fl1))
710 mant1 = -mant1;
711 if (SIGND(fl2))
712 mant2 = -mant2;
714 if (exp1 > exp2)
715 mant2 >>= exp1 - exp2;
716 else {
717 mant1 >>= exp2 - exp1;
718 exp1 = exp2;
720 mant1 += mant2;
722 if (mant1 < 0) {
723 mant1 = -mant1;
724 sign = SIGNBIT;
725 } else if (!mant1) {
726 fl1.d = 0;
727 goto test_done;
730 /* normalize up */
731 while (!(mant1 & ((long long)7<<61))) {
732 mant1 <<= 1;
733 exp1--;
736 /* normalize down? */
737 if (mant1 & ((long long)3<<62)) {
738 mant1 >>= 1;
739 exp1++;
742 /* round to even */
743 mant1 += (mant1 & (1<<9)) ? (1<<8) : ((1<<8)-1);
745 /* normalize down? */
746 if (mant1 & ((long long)3<<62)) {
747 mant1 >>= 1;
748 exp1++;
751 /* lose extra precision */
752 mant1 >>= 9;
754 /* turn off hidden bit */
755 mant1 &= ~HIDDEND_LL;
757 /* pack up and go home */
758 fl1.ll = PACKD_LL(sign,exp1,mant1);
760 test_done:
761 return (fl1.d);
764 /* subtract two doubles */
765 double
766 __subdf3 (double a1, double a2)
768 union double_long fl1, fl2;
770 fl1.d = a1;
771 fl2.d = a2;
773 /* check for zero args */
774 if (!fl2.ll)
775 return (fl1.d);
776 /* twiddle sign bit and add */
777 fl2.l.upper ^= SIGNBIT;
778 if (!fl1.ll)
779 return (fl2.d);
780 return __adddf3 (a1, fl2.d);
783 /* multiply two doubles */
784 double
785 __muldf3 (double a1, double a2)
787 union double_long fl1, fl2;
788 unsigned long long result;
789 int exp;
790 int sign;
792 fl1.d = a1;
793 fl2.d = a2;
795 if (!fl1.ll || !fl2.ll) {
796 fl1.d = 0;
797 goto test_done;
800 /* compute sign and exponent */
801 sign = SIGND(fl1) ^ SIGND(fl2);
802 exp = EXPD(fl1) - EXCESSD;
803 exp += EXPD(fl2);
805 fl1.ll = MANTD_LL(fl1);
806 fl2.ll = MANTD_LL(fl2);
808 /* the multiply is done as one 31x31 multiply and two 31x21 multiples */
809 result = (fl1.ll >> 21) * (fl2.ll >> 21);
810 result += ((fl1.ll & 0x1FFFFF) * (fl2.ll >> 21)) >> 21;
811 result += ((fl2.ll & 0x1FFFFF) * (fl1.ll >> 21)) >> 21;
813 result >>= 2;
814 if (result & ((long long)1<<61)) {
815 /* round */
816 result += 1<<8;
817 result >>= 9;
818 } else {
819 /* round */
820 result += 1<<7;
821 result >>= 8;
822 exp--;
824 if (result & (HIDDEND_LL<<1)) {
825 result >>= 1;
826 exp++;
829 result &= ~HIDDEND_LL;
831 /* pack up and go home */
832 fl1.ll = PACKD_LL(sign,exp,result);
833 test_done:
834 return (fl1.d);
837 /* divide two doubles */
838 double
839 __divdf3 (double a1, double a2)
841 union double_long fl1, fl2;
842 long long mask,result;
843 int exp, sign;
845 fl1.d = a1;
846 fl2.d = a2;
848 /* subtract exponents */
849 exp = EXPD(fl1) - EXPD(fl2) + EXCESSD;
851 /* compute sign */
852 sign = SIGND(fl1) ^ SIGND(fl2);
854 /* numerator zero??? */
855 if (fl1.ll == 0) {
856 /* divide by zero??? */
857 if (fl2.ll == 0)
858 fl1.ll = ((unsigned long long)1<<63)-1; /* NaN */
859 else
860 fl1.ll = 0;
861 goto test_done;
864 /* return +Inf or -Inf */
865 if (fl2.ll == 0) {
866 fl1.ll = PACKD_LL(SIGND(fl1),2047,0);
867 goto test_done;
871 /* now get mantissas */
872 fl1.ll = MANTD_LL(fl1);
873 fl2.ll = MANTD_LL(fl2);
875 /* this assures we have 54 bits of precision in the end */
876 if (fl1.ll < fl2.ll) {
877 fl1.ll <<= 1;
878 exp--;
881 /* now we perform repeated subtraction of fl2.ll from fl1.ll */
882 mask = (long long)1<<53;
883 result = 0;
884 while (mask) {
885 if (fl1.ll >= fl2.ll)
887 result |= mask;
888 fl1.ll -= fl2.ll;
890 fl1.ll <<= 1;
891 mask >>= 1;
894 /* round */
895 result += 1;
897 /* normalize down */
898 exp++;
899 result >>= 1;
901 result &= ~HIDDEND_LL;
903 /* pack up and go home */
904 fl1.ll = PACKD_LL(sign, exp, result);
906 test_done:
907 return (fl1.d);
911 __gtdf2 (double a1, double a2)
913 return __cmpdf2 ((float) a1, (float) a2) > 0;
917 __gedf2 (double a1, double a2)
919 return (__cmpdf2 ((float) a1, (float) a2) >= 0) - 1;
923 __ltdf2 (double a1, double a2)
925 return - (__cmpdf2 ((float) a1, (float) a2) < 0);
929 __ledf2 (double a1, double a2)
931 return __cmpdf2 ((float) a1, (float) a2) > 0;
935 __eqdf2 (double a1, double a2)
937 return *(long long *) &a1 == *(long long *) &a2;
941 __nedf2 (double a1, double a2)
943 return *(long long *) &a1 != *(long long *) &a2;