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