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
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.
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 */
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 */
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) | \
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 */
101 __addsf3 (float a1, float a2)
104 union float_long fl1, fl2;
111 /* check for zero args */
122 if (exp1 > exp2 + 25)
124 if (exp2 > exp1 + 25) {
129 /* do everything in excess precision so's we can round later */
130 mant1 = MANT (fl1.l) << 6;
131 mant2 = MANT (fl2.l) << 6;
140 mant2 >>= exp1 - exp2;
144 mant1 >>= exp2 - exp1;
160 while (!(mant1 & 0xE0000000))
166 /* normalize down? */
167 if (mant1 & (1 << 30))
174 mant1 += (mant1 & 0x40) ? 0x20 : 0x1F;
176 /* normalize down? */
177 if (mant1 & (1 << 30))
183 /* lose extra precision */
186 /* turn off hidden bit */
189 /* pack up and go home */
190 fl1.l = PACK (sign, exp1, mant1);
195 /* subtract two floats */
197 __subsf3 (float a1, float a2)
199 union float_long fl1, fl2;
204 /* check for zero args */
210 /* twiddle sign bit and add */
212 return __addsf3 (a1, fl2.f);
215 /* compare two floats */
217 __cmpsf2 (float a1, float a2)
219 union float_long fl1, fl2;
224 if (SIGN (fl1.l) && SIGN (fl2.l))
236 /* multiply two floats */
238 __mulsf3 (float a1, float a2)
240 union float_long fl1, fl2;
241 unsigned long result;
248 if (!fl1.l || !fl2.l) {
253 /* compute sign and exponent */
254 sign = SIGN (fl1.l) ^ SIGN (fl2.l);
255 exp = EXP (fl1.l) - EXCESS;
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;
267 if (result & 0x20000000)
280 if (result & (HIDDEN<<1)) {
287 /* pack up and go home */
288 fl1.l = PACK (sign, exp, result);
293 /* divide two floats */
295 __divsf3 (float a1, float a2)
297 union float_long fl1, fl2;
305 /* subtract exponents */
306 exp = EXP (fl1.l) - EXP (fl2.l) + EXCESS;
309 sign = SIGN (fl1.l) ^ SIGN (fl2.l);
311 /* divide by zero??? */
313 /* return NaN or -NaN */
314 return (sign ? 0xFFFFFFFF : 0x7FFFFFFF);
316 /* numerator zero??? */
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 */
331 /* now we perform repeated subtraction of fl2.l from fl1.l */
354 /* pack up and go home */
355 fl1.l = PACK (sign, exp, result);
359 /* convert int to double */
361 __floatsidf (long a1)
363 int sign = 0, exp = 31 + EXCESSD;
364 union double_long dl;
368 dl.l.upper = dl.l.lower = 0;
378 while (a1 < 0x1000000)
384 while (a1 < 0x40000000)
390 /* pack up and go home */
392 dl.l.upper |= exp << 20;
393 dl.l.upper |= (a1 >> 10) & ~HIDDEND;
394 dl.l.lower = a1 << 22;
400 __floatdidf (long long a1)
402 int exp = 63 + EXCESSD;
403 union double_long dl;
405 dl.l.upper = dl.l.lower = 0;
410 dl.l.upper = SIGNBIT;
414 while (a1 < (long long)1<<54) {
418 while (a1 < (long long)1<<62) {
423 /* pack up and go home */
424 dl.ll |= (a1 >> 10) & ~HIDDEND_LL;
425 dl.l.upper |= exp << 20;
431 __floatsisf (long a1)
433 (float)__floatsidf(a1);
437 __floatdisf (long long a1)
439 (float)__floatdidf(a1);
446 union float_long fl1;
456 /* negate a double */
460 union double_long dl1;
464 if (!dl1.l.upper && !dl1.l.lower)
467 dl1.l.upper ^= SIGNBIT;
471 /* convert float to double */
473 __extendsfdf2 (float a1)
475 union float_long fl1;
476 union double_long dl;
483 dl.l.upper = dl.l.lower = 0;
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;
496 /* convert double to float */
498 __truncdfsf2 (double a1)
503 union double_long dl1;
507 if (!dl1.l.upper && !dl1.l.lower)
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 */
519 /* did the round overflow? */
520 if (mant & 0xFE000000)
528 /* pack up and go home */
529 fl.l = PACK (SIGND (dl1), exp, mant);
533 /* compare two doubles */
535 __cmpdf2 (double a1, double a2)
537 union double_long dl1, dl2;
542 if (SIGND (dl1) && SIGND (dl2))
544 dl1.l.upper ^= SIGNBIT;
545 dl2.l.upper ^= SIGNBIT;
547 if (dl1.l.upper < dl2.l.upper)
549 if (dl1.l.upper > dl2.l.upper)
551 if (dl1.l.lower < dl2.l.lower)
553 if (dl1.l.lower > dl2.l.lower)
558 /* convert double to int */
560 __fixdfsi (double a1)
562 union double_long dl1;
568 if (!dl1.l.upper && !dl1.l.lower)
571 exp = EXPD (dl1) - EXCESSD - 31;
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)
583 return (SIGND (dl1) ? -l : l);
586 /* convert double to int */
588 __fixdfdi (double a1)
590 union double_long dl1;
596 if (!dl1.l.upper && !dl1.l.lower)
599 exp = EXPD (dl1) - EXCESSD - 64;
603 l = (long long)1<<63;
609 /* shift down until exp = 0 or l = 0 */
610 if (exp < 0 && exp > -64 && l)
615 return (SIGND (dl1) ? -l : l);
618 /* convert double to unsigned int */
620 __fixunsdfsi (double a1)
622 union double_long dl1;
628 if (!dl1.l.upper && !dl1.l.lower)
631 exp = EXPD (dl1) - EXCESSD - 32;
632 l = (((((dl1.l.upper) & 0xFFFFF) | HIDDEND) << 11) | (dl1.l.lower >> 21));
635 return (0xFFFFFFFFul); /* largest integer */
637 /* shift down until exp = 0 or l = 0 */
638 if (exp < 0 && exp > -32 && l)
646 /* convert double to unsigned int */
648 __fixunsdfdi (double a1)
650 union double_long dl1;
652 unsigned long long l;
659 exp = EXPD (dl1) - EXCESSD - 64;
664 return (unsigned long long)-1;
666 /* shift down until exp = 0 or l = 0 */
667 if (exp < 0 && exp > -64 && l)
677 __adddf3 (double a1, double a2)
679 long long mant1, mant2;
680 union double_long fl1, fl2;
687 /* check for zero args */
698 if (exp1 > exp2 + 54)
700 if (exp2 > exp1 + 54) {
705 /* do everything in excess precision so's we can round later */
706 mant1 = MANTD_LL(fl1) << 9;
707 mant2 = MANTD_LL(fl2) << 9;
715 mant2 >>= exp1 - exp2;
717 mant1 >>= exp2 - exp1;
731 while (!(mant1 & ((long long)7<<61))) {
736 /* normalize down? */
737 if (mant1 & ((long long)3<<62)) {
743 mant1 += (mant1 & (1<<9)) ? (1<<8) : ((1<<8)-1);
745 /* normalize down? */
746 if (mant1 & ((long long)3<<62)) {
751 /* lose extra precision */
754 /* turn off hidden bit */
755 mant1 &= ~HIDDEND_LL;
757 /* pack up and go home */
758 fl1.ll = PACKD_LL(sign,exp1,mant1);
764 /* subtract two doubles */
766 __subdf3 (double a1, double a2)
768 union double_long fl1, fl2;
773 /* check for zero args */
776 /* twiddle sign bit and add */
777 fl2.l.upper ^= SIGNBIT;
780 return __adddf3 (a1, fl2.d);
783 /* multiply two doubles */
785 __muldf3 (double a1, double a2)
787 union double_long fl1, fl2;
788 unsigned long long result;
795 if (!fl1.ll || !fl2.ll) {
800 /* compute sign and exponent */
801 sign = SIGND(fl1) ^ SIGND(fl2);
802 exp = EXPD(fl1) - EXCESSD;
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;
814 if (result & ((long long)1<<61)) {
824 if (result & (HIDDEND_LL<<1)) {
829 result &= ~HIDDEND_LL;
831 /* pack up and go home */
832 fl1.ll = PACKD_LL(sign,exp,result);
837 /* divide two doubles */
839 __divdf3 (double a1, double a2)
841 union double_long fl1, fl2;
842 long long mask,result;
848 /* subtract exponents */
849 exp = EXPD(fl1) - EXPD(fl2) + EXCESSD;
852 sign = SIGND(fl1) ^ SIGND(fl2);
854 /* numerator zero??? */
856 /* divide by zero??? */
858 fl1.ll = ((unsigned long long)1<<63)-1; /* NaN */
864 /* return +Inf or -Inf */
866 fl1.ll = PACKD_LL(SIGND(fl1),2047,0);
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) {
881 /* now we perform repeated subtraction of fl2.ll from fl1.ll */
882 mask = (long long)1<<53;
885 if (fl1.ll >= fl2.ll)
901 result &= ~HIDDEND_LL;
903 /* pack up and go home */
904 fl1.ll = PACKD_LL(sign, exp, result);
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;