kit

kit
git clone https://git.ryansepassi.com/git/kit.git
Log | Files | Refs | README

fp_mul_impl.inc (4916B)

      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 #define __mulXf3__ _FP_NAME(__mulXf3__)
     17 
     18 #if defined SINGLE_PRECISION && !defined FP_MUL_SF_EMITTED
     19 #define FP_MUL_SF_EMITTED
     20 #define _FP_MUL_EMIT 1
     21 #elif defined DOUBLE_PRECISION && !defined FP_MUL_DF_EMITTED
     22 #define FP_MUL_DF_EMITTED
     23 #define _FP_MUL_EMIT 1
     24 #elif defined QUAD_PRECISION && !defined FP_MUL_TF_EMITTED
     25 #define FP_MUL_TF_EMITTED
     26 #define _FP_MUL_EMIT 1
     27 #endif
     28 
     29 #ifdef _FP_MUL_EMIT
     30 #undef _FP_MUL_EMIT
     31 
     32 static inline fp_t __mulXf3__(fp_t a, fp_t b) {
     33   const unsigned int aExponent = toRep(a) >> significandBits & maxExponent;
     34   const unsigned int bExponent = toRep(b) >> significandBits & maxExponent;
     35   const rep_t productSign = (toRep(a) ^ toRep(b)) & signBit;
     36 
     37   rep_t aSignificand = toRep(a) & significandMask;
     38   rep_t bSignificand = toRep(b) & significandMask;
     39   int scale = 0;
     40 
     41   // Detect if a or b is zero, denormal, infinity, or NaN.
     42   if (aExponent - 1U >= maxExponent - 1U ||
     43       bExponent - 1U >= maxExponent - 1U) {
     44 
     45     const rep_t aAbs = toRep(a) & absMask;
     46     const rep_t bAbs = toRep(b) & absMask;
     47 
     48     // NaN * anything = qNaN
     49     if (aAbs > infRep)
     50       return fromRep(toRep(a) | quietBit);
     51     // anything * NaN = qNaN
     52     if (bAbs > infRep)
     53       return fromRep(toRep(b) | quietBit);
     54 
     55     if (aAbs == infRep) {
     56       // infinity * non-zero = +/- infinity
     57       if (bAbs)
     58         return fromRep(aAbs | productSign);
     59       // infinity * zero = NaN
     60       else
     61         return fromRep(qnanRep);
     62     }
     63 
     64     if (bAbs == infRep) {
     65       // non-zero * infinity = +/- infinity
     66       if (aAbs)
     67         return fromRep(bAbs | productSign);
     68       // zero * infinity = NaN
     69       else
     70         return fromRep(qnanRep);
     71     }
     72 
     73     // zero * anything = +/- zero
     74     if (!aAbs)
     75       return fromRep(productSign);
     76     // anything * zero = +/- zero
     77     if (!bAbs)
     78       return fromRep(productSign);
     79 
     80     // One or both of a or b is denormal.  The other (if applicable) is a
     81     // normal number.  Renormalize one or both of a and b, and set scale to
     82     // include the necessary exponent adjustment.
     83     if (aAbs < implicitBit)
     84       scale += normalize(&aSignificand);
     85     if (bAbs < implicitBit)
     86       scale += normalize(&bSignificand);
     87   }
     88 
     89   // Set the implicit significand bit.  If we fell through from the
     90   // denormal path it was already set by normalize( ), but setting it twice
     91   // won't hurt anything.
     92   aSignificand |= implicitBit;
     93   bSignificand |= implicitBit;
     94 
     95   // Perform a basic multiplication on the significands.  One of them must be
     96   // shifted beforehand to be aligned with the exponent.
     97   rep_t productHi, productLo;
     98   wideMultiply(aSignificand, bSignificand << exponentBits, &productHi,
     99                &productLo);
    100 
    101   int productExponent = aExponent + bExponent - exponentBias + scale;
    102 
    103   // Normalize the significand and adjust the exponent if needed.
    104   if (productHi & implicitBit)
    105     productExponent++;
    106   else
    107     wideLeftShift(&productHi, &productLo, 1);
    108 
    109   // If we have overflowed the type, return +/- infinity.
    110   if (productExponent >= maxExponent)
    111     return fromRep(infRep | productSign);
    112 
    113   if (productExponent <= 0) {
    114     // The result is denormal before rounding.
    115     //
    116     // If the result is so small that it just underflows to zero, return
    117     // zero with the appropriate sign.  Mathematically, there is no need to
    118     // handle this case separately, but we make it a special case to
    119     // simplify the shift logic.
    120     const unsigned int shift = REP_C(1) - (unsigned int)productExponent;
    121     if (shift >= typeWidth)
    122       return fromRep(productSign);
    123 
    124     // Otherwise, shift the significand of the result so that the round
    125     // bit is the high bit of productLo.
    126     wideRightShiftWithSticky(&productHi, &productLo, shift);
    127   } else {
    128     // The result is normal before rounding.  Insert the exponent.
    129     productHi &= significandMask;
    130     productHi |= (rep_t)productExponent << significandBits;
    131   }
    132 
    133   // Insert the sign of the result.
    134   productHi |= productSign;
    135 
    136   // Perform the final rounding.  The final result may overflow to infinity,
    137   // or underflow to zero, but those are the correct results in those cases.
    138   // We use the default IEEE-754 round-to-nearest, ties-to-even rounding mode.
    139   if (productLo > signBit)
    140     productHi++;
    141   if (productLo == signBit)
    142     productHi += productHi & 1;
    143   return fromRep(productHi);
    144 }
    145 
    146 #endif // _FP_MUL_EMIT