1 files changed, 128 insertions, 0 deletions
diff --git a/compiler-rt/lib/builtins/fp_mul_impl.inc b/compiler-rt/lib/builtins/fp_mul_impl.inc
new file mode 100644
index 000000000000..a93f2d78ad6b
--- /dev/null
+++ b/compiler-rt/lib/builtins/fp_mul_impl.inc
@@ -0,0 +1,128 @@
+//===---- lib/fp_mul_impl.inc - floating point multiplication -----*- C -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+//
+// This file implements soft-float multiplication with the IEEE-754 default
+// rounding (to nearest, ties to even).
+//
+//===----------------------------------------------------------------------===//
+
+#include "fp_lib.h"
+
+static __inline fp_t __mulXf3__(fp_t a, fp_t b) {
+  const unsigned int aExponent = toRep(a) >> significandBits & maxExponent;
+  const unsigned int bExponent = toRep(b) >> significandBits & maxExponent;
+  const rep_t productSign = (toRep(a) ^ toRep(b)) & signBit;
+
+  rep_t aSignificand = toRep(a) & significandMask;
+  rep_t bSignificand = toRep(b) & significandMask;
+  int scale = 0;
+
+  // Detect if a or b is zero, denormal, infinity, or NaN.
+  if (aExponent - 1U >= maxExponent - 1U ||
+      bExponent - 1U >= maxExponent - 1U) {
+
+    const rep_t aAbs = toRep(a) & absMask;
+    const rep_t bAbs = toRep(b) & absMask;
+
+    // NaN * anything = qNaN
+    if (aAbs > infRep)
+      return fromRep(toRep(a) | quietBit);
+    // anything * NaN = qNaN
+    if (bAbs > infRep)
+      return fromRep(toRep(b) | quietBit);
+
+    if (aAbs == infRep) {
+      // infinity * non-zero = +/- infinity
+      if (bAbs)
+        return fromRep(aAbs | productSign);
+      // infinity * zero = NaN
+      else
+        return fromRep(qnanRep);
+    }
+
+    if (bAbs == infRep) {
+      // non-zero * infinity = +/- infinity
+      if (aAbs)
+        return fromRep(bAbs | productSign);
+      // zero * infinity = NaN
+      else
+        return fromRep(qnanRep);
+    }
+
+    // zero * anything = +/- zero
+    if (!aAbs)
+      return fromRep(productSign);
+    // anything * zero = +/- zero
+    if (!bAbs)
+      return fromRep(productSign);
+
+    // One or both of a or b is denormal.  The other (if applicable) is a
+    // normal number.  Renormalize one or both of a and b, and set scale to
+    // include the necessary exponent adjustment.
+    if (aAbs < implicitBit)
+      scale += normalize(&aSignificand);
+    if (bAbs < implicitBit)
+      scale += normalize(&bSignificand);
+  }
+
+  // Set the implicit significand bit.  If we fell through from the
+  // denormal path it was already set by normalize( ), but setting it twice
+  // won't hurt anything.
+  aSignificand |= implicitBit;
+  bSignificand |= implicitBit;
+
+  // Perform a basic multiplication on the significands.  One of them must be
+  // shifted beforehand to be aligned with the exponent.
+  rep_t productHi, productLo;
+  wideMultiply(aSignificand, bSignificand << exponentBits, &productHi,
+               &productLo);
+
+  int productExponent = aExponent + bExponent - exponentBias + scale;
+
+  // Normalize the significand and adjust the exponent if needed.
+  if (productHi & implicitBit)
+    productExponent++;
+  else
+    wideLeftShift(&productHi, &productLo, 1);
+
+  // If we have overflowed the type, return +/- infinity.
+  if (productExponent >= maxExponent)
+    return fromRep(infRep | productSign);
+
+  if (productExponent <= 0) {
+    // The result is denormal before rounding.
+    //
+    // If the result is so small that it just underflows to zero, return
+    // zero with the appropriate sign.  Mathematically, there is no need to
+    // handle this case separately, but we make it a special case to
+    // simplify the shift logic.
+    const unsigned int shift = REP_C(1) - (unsigned int)productExponent;
+    if (shift >= typeWidth)
+      return fromRep(productSign);
+
+    // Otherwise, shift the significand of the result so that the round
+    // bit is the high bit of productLo.
+    wideRightShiftWithSticky(&productHi, &productLo, shift);
+  } else {
+    // The result is normal before rounding.  Insert the exponent.
+    productHi &= significandMask;
+    productHi |= (rep_t)productExponent << significandBits;
+  }
+
+  // Insert the sign of the result.
+  productHi |= productSign;
+
+  // Perform the final rounding.  The final result may overflow to infinity,
+  // or underflow to zero, but those are the correct results in those cases.
+  // We use the default IEEE-754 round-to-nearest, ties-to-even rounding mode.
+  if (productLo > signBit)
+    productHi++;
+  if (productLo == signBit)
+    productHi += productHi & 1;
+  return fromRep(productHi);
+}