1 //===-- High Precision Decimal ----------------------------------*- C++ -*-===//
2 //
3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4 // See httpss//llvm.org/LICENSE.txt for license information.
5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6 //
7 //===----------------------------------------------------------------------===//
8 
9 #ifndef LIBC_SRC_SUPPORT_HIGH_PRECISION_DECIMAL_H
10 #define LIBC_SRC_SUPPORT_HIGH_PRECISION_DECIMAL_H
11 
12 #include "src/__support/ctype_utils.h"
13 #include "src/__support/str_to_integer.h"
14 #include <stdint.h>
15 
16 namespace __llvm_libc {
17 namespace internal {
18 
19 struct LShiftTableEntry {
20   uint32_t new_digits;
21   char const *power_of_five;
22 };
23 
24 // This is based on the HPD data structure described as part of the Simple
25 // Decimal Conversion algorithm by Nigel Tao, described at this link:
26 // https://nigeltao.github.io/blog/2020/parse-number-f64-simple.html
27 class HighPrecisionDecimal {
28 
29   // This precomputed table speeds up left shifts by having the number of new
30   // digits that will be added by multiplying 5^i by 2^i. If the number is less
31   // than 5^i then it will add one fewer digit. There are only 60 entries since
32   // that's the max shift amount.
33   // This table was generated by the script at
34   // libc/utils/mathtools/GenerateHPDConstants.py
35   static constexpr LShiftTableEntry LEFT_SHIFT_DIGIT_TABLE[] = {
36       {0, ""},
37       {1, "5"},
38       {1, "25"},
39       {1, "125"},
40       {2, "625"},
41       {2, "3125"},
42       {2, "15625"},
43       {3, "78125"},
44       {3, "390625"},
45       {3, "1953125"},
46       {4, "9765625"},
47       {4, "48828125"},
48       {4, "244140625"},
49       {4, "1220703125"},
50       {5, "6103515625"},
51       {5, "30517578125"},
52       {5, "152587890625"},
53       {6, "762939453125"},
54       {6, "3814697265625"},
55       {6, "19073486328125"},
56       {7, "95367431640625"},
57       {7, "476837158203125"},
58       {7, "2384185791015625"},
59       {7, "11920928955078125"},
60       {8, "59604644775390625"},
61       {8, "298023223876953125"},
62       {8, "1490116119384765625"},
63       {9, "7450580596923828125"},
64       {9, "37252902984619140625"},
65       {9, "186264514923095703125"},
66       {10, "931322574615478515625"},
67       {10, "4656612873077392578125"},
68       {10, "23283064365386962890625"},
69       {10, "116415321826934814453125"},
70       {11, "582076609134674072265625"},
71       {11, "2910383045673370361328125"},
72       {11, "14551915228366851806640625"},
73       {12, "72759576141834259033203125"},
74       {12, "363797880709171295166015625"},
75       {12, "1818989403545856475830078125"},
76       {13, "9094947017729282379150390625"},
77       {13, "45474735088646411895751953125"},
78       {13, "227373675443232059478759765625"},
79       {13, "1136868377216160297393798828125"},
80       {14, "5684341886080801486968994140625"},
81       {14, "28421709430404007434844970703125"},
82       {14, "142108547152020037174224853515625"},
83       {15, "710542735760100185871124267578125"},
84       {15, "3552713678800500929355621337890625"},
85       {15, "17763568394002504646778106689453125"},
86       {16, "88817841970012523233890533447265625"},
87       {16, "444089209850062616169452667236328125"},
88       {16, "2220446049250313080847263336181640625"},
89       {16, "11102230246251565404236316680908203125"},
90       {17, "55511151231257827021181583404541015625"},
91       {17, "277555756156289135105907917022705078125"},
92       {17, "1387778780781445675529539585113525390625"},
93       {18, "6938893903907228377647697925567626953125"},
94       {18, "34694469519536141888238489627838134765625"},
95       {18, "173472347597680709441192448139190673828125"},
96       {19, "867361737988403547205962240695953369140625"},
97   };
98 
99   // The maximum amount we can shift is the number of bits used in the
100   // accumulator, minus the number of bits needed to represent the base (in this
101   // case 4).
102   static constexpr uint32_t MAX_SHIFT_AMOUNT = sizeof(uint64_t) - 4;
103 
104   // 800 is an arbitrary number of digits, but should be
105   // large enough for any practical number.
106   static constexpr uint32_t MAX_NUM_DIGITS = 800;
107 
108   uint32_t num_digits = 0;
109   int32_t decimal_point = 0;
110   bool truncated = false;
111   uint8_t digits[MAX_NUM_DIGITS];
112 
113 private:
should_round_up(int32_t roundToDigit)114   bool should_round_up(int32_t roundToDigit) {
115     if (roundToDigit < 0 ||
116         static_cast<uint32_t>(roundToDigit) >= this->num_digits) {
117       return false;
118     }
119 
120     // If we're right in the middle and there are no extra digits
121     if (this->digits[roundToDigit] == 5 &&
122         static_cast<uint32_t>(roundToDigit + 1) == this->num_digits) {
123 
124       // Round up if we've truncated (since that means the result is slightly
125       // higher than what's represented.)
126       if (this->truncated) {
127         return true;
128       }
129 
130       // If this exactly halfway, round to even.
131       if (roundToDigit == 0)
132         // When the input is ".5".
133         return false;
134       return this->digits[roundToDigit - 1] % 2 != 0;
135     }
136     // If there are digits after roundToDigit, they must be non-zero since we
137     // trim trailing zeroes after all operations that change digits.
138     return this->digits[roundToDigit] >= 5;
139   }
140 
141   // Takes an amount to left shift and returns the number of new digits needed
142   // to store the result based on LEFT_SHIFT_DIGIT_TABLE.
get_num_new_digits(uint32_t lShiftAmount)143   uint32_t get_num_new_digits(uint32_t lShiftAmount) {
144     const char *power_of_five =
145         LEFT_SHIFT_DIGIT_TABLE[lShiftAmount].power_of_five;
146     uint32_t new_digits = LEFT_SHIFT_DIGIT_TABLE[lShiftAmount].new_digits;
147     uint32_t digit_index = 0;
148     while (power_of_five[digit_index] != 0) {
149       if (digit_index >= this->num_digits) {
150         return new_digits - 1;
151       }
152       if (this->digits[digit_index] != power_of_five[digit_index] - '0') {
153         return new_digits -
154                ((this->digits[digit_index] < power_of_five[digit_index] - '0')
155                     ? 1
156                     : 0);
157       }
158       ++digit_index;
159     }
160     return new_digits;
161   }
162 
163   // Trim all trailing 0s
trim_trailing_zeroes()164   void trim_trailing_zeroes() {
165     while (this->num_digits > 0 && this->digits[this->num_digits - 1] == 0) {
166       --this->num_digits;
167     }
168     if (this->num_digits == 0) {
169       this->decimal_point = 0;
170     }
171   }
172 
173   // Perform a digitwise binary non-rounding right shift on this value by
174   // shiftAmount. The shiftAmount can't be more than MAX_SHIFT_AMOUNT to prevent
175   // overflow.
right_shift(uint32_t shiftAmount)176   void right_shift(uint32_t shiftAmount) {
177     uint32_t read_index = 0;
178     uint32_t write_index = 0;
179 
180     uint64_t accumulator = 0;
181 
182     const uint64_t shift_mask = (uint64_t(1) << shiftAmount) - 1;
183 
184     // Warm Up phase: we don't have enough digits to start writing, so just
185     // read them into the accumulator.
186     while (accumulator >> shiftAmount == 0) {
187       uint64_t read_digit = 0;
188       // If there are still digits to read, read the next one, else the digit is
189       // assumed to be 0.
190       if (read_index < this->num_digits) {
191         read_digit = this->digits[read_index];
192       }
193       accumulator = accumulator * 10 + read_digit;
194       ++read_index;
195     }
196 
197     // Shift the decimal point by the number of digits it took to fill the
198     // accumulator.
199     this->decimal_point -= read_index - 1;
200 
201     // Middle phase: we have enough digits to write, as well as more digits to
202     // read. Keep reading until we run out of digits.
203     while (read_index < this->num_digits) {
204       uint64_t read_digit = this->digits[read_index];
205       uint64_t write_digit = accumulator >> shiftAmount;
206       accumulator &= shift_mask;
207       this->digits[write_index] = static_cast<uint8_t>(write_digit);
208       accumulator = accumulator * 10 + read_digit;
209       ++read_index;
210       ++write_index;
211     }
212 
213     // Cool Down phase: All of the readable digits have been read, so just write
214     // the remainder, while treating any more digits as 0.
215     while (accumulator > 0) {
216       uint64_t write_digit = accumulator >> shiftAmount;
217       accumulator &= shift_mask;
218       if (write_index < MAX_NUM_DIGITS) {
219         this->digits[write_index] = static_cast<uint8_t>(write_digit);
220         ++write_index;
221       } else if (write_digit > 0) {
222         this->truncated = true;
223       }
224       accumulator = accumulator * 10;
225     }
226     this->num_digits = write_index;
227     this->trim_trailing_zeroes();
228   }
229 
230   // Perform a digitwise binary non-rounding left shift on this value by
231   // shiftAmount. The shiftAmount can't be more than MAX_SHIFT_AMOUNT to prevent
232   // overflow.
left_shift(uint32_t shiftAmount)233   void left_shift(uint32_t shiftAmount) {
234     uint32_t new_digits = this->get_num_new_digits(shiftAmount);
235 
236     int32_t read_index = this->num_digits - 1;
237     uint32_t write_index = this->num_digits + new_digits;
238 
239     uint64_t accumulator = 0;
240 
241     // No Warm Up phase. Since we're putting digits in at the top and taking
242     // digits from the bottom we don't have to wait for the accumulator to fill.
243 
244     // Middle phase: while we have more digits to read, keep reading as well as
245     // writing.
246     while (read_index >= 0) {
247       accumulator += static_cast<uint64_t>(this->digits[read_index])
248                      << shiftAmount;
249       uint64_t next_accumulator = accumulator / 10;
250       uint64_t write_digit = accumulator - (10 * next_accumulator);
251       --write_index;
252       if (write_index < MAX_NUM_DIGITS) {
253         this->digits[write_index] = static_cast<uint8_t>(write_digit);
254       } else if (write_digit != 0) {
255         this->truncated = true;
256       }
257       accumulator = next_accumulator;
258       --read_index;
259     }
260 
261     // Cool Down phase: there are no more digits to read, so just write the
262     // remaining digits in the accumulator.
263     while (accumulator > 0) {
264       uint64_t next_accumulator = accumulator / 10;
265       uint64_t write_digit = accumulator - (10 * next_accumulator);
266       --write_index;
267       if (write_index < MAX_NUM_DIGITS) {
268         this->digits[write_index] = static_cast<uint8_t>(write_digit);
269       } else if (write_digit != 0) {
270         this->truncated = true;
271       }
272       accumulator = next_accumulator;
273     }
274 
275     this->num_digits += new_digits;
276     if (this->num_digits > MAX_NUM_DIGITS) {
277       this->num_digits = MAX_NUM_DIGITS;
278     }
279     this->decimal_point += new_digits;
280     this->trim_trailing_zeroes();
281   }
282 
283 public:
284   // numString is assumed to be a string of numeric characters. It doesn't
285   // handle leading spaces.
HighPrecisionDecimal(const char * __restrict numString)286   HighPrecisionDecimal(const char *__restrict numString) {
287     bool saw_dot = false;
288     // This counts the digits in the number, even if there isn't space to store
289     // them all.
290     uint32_t total_digits = 0;
291     while (isdigit(*numString) || *numString == '.') {
292       if (*numString == '.') {
293         if (saw_dot) {
294           break;
295         }
296         this->decimal_point = total_digits;
297         saw_dot = true;
298       } else {
299         if (*numString == '0' && this->num_digits == 0) {
300           --this->decimal_point;
301           ++numString;
302           continue;
303         }
304         ++total_digits;
305         if (this->num_digits < MAX_NUM_DIGITS) {
306           this->digits[this->num_digits] =
307               static_cast<uint8_t>(*numString - '0');
308           ++this->num_digits;
309         } else if (*numString != '0') {
310           this->truncated = true;
311         }
312       }
313       ++numString;
314     }
315 
316     if (!saw_dot)
317       this->decimal_point = total_digits;
318 
319     if ((*numString | 32) == 'e') {
320       ++numString;
321       if (isdigit(*numString) || *numString == '+' || *numString == '-') {
322         int32_t add_to_exp = strtointeger<int32_t>(numString, nullptr, 10);
323         if (add_to_exp > 100000) {
324           add_to_exp = 100000;
325         } else if (add_to_exp < -100000) {
326           add_to_exp = -100000;
327         }
328         this->decimal_point += add_to_exp;
329       }
330     }
331 
332     this->trim_trailing_zeroes();
333   }
334 
335   // Binary shift left (shiftAmount > 0) or right (shiftAmount < 0)
shift(int shiftAmount)336   void shift(int shiftAmount) {
337     if (shiftAmount == 0) {
338       return;
339     }
340     // Left
341     else if (shiftAmount > 0) {
342       while (static_cast<uint32_t>(shiftAmount) > MAX_SHIFT_AMOUNT) {
343         this->left_shift(MAX_SHIFT_AMOUNT);
344         shiftAmount -= MAX_SHIFT_AMOUNT;
345       }
346       this->left_shift(shiftAmount);
347     }
348     // Right
349     else {
350       while (static_cast<uint32_t>(shiftAmount) < -MAX_SHIFT_AMOUNT) {
351         this->right_shift(MAX_SHIFT_AMOUNT);
352         shiftAmount += MAX_SHIFT_AMOUNT;
353       }
354       this->right_shift(-shiftAmount);
355     }
356   }
357 
358   // Round the number represented to the closest value of unsigned int type T.
359   // This is done ignoring overflow.
round_to_integer_type()360   template <class T> T round_to_integer_type() {
361     T result = 0;
362     uint32_t cur_digit = 0;
363 
364     while (static_cast<int32_t>(cur_digit) < this->decimal_point &&
365            cur_digit < this->num_digits) {
366       result = result * 10 + (this->digits[cur_digit]);
367       ++cur_digit;
368     }
369 
370     // If there are implicit 0s at the end of the number, include those.
371     while (static_cast<int32_t>(cur_digit) < this->decimal_point) {
372       result *= 10;
373       ++cur_digit;
374     }
375     if (this->should_round_up(this->decimal_point)) {
376       ++result;
377     }
378     return result;
379   }
380 
381   // Extra functions for testing.
382 
get_digits()383   uint8_t *get_digits() { return this->digits; }
get_num_digits()384   uint32_t get_num_digits() { return this->num_digits; }
get_decimal_point()385   int32_t get_decimal_point() { return this->decimal_point; }
set_truncated(bool trunc)386   void set_truncated(bool trunc) { this->truncated = trunc; }
387 };
388 
389 } // namespace internal
390 } // namespace __llvm_libc
391 
392 #endif // LIBC_SRC_SUPPORT_HIGH_PRECISION_DECIMAL_H
393