lib/muldf3.c

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