diff --git a/libc/src/__support/CMakeLists.txt b/libc/src/__support/CMakeLists.txt --- a/libc/src/__support/CMakeLists.txt +++ b/libc/src/__support/CMakeLists.txt @@ -12,12 +12,20 @@ ctype_utils.h ) +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 --- /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 + +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(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(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(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(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(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(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(shiftAmount) > MAX_SHIFT_AMOUNT) { + this->leftShift(MAX_SHIFT_AMOUNT); + shiftAmount -= MAX_SHIFT_AMOUNT; + } + this->leftShift(shiftAmount); + } + // Right + else { + while (static_cast(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 T roundToIntegerType() { + T result = 0; + uint32_t curDigit = 0; + + while (static_cast(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(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 --- a/libc/test/src/__support/CMakeLists.txt +++ b/libc/test/src/__support/CMakeLists.txt @@ -9,3 +9,13 @@ 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 --- /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(1)); + EXPECT_EQ(hpd.roundToIntegerType(), 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(2)); + EXPECT_EQ(hpd.roundToIntegerType(), 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(5)); + EXPECT_EQ(hpd.roundToIntegerType(), 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(2)); + EXPECT_EQ(hpd.roundToIntegerType(), 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(3)); + EXPECT_EQ(hpd.roundToIntegerType(), 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(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 --- /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)