libc/src/__support/high_precision_decimal.h

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