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 */
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) | \
73 /* define SWAP for 386/960 reverse-byte-order brain-damaged CPUs */
98 __addsf3 (float a1, float a2)
100 register long mant1, mant2;
101 register union float_long fl1, fl2;
102 register int exp1, exp2;
108 /* check for zero args */
117 if (exp1 > exp2 + 25)
119 if (exp2 > exp1 + 25)
122 /* do everything in excess precision so's we can round later */
123 mant1 = MANT (fl1.l) << 6;
124 mant2 = MANT (fl2.l) << 6;
133 mant2 >>= exp1 - exp2;
137 mant1 >>= exp2 - exp1;
151 while (!(mant1 & 0xE0000000))
157 /* normalize down? */
158 if (mant1 & (1 << 30))
165 mant1 += (mant1 & 0x40) ? 0x20 : 0x1F;
167 /* normalize down? */
168 if (mant1 & (1 << 30))
174 /* lose extra precision */
177 /* turn off hidden bit */
180 /* pack up and go home */
181 fl1.l = PACK (sign, exp1, mant1);
185 /* subtract two floats */
187 __subsf3 (float a1, float a2)
189 register union float_long fl1, fl2;
194 /* check for zero args */
200 /* twiddle sign bit and add */
202 return __addsf3 (a1, fl2.f);
205 /* compare two floats */
207 __cmpsf2 (float a1, float a2)
209 register union float_long fl1, fl2;
214 if (SIGN (fl1.l) && SIGN (fl2.l))
226 /* multiply two floats */
228 __mulsf3 (float a1, float a2)
230 register union float_long fl1, fl2;
231 register unsigned long result;
238 if (!fl1.l || !fl2.l)
241 /* compute sign and exponent */
242 sign = SIGN (fl1.l) ^ SIGN (fl2.l);
243 exp = EXP (fl1.l) - EXCESS;
246 fl1.l = MANT (fl1.l);
247 fl2.l = MANT (fl2.l);
249 /* the multiply is done as one 16x16 multiply and two 16x8 multiples */
250 result = (fl1.l >> 8) * (fl2.l >> 8);
251 result += ((fl1.l & 0xFF) * (fl2.l >> 8)) >> 8;
252 result += ((fl2.l & 0xFF) * (fl1.l >> 8)) >> 8;
254 if (result & 0x80000000)
270 /* pack up and go home */
271 fl1.l = PACK (sign, exp, result);
275 /* divide two floats */
277 __divsf3 (float a1, float a2)
279 register union float_long fl1, fl2;
282 register int exp, sign;
287 /* subtract exponents */
288 exp = EXP (fl1.l) - EXP (fl2.l) + EXCESS;
291 sign = SIGN (fl1.l) ^ SIGN (fl2.l);
293 /* divide by zero??? */
295 /* return NaN or -NaN */
296 return (sign ? 0xFFFFFFFF : 0x7FFFFFFF);
298 /* numerator zero??? */
302 /* now get mantissas */
303 fl1.l = MANT (fl1.l);
304 fl2.l = MANT (fl2.l);
306 /* this assures we have 25 bits of precision in the end */
313 /* now we perform repeated subtraction of fl2.l from fl1.l */
336 /* pack up and go home */
337 fl1.l = PACK (sign, exp, result);
341 /* convert int to double */
343 __floatsidf (register long a1)
345 register int sign = 0, exp = 31 + EXCESSD;
346 union double_long dl;
350 dl.l.upper = dl.l.lower = 0;
360 while (a1 < 0x1000000)
366 while (a1 < 0x40000000)
372 /* pack up and go home */
374 dl.l.upper |= exp << 20;
375 dl.l.upper |= (a1 >> 10) & ~HIDDEND;
376 dl.l.lower = a1 << 22;
385 register union float_long fl1;
395 /* negate a double */
399 register union double_long dl1;
403 if (!dl1.l.upper && !dl1.l.lower)
406 dl1.l.upper ^= SIGNBIT;
410 /* convert float to double */
412 __extendsfdf2 (float a1)
414 register union float_long fl1;
415 register union double_long dl;
422 dl.l.upper = dl.l.lower = 0;
426 dl.l.upper = SIGN (fl1.l);
427 exp = EXP (fl1.l) - EXCESS + EXCESSD;
428 dl.l.upper |= exp << 20;
429 dl.l.upper |= (MANT (fl1.l) & ~HIDDEN) >> 3;
430 dl.l.lower = MANT (fl1.l) << 29;
435 /* convert double to float */
437 __truncdfsf2 (double a1)
441 register union float_long fl;
442 register union double_long dl1;
446 if (!dl1.l.upper && !dl1.l.lower)
449 exp = EXPD (dl1) - EXCESSD + EXCESS;
451 /* shift double mantissa 6 bits so we can round */
452 mant = MANTD (dl1) >> 6;
454 /* now round and shift down */
458 /* did the round overflow? */
459 if (mant & 0xFF000000)
467 /* pack up and go home */
468 fl.l = PACK (SIGND (dl1), exp, mant);
472 /* compare two doubles */
474 __cmpdf2 (double a1, double a2)
476 register union double_long dl1, dl2;
481 if (SIGND (dl1) && SIGND (dl2))
483 dl1.l.upper ^= SIGNBIT;
484 dl2.l.upper ^= SIGNBIT;
486 if (dl1.l.upper < dl2.l.upper)
488 if (dl1.l.upper > dl2.l.upper)
490 if (dl1.l.lower < dl2.l.lower)
492 if (dl1.l.lower > dl2.l.lower)
497 /* convert double to int */
499 __fixdfsi (double a1)
501 register union double_long dl1;
507 if (!dl1.l.upper && !dl1.l.lower)
510 exp = EXPD (dl1) - EXCESSD - 31;
514 return (0x7FFFFFFF | SIGND (dl1)); /* largest integer */
516 /* shift down until exp = 0 or l = 0 */
517 if (exp < 0 && exp > -32 && l)
522 return (SIGND (dl1) ? -l : l);
525 /* convert double to unsigned int */
527 long __fixunsdfsi (double a1)
529 register union double_long dl1;
531 register unsigned long l;
535 if (!dl1.l.upper && !dl1.l.lower)
538 exp = EXPD (dl1) - EXCESSD - 32;
539 l = (((((dl1.l.upper) & 0xFFFFF) | HIDDEND) << 11) | (dl1.l.lower >> 21));
542 return (0xFFFFFFFF); /* largest integer */
544 /* shift down until exp = 0 or l = 0 */
545 if (exp < 0 && exp > -32 && l)
553 /* For now, the hard double-precision routines simply
554 punt and do it in single */
557 __adddf3 (double a1, double a2)
559 return ((float) a1 + (float) a2);
562 /* subtract two doubles */
564 __subdf3 (double a1, double a2)
566 return ((float) a1 - (float) a2);
569 /* multiply two doubles */
571 __muldf3 (double a1, double a2)
573 return ((float) a1 * (float) a2);
576 /* divide two doubles */
578 __divdf3 (double a1, double a2)
580 return ((float) a1 / (float) a2);