1//===---- lib/fp_mul_impl.inc - floating point multiplication -----*- C -*-===//
2//
3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4// See https://llvm.org/LICENSE.txt for license information.
5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6//
7//===----------------------------------------------------------------------===//
8//
9// This file implements soft-float multiplication with the IEEE-754 default
10// rounding (to nearest, ties to even).
11//
12//===----------------------------------------------------------------------===//
13
14#include "fp_lib.h"
15
16static __inline fp_t __mulXf3__(fp_t a, fp_t b) {
17 const unsigned int aExponent = toRep(x: a) >> significandBits & maxExponent;
18 const unsigned int bExponent = toRep(x: b) >> significandBits & maxExponent;
19 const rep_t productSign = (toRep(x: a) ^ toRep(x: b)) & signBit;
20
21 rep_t aSignificand = toRep(x: a) & significandMask;
22 rep_t bSignificand = toRep(x: b) & significandMask;
23 int scale = 0;
24
25 // Detect if a or b is zero, denormal, infinity, or NaN.
26 if (aExponent - 1U >= maxExponent - 1U ||
27 bExponent - 1U >= maxExponent - 1U) {
28
29 const rep_t aAbs = toRep(x: a) & absMask;
30 const rep_t bAbs = toRep(x: b) & absMask;
31
32 // NaN * anything = qNaN
33 if (aAbs > infRep)
34 return fromRep(x: toRep(x: a) | quietBit);
35 // anything * NaN = qNaN
36 if (bAbs > infRep)
37 return fromRep(x: toRep(x: b) | quietBit);
38
39 if (aAbs == infRep) {
40 // infinity * non-zero = +/- infinity
41 if (bAbs)
42 return fromRep(x: aAbs | productSign);
43 // infinity * zero = NaN
44 else
45 return fromRep(qnanRep);
46 }
47
48 if (bAbs == infRep) {
49 // non-zero * infinity = +/- infinity
50 if (aAbs)
51 return fromRep(x: bAbs | productSign);
52 // zero * infinity = NaN
53 else
54 return fromRep(qnanRep);
55 }
56
57 // zero * anything = +/- zero
58 if (!aAbs)
59 return fromRep(x: productSign);
60 // anything * zero = +/- zero
61 if (!bAbs)
62 return fromRep(x: productSign);
63
64 // One or both of a or b is denormal. The other (if applicable) is a
65 // normal number. Renormalize one or both of a and b, and set scale to
66 // include the necessary exponent adjustment.
67 if (aAbs < implicitBit)
68 scale += normalize(significand: &aSignificand);
69 if (bAbs < implicitBit)
70 scale += normalize(significand: &bSignificand);
71 }
72
73 // Set the implicit significand bit. If we fell through from the
74 // denormal path it was already set by normalize( ), but setting it twice
75 // won't hurt anything.
76 aSignificand |= implicitBit;
77 bSignificand |= implicitBit;
78
79 // Perform a basic multiplication on the significands. One of them must be
80 // shifted beforehand to be aligned with the exponent.
81 rep_t productHi, productLo;
82 wideMultiply(a: aSignificand, b: bSignificand << exponentBits, hi: &productHi,
83 lo: &productLo);
84
85 int productExponent = aExponent + bExponent - exponentBias + scale;
86
87 // Normalize the significand and adjust the exponent if needed.
88 if (productHi & implicitBit)
89 productExponent++;
90 else
91 wideLeftShift(hi: &productHi, lo: &productLo, count: 1);
92
93 // If we have overflowed the type, return +/- infinity.
94 if (productExponent >= maxExponent)
95 return fromRep(infRep | productSign);
96
97 if (productExponent <= 0) {
98 // The result is denormal before rounding.
99 //
100 // If the result is so small that it just underflows to zero, return
101 // zero with the appropriate sign. Mathematically, there is no need to
102 // handle this case separately, but we make it a special case to
103 // simplify the shift logic.
104 const unsigned int shift = REP_C(1) - (unsigned int)productExponent;
105 if (shift >= typeWidth)
106 return fromRep(x: productSign);
107
108 // Otherwise, shift the significand of the result so that the round
109 // bit is the high bit of productLo.
110 wideRightShiftWithSticky(hi: &productHi, lo: &productLo, count: shift);
111 } else {
112 // The result is normal before rounding. Insert the exponent.
113 productHi &= significandMask;
114 productHi |= (rep_t)productExponent << significandBits;
115 }
116
117 // Insert the sign of the result.
118 productHi |= productSign;
119
120 // Perform the final rounding. The final result may overflow to infinity,
121 // or underflow to zero, but those are the correct results in those cases.
122 // We use the default IEEE-754 round-to-nearest, ties-to-even rounding mode.
123 if (productLo > signBit)
124 productHi++;
125 if (productLo == signBit)
126 productHi += productHi & 1;
127 return fromRep(x: productHi);
128}
129

source code of compiler-rt/lib/builtins/fp_mul_impl.inc