boost/charconv
2
0
Fork 0
mirror of https://github.com/boostorg/charconv.git synced 2026-02-03 09:02:11 +00:00
Code Issues Projects Releases Wiki Activity

port ryu_generic_128

Browse Source
This commit is contained in:
Matt Borland
2023-06-01 14:23:36 +02:00
parent 77655ab9fa
commit ccfcc0ff1f
3 changed files with 431 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
include/boost/charconv/detail/from_chars_float_impl.hpp
View File

@@ -16,6 +16,10 @@
#include <cstdlib>
#include <cmath>
#ifdef BOOST_CHARCONV_HAS_INT128
# include <boost/charconv/detail/ryu/ryu_generic_128.hpp>
#endif
namespace boost { namespace charconv { namespace detail {
#ifdef BOOST_MSVC

29
include/boost/charconv/detail/ryu/generic_128.hpp
View File

@@ -9,6 +9,10 @@
#include <boost/charconv/detail/config.hpp>
#include <cstdint>
#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<uint32_t>(((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

408
include/boost/charconv/detail/ryu/ryu_generic_128.hpp Normal file
View File

@@ -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 <boost/charconv/detail/ryu/generic_128.hpp>
#include <boost/charconv/detail/config.hpp>
#include <boost/charconv/detail/bit_layouts.hpp>
#include <cinttypes>
#include <cstdio>
#include <cstdint>
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