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(a) >> significandBits & maxExponent;
18  const unsigned int bExponent = toRep(b) >> significandBits & maxExponent;
19  const rep_t productSign = (toRep(a) ^ toRep(b)) & signBit;
20
21  rep_t aSignificand = toRep(a) & significandMask;
22  rep_t bSignificand = toRep(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(a) & absMask;
30    const rep_t bAbs = toRep(b) & absMask;
31
32    // NaN * anything = qNaN
33    if (aAbs > infRep)
34      return fromRep(toRep(a) | quietBit);
35    // anything * NaN = qNaN
36    if (bAbs > infRep)
37      return fromRep(toRep(b) | quietBit);
38
39    if (aAbs == infRep) {
40      // infinity * non-zero = +/- infinity
41      if (bAbs)
42        return fromRep(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(bAbs | productSign);
52      // zero * infinity = NaN
53      else
54        return fromRep(qnanRep);
55    }
56
57    // zero * anything = +/- zero
58    if (!aAbs)
59      return fromRep(productSign);
60    // anything * zero = +/- zero
61    if (!bAbs)
62      return fromRep(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(&aSignificand);
69    if (bAbs < implicitBit)
70      scale += normalize(&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(aSignificand, bSignificand << exponentBits, &productHi,
83               &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(&productHi, &productLo, 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(productSign);
107
108    // Otherwise, shift the significand of the result so that the round
109    // bit is the high bit of productLo.
110    wideRightShiftWithSticky(&productHi, &productLo, 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(productHi);
128}
129