Write doubles with the shortest digits (Żmij), digits in registers

Write doubles with the conversion of Zmij by Victor Zverovich (MIT),
ported to C++11 (detail/conversions/zmij.hpp). It finds the
shortest decimal that reads back as the same double, and the
closest one if there are several. Grisu2, used until now, is fast
but not always shortest: it sometimes writes a 17th digit where 16
suffice, or a last digit that is not the closest. The layout is
unchanged (1.5, 100.0, 1e+100, -0.0); float keeps Grisu2.

Digits are converted eight at a time with the BCD conversion of
Xiang JunBo, as in Zmij, and written with one byte swap per eight
digits and fixed-size moves instead of per-digit loops. Leading and
trailing zeros are counted from those bytes. to_chars() uses a
local buffer when the caller's is shorter than the 41 bytes this
may write. The powers of ten come from the number-parsing table,
adjusted where it holds values rounded up, and extended with Zmij's
compressed tables beyond 10^308.

write_shortest() converts its 16 digits in one vector register
(SSE2 on x86-64, NEON on 64-bit Arm, both baseline) and inserts the
decimal point inside the register, avoiding a store-forwarding
stall that cost about 25% of the time to write a double. dump()
writes floats and integers straight into the serializer's write
buffer instead of copying them from a member buffer, and small
integers eight digits at a time. read_eight_bytes() and
parse_eight_digits() are marked always-inline, which GCC had been
calling out of line in the number-parsing loops.

Of one million random doubles, about 0.14% are now written with
different digits, always to a value that still reads back as the
same double.

Signed-off-by: Niels Lohmann <mail@nlohmann.me>
This commit is contained in:
Niels Lohmann committed 2026-10-06 10:48:06 +02:00
1 parent 953d74ddcd
commit 17842694f6
15 files changed
+1823 -107

No files matched your search

+11 -11
View File
@@ -38,7 +38,7 @@ TEST_CASE("Binary Formats" * doctest::skip())
const auto ubjson_2_size = json::to_ubjson(j, true).size();
const auto ubjson_3_size = json::to_ubjson(j, true, true).size();
CHECK(json_size == 2090303);
CHECK(json_size == 2090234);
CHECK(bjdata_1_size == 1112030);
CHECK(bjdata_2_size == 1224148);
CHECK(bjdata_3_size == 1224148);
@@ -51,16 +51,16 @@ TEST_CASE("Binary Formats" * doctest::skip())
CHECK(ubjson_3_size == 1169069);
CHECK((100.0 * double(json_size) / double(json_size)) == Approx(100.0));
CHECK((100.0 * double(bjdata_1_size) / double(json_size)) == Approx(53.199));
CHECK((100.0 * double(bjdata_2_size) / double(json_size)) == Approx(58.563));
CHECK((100.0 * double(bjdata_3_size) / double(json_size)) == Approx(58.563));
CHECK((100.0 * double(bon8_size) / double(json_size)) == Approx(50.509));
CHECK((100.0 * double(bson_size) / double(json_size)) == Approx(85.849));
CHECK((100.0 * double(cbor_size) / double(json_size)) == Approx(50.497));
CHECK((100.0 * double(msgpack_size) / double(json_size)) == Approx(50.526));
CHECK((100.0 * double(ubjson_1_size) / double(json_size)) == Approx(53.199));
CHECK((100.0 * double(ubjson_2_size) / double(json_size)) == Approx(58.563));
CHECK((100.0 * double(ubjson_3_size) / double(json_size)) == Approx(55.928));
CHECK((100.0 * double(bjdata_1_size) / double(json_size)) == Approx(53.201));
CHECK((100.0 * double(bjdata_2_size) / double(json_size)) == Approx(58.565));
CHECK((100.0 * double(bjdata_3_size) / double(json_size)) == Approx(58.565));
CHECK((100.0 * double(bon8_size) / double(json_size)) == Approx(50.511));
CHECK((100.0 * double(bson_size) / double(json_size)) == Approx(85.853));
CHECK((100.0 * double(cbor_size) / double(json_size)) == Approx(50.499));
CHECK((100.0 * double(msgpack_size) / double(json_size)) == Approx(50.528));
CHECK((100.0 * double(ubjson_1_size) / double(json_size)) == Approx(53.201));
CHECK((100.0 * double(ubjson_2_size) / double(json_size)) == Approx(58.565));
CHECK((100.0 * double(ubjson_3_size) / double(json_size)) == Approx(55.930));
}
SECTION("twitter.json")
+261 -1
View File
@@ -15,6 +15,23 @@
#include <nlohmann/json.hpp>
using nlohmann::detail::dtoa_impl::reinterpret_bits;
#include <array>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <iomanip>
#include <limits>
#include <locale>
#include <random>
#include <sstream>
#include <string>
#include <utility>
#include <vector>
#if defined(JSON_HAS_CPP_17)
#include <charconv>
#endif
namespace
{
float make_float(uint32_t sign_bit, uint32_t biased_exponent, uint32_t significand)
@@ -450,7 +467,7 @@ TEST_CASE("formatting")
check_double( 1.2345e+18, "1.2345e+18" ); // 1.2345e+18 1.2345e+18 1.2345e18
check_double( 1.2345e+19, "1.2345e+19" ); // 1.2345e+19 1.2345e+19 1.2345e19
check_double( 1.2345e+20, "1.2345e+20" ); // 1.2345e+20 1.2345e+20 1.2345e20
check_double( 1.2345e+21, "1.2344999999999999e+21" ); // 1.2345e+21 1.2344999999999999e+21 1.2345e21
check_double( 1.2345e+21, "1.2345e+21" ); // 1.2345e+21 1.2344999999999999e+21 1.2345e21
check_double( 1.2345e+22, "1.2345e+22" ); // 1.2345e+22 1.2345e+22 1.2345e22
}
@@ -514,3 +531,246 @@ TEST_CASE("formatting")
check_integer(1000000000000000000LL, "1000000000000000000");
}
}
namespace
{
// a small unsigned big integer (32-bit limbs, least significant first), to
// recompute the powers of ten of the shortest double conversion
using big = std::vector<std::uint32_t>;
void big_mul_small(big& x, std::uint32_t m)
{
std::uint64_t carry = 0;
for (auto& limb : x)
{
const std::uint64_t v = (static_cast<std::uint64_t>(limb) * m) + carry;
limb = static_cast<std::uint32_t>(v);
carry = v >> 32u;
}
if (carry != 0)
{
x.push_back(static_cast<std::uint32_t>(carry));
}
}
void big_div_small(big& x, std::uint32_t d)
{
std::uint64_t rest = 0;
for (std::size_t i = x.size(); i-- > 0;)
{
const std::uint64_t v = (rest << 32u) | x[i];
x[i] = static_cast<std::uint32_t>(v / d);
rest = v % d;
}
while (!x.empty() && x.back() == 0)
{
x.pop_back();
}
}
std::size_t big_bit_length(const big& x)
{
std::size_t n = 32 * x.size();
for (std::uint32_t top = x.back(); (top & 0x80000000u) == 0; top <<= 1u)
{
--n;
}
return n;
}
bool big_bit(const big& x, std::size_t i)
{
return ((x[i / 32] >> (i % 32)) & 1u) != 0;
}
/// the 128 most significant bits of x (floor), shifted left if x has fewer bits
std::pair<std::uint64_t, std::uint64_t> big_top128(const big& x)
{
const std::size_t n = big_bit_length(x);
std::uint64_t high = 0;
std::uint64_t low = 0;
for (std::size_t k = 0; k < 128; ++k)
{
const bool bit = k < n && big_bit(x, n - 1 - k);
if (k < 64)
{
high = (high << 1u) | (bit ? 1u : 0u);
}
else
{
low = (low << 1u) | (bit ? 1u : 0u);
}
}
return {high, low};
}
/// the digits (without trailing zeros) and the decimal exponent of a
/// representation "[-]d[.ddd][e[+-]x]"
std::pair<std::string, int> digits_and_exponent(const std::string& s)
{
std::string digits;
int point = -1;
int exponent = 0;
for (std::size_t i = 0; i < s.size(); ++i)
{
const char c = s[i];
if (c >= '0' && c <= '9')
{
digits += c;
}
else if (c == '.')
{
point = static_cast<int>(digits.size());
}
else if (c == 'e' || c == 'E')
{
exponent = std::stoi(s.substr(i + 1));
break;
}
}
int e = exponent + (point < 0 ? static_cast<int>(digits.size()) : point) - static_cast<int>(digits.size());
const std::size_t first = digits.find_first_not_of('0');
digits = first == std::string::npos ? "0" : digits.substr(first);
while (digits.size() > 1 && digits.back() == '0')
{
digits.pop_back();
++e;
}
return {digits, e};
}
/// whether the decimal digits * 10^e reads back as v
bool reads_back(const std::string& digits, int e, double v)
{
const std::string text = digits + "e" + std::to_string(e);
// (compared bit for bit: v is positive and finite, and -Wfloat-equal)
return reinterpret_bits<std::uint64_t>(std::strtod(text.c_str(), nullptr)) == reinterpret_bits<std::uint64_t>(v);
}
/// Check the representation of a positive finite double: it reads back as
/// the same value, and no representation with fewer digits does.
void check_shortest(double v)
{
std::array<char, 33> buf{};
char* end = nlohmann::detail::to_chars(buf.data(), buf.data() + 32, v);
const std::string text(buf.data(), end);
CAPTURE(text)
CHECK(std::strtod(text.c_str(), nullptr) == v);
// the layout is that of format_buffer() for the same digits
std::array<char, 64> reference{};
int len = 0;
int exponent = 0;
nlohmann::detail::dtoa_impl::shortest_digits(reference.data(), len, exponent, v);
const char* const reference_end = nlohmann::detail::dtoa_impl::format_buffer(reference.data(), len, exponent, -4, 15);
CHECK(text == std::string(reference.data(), static_cast<std::size_t>(reference_end - reference.data())));
const auto de = digits_and_exponent(text);
const std::string& digits = de.first;
if (digits.size() > 1)
{
// the decimals of one digit fewer next to the value
// (a stream rather than snprintf("%.*e"), whose output GCC cannot bound)
std::ostringstream shorter;
shorter.imbue(std::locale::classic());
shorter << std::scientific << std::setprecision(static_cast<int>(digits.size()) - 2) << v;
const auto near = digits_and_exponent(shorter.str());
// as an integer with digits.size() - 1 digits
std::string m = near.first;
int e = near.second;
while (m.size() < digits.size() - 1)
{
m += '0';
--e;
}
const std::uint64_t mid = std::stoull(m);
for (const std::uint64_t candidate :
{
mid - 1, mid, mid + 1
})
{
CAPTURE(candidate)
CHECK(!reads_back(std::to_string(candidate), e, v));
}
}
#if defined(JSON_HAS_CPP_17) && defined(__cpp_lib_to_chars)
// the closest of the shortest representations, as std::to_chars finds it
std::array<char, 64> std_text{};
const auto r = std::to_chars(std_text.data(), std_text.data() + std_text.size(), v, std::chars_format::scientific);
CHECK(digits_and_exponent(std::string(std_text.data(), r.ptr)) == de);
#endif
}
} // namespace
TEST_CASE("shortest digits of doubles")
{
SECTION("powers of ten")
{
// the 128-bit significands of 10^k, rounded down, recomputed
for (int k = -342; k <= 341; ++k)
{
CAPTURE(k)
big x{1};
if (k >= 0)
{
for (int i = 0; i < k; ++i)
{
big_mul_small(x, 10);
}
}
else
{
// floor(2^b / 10^-k) for a b that leaves more than 128 bits
const int b = 128 + 64 + (4 * -k);
x.assign(static_cast<std::size_t>(b / 32) + 1, 0);
x.back() = 1u << (b % 32);
for (int i = 0; i < -k; ++i)
{
big_div_small(x, 10);
}
}
const auto expected = big_top128(x);
const auto actual = nlohmann::detail::zmij::pow10(k);
CHECK(actual.high == expected.first);
CHECK(actual.low == expected.second);
}
}
SECTION("boundary values")
{
for (const double v :
{
std::numeric_limits<double>::min(), std::numeric_limits<double>::max(), std::numeric_limits<double>::denorm_min(),
std::nextafter(std::numeric_limits<double>::min(), 0.0), 1.0, 2.0, 0.1, 0.3, 1e21, 1e22, 1e23, 5e-324, 9007199254740993.0,
1.2345e+21, 2.2250738585072014e-308, 1.7976931348623157e308, 4.9406564584124654e-324, 123456789012345680.0
})
{
check_shortest(v);
}
// all powers of two (their rounding interval is narrower below)
for (int e = -1074; e <= 1023; ++e)
{
check_shortest(std::ldexp(1.0, e));
}
// powers of ten and their neighbors
for (int e = -323; e <= 308; ++e)
{
const double p = std::strtod(("1e" + std::to_string(e)).c_str(), nullptr);
check_shortest(p);
check_shortest(std::nextafter(p, 0.0));
check_shortest(std::nextafter(p, std::numeric_limits<double>::infinity()));
}
}
SECTION("random doubles")
{
std::mt19937_64 rng(5295); // NOLINT(cert-msc32-c,cert-msc51-cpp,bugprone-random-generator-seed): reproducible
for (int i = 0; i < 100000; ++i)
{
const std::uint64_t bits = rng() & 0x7FFFFFFFFFFFFFFFu;
const auto v = reinterpret_bits<double>(bits);
if (std::isfinite(v) && bits != 0)
{
check_shortest(v);
}
}
}
}