1 //==- lib/Support/ScaledNumber.cpp - Support for scaled numbers -*- 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 // Implementation of some scaled number algorithms. 10 // 11 //===----------------------------------------------------------------------===// 12 13 #include "llvm/Support/ScaledNumber.h" 14 #include "llvm/ADT/APFloat.h" 15 #include "llvm/ADT/ArrayRef.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