Home | History | Annotate | Download | only in Support 
      1 //==- lib/Support/ScaledNumber.cpp - Support for scaled numbers -*- C++ -*-===//
      2 //
      3 //                     The LLVM Compiler Infrastructure
      4 //
      5 // This file is distributed under the University of Illinois Open Source
      6 // License. See LICENSE.TXT for details.
      7 //
      8 //===----------------------------------------------------------------------===//
      9 //
     10 // Implementation of some scaled number algorithms.
     11 //
     12 //===----------------------------------------------------------------------===//
     13 
     14 #include "llvm/Support/ScaledNumber.h"
     15 #include "llvm/ADT/APFloat.h"
     16 #include "llvm/Support/Debug.h"
     17 #include "llvm/Support/raw_ostream.h"
     18 
     19 using namespace llvm;
     20 using namespace llvm::ScaledNumbers;
     21 
     22 std::pair<uint64_t, int16_t> ScaledNumbers::multiply64(uint64_t LHS,
     23                                                        uint64_t RHS) {
     24   // Separate into two 32-bit digits (U.L).
     25   auto getU = [](uint64_t N) { return N >> 32; };
     26   auto getL = [](uint64_t N) { return N & UINT32_MAX; };
     27   uint64_t UL = getU(LHS), LL = getL(LHS), UR = getU(RHS), LR = getL(RHS);
     28 
     29   // Compute cross products.
     30   uint64_t P1 = UL * UR, P2 = UL * LR, P3 = LL * UR, P4 = LL * LR;
     31 
     32   // Sum into two 64-bit digits.
     33   uint64_t Upper = P1, Lower = P4;
     34   auto addWithCarry = [&](uint64_t N) {
     35     uint64_t NewLower = Lower + (getL(N) << 32);
     36     Upper += getU(N) + (NewLower < Lower);
     37     Lower = NewLower;
     38   };
     39   addWithCarry(P2);
     40   addWithCarry(P3);
     41 
     42   // Check whether the upper digit is empty.
     43   if (!Upper)
     44     return std::make_pair(Lower, 0);
     45 
     46   // Shift as little as possible to maximize precision.
     47   unsigned LeadingZeros = countLeadingZeros(Upper);
     48   int Shift = 64 - LeadingZeros;
     49   if (LeadingZeros)
     50     Upper = Upper << LeadingZeros | Lower >> Shift;
     51   return getRounded(Upper, Shift,
     52                     Shift && (Lower & UINT64_C(1) << (Shift - 1)));
     53 }
     54 
     55 static uint64_t getHalf(uint64_t N) { return (N >> 1) + (N & 1); }
     56 
     57 std::pair<uint32_t, int16_t> ScaledNumbers::divide32(uint32_t Dividend,
     58                                                      uint32_t Divisor) {
     59   assert(Dividend && "expected non-zero dividend");
     60   assert(Divisor && "expected non-zero divisor");
     61 
     62   // Use 64-bit math and canonicalize the dividend to gain precision.
     63   uint64_t Dividend64 = Dividend;
     64   int Shift = 0;
     65   if (int Zeros = countLeadingZeros(Dividend64)) {
     66     Shift -= Zeros;
     67     Dividend64 <<= Zeros;
     68   }
     69   uint64_t Quotient = Dividend64 / Divisor;
     70   uint64_t Remainder = Dividend64 % Divisor;
     71 
     72   // If Quotient needs to be shifted, leave the rounding to getAdjusted().
     73   if (Quotient > UINT32_MAX)
     74     return getAdjusted<uint32_t>(Quotient, Shift);
     75 
     76   // Round based on the value of the next bit.
     77   return getRounded<uint32_t>(Quotient, Shift, Remainder >= getHalf(Divisor));
     78 }
     79 
     80 std::pair<uint64_t, int16_t> ScaledNumbers::divide64(uint64_t Dividend,
     81                                                      uint64_t Divisor) {
     82   assert(Dividend && "expected non-zero dividend");
     83   assert(Divisor && "expected non-zero divisor");
     84 
     85   // Minimize size of divisor.
     86   int Shift = 0;
     87   if (int Zeros = countTrailingZeros(Divisor)) {
     88     Shift -= Zeros;
     89     Divisor >>= Zeros;
     90   }
     91 
     92   // Check for powers of two.
     93   if (Divisor == 1)
     94     return std::make_pair(Dividend, Shift);
     95 
     96   // Maximize size of dividend.
     97   if (int Zeros = countLeadingZeros(Dividend)) {
     98     Shift -= Zeros;
     99     Dividend <<= Zeros;
    100   }
    101 
    102   // Start with the result of a divide.
    103   uint64_t Quotient = Dividend / Divisor;
    104   Dividend %= Divisor;
    105 
    106   // Continue building the quotient with long division.
    107   while (!(Quotient >> 63) && Dividend) {
    108     // Shift Dividend and check for overflow.
    109     bool IsOverflow = Dividend >> 63;
    110     Dividend <<= 1;
    111     --Shift;
    112 
    113     // Get the next bit of Quotient.
    114     Quotient <<= 1;
    115     if (IsOverflow || Divisor <= Dividend) {
    116       Quotient |= 1;
    117       Dividend -= Divisor;
    118     }
    119   }
    120 
    121   return getRounded(Quotient, Shift, Dividend >= getHalf(Divisor));
    122 }
    123 
    124 int ScaledNumbers::compareImpl(uint64_t L, uint64_t R, int ScaleDiff) {
    125   assert(ScaleDiff >= 0 && "wrong argument order");
    126   assert(ScaleDiff < 64 && "numbers too far apart");
    127 
    128   uint64_t L_adjusted = L >> ScaleDiff;
    129   if (L_adjusted < R)
    130     return -1;
    131   if (L_adjusted > R)
    132     return 1;
    133 
    134   return L > L_adjusted << ScaleDiff ? 1 : 0;
    135 }
    136 
    137 static void appendDigit(std::string &Str, unsigned D) {
    138   assert(D < 10);
    139   Str += '0' + D % 10;
    140 }
    141 
    142 static void appendNumber(std::string &Str, uint64_t N) {
    143   while (N) {
    144     appendDigit(Str, N % 10);
    145     N /= 10;
    146   }
    147 }
    148 
    149 static bool doesRoundUp(char Digit) {
    150   switch (Digit) {
    151   case '5':
    152   case '6':
    153   case '7':
    154   case '8':
    155   case '9':
    156     return true;
    157   default:
    158     return false;
    159   }
    160 }
    161 
    162 static std::string toStringAPFloat(uint64_t D, int E, unsigned Precision) {
    163   assert(E >= ScaledNumbers::MinScale);
    164   assert(E <= ScaledNumbers::MaxScale);
    165 
    166   // Find a new E, but don't let it increase past MaxScale.
    167   int LeadingZeros = ScaledNumberBase::countLeadingZeros64(D);
    168   int NewE = std::min(ScaledNumbers::MaxScale, E + 63 - LeadingZeros);
    169   int Shift = 63 - (NewE - E);
    170   assert(Shift <= LeadingZeros);
    171   assert(Shift == LeadingZeros || NewE == ScaledNumbers::MaxScale);
    172   assert(Shift >= 0 && Shift < 64 && "undefined behavior");
    173   D <<= Shift;
    174   E = NewE;
    175 
    176   // Check for a denormal.
    177   unsigned AdjustedE = E + 16383;
    178   if (!(D >> 63)) {
    179     assert(E == ScaledNumbers::MaxScale);
    180     AdjustedE = 0;
    181   }
    182 
    183   // Build the float and print it.
    184   uint64_t RawBits[2] = {D, AdjustedE};
    185   APFloat Float(APFloat::x87DoubleExtended, APInt(80, RawBits));
    186   SmallVector<char, 24> Chars;
    187   Float.toString(Chars, Precision, 0);
    188   return std::string(Chars.begin(), Chars.end());
    189 }
    190 
    191 static std::string stripTrailingZeros(const std::string &Float) {
    192   size_t NonZero = Float.find_last_not_of('0');
    193   assert(NonZero != std::string::npos && "no . in floating point string");
    194 
    195   if (Float[NonZero] == '.')
    196     ++NonZero;
    197 
    198   return Float.substr(0, NonZero + 1);
    199 }
    200 
    201 std::string ScaledNumberBase::toString(uint64_t D, int16_t E, int Width,
    202                                        unsigned Precision) {
    203   if (!D)
    204     return "0.0";
    205 
    206   // Canonicalize exponent and digits.
    207   uint64_t Above0 = 0;
    208   uint64_t Below0 = 0;
    209   uint64_t Extra = 0;
    210   int ExtraShift = 0;
    211   if (E == 0) {
    212     Above0 = D;
    213   } else if (E > 0) {
    214     if (int Shift = std::min(int16_t(countLeadingZeros64(D)), E)) {
    215       D <<= Shift;
    216       E -= Shift;
    217 
    218       if (!E)
    219         Above0 = D;
    220     }
    221   } else if (E > -64) {
    222     Above0 = D >> -E;
    223     Below0 = D << (64 + E);
    224   } else if (E == -64) {
    225     // Special case: shift by 64 bits is undefined behavior.
    226     Below0 = D;
    227   } else if (E > -120) {
    228     Below0 = D >> (-E - 64);
    229     Extra = D << (128 + E);
    230     ExtraShift = -64 - E;
    231   }
    232 
    233   // Fall back on APFloat for very small and very large numbers.
    234   if (!Above0 && !Below0)
    235     return toStringAPFloat(D, E, Precision);
    236 
    237   // Append the digits before the decimal.
    238   std::string Str;
    239   size_t DigitsOut = 0;
    240   if (Above0) {
    241     appendNumber(Str, Above0);
    242     DigitsOut = Str.size();
    243   } else
    244     appendDigit(Str, 0);
    245   std::reverse(Str.begin(), Str.end());
    246 
    247   // Return early if there's nothing after the decimal.
    248   if (!Below0)
    249     return Str + ".0";
    250 
    251   // Append the decimal and beyond.
    252   Str += '.';
    253   uint64_t Error = UINT64_C(1) << (64 - Width);
    254 
    255   // We need to shift Below0 to the right to make space for calculating
    256   // digits.  Save the precision we're losing in Extra.
    257   Extra = (Below0 & 0xf) << 56 | (Extra >> 8);
    258   Below0 >>= 4;
    259   size_t SinceDot = 0;
    260   size_t AfterDot = Str.size();
    261   do {
    262     if (ExtraShift) {
    263       --ExtraShift;
    264       Error *= 5;
    265     } else
    266       Error *= 10;
    267 
    268     Below0 *= 10;
    269     Extra *= 10;
    270     Below0 += (Extra >> 60);
    271     Extra = Extra & (UINT64_MAX >> 4);
    272     appendDigit(Str, Below0 >> 60);
    273     Below0 = Below0 & (UINT64_MAX >> 4);
    274     if (DigitsOut || Str.back() != '0')
    275       ++DigitsOut;
    276     ++SinceDot;
    277   } while (Error && (Below0 << 4 | Extra >> 60) >= Error / 2 &&
    278            (!Precision || DigitsOut <= Precision || SinceDot < 2));
    279 
    280   // Return early for maximum precision.
    281   if (!Precision || DigitsOut <= Precision)
    282     return stripTrailingZeros(Str);
    283 
    284   // Find where to truncate.
    285   size_t Truncate =
    286       std::max(Str.size() - (DigitsOut - Precision), AfterDot + 1);
    287 
    288   // Check if there's anything to truncate.
    289   if (Truncate >= Str.size())
    290     return stripTrailingZeros(Str);
    291 
    292   bool Carry = doesRoundUp(Str[Truncate]);
    293   if (!Carry)
    294     return stripTrailingZeros(Str.substr(0, Truncate));
    295 
    296   // Round with the first truncated digit.
    297   for (std::string::reverse_iterator I(Str.begin() + Truncate), E = Str.rend();
    298        I != E; ++I) {
    299     if (*I == '.')
    300       continue;
    301     if (*I == '9') {
    302       *I = '0';
    303       continue;
    304     }
    305 
    306     ++*I;
    307     Carry = false;
    308     break;
    309   }
    310 
    311   // Add "1" in front if we still need to carry.
    312   return stripTrailingZeros(std::string(Carry, '1') + Str.substr(0, Truncate));
    313 }
    314 
    315 raw_ostream &ScaledNumberBase::print(raw_ostream &OS, uint64_t D, int16_t E,
    316                                      int Width, unsigned Precision) {
    317   return OS << toString(D, E, Width, Precision);
    318 }
    319 
    320 void ScaledNumberBase::dump(uint64_t D, int16_t E, int Width) {
    321   print(dbgs(), D, E, Width, 0) << "[" << Width << ":" << D << "*2^" << E
    322                                 << "]";
    323 }
    324