1//===-- High Precision Decimal ----------------------------------*- 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// -----------------------------------------------------------------------------
10// **** WARNING ****
11// This file is shared with libc++. You should also be careful when adding
12// dependencies to this file, since it needs to build for all libc++ targets.
13// -----------------------------------------------------------------------------
14
15#ifndef LLVM_LIBC_SRC___SUPPORT_HIGH_PRECISION_DECIMAL_H
16#define LLVM_LIBC_SRC___SUPPORT_HIGH_PRECISION_DECIMAL_H
17
18#include "hdr/stdint_proxy.h"
19#include "src/__support/CPP/limits.h"
20#include "src/__support/ctype_utils.h"
21#include "src/__support/macros/config.h"
22#include "src/__support/str_to_integer.h"
23#include "src/__support/wctype_utils.h"
24
25namespace LIBC_NAMESPACE_DECL {
26namespace internal {
27
28struct LShiftTableEntry {
29 uint32_t new_digits;
30 char const *power_of_five;
31};
32
33// -----------------------------------------------------------------------------
34// **** WARNING ****
35// This interface is shared with libc++, if you change this interface you need
36// to update it in both libc and libc++.
37// -----------------------------------------------------------------------------
38// This is used in both this file and in the main str_to_float.h.
39// TODO: Figure out where to put this.
40enum class RoundDirection { Up, Down, Nearest };
41
42// These constants are used in both this file and in the main str_to_float.h.
43// TODO: Figure out where to put this.
44template <typename CharType> struct constants;
45template <> struct constants<char> {
46 static constexpr char DECIMAL_POINT = '.';
47 static constexpr char DECIMAL_EXPONENT_MARKER = 'e';
48 static constexpr char HEX_EXPONENT_MARKER = 'p';
49 static constexpr char INF_STRING[] = "infinity";
50 static constexpr char NAN_STRING[] = "nan";
51};
52template <> struct constants<wchar_t> {
53 static constexpr wchar_t DECIMAL_POINT = L'.';
54 static constexpr wchar_t DECIMAL_EXPONENT_MARKER = L'e';
55 static constexpr wchar_t HEX_EXPONENT_MARKER = L'p';
56 static constexpr wchar_t INF_STRING[] = L"infinity";
57 static constexpr wchar_t NAN_STRING[] = L"nan";
58};
59
60// This is based on the HPD data structure described as part of the Simple
61// Decimal Conversion algorithm by Nigel Tao, described at this link:
62// https://nigeltao.github.io/blog/2020/parse-number-f64-simple.html
63class HighPrecisionDecimal {
64
65 // This precomputed table speeds up left shifts by having the number of new
66 // digits that will be added by multiplying 5^i by 2^i. If the number is less
67 // than 5^i then it will add one fewer digit. There are only 60 entries since
68 // that's the max shift amount.
69 // This table was generated by the script at
70 // libc/utils/mathtools/GenerateHPDConstants.py
71 static constexpr LShiftTableEntry LEFT_SHIFT_DIGIT_TABLE[] = {
72 {.new_digits: 0, .power_of_five: ""},
73 {.new_digits: 1, .power_of_five: "5"},
74 {.new_digits: 1, .power_of_five: "25"},
75 {.new_digits: 1, .power_of_five: "125"},
76 {.new_digits: 2, .power_of_five: "625"},
77 {.new_digits: 2, .power_of_five: "3125"},
78 {.new_digits: 2, .power_of_five: "15625"},
79 {.new_digits: 3, .power_of_five: "78125"},
80 {.new_digits: 3, .power_of_five: "390625"},
81 {.new_digits: 3, .power_of_five: "1953125"},
82 {.new_digits: 4, .power_of_five: "9765625"},
83 {.new_digits: 4, .power_of_five: "48828125"},
84 {.new_digits: 4, .power_of_five: "244140625"},
85 {.new_digits: 4, .power_of_five: "1220703125"},
86 {.new_digits: 5, .power_of_five: "6103515625"},
87 {.new_digits: 5, .power_of_five: "30517578125"},
88 {.new_digits: 5, .power_of_five: "152587890625"},
89 {.new_digits: 6, .power_of_five: "762939453125"},
90 {.new_digits: 6, .power_of_five: "3814697265625"},
91 {.new_digits: 6, .power_of_five: "19073486328125"},
92 {.new_digits: 7, .power_of_five: "95367431640625"},
93 {.new_digits: 7, .power_of_five: "476837158203125"},
94 {.new_digits: 7, .power_of_five: "2384185791015625"},
95 {.new_digits: 7, .power_of_five: "11920928955078125"},
96 {.new_digits: 8, .power_of_five: "59604644775390625"},
97 {.new_digits: 8, .power_of_five: "298023223876953125"},
98 {.new_digits: 8, .power_of_five: "1490116119384765625"},
99 {.new_digits: 9, .power_of_five: "7450580596923828125"},
100 {.new_digits: 9, .power_of_five: "37252902984619140625"},
101 {.new_digits: 9, .power_of_five: "186264514923095703125"},
102 {.new_digits: 10, .power_of_five: "931322574615478515625"},
103 {.new_digits: 10, .power_of_five: "4656612873077392578125"},
104 {.new_digits: 10, .power_of_five: "23283064365386962890625"},
105 {.new_digits: 10, .power_of_five: "116415321826934814453125"},
106 {.new_digits: 11, .power_of_five: "582076609134674072265625"},
107 {.new_digits: 11, .power_of_five: "2910383045673370361328125"},
108 {.new_digits: 11, .power_of_five: "14551915228366851806640625"},
109 {.new_digits: 12, .power_of_five: "72759576141834259033203125"},
110 {.new_digits: 12, .power_of_five: "363797880709171295166015625"},
111 {.new_digits: 12, .power_of_five: "1818989403545856475830078125"},
112 {.new_digits: 13, .power_of_five: "9094947017729282379150390625"},
113 {.new_digits: 13, .power_of_five: "45474735088646411895751953125"},
114 {.new_digits: 13, .power_of_five: "227373675443232059478759765625"},
115 {.new_digits: 13, .power_of_five: "1136868377216160297393798828125"},
116 {.new_digits: 14, .power_of_five: "5684341886080801486968994140625"},
117 {.new_digits: 14, .power_of_five: "28421709430404007434844970703125"},
118 {.new_digits: 14, .power_of_five: "142108547152020037174224853515625"},
119 {.new_digits: 15, .power_of_five: "710542735760100185871124267578125"},
120 {.new_digits: 15, .power_of_five: "3552713678800500929355621337890625"},
121 {.new_digits: 15, .power_of_five: "17763568394002504646778106689453125"},
122 {.new_digits: 16, .power_of_five: "88817841970012523233890533447265625"},
123 {.new_digits: 16, .power_of_five: "444089209850062616169452667236328125"},
124 {.new_digits: 16, .power_of_five: "2220446049250313080847263336181640625"},
125 {.new_digits: 16, .power_of_five: "11102230246251565404236316680908203125"},
126 {.new_digits: 17, .power_of_five: "55511151231257827021181583404541015625"},
127 {.new_digits: 17, .power_of_five: "277555756156289135105907917022705078125"},
128 {.new_digits: 17, .power_of_five: "1387778780781445675529539585113525390625"},
129 {.new_digits: 18, .power_of_five: "6938893903907228377647697925567626953125"},
130 {.new_digits: 18, .power_of_five: "34694469519536141888238489627838134765625"},
131 {.new_digits: 18, .power_of_five: "173472347597680709441192448139190673828125"},
132 {.new_digits: 19, .power_of_five: "867361737988403547205962240695953369140625"},
133 };
134
135 // The maximum amount we can shift is the number of bits used in the
136 // accumulator, minus the number of bits needed to represent the base (in this
137 // case 4).
138 static constexpr uint32_t MAX_SHIFT_AMOUNT =
139 cpp::numeric_limits<uint64_t>::digits - 4;
140
141 // 800 is an arbitrary number of digits, but should be
142 // large enough for any practical number.
143 static constexpr uint32_t MAX_NUM_DIGITS = 800;
144
145 uint32_t num_digits = 0;
146 int32_t decimal_point = 0;
147 bool truncated = false;
148 uint8_t digits[MAX_NUM_DIGITS];
149
150private:
151 LIBC_INLINE bool should_round_up(int32_t round_to_digit,
152 RoundDirection round) {
153 if (round_to_digit < 0 ||
154 static_cast<uint32_t>(round_to_digit) >= this->num_digits) {
155 return false;
156 }
157
158 // The above condition handles all cases where all of the trailing digits
159 // are zero. In that case, if the rounding mode is up, then this number
160 // should be rounded up. Similarly, if the rounding mode is down, then it
161 // should always round down.
162 if (round == RoundDirection::Up) {
163 return true;
164 } else if (round == RoundDirection::Down) {
165 return false;
166 }
167 // Else round to nearest.
168
169 // If we're right in the middle and there are no extra digits
170 if (this->digits[round_to_digit] == 5 &&
171 static_cast<uint32_t>(round_to_digit + 1) == this->num_digits) {
172
173 // Round up if we've truncated (since that means the result is slightly
174 // higher than what's represented.)
175 if (this->truncated) {
176 return true;
177 }
178
179 // If this exactly halfway, round to even.
180 if (round_to_digit == 0)
181 // When the input is ".5".
182 return false;
183 return this->digits[round_to_digit - 1] % 2 != 0;
184 }
185 // If there are digits after round_to_digit, they must be non-zero since we
186 // trim trailing zeroes after all operations that change digits.
187 return this->digits[round_to_digit] >= 5;
188 }
189
190 // Takes an amount to left shift and returns the number of new digits needed
191 // to store the result based on LEFT_SHIFT_DIGIT_TABLE.
192 LIBC_INLINE uint32_t get_num_new_digits(uint32_t lshift_amount) {
193 const char *power_of_five =
194 LEFT_SHIFT_DIGIT_TABLE[lshift_amount].power_of_five;
195 uint32_t new_digits = LEFT_SHIFT_DIGIT_TABLE[lshift_amount].new_digits;
196 uint32_t digit_index = 0;
197 while (power_of_five[digit_index] != 0) {
198 if (digit_index >= this->num_digits) {
199 return new_digits - 1;
200 }
201 if (this->digits[digit_index] !=
202 internal::b36_char_to_int(ch: power_of_five[digit_index])) {
203 return new_digits -
204 ((this->digits[digit_index] <
205 internal::b36_char_to_int(ch: power_of_five[digit_index]))
206 ? 1
207 : 0);
208 }
209 ++digit_index;
210 }
211 return new_digits;
212 }
213
214 // Trim all trailing 0s
215 LIBC_INLINE void trim_trailing_zeroes() {
216 while (this->num_digits > 0 && this->digits[this->num_digits - 1] == 0) {
217 --this->num_digits;
218 }
219 if (this->num_digits == 0) {
220 this->decimal_point = 0;
221 }
222 }
223
224 // Perform a digitwise binary non-rounding right shift on this value by
225 // shift_amount. The shift_amount can't be more than MAX_SHIFT_AMOUNT to
226 // prevent overflow.
227 LIBC_INLINE void right_shift(uint32_t shift_amount) {
228 uint32_t read_index = 0;
229 uint32_t write_index = 0;
230
231 uint64_t accumulator = 0;
232
233 const uint64_t shift_mask = (uint64_t(1) << shift_amount) - 1;
234
235 // Warm Up phase: we don't have enough digits to start writing, so just
236 // read them into the accumulator.
237 while (accumulator >> shift_amount == 0) {
238 uint64_t read_digit = 0;
239 // If there are still digits to read, read the next one, else the digit is
240 // assumed to be 0.
241 if (read_index < this->num_digits) {
242 read_digit = this->digits[read_index];
243 }
244 accumulator = accumulator * 10 + read_digit;
245 ++read_index;
246 }
247
248 // Shift the decimal point by the number of digits it took to fill the
249 // accumulator.
250 this->decimal_point -= read_index - 1;
251
252 // Middle phase: we have enough digits to write, as well as more digits to
253 // read. Keep reading until we run out of digits.
254 while (read_index < this->num_digits) {
255 uint64_t read_digit = this->digits[read_index];
256 uint64_t write_digit = accumulator >> shift_amount;
257 accumulator &= shift_mask;
258 this->digits[write_index] = static_cast<uint8_t>(write_digit);
259 accumulator = accumulator * 10 + read_digit;
260 ++read_index;
261 ++write_index;
262 }
263
264 // Cool Down phase: All of the readable digits have been read, so just write
265 // the remainder, while treating any more digits as 0.
266 while (accumulator > 0) {
267 uint64_t write_digit = accumulator >> shift_amount;
268 accumulator &= shift_mask;
269 if (write_index < MAX_NUM_DIGITS) {
270 this->digits[write_index] = static_cast<uint8_t>(write_digit);
271 ++write_index;
272 } else if (write_digit > 0) {
273 this->truncated = true;
274 }
275 accumulator = accumulator * 10;
276 }
277 this->num_digits = write_index;
278 this->trim_trailing_zeroes();
279 }
280
281 // Perform a digitwise binary non-rounding left shift on this value by
282 // shift_amount. The shift_amount can't be more than MAX_SHIFT_AMOUNT to
283 // prevent overflow.
284 LIBC_INLINE void left_shift(uint32_t shift_amount) {
285 uint32_t new_digits = this->get_num_new_digits(lshift_amount: shift_amount);
286
287 int32_t read_index = static_cast<int32_t>(this->num_digits - 1);
288 uint32_t write_index = this->num_digits + new_digits;
289
290 uint64_t accumulator = 0;
291
292 // No Warm Up phase. Since we're putting digits in at the top and taking
293 // digits from the bottom we don't have to wait for the accumulator to fill.
294
295 // Middle phase: while we have more digits to read, keep reading as well as
296 // writing.
297 while (read_index >= 0) {
298 accumulator += static_cast<uint64_t>(this->digits[read_index])
299 << shift_amount;
300 uint64_t next_accumulator = accumulator / 10;
301 uint64_t write_digit = accumulator - (10 * next_accumulator);
302 --write_index;
303 if (write_index < MAX_NUM_DIGITS) {
304 this->digits[write_index] = static_cast<uint8_t>(write_digit);
305 } else if (write_digit != 0) {
306 this->truncated = true;
307 }
308 accumulator = next_accumulator;
309 --read_index;
310 }
311
312 // Cool Down phase: there are no more digits to read, so just write the
313 // remaining digits in the accumulator.
314 while (accumulator > 0) {
315 uint64_t next_accumulator = accumulator / 10;
316 uint64_t write_digit = accumulator - (10 * next_accumulator);
317 --write_index;
318 if (write_index < MAX_NUM_DIGITS) {
319 this->digits[write_index] = static_cast<uint8_t>(write_digit);
320 } else if (write_digit != 0) {
321 this->truncated = true;
322 }
323 accumulator = next_accumulator;
324 }
325
326 this->num_digits += new_digits;
327 if (this->num_digits > MAX_NUM_DIGITS) {
328 this->num_digits = MAX_NUM_DIGITS;
329 }
330 this->decimal_point += new_digits;
331 this->trim_trailing_zeroes();
332 }
333
334public:
335 // num_string is assumed to be a string of numeric characters. It doesn't
336 // handle leading spaces.
337 template <typename CharType>
338 LIBC_INLINE HighPrecisionDecimal(
339 const CharType *__restrict num_string,
340 const size_t num_len = cpp::numeric_limits<size_t>::max()) {
341 bool saw_dot = false;
342 size_t num_cur = 0;
343 // This counts the digits in the number, even if there isn't space to store
344 // them all.
345 uint32_t total_digits = 0;
346 while (num_cur < num_len &&
347 (isdigit(num_string[num_cur]) ||
348 num_string[num_cur] == constants<CharType>::DECIMAL_POINT)) {
349 if (num_string[num_cur] == constants<CharType>::DECIMAL_POINT) {
350 if (saw_dot) {
351 break;
352 }
353 this->decimal_point = static_cast<int32_t>(total_digits);
354 saw_dot = true;
355 } else {
356 int digit = b36_char_to_int(num_string[num_cur]);
357 if (digit == 0 && this->num_digits == 0) {
358 --this->decimal_point;
359 ++num_cur;
360 continue;
361 }
362 ++total_digits;
363 if (this->num_digits < MAX_NUM_DIGITS) {
364 this->digits[this->num_digits] = static_cast<uint8_t>(digit);
365 ++this->num_digits;
366 } else if (digit != 0) {
367 this->truncated = true;
368 }
369 }
370 ++num_cur;
371 }
372
373 if (!saw_dot)
374 this->decimal_point = static_cast<int32_t>(total_digits);
375
376 if (num_cur < num_len && tolower(num_string[num_cur]) ==
377 constants<CharType>::DECIMAL_EXPONENT_MARKER) {
378 ++num_cur;
379 if (isdigit(num_string[num_cur]) || get_sign(num_string + num_cur) != 0) {
380 auto result =
381 strtointeger<int32_t>(num_string + num_cur, 10, num_len - num_cur);
382 if (result.has_error()) {
383 // TODO: handle error
384 }
385 int32_t add_to_exponent = result.value;
386
387 // Here we do this operation as int64 to avoid overflow.
388 int64_t temp_exponent = static_cast<int64_t>(this->decimal_point) +
389 static_cast<int64_t>(add_to_exponent);
390
391 // Theoretically these numbers should be MAX_BIASED_EXPONENT for long
392 // double, but that should be ~16,000 which is much less than 1 << 30.
393 if (temp_exponent > (1 << 30)) {
394 temp_exponent = (1 << 30);
395 } else if (temp_exponent < -(1 << 30)) {
396 temp_exponent = -(1 << 30);
397 }
398 this->decimal_point = static_cast<int32_t>(temp_exponent);
399 }
400 }
401
402 this->trim_trailing_zeroes();
403 }
404
405 // Binary shift left (shift_amount > 0) or right (shift_amount < 0)
406 LIBC_INLINE void shift(int shift_amount) {
407 if (shift_amount == 0) {
408 return;
409 }
410 // Left
411 else if (shift_amount > 0) {
412 while (static_cast<uint32_t>(shift_amount) > MAX_SHIFT_AMOUNT) {
413 this->left_shift(shift_amount: MAX_SHIFT_AMOUNT);
414 shift_amount -= MAX_SHIFT_AMOUNT;
415 }
416 this->left_shift(shift_amount: static_cast<uint32_t>(shift_amount));
417 }
418 // Right
419 else {
420 while (static_cast<uint32_t>(shift_amount) < -MAX_SHIFT_AMOUNT) {
421 this->right_shift(shift_amount: MAX_SHIFT_AMOUNT);
422 shift_amount += MAX_SHIFT_AMOUNT;
423 }
424 this->right_shift(shift_amount: static_cast<uint32_t>(-shift_amount));
425 }
426 }
427
428 // Round the number represented to the closest value of unsigned int type T.
429 // This is done ignoring overflow.
430 template <class T>
431 LIBC_INLINE T
432 round_to_integer_type(RoundDirection round = RoundDirection::Nearest) {
433 T result = 0;
434 uint32_t cur_digit = 0;
435
436 while (static_cast<int32_t>(cur_digit) < this->decimal_point &&
437 cur_digit < this->num_digits) {
438 result = result * 10 + (this->digits[cur_digit]);
439 ++cur_digit;
440 }
441
442 // If there are implicit 0s at the end of the number, include those.
443 while (static_cast<int32_t>(cur_digit) < this->decimal_point) {
444 result *= 10;
445 ++cur_digit;
446 }
447 return result +
448 static_cast<T>(this->should_round_up(round_to_digit: this->decimal_point, round));
449 }
450
451 // Extra functions for testing.
452
453 LIBC_INLINE uint8_t *get_digits() { return this->digits; }
454 LIBC_INLINE uint32_t get_num_digits() { return this->num_digits; }
455 LIBC_INLINE int32_t get_decimal_point() { return this->decimal_point; }
456 LIBC_INLINE void set_truncated(bool trunc) { this->truncated = trunc; }
457};
458
459} // namespace internal
460} // namespace LIBC_NAMESPACE_DECL
461
462#endif // LLVM_LIBC_SRC___SUPPORT_HIGH_PRECISION_DECIMAL_H
463