compiler-rt/lib/builtins/fp_mul_impl.inc

   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
  16 static __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 }