]>
Commit | Line | Data |
---|---|---|
92a42be0 SL |
1 | //===---- lib/fp_mul_impl.inc - floating point multiplication -----*- C -*-===// |
2 | // | |
3 | // The LLVM Compiler Infrastructure | |
4 | // | |
5 | // This file is dual licensed under the MIT and the University of Illinois Open | |
6 | // Source Licenses. See LICENSE.TXT for details. | |
7 | // | |
8 | //===----------------------------------------------------------------------===// | |
9 | // | |
10 | // This file implements soft-float multiplication with the IEEE-754 default | |
11 | // rounding (to nearest, ties to even). | |
12 | // | |
13 | //===----------------------------------------------------------------------===// | |
14 | ||
15 | #include "fp_lib.h" | |
16 | ||
17 | static __inline fp_t __mulXf3__(fp_t a, fp_t b) { | |
18 | const unsigned int aExponent = toRep(a) >> significandBits & maxExponent; | |
19 | const unsigned int bExponent = toRep(b) >> significandBits & maxExponent; | |
20 | const rep_t productSign = (toRep(a) ^ toRep(b)) & signBit; | |
21 | ||
22 | rep_t aSignificand = toRep(a) & significandMask; | |
23 | rep_t bSignificand = toRep(b) & significandMask; | |
24 | int scale = 0; | |
25 | ||
26 | // Detect if a or b is zero, denormal, infinity, or NaN. | |
27 | if (aExponent-1U >= maxExponent-1U || bExponent-1U >= maxExponent-1U) { | |
28 | ||
29 | const rep_t aAbs = toRep(a) & absMask; | |
30 | const rep_t bAbs = toRep(b) & absMask; | |
31 | ||
32 | // NaN * anything = qNaN | |
33 | if (aAbs > infRep) return fromRep(toRep(a) | quietBit); | |
34 | // anything * NaN = qNaN | |
35 | if (bAbs > infRep) return fromRep(toRep(b) | quietBit); | |
36 | ||
37 | if (aAbs == infRep) { | |
38 | // infinity * non-zero = +/- infinity | |
39 | if (bAbs) return fromRep(aAbs | productSign); | |
40 | // infinity * zero = NaN | |
41 | else return fromRep(qnanRep); | |
42 | } | |
43 | ||
44 | if (bAbs == infRep) { | |
45 | //? non-zero * infinity = +/- infinity | |
46 | if (aAbs) return fromRep(bAbs | productSign); | |
47 | // zero * infinity = NaN | |
48 | else return fromRep(qnanRep); | |
49 | } | |
50 | ||
51 | // zero * anything = +/- zero | |
52 | if (!aAbs) return fromRep(productSign); | |
53 | // anything * zero = +/- zero | |
54 | if (!bAbs) return fromRep(productSign); | |
55 | ||
56 | // one or both of a or b is denormal, the other (if applicable) is a | |
57 | // normal number. Renormalize one or both of a and b, and set scale to | |
58 | // include the necessary exponent adjustment. | |
59 | if (aAbs < implicitBit) scale += normalize(&aSignificand); | |
60 | if (bAbs < implicitBit) scale += normalize(&bSignificand); | |
61 | } | |
62 | ||
63 | // Or in the implicit significand bit. (If we fell through from the | |
64 | // denormal path it was already set by normalize( ), but setting it twice | |
65 | // won't hurt anything.) | |
66 | aSignificand |= implicitBit; | |
67 | bSignificand |= implicitBit; | |
68 | ||
69 | // Get the significand of a*b. Before multiplying the significands, shift | |
70 | // one of them left to left-align it in the field. Thus, the product will | |
71 | // have (exponentBits + 2) integral digits, all but two of which must be | |
72 | // zero. Normalizing this result is just a conditional left-shift by one | |
73 | // and bumping the exponent accordingly. | |
74 | rep_t productHi, productLo; | |
75 | wideMultiply(aSignificand, bSignificand << exponentBits, | |
76 | &productHi, &productLo); | |
77 | ||
78 | int productExponent = aExponent + bExponent - exponentBias + scale; | |
79 | ||
80 | // Normalize the significand, adjust exponent if needed. | |
81 | if (productHi & implicitBit) productExponent++; | |
82 | else wideLeftShift(&productHi, &productLo, 1); | |
83 | ||
84 | // If we have overflowed the type, return +/- infinity. | |
85 | if (productExponent >= maxExponent) return fromRep(infRep | productSign); | |
86 | ||
87 | if (productExponent <= 0) { | |
88 | // Result is denormal before rounding | |
89 | // | |
90 | // If the result is so small that it just underflows to zero, return | |
91 | // a zero of the appropriate sign. Mathematically there is no need to | |
92 | // handle this case separately, but we make it a special case to | |
93 | // simplify the shift logic. | |
94 | const unsigned int shift = REP_C(1) - (unsigned int)productExponent; | |
95 | if (shift >= typeWidth) return fromRep(productSign); | |
96 | ||
97 | // Otherwise, shift the significand of the result so that the round | |
98 | // bit is the high bit of productLo. | |
99 | wideRightShiftWithSticky(&productHi, &productLo, shift); | |
100 | } | |
101 | else { | |
102 | // Result is normal before rounding; insert the exponent. | |
103 | productHi &= significandMask; | |
104 | productHi |= (rep_t)productExponent << significandBits; | |
105 | } | |
106 | ||
107 | // Insert the sign of the result: | |
108 | productHi |= productSign; | |
109 | ||
110 | // Final rounding. The final result may overflow to infinity, or underflow | |
111 | // to zero, but those are the correct results in those cases. We use the | |
112 | // default IEEE-754 round-to-nearest, ties-to-even rounding mode. | |
113 | if (productLo > signBit) productHi++; | |
114 | if (productLo == signBit) productHi += productHi & 1; | |
115 | return fromRep(productHi); | |
116 | } |