contrib/compiler-rt/lib/mulsf3.c

   1 //===-- lib/mulsf3.c - Single-precision 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 single-precision soft-float multiplication
  11 // with the IEEE-754 default rounding (to nearest, ties to even).
  12 //
  13 //===----------------------------------------------------------------------===//
  14 #include "abi.h"
  15
  16 #define SINGLE_PRECISION
  17 #include "fp_lib.h"
  18
  19 ARM_EABI_FNALIAS(fmul, mulsf3);
  20
  21 COMPILER_RT_ABI fp_t
  22 __mulsf3(fp_t a, fp_t b) {
  23
  24     const unsigned int aExponent = toRep(a) >> significandBits & maxExponent;
  25     const unsigned int bExponent = toRep(b) >> significandBits & maxExponent;
  26     const rep_t productSign = (toRep(a) ^ toRep(b)) & signBit;
  27
  28     rep_t aSignificand = toRep(a) & significandMask;
  29     rep_t bSignificand = toRep(b) & significandMask;
  30     int scale = 0;
  31
  32     // Detect if a or b is zero, denormal, infinity, or NaN.
  33     if (aExponent-1U >= maxExponent-1U || bExponent-1U >= maxExponent-1U) {
  34
  35         const rep_t aAbs = toRep(a) & absMask;
  36         const rep_t bAbs = toRep(b) & absMask;
  37
  38         // NaN * anything = qNaN
  39         if (aAbs > infRep) return fromRep(toRep(a) | quietBit);
  40         // anything * NaN = qNaN
  41         if (bAbs > infRep) return fromRep(toRep(b) | quietBit);
  42
  43         if (aAbs == infRep) {
  44             // infinity * non-zero = +/- infinity
  45             if (bAbs) return fromRep(aAbs | productSign);
  46             // infinity * zero = NaN
  47             else return fromRep(qnanRep);
  48         }
  49
  50         if (bAbs == infRep) {
  51             // non-zero * infinity = +/- infinity
  52             if (aAbs) return fromRep(bAbs | productSign);
  53             // zero * infinity = NaN
  54             else return fromRep(qnanRep);
  55         }
  56
  57         // zero * anything = +/- zero
  58         if (!aAbs) return fromRep(productSign);
  59         // anything * zero = +/- zero
  60         if (!bAbs) return fromRep(productSign);
  61
  62         // one or both of a or b is denormal, the other (if applicable) is a
  63         // normal number.  Renormalize one or both of a and b, and set scale to
  64         // include the necessary exponent adjustment.
  65         if (aAbs < implicitBit) scale += normalize(&aSignificand);
  66         if (bAbs < implicitBit) scale += normalize(&bSignificand);
  67     }
  68
  69     // Or in the implicit significand bit.  (If we fell through from the
  70     // denormal path it was already set by normalize( ), but setting it twice
  71     // won't hurt anything.)
  72     aSignificand |= implicitBit;
  73     bSignificand |= implicitBit;
  74
  75     // Get the significand of a*b.  Before multiplying the significands, shift
  76     // one of them left to left-align it in the field.  Thus, the product will
  77     // have (exponentBits + 2) integral digits, all but two of which must be
  78     // zero.  Normalizing this result is just a conditional left-shift by one
  79     // and bumping the exponent accordingly.
  80     rep_t productHi, productLo;
  81     wideMultiply(aSignificand, bSignificand << exponentBits,
  82                  &productHi, &productLo);
  83
  84     int productExponent = aExponent + bExponent - exponentBias + scale;
  85
  86     // Normalize the significand, adjust exponent if needed.
  87     if (productHi & implicitBit) productExponent++;
  88     else wideLeftShift(&productHi, &productLo, 1);
  89
  90     // If we have overflowed the type, return +/- infinity.
  91     if (productExponent >= maxExponent) return fromRep(infRep | productSign);
  92
  93     if (productExponent <= 0) {
  94         // Result is denormal before rounding, the exponent is zero and we
  95         // need to shift the significand.
  96         wideRightShiftWithSticky(&productHi, &productLo, 1 - productExponent);
  97     }
  98
  99     else {
 100         // Result is normal before rounding; insert the exponent.
 101         productHi &= significandMask;
 102         productHi |= (rep_t)productExponent << significandBits;
 103     }
 104
 105     // Insert the sign of the result:
 106     productHi |= productSign;
 107
 108     // Final rounding.  The final result may overflow to infinity, or underflow
 109     // to zero, but those are the correct results in those cases.
 110     if (productLo > signBit) productHi++;
 111     if (productLo == signBit) productHi += productHi & 1;
 112     return fromRep(productHi);
 113 }