[libc][NFC] Add supporting class for atof implementation

This change adds the High Precision Decimal described here:
https://nigeltao.github.io/blog/2020/parse-number-f64-simple.html
It will be used for the atof implementation later, but is complete and
tested now.

The code is inspired by the golang implmentation of the HPD class, which
can be found here: https://github.com/golang/go/blob/release-branch.go1.16/src/strconv/decimal.go

Reviewed By: sivachandra

Differential Revision: https://reviews.llvm.org/D110454

diff --git a/libc/src/__support/CMakeLists.txt b/libc/src/__support/CMakeLists.txt
index a09ca4b..2699b55 100644 (file)
--- a/libc/src/__support/CMakeLists.txt
+++ b/libc/src/__support/CMakeLists.txt
@@ -13,11 +13,19 @@ add_header_library(
 )
 
 add_header_library(
+  high_precision_decimal
+  HDRS
+    high_precision_decimal.h
+
+)
+
+add_header_library(
   str_conv_utils
   HDRS
     str_conv_utils.h
   DEPENDS
     .ctype_utils
+    .high_precision_decimal
     libc.include.errno
     libc.src.errno.__errno_location
     libc.utils.CPP.standalone_cpp

diff --git a/libc/src/__support/high_precision_decimal.h b/libc/src/__support/high_precision_decimal.h
new file mode 100644 (file)
index 0000000..d61ba80
--- /dev/null
+++ b/libc/src/__support/high_precision_decimal.h
@@ -0,0 +1,378 @@
+//===-- High Precision Decimal ----------------------------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See httpss//llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LIBC_SRC_SUPPORT_HIGH_PRECISION_DECIMAL_H
+#define LIBC_SRC_SUPPORT_HIGH_PRECISION_DECIMAL_H
+
+#include "src/__support/ctype_utils.h"
+#include "src/__support/str_conv_utils.h"
+#include <stdint.h>
+
+namespace __llvm_libc {
+namespace internal {
+
+struct LShiftTableEntry {
+  uint32_t newDigits;
+  char const *powerOfFive;
+};
+
+// This is based on the HPD data structure described as part of the Simple
+// Decimal Conversion algorithm by Nigel Tao, described at this link:
+// https://nigeltao.github.io/blog/2020/parse-number-f64-simple.html
+class HighPrecsisionDecimal {
+
+  // This precomputed table speeds up left shifts by having the number of new
+  // digits that will be added by multiplying 5^i by 2^i. If the number is less
+  // than 5^i then it will add one fewer digit. There are only 60 entries since
+  // that's the max shift amount.
+  // This table was generated by the script at
+  // libc/utils/mathtools/GenerateHPDConstants.py
+  static constexpr LShiftTableEntry LEFT_SHIFT_DIGIT_TABLE[] = {
+      {0, ""},
+      {1, "5"},
+      {1, "25"},
+      {1, "125"},
+      {2, "625"},
+      {2, "3125"},
+      {2, "15625"},
+      {3, "78125"},
+      {3, "390625"},
+      {3, "1953125"},
+      {4, "9765625"},
+      {4, "48828125"},
+      {4, "244140625"},
+      {4, "1220703125"},
+      {5, "6103515625"},
+      {5, "30517578125"},
+      {5, "152587890625"},
+      {6, "762939453125"},
+      {6, "3814697265625"},
+      {6, "19073486328125"},
+      {7, "95367431640625"},
+      {7, "476837158203125"},
+      {7, "2384185791015625"},
+      {7, "11920928955078125"},
+      {8, "59604644775390625"},
+      {8, "298023223876953125"},
+      {8, "1490116119384765625"},
+      {9, "7450580596923828125"},
+      {9, "37252902984619140625"},
+      {9, "186264514923095703125"},
+      {10, "931322574615478515625"},
+      {10, "4656612873077392578125"},
+      {10, "23283064365386962890625"},
+      {10, "116415321826934814453125"},
+      {11, "582076609134674072265625"},
+      {11, "2910383045673370361328125"},
+      {11, "14551915228366851806640625"},
+      {12, "72759576141834259033203125"},
+      {12, "363797880709171295166015625"},
+      {12, "1818989403545856475830078125"},
+      {13, "9094947017729282379150390625"},
+      {13, "45474735088646411895751953125"},
+      {13, "227373675443232059478759765625"},
+      {13, "1136868377216160297393798828125"},
+      {14, "5684341886080801486968994140625"},
+      {14, "28421709430404007434844970703125"},
+      {14, "142108547152020037174224853515625"},
+      {15, "710542735760100185871124267578125"},
+      {15, "3552713678800500929355621337890625"},
+      {15, "17763568394002504646778106689453125"},
+      {16, "88817841970012523233890533447265625"},
+      {16, "444089209850062616169452667236328125"},
+      {16, "2220446049250313080847263336181640625"},
+      {16, "11102230246251565404236316680908203125"},
+      {17, "55511151231257827021181583404541015625"},
+      {17, "277555756156289135105907917022705078125"},
+      {17, "1387778780781445675529539585113525390625"},
+      {18, "6938893903907228377647697925567626953125"},
+      {18, "34694469519536141888238489627838134765625"},
+      {18, "173472347597680709441192448139190673828125"},
+      {19, "867361737988403547205962240695953369140625"},
+  };
+
+  // The maximum amount we can shift is the number of bits used in the
+  // accumulator, minus the number of bits needed to represent the base (in this
+  // case 4).
+  static constexpr uint32_t MAX_SHIFT_AMOUNT = sizeof(uint64_t) - 4;
+
+  // 800 is an arbitrary number of digits, but should be
+  // large enough for any practical number.
+  static constexpr uint32_t MAX_NUM_DIGITS = 800;
+
+  uint32_t numDigits = 0;
+  int32_t decimalPoint = 0;
+  bool truncated = false;
+  uint8_t digits[MAX_NUM_DIGITS];
+
+private:
+  bool shouldRoundUp(uint32_t roundToDigit) {
+    if (roundToDigit < 0 || roundToDigit >= this->numDigits) {
+      return false;
+    }
+
+    // If we're right in the middle and there are no extra digits
+    if (this->digits[roundToDigit] == 5 &&
+        roundToDigit + 1 == this->numDigits) {
+
+      // Round up if we've truncated (since that means the result is slightly
+      // higher than what's represented.)
+      if (this->truncated) {
+        return true;
+      }
+
+      // If this exactly halfway, round to even.
+      return this->digits[roundToDigit - 1] % 2 != 0;
+    }
+    // If there are digits after roundToDigit, they must be non-zero since we
+    // trim trailing zeroes after all operations that change digits.
+    return this->digits[roundToDigit] >= 5;
+  }
+
+  // Takes an amount to left shift and returns the number of new digits needed
+  // to store the result based on LEFT_SHIFT_DIGIT_TABLE.
+  uint32_t getNumNewDigits(uint32_t lShiftAmount) {
+    const char *powerOfFive = LEFT_SHIFT_DIGIT_TABLE[lShiftAmount].powerOfFive;
+    uint32_t newDigits = LEFT_SHIFT_DIGIT_TABLE[lShiftAmount].newDigits;
+    uint32_t digitIndex = 0;
+    while (powerOfFive[digitIndex] != 0) {
+      if (digitIndex >= this->numDigits) {
+        return newDigits - 1;
+      }
+      if (this->digits[digitIndex] != powerOfFive[digitIndex] - '0') {
+        return newDigits -
+               ((this->digits[digitIndex] < powerOfFive[digitIndex] - '0') ? 1
+                                                                           : 0);
+      }
+      ++digitIndex;
+    }
+    return newDigits;
+  }
+
+  // Trim all trailing 0s
+  void trimTrailingZeroes() {
+    while (this->numDigits > 0 && this->digits[this->numDigits - 1] == 0) {
+      --this->numDigits;
+    }
+    if (this->numDigits == 0) {
+      this->decimalPoint = 0;
+    }
+  }
+
+  // Perform a digitwise binary non-rounding right shift on this value by
+  // shiftAmount. The shiftAmount can't be more than MAX_SHIFT_AMOUNT to prevent
+  // overflow.
+  void rightShift(uint32_t shiftAmount) {
+    uint32_t readIndex = 0;
+    uint32_t writeIndex = 0;
+
+    uint64_t accumulator = 0;
+
+    const uint64_t shiftMask = (uint64_t(1) << shiftAmount) - 1;
+
+    // Warm Up phase: we don't have enough digits to start writing, so just
+    // read them into the accumulator.
+    while (accumulator >> shiftAmount == 0) {
+      uint64_t readDigit = 0;
+      // If there are still digits to read, read the next one, else the digit is
+      // assumed to be 0.
+      if (readIndex < this->numDigits) {
+        readDigit = this->digits[readIndex];
+      }
+      accumulator = accumulator * 10 + readDigit;
+      ++readIndex;
+    }
+
+    // Shift the decimal point by the number of digits it took to fill the
+    // accumulator.
+    this->decimalPoint -= readIndex - 1;
+
+    // Middle phase: we have enough digits to write, as well as more digits to
+    // read. Keep reading until we run out of digits.
+    while (readIndex < this->numDigits) {
+      uint64_t readDigit = this->digits[readIndex];
+      uint64_t writeDigit = accumulator >> shiftAmount;
+      accumulator &= shiftMask;
+      this->digits[writeIndex] = static_cast<uint8_t>(writeDigit);
+      accumulator = accumulator * 10 + readDigit;
+      ++readIndex;
+      ++writeIndex;
+    }
+
+    // Cool Down phase: All of the readable digits have been read, so just write
+    // the remainder, while treating any more digits as 0.
+    while (accumulator > 0) {
+      uint64_t writeDigit = accumulator >> shiftAmount;
+      accumulator &= shiftMask;
+      if (writeIndex < MAX_NUM_DIGITS) {
+        this->digits[writeIndex] = static_cast<uint8_t>(writeDigit);
+        ++writeIndex;
+      } else if (writeDigit > 0) {
+        this->truncated = true;
+      }
+      accumulator = accumulator * 10;
+    }
+    this->numDigits = writeIndex;
+    this->trimTrailingZeroes();
+  }
+
+  // Perform a digitwise binary non-rounding left shift on this value by
+  // shiftAmount. The shiftAmount can't be more than MAX_SHIFT_AMOUNT to prevent
+  // overflow.
+  void leftShift(uint32_t shiftAmount) {
+    uint32_t newDigits = this->getNumNewDigits(shiftAmount);
+
+    int32_t readIndex = this->numDigits - 1;
+    uint32_t writeIndex = this->numDigits + newDigits;
+
+    uint64_t accumulator = 0;
+
+    // No Warm Up phase. Since we're putting digits in at the top and taking
+    // digits from the bottom we don't have to wait for the accumulator to fill.
+
+    // Middle phase: while we have more digits to read, keep reading as well as
+    // writing.
+    while (readIndex >= 0) {
+      accumulator += static_cast<uint64_t>(this->digits[readIndex])
+                     << shiftAmount;
+      uint64_t nextAccumulator = accumulator / 10;
+      uint64_t writeDigit = accumulator - (10 * nextAccumulator);
+      --writeIndex;
+      if (writeIndex < MAX_NUM_DIGITS) {
+        this->digits[writeIndex] = static_cast<uint8_t>(writeDigit);
+      } else if (writeDigit != 0) {
+        this->truncated = true;
+      }
+      accumulator = nextAccumulator;
+      --readIndex;
+    }
+
+    // Cool Down phase: there are no more digits to read, so just write the
+    // remaining digits in the accumulator.
+    while (accumulator > 0) {
+      uint64_t nextAccumulator = accumulator / 10;
+      uint64_t writeDigit = accumulator - (10 * nextAccumulator);
+      --writeIndex;
+      if (writeIndex < MAX_NUM_DIGITS) {
+        this->digits[writeIndex] = static_cast<uint8_t>(writeDigit);
+      } else if (writeDigit != 0) {
+        this->truncated = true;
+      }
+      accumulator = nextAccumulator;
+    }
+
+    this->numDigits += newDigits;
+    if (this->numDigits > MAX_NUM_DIGITS) {
+      this->numDigits = MAX_NUM_DIGITS;
+    }
+    this->decimalPoint += newDigits;
+    this->trimTrailingZeroes();
+  }
+
+public:
+  // numString is assumed to be a string of numeric characters. It doesn't
+  // handle leading spaces.
+  HighPrecsisionDecimal(const char *__restrict numString) {
+    bool sawDot = false;
+    bool sawDigit = false;
+    while (isdigit(*numString) || *numString == '.') {
+      if (*numString == '.') {
+        if (sawDot) {
+          break;
+        }
+        this->decimalPoint = this->numDigits;
+        sawDot = true;
+      } else {
+        sawDigit = true;
+        if (*numString == '0' && this->numDigits == 0) {
+          --this->decimalPoint;
+          continue;
+        }
+        if (this->numDigits < MAX_NUM_DIGITS) {
+          this->digits[this->numDigits] = *numString - '0';
+          ++this->numDigits;
+        } else if (*numString != '0') {
+          this->truncated = true;
+        }
+      }
+      ++numString;
+    }
+
+    if (!sawDot) {
+      this->decimalPoint = this->numDigits;
+    }
+
+    if ((*numString | 32) == 'e') {
+      ++numString;
+      if (isdigit(*numString) || *numString == '+' || *numString == '-') {
+        int32_t addToExp = strtointeger<int32_t>(numString, nullptr, 10);
+        this->decimalPoint += addToExp;
+      }
+    }
+
+    this->trimTrailingZeroes();
+  }
+
+  // Binary shift left (shiftAmount > 0) or right (shiftAmount < 0)
+  void shift(int shiftAmount) {
+    if (shiftAmount == 0) {
+      return;
+    }
+    // Left
+    else if (shiftAmount > 0) {
+      while (static_cast<uint32_t>(shiftAmount) > MAX_SHIFT_AMOUNT) {
+        this->leftShift(MAX_SHIFT_AMOUNT);
+        shiftAmount -= MAX_SHIFT_AMOUNT;
+      }
+      this->leftShift(shiftAmount);
+    }
+    // Right
+    else {
+      while (static_cast<uint32_t>(shiftAmount) < -MAX_SHIFT_AMOUNT) {
+        this->rightShift(MAX_SHIFT_AMOUNT);
+        shiftAmount += MAX_SHIFT_AMOUNT;
+      }
+      this->rightShift(-shiftAmount);
+    }
+  }
+
+  // Round the number represented to the closest value of unsigned int type T.
+  // This is done ignoring overflow.
+  template <class T> T roundToIntegerType() {
+    T result = 0;
+    uint32_t curDigit = 0;
+
+    while (static_cast<int32_t>(curDigit) < this->decimalPoint &&
+           curDigit < this->numDigits) {
+      result = result * 10 + (this->digits[curDigit]);
+      ++curDigit;
+    }
+
+    // If there are implicit 0s at the end of the number, include those.
+    while (static_cast<int32_t>(curDigit) < this->decimalPoint) {
+      result *= 10;
+      ++curDigit;
+    }
+    if (this->shouldRoundUp(this->decimalPoint)) {
+      ++result;
+    }
+    return result;
+  }
+
+  // Extra functions for testing.
+
+  uint8_t *getDigits() { return this->digits; }
+  uint32_t getNumDigits() { return this->numDigits; }
+  int32_t getDecimalPoint() { return this->decimalPoint; }
+  void setTruncated(bool trunc) { this->truncated = trunc; }
+};
+
+} // namespace internal
+} // namespace __llvm_libc
+
+#endif // LIBC_SRC_SUPPORT_HIGH_PRECISION_DECIMAL_H

diff --git a/libc/test/src/__support/CMakeLists.txt b/libc/test/src/__support/CMakeLists.txt
index 813e413..7ae4d90 100644 (file)
--- a/libc/test/src/__support/CMakeLists.txt
+++ b/libc/test/src/__support/CMakeLists.txt
@@ -9,3 +9,13 @@ add_libc_unittest(
   DEPENDS
     libc.src.__support.common
 )
+
+add_libc_unittest(
+  high_precision_decimal_test
+  SUITE
+    libc_support_unittests
+  SRCS
+  high_precision_decimal_test.cpp
+  DEPENDS
+    libc.src.__support.high_precision_decimal
+)

diff --git a/libc/test/src/__support/high_precision_decimal_test.cpp b/libc/test/src/__support/high_precision_decimal_test.cpp
new file mode 100644 (file)
index 0000000..0234426
--- /dev/null
+++ b/libc/test/src/__support/high_precision_decimal_test.cpp
@@ -0,0 +1,381 @@
+//===-- Unittests for high_precision_decimal ------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/__support/high_precision_decimal.h"
+
+#include "utils/UnitTest/Test.h"
+
+TEST(LlvmLibcHighPrecisionDecimalTest, BasicInit) {
+  __llvm_libc::internal::HighPrecsisionDecimal hpd =
+      __llvm_libc::internal::HighPrecsisionDecimal("1.2345");
+  uint8_t *digits = hpd.getDigits();
+
+  EXPECT_EQ(digits[0], uint8_t(1));
+  EXPECT_EQ(digits[1], uint8_t(2));
+  EXPECT_EQ(digits[2], uint8_t(3));
+  EXPECT_EQ(digits[3], uint8_t(4));
+  EXPECT_EQ(digits[4], uint8_t(5));
+  EXPECT_EQ(hpd.getNumDigits(), 5u);
+  EXPECT_EQ(hpd.getDecimalPoint(), 1);
+}
+
+TEST(LlvmLibcHighPrecisionDecimalTest, BasicShift) {
+  __llvm_libc::internal::HighPrecsisionDecimal hpd =
+      __llvm_libc::internal::HighPrecsisionDecimal("1");
+  uint8_t *digits = hpd.getDigits();
+
+  hpd.shift(1); // shift left 1, equal to multiplying by 2.
+
+  EXPECT_EQ(digits[0], uint8_t(2));
+  EXPECT_EQ(hpd.getNumDigits(), 1u);
+  EXPECT_EQ(hpd.getDecimalPoint(), 1);
+}
+
+TEST(LlvmLibcHighPrecisionDecimalTest, SmallShift) {
+  __llvm_libc::internal::HighPrecsisionDecimal hpd =
+      __llvm_libc::internal::HighPrecsisionDecimal("1.2345");
+  uint8_t *digits = hpd.getDigits();
+
+  hpd.shift(-1); // shift right one, equal to dividing by 2
+  // result should be 0.61725
+
+  EXPECT_EQ(digits[0], uint8_t(6));
+  EXPECT_EQ(digits[1], uint8_t(1));
+  EXPECT_EQ(digits[2], uint8_t(7));
+  EXPECT_EQ(digits[3], uint8_t(2));
+  EXPECT_EQ(digits[4], uint8_t(5));
+  EXPECT_EQ(hpd.getNumDigits(), 5u);
+  EXPECT_EQ(hpd.getDecimalPoint(), 0);
+
+  hpd.shift(1); // shift left one, equal to multiplying by 2
+  // result should be 1.2345 again
+
+  EXPECT_EQ(digits[0], uint8_t(1));
+  EXPECT_EQ(digits[1], uint8_t(2));
+  EXPECT_EQ(digits[2], uint8_t(3));
+  EXPECT_EQ(digits[3], uint8_t(4));
+  EXPECT_EQ(digits[4], uint8_t(5));
+  EXPECT_EQ(hpd.getNumDigits(), 5u);
+  EXPECT_EQ(hpd.getDecimalPoint(), 1);
+
+  hpd.shift(1); // shift left one again
+  // result should be 2.469
+
+  EXPECT_EQ(digits[0], uint8_t(2));
+  EXPECT_EQ(digits[1], uint8_t(4));
+  EXPECT_EQ(digits[2], uint8_t(6));
+  EXPECT_EQ(digits[3], uint8_t(9));
+  EXPECT_EQ(hpd.getNumDigits(), 4u);
+  EXPECT_EQ(hpd.getDecimalPoint(), 1);
+
+  hpd.shift(-1); // shift right one again
+  // result should be 1.2345 again
+
+  EXPECT_EQ(digits[0], uint8_t(1));
+  EXPECT_EQ(digits[1], uint8_t(2));
+  EXPECT_EQ(digits[2], uint8_t(3));
+  EXPECT_EQ(digits[3], uint8_t(4));
+  EXPECT_EQ(digits[4], uint8_t(5));
+  EXPECT_EQ(hpd.getNumDigits(), 5u);
+  EXPECT_EQ(hpd.getDecimalPoint(), 1);
+}
+
+TEST(LlvmLibcHighPrecisionDecimalTest, MediumShift) {
+  __llvm_libc::internal::HighPrecsisionDecimal hpd =
+      __llvm_libc::internal::HighPrecsisionDecimal(".299792458");
+  uint8_t *digits = hpd.getDigits();
+
+  hpd.shift(-3); // shift right three, equal to dividing by 8
+  // result should be 0.03747405725
+
+  EXPECT_EQ(digits[0], uint8_t(3));
+  EXPECT_EQ(digits[1], uint8_t(7));
+  EXPECT_EQ(digits[2], uint8_t(4));
+  EXPECT_EQ(digits[3], uint8_t(7));
+  EXPECT_EQ(digits[4], uint8_t(4));
+  EXPECT_EQ(digits[5], uint8_t(0));
+  EXPECT_EQ(digits[6], uint8_t(5));
+  EXPECT_EQ(digits[7], uint8_t(7));
+  EXPECT_EQ(digits[8], uint8_t(2));
+  EXPECT_EQ(digits[9], uint8_t(5));
+  EXPECT_EQ(hpd.getNumDigits(), 10u);
+  EXPECT_EQ(hpd.getDecimalPoint(), -1);
+
+  hpd.shift(3); // shift left three, equal to multiplying by 8
+  // result should be 0.299792458 again
+
+  EXPECT_EQ(digits[0], uint8_t(2));
+  EXPECT_EQ(digits[1], uint8_t(9));
+  EXPECT_EQ(digits[2], uint8_t(9));
+  EXPECT_EQ(digits[3], uint8_t(7));
+  EXPECT_EQ(digits[4], uint8_t(9));
+  EXPECT_EQ(digits[5], uint8_t(2));
+  EXPECT_EQ(digits[6], uint8_t(4));
+  EXPECT_EQ(digits[7], uint8_t(5));
+  EXPECT_EQ(digits[8], uint8_t(8));
+  EXPECT_EQ(hpd.getNumDigits(), 9u);
+  EXPECT_EQ(hpd.getDecimalPoint(), 0);
+}
+
+TEST(LlvmLibcHighPrecisionDecimalTest, BigShift) {
+  __llvm_libc::internal::HighPrecsisionDecimal hpd =
+      __llvm_libc::internal::HighPrecsisionDecimal(".299792458");
+  uint8_t *digits = hpd.getDigits();
+
+  hpd.shift(-29); // shift right 29, equal to dividing by 536,870,912
+  // result should be 0.0000000005584069676697254180908203125
+
+  EXPECT_EQ(digits[0], uint8_t(5));
+  EXPECT_EQ(digits[1], uint8_t(5));
+  EXPECT_EQ(digits[2], uint8_t(8));
+  EXPECT_EQ(digits[3], uint8_t(4));
+  EXPECT_EQ(digits[4], uint8_t(0));
+  EXPECT_EQ(digits[5], uint8_t(6));
+  EXPECT_EQ(digits[6], uint8_t(9));
+  EXPECT_EQ(digits[7], uint8_t(6));
+  EXPECT_EQ(digits[8], uint8_t(7));
+  EXPECT_EQ(digits[9], uint8_t(6));
+  EXPECT_EQ(digits[10], uint8_t(6));
+  EXPECT_EQ(digits[11], uint8_t(9));
+  EXPECT_EQ(digits[12], uint8_t(7));
+  EXPECT_EQ(digits[13], uint8_t(2));
+  EXPECT_EQ(digits[14], uint8_t(5));
+  EXPECT_EQ(digits[15], uint8_t(4));
+  EXPECT_EQ(digits[16], uint8_t(1));
+  EXPECT_EQ(digits[17], uint8_t(8));
+  EXPECT_EQ(digits[18], uint8_t(0));
+  EXPECT_EQ(digits[19], uint8_t(9));
+  EXPECT_EQ(digits[20], uint8_t(0));
+  EXPECT_EQ(digits[21], uint8_t(8));
+  EXPECT_EQ(digits[22], uint8_t(2));
+  EXPECT_EQ(digits[23], uint8_t(0));
+  EXPECT_EQ(digits[24], uint8_t(3));
+  EXPECT_EQ(digits[25], uint8_t(1));
+  EXPECT_EQ(digits[26], uint8_t(2));
+  EXPECT_EQ(digits[27], uint8_t(5));
+  EXPECT_EQ(hpd.getNumDigits(), 28u);
+  EXPECT_EQ(hpd.getDecimalPoint(), -9);
+
+  hpd.shift(29); // shift left 29, equal to multiplying by 536,870,912
+  // result should be 0.299792458 again
+
+  EXPECT_EQ(digits[0], uint8_t(2));
+  EXPECT_EQ(digits[1], uint8_t(9));
+  EXPECT_EQ(digits[2], uint8_t(9));
+  EXPECT_EQ(digits[3], uint8_t(7));
+  EXPECT_EQ(digits[4], uint8_t(9));
+  EXPECT_EQ(digits[5], uint8_t(2));
+  EXPECT_EQ(digits[6], uint8_t(4));
+  EXPECT_EQ(digits[7], uint8_t(5));
+  EXPECT_EQ(digits[8], uint8_t(8));
+  EXPECT_EQ(hpd.getNumDigits(), 9u);
+  EXPECT_EQ(hpd.getDecimalPoint(), 0);
+}
+
+TEST(LlvmLibcHighPrecisionDecimalTest, BigShiftInSteps) {
+  __llvm_libc::internal::HighPrecsisionDecimal hpd =
+      __llvm_libc::internal::HighPrecsisionDecimal("1");
+  uint8_t *digits = hpd.getDigits();
+
+  hpd.shift(60); // shift left 60, equal to multiplying by
+                 // 1152921504606846976.
+
+  EXPECT_EQ(digits[0], uint8_t(1));
+  EXPECT_EQ(digits[1], uint8_t(1));
+  EXPECT_EQ(digits[2], uint8_t(5));
+  EXPECT_EQ(digits[3], uint8_t(2));
+  EXPECT_EQ(digits[4], uint8_t(9));
+  EXPECT_EQ(digits[5], uint8_t(2));
+  EXPECT_EQ(digits[6], uint8_t(1));
+  EXPECT_EQ(digits[7], uint8_t(5));
+  EXPECT_EQ(digits[8], uint8_t(0));
+  EXPECT_EQ(digits[9], uint8_t(4));
+  EXPECT_EQ(digits[10], uint8_t(6));
+  EXPECT_EQ(digits[11], uint8_t(0));
+  EXPECT_EQ(digits[12], uint8_t(6));
+  EXPECT_EQ(digits[13], uint8_t(8));
+  EXPECT_EQ(digits[14], uint8_t(4));
+  EXPECT_EQ(digits[15], uint8_t(6));
+  EXPECT_EQ(digits[16], uint8_t(9));
+  EXPECT_EQ(digits[17], uint8_t(7));
+  EXPECT_EQ(digits[18], uint8_t(6));
+  EXPECT_EQ(hpd.getNumDigits(), 19u);
+  EXPECT_EQ(hpd.getDecimalPoint(), 19);
+
+  hpd.shift(40); // shift left 40, equal to multiplying by
+                 // 1099511627776. Result should be 2^100
+
+  EXPECT_EQ(digits[0], uint8_t(1));
+  EXPECT_EQ(digits[1], uint8_t(2));
+  EXPECT_EQ(digits[2], uint8_t(6));
+  EXPECT_EQ(digits[3], uint8_t(7));
+  EXPECT_EQ(digits[4], uint8_t(6));
+  EXPECT_EQ(digits[5], uint8_t(5));
+  EXPECT_EQ(digits[6], uint8_t(0));
+  EXPECT_EQ(digits[7], uint8_t(6));
+  EXPECT_EQ(digits[8], uint8_t(0));
+  EXPECT_EQ(digits[9], uint8_t(0));
+  EXPECT_EQ(digits[10], uint8_t(2));
+  EXPECT_EQ(digits[11], uint8_t(2));
+  EXPECT_EQ(digits[12], uint8_t(8));
+  EXPECT_EQ(digits[13], uint8_t(2));
+  EXPECT_EQ(digits[14], uint8_t(2));
+  EXPECT_EQ(digits[15], uint8_t(9));
+  EXPECT_EQ(digits[16], uint8_t(4));
+  EXPECT_EQ(digits[17], uint8_t(0));
+  EXPECT_EQ(digits[18], uint8_t(1));
+  EXPECT_EQ(digits[19], uint8_t(4));
+  EXPECT_EQ(digits[20], uint8_t(9));
+  EXPECT_EQ(digits[21], uint8_t(6));
+  EXPECT_EQ(digits[22], uint8_t(7));
+  EXPECT_EQ(digits[23], uint8_t(0));
+  EXPECT_EQ(digits[24], uint8_t(3));
+  EXPECT_EQ(digits[25], uint8_t(2));
+  EXPECT_EQ(digits[26], uint8_t(0));
+  EXPECT_EQ(digits[27], uint8_t(5));
+  EXPECT_EQ(digits[28], uint8_t(3));
+  EXPECT_EQ(digits[29], uint8_t(7));
+  EXPECT_EQ(digits[30], uint8_t(6));
+
+  EXPECT_EQ(hpd.getNumDigits(), 31u);
+  EXPECT_EQ(hpd.getDecimalPoint(), 31);
+
+  hpd.shift(-60); // shift right 60, equal to dividing by
+                  // 1152921504606846976. Result should be 2^40
+
+  EXPECT_EQ(digits[0], uint8_t(1));
+  EXPECT_EQ(digits[1], uint8_t(0));
+  EXPECT_EQ(digits[2], uint8_t(9));
+  EXPECT_EQ(digits[3], uint8_t(9));
+  EXPECT_EQ(digits[4], uint8_t(5));
+  EXPECT_EQ(digits[5], uint8_t(1));
+  EXPECT_EQ(digits[6], uint8_t(1));
+  EXPECT_EQ(digits[7], uint8_t(6));
+  EXPECT_EQ(digits[8], uint8_t(2));
+  EXPECT_EQ(digits[9], uint8_t(7));
+  EXPECT_EQ(digits[10], uint8_t(7));
+  EXPECT_EQ(digits[11], uint8_t(7));
+  EXPECT_EQ(digits[12], uint8_t(6));
+
+  EXPECT_EQ(hpd.getNumDigits(), 13u);
+  EXPECT_EQ(hpd.getDecimalPoint(), 13);
+
+  hpd.shift(-40); // shift right 40, equal to dividing by
+                  // 1099511627776. Result should be 1
+
+  EXPECT_EQ(digits[0], uint8_t(1));
+
+  EXPECT_EQ(hpd.getNumDigits(), 1u);
+  EXPECT_EQ(hpd.getDecimalPoint(), 1);
+}
+
+TEST(LlvmLibcHighPrecisionDecimalTest, VeryBigShift) {
+  __llvm_libc::internal::HighPrecsisionDecimal hpd =
+      __llvm_libc::internal::HighPrecsisionDecimal("1");
+  uint8_t *digits = hpd.getDigits();
+
+  hpd.shift(100); // shift left 100, equal to multiplying by
+                  // 1267650600228229401496703205376.
+  // result should be 2^100
+
+  EXPECT_EQ(digits[0], uint8_t(1));
+  EXPECT_EQ(digits[1], uint8_t(2));
+  EXPECT_EQ(digits[2], uint8_t(6));
+  EXPECT_EQ(digits[3], uint8_t(7));
+  EXPECT_EQ(digits[4], uint8_t(6));
+  EXPECT_EQ(digits[5], uint8_t(5));
+  EXPECT_EQ(digits[6], uint8_t(0));
+  EXPECT_EQ(digits[7], uint8_t(6));
+  EXPECT_EQ(digits[8], uint8_t(0));
+  EXPECT_EQ(digits[9], uint8_t(0));
+  EXPECT_EQ(digits[10], uint8_t(2));
+  EXPECT_EQ(digits[11], uint8_t(2));
+  EXPECT_EQ(digits[12], uint8_t(8));
+  EXPECT_EQ(digits[13], uint8_t(2));
+  EXPECT_EQ(digits[14], uint8_t(2));
+  EXPECT_EQ(digits[15], uint8_t(9));
+  EXPECT_EQ(digits[16], uint8_t(4));
+  EXPECT_EQ(digits[17], uint8_t(0));
+  EXPECT_EQ(digits[18], uint8_t(1));
+  EXPECT_EQ(digits[19], uint8_t(4));
+  EXPECT_EQ(digits[20], uint8_t(9));
+  EXPECT_EQ(digits[21], uint8_t(6));
+  EXPECT_EQ(digits[22], uint8_t(7));
+  EXPECT_EQ(digits[23], uint8_t(0));
+  EXPECT_EQ(digits[24], uint8_t(3));
+  EXPECT_EQ(digits[25], uint8_t(2));
+  EXPECT_EQ(digits[26], uint8_t(0));
+  EXPECT_EQ(digits[27], uint8_t(5));
+  EXPECT_EQ(digits[28], uint8_t(3));
+  EXPECT_EQ(digits[29], uint8_t(7));
+  EXPECT_EQ(digits[30], uint8_t(6));
+
+  EXPECT_EQ(hpd.getNumDigits(), 31u);
+  EXPECT_EQ(hpd.getDecimalPoint(), 31);
+
+  hpd.shift(-100); // shift right 100, equal to dividing by
+                   // 1267650600228229401496703205376.
+  // result should be 1
+
+  EXPECT_EQ(digits[0], uint8_t(1));
+  EXPECT_EQ(hpd.getNumDigits(), 1u);
+  EXPECT_EQ(hpd.getDecimalPoint(), 1);
+}
+
+TEST(LlvmLibcHighPrecisionDecimalTest, RoundingTest) {
+  __llvm_libc::internal::HighPrecsisionDecimal hpd =
+      __llvm_libc::internal::HighPrecsisionDecimal("1.2345");
+
+  EXPECT_EQ(hpd.roundToIntegerType<uint32_t>(), uint32_t(1));
+  EXPECT_EQ(hpd.roundToIntegerType<uint64_t>(), uint64_t(1));
+  EXPECT_EQ(hpd.roundToIntegerType<__uint128_t>(), __uint128_t(1));
+
+  hpd.shift(1); // shift left 1 to get 2.469 (rounds to 2)
+
+  EXPECT_EQ(hpd.roundToIntegerType<uint32_t>(), uint32_t(2));
+  EXPECT_EQ(hpd.roundToIntegerType<uint64_t>(), uint64_t(2));
+  EXPECT_EQ(hpd.roundToIntegerType<__uint128_t>(), __uint128_t(2));
+
+  hpd.shift(1); // shift left 1 to get 4.938 (rounds to 5)
+
+  EXPECT_EQ(hpd.roundToIntegerType<uint32_t>(), uint32_t(5));
+  EXPECT_EQ(hpd.roundToIntegerType<uint64_t>(), uint64_t(5));
+  EXPECT_EQ(hpd.roundToIntegerType<__uint128_t>(), __uint128_t(5));
+
+  // 2.5 is right between two integers, so we round to even (2)
+  hpd = __llvm_libc::internal::HighPrecsisionDecimal("2.5");
+
+  EXPECT_EQ(hpd.roundToIntegerType<uint32_t>(), uint32_t(2));
+  EXPECT_EQ(hpd.roundToIntegerType<uint64_t>(), uint64_t(2));
+  EXPECT_EQ(hpd.roundToIntegerType<__uint128_t>(), __uint128_t(2));
+
+  // unless it's marked as having truncated, which means it's actually slightly
+  // higher, forcing a round up (3)
+  hpd.setTruncated(true);
+
+  EXPECT_EQ(hpd.roundToIntegerType<uint32_t>(), uint32_t(3));
+  EXPECT_EQ(hpd.roundToIntegerType<uint64_t>(), uint64_t(3));
+  EXPECT_EQ(hpd.roundToIntegerType<__uint128_t>(), __uint128_t(3));
+
+  // Check that the larger int types are being handled properly (overflow is not
+  // handled, so int types that are too small are ignored for this test.)
+
+  // 1099511627776 = 2^40
+  hpd = __llvm_libc::internal::HighPrecsisionDecimal("1099511627776");
+
+  EXPECT_EQ(hpd.roundToIntegerType<uint64_t>(), uint64_t(1099511627776));
+  EXPECT_EQ(hpd.roundToIntegerType<__uint128_t>(), __uint128_t(1099511627776));
+
+  // 1267650600228229401496703205376 = 2^100
+  hpd = __llvm_libc::internal::HighPrecsisionDecimal(
+      "1267650600228229401496703205376");
+
+  __uint128_t result = __uint128_t(1) << 100;
+
+  EXPECT_EQ(hpd.roundToIntegerType<__uint128_t>(), result);
+}

diff --git a/libc/utils/mathtools/GenerateHPDConstants.py b/libc/utils/mathtools/GenerateHPDConstants.py
new file mode 100644 (file)
index 0000000..4c87d45
--- /dev/null
+++ b/libc/utils/mathtools/GenerateHPDConstants.py
@@ -0,0 +1,65 @@
+from math import *
+
+'''
+This script is used to generate a table used by 
+libc/src/__support/high_precision_decimal.h.
+
+For the ith entry in the table there are two values (indexed starting at 0).
+The first value is the number of digits longer the second value would be if
+multiplied by 2^i.
+The second value is the smallest number that would create that number of 
+additional digits (which in base ten is always 5^i). Anything less creates one 
+fewer digit.
+
+As an example, the 3rd entry in the table is {1, "125"}. This means that if 
+125 is multiplied by 2^3 = 8, it will have exactly one more digit.
+Multiplying it out we get 125 * 8 = 1000. 125 is the smallest number that gives
+that extra digit, for example 124 * 8 = 992, and all larger 3 digit numbers
+also give only one extra digit when multiplied by 8, for example 8 * 999 = 7992.
+This makes sense because 5^3 * 2^3 = 10^3, the smallest 4 digit number.
+
+For numbers with more digits we can ignore the digits past what's in the second
+value, since the most significant digits determine how many extra digits there 
+will be. Looking at the previous example, if we have 1000, and we look at just 
+the first 3 digits (since 125 has 3 digits), we see that 100 < 125, so we get 
+one fewer than 1 extra digits, which is 0. 
+Multiplying it out we get 1000 * 8 = 8000, which fits the expectation. 
+Another few quick examples: 
+For 1255, 125 !< 125, so 1 digit more: 1255 * 8 = 10040
+For 9999, 999 !< 125, so 1 digit more: 9999 * 8 = 79992
+
+Now let's try an example with the 10th entry: {4, "9765625"}. This one means 
+that 9765625 * 2^10 will have 4 extra digits. 
+Let's skip straight to the examples:
+For 1, 1 < 9765625, so 4-1=3 extra digits: 1 * 2^10 = 1024, 1 digit to 4 digits is a difference of 3.
+For 9765624, 9765624 < 9765625 so 3 extra digits: 9765624 * 1024 = 9999998976, 7 digits to 10 digits is a difference of 3.
+For 9765625, 9765625 !< 9765625 so 4 extra digits: 9765625 * 1024 = 10000000000, 7 digits to 11 digits is a difference of 4.
+For 9999999, 9999999 !< 9765625 so 4 extra digits: 9999999 * 1024 = 10239998976, 7 digits to 11 digits is a difference of 4.
+For 12345678, 1234567 < 9765625 so 3 extra digits: 12345678 * 1024 = 12641974272, 8 digits to 11 digits is a difference of 3.
+
+
+This is important when left shifting in the HPD class because it reads and
+writes the number backwards, and to shift in place it needs to know where the
+last digit should go. Since a binary left shift by i bits is the same as
+multiplying by 2^i we know that looking up the ith element in the table will
+tell us the number of additional digits. If the first digits of the number
+being shifted are greater than or equal to the digits of 5^i (the second value
+of each entry) then it is just the first value in the entry, else it is one
+fewer.
+'''
+
+
+# Generate Left Shift Table
+outStr = ""
+for i in range(61):
+  tenToTheI = 10**i
+  fiveToTheI = 5**i
+  outStr += "{"
+  # The number of new digits that would be created by multiplying 5**i by 2**i
+  outStr += str(ceil(log10(tenToTheI) - log10(fiveToTheI)))
+  outStr += ', "'
+  if not i == 0:
+    outStr += str(fiveToTheI)
+  outStr += '"},\n'
+
+print(outStr)