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 eventually become a
14 replacement for libgcc1.c. It might also make possible
15 cross-compilation for an IEEE target machine from a non-IEEE
18 If you'd like to work on completing this, please talk to rms@gnu.ai.mit.edu.
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 D_PHANTOM_BIT 0x00100000
59 #define SIGNBIT 0x80000000
60 #define HIDDEN (1 << 23)
61 #define SIGN(fp) ((fp) & SIGNBIT)
62 #define EXP(fp) (((fp) >> 23) & 0xFF)
63 #define MANT(fp) (((fp) & 0x7FFFFF) | HIDDEN)
64 #define PACK(s,e,m) ((s) | ((e) << 23) | (m))
66 /* the following deal with IEEE double-precision numbers */
68 #define HIDDEND (1 << 20)
69 #define EXPD(fp) (((fp.l.upper) >> 20) & 0x7FF)
70 #define SIGND(fp) ((fp.l.upper) & SIGNBIT)
71 #define MANTD(fp) (((((fp.l.upper) & 0xFFFFF) | HIDDEND) << 10) | \
74 /* define SWAP for 386/960 reverse-byte-order brain-damaged CPUs */
99 unsigned mantissa2 : 32;
100 unsigned mantissa1 : 20;
101 unsigned exponent : 11;
104 unsigned exponent : 11;
106 unsigned mantissa2 : 32;
107 unsigned mantissa1 : 20;
130 __addsf3 (float a1, float a2)
132 register long mant1, mant2;
133 register union float_long fl1, fl2;
134 register int exp1, exp2;
140 /* check for zero args */
149 if (exp1 > exp2 + 25)
151 if (exp2 > exp1 + 25)
154 /* do everything in excess precision so's we can round later */
155 mant1 = MANT (fl1.l) << 6;
156 mant2 = MANT (fl2.l) << 6;
165 mant2 >>= exp1 - exp2;
169 mant1 >>= exp2 - exp1;
183 while (!(mant1 & 0xE0000000))
189 /* normalize down? */
190 if (mant1 & (1 << 30))
197 mant1 += (mant1 & 0x40) ? 0x20 : 0x1F;
199 /* normalize down? */
200 if (mant1 & (1 << 30))
206 /* lose extra precision */
209 /* turn off hidden bit */
212 /* pack up and go home */
213 fl1.l = PACK (sign, exp1, mant1);
217 /* subtract two floats */
220 __subsf3 (float a1, float a2)
222 register union float_long fl1, fl2;
227 /* check for zero args */
233 /* twiddle sign bit and add */
235 return __addsf3 (a1, fl2.f);
238 /* compare two floats */
241 __cmpsf2 (float a1, float a2)
243 register union float_long fl1, fl2;
248 if (SIGN (fl1.l) && SIGN (fl2.l))
260 /* multiply two floats */
263 __mulsf3 (float a1, float a2)
265 register union float_long fl1, fl2;
266 register unsigned long result;
273 if (!fl1.l || !fl2.l)
276 /* compute sign and exponent */
277 sign = SIGN (fl1.l) ^ SIGN (fl2.l);
278 exp = EXP (fl1.l) - EXCESS;
281 fl1.l = MANT (fl1.l);
282 fl2.l = MANT (fl2.l);
284 /* the multiply is done as one 16x16 multiply and two 16x8 multiples */
285 result = (fl1.l >> 8) * (fl2.l >> 8);
286 result += ((fl1.l & 0xFF) * (fl2.l >> 8)) >> 8;
287 result += ((fl2.l & 0xFF) * (fl1.l >> 8)) >> 8;
289 if (result & 0x80000000)
305 /* pack up and go home */
306 fl1.l = PACK (sign, exp, result);
310 /* divide two floats */
313 __divsf3 (float a1, float a2)
315 register union float_long fl1, fl2;
318 register int exp, sign;
323 /* subtract exponents */
324 exp = EXP (fl1.l) - EXP (fl2.l) + EXCESS;
327 sign = SIGN (fl1.l) ^ SIGN (fl2.l);
329 /* divide by zero??? */
331 /* return NaN or -NaN */
332 return (sign ? 0xFFFFFFFF : 0x7FFFFFFF);
334 /* numerator zero??? */
338 /* now get mantissas */
339 fl1.l = MANT (fl1.l);
340 fl2.l = MANT (fl2.l);
342 /* this assures we have 25 bits of precision in the end */
349 /* now we perform repeated subtraction of fl2.l from fl1.l */
372 /* pack up and go home */
373 fl1.l = PACK (sign, exp, result);
377 /* convert int to double */
380 __floatsidf (register long a1)
382 register int sign = 0, exp = 31 + EXCESSD;
383 union double_long dl;
387 dl.l.upper = dl.l.lower = 0;
397 while (a1 < 0x1000000)
403 while (a1 < 0x40000000)
409 /* pack up and go home */
411 dl.l.upper |= exp << 20;
412 dl.l.upper |= (a1 >> 10) & ~HIDDEND;
413 dl.l.lower = a1 << 22;
423 register union float_long fl1;
433 /* negate a double */
438 register union double_long dl1;
442 if (!dl1.l.upper && !dl1.l.lower)
445 dl1.l.upper ^= SIGNBIT;
449 /* convert float to double */
452 __extendsfdf2 (float a1)
454 register union float_long fl1;
455 register union double_long dl;
462 dl.l.upper = dl.l.lower = 0;
466 dl.l.upper = SIGN (fl1.l);
467 exp = EXP (fl1.l) - EXCESS + EXCESSD;
468 dl.l.upper |= exp << 20;
469 dl.l.upper |= (MANT (fl1.l) & ~HIDDEN) >> 3;
470 dl.l.lower = MANT (fl1.l) << 29;
475 /* convert double to float */
478 __truncdfsf2 (double a1)
482 register union float_long fl;
483 register union double_long dl1;
487 if (!dl1.l.upper && !dl1.l.lower)
490 exp = EXPD (dl1) - EXCESSD + EXCESS;
492 /* shift double mantissa 6 bits so we can round */
493 mant = MANTD (dl1) >> 6;
495 /* now round and shift down */
499 /* did the round overflow? */
500 if (mant & 0xFF000000)
508 /* pack up and go home */
509 fl.l = PACK (SIGND (dl1), exp, mant);
513 /* compare two doubles */
516 __cmpdf2 (double a1, double a2)
518 register union double_long dl1, dl2;
523 if (SIGND (dl1) && SIGND (dl2))
525 dl1.l.upper ^= SIGNBIT;
526 dl2.l.upper ^= SIGNBIT;
528 if (dl1.l.upper < dl2.l.upper)
530 if (dl1.l.upper > dl2.l.upper)
532 if (dl1.l.lower < dl2.l.lower)
534 if (dl1.l.lower > dl2.l.lower)
539 /* convert double to int */
542 __fixdfsi (double a1)
544 register union double_long dl1;
550 if (!dl1.l.upper && !dl1.l.lower)
553 exp = EXPD (dl1) - EXCESSD - 31;
557 return (0x7FFFFFFF | SIGND (dl1)); /* largest integer */
559 /* shift down until exp = 0 or l = 0 */
560 if (exp < 0 && exp > -32 && l)
565 return (SIGND (dl1) ? -l : l);
568 /* convert double to unsigned int */
571 long __fixunsdfsi (double a1)
573 register union double_long dl1;
575 register unsigned long l;
579 if (!dl1.l.upper && !dl1.l.lower)
582 exp = EXPD (dl1) - EXCESSD - 32;
583 l = (((((dl1.l.upper) & 0xFFFFF) | HIDDEND) << 11) | (dl1.l.lower >> 21));
586 return (0xFFFFFFFF); /* largest integer */
588 /* shift down until exp = 0 or l = 0 */
589 if (exp < 0 && exp > -32 && l)
597 /* For now, the hard double-precision routines simply
598 punt and do it in single */
602 __adddf3 (double a1, double a2)
604 return ((float) a1 + (float) a2);
607 /* subtract two doubles */
610 __subdf3 (double a1, double a2)
612 return ((float) a1 - (float) a2);
615 /* multiply two doubles */
618 __muldf3 (double a1, double a2)
620 return ((float) a1 * (float) a2);
625 * Name: Barrett Richardson
626 * E-mail: barrett@iglou.com
627 * When: Thu Dec 15 10:31:11 EST 1994
631 * double __divdf3(double a1, double a2);
633 * Does software divide of a1 / a2.
635 * Based largely on __divsf3() in floatlib.c in the gcc
638 * Purpose: To be used in conjunction with the -msoft-float
639 * option of gcc. You should be able to tack it to the
640 * end of floatlib.c included in the gcc distribution,
641 * and delete the __divdf3() already there which just
642 * calls the single precision function (or may just
643 * use the floating point processor with some configurations).
645 * You may use this code for whatever your heart desires.
652 * Compare the mantissas of two doubles.
653 * Each mantissa is in two longs.
655 * return 1 if x1's mantissa is greater than x2's
656 * -1 if x1's mantissa is less than x2's
657 * 0 if the two mantissa's are equal.
659 * The Mantissas won't fit into a 4 byte word, so they are
660 * broken up into two parts.
662 * This function is used internally by __divdf3()
666 __dcmp (long x1m1, long x1m2, long x2m1, long x2m2)
674 /* If the first word in the two mantissas were equal check the second word */
686 /* divide two doubles */
689 __divdf3 (double a1, double a2)
696 register unsigned long mantissa1,
712 exponent = x1.ieee.exponent - x2.ieee.exponent + EXCESSD;
714 sign = x1.ieee.sign ^ x2.ieee.sign;
716 x2.ieee.sign = 0; /* don't want the sign bit to affect any zero */
717 /* comparisons when checking for zero divide */
719 if (!x2.l.lower && !x2.l.upper) { /* check for zero divide */
720 result.l.lower = 0x0;
722 result.l.upper = 0xFFF00000; /* negative infinity */
724 result.l.upper = 0x7FF00000; /* positive infinity */
728 if (!x1.l.upper && !x1.l.lower) /* check for 0.0 numerator */
731 x1m1 = x1.ieee.mantissa1 | D_PHANTOM_BIT; /* turn on phantom bit */
732 x1m2 = x1.ieee.mantissa2;
734 x2m1 = x2.ieee.mantissa1 | D_PHANTOM_BIT; /* turn on phantom bit */
735 x2m2 = x2.ieee.mantissa2;
737 if (__dcmp(x1m1,x1m2,x2m1,x2m2) < 0) {
739 /* if x1's mantissa is less than x2's shift it left one and decrement */
740 /* the exponent to accommodate the change in the mantissa */
743 bit_bucket = x1m2 >> 31; /* Shift mantissa left one */
744 x1m1 |= bit_bucket; /* */
755 /* Get the first part of the results mantissa using successive */
761 if (__dcmp(x1m1,x1m2,x2m1,x2m2) >= 0) {
763 /* subtract x2's mantissa from x1's */
765 mantissa1 |= mask; /* turn on a bit in the result */
774 bit_bucket = x1m2 >> 31; /* Shift mantissa left one */
775 x1m1 |= bit_bucket; /* */
781 /* Get the second part of the results mantissa using successive */
787 if (__dcmp(x1m1,x1m2,x2m1,x2m2) >= 0) {
789 /* subtract x2's mantissa from x1's */
791 mantissa2 |= mask; /* turn on a bit in the result */
799 bit_bucket = x1m2 >> 31; /* Shift mantissa left one */
800 x1m1 |= bit_bucket; /* */
806 /* round up by adding 1 to mantissa */
808 if (mantissa2 == 0xFFFFFFFF) { /* check for over flow */
810 /* spill if overflow */
818 exponent++; /* increment exponent (mantissa must be shifted right */
821 /* shift mantissa right one and assume a phantom bit (which really gives */
822 /* 53 bits of precision in the mantissa) */
825 bit_bucket = mantissa1 & 1;
826 mantissa2 |= (bit_bucket << 31);
829 /* put all the info into the result */
831 result.ieee.exponent = exponent;
832 result.ieee.sign = sign;
833 result.ieee.mantissa1 = mantissa1;
834 result.ieee.mantissa2 = mantissa2;