diff options
Diffstat (limited to 'compiler-rt/lib/builtins/fp_mul_impl.inc')
-rw-r--r-- | compiler-rt/lib/builtins/fp_mul_impl.inc | 128 |
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); +} |