From ccfcc0ff1f27dfebfda93d0e79e10e4261cd850d Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Thu, 1 Jun 2023 14:23:36 +0200 Subject: [PATCH] port ryu_generic_128 --- .../charconv/detail/from_chars_float_impl.hpp | 4 + .../boost/charconv/detail/ryu/generic_128.hpp | 29 +- .../charconv/detail/ryu/ryu_generic_128.hpp | 408 ++++++++++++++++++ 3 files changed, 431 insertions(+), 10 deletions(-) create mode 100644 include/boost/charconv/detail/ryu/ryu_generic_128.hpp diff --git a/include/boost/charconv/detail/from_chars_float_impl.hpp b/include/boost/charconv/detail/from_chars_float_impl.hpp index d1e2d02..053e96b 100644 --- a/include/boost/charconv/detail/from_chars_float_impl.hpp +++ b/include/boost/charconv/detail/from_chars_float_impl.hpp @@ -16,6 +16,10 @@ #include #include +#ifdef BOOST_CHARCONV_HAS_INT128 +# include +#endif + namespace boost { namespace charconv { namespace detail { #ifdef BOOST_MSVC diff --git a/include/boost/charconv/detail/ryu/generic_128.hpp b/include/boost/charconv/detail/ryu/generic_128.hpp index ee0d91f..80a3ec0 100644 --- a/include/boost/charconv/detail/ryu/generic_128.hpp +++ b/include/boost/charconv/detail/ryu/generic_128.hpp @@ -9,6 +9,10 @@ #include #include +#define BOOST_CHARCONV_POW5_TABLE_SIZE 56 +#define BOOST_CHARCONV_POW5_BITCOUNT 249 +#define BOOST_CHARCONV_POW5_INV_BITCOUNT 249 + namespace boost { namespace charconv { namespace detail { namespace ryu { // These tables are ~4.5 kByte total, compared to ~160 kByte for the full tables. @@ -16,11 +20,7 @@ namespace boost { namespace charconv { namespace detail { namespace ryu { // There's no way to define 128-bit constants in C, so we use little-endian // pairs of 64-bit constants. -static constexpr size_t pow5_table_size = 56; -static constexpr uint32_t pow5_bitcount = 249; -static constexpr uint32_t pow5_inv_bitcount = 249; - -static constexpr uint64_t GENERIC_POW5_TABLE[pow5_table_size][2] = { +static constexpr uint64_t GENERIC_POW5_TABLE[BOOST_CHARCONV_POW5_TABLE_SIZE][2] = { { 1u, 0u }, { 5u, 0u }, { 25u, 0u }, @@ -350,7 +350,7 @@ static constexpr uint64_t POW5_INV_ERRORS[154] = { }; // Returns e == 0 ? 1 : ceil(log_2(5^e)); requires 0 <= e <= 32768. -static constexpr uint32_t pow5bits(const uint32_t e) noexcept +static BOOST_CHARCONV_CXX14_CONSTEXPR uint32_t pow5bits(const uint32_t e) noexcept { BOOST_CHARCONV_ASSERT(e <= 1 << 15); return static_cast(((e * UINT64_C(163391164108059)) >> 46) + 1); @@ -421,8 +421,8 @@ void mul_128_256_shift( // Computes 5^i in the form required by Ryu, and stores it in the given pointer. static BOOST_CXX14_CONSTEXPR void generic_computePow5(const uint32_t i, uint64_t* const result) noexcept { - const uint32_t base = i / pow5_table_size; - const uint32_t base2 = base * pow5_table_size; + const uint32_t base = i / BOOST_CHARCONV_POW5_TABLE_SIZE; + const uint32_t base2 = base * BOOST_CHARCONV_POW5_TABLE_SIZE; const uint64_t* const mul = GENERIC_POW5_SPLIT[base]; if (i == base2) { @@ -444,8 +444,8 @@ static BOOST_CXX14_CONSTEXPR void generic_computePow5(const uint32_t i, uint64_t // Computes 5^-i in the form required by Ryu, and stores it in the given pointer. static BOOST_CXX14_CONSTEXPR void generic_computeInvPow5(const uint32_t i, uint64_t* const result) noexcept { - const uint32_t base = (i + pow5_table_size - 1) / pow5_table_size; - const uint32_t base2 = base * pow5_table_size; + const uint32_t base = (i + BOOST_CHARCONV_POW5_TABLE_SIZE - 1) / BOOST_CHARCONV_POW5_TABLE_SIZE; + const uint32_t base2 = base * BOOST_CHARCONV_POW5_TABLE_SIZE; const uint64_t* const mul = GENERIC_POW5_INV_SPLIT[base]; // 1/5^base2 if (i == base2) { @@ -526,6 +526,15 @@ static BOOST_CHARCONV_CXX14_CONSTEXPR uint32_t log10Pow2(const int32_t e) noexce return (uint32_t) ((((uint64_t) e) * UINT64_C(169464822037455)) >> 49); } +// Returns floor(log_10(5^e)). +static BOOST_CHARCONV_CXX14_CONSTEXPR uint32_t log10Pow5(const int32_t e) noexcept +{ + // The first value this approximation fails for is 5^2621 which is just greater than 10^1832. + assert(e >= 0); + assert(e <= 1 << 15); + return (uint32_t) ((((uint64_t) e) * UINT64_C(196742565691928)) >> 48); +} + }}}} // Namespaces #endif // BOOST_CHARCONV_DETAIL_RYU_GENERIC_128_HPP diff --git a/include/boost/charconv/detail/ryu/ryu_generic_128.hpp b/include/boost/charconv/detail/ryu/ryu_generic_128.hpp new file mode 100644 index 0000000..0abcc90 --- /dev/null +++ b/include/boost/charconv/detail/ryu/ryu_generic_128.hpp @@ -0,0 +1,408 @@ +// Copyright 2018 - 2023 Ulf Adams +// Copyright 2023 Matt Borland +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt + +#ifndef BOOST_CHARCONV_DETAIL_RYU_RYU_GENERIC_128_HPP +#define BOOST_CHARCONV_DETAIL_RYU_RYU_GENERIC_128_HPP// + +#include +#include +#include +#include +#include +#include + +namespace boost { namespace charconv { namespace detail { namespace ryu { + +static constexpr int32_t fd128_exceptional_exponent = 0x7FFFFFFF; +static constexpr boost::uint128_type one = 1; + +struct floating_decimal_128 +{ + boost::uint128_type mantissa; + int32_t exponent; + bool sign; +}; + +#ifdef BOOST_CHARCONV_DEBUG +static char* s(boost::uint128_type v) { + int len = decimalLength(v); + char* b = (char*) malloc((len + 1) * sizeof(char)); + for (int i = 0; i < len; i++) { + const uint32_t c = (uint32_t) (v % 10); + v /= 10; + b[len - 1 - i] = (char) ('0' + c); + } + b[len] = 0; + return b; +} +#endif + +struct floating_decimal_128 generic_binary_to_decimal( + const boost::uint128_type bits, + const uint32_t mantissaBits, const uint32_t exponentBits, const bool explicitLeadingBit) noexcept +{ + #ifdef BOOST_CHARCONV_DEBUG + printf("IN="); + for (int32_t bit = 127; bit >= 0; --bit) + { + printf("%u", (uint32_t) ((bits >> bit) & 1)); + } + printf("\n"); + #endif + + const uint32_t bias = (1u << (exponentBits - 1)) - 1; + const bool ieeeSign = ((bits >> (mantissaBits + exponentBits)) & 1) != 0; + const boost::uint128_type ieeeMantissa = bits & ((one << mantissaBits) - 1); + const uint32_t ieeeExponent = (uint32_t) ((bits >> mantissaBits) & ((one << exponentBits) - 1u)); + + if (ieeeExponent == 0 && ieeeMantissa == 0) + { + struct floating_decimal_128 fd {0, 0, ieeeSign}; + return fd; + } + if (ieeeExponent == ((1u << exponentBits) - 1u)) + { + struct floating_decimal_128 fd; + fd.mantissa = explicitLeadingBit ? ieeeMantissa & ((one << (mantissaBits - 1)) - 1) : ieeeMantissa; + fd.exponent = fd128_exceptional_exponent; + fd.sign = ieeeSign; + return fd; + } + + int32_t e2; + boost::uint128_type m2; + // We subtract 2 in all cases so that the bounds computation has 2 additional bits. + if (explicitLeadingBit) + { + // mantissaBits includes the explicit leading bit, so we need to correct for that here. + if (ieeeExponent == 0) + { + e2 = (int32_t)(1 - bias - mantissaBits + 1 - 2); + } + else + { + e2 = (int32_t)(ieeeExponent - bias - mantissaBits + 1 - 2); + } + m2 = ieeeMantissa; + } + else + { + if (ieeeExponent == 0) + { + e2 = (int32_t)(1 - bias - mantissaBits - 2); + m2 = ieeeMantissa; + } else + { + e2 = ieeeExponent - bias - mantissaBits - 2; + m2 = (one << mantissaBits) | ieeeMantissa; + } + } + const bool even = (m2 & 1) == 0; + const bool acceptBounds = even; + + #ifdef BOOST_CHARCONV_DEBUG + printf("-> %s %s * 2^%d\n", ieeeSign ? "-" : "+", s(m2), e2 + 2); + #endif + + // Step 2: Determine the interval of legal decimal representations. + const boost::uint128_type mv = 4 * m2; + // Implicit bool -> int conversion. True is 1, false is 0. + const uint32_t mmShift = + (ieeeMantissa != (explicitLeadingBit ? one << (mantissaBits - 1) : 0)) + || (ieeeExponent == 0); + + // Step 3: Convert to a decimal power base using 128-bit arithmetic. + boost::uint128_type vr; + boost::uint128_type vp; + boost::uint128_type vm; + int32_t e10; + bool vmIsTrailingZeros = false; + bool vrIsTrailingZeros = false; + if (e2 >= 0) + { + // I tried special-casing q == 0, but there was no effect on performance. + // This expression is slightly faster than max(0, log10Pow2(e2) - 1). + const uint32_t q = log10Pow2(e2) - (e2 > 3); + e10 = (int32_t)q; + const int32_t k = BOOST_CHARCONV_POW5_INV_BITCOUNT + pow5bits(q) - 1; + const int32_t i = (int32_t)(-e2 + q + k); + uint64_t pow5[4]; + generic_computeInvPow5(q, pow5); + vr = mulShift(4 * m2, pow5, i); + vp = mulShift(4 * m2 + 2, pow5, i); + vm = mulShift(4 * m2 - 1 - mmShift, pow5, i); +#ifdef BOOST_CHARCONV_DEBUG + printf("%s * 2^%d / 10^%d\n", s(mv), e2, q); + printf("V+=%s\nV =%s\nV-=%s\n", s(vp), s(vr), s(vm)); +#endif + // floor(log_5(2^128)) = 55, this is very conservative + if (q <= 55) + { + // Only one of mp, mv, and mm can be a multiple of 5, if any. + if (mv % 5 == 0) + { + vrIsTrailingZeros = multipleOfPowerOf5(mv, q - 1); + } + else if (acceptBounds) + { + // Same as min(e2 + (~mm & 1), pow5Factor(mm)) >= q + // <=> e2 + (~mm & 1) >= q && pow5Factor(mm) >= q + // <=> true && pow5Factor(mm) >= q, since e2 >= q. + vmIsTrailingZeros = multipleOfPowerOf5(mv - 1 - mmShift, q); + } + else + { + // Same as min(e2 + 1, pow5Factor(mp)) >= q. + vp -= multipleOfPowerOf5(mv + 2, q); + } + } + } + else + { + // This expression is slightly faster than max(0, log10Pow5(-e2) - 1). + const uint32_t q = log10Pow5(-e2) - (int32_t)(-e2 > 1); + e10 = (int32_t)q + e2; + const int32_t i = (int32_t)(-e2 - q); + const int32_t k = (int32_t)pow5bits(i) - BOOST_CHARCONV_POW5_BITCOUNT; + const int32_t j = (int32_t)q - k; + uint64_t pow5[4]; + generic_computePow5(i, pow5); + vr = mulShift(4 * m2, pow5, j); + vp = mulShift(4 * m2 + 2, pow5, j); + vm = mulShift(4 * m2 - 1 - mmShift, pow5, j); +#ifdef BOOST_CHARCONV_DEBUG + printf("%s * 5^%d / 10^%d\n", s(mv), -e2, q); + printf("%d %d %d %d\n", q, i, k, j); + printf("V+=%s\nV =%s\nV-=%s\n", s(vp), s(vr), s(vm)); +#endif + if (q <= 1) + { + // {vr,vp,vm} is trailing zeros if {mv,mp,mm} has at least q trailing 0 bits. + // mv = 4 m2, so it always has at least two trailing 0 bits. + vrIsTrailingZeros = true; + if (acceptBounds) + { + // mm = mv - 1 - mmShift, so it has 1 trailing 0 bit iff mmShift == 1. + vmIsTrailingZeros = mmShift == 1; + } + else + { + // mp = mv + 2, so it always has at least one trailing 0 bit. + --vp; + } + } + else if (q < 127) + { + // We need to compute min(ntz(mv), pow5Factor(mv) - e2) >= q-1 + // <=> ntz(mv) >= q-1 && pow5Factor(mv) - e2 >= q-1 + // <=> ntz(mv) >= q-1 (e2 is negative and -e2 >= q) + // <=> (mv & ((1 << (q-1)) - 1)) == 0 + // We also need to make sure that the left shift does not overflow. + vrIsTrailingZeros = multipleOfPowerOf2(mv, q - 1); + + #ifdef BOOST_CHARCONV_DEBUG + printf("vr is trailing zeros=%s\n", vrIsTrailingZeros ? "true" : "false"); + #endif + } + } + + #ifdef BOOST_CHARCONV_DEBUG + printf("e10=%d\n", e10); + printf("V+=%s\nV =%s\nV-=%s\n", s(vp), s(vr), s(vm)); + printf("vm is trailing zeros=%s\n", vmIsTrailingZeros ? "true" : "false"); + printf("vr is trailing zeros=%s\n", vrIsTrailingZeros ? "true" : "false"); + #endif + + // Step 4: Find the shortest decimal representation in the interval of legal representations. + uint32_t removed = 0; + uint8_t lastRemovedDigit = 0; + boost::uint128_type output; + + while (vp / 10 > vm / 10) + { + vmIsTrailingZeros &= vm % 10 == 0; + vrIsTrailingZeros &= lastRemovedDigit == 0; + lastRemovedDigit = (uint8_t) (vr % 10); + vr /= 10; + vp /= 10; + vm /= 10; + ++removed; + } + + #ifdef BOOST_CHARCONV_DEBUG + printf("V+=%s\nV =%s\nV-=%s\n", s(vp), s(vr), s(vm)); + printf("d-10=%s\n", vmIsTrailingZeros ? "true" : "false"); + #endif + + if (vmIsTrailingZeros) + { + while (vm % 10 == 0) + { + vrIsTrailingZeros &= lastRemovedDigit == 0; + lastRemovedDigit = (uint8_t) (vr % 10); + vr /= 10; + vp /= 10; + vm /= 10; + ++removed; + } + } + + #ifdef BOOST_CHARCONV_DEBUG + printf("%s %d\n", s(vr), lastRemovedDigit); + printf("vr is trailing zeros=%s\n", vrIsTrailingZeros ? "true" : "false"); + #endif + + if (vrIsTrailingZeros && (lastRemovedDigit == 5) && (vr % 2 == 0)) + { + // Round even if the exact numbers is .....50..0. + lastRemovedDigit = 4; + } + // We need to take vr+1 if vr is outside bounds or we need to round up. + output = vr + (boost::uint128_type)((vr == vm && (!acceptBounds || !vmIsTrailingZeros)) || (lastRemovedDigit >= 5)); + const int32_t exp = e10 + removed; + + #ifdef BOOST_CHARCONV_DEBUG + printf("V+=%s\nV =%s\nV-=%s\n", s(vp), s(vr), s(vm)); + printf("O=%s\n", s(output)); + printf("EXP=%d\n", exp); + #endif + + return {output, exp, ieeeSign}; +} + +static inline int copy_special_str(char * const result, const struct floating_decimal_128 fd) noexcept +{ + // TODO(mborland): match the handling of the other NaNs and infinities in the library + if (fd.mantissa) + { + memcpy(result, "NaN", 3); + return 3; + } + if (fd.sign) + { + result[0] = '-'; + } + + memcpy(result + fd.sign, "Infinity", 8); + return fd.sign + 8; +} + +// Converts the given decimal floating point number to a string, writing to result, and returning +// the number characters written. Does not terminate the buffer with a 0. In the worst case, this +// function can write up to 53 characters. +// +// Maximal char buffer requirement: +// sign + mantissa digits + decimal dot + 'E' + exponent sign + exponent digits +// = 1 + 39 + 1 + 1 + 1 + 10 = 53 +int generic_to_chars(const struct floating_decimal_128 v, char* const result) noexcept +{ + if (v.exponent == fd128_exceptional_exponent) + { + return copy_special_str(result, v); + } + + // Step 5: Print the decimal representation. + size_t index = 0; + if (v.sign) + { + result[index++] = '-'; + } + + boost::uint128_type output = v.mantissa; + const uint32_t olength = decimalLength(output); + + #ifdef BOOST_CHARCONV_DEBUG + printf("DIGITS=%s\n", s(v.mantissa)); + printf("OLEN=%u\n", olength); + printf("EXP=%u\n", v.exponent + olength); + #endif + + for (uint32_t i = 0; i < olength - 1; ++i) + { + const uint32_t c = (uint32_t) (output % 10); + output /= 10; + result[index + olength - i] = (char) ('0' + c); + } + result[index] = '0' + (uint32_t)(output % 10); // output should be < 10 by now. + + // Print decimal point if needed. + if (olength > 1) + { + result[index + 1] = '.'; + index += olength + 1; + } + else + { + ++index; + } + + // Print the exponent. + result[index++] = 'E'; + int32_t exp = v.exponent + olength - 1; + if (exp < 0) + { + result[index++] = '-'; + exp = -exp; + } + + uint32_t elength = decimalLength(exp); + for (uint32_t i = 0; i < elength; ++i) + { + const uint32_t c = exp % 10; + exp /= 10; + result[index + elength - 1 - i] = (char) ('0' + c); + } + index += elength; + return index; +} + +struct floating_decimal_128 float_to_fd128(float f) noexcept +{ + static_assert(sizeof(float) == sizeof(uint32_t), "Float is not 32 bits"); + uint32_t bits = 0; + memcpy(&bits, &f, sizeof(float)); + return generic_binary_to_decimal(bits, 23, 8, false); +} + +struct floating_decimal_128 double_to_fd128(double d) noexcept +{ + static_assert(sizeof(double) == sizeof(uint64_t), "Float is not 64 bits"); + uint64_t bits = 0; + memcpy(&bits, &d, sizeof(double)); + return generic_binary_to_decimal(bits, 52, 11, false); +} + +#if BOOST_CHARCONV_LDBL_BITS == 80 + +struct floating_decimal_128 long_double_to_fd128(long double d) noexcept +{ + boost::uint128_type bits = 0; + memcpy(&bits, &d, sizeof(long double)); + + #ifdef BOOST_CHARCONV_DEBUG + // For some odd reason, this ends up with noise in the top 48 bits. We can + // clear out those bits with the following line; this is not required, the + // conversion routine should ignore those bits, but the debug output can be + // confusing if they aren't 0s. + bits &= (one << 80) - 1; + #endif + + return generic_binary_to_decimal(bits, 64, 15, true); +} + +#elif BOOST_CHARCONV_LDBL_BITS == 128 + +struct floating_decimal_128 long_double_to_fd128(long double d) noexcept +{ + boost::uint128_type bits = 0; + memcpy(&bits, &d, sizeof(long double)); + return generic_binary_to_decimal(bits, 113, 15, true); +} + +#endif + +}}}} // Namespaces + +#endif //BOOST_RYU_GENERIC_128_HPP