• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 //===-- lib/mulsf3.c - Single-precision 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 single-precision soft-float multiplication
11 // with the IEEE-754 default rounding (to nearest, ties to even).
12 //
13 //===----------------------------------------------------------------------===//
14 
15 #define SINGLE_PRECISION
16 #include "fp_lib.h"
17 
ARM_EABI_FNALIAS(fmul,mulsf3)18 ARM_EABI_FNALIAS(fmul, mulsf3)
19 
20 COMPILER_RT_ABI fp_t
21 __mulsf3(fp_t a, fp_t b) {
22 
23     const unsigned int aExponent = toRep(a) >> significandBits & maxExponent;
24     const unsigned int bExponent = toRep(b) >> significandBits & maxExponent;
25     const rep_t productSign = (toRep(a) ^ toRep(b)) & signBit;
26 
27     rep_t aSignificand = toRep(a) & significandMask;
28     rep_t bSignificand = toRep(b) & significandMask;
29     int scale = 0;
30 
31     // Detect if a or b is zero, denormal, infinity, or NaN.
32     if (aExponent-1U >= maxExponent-1U || bExponent-1U >= maxExponent-1U) {
33 
34         const rep_t aAbs = toRep(a) & absMask;
35         const rep_t bAbs = toRep(b) & absMask;
36 
37         // NaN * anything = qNaN
38         if (aAbs > infRep) return fromRep(toRep(a) | quietBit);
39         // anything * NaN = qNaN
40         if (bAbs > infRep) return fromRep(toRep(b) | quietBit);
41 
42         if (aAbs == infRep) {
43             // infinity * non-zero = +/- infinity
44             if (bAbs) return fromRep(aAbs | productSign);
45             // infinity * zero = NaN
46             else return fromRep(qnanRep);
47         }
48 
49         if (bAbs == infRep) {
50             // non-zero * infinity = +/- infinity
51             if (aAbs) return fromRep(bAbs | productSign);
52             // zero * infinity = NaN
53             else return fromRep(qnanRep);
54         }
55 
56         // zero * anything = +/- zero
57         if (!aAbs) return fromRep(productSign);
58         // anything * zero = +/- zero
59         if (!bAbs) return fromRep(productSign);
60 
61         // one or both of a or b is denormal, the other (if applicable) is a
62         // normal number.  Renormalize one or both of a and b, and set scale to
63         // include the necessary exponent adjustment.
64         if (aAbs < implicitBit) scale += normalize(&aSignificand);
65         if (bAbs < implicitBit) scale += normalize(&bSignificand);
66     }
67 
68     // Or in the implicit significand bit.  (If we fell through from the
69     // denormal path it was already set by normalize( ), but setting it twice
70     // won't hurt anything.)
71     aSignificand |= implicitBit;
72     bSignificand |= implicitBit;
73 
74     // Get the significand of a*b.  Before multiplying the significands, shift
75     // one of them left to left-align it in the field.  Thus, the product will
76     // have (exponentBits + 2) integral digits, all but two of which must be
77     // zero.  Normalizing this result is just a conditional left-shift by one
78     // and bumping the exponent accordingly.
79     rep_t productHi, productLo;
80     wideMultiply(aSignificand, bSignificand << exponentBits,
81                  &productHi, &productLo);
82 
83     int productExponent = aExponent + bExponent - exponentBias + scale;
84 
85     // Normalize the significand, adjust exponent if needed.
86     if (productHi & implicitBit) productExponent++;
87     else wideLeftShift(&productHi, &productLo, 1);
88 
89     // If we have overflowed the type, return +/- infinity.
90     if (productExponent >= maxExponent) return fromRep(infRep | productSign);
91 
92     if (productExponent <= 0) {
93         // Result is denormal before rounding, the exponent is zero and we
94         // need to shift the significand.
95         wideRightShiftWithSticky(&productHi, &productLo, 1U - (unsigned)productExponent);
96     }
97 
98     else {
99         // Result is normal before rounding; insert the exponent.
100         productHi &= significandMask;
101         productHi |= (rep_t)productExponent << significandBits;
102     }
103 
104     // Insert the sign of the result:
105     productHi |= productSign;
106 
107     // Final rounding.  The final result may overflow to infinity, or underflow
108     // to zero, but those are the correct results in those cases.
109     if (productLo > signBit) productHi++;
110     if (productLo == signBit) productHi += productHi & 1;
111     return fromRep(productHi);
112 }
113