Write doubles with the shortest digits (Zmij)

dump() writes doubles with the conversion of Zmij by Victor Zverovich (https://github.com/vitaut/zmij, MIT), ported to C++11 in detail/conversions/zmij.hpp: the shortest decimal in the rounding interval, the closest one if there are several. Grisu2 does not always find the shortest digits; about 0.14% of random doubles are now written differently (0.08% with fewer digits, 0.06% with the closest last digit); short decimals such as 0.1 or 2555.56 are not affected. float keeps Grisu2. The layout of doubles is unchanged, but written differently: the digits are converted eight at a time (the BCD conversion of Xiang JunBo, as in Zmij) and stored with one byte swap per eight digits; leading and trailing zeros are counted from those bytes; and the layouts of format_buffer() are written with fixed-size moves instead of per-digit loops and moves of the buffer (to_chars() uses a local buffer if the caller's is shorter than the 41 bytes this may write). The powers of ten come from the table for number parsing, adjusted where it holds them rounded up, and from the compressed tables of Zmij beyond 10^308. json::dump() gets faster on floats: canada -53%, numbers -46%, mesh -37%, marine_ik -30%. Tests: the powers of ten recomputed with a small big-integer; for random doubles, all powers of two and of ten and their neighbors, and boundary values: the output reads back as the same value, no decimal with one digit fewer does, the layout equals that of format_buffer() for the same digits, and (C++17) the digits equal those of std::to_chars. The size ratios of canada.json in unit-binary_formats.cpp and one expectation in unit-to_chars.cpp change with the shorter output. Signed-off-by: Niels Lohmann <mail@nlohmann.me>
2026-09-29 19:20:30 +00:00 · 2026-09-29 14:03:25 +02:00
11 changed files with 1243 additions and 65 deletions
@@ -26,6 +26,7 @@ cc_library(
        "include/nlohmann/detail/conversions/from_json.hpp",
        "include/nlohmann/detail/conversions/to_chars.hpp",
        "include/nlohmann/detail/conversions/to_json.hpp",
+        "include/nlohmann/detail/conversions/zmij.hpp",
        "include/nlohmann/detail/exceptions.hpp",
        "include/nlohmann/detail/hash.hpp",
        "include/nlohmann/detail/input/binary_reader.hpp",
@@ -1401,6 +1401,7 @@ THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR I

 - The class contains the UTF-8 Decoder from Bjoern Hoehrmann which is licensed under the [MIT License](https://opensource.org/licenses/MIT) (see above). Copyright &copy; 2008-2009 [Björn Hoehrmann](https://bjoern.hoehrmann.de/) <bjoern@hoehrmann.de>
 - The class contains a slightly modified version of the Grisu2 algorithm from Florian Loitsch which is licensed under the [MIT License](https://opensource.org/licenses/MIT) (see above). Copyright &copy; 2009 [Florian Loitsch](https://florian.loitsch.com/)
+- The class contains a port of the shortest double-to-decimal conversion of [Żmij](https://github.com/vitaut/zmij) by Victor Zverovich, which is licensed under the [MIT License](https://opensource.org/licenses/MIT) (see above). Copyright &copy; 2025 [Victor Zverovich](https://github.com/vitaut)
 - The class contains a copy of [Hedley](https://nemequ.github.io/hedley/) from Evan Nemerson which is licensed as [CC0-1.0](https://creativecommons.org/publicdomain/zero/1.0/).
 - The class contains parts of [Google Abseil](https://github.com/abseil/abseil-cpp) which is licensed under the [Apache 2.0 License](https://opensource.org/licenses/Apache-2.0).
 - The class contains an adapted version of the Eisel-Lemire algorithm and its table of powers of five from [fast_float](https://github.com/fastfloat/fast_float) by Daniel Lemire and contributors, which is available under the [MIT License](https://opensource.org/licenses/MIT) (used here), the Apache 2.0 License, and the Boost Software License. Copyright &copy; 2021 The fast_float authors
@@ -60,6 +60,9 @@ Linear.

 ## Notes

+Floating-point numbers are written with the fewest digits that read back as the same value (for `#!cpp double`; see
+[number handling](../../features/types/number_handling.md#number-serialization)).
+
 Binary values are serialized as an object containing two keys:

 - "bytes": an array of bytes as integers
@@ -96,3 +99,5 @@ Binary values are serialized as an object containing two keys:
 - Indentation character `indent_char`, option `ensure_ascii` and exceptions added in version 3.0.0.
 - Error handlers added in version 3.4.0.
 - Serialization of binary values added in version 3.8.0.
+- Doubles are written with the shortest digits (Żmij instead of Grisu2) since version 3.13.0; about 0.1% of doubles are
+  written differently, most of them with fewer digits.
@@ -118,9 +118,10 @@ That is, `-0` is stored as a signed integer, but the serialization does not repr
 ### Number serialization

 - Integer numbers are serialized as is; that is, no scientific notation is used.
- Floating-point numbers are serialized as specified by the `#!c %g` printf modifier with 
-  [`std::numeric_limits<double>::max_digits10`](https://en.cppreference.com/w/cpp/types/numeric_limits/max_digits10)
-  significant digits. The rationale is to use the shortest representation while still allowing round-tripping.
+- Floating-point numbers are serialized with the fewest digits that read back as the same value (the closest such
+  digits if there are several), in the layout of the `#!c %g` printf modifier: `#!c 1.5`, `#!c 100.0`, `#!c 1e+100`.
+  Doubles are converted with the algorithm of [Żmij](https://github.com/vitaut/zmij), floats with Grisu2, which
+  can write more digits than necessary.

 !!! hint "Notes regarding precision of floating-point numbers"

@@ -540,9 +540,10 @@ therefore silently changes parse results rather than raising an error. See
  specifiers, for which the library likewise provides only `#!cpp double` and `#!cpp long double` overloads
  (`#!cpp float` is promoted to `#!cpp double`).

-If `#!cpp std::numeric_limits<NumberFloatType>` describes an IEEE 754 binary32 or binary64 number, `dump` uses the
-Grisu2 algorithm, which produces the shortest representation that round-trips. Otherwise the `snprintf` fallback with
-`max_digits10` digits is used.
+If `#!cpp std::numeric_limits<NumberFloatType>` describes an IEEE 754 binary64 number, `dump` uses the algorithm of
+Żmij, which produces the shortest representation that round-trips. For IEEE 754 binary32 numbers, it uses Grisu2,
+which produces a short representation that round-trips. Otherwise the `snprintf` fallback with `max_digits10` digits is
+used.

 ### Required for the binary formats

@@ -554,7 +555,7 @@ binary32 or binary64 field and have no encoding for `#!cpp long double`.

 | Type                     | Support                                                                                                               |
 |--------------------------|-----------------------------------------------------------------------------------------------------------------------|
-| `#!cpp double` (default) | full; short round-trip output through Grisu2                                                                          |
+| `#!cpp double` (default) | full; shortest round-trip output through Żmij                                                                         |
 | `#!cpp float`            | full; short round-trip output through Grisu2                                                                          |
 | `#!cpp long double`      | `dump` and `parse` only; the binary format writers do not compile, as they only handle IEEE 754 binary32 and binary64 |
 | any other type           | not usable                                                                                                            |
@@ -18,6 +18,8 @@ The class contains the UTF-8 Decoder from Bjoern Hoehrmann which is licensed und

 The class contains a slightly modified version of the Grisu2 algorithm from Florian Loitsch which is licensed under the [MIT License](https://opensource.org/licenses/MIT) (see above). Copyright &copy; 2009 [Florian Loitsch](https://florian.loitsch.com/)

+The class contains a port of the shortest double-to-decimal conversion of [Żmij](https://github.com/vitaut/zmij) by Victor Zverovich, which is licensed under the [MIT License](https://opensource.org/licenses/MIT) (see above). Copyright &copy; 2025 [Victor Zverovich](https://github.com/vitaut)
+
 The class contains a copy of [Hedley](https://nemequ.github.io/hedley/) from Evan Nemerson which is licensed as [CC0-1.0](https://creativecommons.org/publicdomain/zero/1.0/).

 The class contains an adapted version of the Eisel-Lemire algorithm and its table of powers of five from [fast_float](https://github.com/fastfloat/fast_float) by Daniel Lemire and contributors, which is available under the [MIT License](https://opensource.org/licenses/MIT) (used here), the Apache 2.0 License, and the Boost Software License. Copyright &copy; 2021 The fast_float authors
@@ -11,11 +11,17 @@

 #include <array> // array
 #include <cmath>   // signbit, isfinite
+#include <cstddef> // size_t
 #include <cstdint> // intN_t, uintN_t
 #include <cstring> // memcpy, memmove
 #include <limits> // numeric_limits
 #include <type_traits> // conditional

+#ifdef _MSC_VER
+    #include <cstdlib> // _byteswap_uint64
+#endif
+
+#include <nlohmann/detail/conversions/zmij.hpp>
 #include <nlohmann/detail/macro_scope.hpp>

 NLOHMANN_JSON_NAMESPACE_BEGIN
@@ -918,6 +924,87 @@ void grisu2(char* buf, int& len, int& decimal_exponent, FloatType value)
    grisu2(buf, len, decimal_exponent, w.minus, w.w, w.plus);
 }

+/*!
+@brief the shortest digits of a positive finite float (other than double): Grisu2
+*/
+template<typename FloatType>
+JSON_HEDLEY_NON_NULL(1)
+void shortest_digits(char* buf, int& len, int& decimal_exponent, FloatType value)
+{
+    grisu2(buf, len, decimal_exponent, value);
+}
+
+/*!
+@brief the shortest digits of a positive finite double: the conversion of
+Zmij (see zmij.hpp), which always finds the shortest digits that read back as
+the same value (Grisu2 does not for about one double in a thousand), and the
+closest of them if there are several
+
+v = buf * 10^decimal_exponent, as for grisu2()
+*/
+JSON_HEDLEY_NON_NULL(1)
+inline void shortest_digits(char* buf, int& len, int& decimal_exponent, double value)
+{
+    static_assert(std::numeric_limits<double>::is_iec559 && std::numeric_limits<double>::digits == 53,
+                  "internal error: the conversion of Zmij needs IEEE 754 binary64 doubles");
+    JSON_ASSERT(std::isfinite(value));
+    JSON_ASSERT(value > 0);
+
+    std::uint64_t bits = 0;
+    std::memcpy(&bits, &value, sizeof(bits));
+    zmij::decimal d = zmij::to_decimal(bits);
+    // without trailing zeros (up to 16): 8, 4, 2, 1 at a time
+    while (d.significand % 100000000 == 0)
+    {
+        d.significand /= 100000000;
+        d.exponent += 8;
+    }
+    if (d.significand % 10000 == 0)
+    {
+        d.significand /= 10000;
+        d.exponent += 4;
+    }
+    if (d.significand % 100 == 0)
+    {
+        d.significand /= 100;
+        d.exponent += 2;
+    }
+    if (d.significand % 10 == 0)
+    {
+        d.significand /= 10;
+        d.exponent += 1;
+    }
+    // at most 17 digits, written from the back two at a time
+    static constexpr const char* pairs =
+        "00010203040506070809101112131415161718192021222324252627282930313233343536373839"
+        "40414243444546474849505152535455565758596061626364656667686970717273747576777879"
+        "8081828384858687888990919293949596979899";
+    std::array<char, 20> digits{};
+    std::size_t n = digits.size();
+    while (d.significand >= 100)
+    {
+        const auto i = static_cast<std::size_t>(d.significand % 100) * 2;
+        d.significand /= 100;
+        n -= 2;
+        digits[n] = pairs[i];
+        digits[n + 1] = pairs[i + 1];
+    }
+    if (d.significand >= 10)
+    {
+        const auto i = static_cast<std::size_t>(d.significand) * 2;
+        n -= 2;
+        digits[n] = pairs[i];
+        digits[n + 1] = pairs[i + 1];
+    }
+    else
+    {
+        digits[--n] = static_cast<char>('0' + d.significand);
+    }
+    len = static_cast<int>(digits.size() - n);
+    std::memcpy(buf, digits.data() + n, static_cast<std::size_t>(len));
+    decimal_exponent = d.exponent;
+}
+
 /*!
@brief appends a decimal representation of e to buf
@return a pointer to the element following the exponent.
@@ -1047,6 +1134,177 @@ inline char* format_buffer(char* buf, int len, int decimal_exponent,
    return append_exponent(buf, n - 1);
 }

+/// eight decimal digits (a value below 10^8) as bytes 0..9, the first digit
+/// in the most significant byte: three steps that divide all lanes at once
+/// by a multiplication (the conversion of Xiang JunBo, as in Zmij)
+inline std::uint64_t eight_digit_bytes(std::uint64_t abcdefgh) noexcept
+{
+    const std::uint64_t abcd_efgh = abcdefgh + (((std::uint64_t{1} << 32u) - 10000u) * ((abcdefgh * (((std::uint64_t{1} << 40u) / 10000u) + 1u)) >> 40u));
+    const std::uint64_t ab_cd_ef_gh = abcd_efgh + (((std::uint64_t{1} << 16u) - 100u) * (((abcd_efgh * (((std::uint64_t{1} << 19u) / 100u) + 1u)) >> 19u) & 0x7F0000007Fu));
+    return ab_cd_ef_gh + (((std::uint64_t{1} << 8u) - 10u) * (((ab_cd_ef_gh * (((std::uint64_t{1} << 10u) / 10u) + 1u)) >> 10u) & 0x000F000F000F000Fu));
+}
+
+/// store the bytes of v, the most significant one first (one byte swap and
+/// one store where the byte order is known: compilers do not reliably merge
+/// the byte stores once this is inlined)
+inline void store_msb_first(char* p, std::uint64_t v) noexcept
+{
+#if defined(__BYTE_ORDER__) && defined(__ORDER_LITTLE_ENDIAN__) && __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
+    v = __builtin_bswap64(v);
+    std::memcpy(p, &v, sizeof(v));
+#elif defined(__BYTE_ORDER__) && defined(__ORDER_BIG_ENDIAN__) && __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__
+    std::memcpy(p, &v, sizeof(v));
+#elif defined(_MSC_VER) // (little-endian on all its targets)
+    v = _byteswap_uint64(v);
+    std::memcpy(p, &v, sizeof(v));
+#else
+    for (unsigned i = 0; i < 8; ++i)
+    {
+        p[i] = static_cast<char>(v >> (56u - (8u * i)));
+    }
+#endif
+}
+
+/*!
+@brief digits * 10^exp for a double, in the layout of format_buffer()
+
+The layout is that of format_buffer() with min_exp -4 and max_exp 15 (the
+digits10 of double). The digits are converted eight at a time and placed
+with fixed-size moves instead of per-digit loops and moves of the buffer.
+
+@param[in] digits  the digits (not 0, at most 17 digits; trailing zeros allowed)
+@param[in] exp     the decimal exponent of the last digit
+@return a pointer past the text; up to 41 bytes at @a first are written
+        (some beyond the returned end)
+*/
+JSON_HEDLEY_NON_NULL(1)
+JSON_HEDLEY_RETURNS_NON_NULL
+inline char* write_decimal(char* first, std::uint64_t digits, int exp) noexcept
+{
+    JSON_ASSERT(digits != 0 && digits < 100000000000000000u);
+    const std::uint64_t upper = digits / 100000000u;
+    const std::uint64_t b0 = upper / 100000000u; // (one digit: it is its own byte)
+    const std::uint64_t b1 = eight_digit_bytes(upper % 100000000u);
+    const std::uint64_t b2 = eight_digit_bytes(digits % 100000000u);
+    // leading and trailing zero digits: zero bytes, counted without division
+    int leading = 16;
+    int zeros = 16;
+    if (b0 != 0)
+    {
+        leading = count_leading_zeros(b0) / 8;
+    }
+    else if (b1 != 0)
+    {
+        leading = 8 + (count_leading_zeros(b1) / 8);
+    }
+    else
+    {
+        leading += count_leading_zeros(b2) / 8;
+    }
+    if (b2 != 0)
+    {
+        zeros = count_trailing_zeros(b2) / 8;
+    }
+    else if (b1 != 0)
+    {
+        zeros = 8 + (count_trailing_zeros(b1) / 8);
+    }
+    // (else: 16, b0 is the one digit that is not 0)
+    // the digits as text at text + leading, then '0's, so that fixed-size
+    // moves need not check how many digits there are
+    std::array<char, 64> text; // NOLINT(cppcoreguidelines-pro-type-member-init,hicpp-member-init): written before read
+    store_msb_first(text.data(), b0 + 0x3030303030303030u);
+    store_msb_first(text.data() + 8, b1 + 0x3030303030303030u);
+    store_msb_first(text.data() + 16, b2 + 0x3030303030303030u);
+    std::memset(text.data() + 24, '0', 40);
+    const int k = 24 - leading - zeros; // significant digits
+    const int n = k + exp + zeros;      // position of the decimal point after the first digit
+    const char* const s0 = text.data() + leading;
+
+    if (-4 < n && n <= 15)
+    {
+        // "0.[000]digits" (n <= 0) is the digits after 1 - n leading '0's
+        // with the point after the first; "digits[000].0" (n >= k) and
+        // "dig.its" put the point after n characters
+        const int pad = n <= 0 ? 1 - n : 0;
+        const char* const s = s0 - pad;
+        const int len = k + pad;
+        const int point = n + pad;
+        std::memcpy(first, s, 16);
+        std::memcpy(first + point + 1, s + point, 24);
+        first[point] = '.';
+        return first + (point >= len ? point + 2 : len + 1);
+    }
+
+    // d.igitse+XX, with at least two exponent digits (as append_exponent())
+    std::memcpy(first, s0, 16);
+    std::memcpy(first + 2, s0 + 1, 16);
+    first[1] = '.';
+    char* const end = first + (k == 1 ? 1 : k + 1);
+    const int e = n - 1;
+    const auto ea = static_cast<unsigned>(e < 0 ? -e : e);
+    const bool three = ea >= 100;
+    end[0] = 'e';
+    end[1] = e < 0 ? '-' : '+';
+    end[2] = static_cast<char>('0' + (three ? ea / 100 : (ea / 10) % 10));
+    end[3] = static_cast<char>('0' + (three ? (ea / 10) % 10 : ea % 10));
+    end[4] = static_cast<char>('0' + (ea % 10));
+    return end + (three ? 5 : 4);
+}
+
+/// a positive finite float (other than double): Grisu2 and format_buffer()
+template<typename FloatType>
+JSON_HEDLEY_NON_NULL(1, 2)
+JSON_HEDLEY_RETURNS_NON_NULL
+char* write_positive(char* first, const char* last, FloatType value)
+{
+    JSON_ASSERT(last - first >= std::numeric_limits<FloatType>::max_digits10);
+
+    // Compute v = buffer * 10^decimal_exponent.
+    // The decimal digits are stored in the buffer, which needs to be interpreted
+    // as an unsigned decimal integer.
+    // len is the length of the buffer, i.e., the number of decimal digits.
+    int len = 0;
+    int decimal_exponent = 0;
+    shortest_digits(first, len, decimal_exponent, value);
+
+    JSON_ASSERT(len <= std::numeric_limits<FloatType>::max_digits10);
+
+    // Format the buffer like printf("%.*g", prec, value)
+    constexpr int kMinExp = -4;
+    // Use digits10 here to increase compatibility with version 2.
+    constexpr int kMaxExp = std::numeric_limits<FloatType>::digits10;
+
+    JSON_ASSERT(last - first >= kMaxExp + 2);
+    JSON_ASSERT(last - first >= 2 + (-kMinExp - 1) + std::numeric_limits<FloatType>::max_digits10);
+    JSON_ASSERT(last - first >= std::numeric_limits<FloatType>::max_digits10 + 6);
+
+    return format_buffer(first, len, decimal_exponent, kMinExp, kMaxExp);
+}
+
+/// a positive finite double: the shortest digits (Zmij), laid out by
+/// write_decimal() (through a local buffer if [first, last) is shorter than
+/// the 41 bytes it may write)
+JSON_HEDLEY_NON_NULL(1, 2)
+JSON_HEDLEY_RETURNS_NON_NULL
+inline char* write_positive(char* first, const char* last, double value)
+{
+    static_assert(std::numeric_limits<double>::is_iec559 && std::numeric_limits<double>::digits == 53,
+                  "internal error: the conversion of Zmij needs IEEE 754 binary64 doubles");
+    std::uint64_t bits = 0;
+    std::memcpy(&bits, &value, sizeof(bits));
+    const zmij::decimal d = zmij::to_decimal(bits);
+    if (JSON_HEDLEY_LIKELY(last - first >= 41))
+    {
+        return write_decimal(first, d.significand, d.exponent);
+    }
+    std::array<char, 64> buf; // NOLINT(cppcoreguidelines-pro-type-member-init,hicpp-member-init): written before read
+    const auto len = static_cast<std::size_t>(write_decimal(buf.data(), d.significand, d.exponent) - buf.data());
+    JSON_ASSERT(static_cast<std::size_t>(last - first) >= len);
+    std::memcpy(first, buf.data(), len);
+    return first + len;
+}
+
 }  // namespace dtoa_impl

 /*!
@@ -1064,7 +1322,6 @@ JSON_HEDLEY_NON_NULL(1, 2)
 JSON_HEDLEY_RETURNS_NON_NULL
 char* to_chars(char* first, const char* last, FloatType value)
 {
-    static_cast<void>(last); // maybe unused - fix warning
    JSON_ASSERT(std::isfinite(value));

    // Use signbit(value) instead of (value < 0) since signbit works for -0.
@@ -1090,28 +1347,7 @@ char* to_chars(char* first, const char* last, FloatType value)
    JSON_HEDLEY_DIAGNOSTIC_POP
 #endif

-    JSON_ASSERT(last - first >= std::numeric_limits<FloatType>::max_digits10);
-
-    // Compute v = buffer * 10^decimal_exponent.
-    // The decimal digits are stored in the buffer, which needs to be interpreted
-    // as an unsigned decimal integer.
-    // len is the length of the buffer, i.e., the number of decimal digits.
-    int len = 0;
-    int decimal_exponent = 0;
-    dtoa_impl::grisu2(first, len, decimal_exponent, value);
-
-    JSON_ASSERT(len <= std::numeric_limits<FloatType>::max_digits10);
-
-    // Format the buffer like printf("%.*g", prec, value)
-    constexpr int kMinExp = -4;
-    // Use digits10 here to increase compatibility with version 2.
-    constexpr int kMaxExp = std::numeric_limits<FloatType>::digits10;
-
-    JSON_ASSERT(last - first >= kMaxExp + 2);
-    JSON_ASSERT(last - first >= 2 + (-kMinExp - 1) + std::numeric_limits<FloatType>::max_digits10);
-    JSON_ASSERT(last - first >= std::numeric_limits<FloatType>::max_digits10 + 6);
-
-    return dtoa_impl::format_buffer(first, len, decimal_exponent, kMinExp, kMaxExp);
+    return dtoa_impl::write_positive(first, last, value);
 }

 }  // namespace detail
@@ -0,0 +1,218 @@
+//     __ _____ _____ _____
+//  __|  |   __|     |   | |  JSON for Modern C++
+// |  |  |__   |  |  | | | |  version 3.12.0
+// |_____|_____|_____|_|___|  https://github.com/nlohmann/json
+//
+// SPDX-FileCopyrightText: 2025 Victor Zverovich <https://github.com/vitaut/zmij>
+// SPDX-FileCopyrightText: 2013-2026 Niels Lohmann <https://nlohmann.me>
+// SPDX-License-Identifier: MIT
+
+#pragma once
+
+#include <array> // array
+#include <cstddef> // size_t
+#include <cstdint> // uint32_t, uint64_t
+
+#include <nlohmann/detail/abi_macros.hpp>
+#include <nlohmann/detail/bit_ops.hpp>
+#include <nlohmann/detail/input/pow5_table.hpp>
+#include <nlohmann/detail/macro_scope.hpp>
+
+NLOHMANN_JSON_NAMESPACE_BEGIN
+namespace detail
+{
+
+/*!
+@brief the shortest decimal representation of a double
+
+A C++11 port of the conversion of Zmij by Victor Zverovich
+(https://github.com/vitaut/zmij, MIT license): the shortest decimal in the
+rounding interval of a double, the closest one if there are several. Zmij
+credits Xiang JunBo (producing the shorter candidate without a division) and
+Dougall Johnson (the compressed powers of ten). The powers of ten are taken
+from the table for number parsing (pow5_table.hpp) where it holds them, and
+computed from the compressed tables of Zmij beyond it.
+*/
+namespace zmij
+{
+
+/// significand * 10^exponent
+struct decimal
+{
+    std::uint64_t significand;
+    int exponent;
+};
+
+/// the compressed powers of ten of Zmij
+inline const std::array<std::uint64_t, 28>& pow10_minor() noexcept
+{
+    static const std::array<std::uint64_t, 28> table =
+    {
+        {
+            0x8000000000000000u, 0xa000000000000000u, 0xc800000000000000u, 0xfa00000000000000u, 0x9c40000000000000u,
+            0xc350000000000000u, 0xf424000000000000u, 0x9896800000000000u, 0xbebc200000000000u, 0xee6b280000000000u,
+            0x9502f90000000000u, 0xba43b74000000000u, 0xe8d4a51000000000u, 0x9184e72a00000000u, 0xb5e620f480000000u,
+            0xe35fa931a0000000u, 0x8e1bc9bf04000000u, 0xb1a2bc2ec5000000u, 0xde0b6b3a76400000u, 0x8ac7230489e80000u,
+            0xad78ebc5ac620000u, 0xd8d726b7177a8000u, 0x878678326eac9000u, 0xa968163f0a57b400u, 0xd3c21bcecceda100u,
+            0x84595161401484a0u, 0xa56fa5b99019a5c8u, 0xcecb8f27f4200f3au
+        }
+    };
+    return table;
+}
+
+/// (high, low) pairs
+inline const std::array<std::uint64_t, 50>& pow10_major() noexcept
+{
+    static const std::array<std::uint64_t, 50> table =
+    {
+        {
+            0xaddcb9e83c6b1793u, 0xdf4abe242a1bbf3eu, 0xaf8e5410288e1b6fu, 0x07ecf0ae5ee44ddau, 0xb1442798f49ffb4au, 0x99cd11cfdf41779du,
+            0xb2fe3f0b8599ef07u, 0x861fa7e6dcb4aa15u, 0xb4bca50b065abe63u, 0x0fed077a756b53aau, 0xb67f6455292cbf08u, 0x1a3bc84c17b1d543u,
+            0xb84687c269ef3bfbu, 0x3d5d514f40eea742u, 0xba121a4650e4ddebu, 0x92f34d62616ce413u, 0xbbe226efb628afeau, 0x890489f70a55368cu,
+            0xbdb6b8e905cb600fu, 0x5400e987bbc1c921u, 0xbf8fdb78849a5f96u, 0xde98520472bdd034u, 0xc16d9a0095928a27u, 0x75b7053c0f178294u,
+            0xc350000000000000u, 0x0000000000000000u, 0xc5371912364ce305u, 0x6c28000000000000u, 0xc722f0ef9d80aad6u, 0x424d3ad2b7b97ef6u,
+            0xc913936dd571c84cu, 0x03bc3a19cd1e38eau, 0xcb090c8001ab551cu, 0x5cadf5bfd3072cc6u, 0xcd036837130890a1u, 0x36dba887c37a8c10u,
+            0xcf02b2c21207ef2eu, 0x94f967e45e03f4bcu, 0xd106f86e69d785c7u, 0xe13336d701beba52u, 0xd31045a8341ca07cu, 0x1ede48111209a051u,
+            0xd51ea6fa85785631u, 0x552a74227f3ea566u, 0xd732290fbacaf133u, 0xa97c177947ad4096u, 0xd94ad8b1c7380874u, 0x18375281ae7822bdu,
+            0xdb68c2ca82ed2a05u, 0xa67398db9f6820e1u
+        }
+    };
+    return table;
+}
+
+/// one bit per power: whether the computed value is one unit too large
+inline const std::array<std::uint32_t, 21>& pow10_fixups() noexcept
+{
+    static const std::array<std::uint32_t, 21> table =
+    {
+        {
+            0x8d8fc810u, 0x06100293u, 0x19000000u, 0x00100000u, 0x00000908u, 0x00000000u, 0x04e00300u, 0x3807e0b2u, 0x3d83d793u, 0x0006f5ccu,
+            0x00000000u, 0xffff0000u, 0x8076337du, 0x4ff45ba0u, 0x09405033u, 0x034376d9u, 0x09000000u, 0x4e100501u, 0x076d14dcu, 0xf964f45eu,
+            0x0000003du
+        }
+    };
+    return table;
+}
+
+/// the 128-bit significand of 10^k, rounded down, for k in [-307, 341]
+/// (compute_pow10 of Zmij)
+inline uint128_parts compute_pow10(int k) noexcept
+{
+    const auto i = static_cast<unsigned>(k + 307);
+    const std::uint64_t m = pow10_minor()[(i + 24) % 28];
+    const std::size_t j = 2 * static_cast<std::size_t>((i + 24) / 28);
+    const std::uint64_t h_hi = pow10_major()[j];
+    const std::uint64_t h_lo = pow10_major()[j + 1];
+    const std::uint64_t h1 = full_multiplication(h_lo, m).high;
+    const std::uint64_t c0 = h_lo * m;
+    const std::uint64_t c1 = h1 + (h_hi * m);
+    const std::uint64_t c2 = (c1 < h1 ? 1u : 0u) + full_multiplication(h_hi, m).high;
+    uint128_parts r{};
+    if ((c2 >> 63u) != 0)
+    {
+        r.high = c2;
+        r.low = c1;
+    }
+    else
+    {
+        r.high = (c2 << 1u) | (c1 >> 63u);
+        r.low = (c1 << 1u) | (c0 >> 63u);
+    }
+    r.low -= (pow10_fixups()[i >> 5u] >> (i & 31u)) & 1u;
+    return r;
+}
+
+/// The 128-bit significand of 10^k, rounded down, for k in [-342, 341].
+/// Up to 10^308, the table for number parsing holds the same significands
+/// (those of 5^k), except for k in [-27, -1], where it holds them one unit
+/// larger (as the Eisel-Lemire algorithm needs them).
+inline uint128_parts pow10(int k) noexcept
+{
+    if (k > pow5_128_largest_power)
+    {
+        return compute_pow10(k); // (only for the smallest doubles)
+    }
+    const auto i = 2 * static_cast<std::size_t>(k - pow5_128_smallest_power);
+    uint128_parts r{pow5_128()[i + 1], pow5_128()[i]};
+    const std::uint64_t adjust = static_cast<unsigned>(k + 27) < 27u ? 1u : 0u;
+    r.high -= r.low < adjust ? 1u : 0u;
+    r.low -= adjust;
+    return r;
+}
+
+/// (x_hi * 2^64 + x_lo) * y >> 64, as 128 bits
+inline uint128_parts umul192_hi128(std::uint64_t x_hi, std::uint64_t x_lo, std::uint64_t y) noexcept
+{
+    const uint128_parts p = full_multiplication(x_hi, y);
+    uint128_parts r{};
+    r.low = p.low + full_multiplication(x_lo, y).high;
+    r.high = p.high + (r.low < p.low ? 1u : 0u);
+    return r;
+}
+
+/// (x * y + c) >> 64
+inline std::uint64_t umul128_add_hi64(std::uint64_t x, std::uint64_t y, std::uint64_t c) noexcept
+{
+    const uint128_parts p = full_multiplication(x, y);
+    return p.high + (p.low + c < p.low ? 1u : 0u);
+}
+
+/// The shortest decimal in the rounding interval of a positive finite double
+/// given by its bits, the closest one if there are several (to_decimal of
+/// Zmij). The significand can end in zeros.
+inline decimal to_decimal(std::uint64_t bits) noexcept
+{
+    constexpr int extra_shift = 9;
+    const auto raw_exp = static_cast<int>((bits >> 52u) & 0x7FFu);
+    std::uint64_t bin_sig = bits & ((std::uint64_t{1} << 52u) - 1);
+    // a power of two has a narrower interval below (except the smallest normal)
+    const bool regular = bin_sig != 0 || raw_exp <= 1;
+    const int bin_exp = (raw_exp == 0 ? 1 : raw_exp) - 1075;
+    if (raw_exp != 0)
+    {
+        bin_sig |= std::uint64_t{1} << 52u;
+    }
+    // floor(log10(2^bin_exp)), or floor(log10(3/4 * 2^bin_exp)) for the irregular case
+    const int dec_exp = ((bin_exp * 315653) - (regular ? 0 : 131072)) >> 20;
+    // scaled by 10^(-dec_exp - 1): the integral part is the shorter candidate
+    const int shift = bin_exp + ((-(dec_exp + 1) * 217707) >> 16) + 1 + extra_shift;
+    const uint128_parts p10 = pow10(-dec_exp - 1);
+    const uint128_parts p = umul192_hi128(p10.high, p10.low, bin_sig << static_cast<unsigned>(shift));
+    std::uint64_t integral = p.high >> static_cast<unsigned>(extra_shift);
+    const std::uint64_t fractional = (p.high << static_cast<unsigned>(64 - extra_shift)) | (p.low >> static_cast<unsigned>(extra_shift));
+    std::uint64_t digit = 0;
+    bool round_up = false;
+    bool round_down = false;
+    if (JSON_HEDLEY_LIKELY(regular))
+    {
+        const std::uint64_t half_ulp = (p10.high >> static_cast<unsigned>(extra_shift + 1 - shift)) + (1 - (bin_sig & 1u));
+        round_up = fractional + half_ulp < fractional;
+        round_down = half_ulp > fractional;
+        // the last digit of the longer candidate, rounded to nearest
+        digit = umul128_add_hi64(fractional, 10, (std::uint64_t{1} << 63u) + 6);
+        if (fractional == (std::uint64_t{1} << 62u))
+        {
+            digit = 2; // 2.5 rounds to 2
+        }
+    }
+    else
+    {
+        const std::uint64_t half_ulp = p10.high >> static_cast<unsigned>(extra_shift + 1 - shift);
+        round_up = half_ulp > ~std::uint64_t{0} - fractional;
+        round_down = (half_ulp >> 1u) > fractional;
+        digit = umul128_add_hi64(fractional, 10, (std::uint64_t{1} << 63u) - 1);
+        const std::uint64_t lowest = umul128_add_hi64(fractional - (half_ulp >> 1u), 10, ~std::uint64_t{0});
+        digit = digit < lowest ? lowest : digit;
+    }
+    integral += round_up ? 1u : 0u;
+    if (!round_up && !round_down)
+    {
+        // the shorter candidate is outside the rounding interval: one digit more
+        return decimal{(integral * 10) + digit, dec_exp};
+    }
+    return decimal{integral, dec_exp + 1};
+}
+
+}  // namespace zmij
+}  // namespace detail
+NLOHMANN_JSON_NAMESPACE_END
@@ -23918,11 +23918,240 @@ NLOHMANN_JSON_NAMESPACE_END

 #include <array> // array
 #include <cmath>   // signbit, isfinite
+#include <cstddef> // size_t
 #include <cstdint> // intN_t, uintN_t
 #include <cstring> // memcpy, memmove
 #include <limits> // numeric_limits
 #include <type_traits> // conditional

+#ifdef _MSC_VER
+    #include <cstdlib> // _byteswap_uint64
+#endif
+
+// #include <nlohmann/detail/conversions/zmij.hpp>
+//     __ _____ _____ _____
+//  __|  |   __|     |   | |  JSON for Modern C++
+// |  |  |__   |  |  | | | |  version 3.12.0
+// |_____|_____|_____|_|___|  https://github.com/nlohmann/json
+//
+// SPDX-FileCopyrightText: 2025 Victor Zverovich <https://github.com/vitaut/zmij>
+// SPDX-FileCopyrightText: 2013-2026 Niels Lohmann <https://nlohmann.me>
+// SPDX-License-Identifier: MIT
+
+
+
+#include <array> // array
+#include <cstddef> // size_t
+#include <cstdint> // uint32_t, uint64_t
+
+// #include <nlohmann/detail/abi_macros.hpp>
+
+// #include <nlohmann/detail/bit_ops.hpp>
+
+// #include <nlohmann/detail/input/pow5_table.hpp>
+
+// #include <nlohmann/detail/macro_scope.hpp>
+
+
+NLOHMANN_JSON_NAMESPACE_BEGIN
+namespace detail
+{
+
+/*!
+@brief the shortest decimal representation of a double
+
+A C++11 port of the conversion of Zmij by Victor Zverovich
+(https://github.com/vitaut/zmij, MIT license): the shortest decimal in the
+rounding interval of a double, the closest one if there are several. Zmij
+credits Xiang JunBo (producing the shorter candidate without a division) and
+Dougall Johnson (the compressed powers of ten). The powers of ten are taken
+from the table for number parsing (pow5_table.hpp) where it holds them, and
+computed from the compressed tables of Zmij beyond it.
+*/
+namespace zmij
+{
+
+/// significand * 10^exponent
+struct decimal
+{
+    std::uint64_t significand;
+    int exponent;
+};
+
+/// the compressed powers of ten of Zmij
+inline const std::array<std::uint64_t, 28>& pow10_minor() noexcept
+{
+    static const std::array<std::uint64_t, 28> table =
+    {
+        {
+            0x8000000000000000u, 0xa000000000000000u, 0xc800000000000000u, 0xfa00000000000000u, 0x9c40000000000000u,
+            0xc350000000000000u, 0xf424000000000000u, 0x9896800000000000u, 0xbebc200000000000u, 0xee6b280000000000u,
+            0x9502f90000000000u, 0xba43b74000000000u, 0xe8d4a51000000000u, 0x9184e72a00000000u, 0xb5e620f480000000u,
+            0xe35fa931a0000000u, 0x8e1bc9bf04000000u, 0xb1a2bc2ec5000000u, 0xde0b6b3a76400000u, 0x8ac7230489e80000u,
+            0xad78ebc5ac620000u, 0xd8d726b7177a8000u, 0x878678326eac9000u, 0xa968163f0a57b400u, 0xd3c21bcecceda100u,
+            0x84595161401484a0u, 0xa56fa5b99019a5c8u, 0xcecb8f27f4200f3au
+        }
+    };
+    return table;
+}
+
+/// (high, low) pairs
+inline const std::array<std::uint64_t, 50>& pow10_major() noexcept
+{
+    static const std::array<std::uint64_t, 50> table =
+    {
+        {
+            0xaddcb9e83c6b1793u, 0xdf4abe242a1bbf3eu, 0xaf8e5410288e1b6fu, 0x07ecf0ae5ee44ddau, 0xb1442798f49ffb4au, 0x99cd11cfdf41779du,
+            0xb2fe3f0b8599ef07u, 0x861fa7e6dcb4aa15u, 0xb4bca50b065abe63u, 0x0fed077a756b53aau, 0xb67f6455292cbf08u, 0x1a3bc84c17b1d543u,
+            0xb84687c269ef3bfbu, 0x3d5d514f40eea742u, 0xba121a4650e4ddebu, 0x92f34d62616ce413u, 0xbbe226efb628afeau, 0x890489f70a55368cu,
+            0xbdb6b8e905cb600fu, 0x5400e987bbc1c921u, 0xbf8fdb78849a5f96u, 0xde98520472bdd034u, 0xc16d9a0095928a27u, 0x75b7053c0f178294u,
+            0xc350000000000000u, 0x0000000000000000u, 0xc5371912364ce305u, 0x6c28000000000000u, 0xc722f0ef9d80aad6u, 0x424d3ad2b7b97ef6u,
+            0xc913936dd571c84cu, 0x03bc3a19cd1e38eau, 0xcb090c8001ab551cu, 0x5cadf5bfd3072cc6u, 0xcd036837130890a1u, 0x36dba887c37a8c10u,
+            0xcf02b2c21207ef2eu, 0x94f967e45e03f4bcu, 0xd106f86e69d785c7u, 0xe13336d701beba52u, 0xd31045a8341ca07cu, 0x1ede48111209a051u,
+            0xd51ea6fa85785631u, 0x552a74227f3ea566u, 0xd732290fbacaf133u, 0xa97c177947ad4096u, 0xd94ad8b1c7380874u, 0x18375281ae7822bdu,
+            0xdb68c2ca82ed2a05u, 0xa67398db9f6820e1u
+        }
+    };
+    return table;
+}
+
+/// one bit per power: whether the computed value is one unit too large
+inline const std::array<std::uint32_t, 21>& pow10_fixups() noexcept
+{
+    static const std::array<std::uint32_t, 21> table =
+    {
+        {
+            0x8d8fc810u, 0x06100293u, 0x19000000u, 0x00100000u, 0x00000908u, 0x00000000u, 0x04e00300u, 0x3807e0b2u, 0x3d83d793u, 0x0006f5ccu,
+            0x00000000u, 0xffff0000u, 0x8076337du, 0x4ff45ba0u, 0x09405033u, 0x034376d9u, 0x09000000u, 0x4e100501u, 0x076d14dcu, 0xf964f45eu,
+            0x0000003du
+        }
+    };
+    return table;
+}
+
+/// the 128-bit significand of 10^k, rounded down, for k in [-307, 341]
+/// (compute_pow10 of Zmij)
+inline uint128_parts compute_pow10(int k) noexcept
+{
+    const auto i = static_cast<unsigned>(k + 307);
+    const std::uint64_t m = pow10_minor()[(i + 24) % 28];
+    const std::size_t j = 2 * static_cast<std::size_t>((i + 24) / 28);
+    const std::uint64_t h_hi = pow10_major()[j];
+    const std::uint64_t h_lo = pow10_major()[j + 1];
+    const std::uint64_t h1 = full_multiplication(h_lo, m).high;
+    const std::uint64_t c0 = h_lo * m;
+    const std::uint64_t c1 = h1 + (h_hi * m);
+    const std::uint64_t c2 = (c1 < h1 ? 1u : 0u) + full_multiplication(h_hi, m).high;
+    uint128_parts r{};
+    if ((c2 >> 63u) != 0)
+    {
+        r.high = c2;
+        r.low = c1;
+    }
+    else
+    {
+        r.high = (c2 << 1u) | (c1 >> 63u);
+        r.low = (c1 << 1u) | (c0 >> 63u);
+    }
+    r.low -= (pow10_fixups()[i >> 5u] >> (i & 31u)) & 1u;
+    return r;
+}
+
+/// The 128-bit significand of 10^k, rounded down, for k in [-342, 341].
+/// Up to 10^308, the table for number parsing holds the same significands
+/// (those of 5^k), except for k in [-27, -1], where it holds them one unit
+/// larger (as the Eisel-Lemire algorithm needs them).
+inline uint128_parts pow10(int k) noexcept
+{
+    if (k > pow5_128_largest_power)
+    {
+        return compute_pow10(k); // (only for the smallest doubles)
+    }
+    const auto i = 2 * static_cast<std::size_t>(k - pow5_128_smallest_power);
+    uint128_parts r{pow5_128()[i + 1], pow5_128()[i]};
+    const std::uint64_t adjust = static_cast<unsigned>(k + 27) < 27u ? 1u : 0u;
+    r.high -= r.low < adjust ? 1u : 0u;
+    r.low -= adjust;
+    return r;
+}
+
+/// (x_hi * 2^64 + x_lo) * y >> 64, as 128 bits
+inline uint128_parts umul192_hi128(std::uint64_t x_hi, std::uint64_t x_lo, std::uint64_t y) noexcept
+{
+    const uint128_parts p = full_multiplication(x_hi, y);
+    uint128_parts r{};
+    r.low = p.low + full_multiplication(x_lo, y).high;
+    r.high = p.high + (r.low < p.low ? 1u : 0u);
+    return r;
+}
+
+/// (x * y + c) >> 64
+inline std::uint64_t umul128_add_hi64(std::uint64_t x, std::uint64_t y, std::uint64_t c) noexcept
+{
+    const uint128_parts p = full_multiplication(x, y);
+    return p.high + (p.low + c < p.low ? 1u : 0u);
+}
+
+/// The shortest decimal in the rounding interval of a positive finite double
+/// given by its bits, the closest one if there are several (to_decimal of
+/// Zmij). The significand can end in zeros.
+inline decimal to_decimal(std::uint64_t bits) noexcept
+{
+    constexpr int extra_shift = 9;
+    const auto raw_exp = static_cast<int>((bits >> 52u) & 0x7FFu);
+    std::uint64_t bin_sig = bits & ((std::uint64_t{1} << 52u) - 1);
+    // a power of two has a narrower interval below (except the smallest normal)
+    const bool regular = bin_sig != 0 || raw_exp <= 1;
+    const int bin_exp = (raw_exp == 0 ? 1 : raw_exp) - 1075;
+    if (raw_exp != 0)
+    {
+        bin_sig |= std::uint64_t{1} << 52u;
+    }
+    // floor(log10(2^bin_exp)), or floor(log10(3/4 * 2^bin_exp)) for the irregular case
+    const int dec_exp = ((bin_exp * 315653) - (regular ? 0 : 131072)) >> 20;
+    // scaled by 10^(-dec_exp - 1): the integral part is the shorter candidate
+    const int shift = bin_exp + ((-(dec_exp + 1) * 217707) >> 16) + 1 + extra_shift;
+    const uint128_parts p10 = pow10(-dec_exp - 1);
+    const uint128_parts p = umul192_hi128(p10.high, p10.low, bin_sig << static_cast<unsigned>(shift));
+    std::uint64_t integral = p.high >> static_cast<unsigned>(extra_shift);
+    const std::uint64_t fractional = (p.high << static_cast<unsigned>(64 - extra_shift)) | (p.low >> static_cast<unsigned>(extra_shift));
+    std::uint64_t digit = 0;
+    bool round_up = false;
+    bool round_down = false;
+    if (JSON_HEDLEY_LIKELY(regular))
+    {
+        const std::uint64_t half_ulp = (p10.high >> static_cast<unsigned>(extra_shift + 1 - shift)) + (1 - (bin_sig & 1u));
+        round_up = fractional + half_ulp < fractional;
+        round_down = half_ulp > fractional;
+        // the last digit of the longer candidate, rounded to nearest
+        digit = umul128_add_hi64(fractional, 10, (std::uint64_t{1} << 63u) + 6);
+        if (fractional == (std::uint64_t{1} << 62u))
+        {
+            digit = 2; // 2.5 rounds to 2
+        }
+    }
+    else
+    {
+        const std::uint64_t half_ulp = p10.high >> static_cast<unsigned>(extra_shift + 1 - shift);
+        round_up = half_ulp > ~std::uint64_t{0} - fractional;
+        round_down = (half_ulp >> 1u) > fractional;
+        digit = umul128_add_hi64(fractional, 10, (std::uint64_t{1} << 63u) - 1);
+        const std::uint64_t lowest = umul128_add_hi64(fractional - (half_ulp >> 1u), 10, ~std::uint64_t{0});
+        digit = digit < lowest ? lowest : digit;
+    }
+    integral += round_up ? 1u : 0u;
+    if (!round_up && !round_down)
+    {
+        // the shorter candidate is outside the rounding interval: one digit more
+        return decimal{(integral * 10) + digit, dec_exp};
+    }
+    return decimal{integral, dec_exp + 1};
+}
+
+}  // namespace zmij
+}  // namespace detail
+NLOHMANN_JSON_NAMESPACE_END
+
 // #include <nlohmann/detail/macro_scope.hpp>


@@ -24826,6 +25055,87 @@ void grisu2(char* buf, int& len, int& decimal_exponent, FloatType value)
    grisu2(buf, len, decimal_exponent, w.minus, w.w, w.plus);
 }

+/*!
+@brief the shortest digits of a positive finite float (other than double): Grisu2
+*/
+template<typename FloatType>
+JSON_HEDLEY_NON_NULL(1)
+void shortest_digits(char* buf, int& len, int& decimal_exponent, FloatType value)
+{
+    grisu2(buf, len, decimal_exponent, value);
+}
+
+/*!
+@brief the shortest digits of a positive finite double: the conversion of
+Zmij (see zmij.hpp), which always finds the shortest digits that read back as
+the same value (Grisu2 does not for about one double in a thousand), and the
+closest of them if there are several
+
+v = buf * 10^decimal_exponent, as for grisu2()
+*/
+JSON_HEDLEY_NON_NULL(1)
+inline void shortest_digits(char* buf, int& len, int& decimal_exponent, double value)
+{
+    static_assert(std::numeric_limits<double>::is_iec559 && std::numeric_limits<double>::digits == 53,
+                  "internal error: the conversion of Zmij needs IEEE 754 binary64 doubles");
+    JSON_ASSERT(std::isfinite(value));
+    JSON_ASSERT(value > 0);
+
+    std::uint64_t bits = 0;
+    std::memcpy(&bits, &value, sizeof(bits));
+    zmij::decimal d = zmij::to_decimal(bits);
+    // without trailing zeros (up to 16): 8, 4, 2, 1 at a time
+    while (d.significand % 100000000 == 0)
+    {
+        d.significand /= 100000000;
+        d.exponent += 8;
+    }
+    if (d.significand % 10000 == 0)
+    {
+        d.significand /= 10000;
+        d.exponent += 4;
+    }
+    if (d.significand % 100 == 0)
+    {
+        d.significand /= 100;
+        d.exponent += 2;
+    }
+    if (d.significand % 10 == 0)
+    {
+        d.significand /= 10;
+        d.exponent += 1;
+    }
+    // at most 17 digits, written from the back two at a time
+    static constexpr const char* pairs =
+        "00010203040506070809101112131415161718192021222324252627282930313233343536373839"
+        "40414243444546474849505152535455565758596061626364656667686970717273747576777879"
+        "8081828384858687888990919293949596979899";
+    std::array<char, 20> digits{};
+    std::size_t n = digits.size();
+    while (d.significand >= 100)
+    {
+        const auto i = static_cast<std::size_t>(d.significand % 100) * 2;
+        d.significand /= 100;
+        n -= 2;
+        digits[n] = pairs[i];
+        digits[n + 1] = pairs[i + 1];
+    }
+    if (d.significand >= 10)
+    {
+        const auto i = static_cast<std::size_t>(d.significand) * 2;
+        n -= 2;
+        digits[n] = pairs[i];
+        digits[n + 1] = pairs[i + 1];
+    }
+    else
+    {
+        digits[--n] = static_cast<char>('0' + d.significand);
+    }
+    len = static_cast<int>(digits.size() - n);
+    std::memcpy(buf, digits.data() + n, static_cast<std::size_t>(len));
+    decimal_exponent = d.exponent;
+}
+
 /*!
@brief appends a decimal representation of e to buf
@return a pointer to the element following the exponent.
@@ -24955,6 +25265,177 @@ inline char* format_buffer(char* buf, int len, int decimal_exponent,
    return append_exponent(buf, n - 1);
 }

+/// eight decimal digits (a value below 10^8) as bytes 0..9, the first digit
+/// in the most significant byte: three steps that divide all lanes at once
+/// by a multiplication (the conversion of Xiang JunBo, as in Zmij)
+inline std::uint64_t eight_digit_bytes(std::uint64_t abcdefgh) noexcept
+{
+    const std::uint64_t abcd_efgh = abcdefgh + (((std::uint64_t{1} << 32u) - 10000u) * ((abcdefgh * (((std::uint64_t{1} << 40u) / 10000u) + 1u)) >> 40u));
+    const std::uint64_t ab_cd_ef_gh = abcd_efgh + (((std::uint64_t{1} << 16u) - 100u) * (((abcd_efgh * (((std::uint64_t{1} << 19u) / 100u) + 1u)) >> 19u) & 0x7F0000007Fu));
+    return ab_cd_ef_gh + (((std::uint64_t{1} << 8u) - 10u) * (((ab_cd_ef_gh * (((std::uint64_t{1} << 10u) / 10u) + 1u)) >> 10u) & 0x000F000F000F000Fu));
+}
+
+/// store the bytes of v, the most significant one first (one byte swap and
+/// one store where the byte order is known: compilers do not reliably merge
+/// the byte stores once this is inlined)
+inline void store_msb_first(char* p, std::uint64_t v) noexcept
+{
+#if defined(__BYTE_ORDER__) && defined(__ORDER_LITTLE_ENDIAN__) && __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
+    v = __builtin_bswap64(v);
+    std::memcpy(p, &v, sizeof(v));
+#elif defined(__BYTE_ORDER__) && defined(__ORDER_BIG_ENDIAN__) && __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__
+    std::memcpy(p, &v, sizeof(v));
+#elif defined(_MSC_VER) // (little-endian on all its targets)
+    v = _byteswap_uint64(v);
+    std::memcpy(p, &v, sizeof(v));
+#else
+    for (unsigned i = 0; i < 8; ++i)
+    {
+        p[i] = static_cast<char>(v >> (56u - (8u * i)));
+    }
+#endif
+}
+
+/*!
+@brief digits * 10^exp for a double, in the layout of format_buffer()
+
+The layout is that of format_buffer() with min_exp -4 and max_exp 15 (the
+digits10 of double). The digits are converted eight at a time and placed
+with fixed-size moves instead of per-digit loops and moves of the buffer.
+
+@param[in] digits  the digits (not 0, at most 17 digits; trailing zeros allowed)
+@param[in] exp     the decimal exponent of the last digit
+@return a pointer past the text; up to 41 bytes at @a first are written
+        (some beyond the returned end)
+*/
+JSON_HEDLEY_NON_NULL(1)
+JSON_HEDLEY_RETURNS_NON_NULL
+inline char* write_decimal(char* first, std::uint64_t digits, int exp) noexcept
+{
+    JSON_ASSERT(digits != 0 && digits < 100000000000000000u);
+    const std::uint64_t upper = digits / 100000000u;
+    const std::uint64_t b0 = upper / 100000000u; // (one digit: it is its own byte)
+    const std::uint64_t b1 = eight_digit_bytes(upper % 100000000u);
+    const std::uint64_t b2 = eight_digit_bytes(digits % 100000000u);
+    // leading and trailing zero digits: zero bytes, counted without division
+    int leading = 16;
+    int zeros = 16;
+    if (b0 != 0)
+    {
+        leading = count_leading_zeros(b0) / 8;
+    }
+    else if (b1 != 0)
+    {
+        leading = 8 + (count_leading_zeros(b1) / 8);
+    }
+    else
+    {
+        leading += count_leading_zeros(b2) / 8;
+    }
+    if (b2 != 0)
+    {
+        zeros = count_trailing_zeros(b2) / 8;
+    }
+    else if (b1 != 0)
+    {
+        zeros = 8 + (count_trailing_zeros(b1) / 8);
+    }
+    // (else: 16, b0 is the one digit that is not 0)
+    // the digits as text at text + leading, then '0's, so that fixed-size
+    // moves need not check how many digits there are
+    std::array<char, 64> text; // NOLINT(cppcoreguidelines-pro-type-member-init,hicpp-member-init): written before read
+    store_msb_first(text.data(), b0 + 0x3030303030303030u);
+    store_msb_first(text.data() + 8, b1 + 0x3030303030303030u);
+    store_msb_first(text.data() + 16, b2 + 0x3030303030303030u);
+    std::memset(text.data() + 24, '0', 40);
+    const int k = 24 - leading - zeros; // significant digits
+    const int n = k + exp + zeros;      // position of the decimal point after the first digit
+    const char* const s0 = text.data() + leading;
+
+    if (-4 < n && n <= 15)
+    {
+        // "0.[000]digits" (n <= 0) is the digits after 1 - n leading '0's
+        // with the point after the first; "digits[000].0" (n >= k) and
+        // "dig.its" put the point after n characters
+        const int pad = n <= 0 ? 1 - n : 0;
+        const char* const s = s0 - pad;
+        const int len = k + pad;
+        const int point = n + pad;
+        std::memcpy(first, s, 16);
+        std::memcpy(first + point + 1, s + point, 24);
+        first[point] = '.';
+        return first + (point >= len ? point + 2 : len + 1);
+    }
+
+    // d.igitse+XX, with at least two exponent digits (as append_exponent())
+    std::memcpy(first, s0, 16);
+    std::memcpy(first + 2, s0 + 1, 16);
+    first[1] = '.';
+    char* const end = first + (k == 1 ? 1 : k + 1);
+    const int e = n - 1;
+    const auto ea = static_cast<unsigned>(e < 0 ? -e : e);
+    const bool three = ea >= 100;
+    end[0] = 'e';
+    end[1] = e < 0 ? '-' : '+';
+    end[2] = static_cast<char>('0' + (three ? ea / 100 : (ea / 10) % 10));
+    end[3] = static_cast<char>('0' + (three ? (ea / 10) % 10 : ea % 10));
+    end[4] = static_cast<char>('0' + (ea % 10));
+    return end + (three ? 5 : 4);
+}
+
+/// a positive finite float (other than double): Grisu2 and format_buffer()
+template<typename FloatType>
+JSON_HEDLEY_NON_NULL(1, 2)
+JSON_HEDLEY_RETURNS_NON_NULL
+char* write_positive(char* first, const char* last, FloatType value)
+{
+    JSON_ASSERT(last - first >= std::numeric_limits<FloatType>::max_digits10);
+
+    // Compute v = buffer * 10^decimal_exponent.
+    // The decimal digits are stored in the buffer, which needs to be interpreted
+    // as an unsigned decimal integer.
+    // len is the length of the buffer, i.e., the number of decimal digits.
+    int len = 0;
+    int decimal_exponent = 0;
+    shortest_digits(first, len, decimal_exponent, value);
+
+    JSON_ASSERT(len <= std::numeric_limits<FloatType>::max_digits10);
+
+    // Format the buffer like printf("%.*g", prec, value)
+    constexpr int kMinExp = -4;
+    // Use digits10 here to increase compatibility with version 2.
+    constexpr int kMaxExp = std::numeric_limits<FloatType>::digits10;
+
+    JSON_ASSERT(last - first >= kMaxExp + 2);
+    JSON_ASSERT(last - first >= 2 + (-kMinExp - 1) + std::numeric_limits<FloatType>::max_digits10);
+    JSON_ASSERT(last - first >= std::numeric_limits<FloatType>::max_digits10 + 6);
+
+    return format_buffer(first, len, decimal_exponent, kMinExp, kMaxExp);
+}
+
+/// a positive finite double: the shortest digits (Zmij), laid out by
+/// write_decimal() (through a local buffer if [first, last) is shorter than
+/// the 41 bytes it may write)
+JSON_HEDLEY_NON_NULL(1, 2)
+JSON_HEDLEY_RETURNS_NON_NULL
+inline char* write_positive(char* first, const char* last, double value)
+{
+    static_assert(std::numeric_limits<double>::is_iec559 && std::numeric_limits<double>::digits == 53,
+                  "internal error: the conversion of Zmij needs IEEE 754 binary64 doubles");
+    std::uint64_t bits = 0;
+    std::memcpy(&bits, &value, sizeof(bits));
+    const zmij::decimal d = zmij::to_decimal(bits);
+    if (JSON_HEDLEY_LIKELY(last - first >= 41))
+    {
+        return write_decimal(first, d.significand, d.exponent);
+    }
+    std::array<char, 64> buf; // NOLINT(cppcoreguidelines-pro-type-member-init,hicpp-member-init): written before read
+    const auto len = static_cast<std::size_t>(write_decimal(buf.data(), d.significand, d.exponent) - buf.data());
+    JSON_ASSERT(static_cast<std::size_t>(last - first) >= len);
+    std::memcpy(first, buf.data(), len);
+    return first + len;
+}
+
 }  // namespace dtoa_impl

 /*!
@@ -24972,7 +25453,6 @@ JSON_HEDLEY_NON_NULL(1, 2)
 JSON_HEDLEY_RETURNS_NON_NULL
 char* to_chars(char* first, const char* last, FloatType value)
 {
-    static_cast<void>(last); // maybe unused - fix warning
    JSON_ASSERT(std::isfinite(value));

    // Use signbit(value) instead of (value < 0) since signbit works for -0.
@@ -24998,28 +25478,7 @@ char* to_chars(char* first, const char* last, FloatType value)
    JSON_HEDLEY_DIAGNOSTIC_POP
 #endif

-    JSON_ASSERT(last - first >= std::numeric_limits<FloatType>::max_digits10);
-
-    // Compute v = buffer * 10^decimal_exponent.
-    // The decimal digits are stored in the buffer, which needs to be interpreted
-    // as an unsigned decimal integer.
-    // len is the length of the buffer, i.e., the number of decimal digits.
-    int len = 0;
-    int decimal_exponent = 0;
-    dtoa_impl::grisu2(first, len, decimal_exponent, value);
-
-    JSON_ASSERT(len <= std::numeric_limits<FloatType>::max_digits10);
-
-    // Format the buffer like printf("%.*g", prec, value)
-    constexpr int kMinExp = -4;
-    // Use digits10 here to increase compatibility with version 2.
-    constexpr int kMaxExp = std::numeric_limits<FloatType>::digits10;
-
-    JSON_ASSERT(last - first >= kMaxExp + 2);
-    JSON_ASSERT(last - first >= 2 + (-kMinExp - 1) + std::numeric_limits<FloatType>::max_digits10);
-    JSON_ASSERT(last - first >= std::numeric_limits<FloatType>::max_digits10 + 6);
-
-    return dtoa_impl::format_buffer(first, len, decimal_exponent, kMinExp, kMaxExp);
+    return dtoa_impl::write_positive(first, last, value);
 }

 }  // namespace detail
@@ -33,7 +33,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);
@@ -46,16 +46,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")
@@ -15,6 +15,20 @@
 #include <nlohmann/json.hpp>
 using nlohmann::detail::dtoa_impl::reinterpret_bits;

+#include <array>
+#include <cmath>
+#include <cstdint>
+#include <cstdio>
+#include <cstdlib>
+#include <limits>
+#include <random>
+#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 +464,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 +528,243 @@ 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);
+    return std::strtod(text.c_str(), nullptr) == 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
+        std::array<char, 64> shorter{};
+        const int n = std::snprintf(shorter.data(), shorter.size(), "%.*e", static_cast<int>(digits.size()) - 2, v); // NOLINT(cppcoreguidelines-pro-type-vararg,hicpp-vararg)
+        const auto near = digits_and_exponent(std::string(shorter.data(), static_cast<std::size_t>(n)));
+        // 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) && v != 0)
+            {
+                check_shortest(v);
+            }
+        }
+    }
+}