lib/builtins/fp_mul_impl.inc

   1 //===---- lib/fp_mul_impl.inc - floating point multiplication -----*- C -*-===//
   2 //
   3 //                     The LLVM Compiler Infrastructure
   4 //
   5 // This file is dual licensed under the MIT and the University of Illinois Open
   6 // Source Licenses. See LICENSE.TXT for details.
   7 //
   8 //===----------------------------------------------------------------------===//
   9 //
  10 // This file implements soft-float multiplication with the IEEE-754 default
  11 // rounding (to nearest, ties to even).
  12 //
  13 //===----------------------------------------------------------------------===//
  14
  15 #include "fp_lib.h"
  16
  17 static inline fp_t __mulXf3__(fp_t a, fp_t b) {
  18     const unsigned int aExponent = toRep(a) >> significandBits & maxExponent;
  19     const unsigned int bExponent = toRep(b) >> significandBits & maxExponent;
  20     const rep_t productSign = (toRep(a) ^ toRep(b)) & signBit;
  21
  22     rep_t aSignificand = toRep(a) & significandMask;
  23     rep_t bSignificand = toRep(b) & significandMask;
  24     int scale = 0;
  25
  26     // Detect if a or b is zero, denormal, infinity, or NaN.
  27     if (aExponent-1U >= maxExponent-1U || 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) return fromRep(toRep(a) | quietBit);
  34         // anything * NaN = qNaN
  35         if (bAbs > infRep) return fromRep(toRep(b) | quietBit);
  36
  37         if (aAbs == infRep) {
  38             // infinity * non-zero = +/- infinity
  39             if (bAbs) return fromRep(aAbs | productSign);
  40             // infinity * zero = NaN
  41             else return fromRep(qnanRep);
  42         }
  43
  44         if (bAbs == infRep) {
  45             //? non-zero * infinity = +/- infinity
  46             if (aAbs) return fromRep(bAbs | productSign);
  47             // zero * infinity = NaN
  48             else return fromRep(qnanRep);
  49         }
  50
  51         // zero * anything = +/- zero
  52         if (!aAbs) return fromRep(productSign);
  53         // anything * zero = +/- zero
  54         if (!bAbs) return fromRep(productSign);
  55
  56         // one or both of a or b is denormal, the other (if applicable) is a
  57         // normal number.  Renormalize one or both of a and b, and set scale to
  58         // include the necessary exponent adjustment.
  59         if (aAbs < implicitBit) scale += normalize(&aSignificand);
  60         if (bAbs < implicitBit) scale += normalize(&bSignificand);
  61     }
  62
  63     // Or in the implicit significand bit.  (If we fell through from the
  64     // denormal path it was already set by normalize( ), but setting it twice
  65     // won't hurt anything.)
  66     aSignificand |= implicitBit;
  67     bSignificand |= implicitBit;
  68
  69     // Get the significand of a*b.  Before multiplying the significands, shift
  70     // one of them left to left-align it in the field.  Thus, the product will
  71     // have (exponentBits + 2) integral digits, all but two of which must be
  72     // zero.  Normalizing this result is just a conditional left-shift by one
  73     // and bumping the exponent accordingly.
  74     rep_t productHi, productLo;
  75     wideMultiply(aSignificand, bSignificand << exponentBits,
  76                  &productHi, &productLo);
  77
  78     int productExponent = aExponent + bExponent - exponentBias + scale;
  79
  80     // Normalize the significand, adjust exponent if needed.
  81     if (productHi & implicitBit) productExponent++;
  82     else wideLeftShift(&productHi, &productLo, 1);
  83
  84     // If we have overflowed the type, return +/- infinity.
  85     if (productExponent >= maxExponent) return fromRep(infRep | productSign);
  86
  87     if (productExponent <= 0) {
  88         // Result is denormal before rounding
  89         //
  90         // If the result is so small that it just underflows to zero, return
  91         // a zero of the appropriate sign.  Mathematically there is no need to
  92         // handle this case separately, but we make it a special case to
  93         // simplify the shift logic.
  94         const unsigned int shift = REP_C(1) - (unsigned int)productExponent;
  95         if (shift >= typeWidth) return fromRep(productSign);
  96
  97         // Otherwise, shift the significand of the result so that the round
  98         // bit is the high bit of productLo.
  99         wideRightShiftWithSticky(&productHi, &productLo, shift);
 100     }
 101     else {
 102         // Result is normal before rounding; insert the exponent.
 103         productHi &= significandMask;
 104         productHi |= (rep_t)productExponent << significandBits;
 105     }
 106
 107     // Insert the sign of the result:
 108     productHi |= productSign;
 109
 110     // Final rounding.  The final result may overflow to infinity, or underflow
 111     // to zero, but those are the correct results in those cases.  We use the
 112     // default IEEE-754 round-to-nearest, ties-to-even rounding mode.
 113     if (productLo > signBit) productHi++;
 114     if (productLo == signBit) productHi += productHi & 1;
 115     return fromRep(productHi);
 116 }