326 lines · cpp
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-exception6//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 19using namespace llvm;20using namespace llvm::ScaledNumbers;21 22std::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 {Lower, 0};45 46 // Shift as little as possible to maximize precision.47 unsigned LeadingZeros = llvm::countl_zero(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 55static uint64_t getHalf(uint64_t N) { return (N >> 1) + (N & 1); }56 57std::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 = llvm::countl_zero(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 80std::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 = llvm::countr_zero(Divisor)) {88 Shift -= Zeros;89 Divisor >>= Zeros;90 }91 92 // Check for powers of two.93 if (Divisor == 1)94 return {Dividend, Shift};95 96 // Maximize size of dividend.97 if (int Zeros = llvm::countl_zero(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 124int 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 137static void appendDigit(std::string &Str, unsigned D) {138 assert(D < 10);139 Str += '0' + D % 10;140}141 142static void appendNumber(std::string &Str, uint64_t N) {143 while (N) {144 appendDigit(Str, N % 10);145 N /= 10;146 }147}148 149static 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 162static 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 191static 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 201std::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 }246 std::reverse(Str.begin(), Str.end());247 248 // Return early if there's nothing after the decimal.249 if (!Below0)250 return Str + ".0";251 252 // Append the decimal and beyond.253 Str += '.';254 uint64_t Error = UINT64_C(1) << (64 - Width);255 256 // We need to shift Below0 to the right to make space for calculating257 // digits. Save the precision we're losing in Extra.258 Extra = (Below0 & 0xf) << 56 | (Extra >> 8);259 Below0 >>= 4;260 size_t SinceDot = 0;261 size_t AfterDot = Str.size();262 do {263 if (ExtraShift) {264 --ExtraShift;265 Error *= 5;266 } else {267 Error *= 10;268 }269 270 Below0 *= 10;271 Extra *= 10;272 Below0 += (Extra >> 60);273 Extra = Extra & (UINT64_MAX >> 4);274 appendDigit(Str, Below0 >> 60);275 Below0 = Below0 & (UINT64_MAX >> 4);276 if (DigitsOut || Str.back() != '0')277 ++DigitsOut;278 ++SinceDot;279 } while (Error && (Below0 << 4 | Extra >> 60) >= Error / 2 &&280 (!Precision || DigitsOut <= Precision || SinceDot < 2));281 282 // Return early for maximum precision.283 if (!Precision || DigitsOut <= Precision)284 return stripTrailingZeros(Str);285 286 // Find where to truncate.287 size_t Truncate =288 std::max(Str.size() - (DigitsOut - Precision), AfterDot + 1);289 290 // Check if there's anything to truncate.291 if (Truncate >= Str.size())292 return stripTrailingZeros(Str);293 294 bool Carry = doesRoundUp(Str[Truncate]);295 if (!Carry)296 return stripTrailingZeros(Str.substr(0, Truncate));297 298 // Round with the first truncated digit.299 for (std::string::reverse_iterator I(Str.begin() + Truncate), E = Str.rend();300 I != E; ++I) {301 if (*I == '.')302 continue;303 if (*I == '9') {304 *I = '0';305 continue;306 }307 308 ++*I;309 Carry = false;310 break;311 }312 313 // Add "1" in front if we still need to carry.314 return stripTrailingZeros(std::string(Carry, '1') + Str.substr(0, Truncate));315}316 317raw_ostream &ScaledNumberBase::print(raw_ostream &OS, uint64_t D, int16_t E,318 int Width, unsigned Precision) {319 return OS << toString(D, E, Width, Precision);320}321 322void ScaledNumberBase::dump(uint64_t D, int16_t E, int Width) {323 print(dbgs(), D, E, Width, 0) << "[" << Width << ":" << D << "*2^" << E324 << "]";325}326