Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 10:02:17

0001 //     __ _____ _____ _____
0002 //  __|  |   __|     |   | |  JSON for Modern C++
0003 // |  |  |__   |  |  | | | |  version 3.11.2
0004 // |_____|_____|_____|_|___|  https://github.com/nlohmann/json
0005 //
0006 // SPDX-FileCopyrightText: 2009 Florian Loitsch <https://florian.loitsch.com/>
0007 // SPDX-FileCopyrightText: 2013-2022 Niels Lohmann <https://nlohmann.me>
0008 // SPDX-License-Identifier: MIT
0009 
0010 #pragma once
0011 
0012 #include <array> // array
0013 #include <cmath>   // signbit, isfinite
0014 #include <cstdint> // intN_t, uintN_t
0015 #include <cstring> // memcpy, memmove
0016 #include <limits> // numeric_limits
0017 #include <type_traits> // conditional
0018 
0019 #include <nlohmann/detail/macro_scope.hpp>
0020 
0021 NLOHMANN_JSON_NAMESPACE_BEGIN
0022 namespace detail
0023 {
0024 
0025 /*!
0026 @brief implements the Grisu2 algorithm for binary to decimal floating-point
0027 conversion.
0028 
0029 This implementation is a slightly modified version of the reference
0030 implementation which may be obtained from
0031 http://florian.loitsch.com/publications (bench.tar.gz).
0032 
0033 The code is distributed under the MIT license, Copyright (c) 2009 Florian Loitsch.
0034 
0035 For a detailed description of the algorithm see:
0036 
0037 [1] Loitsch, "Printing Floating-Point Numbers Quickly and Accurately with
0038     Integers", Proceedings of the ACM SIGPLAN 2010 Conference on Programming
0039     Language Design and Implementation, PLDI 2010
0040 [2] Burger, Dybvig, "Printing Floating-Point Numbers Quickly and Accurately",
0041     Proceedings of the ACM SIGPLAN 1996 Conference on Programming Language
0042     Design and Implementation, PLDI 1996
0043 */
0044 namespace dtoa_impl
0045 {
0046 
0047 template<typename Target, typename Source>
0048 Target reinterpret_bits(const Source source)
0049 {
0050     static_assert(sizeof(Target) == sizeof(Source), "size mismatch");
0051 
0052     Target target;
0053     std::memcpy(&target, &source, sizeof(Source));
0054     return target;
0055 }
0056 
0057 struct diyfp // f * 2^e
0058 {
0059     static constexpr int kPrecision = 64; // = q
0060 
0061     std::uint64_t f = 0;
0062     int e = 0;
0063 
0064     constexpr diyfp(std::uint64_t f_, int e_) noexcept : f(f_), e(e_) {}
0065 
0066     /*!
0067     @brief returns x - y
0068     @pre x.e == y.e and x.f >= y.f
0069     */
0070     static diyfp sub(const diyfp& x, const diyfp& y) noexcept
0071     {
0072         JSON_ASSERT(x.e == y.e);
0073         JSON_ASSERT(x.f >= y.f);
0074 
0075         return {x.f - y.f, x.e};
0076     }
0077 
0078     /*!
0079     @brief returns x * y
0080     @note The result is rounded. (Only the upper q bits are returned.)
0081     */
0082     static diyfp mul(const diyfp& x, const diyfp& y) noexcept
0083     {
0084         static_assert(kPrecision == 64, "internal error");
0085 
0086         // Computes:
0087         //  f = round((x.f * y.f) / 2^q)
0088         //  e = x.e + y.e + q
0089 
0090         // Emulate the 64-bit * 64-bit multiplication:
0091         //
0092         // p = u * v
0093         //   = (u_lo + 2^32 u_hi) (v_lo + 2^32 v_hi)
0094         //   = (u_lo v_lo         ) + 2^32 ((u_lo v_hi         ) + (u_hi v_lo         )) + 2^64 (u_hi v_hi         )
0095         //   = (p0                ) + 2^32 ((p1                ) + (p2                )) + 2^64 (p3                )
0096         //   = (p0_lo + 2^32 p0_hi) + 2^32 ((p1_lo + 2^32 p1_hi) + (p2_lo + 2^32 p2_hi)) + 2^64 (p3                )
0097         //   = (p0_lo             ) + 2^32 (p0_hi + p1_lo + p2_lo                      ) + 2^64 (p1_hi + p2_hi + p3)
0098         //   = (p0_lo             ) + 2^32 (Q                                          ) + 2^64 (H                 )
0099         //   = (p0_lo             ) + 2^32 (Q_lo + 2^32 Q_hi                           ) + 2^64 (H                 )
0100         //
0101         // (Since Q might be larger than 2^32 - 1)
0102         //
0103         //   = (p0_lo + 2^32 Q_lo) + 2^64 (Q_hi + H)
0104         //
0105         // (Q_hi + H does not overflow a 64-bit int)
0106         //
0107         //   = p_lo + 2^64 p_hi
0108 
0109         const std::uint64_t u_lo = x.f & 0xFFFFFFFFu;
0110         const std::uint64_t u_hi = x.f >> 32u;
0111         const std::uint64_t v_lo = y.f & 0xFFFFFFFFu;
0112         const std::uint64_t v_hi = y.f >> 32u;
0113 
0114         const std::uint64_t p0 = u_lo * v_lo;
0115         const std::uint64_t p1 = u_lo * v_hi;
0116         const std::uint64_t p2 = u_hi * v_lo;
0117         const std::uint64_t p3 = u_hi * v_hi;
0118 
0119         const std::uint64_t p0_hi = p0 >> 32u;
0120         const std::uint64_t p1_lo = p1 & 0xFFFFFFFFu;
0121         const std::uint64_t p1_hi = p1 >> 32u;
0122         const std::uint64_t p2_lo = p2 & 0xFFFFFFFFu;
0123         const std::uint64_t p2_hi = p2 >> 32u;
0124 
0125         std::uint64_t Q = p0_hi + p1_lo + p2_lo;
0126 
0127         // The full product might now be computed as
0128         //
0129         // p_hi = p3 + p2_hi + p1_hi + (Q >> 32)
0130         // p_lo = p0_lo + (Q << 32)
0131         //
0132         // But in this particular case here, the full p_lo is not required.
0133         // Effectively we only need to add the highest bit in p_lo to p_hi (and
0134         // Q_hi + 1 does not overflow).
0135 
0136         Q += std::uint64_t{1} << (64u - 32u - 1u); // round, ties up
0137 
0138         const std::uint64_t h = p3 + p2_hi + p1_hi + (Q >> 32u);
0139 
0140         return {h, x.e + y.e + 64};
0141     }
0142 
0143     /*!
0144     @brief normalize x such that the significand is >= 2^(q-1)
0145     @pre x.f != 0
0146     */
0147     static diyfp normalize(diyfp x) noexcept
0148     {
0149         JSON_ASSERT(x.f != 0);
0150 
0151         while ((x.f >> 63u) == 0)
0152         {
0153             x.f <<= 1u;
0154             x.e--;
0155         }
0156 
0157         return x;
0158     }
0159 
0160     /*!
0161     @brief normalize x such that the result has the exponent E
0162     @pre e >= x.e and the upper e - x.e bits of x.f must be zero.
0163     */
0164     static diyfp normalize_to(const diyfp& x, const int target_exponent) noexcept
0165     {
0166         const int delta = x.e - target_exponent;
0167 
0168         JSON_ASSERT(delta >= 0);
0169         JSON_ASSERT(((x.f << delta) >> delta) == x.f);
0170 
0171         return {x.f << delta, target_exponent};
0172     }
0173 };
0174 
0175 struct boundaries
0176 {
0177     diyfp w;
0178     diyfp minus;
0179     diyfp plus;
0180 };
0181 
0182 /*!
0183 Compute the (normalized) diyfp representing the input number 'value' and its
0184 boundaries.
0185 
0186 @pre value must be finite and positive
0187 */
0188 template<typename FloatType>
0189 boundaries compute_boundaries(FloatType value)
0190 {
0191     JSON_ASSERT(std::isfinite(value));
0192     JSON_ASSERT(value > 0);
0193 
0194     // Convert the IEEE representation into a diyfp.
0195     //
0196     // If v is denormal:
0197     //      value = 0.F * 2^(1 - bias) = (          F) * 2^(1 - bias - (p-1))
0198     // If v is normalized:
0199     //      value = 1.F * 2^(E - bias) = (2^(p-1) + F) * 2^(E - bias - (p-1))
0200 
0201     static_assert(std::numeric_limits<FloatType>::is_iec559,
0202                   "internal error: dtoa_short requires an IEEE-754 floating-point implementation");
0203 
0204     constexpr int      kPrecision = std::numeric_limits<FloatType>::digits; // = p (includes the hidden bit)
0205     constexpr int      kBias      = std::numeric_limits<FloatType>::max_exponent - 1 + (kPrecision - 1);
0206     constexpr int      kMinExp    = 1 - kBias;
0207     constexpr std::uint64_t kHiddenBit = std::uint64_t{1} << (kPrecision - 1); // = 2^(p-1)
0208 
0209     using bits_type = typename std::conditional<kPrecision == 24, std::uint32_t, std::uint64_t >::type;
0210 
0211     const auto bits = static_cast<std::uint64_t>(reinterpret_bits<bits_type>(value));
0212     const std::uint64_t E = bits >> (kPrecision - 1);
0213     const std::uint64_t F = bits & (kHiddenBit - 1);
0214 
0215     const bool is_denormal = E == 0;
0216     const diyfp v = is_denormal
0217                     ? diyfp(F, kMinExp)
0218                     : diyfp(F + kHiddenBit, static_cast<int>(E) - kBias);
0219 
0220     // Compute the boundaries m- and m+ of the floating-point value
0221     // v = f * 2^e.
0222     //
0223     // Determine v- and v+, the floating-point predecessor and successor if v,
0224     // respectively.
0225     //
0226     //      v- = v - 2^e        if f != 2^(p-1) or e == e_min                (A)
0227     //         = v - 2^(e-1)    if f == 2^(p-1) and e > e_min                (B)
0228     //
0229     //      v+ = v + 2^e
0230     //
0231     // Let m- = (v- + v) / 2 and m+ = (v + v+) / 2. All real numbers _strictly_
0232     // between m- and m+ round to v, regardless of how the input rounding
0233     // algorithm breaks ties.
0234     //
0235     //      ---+-------------+-------------+-------------+-------------+---  (A)
0236     //         v-            m-            v             m+            v+
0237     //
0238     //      -----------------+------+------+-------------+-------------+---  (B)
0239     //                       v-     m-     v             m+            v+
0240 
0241     const bool lower_boundary_is_closer = F == 0 && E > 1;
0242     const diyfp m_plus = diyfp(2 * v.f + 1, v.e - 1);
0243     const diyfp m_minus = lower_boundary_is_closer
0244                           ? diyfp(4 * v.f - 1, v.e - 2)  // (B)
0245                           : diyfp(2 * v.f - 1, v.e - 1); // (A)
0246 
0247     // Determine the normalized w+ = m+.
0248     const diyfp w_plus = diyfp::normalize(m_plus);
0249 
0250     // Determine w- = m- such that e_(w-) = e_(w+).
0251     const diyfp w_minus = diyfp::normalize_to(m_minus, w_plus.e);
0252 
0253     return {diyfp::normalize(v), w_minus, w_plus};
0254 }
0255 
0256 // Given normalized diyfp w, Grisu needs to find a (normalized) cached
0257 // power-of-ten c, such that the exponent of the product c * w = f * 2^e lies
0258 // within a certain range [alpha, gamma] (Definition 3.2 from [1])
0259 //
0260 //      alpha <= e = e_c + e_w + q <= gamma
0261 //
0262 // or
0263 //
0264 //      f_c * f_w * 2^alpha <= f_c 2^(e_c) * f_w 2^(e_w) * 2^q
0265 //                          <= f_c * f_w * 2^gamma
0266 //
0267 // Since c and w are normalized, i.e. 2^(q-1) <= f < 2^q, this implies
0268 //
0269 //      2^(q-1) * 2^(q-1) * 2^alpha <= c * w * 2^q < 2^q * 2^q * 2^gamma
0270 //
0271 // or
0272 //
0273 //      2^(q - 2 + alpha) <= c * w < 2^(q + gamma)
0274 //
0275 // The choice of (alpha,gamma) determines the size of the table and the form of
0276 // the digit generation procedure. Using (alpha,gamma)=(-60,-32) works out well
0277 // in practice:
0278 //
0279 // The idea is to cut the number c * w = f * 2^e into two parts, which can be
0280 // processed independently: An integral part p1, and a fractional part p2:
0281 //
0282 //      f * 2^e = ( (f div 2^-e) * 2^-e + (f mod 2^-e) ) * 2^e
0283 //              = (f div 2^-e) + (f mod 2^-e) * 2^e
0284 //              = p1 + p2 * 2^e
0285 //
0286 // The conversion of p1 into decimal form requires a series of divisions and
0287 // modulos by (a power of) 10. These operations are faster for 32-bit than for
0288 // 64-bit integers, so p1 should ideally fit into a 32-bit integer. This can be
0289 // achieved by choosing
0290 //
0291 //      -e >= 32   or   e <= -32 := gamma
0292 //
0293 // In order to convert the fractional part
0294 //
0295 //      p2 * 2^e = p2 / 2^-e = d[-1] / 10^1 + d[-2] / 10^2 + ...
0296 //
0297 // into decimal form, the fraction is repeatedly multiplied by 10 and the digits
0298 // d[-i] are extracted in order:
0299 //
0300 //      (10 * p2) div 2^-e = d[-1]
0301 //      (10 * p2) mod 2^-e = d[-2] / 10^1 + ...
0302 //
0303 // The multiplication by 10 must not overflow. It is sufficient to choose
0304 //
0305 //      10 * p2 < 16 * p2 = 2^4 * p2 <= 2^64.
0306 //
0307 // Since p2 = f mod 2^-e < 2^-e,
0308 //
0309 //      -e <= 60   or   e >= -60 := alpha
0310 
0311 constexpr int kAlpha = -60;
0312 constexpr int kGamma = -32;
0313 
0314 struct cached_power // c = f * 2^e ~= 10^k
0315 {
0316     std::uint64_t f;
0317     int e;
0318     int k;
0319 };
0320 
0321 /*!
0322 For a normalized diyfp w = f * 2^e, this function returns a (normalized) cached
0323 power-of-ten c = f_c * 2^e_c, such that the exponent of the product w * c
0324 satisfies (Definition 3.2 from [1])
0325 
0326      alpha <= e_c + e + q <= gamma.
0327 */
0328 inline cached_power get_cached_power_for_binary_exponent(int e)
0329 {
0330     // Now
0331     //
0332     //      alpha <= e_c + e + q <= gamma                                    (1)
0333     //      ==> f_c * 2^alpha <= c * 2^e * 2^q
0334     //
0335     // and since the c's are normalized, 2^(q-1) <= f_c,
0336     //
0337     //      ==> 2^(q - 1 + alpha) <= c * 2^(e + q)
0338     //      ==> 2^(alpha - e - 1) <= c
0339     //
0340     // If c were an exact power of ten, i.e. c = 10^k, one may determine k as
0341     //
0342     //      k = ceil( log_10( 2^(alpha - e - 1) ) )
0343     //        = ceil( (alpha - e - 1) * log_10(2) )
0344     //
0345     // From the paper:
0346     // "In theory the result of the procedure could be wrong since c is rounded,
0347     //  and the computation itself is approximated [...]. In practice, however,
0348     //  this simple function is sufficient."
0349     //
0350     // For IEEE double precision floating-point numbers converted into
0351     // normalized diyfp's w = f * 2^e, with q = 64,
0352     //
0353     //      e >= -1022      (min IEEE exponent)
0354     //           -52        (p - 1)
0355     //           -52        (p - 1, possibly normalize denormal IEEE numbers)
0356     //           -11        (normalize the diyfp)
0357     //         = -1137
0358     //
0359     // and
0360     //
0361     //      e <= +1023      (max IEEE exponent)
0362     //           -52        (p - 1)
0363     //           -11        (normalize the diyfp)
0364     //         = 960
0365     //
0366     // This binary exponent range [-1137,960] results in a decimal exponent
0367     // range [-307,324]. One does not need to store a cached power for each
0368     // k in this range. For each such k it suffices to find a cached power
0369     // such that the exponent of the product lies in [alpha,gamma].
0370     // This implies that the difference of the decimal exponents of adjacent
0371     // table entries must be less than or equal to
0372     //
0373     //      floor( (gamma - alpha) * log_10(2) ) = 8.
0374     //
0375     // (A smaller distance gamma-alpha would require a larger table.)
0376 
0377     // NB:
0378     // Actually this function returns c, such that -60 <= e_c + e + 64 <= -34.
0379 
0380     constexpr int kCachedPowersMinDecExp = -300;
0381     constexpr int kCachedPowersDecStep = 8;
0382 
0383     static constexpr std::array<cached_power, 79> kCachedPowers =
0384     {
0385         {
0386             { 0xAB70FE17C79AC6CA, -1060, -300 },
0387             { 0xFF77B1FCBEBCDC4F, -1034, -292 },
0388             { 0xBE5691EF416BD60C, -1007, -284 },
0389             { 0x8DD01FAD907FFC3C,  -980, -276 },
0390             { 0xD3515C2831559A83,  -954, -268 },
0391             { 0x9D71AC8FADA6C9B5,  -927, -260 },
0392             { 0xEA9C227723EE8BCB,  -901, -252 },
0393             { 0xAECC49914078536D,  -874, -244 },
0394             { 0x823C12795DB6CE57,  -847, -236 },
0395             { 0xC21094364DFB5637,  -821, -228 },
0396             { 0x9096EA6F3848984F,  -794, -220 },
0397             { 0xD77485CB25823AC7,  -768, -212 },
0398             { 0xA086CFCD97BF97F4,  -741, -204 },
0399             { 0xEF340A98172AACE5,  -715, -196 },
0400             { 0xB23867FB2A35B28E,  -688, -188 },
0401             { 0x84C8D4DFD2C63F3B,  -661, -180 },
0402             { 0xC5DD44271AD3CDBA,  -635, -172 },
0403             { 0x936B9FCEBB25C996,  -608, -164 },
0404             { 0xDBAC6C247D62A584,  -582, -156 },
0405             { 0xA3AB66580D5FDAF6,  -555, -148 },
0406             { 0xF3E2F893DEC3F126,  -529, -140 },
0407             { 0xB5B5ADA8AAFF80B8,  -502, -132 },
0408             { 0x87625F056C7C4A8B,  -475, -124 },
0409             { 0xC9BCFF6034C13053,  -449, -116 },
0410             { 0x964E858C91BA2655,  -422, -108 },
0411             { 0xDFF9772470297EBD,  -396, -100 },
0412             { 0xA6DFBD9FB8E5B88F,  -369,  -92 },
0413             { 0xF8A95FCF88747D94,  -343,  -84 },
0414             { 0xB94470938FA89BCF,  -316,  -76 },
0415             { 0x8A08F0F8BF0F156B,  -289,  -68 },
0416             { 0xCDB02555653131B6,  -263,  -60 },
0417             { 0x993FE2C6D07B7FAC,  -236,  -52 },
0418             { 0xE45C10C42A2B3B06,  -210,  -44 },
0419             { 0xAA242499697392D3,  -183,  -36 },
0420             { 0xFD87B5F28300CA0E,  -157,  -28 },
0421             { 0xBCE5086492111AEB,  -130,  -20 },
0422             { 0x8CBCCC096F5088CC,  -103,  -12 },
0423             { 0xD1B71758E219652C,   -77,   -4 },
0424             { 0x9C40000000000000,   -50,    4 },
0425             { 0xE8D4A51000000000,   -24,   12 },
0426             { 0xAD78EBC5AC620000,     3,   20 },
0427             { 0x813F3978F8940984,    30,   28 },
0428             { 0xC097CE7BC90715B3,    56,   36 },
0429             { 0x8F7E32CE7BEA5C70,    83,   44 },
0430             { 0xD5D238A4ABE98068,   109,   52 },
0431             { 0x9F4F2726179A2245,   136,   60 },
0432             { 0xED63A231D4C4FB27,   162,   68 },
0433             { 0xB0DE65388CC8ADA8,   189,   76 },
0434             { 0x83C7088E1AAB65DB,   216,   84 },
0435             { 0xC45D1DF942711D9A,   242,   92 },
0436             { 0x924D692CA61BE758,   269,  100 },
0437             { 0xDA01EE641A708DEA,   295,  108 },
0438             { 0xA26DA3999AEF774A,   322,  116 },
0439             { 0xF209787BB47D6B85,   348,  124 },
0440             { 0xB454E4A179DD1877,   375,  132 },
0441             { 0x865B86925B9BC5C2,   402,  140 },
0442             { 0xC83553C5C8965D3D,   428,  148 },
0443             { 0x952AB45CFA97A0B3,   455,  156 },
0444             { 0xDE469FBD99A05FE3,   481,  164 },
0445             { 0xA59BC234DB398C25,   508,  172 },
0446             { 0xF6C69A72A3989F5C,   534,  180 },
0447             { 0xB7DCBF5354E9BECE,   561,  188 },
0448             { 0x88FCF317F22241E2,   588,  196 },
0449             { 0xCC20CE9BD35C78A5,   614,  204 },
0450             { 0x98165AF37B2153DF,   641,  212 },
0451             { 0xE2A0B5DC971F303A,   667,  220 },
0452             { 0xA8D9D1535CE3B396,   694,  228 },
0453             { 0xFB9B7CD9A4A7443C,   720,  236 },
0454             { 0xBB764C4CA7A44410,   747,  244 },
0455             { 0x8BAB8EEFB6409C1A,   774,  252 },
0456             { 0xD01FEF10A657842C,   800,  260 },
0457             { 0x9B10A4E5E9913129,   827,  268 },
0458             { 0xE7109BFBA19C0C9D,   853,  276 },
0459             { 0xAC2820D9623BF429,   880,  284 },
0460             { 0x80444B5E7AA7CF85,   907,  292 },
0461             { 0xBF21E44003ACDD2D,   933,  300 },
0462             { 0x8E679C2F5E44FF8F,   960,  308 },
0463             { 0xD433179D9C8CB841,   986,  316 },
0464             { 0x9E19DB92B4E31BA9,  1013,  324 },
0465         }
0466     };
0467 
0468     // This computation gives exactly the same results for k as
0469     //      k = ceil((kAlpha - e - 1) * 0.30102999566398114)
0470     // for |e| <= 1500, but doesn't require floating-point operations.
0471     // NB: log_10(2) ~= 78913 / 2^18
0472     JSON_ASSERT(e >= -1500);
0473     JSON_ASSERT(e <=  1500);
0474     const int f = kAlpha - e - 1;
0475     const int k = (f * 78913) / (1 << 18) + static_cast<int>(f > 0);
0476 
0477     const int index = (-kCachedPowersMinDecExp + k + (kCachedPowersDecStep - 1)) / kCachedPowersDecStep;
0478     JSON_ASSERT(index >= 0);
0479     JSON_ASSERT(static_cast<std::size_t>(index) < kCachedPowers.size());
0480 
0481     const cached_power cached = kCachedPowers[static_cast<std::size_t>(index)];
0482     JSON_ASSERT(kAlpha <= cached.e + e + 64);
0483     JSON_ASSERT(kGamma >= cached.e + e + 64);
0484 
0485     return cached;
0486 }
0487 
0488 /*!
0489 For n != 0, returns k, such that pow10 := 10^(k-1) <= n < 10^k.
0490 For n == 0, returns 1 and sets pow10 := 1.
0491 */
0492 inline int find_largest_pow10(const std::uint32_t n, std::uint32_t& pow10)
0493 {
0494     // LCOV_EXCL_START
0495     if (n >= 1000000000)
0496     {
0497         pow10 = 1000000000;
0498         return 10;
0499     }
0500     // LCOV_EXCL_STOP
0501     if (n >= 100000000)
0502     {
0503         pow10 = 100000000;
0504         return  9;
0505     }
0506     if (n >= 10000000)
0507     {
0508         pow10 = 10000000;
0509         return  8;
0510     }
0511     if (n >= 1000000)
0512     {
0513         pow10 = 1000000;
0514         return  7;
0515     }
0516     if (n >= 100000)
0517     {
0518         pow10 = 100000;
0519         return  6;
0520     }
0521     if (n >= 10000)
0522     {
0523         pow10 = 10000;
0524         return  5;
0525     }
0526     if (n >= 1000)
0527     {
0528         pow10 = 1000;
0529         return  4;
0530     }
0531     if (n >= 100)
0532     {
0533         pow10 = 100;
0534         return  3;
0535     }
0536     if (n >= 10)
0537     {
0538         pow10 = 10;
0539         return  2;
0540     }
0541 
0542     pow10 = 1;
0543     return 1;
0544 }
0545 
0546 inline void grisu2_round(char* buf, int len, std::uint64_t dist, std::uint64_t delta,
0547                          std::uint64_t rest, std::uint64_t ten_k)
0548 {
0549     JSON_ASSERT(len >= 1);
0550     JSON_ASSERT(dist <= delta);
0551     JSON_ASSERT(rest <= delta);
0552     JSON_ASSERT(ten_k > 0);
0553 
0554     //               <--------------------------- delta ---->
0555     //                                  <---- dist --------->
0556     // --------------[------------------+-------------------]--------------
0557     //               M-                 w                   M+
0558     //
0559     //                                  ten_k
0560     //                                <------>
0561     //                                       <---- rest ---->
0562     // --------------[------------------+----+--------------]--------------
0563     //                                  w    V
0564     //                                       = buf * 10^k
0565     //
0566     // ten_k represents a unit-in-the-last-place in the decimal representation
0567     // stored in buf.
0568     // Decrement buf by ten_k while this takes buf closer to w.
0569 
0570     // The tests are written in this order to avoid overflow in unsigned
0571     // integer arithmetic.
0572 
0573     while (rest < dist
0574             && delta - rest >= ten_k
0575             && (rest + ten_k < dist || dist - rest > rest + ten_k - dist))
0576     {
0577         JSON_ASSERT(buf[len - 1] != '0');
0578         buf[len - 1]--;
0579         rest += ten_k;
0580     }
0581 }
0582 
0583 /*!
0584 Generates V = buffer * 10^decimal_exponent, such that M- <= V <= M+.
0585 M- and M+ must be normalized and share the same exponent -60 <= e <= -32.
0586 */
0587 inline void grisu2_digit_gen(char* buffer, int& length, int& decimal_exponent,
0588                              diyfp M_minus, diyfp w, diyfp M_plus)
0589 {
0590     static_assert(kAlpha >= -60, "internal error");
0591     static_assert(kGamma <= -32, "internal error");
0592 
0593     // Generates the digits (and the exponent) of a decimal floating-point
0594     // number V = buffer * 10^decimal_exponent in the range [M-, M+]. The diyfp's
0595     // w, M- and M+ share the same exponent e, which satisfies alpha <= e <= gamma.
0596     //
0597     //               <--------------------------- delta ---->
0598     //                                  <---- dist --------->
0599     // --------------[------------------+-------------------]--------------
0600     //               M-                 w                   M+
0601     //
0602     // Grisu2 generates the digits of M+ from left to right and stops as soon as
0603     // V is in [M-,M+].
0604 
0605     JSON_ASSERT(M_plus.e >= kAlpha);
0606     JSON_ASSERT(M_plus.e <= kGamma);
0607 
0608     std::uint64_t delta = diyfp::sub(M_plus, M_minus).f; // (significand of (M+ - M-), implicit exponent is e)
0609     std::uint64_t dist  = diyfp::sub(M_plus, w      ).f; // (significand of (M+ - w ), implicit exponent is e)
0610 
0611     // Split M+ = f * 2^e into two parts p1 and p2 (note: e < 0):
0612     //
0613     //      M+ = f * 2^e
0614     //         = ((f div 2^-e) * 2^-e + (f mod 2^-e)) * 2^e
0615     //         = ((p1        ) * 2^-e + (p2        )) * 2^e
0616     //         = p1 + p2 * 2^e
0617 
0618     const diyfp one(std::uint64_t{1} << -M_plus.e, M_plus.e);
0619 
0620     auto p1 = static_cast<std::uint32_t>(M_plus.f >> -one.e); // p1 = f div 2^-e (Since -e >= 32, p1 fits into a 32-bit int.)
0621     std::uint64_t p2 = M_plus.f & (one.f - 1);                    // p2 = f mod 2^-e
0622 
0623     // 1)
0624     //
0625     // Generate the digits of the integral part p1 = d[n-1]...d[1]d[0]
0626 
0627     JSON_ASSERT(p1 > 0);
0628 
0629     std::uint32_t pow10{};
0630     const int k = find_largest_pow10(p1, pow10);
0631 
0632     //      10^(k-1) <= p1 < 10^k, pow10 = 10^(k-1)
0633     //
0634     //      p1 = (p1 div 10^(k-1)) * 10^(k-1) + (p1 mod 10^(k-1))
0635     //         = (d[k-1]         ) * 10^(k-1) + (p1 mod 10^(k-1))
0636     //
0637     //      M+ = p1                                             + p2 * 2^e
0638     //         = d[k-1] * 10^(k-1) + (p1 mod 10^(k-1))          + p2 * 2^e
0639     //         = d[k-1] * 10^(k-1) + ((p1 mod 10^(k-1)) * 2^-e + p2) * 2^e
0640     //         = d[k-1] * 10^(k-1) + (                         rest) * 2^e
0641     //
0642     // Now generate the digits d[n] of p1 from left to right (n = k-1,...,0)
0643     //
0644     //      p1 = d[k-1]...d[n] * 10^n + d[n-1]...d[0]
0645     //
0646     // but stop as soon as
0647     //
0648     //      rest * 2^e = (d[n-1]...d[0] * 2^-e + p2) * 2^e <= delta * 2^e
0649 
0650     int n = k;
0651     while (n > 0)
0652     {
0653         // Invariants:
0654         //      M+ = buffer * 10^n + (p1 + p2 * 2^e)    (buffer = 0 for n = k)
0655         //      pow10 = 10^(n-1) <= p1 < 10^n
0656         //
0657         const std::uint32_t d = p1 / pow10;  // d = p1 div 10^(n-1)
0658         const std::uint32_t r = p1 % pow10;  // r = p1 mod 10^(n-1)
0659         //
0660         //      M+ = buffer * 10^n + (d * 10^(n-1) + r) + p2 * 2^e
0661         //         = (buffer * 10 + d) * 10^(n-1) + (r + p2 * 2^e)
0662         //
0663         JSON_ASSERT(d <= 9);
0664         buffer[length++] = static_cast<char>('0' + d); // buffer := buffer * 10 + d
0665         //
0666         //      M+ = buffer * 10^(n-1) + (r + p2 * 2^e)
0667         //
0668         p1 = r;
0669         n--;
0670         //
0671         //      M+ = buffer * 10^n + (p1 + p2 * 2^e)
0672         //      pow10 = 10^n
0673         //
0674 
0675         // Now check if enough digits have been generated.
0676         // Compute
0677         //
0678         //      p1 + p2 * 2^e = (p1 * 2^-e + p2) * 2^e = rest * 2^e
0679         //
0680         // Note:
0681         // Since rest and delta share the same exponent e, it suffices to
0682         // compare the significands.
0683         const std::uint64_t rest = (std::uint64_t{p1} << -one.e) + p2;
0684         if (rest <= delta)
0685         {
0686             // V = buffer * 10^n, with M- <= V <= M+.
0687 
0688             decimal_exponent += n;
0689 
0690             // We may now just stop. But instead look if the buffer could be
0691             // decremented to bring V closer to w.
0692             //
0693             // pow10 = 10^n is now 1 ulp in the decimal representation V.
0694             // The rounding procedure works with diyfp's with an implicit
0695             // exponent of e.
0696             //
0697             //      10^n = (10^n * 2^-e) * 2^e = ulp * 2^e
0698             //
0699             const std::uint64_t ten_n = std::uint64_t{pow10} << -one.e;
0700             grisu2_round(buffer, length, dist, delta, rest, ten_n);
0701 
0702             return;
0703         }
0704 
0705         pow10 /= 10;
0706         //
0707         //      pow10 = 10^(n-1) <= p1 < 10^n
0708         // Invariants restored.
0709     }
0710 
0711     // 2)
0712     //
0713     // The digits of the integral part have been generated:
0714     //
0715     //      M+ = d[k-1]...d[1]d[0] + p2 * 2^e
0716     //         = buffer            + p2 * 2^e
0717     //
0718     // Now generate the digits of the fractional part p2 * 2^e.
0719     //
0720     // Note:
0721     // No decimal point is generated: the exponent is adjusted instead.
0722     //
0723     // p2 actually represents the fraction
0724     //
0725     //      p2 * 2^e
0726     //          = p2 / 2^-e
0727     //          = d[-1] / 10^1 + d[-2] / 10^2 + ...
0728     //
0729     // Now generate the digits d[-m] of p1 from left to right (m = 1,2,...)
0730     //
0731     //      p2 * 2^e = d[-1]d[-2]...d[-m] * 10^-m
0732     //                      + 10^-m * (d[-m-1] / 10^1 + d[-m-2] / 10^2 + ...)
0733     //
0734     // using
0735     //
0736     //      10^m * p2 = ((10^m * p2) div 2^-e) * 2^-e + ((10^m * p2) mod 2^-e)
0737     //                = (                   d) * 2^-e + (                   r)
0738     //
0739     // or
0740     //      10^m * p2 * 2^e = d + r * 2^e
0741     //
0742     // i.e.
0743     //
0744     //      M+ = buffer + p2 * 2^e
0745     //         = buffer + 10^-m * (d + r * 2^e)
0746     //         = (buffer * 10^m + d) * 10^-m + 10^-m * r * 2^e
0747     //
0748     // and stop as soon as 10^-m * r * 2^e <= delta * 2^e
0749 
0750     JSON_ASSERT(p2 > delta);
0751 
0752     int m = 0;
0753     for (;;)
0754     {
0755         // Invariant:
0756         //      M+ = buffer * 10^-m + 10^-m * (d[-m-1] / 10 + d[-m-2] / 10^2 + ...) * 2^e
0757         //         = buffer * 10^-m + 10^-m * (p2                                 ) * 2^e
0758         //         = buffer * 10^-m + 10^-m * (1/10 * (10 * p2)                   ) * 2^e
0759         //         = buffer * 10^-m + 10^-m * (1/10 * ((10*p2 div 2^-e) * 2^-e + (10*p2 mod 2^-e)) * 2^e
0760         //
0761         JSON_ASSERT(p2 <= (std::numeric_limits<std::uint64_t>::max)() / 10);
0762         p2 *= 10;
0763         const std::uint64_t d = p2 >> -one.e;     // d = (10 * p2) div 2^-e
0764         const std::uint64_t r = p2 & (one.f - 1); // r = (10 * p2) mod 2^-e
0765         //
0766         //      M+ = buffer * 10^-m + 10^-m * (1/10 * (d * 2^-e + r) * 2^e
0767         //         = buffer * 10^-m + 10^-m * (1/10 * (d + r * 2^e))
0768         //         = (buffer * 10 + d) * 10^(-m-1) + 10^(-m-1) * r * 2^e
0769         //
0770         JSON_ASSERT(d <= 9);
0771         buffer[length++] = static_cast<char>('0' + d); // buffer := buffer * 10 + d
0772         //
0773         //      M+ = buffer * 10^(-m-1) + 10^(-m-1) * r * 2^e
0774         //
0775         p2 = r;
0776         m++;
0777         //
0778         //      M+ = buffer * 10^-m + 10^-m * p2 * 2^e
0779         // Invariant restored.
0780 
0781         // Check if enough digits have been generated.
0782         //
0783         //      10^-m * p2 * 2^e <= delta * 2^e
0784         //              p2 * 2^e <= 10^m * delta * 2^e
0785         //                    p2 <= 10^m * delta
0786         delta *= 10;
0787         dist  *= 10;
0788         if (p2 <= delta)
0789         {
0790             break;
0791         }
0792     }
0793 
0794     // V = buffer * 10^-m, with M- <= V <= M+.
0795 
0796     decimal_exponent -= m;
0797 
0798     // 1 ulp in the decimal representation is now 10^-m.
0799     // Since delta and dist are now scaled by 10^m, we need to do the
0800     // same with ulp in order to keep the units in sync.
0801     //
0802     //      10^m * 10^-m = 1 = 2^-e * 2^e = ten_m * 2^e
0803     //
0804     const std::uint64_t ten_m = one.f;
0805     grisu2_round(buffer, length, dist, delta, p2, ten_m);
0806 
0807     // By construction this algorithm generates the shortest possible decimal
0808     // number (Loitsch, Theorem 6.2) which rounds back to w.
0809     // For an input number of precision p, at least
0810     //
0811     //      N = 1 + ceil(p * log_10(2))
0812     //
0813     // decimal digits are sufficient to identify all binary floating-point
0814     // numbers (Matula, "In-and-Out conversions").
0815     // This implies that the algorithm does not produce more than N decimal
0816     // digits.
0817     //
0818     //      N = 17 for p = 53 (IEEE double precision)
0819     //      N = 9  for p = 24 (IEEE single precision)
0820 }
0821 
0822 /*!
0823 v = buf * 10^decimal_exponent
0824 len is the length of the buffer (number of decimal digits)
0825 The buffer must be large enough, i.e. >= max_digits10.
0826 */
0827 JSON_HEDLEY_NON_NULL(1)
0828 inline void grisu2(char* buf, int& len, int& decimal_exponent,
0829                    diyfp m_minus, diyfp v, diyfp m_plus)
0830 {
0831     JSON_ASSERT(m_plus.e == m_minus.e);
0832     JSON_ASSERT(m_plus.e == v.e);
0833 
0834     //  --------(-----------------------+-----------------------)--------    (A)
0835     //          m-                      v                       m+
0836     //
0837     //  --------------------(-----------+-----------------------)--------    (B)
0838     //                      m-          v                       m+
0839     //
0840     // First scale v (and m- and m+) such that the exponent is in the range
0841     // [alpha, gamma].
0842 
0843     const cached_power cached = get_cached_power_for_binary_exponent(m_plus.e);
0844 
0845     const diyfp c_minus_k(cached.f, cached.e); // = c ~= 10^-k
0846 
0847     // The exponent of the products is = v.e + c_minus_k.e + q and is in the range [alpha,gamma]
0848     const diyfp w       = diyfp::mul(v,       c_minus_k);
0849     const diyfp w_minus = diyfp::mul(m_minus, c_minus_k);
0850     const diyfp w_plus  = diyfp::mul(m_plus,  c_minus_k);
0851 
0852     //  ----(---+---)---------------(---+---)---------------(---+---)----
0853     //          w-                      w                       w+
0854     //          = c*m-                  = c*v                   = c*m+
0855     //
0856     // diyfp::mul rounds its result and c_minus_k is approximated too. w, w- and
0857     // w+ are now off by a small amount.
0858     // In fact:
0859     //
0860     //      w - v * 10^k < 1 ulp
0861     //
0862     // To account for this inaccuracy, add resp. subtract 1 ulp.
0863     //
0864     //  --------+---[---------------(---+---)---------------]---+--------
0865     //          w-  M-                  w                   M+  w+
0866     //
0867     // Now any number in [M-, M+] (bounds included) will round to w when input,
0868     // regardless of how the input rounding algorithm breaks ties.
0869     //
0870     // And digit_gen generates the shortest possible such number in [M-, M+].
0871     // Note that this does not mean that Grisu2 always generates the shortest
0872     // possible number in the interval (m-, m+).
0873     const diyfp M_minus(w_minus.f + 1, w_minus.e);
0874     const diyfp M_plus (w_plus.f  - 1, w_plus.e );
0875 
0876     decimal_exponent = -cached.k; // = -(-k) = k
0877 
0878     grisu2_digit_gen(buf, len, decimal_exponent, M_minus, w, M_plus);
0879 }
0880 
0881 /*!
0882 v = buf * 10^decimal_exponent
0883 len is the length of the buffer (number of decimal digits)
0884 The buffer must be large enough, i.e. >= max_digits10.
0885 */
0886 template<typename FloatType>
0887 JSON_HEDLEY_NON_NULL(1)
0888 void grisu2(char* buf, int& len, int& decimal_exponent, FloatType value)
0889 {
0890     static_assert(diyfp::kPrecision >= std::numeric_limits<FloatType>::digits + 3,
0891                   "internal error: not enough precision");
0892 
0893     JSON_ASSERT(std::isfinite(value));
0894     JSON_ASSERT(value > 0);
0895 
0896     // If the neighbors (and boundaries) of 'value' are always computed for double-precision
0897     // numbers, all float's can be recovered using strtod (and strtof). However, the resulting
0898     // decimal representations are not exactly "short".
0899     //
0900     // The documentation for 'std::to_chars' (https://en.cppreference.com/w/cpp/utility/to_chars)
0901     // says "value is converted to a string as if by std::sprintf in the default ("C") locale"
0902     // and since sprintf promotes floats to doubles, I think this is exactly what 'std::to_chars'
0903     // does.
0904     // On the other hand, the documentation for 'std::to_chars' requires that "parsing the
0905     // representation using the corresponding std::from_chars function recovers value exactly". That
0906     // indicates that single precision floating-point numbers should be recovered using
0907     // 'std::strtof'.
0908     //
0909     // NB: If the neighbors are computed for single-precision numbers, there is a single float
0910     //     (7.0385307e-26f) which can't be recovered using strtod. The resulting double precision
0911     //     value is off by 1 ulp.
0912 #if 0
0913     const boundaries w = compute_boundaries(static_cast<double>(value));
0914 #else
0915     const boundaries w = compute_boundaries(value);
0916 #endif
0917 
0918     grisu2(buf, len, decimal_exponent, w.minus, w.w, w.plus);
0919 }
0920 
0921 /*!
0922 @brief appends a decimal representation of e to buf
0923 @return a pointer to the element following the exponent.
0924 @pre -1000 < e < 1000
0925 */
0926 JSON_HEDLEY_NON_NULL(1)
0927 JSON_HEDLEY_RETURNS_NON_NULL
0928 inline char* append_exponent(char* buf, int e)
0929 {
0930     JSON_ASSERT(e > -1000);
0931     JSON_ASSERT(e <  1000);
0932 
0933     if (e < 0)
0934     {
0935         e = -e;
0936         *buf++ = '-';
0937     }
0938     else
0939     {
0940         *buf++ = '+';
0941     }
0942 
0943     auto k = static_cast<std::uint32_t>(e);
0944     if (k < 10)
0945     {
0946         // Always print at least two digits in the exponent.
0947         // This is for compatibility with printf("%g").
0948         *buf++ = '0';
0949         *buf++ = static_cast<char>('0' + k);
0950     }
0951     else if (k < 100)
0952     {
0953         *buf++ = static_cast<char>('0' + k / 10);
0954         k %= 10;
0955         *buf++ = static_cast<char>('0' + k);
0956     }
0957     else
0958     {
0959         *buf++ = static_cast<char>('0' + k / 100);
0960         k %= 100;
0961         *buf++ = static_cast<char>('0' + k / 10);
0962         k %= 10;
0963         *buf++ = static_cast<char>('0' + k);
0964     }
0965 
0966     return buf;
0967 }
0968 
0969 /*!
0970 @brief prettify v = buf * 10^decimal_exponent
0971 
0972 If v is in the range [10^min_exp, 10^max_exp) it will be printed in fixed-point
0973 notation. Otherwise it will be printed in exponential notation.
0974 
0975 @pre min_exp < 0
0976 @pre max_exp > 0
0977 */
0978 JSON_HEDLEY_NON_NULL(1)
0979 JSON_HEDLEY_RETURNS_NON_NULL
0980 inline char* format_buffer(char* buf, int len, int decimal_exponent,
0981                            int min_exp, int max_exp)
0982 {
0983     JSON_ASSERT(min_exp < 0);
0984     JSON_ASSERT(max_exp > 0);
0985 
0986     const int k = len;
0987     const int n = len + decimal_exponent;
0988 
0989     // v = buf * 10^(n-k)
0990     // k is the length of the buffer (number of decimal digits)
0991     // n is the position of the decimal point relative to the start of the buffer.
0992 
0993     if (k <= n && n <= max_exp)
0994     {
0995         // digits[000]
0996         // len <= max_exp + 2
0997 
0998         std::memset(buf + k, '0', static_cast<size_t>(n) - static_cast<size_t>(k));
0999         // Make it look like a floating-point number (#362, #378)
1000         buf[n + 0] = '.';
1001         buf[n + 1] = '0';
1002         return buf + (static_cast<size_t>(n) + 2);
1003     }
1004 
1005     if (0 < n && n <= max_exp)
1006     {
1007         // dig.its
1008         // len <= max_digits10 + 1
1009 
1010         JSON_ASSERT(k > n);
1011 
1012         std::memmove(buf + (static_cast<size_t>(n) + 1), buf + n, static_cast<size_t>(k) - static_cast<size_t>(n));
1013         buf[n] = '.';
1014         return buf + (static_cast<size_t>(k) + 1U);
1015     }
1016 
1017     if (min_exp < n && n <= 0)
1018     {
1019         // 0.[000]digits
1020         // len <= 2 + (-min_exp - 1) + max_digits10
1021 
1022         std::memmove(buf + (2 + static_cast<size_t>(-n)), buf, static_cast<size_t>(k));
1023         buf[0] = '0';
1024         buf[1] = '.';
1025         std::memset(buf + 2, '0', static_cast<size_t>(-n));
1026         return buf + (2U + static_cast<size_t>(-n) + static_cast<size_t>(k));
1027     }
1028 
1029     if (k == 1)
1030     {
1031         // dE+123
1032         // len <= 1 + 5
1033 
1034         buf += 1;
1035     }
1036     else
1037     {
1038         // d.igitsE+123
1039         // len <= max_digits10 + 1 + 5
1040 
1041         std::memmove(buf + 2, buf + 1, static_cast<size_t>(k) - 1);
1042         buf[1] = '.';
1043         buf += 1 + static_cast<size_t>(k);
1044     }
1045 
1046     *buf++ = 'e';
1047     return append_exponent(buf, n - 1);
1048 }
1049 
1050 }  // namespace dtoa_impl
1051 
1052 /*!
1053 @brief generates a decimal representation of the floating-point number value in [first, last).
1054 
1055 The format of the resulting decimal representation is similar to printf's %g
1056 format. Returns an iterator pointing past-the-end of the decimal representation.
1057 
1058 @note The input number must be finite, i.e. NaN's and Inf's are not supported.
1059 @note The buffer must be large enough.
1060 @note The result is NOT null-terminated.
1061 */
1062 template<typename FloatType>
1063 JSON_HEDLEY_NON_NULL(1, 2)
1064 JSON_HEDLEY_RETURNS_NON_NULL
1065 char* to_chars(char* first, const char* last, FloatType value)
1066 {
1067     static_cast<void>(last); // maybe unused - fix warning
1068     JSON_ASSERT(std::isfinite(value));
1069 
1070     // Use signbit(value) instead of (value < 0) since signbit works for -0.
1071     if (std::signbit(value))
1072     {
1073         value = -value;
1074         *first++ = '-';
1075     }
1076 
1077 #ifdef __GNUC__
1078 #pragma GCC diagnostic push
1079 #pragma GCC diagnostic ignored "-Wfloat-equal"
1080 #endif
1081     if (value == 0) // +-0
1082     {
1083         *first++ = '0';
1084         // Make it look like a floating-point number (#362, #378)
1085         *first++ = '.';
1086         *first++ = '0';
1087         return first;
1088     }
1089 #ifdef __GNUC__
1090 #pragma GCC diagnostic pop
1091 #endif
1092 
1093     JSON_ASSERT(last - first >= std::numeric_limits<FloatType>::max_digits10);
1094 
1095     // Compute v = buffer * 10^decimal_exponent.
1096     // The decimal digits are stored in the buffer, which needs to be interpreted
1097     // as an unsigned decimal integer.
1098     // len is the length of the buffer, i.e. the number of decimal digits.
1099     int len = 0;
1100     int decimal_exponent = 0;
1101     dtoa_impl::grisu2(first, len, decimal_exponent, value);
1102 
1103     JSON_ASSERT(len <= std::numeric_limits<FloatType>::max_digits10);
1104 
1105     // Format the buffer like printf("%.*g", prec, value)
1106     constexpr int kMinExp = -4;
1107     // Use digits10 here to increase compatibility with version 2.
1108     constexpr int kMaxExp = std::numeric_limits<FloatType>::digits10;
1109 
1110     JSON_ASSERT(last - first >= kMaxExp + 2);
1111     JSON_ASSERT(last - first >= 2 + (-kMinExp - 1) + std::numeric_limits<FloatType>::max_digits10);
1112     JSON_ASSERT(last - first >= std::numeric_limits<FloatType>::max_digits10 + 6);
1113 
1114     return dtoa_impl::format_buffer(first, len, decimal_exponent, kMinExp, kMaxExp);
1115 }
1116 
1117 }  // namespace detail
1118 NLOHMANN_JSON_NAMESPACE_END