From 38e678d9f1b6247101cd4b719f3a9a4c42d3726f Mon Sep 17 00:00:00 2001 From: abolz Date: Sat, 9 Dec 2017 13:14:13 +0100 Subject: [PATCH] Use Grisu2 for floating-point formatting This is an attempt to fix #360. The algorithm produces decimal representations which are guaranteed to roundtrip and in ~99.8% actually produces the shortest possible representation. So this is a nice compromise between using a precision of digits10 and max_digits10. Note 1: The implementation only works for IEEE single/double precision numbers. So the old implementation is kept for compatibility with non-IEEE implementations and 'long double'. Note 2: If number_float_t is 'float', not all serialized numbers can be recovered using strtod (strtof works, though). (There is exactly one such number and the result is off by 1 ulp.) This can be avoided by changing the implementation (the fix is trivial), but then the resulting decimal numbers are not exactly short. Note 3: The code should be replaced by std::to_chars once it is available. --- src/json.hpp | 921 ++++++++++++++++++++++++++++++++++++++++- test/src/unit-dtoa.cpp | 470 +++++++++++++++++++++ 2 files changed, 1390 insertions(+), 1 deletion(-) create mode 100644 test/src/unit-dtoa.cpp diff --git a/src/json.hpp b/src/json.hpp index 23d420547..7112c2bab 100644 --- a/src/json.hpp +++ b/src/json.hpp @@ -6096,6 +6096,902 @@ class binary_writer // serialization // /////////////////// +// Implements the Grisu2 algorithm for binary to decimal floating-point +// conversion. +// +// This implementation is a slightly modified version of the reference +// implementation by Florian Loitsch which can be obtained from +// http://florian.loitsch.com/publications (bench.tar.gz) +// +// References: +// +// [1] Loitsch, "Printing Floating-Point Numbers Quickly and Accurately with Integers", +// Proceedings of the ACM SIGPLAN 2010 Conference on Programming Language Design and Implementation, PLDI 2010 +// [2] Burger, Dybvig, "Printing Floating-Point Numbers Quickly and Accurately", +// Proceedings of the ACM SIGPLAN 1996 Conference on Programming Language Design and Implementation, PLDI 1996 +// +// Note: +// dtoa::to_short_string should be replaced with std::to_chars if available... +namespace dtoa { + +template +inline Target reinterpret_bits(Source source) +{ + static_assert(sizeof(Target) == sizeof(Source), "Size mismatch"); + + Target target; + std::memcpy(&target, &source, sizeof(Source)); + return target; +} + +struct diyfp // f * 2^e +{ + static constexpr int kPrecision = 64; // = q + + uint64_t f; + int e; + + constexpr diyfp() : f(0), e(0) {} + constexpr diyfp(uint64_t f_, int e_) : f(f_), e(e_) {} + + // Returns x - y. + // PRE: x.e == y.e and x.f >= y.f + static diyfp sub(diyfp x, diyfp y); + + // Returns x * y. + // The result is rounded. (Only the upper q bits are returned.) + static diyfp mul(diyfp x, diyfp y); + + // Normalize x such that the significand is >= 2^(q-1). + // PRE: x.f != 0 + static diyfp normalize(diyfp x); + + // Normalize x such that the result has the exponent E. + // PRE: e >= x.e and the upper e - x.e bits of x.f must be zero. + static diyfp normalize_to(diyfp x, int e); +}; + +inline diyfp diyfp::sub(diyfp x, diyfp y) +{ + assert(x.e == y.e); + assert(x.f >= y.f); + + return diyfp(x.f - y.f, x.e); +} + +inline diyfp diyfp::mul(diyfp x, diyfp y) +{ + // Computes: + // f = round((x.f * y.f) / 2^q) + // e = x.e + y.e + q + + // Emulate the 64-bit * 64-bit multiplication: + // + // p = u * v + // = (u_lo + 2^32 u_hi) (v_lo + 2^32 v_hi) + // = (u_lo v_lo ) + 2^32 ((u_lo v_hi ) + (u_hi v_lo )) + 2^64 (u_hi v_hi ) + // = (p0 ) + 2^32 ((p1 ) + (p2 )) + 2^64 (p3 ) + // = (p0_lo + 2^32 p0_hi) + 2^32 ((p1_lo + 2^32 p1_hi) + (p2_lo + 2^32 p2_hi)) + 2^64 (p3 ) + // = (p0_lo ) + 2^32 (p0_hi + p1_lo + p2_lo ) + 2^64 (p1_hi + p2_hi + p3) + // = (p0_lo ) + 2^32 (Q ) + 2^64 (H ) + // = (p0_lo ) + 2^32 (Q_lo + 2^32 Q_hi ) + 2^64 (H ) + // + // (Since Q might be larger than 2^32 - 1) + // + // = (p0_lo + 2^32 Q_lo) + 2^64 (Q_hi + H) + // + // (Q_hi + H does not overflow a 64-bit int) + // + // = p_lo + 2^64 p_hi + + const uint64_t u_lo = x.f & 0xFFFFFFFF; + const uint64_t u_hi = x.f >> 32; + const uint64_t v_lo = y.f & 0xFFFFFFFF; + const uint64_t v_hi = y.f >> 32; + + const uint64_t p0 = u_lo * v_lo; + const uint64_t p1 = u_lo * v_hi; + const uint64_t p2 = u_hi * v_lo; + const uint64_t p3 = u_hi * v_hi; + + const uint64_t p0_hi = p0 >> 32; + const uint64_t p1_lo = p1 & 0xFFFFFFFF; + const uint64_t p1_hi = p1 >> 32; + const uint64_t p2_lo = p2 & 0xFFFFFFFF; + const uint64_t p2_hi = p2 >> 32; + + uint64_t Q = p0_hi + p1_lo + p2_lo; + + // The full product might now be computed as + // + // p_hi = p3 + p2_hi + p1_hi + (Q >> 32) + // p_lo = p0_lo + (Q << 32) + // + // But in this particular case here, the full p_lo is not required. + // Effectively we only need to add the highest bit in p_lo to p_hi (and + // Q_hi + 1 does not overflow). + + Q += uint64_t{1} << (63 - 32); // round, ties up + + const uint64_t h = p3 + p2_hi + p1_hi + (Q >> 32); + + return diyfp(h, x.e + y.e + 64); +} + +inline diyfp diyfp::normalize(diyfp x) +{ + assert(x.f != 0); + + while ((x.f >> 63) == 0) + { + x.f <<= 1; + x.e--; + } + + return x; +} + +inline diyfp diyfp::normalize_to(diyfp x, int target_exponent) +{ + const int delta = x.e - target_exponent; + + assert(delta >= 0); + assert(((x.f << delta) >> delta) == x.f); + + return diyfp(x.f << delta, target_exponent); +} + +struct boundaries { + diyfp w; + diyfp minus; + diyfp plus; +}; + +// Compute the (normalized) diyfp representing the input number 'value' and its boundaries. +// PRE: value must be finite and positive +template +inline boundaries compute_boundaries(FloatType value) +{ + assert(std::isfinite(value)); + assert(value > 0); + + // Convert the IEEE representation into a diyfp. + // + // If v is denormal: + // value = 0.F * 2^(1 - bias) = ( F) * 2^(1 - bias - (p-1)) + // If v is normalized: + // value = 1.F * 2^(E - bias) = (2^(p-1) + F) * 2^(E - bias - (p-1)) + + static_assert(std::numeric_limits::is_iec559, "Requires an IEEE-754 floating-point implementation"); + + constexpr int kPrecision = std::numeric_limits::digits; // = p (includes the hidden bit) + constexpr int kBias = std::numeric_limits::max_exponent - 1 + (kPrecision - 1); + constexpr int kMinExp = 1 - kBias; + constexpr uint64_t kHiddenBit = uint64_t{1} << (kPrecision - 1); // = 2^(p-1) + + using bits_type = typename std::conditional::type; + + const uint64_t bits = reinterpret_bits(value); + const uint64_t E = bits >> (kPrecision - 1); + const uint64_t F = bits & (kHiddenBit - 1); + + const bool is_denormal = (E == 0); + + const diyfp v + = is_denormal + ? diyfp(F, 1 - kBias) + : diyfp(F + kHiddenBit, static_cast(E) - kBias); + + // Compute the boundaries m- and m+ of the floating-point value + // v = f * 2^e. + // + // Determine v- and v+, the floating-point predecessor and successor if v, + // respectively. + // + // v- = v - 2^e if f != 2^(p-1) or e == e_min (A) + // = v - 2^(e-1) if f == 2^(p-1) and e > e_min (B) + // + // v+ = v + 2^e + // + // Let m- = (v- + v) / 2 and m+ = (v + v+) / 2. All real numbers _strictly_ + // between m- and m+ round to v, regardless of how the input rounding + // algorithm breaks ties. + // + // ---+-------------+-------------+-------------+-------------+--- (A) + // v- m- v m+ v+ + // + // -----------------+------+------+-------------+-------------+--- (B) + // v- m- v m+ v+ + +// const bool lower_boundary_is_closer = (F == 0 and E > 1); + const bool lower_boundary_is_closer = (v.f == kHiddenBit and v.e > kMinExp); + + const diyfp m_plus = diyfp(2*v.f + 1, v.e - 1); + const diyfp m_minus + = lower_boundary_is_closer + ? diyfp(4*v.f - 1, v.e - 2) // (B) + : diyfp(2*v.f - 1, v.e - 1); // (A) + + // Determine the normalized w+ = m+. + const diyfp w_plus = diyfp::normalize(m_plus); + + // Determine w- = m- such that e_(w-) = e_(w+). + const diyfp w_minus = diyfp::normalize_to(m_minus, w_plus.e); + + return {diyfp::normalize(v), w_minus, w_plus}; +} + +// Given normalized diyfp w, Grisu needs to find a (normalized) cached +// power-of-ten c, such that the exponent of the product c * w = f * 2^e lies +// within a certain range [alpha, gamma] (Definition 3.2 from [1]) +// +// alpha <= e = e_c + e_w + q <= gamma +// +// The choice of (alpha,gamma) determines the size of the table and the digit +// generation procedure. Using (alpha,gamma)=(-60,-32) works out well in +// practice. +constexpr int kAlpha = -60; +constexpr int kGamma = -32; + +struct cached_power { // c = f * 2^e ~= 10^k + uint64_t f; + int e; + int k; +}; + +inline cached_power get_cached_power_for_binary_exponent(int e) +{ + // NB: + // Actually this function returns c, such that -60 <= e_c + e + 64 <= -34. + + constexpr int kCachedPowersSize = 79; + constexpr int kCachedPowersMinDecExp = -300; + constexpr int kCachedPowersDecStep = 8; + + static constexpr cached_power kCachedPowers[] = { + { 0xAB70FE17C79AC6CA, -1060, -300 }, + { 0xFF77B1FCBEBCDC4F, -1034, -292 }, + { 0xBE5691EF416BD60C, -1007, -284 }, + { 0x8DD01FAD907FFC3C, -980, -276 }, + { 0xD3515C2831559A83, -954, -268 }, + { 0x9D71AC8FADA6C9B5, -927, -260 }, + { 0xEA9C227723EE8BCB, -901, -252 }, + { 0xAECC49914078536D, -874, -244 }, + { 0x823C12795DB6CE57, -847, -236 }, + { 0xC21094364DFB5637, -821, -228 }, + { 0x9096EA6F3848984F, -794, -220 }, + { 0xD77485CB25823AC7, -768, -212 }, + { 0xA086CFCD97BF97F4, -741, -204 }, + { 0xEF340A98172AACE5, -715, -196 }, + { 0xB23867FB2A35B28E, -688, -188 }, + { 0x84C8D4DFD2C63F3B, -661, -180 }, + { 0xC5DD44271AD3CDBA, -635, -172 }, + { 0x936B9FCEBB25C996, -608, -164 }, + { 0xDBAC6C247D62A584, -582, -156 }, + { 0xA3AB66580D5FDAF6, -555, -148 }, + { 0xF3E2F893DEC3F126, -529, -140 }, + { 0xB5B5ADA8AAFF80B8, -502, -132 }, + { 0x87625F056C7C4A8B, -475, -124 }, + { 0xC9BCFF6034C13053, -449, -116 }, + { 0x964E858C91BA2655, -422, -108 }, + { 0xDFF9772470297EBD, -396, -100 }, + { 0xA6DFBD9FB8E5B88F, -369, -92 }, + { 0xF8A95FCF88747D94, -343, -84 }, + { 0xB94470938FA89BCF, -316, -76 }, + { 0x8A08F0F8BF0F156B, -289, -68 }, + { 0xCDB02555653131B6, -263, -60 }, + { 0x993FE2C6D07B7FAC, -236, -52 }, + { 0xE45C10C42A2B3B06, -210, -44 }, + { 0xAA242499697392D3, -183, -36 }, + { 0xFD87B5F28300CA0E, -157, -28 }, + { 0xBCE5086492111AEB, -130, -20 }, + { 0x8CBCCC096F5088CC, -103, -12 }, + { 0xD1B71758E219652C, -77, -4 }, + { 0x9C40000000000000, -50, 4 }, + { 0xE8D4A51000000000, -24, 12 }, + { 0xAD78EBC5AC620000, 3, 20 }, + { 0x813F3978F8940984, 30, 28 }, + { 0xC097CE7BC90715B3, 56, 36 }, + { 0x8F7E32CE7BEA5C70, 83, 44 }, + { 0xD5D238A4ABE98068, 109, 52 }, + { 0x9F4F2726179A2245, 136, 60 }, + { 0xED63A231D4C4FB27, 162, 68 }, + { 0xB0DE65388CC8ADA8, 189, 76 }, + { 0x83C7088E1AAB65DB, 216, 84 }, + { 0xC45D1DF942711D9A, 242, 92 }, + { 0x924D692CA61BE758, 269, 100 }, + { 0xDA01EE641A708DEA, 295, 108 }, + { 0xA26DA3999AEF774A, 322, 116 }, + { 0xF209787BB47D6B85, 348, 124 }, + { 0xB454E4A179DD1877, 375, 132 }, + { 0x865B86925B9BC5C2, 402, 140 }, + { 0xC83553C5C8965D3D, 428, 148 }, + { 0x952AB45CFA97A0B3, 455, 156 }, + { 0xDE469FBD99A05FE3, 481, 164 }, + { 0xA59BC234DB398C25, 508, 172 }, + { 0xF6C69A72A3989F5C, 534, 180 }, + { 0xB7DCBF5354E9BECE, 561, 188 }, + { 0x88FCF317F22241E2, 588, 196 }, + { 0xCC20CE9BD35C78A5, 614, 204 }, + { 0x98165AF37B2153DF, 641, 212 }, + { 0xE2A0B5DC971F303A, 667, 220 }, + { 0xA8D9D1535CE3B396, 694, 228 }, + { 0xFB9B7CD9A4A7443C, 720, 236 }, + { 0xBB764C4CA7A44410, 747, 244 }, + { 0x8BAB8EEFB6409C1A, 774, 252 }, + { 0xD01FEF10A657842C, 800, 260 }, + { 0x9B10A4E5E9913129, 827, 268 }, + { 0xE7109BFBA19C0C9D, 853, 276 }, + { 0xAC2820D9623BF429, 880, 284 }, + { 0x80444B5E7AA7CF85, 907, 292 }, + { 0xBF21E44003ACDD2D, 933, 300 }, + { 0x8E679C2F5E44FF8F, 960, 308 }, + { 0xD433179D9C8CB841, 986, 316 }, + { 0x9E19DB92B4E31BA9, 1013, 324 }, + }; + + // This computation gives exactly the same results for k as + // k = ceil((kAlpha - e - 1) * 0.30102999566398114) + // for |e| <= 1500, but doesn't require floating-point operations. + // NB: log_10(2) ~= 78913 / 2^18 + assert(e >= -1500); + assert(e <= 1500); + const int f = kAlpha - e - 1; + const int k = (f * 78913) / (1 << 18) + (f > 0); + + const int index = (-kCachedPowersMinDecExp + k + (kCachedPowersDecStep - 1)) / kCachedPowersDecStep; + assert(index >= 0); + assert(index < kCachedPowersSize); + static_cast(kCachedPowersSize); // Fix warning. + + const cached_power cached = kCachedPowers[index]; + assert(kAlpha <= cached.e + e + 64); + assert(kGamma >= cached.e + e + 64); + + return cached; +} + +// For n != 0, returns k, such that pow10 := 10^(k-1) <= n < 10^k. +// For n == 0, returns 1 and sets pow10 := 1. +inline int find_largest_pow10(uint32_t n, uint32_t& pow10) +{ + if (n >= 1000000000) { pow10 = 1000000000; return 10; } + if (n >= 100000000) { pow10 = 100000000; return 9; } + if (n >= 10000000) { pow10 = 10000000; return 8; } + if (n >= 1000000) { pow10 = 1000000; return 7; } + if (n >= 100000) { pow10 = 100000; return 6; } + if (n >= 10000) { pow10 = 10000; return 5; } + if (n >= 1000) { pow10 = 1000; return 4; } + if (n >= 100) { pow10 = 100; return 3; } + if (n >= 10) { pow10 = 10; return 2; } + + pow10 = 1; return 1; +} + +inline void grisu2_round(char* buf, int len, uint64_t dist, uint64_t delta, uint64_t rest, uint64_t ten_k) +{ + assert(len >= 1); + assert(dist <= delta); + assert(rest <= delta); + assert(ten_k > 0); + + // <--------------------------- delta ----> + // <---- dist ---------> + // --------------[------------------+-------------------]-------------- + // w- w w+ + // + // ten_k + // <------> + // <---- rest ----> + // --------------[------------------+----+--------------]-------------- + // w V + // = buf * 10^k + // + // ten_k represents a unit-in-the-last-place in the decimal representation + // stored in buf. + // Decrement buf by ten_k while this takes buf closer to w. + + // The tests are written in this order to avoid overflow in unsigned + // integer arithmetic. + + while (rest < dist + and delta - rest >= ten_k + and (rest + ten_k < dist or dist - rest > rest + ten_k - dist)) + { + assert(buf[len - 1] != '0'); + buf[len - 1]--; + rest += ten_k; + } +} + +inline void grisu2_digit_gen(char* buffer, int& length, int& decimal_exponent, diyfp M_minus, diyfp w, diyfp M_plus) +{ + static_assert(diyfp::kPrecision == 64, "invalid config"); + static_assert(kAlpha >= -60, "invalid config"); + static_assert(kGamma <= -32, "invalid config"); + + // Generates the digits (and the exponent) of a decimal floating-point + // number V = buffer * 10^decimal_exponent in the range [M-, M+]. The diyfp's + // w, M- and M+ share the same exponent e, which satisfies alpha <= e <= gamma. + // + // <--------------------------- delta ----> + // <---- dist ---------> + // --------------[------------------+-------------------]-------------- + // M- w M+ + // + // Grisu2 generates the digits of M+ from left to right and stops as soon as + // V is in [M-,M+]. + + assert(M_plus.e >= kAlpha); + assert(M_plus.e <= kGamma); + + uint64_t delta = diyfp::sub(M_plus, M_minus).f; // (significand of (M+ - M-), implicit exponent is e) + uint64_t dist = diyfp::sub(M_plus, w ).f; // (significand of (M+ - w ), implicit exponent is e) + + // Split M+ = f * 2^e into two parts p1 and p2 (note: e < 0): + // + // M+ = f * 2^e + // = ((f div 2^-e) * 2^-e + (f mod 2^-e)) * 2^e + // = ((p1 ) * 2^-e + (p2 )) * 2^e + // = p1 + p2 * 2^e + + const diyfp one(uint64_t{1} << -M_plus.e, M_plus.e); + + uint32_t p1 = static_cast(M_plus.f >> -one.e); // p1 = f div 2^-e (Since -e >= 32, p1 fits into a 32-bit int.) + uint64_t p2 = M_plus.f & (one.f - 1); // p2 = f mod 2^-e + + // 1) + // + // Generate the digits of the integral part p1 = d[n-1]...d[1]d[0] + + assert(p1 > 0); + + uint32_t pow10; + const int k = find_largest_pow10(p1, pow10); + + // 10^(k-1) <= p1 < 10^k, pow10 = 10^(k-1) + // + // p1 = (p1 div 10^(k-1)) * 10^(k-1) + (p1 mod 10^(k-1)) + // = (d[k-1] ) * 10^(k-1) + (p1 mod 10^(k-1)) + // + // M+ = p1 + p2 * 2^e + // = d[k-1] * 10^(k-1) + (p1 mod 10^(k-1)) + p2 * 2^e + // = d[k-1] * 10^(k-1) + ((p1 mod 10^(k-1)) * 2^-e + p2) * 2^e + // = d[k-1] * 10^(k-1) + ( rest) * 2^e + // + // Now generate the digits d[n] of p1 from left to right (n = k-1,...,0) + // + // p1 = d[k-1]...d[n] * 10^n + d[n-1]...d[0] + // + // but stop as soon as + // + // rest * 2^e = (d[n-1]...d[0] * 2^-e + p2) * 2^e <= delta * 2^e + + int n = k; + while (n > 0) + { + // Invariants: + // M+ = buffer * 10^n + (p1 + p2 * 2^e) (buffer = 0 for n = k) + // pow10 = 10^(n-1) <= p1 < 10^n + // + const uint32_t d = p1 / pow10; // d = p1 div 10^(n-1) + const uint32_t r = p1 % pow10; // r = p1 mod 10^(n-1) + // + // M+ = buffer * 10^n + (d * 10^(n-1) + r) + p2 * 2^e + // = (buffer * 10 + d) * 10^(n-1) + (r + p2 * 2^e) + // + assert(d <= 9); + buffer[length++] = static_cast('0' + d); // buffer := buffer * 10 + d + // + // M+ = buffer * 10^(n-1) + (r + p2 * 2^e) + // + p1 = r; + n--; + // + // M+ = buffer * 10^n + (p1 + p2 * 2^e) + // pow10 = 10^n + // + + // Now check if enough digits have been generated. + // Compute + // + // p1 + p2 * 2^e = (p1 * 2^-e + p2) * 2^e = rest * 2^e + // + // Note: + // Since rest and delta share the same exponent e, it suffices to + // compare the significands. + const uint64_t rest = (uint64_t{p1} << -one.e) + p2; + if (rest <= delta) + { + // V = buffer * 10^n, with M- <= V <= M+. + + decimal_exponent += n; + + // We may now just stop. But instead look if the buffer could be + // decremented to bring V closer to w. + // + // pow10 = 10^n is now 1 ulp in the decimal representation V. + // The rounding procedure works with diyfp's with an implicit + // exponent of e. + // + // 10^n = (10^n * 2^-e) * 2^e = ulp * 2^e + // + const uint64_t ten_n = uint64_t{pow10} << -one.e; + grisu2_round(buffer, length, dist, delta, rest, ten_n); + + return; + } + + pow10 /= 10; + // + // pow10 = 10^(n-1) <= p1 < 10^n + // Invariants restored. + } + + // 2) + // + // The digits of the integral part have been generated: + // + // M+ = d[k-1]...d[1]d[0] + p2 * 2^e + // = buffer + p2 * 2^e + // + // Now generate the digits of the fractional part p2 * 2^e. + // + // Note: + // No decimal point is generated: the exponent is adjusted instead. + // + // p2 actually represents the fraction + // + // p2 * 2^e + // = p2 / 2^-e + // = d[-1] / 10^1 + d[-2] / 10^2 + ... + // + // Now generate the digits d[-m] of p1 from left to right (m = 1,2,...) + // + // p2 * 2^e = d[-1]d[-2]...d[-m] * 10^-m + // + 10^-m * (d[-m-1] / 10^1 + d[-m-2] / 10^2 + ...) + // + // using + // + // 10^m * p2 = ((10^m * p2) div 2^-e) * 2^-e + ((10^m * p2) mod 2^-e) + // = ( d) * 2^-e + ( r) + // + // or + // 10^m * p2 * 2^e = d + r * 2^e + // + // i.e. + // + // M+ = buffer + p2 * 2^e + // = buffer + 10^-m * (d + r * 2^e) + // = (buffer * 10^m + d) * 10^-m + 10^-m * r * 2^e + // + // and stop as soon as 10^-m * r * 2^e <= delta * 2^e + + assert(p2 > delta); + + int m = 0; + for (;;) + { + // Invariant: + // M+ = buffer * 10^-m + 10^-m * (d[-m-1] / 10 + d[-m-2] / 10^2 + ...) * 2^e + // = buffer * 10^-m + 10^-m * (p2 ) * 2^e + // = buffer * 10^-m + 10^-m * (1/10 * (10 * p2) ) * 2^e + // = buffer * 10^-m + 10^-m * (1/10 * ((10*p2 div 2^-e) * 2^-e + (10*p2 mod 2^-e)) * 2^e + // + assert(p2 <= UINT64_MAX / 10); + p2 *= 10; + const uint64_t d = p2 >> -one.e; // d = (10 * p2) div 2^-e + const uint64_t r = p2 & (one.f - 1); // r = (10 * p2) mod 2^-e + // + // M+ = buffer * 10^-m + 10^-m * (1/10 * (d * 2^-e + r) * 2^e + // = buffer * 10^-m + 10^-m * (1/10 * (d + r * 2^e)) + // = (buffer * 10 + d) * 10^(-m-1) + 10^(-m-1) * r * 2^e + // + assert(d <= 9); + buffer[length++] = static_cast('0' + d); // buffer := buffer * 10 + d + // + // M+ = buffer * 10^(-m-1) + 10^(-m-1) * r * 2^e + // + p2 = r; + m++; + // + // M+ = buffer * 10^-m + 10^-m * p2 * 2^e + // Invariant restored. + + // Check if enough digits have been generated. + // + // 10^-m * p2 * 2^e <= delta * 2^e + // p2 * 2^e <= 10^m * delta * 2^e + // p2 <= 10^m * delta + delta *= 10; + dist *= 10; + if (p2 <= delta) + { + break; + } + } + + // V = buffer * 10^-m, with M- <= V <= M+. + + decimal_exponent -= m; + + // 1 ulp in the decimal representation is now 10^-m. + // Since delta and dist are now scaled by 10^m, we need to do the + // same with ulp in order to keep the units in sync. + // + // 10^m * 10^-m = 1 = 2^-e * 2^e = ten_m * 2^e + // + const uint64_t ten_m = one.f; + grisu2_round(buffer, length, dist, delta, p2, ten_m); + + // By construction this algorithm generates the shortest possible decimal + // number (Loitsch, Theorem 6.2) which rounds back to w. + // For an input number of precision p, at least + // + // N = 1 + ceil(p * log_10(2)) + // + // decimal digits are sufficient to identify all binary floating-point + // numbers (Matula, "In-and-Out conversions"). + // This implies that the algorithm does not produce more than N decimal + // digits. + // + // N = 17 for p = 53 (IEEE double precision) + // N = 9 for p = 24 (IEEE single precision) +} + +// v = buf * 10^decimal_exponent +// len is the length of the buffer (number of decimal digits) +// The buffer must be large enough, i.e. >= max_digits10. +inline void grisu2(char* buf, int& len, int& decimal_exponent, diyfp m_minus, diyfp v, diyfp m_plus) +{ + assert(m_plus.e == m_minus.e); + assert(m_plus.e == v.e); + + // --------(-----------------------+-----------------------)-------- (A) + // m- v m+ + // + // --------------------(-----------+-----------------------)-------- (B) + // m- v m+ + // + // First scale v (and m- and m+) such that the exponent is in the range + // [alpha, gamma]. + + const cached_power cached = get_cached_power_for_binary_exponent(m_plus.e); + + const diyfp c_minus_k(cached.f, cached.e); // = c ~= 10^-k + + // The exponent of the products is = v.e + c_minus_k.e + q and is in the range [alpha,gamma] + const diyfp w = diyfp::mul(v, c_minus_k); + const diyfp w_minus = diyfp::mul(m_minus, c_minus_k); + const diyfp w_plus = diyfp::mul(m_plus, c_minus_k); + + // ----(---+---)---------------(---+---)---------------(---+---)---- + // w- w w+ + // = c*m- = c*v = c*m+ + // + // diyfp::mul rounds its result and c_minus_k is approximated too. w, w- and + // w+ are now off by a small amount. + // In fact: + // + // w - v * 10^k < 1 ulp + // + // To account for this inaccuracy, add resp. subtract 1 ulp. + // + // --------+---[---------------(---+---)---------------]---+-------- + // w- M- w M+ w+ + // + // Now any number in [M-, M+] (bounds included) will round to w when input, + // regardless of how the input rounding algorithm breaks ties. + // + // And digit_gen generates the shortest possible such number in [M-, M+]. + // Note that this does not mean that Grisu2 always generates the shortest + // possible number in the interval (m-, m+). + const diyfp M_minus(w_minus.f + 1, w_minus.e); + const diyfp M_plus (w_plus.f - 1, w_plus.e ); + + decimal_exponent = -cached.k; // = -(-k) = k + + grisu2_digit_gen(buf, len, decimal_exponent, M_minus, w, M_plus); +} + +// v = buf * 10^decimal_exponent +// len is the length of the buffer (number of decimal digits) +// The buffer must be large enough, i.e. >= max_digits10. +template +inline void grisu2(char* buf, int& len, int& decimal_exponent, FloatType value) +{ + assert(std::isfinite(value)); + assert(value > 0); + +#if 0 + // If the neighbors (and boundaries) of 'value' are always computed for + // double-precision numbers, all float's can be recovered using strtod + // (and strtof). However, the resulting decimal representations are not + // exactly "short". + // + // If the neighbors are computed for single-precision numbers, there is a + // single float (7.0385307e-26f) which can't be recovered using strtod. + // (The resulting double precision is off by 1 ulp.) + const boundaries w = compute_boundaries(static_cast(value)); +#else + const boundaries w = compute_boundaries(value); +#endif + + grisu2(buf, len, decimal_exponent, w.minus, w.w, w.plus); +} + +// Appends a decimal representation of e to buf. +// Returns a pointer to the element following the exponent. +// PRE: -1000 < e < 1000 +inline char* append_exponent(char* buf, int e) +{ + assert(e > -1000); + assert(e < 1000); + + if (e < 0) + { + e = -e; + *buf++ = '-'; + } + else + { + *buf++ = '+'; + } + + uint32_t k = static_cast(e); + if (k < 10) + { + *buf++ = static_cast('0' + k); + } + else if (k < 100) + { + *buf++ = static_cast('0' + k / 10); k %= 10; + *buf++ = static_cast('0' + k ); + } + else + { + *buf++ = static_cast('0' + k / 100); k %= 100; + *buf++ = static_cast('0' + k / 10); k %= 10; + *buf++ = static_cast('0' + k ); + } + + return buf; +} + +// Prettify v = buf * 10^decimal_exponent +// If v is in the range [10^min_exp, 10^max_exp) it will be printed in fixed-point notation. +// Otherwise it will be printed in exponential notation. +// PRE: min_exp < 0 +// PRE: max_exp > 0 +inline char* format_buffer(char* buf, int len, int decimal_exponent, int min_exp, int max_exp) +{ + assert(min_exp < 0); + assert(max_exp > 0); + + int k = len; + int n = len + decimal_exponent; + // v = buf * 10^(n-k) + // k is the length of the buffer (number of decimal digits) + // n is the position of the decimal point relative to the start of the buffer. + + if (k <= n and n <= max_exp) + { + // digits[000] + // len <= max_exp + + std::memset(buf + k, '0', static_cast(n - k)); + buf[n++] = '.'; + buf[n++] = '0'; + return buf + n; + } + + if (0 < n and n <= max_exp) + { + // dig.its + // len <= max_digits10 + 1 + + assert(k > n); + + std::memmove(buf + (n + 1), buf + n, static_cast(k - n)); + buf[n] = '.'; + return buf + (k + 1); + } + + if (min_exp < n and n <= 0) + { + // 0.[000]digits + // len <= 2 + (-min_exp - 1) + max_digits10 + + std::memmove(buf + (2 + -n), buf, static_cast(k)); + buf[0] = '0'; + buf[1] = '.'; + std::memset(buf + 2, '0', static_cast(-n)); + return buf + (2 + (-n) + k); + } + + if (k == 1) + { + // dE+123 + // len <= 1 + 5 + + buf += 1; + } + else + { + // d.igitsE+123 + // len <= max_digits10 + 1 + 5 + + std::memmove(buf + 2, buf + 1, static_cast(k - 1)); + buf[1] = '.'; + buf += 1 + k; + } + + *buf++ = 'e'; + return append_exponent(buf, n - 1); +} + +// Generates a decimal representation of the floating-point number 'value' in +// [first, last). +// +// The input number must be finite, i.e. NaN's and Inf's are not supported. +// The buffer must be large enough. +// +// The result is _not_ null-terminated. +template +inline char* to_short_string(char* first, char* last, FloatType value) +{ + static_assert(diyfp::kPrecision >= std::numeric_limits::digits + 3, "Not enough precision"); + + static_cast(last); // maybe unused - fix warning + assert(!std::isnan(value)); + assert(!std::isinf(value)); + + // Use signbit(value) instead of (value < 0) since signbit works for -0. + if (std::signbit(value)) + { + value = -value; + *first++ = '-'; + } + + if (value == 0) + { + *first++ = '0'; + *first++ = '.'; + *first++ = '0'; + return first; + } + + assert(last - first >= std::numeric_limits::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; + grisu2(first, len, decimal_exponent, value); + + assert(len <= std::numeric_limits::max_digits10); + + // Format the buffer like printf("%.*g", prec, value) + constexpr int kMinExp = -4; +#if 0 + constexpr int kMaxExp = std::numeric_limits::digits10; +#else + constexpr int kMaxExp = std::numeric_limits::max_digits10; +#endif + + assert(last - first >= + std::max({ + kMaxExp, + 2 + (-kMinExp - 1) + std::numeric_limits::max_digits10, + std::numeric_limits::max_digits10 + 6})); + + return format_buffer(first, len, decimal_exponent, kMinExp, kMaxExp); +} + +} // namespace dtoa + template class serializer { @@ -6711,12 +7607,35 @@ class serializer void dump_float(number_float_t x) { // NaN / inf - if (not std::isfinite(x) or std::isnan(x)) + if (not std::isfinite(x)) { o->write_characters("null", 4); return; } + // If number_float_t is an IEEE-754 single or double precision number, + // use the Grisu2 algorithm to produce short numbers which are guaranteed + // to round-trip. + // + // NB: The test below works if == . + static constexpr bool is_ieee_single_or_double + = (std::numeric_limits::is_iec559 and std::numeric_limits::digits == 24 and std::numeric_limits::max_exponent == 128) or + (std::numeric_limits::is_iec559 and std::numeric_limits::digits == 53 and std::numeric_limits::max_exponent == 1024); + + dump_float(x, std::integral_constant()); + } + + private: + void dump_float(number_float_t x, std::true_type) + { + char* begin = number_buffer.data(); + char* end = dtoa::to_short_string(begin, begin + number_buffer.size(), x); + + o->write_characters(begin, static_cast(end - begin)); + } + + void dump_float(number_float_t x, std::false_type) + { // get number of digits for a text -> float -> text round-trip static constexpr auto d = std::numeric_limits::digits10; diff --git a/test/src/unit-dtoa.cpp b/test/src/unit-dtoa.cpp new file mode 100644 index 000000000..209915fb2 --- /dev/null +++ b/test/src/unit-dtoa.cpp @@ -0,0 +1,470 @@ +/* + __ _____ _____ _____ + __| | __| | | | JSON for Modern C++ (test suite) +| | |__ | | | | | | version 2.1.1 +|_____|_____|_____|_|___| https://github.com/nlohmann/json + +Licensed under the MIT License . +Copyright (c) 2013-2017 Niels Lohmann . + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +*/ + +// XXX: +// Only compile these tests if 'float' and 'double' are IEEE-754 single- and +// double-precision numbers, resp. + +#include "catch.hpp" + +#include "json.hpp" + +using nlohmann::detail::dtoa::reinterpret_bits; + +static float make_float(uint32_t sign_bit, uint32_t biased_exponent, uint32_t significand) +{ + assert(sign_bit == 0 || sign_bit == 1); + assert(biased_exponent <= 0xFF); + assert(significand <= 0x007FFFFF); + + uint32_t bits = 0; + + bits |= sign_bit << 31; + bits |= biased_exponent << 23; + bits |= significand; + + return reinterpret_bits(bits); +} + +// ldexp -- convert f * 2^e to IEEE single precision +static float make_float(uint64_t f, int e) +{ + constexpr uint64_t kHiddenBit = 0x00800000; + constexpr uint64_t kSignificandMask = 0x007FFFFF; + constexpr int kPhysicalSignificandSize = 23; // Excludes the hidden bit. + constexpr int kExponentBias = 0x7F + kPhysicalSignificandSize; + constexpr int kDenormalExponent = 1 - kExponentBias; + constexpr int kMaxExponent = 0xFF - kExponentBias; + + while (f > kHiddenBit + kSignificandMask) { + f >>= 1; + e++; + } + if (e >= kMaxExponent) { + return std::numeric_limits::infinity(); + } + if (e < kDenormalExponent) { + return 0.0; + } + while (e > kDenormalExponent && (f & kHiddenBit) == 0) { + f <<= 1; + e--; + } + + uint64_t biased_exponent; + if (e == kDenormalExponent && (f & kHiddenBit) == 0) + biased_exponent = 0; + else + biased_exponent = static_cast(e + kExponentBias); + + uint64_t bits = (f & kSignificandMask) | (biased_exponent << kPhysicalSignificandSize); + return reinterpret_bits(static_cast(bits)); +} + +static double make_double(uint64_t sign_bit, uint64_t biased_exponent, uint64_t significand) +{ + assert(sign_bit == 0 || sign_bit == 1); + assert(biased_exponent <= 0x7FF); + assert(significand <= 0x000FFFFFFFFFFFFF); + + uint64_t bits = 0; + + bits |= sign_bit << 63; + bits |= biased_exponent << 52; + bits |= significand; + + return reinterpret_bits(bits); +} + +// ldexp -- convert f * 2^e to IEEE double precision +static double make_double(uint64_t f, int e) +{ + constexpr uint64_t kHiddenBit = 0x0010000000000000; + constexpr uint64_t kSignificandMask = 0x000FFFFFFFFFFFFF; + constexpr int kPhysicalSignificandSize = 52; // Excludes the hidden bit. + constexpr int kExponentBias = 0x3FF + kPhysicalSignificandSize; + constexpr int kDenormalExponent = 1 - kExponentBias; + constexpr int kMaxExponent = 0x7FF - kExponentBias; + + while (f > kHiddenBit + kSignificandMask) { + f >>= 1; + e++; + } + if (e >= kMaxExponent) { + return std::numeric_limits::infinity(); + } + if (e < kDenormalExponent) { + return 0.0; + } + while (e > kDenormalExponent && (f & kHiddenBit) == 0) { + f <<= 1; + e--; + } + + uint64_t biased_exponent; + if (e == kDenormalExponent && (f & kHiddenBit) == 0) + biased_exponent = 0; + else + biased_exponent = static_cast(e + kExponentBias); + + uint64_t bits = (f & kSignificandMask) | (biased_exponent << kPhysicalSignificandSize); + return reinterpret_bits(bits); +} + +TEST_CASE("digit gen") +{ + SECTION("single precision") + { + auto check_float = [](float number, const std::string& digits, int expected_exponent) + { + CAPTURE(number); + CAPTURE(digits); + CAPTURE(expected_exponent); + + char buf[32]; + int len = 0; + int exponent = 0; + nlohmann::detail::dtoa::grisu2(buf, len, exponent, number); + + CHECK(digits == std::string(buf, buf + len)); + CHECK(expected_exponent == exponent); + }; + + check_float(make_float(0, 0, 0x00000001), "1", -45); // min denormal + check_float(make_float(0, 0, 0x007FFFFF), "11754942", -45); // max denormal + check_float(make_float(0, 1, 0x00000000), "11754944", -45); // min normal + check_float(make_float(0, 1, 0x00000001), "11754945", -45); + check_float(make_float(0, 1, 0x007FFFFF), "23509886", -45); + check_float(make_float(0, 2, 0x00000000), "23509887", -45); + check_float(make_float(0, 2, 0x00000001), "2350989", -44); + check_float(make_float(0, 24, 0x00000000), "98607613", -39); // fail if no special case in normalized boundaries + check_float(make_float(0, 30, 0x00000000), "63108872", -37); // fail if no special case in normalized boundaries + check_float(make_float(0, 31, 0x00000000), "12621775", -36); // fail if no special case in normalized boundaries + check_float(make_float(0, 57, 0x00000000), "84703295", -29); // fail if no special case in normalized boundaries + check_float(make_float(0, 254, 0x007FFFFE), "34028233", 31); + check_float(make_float(0, 254, 0x007FFFFF), "34028235", 31); // max normal + + // V. Paxson and W. Kahan, "A Program for Testing IEEE Binary-Decimal Conversion", manuscript, May 1991, + // ftp://ftp.ee.lbl.gov/testbase-report.ps.Z (report) + // ftp://ftp.ee.lbl.gov/testbase.tar.Z (program) + + // Table 16: Stress Inputs for Converting 24-bit Binary to Decimal, < 1/2 ULP + check_float(make_float(12676506, -102), "25", -25); + check_float(make_float(12676506, -103), "125", -26); + check_float(make_float(15445013, 86), "1195", 30); + check_float(make_float(13734123, -138), "39415", -39); + check_float(make_float(12428269, -130), "913085", -38); + check_float(make_float(15334037, -146), "1719005", -43); + check_float(make_float(11518287, -41), "52379105", -13); + check_float(make_float(12584953, -145), "2821644", -43); + check_float(make_float(15961084, -125), "37524328", -38); + check_float(make_float(14915817, -146), "16721209", -44); + check_float(make_float(10845484, -102), "21388946", -31); + check_float(make_float(16431059, -61), "7125836", -18); + + // Table 17: Stress Inputs for Converting 24-bit Binary to Decimal, > 1/2 ULP + check_float(make_float(16093626, 69), "95", 26); + check_float(make_float( 9983778, 25), "335", 12); + check_float(make_float(12745034, 104), "2585", 35); + check_float(make_float(12706553, 72), "60005", 24); + check_float(make_float(11005028, 45), "387205", 15); + check_float(make_float(15059547, 71), "3555835", 22); + check_float(make_float(16015691, -99), "25268305", -30); + check_float(make_float( 8667859, 56), "6245851", 17); + check_float(make_float(14855922, -82), "30721327", -25); + check_float(make_float(14855922, -83), "15360663", -25); + check_float(make_float(10144164, -110), "781478", -32); + check_float(make_float(13248074, 95), "52481028", 28); + } + + SECTION("double precision") + { + auto check_double = [](double number, const std::string& digits, int expected_exponent) + { + CAPTURE(number); + CAPTURE(digits); + CAPTURE(expected_exponent); + + char buf[32]; + int len = 0; + int exponent = 0; + nlohmann::detail::dtoa::grisu2(buf, len, exponent, number); + + CHECK(digits == std::string(buf, buf + len)); + CHECK(expected_exponent == exponent); + }; + + check_double(make_double(0, 0, 0x0000000000000001), "5", -324); // min denormal + check_double(make_double(0, 0, 0x000FFFFFFFFFFFFF), "2225073858507201", -323); // max denormal + check_double(make_double(0, 1, 0x0000000000000000), "22250738585072014", -324); // min normal + check_double(make_double(0, 1, 0x0000000000000001), "2225073858507202", -323); + check_double(make_double(0, 1, 0x000FFFFFFFFFFFFF), "44501477170144023", -324); + check_double(make_double(0, 2, 0x0000000000000000), "4450147717014403", -323); + check_double(make_double(0, 2, 0x0000000000000001), "4450147717014404", -323); + check_double(make_double(0, 4, 0x0000000000000000), "17800590868057611", -323); // fail if no special case in normalized boundaries + check_double(make_double(0, 5, 0x0000000000000000), "35601181736115222", -323); // fail if no special case in normalized boundaries + check_double(make_double(0, 6, 0x0000000000000000), "7120236347223045", -322); // fail if no special case in normalized boundaries + check_double(make_double(0, 10, 0x0000000000000000), "11392378155556871", -321); // fail if no special case in normalized boundaries + check_double(make_double(0, 2046, 0x000FFFFFFFFFFFFE), "17976931348623155", 292); + check_double(make_double(0, 2046, 0x000FFFFFFFFFFFFF), "17976931348623157", 292); // max normal + + // Test different paths in DigitGen + check_double( 10000, "1", 4); + check_double( 1200000, "12", 5); + check_double(4.9406564584124654e-324, "5", -324); // exit integral loop + check_double(2.2250738585072009e-308, "2225073858507201", -323); // exit fractional loop + check_double( 1.82877982605164e-99, "182877982605164", -113); + check_double( 1.1505466208671903e-09, "11505466208671903", -25); + check_double( 5.5645893133766722e+20, "5564589313376672", 5); + check_double( 53.034830388866226, "53034830388866226", -15); + check_double( 0.0021066531670178605, "21066531670178605", -19); + + // V. Paxson and W. Kahan, "A Program for Testing IEEE Binary-Decimal Conversion", manuscript, May 1991, + // ftp://ftp.ee.lbl.gov/testbase-report.ps.Z (report) + // ftp://ftp.ee.lbl.gov/testbase.tar.Z (program) + + // Table 3: Stress Inputs for Converting 53-bit Binary to Decimal, < 1/2 ULP + check_double(make_double(8511030020275656, -342) /* 9.5e-088 */, "95", -89); + check_double(make_double(5201988407066741, -824) /* 4.65e-233 */, "465", -235); + check_double(make_double(6406892948269899, +237) /* 1.415e+087 */, "1415", 84); + check_double(make_double(8431154198732492, +72) /* 3.9815e+037 */, "39815", 33); + check_double(make_double(6475049196144587, +99) /* 4.10405e+045 */, "410405", 40); + check_double(make_double(8274307542972842, +726) /* 2.920845e+234 */, "2920845", 228); + check_double(make_double(5381065484265332, -456) /* 2.8919465e-122 */, "28919465", -129); + check_double(make_double(6761728585499734, -1057) /* 4.37877185e-303 */, "437877185", -311); + check_double(make_double(7976538478610756, +376) /* 1.227701635e+129 */, "1227701635", 120); + check_double(make_double(5982403858958067, +377) /* 1.8415524525e+129 */, "18415524525", 119); + check_double(make_double(5536995190630837, +93) /* 5.48357443505e+043 */, "548357443505", 32); + check_double(make_double(7225450889282194, +710) /* 3.891901811465e+229 */, "3891901811465", 217); + check_double(make_double(7225450889282194, +709) /* 1.9459509057325e+229 */, "19459509057325", 216); + check_double(make_double(8703372741147379, +117) /* 1.44609583816055e+051 */, "144609583816055", 37); + check_double(make_double(8944262675275217, -1001) /* 4.173677474585315e-286 */, "4173677474585315", -301); + check_double(make_double(7459803696087692, -707) /* 1.1079507728788885e-197 */, "11079507728788885", -213); + check_double(make_double(6080469016670379, -381) /* 1.234550136632744e-099 */, "1234550136632744", -114); + check_double(make_double(8385515147034757, +721) /* 9.2503171196036502e+232 */, "925031711960365", 218); + check_double(make_double(7514216811389786, -828) /* 4.1980471502848898e-234 */, "419804715028489", -248); + check_double(make_double(8397297803260511, -345) /* 1.1716315319786511e-088 */, "11716315319786511", -104); + check_double(make_double(6733459239310543, +202) /* 4.3281007284461249e+076 */, "4328100728446125", 61); + check_double(make_double(8091450587292794, -473) /* 3.3177101181600311e-127 */, "3317710118160031", -142); + + // Table 4: Stress Inputs for Converting 53-bit Binary to Decimal, > 1/2 ULP + check_double(make_double(6567258882077402, +952) /* 2.5e+302 */, "25", 301); + check_double(make_double(6712731423444934, +535) /* 7.55e+176 */, "755", 174); + check_double(make_double(6712731423444934, +534) /* 3.775e+176 */, "3775", 173); + check_double(make_double(5298405411573037, -957) /* 4.3495e-273 */, "43495", -277); + check_double(make_double(5137311167659507, -144) /* 2.30365e-028 */, "230365", -33); + check_double(make_double(6722280709661868, +363) /* 1.263005e+125 */, "1263005", 119); + check_double(make_double(5344436398034927, -169) /* 7.1422105e-036 */, "71422105", -43); + check_double(make_double(8369123604277281, -853) /* 1.39345735e-241 */, "139345735", -249); + check_double(make_double(8995822108487663, -780) /* 1.414634485e-219 */, "1414634485", -228); + check_double(make_double(8942832835564782, -383) /* 4.5392779195e-100 */, "45392779195", -110); + check_double(make_double(8942832835564782, -384) /* 2.26963895975e-100 */, "226963895975", -111); + check_double(make_double(8942832835564782, -385) /* 1.134819479875e-100 */, "1134819479875", -112); + check_double(make_double(6965949469487146, -249) /* 7.7003665618895e-060 */, "77003665618895", -73); + check_double(make_double(6965949469487146, -250) /* 3.85018328094475e-060 */, "385018328094475", -74); + check_double(make_double(6965949469487146, -251) /* 1.925091640472375e-060 */, "1925091640472375", -75); + check_double(make_double(7487252720986826, +548) /* 6.8985865317742005e+180 */, "68985865317742005", 164); + check_double(make_double(5592117679628511, +164) /* 1.3076622631878654e+065 */, "13076622631878654", 49); + check_double(make_double(8887055249355788, +665) /* 1.3605202075612124e+216 */, "13605202075612124", 200); + check_double(make_double(6994187472632449, +690) /* 3.5928102174759597e+223 */, "35928102174759597", 207); + check_double(make_double(8797576579012143, +588) /* 8.9125197712484552e+192 */, "8912519771248455", 177); + check_double(make_double(7363326733505337, +272) /* 5.5876975736230114e+097 */, "55876975736230114", 81); + check_double(make_double(8549497411294502, -448) /* 1.1762578307285404e-119 */, "11762578307285404", -135); + + // Table 20: Stress Inputs for Converting 56-bit Binary to Decimal, < 1/2 ULP + check_double(make_double(50883641005312716, -172) /* 8.4999999999999993e-036 */, "8499999999999999", -51); + check_double(make_double(38162730753984537, -170) /* 2.5499999999999999e-035 */, "255", -37); + check_double(make_double(50832789069151999, -101) /* 2.0049999999999997e-014 */, "20049999999999997", -30); + check_double(make_double(51822367833714164, -109) /* 7.9844999999999994e-017 */, "7984499999999999", -32); + check_double(make_double(66840152193508133, -172) /* 1.1165499999999999e-035 */, "11165499999999999", -51); + check_double(make_double(55111239245584393, -138) /* 1.581615e-025 */, "1581615", -31); + check_double(make_double(71704866733321482, -112) /* 1.3809855e-017 */, "13809855", -24); + check_double(make_double(67160949328233173, -142) /* 1.2046404499999999e-026 */, "12046404499999999", -42); + check_double(make_double(53237141308040189, -152) /* 9.3251405449999991e-030 */, "9325140544999999", -45); + check_double(make_double(62785329394975786, -112) /* 1.2092014595e-017 */, "12092014595", -27); + check_double(make_double(48367680154689523, -77) /* 3.2007045838499998e-007 */, "320070458385", -18); + check_double(make_double(42552223180606797, -102) /* 8.391946324354999e-015 */, "8391946324354999", -30); + check_double(make_double(63626356173011241, -112) /* 1.2253990460585e-017 */, "12253990460585", -30); + check_double(make_double(43566388595783643, -99) /* 6.8735641489760495e-014 */, "687356414897605", -28); + check_double(make_double(54512669636675272, -159) /* 7.459816430480385e-032 */, "7459816430480385", -47); + check_double(make_double(52306490527514614, -167) /* 2.7960588398142552e-034 */, "2796058839814255", -49); + check_double(make_double(52306490527514614, -168) /* 1.3980294199071276e-034 */, "13980294199071276", -50); + check_double(make_double(41024721590449423, -89) /* 6.6279012373057359e-011 */, "6627901237305736", -26); + check_double(make_double(37664020415894738, -132) /* 6.9177880043968072e-024 */, "6917788004396807", -39); + check_double(make_double(37549883692866294, -93) /* 3.7915693108349708e-012 */, "3791569310834971", -27); + check_double(make_double(69124110374399839, -104) /* 3.4080817676591365e-015 */, "34080817676591365", -31); + check_double(make_double(69124110374399839, -105) /* 1.7040408838295683e-015 */, "17040408838295683", -31); + + // Table 21: Stress Inputs for Converting 56-bit Binary to Decimal, > 1/2 ULP + check_double(make_double(49517601571415211, -94) /* 2.4999999999999998e-012 */, "25", -13); + check_double(make_double(49517601571415211, -95) /* 1.2499999999999999e-012 */, "125", -14); + check_double(make_double(54390733528642804, -133) /* 4.9949999999999996e-024 */, "49949999999999996", -40); // shortest: 4995e-27 + check_double(make_double(71805402319113924, -157) /* 3.9304999999999998e-031 */, "39304999999999998", -47); // shortest: 39305e-35 + check_double(make_double(40435277969631694, -179) /* 5.2770499999999992e-038 */, "5277049999999999", -53); + check_double(make_double(57241991568619049, -165) /* 1.223955e-033 */, "1223955", -39); + check_double(make_double(65224162876242886, +58) /* 1.8799584999999998e+034 */, "18799584999999998", 18); + check_double(make_double(70173376848895368, -138) /* 2.01387715e-025 */, "201387715", -33); + check_double(make_double(37072848117383207, -99) /* 5.8490641049999989e-014 */, "5849064104999999", -29); + check_double(make_double(56845051585389697, -176) /* 5.9349003054999999e-037 */, "59349003055", -47); + check_double(make_double(54791673366936431, -145) /* 1.2284718039499998e-027 */, "12284718039499998", -43); + check_double(make_double(66800318669106231, -169) /* 8.9270767180849991e-035 */, "8927076718084999", -50); + check_double(make_double(66800318669106231, -170) /* 4.4635383590424995e-035 */, "44635383590424995", -51); + check_double(make_double(66574323440112438, -119) /* 1.0016990862549499e-019 */, "10016990862549499", -35); + check_double(make_double(65645179969330963, -173) /* 5.4829412628024647e-036 */, "5482941262802465", -51); + check_double(make_double(61847254334681076, -109) /* 9.5290783281036439e-017 */, "9529078328103644", -32); + check_double(make_double(39990712921393606, -145) /* 8.9662279366405553e-028 */, "8966227936640555", -43); + check_double(make_double(59292318184400283, -149) /* 8.3086234418058538e-029 */, "8308623441805854", -44); + check_double(make_double(69116558615326153, -143) /* 6.1985873566126555e-027 */, "61985873566126555", -43); + check_double(make_double(69116558615326153, -144) /* 3.0992936783063277e-027 */, "30992936783063277", -43); + check_double(make_double(39462549494468513, -152) /* 6.9123512506176015e-030 */, "6912351250617602", -45); + check_double(make_double(39462549494468513, -153) /* 3.4561756253088008e-030 */, "3456175625308801", -45); + } +} + +TEST_CASE("formatting") +{ + SECTION("single precision") + { + auto check_float = [](float number, const std::string& expected) + { + char buf[32]; + char* end = nlohmann::detail::dtoa::to_short_string(buf, buf + 32, number); + std::string actual(buf, end); + + CHECK(actual == expected); + }; + // %.9g + check_float( -1.2345e-22f, "-1.2345e-22" ); // -1.23450004e-22 + check_float( -1.2345e-21f, "-1.2345e-21" ); // -1.23450002e-21 + check_float( -1.2345e-20f, "-1.2345e-20" ); // -1.23450002e-20 + check_float( -1.2345e-19f, "-1.2345e-19" ); // -1.23449999e-19 + check_float( -1.2345e-18f, "-1.2345e-18" ); // -1.23449996e-18 + check_float( -1.2345e-17f, "-1.2345e-17" ); // -1.23449998e-17 + check_float( -1.2345e-16f, "-1.2345e-16" ); // -1.23449996e-16 + check_float( -1.2345e-15f, "-1.2345e-15" ); // -1.23450002e-15 + check_float( -1.2345e-14f, "-1.2345e-14" ); // -1.23450004e-14 + check_float( -1.2345e-13f, "-1.2345e-13" ); // -1.23449997e-13 + check_float( -1.2345e-12f, "-1.2345e-12" ); // -1.23450002e-12 + check_float( -1.2345e-11f, "-1.2345e-11" ); // -1.2345e-11 + check_float( -1.2345e-10f, "-1.2345e-10" ); // -1.2345e-10 + check_float( -1.2345e-9f, "-1.2345e-9" ); // -1.23449995e-09 + check_float( -1.2345e-8f, "-1.2345e-8" ); // -1.23449997e-08 + check_float( -1.2345e-7f, "-1.2345e-7" ); // -1.23449993e-07 + check_float( -1.2345e-6f, "-1.2345e-6" ); // -1.23450002e-06 + check_float( -1.2345e-5f, "-1.2345e-5" ); // -1.2345e-05 + check_float( -1.2345e-4f, "-0.00012345" ); // -0.000123449994 + check_float( -1.2345e-3f, "-0.0012345" ); // -0.00123449997 + check_float( -1.2345e-2f, "-0.012345" ); // -0.0123450002 + check_float( -1.2345e-1f, "-0.12345" ); // -0.123450004 + check_float( -0.0f, "-0.0" ); // -0 + check_float( 0.0f, "0.0" ); // 0 + check_float( 1.2345e+0f, "1.2345" ); // 1.23450005 + check_float( 1.2345e+1f, "12.345" ); // 12.3450003 + check_float( 1.2345e+2f, "123.45" ); // 123.449997 + check_float( 1.2345e+3f, "1234.5" ); // 1234.5 + check_float( 1.2345e+4f, "12345.0" ); // 12345 + check_float( 1.2345e+5f, "123450.0" ); // 123450 + check_float( 1.2345e+6f, "1234500.0" ); // 1234500 + check_float( 1.2345e+7f, "12345000.0" ); // 12345000 + check_float( 1.2345e+8f, "123450000.0" ); // 123450000 + check_float( 1.2345e+9f, "1.2345e+9" ); // 1.23449997e+09 + check_float( 1.2345e+10f, "1.2345e+10" ); // 1.23449999e+10 + check_float( 1.2345e+11f, "1.2345e+11" ); // 1.23449999e+11 + check_float( 1.2345e+12f, "1.2345e+12" ); // 1.23450006e+12 + check_float( 1.2345e+13f, "1.2345e+13" ); // 1.23449995e+13 + check_float( 1.2345e+14f, "1.2345e+14" ); // 1.23450002e+14 + check_float( 1.2345e+15f, "1.2345e+15" ); // 1.23450003e+15 + check_float( 1.2345e+16f, "1.2345e+16" ); // 1.23449998e+16 + check_float( 1.2345e+17f, "1.2345e+17" ); // 1.23449996e+17 + check_float( 1.2345e+18f, "1.2345e+18" ); // 1.23450004e+18 + check_float( 1.2345e+19f, "1.2345e+19" ); // 1.23449999e+19 + check_float( 1.2345e+20f, "1.2345e+20" ); // 1.23449999e+20 + check_float( 1.2345e+21f, "1.2345e+21" ); // 1.23449999e+21 + check_float( 1.2345e+22f, "1.2345e+22" ); // 1.23450005e+22 + } + + SECTION("double precision") + { + auto check_double = [](double number, const std::string& expected) + { + char buf[32]; + char* end = nlohmann::detail::dtoa::to_short_string(buf, buf + 32, number); + std::string actual(buf, end); + + CHECK(actual == expected); + }; + // %.17g + check_double( -1.2345e-22, "-1.2345e-22" ); // -1.2345000000000001e-22 + check_double( -1.2345e-21, "-1.2345e-21" ); // -1.2345000000000001e-21 + check_double( -1.2345e-20, "-1.2345e-20" ); // -1.2345e-20 + check_double( -1.2345e-19, "-1.2345e-19" ); // -1.2345000000000001e-19 + check_double( -1.2345e-18, "-1.2345e-18" ); // -1.2345000000000001e-18 + check_double( -1.2345e-17, "-1.2345e-17" ); // -1.2345e-17 + check_double( -1.2345e-16, "-1.2345e-16" ); // -1.2344999999999999e-16 + check_double( -1.2345e-15, "-1.2345e-15" ); // -1.2345e-15 + check_double( -1.2345e-14, "-1.2345e-14" ); // -1.2345e-14 + check_double( -1.2345e-13, "-1.2345e-13" ); // -1.2344999999999999e-13 + check_double( -1.2345e-12, "-1.2345e-12" ); // -1.2345e-12 + check_double( -1.2345e-11, "-1.2345e-11" ); // -1.2345e-11 + check_double( -1.2345e-10, "-1.2345e-10" ); // -1.2345e-10 + check_double( -1.2345e-9, "-1.2345e-9" ); // -1.2345e-09 + check_double( -1.2345e-8, "-1.2345e-8" ); // -1.2345000000000001e-08 + check_double( -1.2345e-7, "-1.2345e-7" ); // -1.2345000000000001e-07 + check_double( -1.2345e-6, "-1.2345e-6" ); // -1.2345e-06 + check_double( -1.2345e-5, "-1.2345e-5" ); // -1.2345e-05 + check_double( -1.2345e-4, "-0.00012345" ); // -0.00012344999999999999 + check_double( -1.2345e-3, "-0.0012345" ); // -0.0012344999999999999 + check_double( -1.2345e-2, "-0.012345" ); // -0.012345 + check_double( -1.2345e-1, "-0.12345" ); // -0.12345 + check_double( -0.0, "-0.0" ); // -0 + check_double( 0.0, "0.0" ); // 0 + check_double( 1.2345e+0, "1.2345" ); // 1.2344999999999999 + check_double( 1.2345e+1, "12.345" ); // 12.345000000000001 + check_double( 1.2345e+2, "123.45" ); // 123.45 + check_double( 1.2345e+3, "1234.5" ); // 1234.5 + check_double( 1.2345e+4, "12345.0" ); // 12345 + check_double( 1.2345e+5, "123450.0" ); // 123450 + check_double( 1.2345e+6, "1234500.0" ); // 1234500 + check_double( 1.2345e+7, "12345000.0" ); // 12345000 + check_double( 1.2345e+8, "123450000.0" ); // 123450000 + check_double( 1.2345e+9, "1234500000.0" ); // 1234500000 + check_double( 1.2345e+10, "12345000000.0" ); // 12345000000 + check_double( 1.2345e+11, "123450000000.0" ); // 123450000000 + check_double( 1.2345e+12, "1234500000000.0" ); // 1234500000000 + check_double( 1.2345e+13, "12345000000000.0" ); // 12345000000000 + check_double( 1.2345e+14, "123450000000000.0" ); // 123450000000000 + check_double( 1.2345e+15, "1234500000000000.0" ); // 1234500000000000 + check_double( 1.2345e+16, "12345000000000000.0" ); // 12345000000000000 + check_double( 1.2345e+17, "1.2345e+17" ); // 1.2345e+17 + check_double( 1.2345e+18, "1.2345e+18" ); // 1.2345e+18 + check_double( 1.2345e+19, "1.2345e+19" ); // 1.2345e+19 + check_double( 1.2345e+20, "1.2345e+20" ); // 1.2345e+20 + check_double( 1.2345e+21, "1.2344999999999999e+21" ); // 1.2344999999999999e+21 (shortest: 1.2345e+21) + check_double( 1.2345e+22, "1.2345e+22" ); // 1.2345e+22 + } +}