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 |
|
|
} |