1#pragma once
2
3#include <array> // array
4#include <cassert> // assert
5#include <ciso646> // or, and, not
6#include <cmath> // signbit, isfinite
7#include <cstdint> // intN_t, uintN_t
8#include <cstring> // memcpy, memmove
9#include <limits> // numeric_limits
10#include <type_traits> // conditional
11#include <nlohmann/detail/macro_scope.hpp>
12
13namespace nlohmann
14{
15namespace detail
16{
17
18/*!
19@brief implements the Grisu2 algorithm for binary to decimal floating-point
20conversion.
21
22This implementation is a slightly modified version of the reference
23implementation which may be obtained from
24http://florian.loitsch.com/publications (bench.tar.gz).
25
26The code is distributed under the MIT license, Copyright (c) 2009 Florian Loitsch.
27
28For a detailed description of the algorithm see:
29
30[1] Loitsch, "Printing Floating-Point Numbers Quickly and Accurately with
31 Integers", Proceedings of the ACM SIGPLAN 2010 Conference on Programming
32 Language Design and Implementation, PLDI 2010
33[2] Burger, Dybvig, "Printing Floating-Point Numbers Quickly and Accurately",
34 Proceedings of the ACM SIGPLAN 1996 Conference on Programming Language
35 Design and Implementation, PLDI 1996
36*/
37namespace dtoa_impl
38{
39
40template <typename Target, typename Source>
41Target reinterpret_bits(const Source source)
42{
43 static_assert(sizeof(Target) == sizeof(Source), "size mismatch");
44
45 Target target;
46 std::memcpy(&target, &source, sizeof(Source));
47 return target;
48}
49
50struct diyfp // f * 2^e
51{
52 static constexpr int kPrecision = 64; // = q
53
54 std::uint64_t f = 0;
55 int e = 0;
56
57 constexpr diyfp(std::uint64_t f_, int e_) noexcept : f(f_), e(e_) {}
58
59 /*!
60 @brief returns x - y
61 @pre x.e == y.e and x.f >= y.f
62 */
63 static diyfp sub(const diyfp& x, const diyfp& y) noexcept
64 {
65 assert(x.e == y.e);
66 assert(x.f >= y.f);
67
68 return {x.f - y.f, x.e};
69 }
70
71 /*!
72 @brief returns x * y
73 @note The result is rounded. (Only the upper q bits are returned.)
74 */
75 static diyfp mul(const diyfp& x, const diyfp& y) noexcept
76 {
77 static_assert(kPrecision == 64, "internal error");
78
79 // Computes:
80 // f = round((x.f * y.f) / 2^q)
81 // e = x.e + y.e + q
82
83 // Emulate the 64-bit * 64-bit multiplication:
84 //
85 // p = u * v
86 // = (u_lo + 2^32 u_hi) (v_lo + 2^32 v_hi)
87 // = (u_lo v_lo ) + 2^32 ((u_lo v_hi ) + (u_hi v_lo )) + 2^64 (u_hi v_hi )
88 // = (p0 ) + 2^32 ((p1 ) + (p2 )) + 2^64 (p3 )
89 // = (p0_lo + 2^32 p0_hi) + 2^32 ((p1_lo + 2^32 p1_hi) + (p2_lo + 2^32 p2_hi)) + 2^64 (p3 )
90 // = (p0_lo ) + 2^32 (p0_hi + p1_lo + p2_lo ) + 2^64 (p1_hi + p2_hi + p3)
91 // = (p0_lo ) + 2^32 (Q ) + 2^64 (H )
92 // = (p0_lo ) + 2^32 (Q_lo + 2^32 Q_hi ) + 2^64 (H )
93 //
94 // (Since Q might be larger than 2^32 - 1)
95 //
96 // = (p0_lo + 2^32 Q_lo) + 2^64 (Q_hi + H)
97 //
98 // (Q_hi + H does not overflow a 64-bit int)
99 //
100 // = p_lo + 2^64 p_hi
101
102 const std::uint64_t u_lo = x.f & 0xFFFFFFFFu;
103 const std::uint64_t u_hi = x.f >> 32u;
104 const std::uint64_t v_lo = y.f & 0xFFFFFFFFu;
105 const std::uint64_t v_hi = y.f >> 32u;
106
107 const std::uint64_t p0 = u_lo * v_lo;
108 const std::uint64_t p1 = u_lo * v_hi;
109 const std::uint64_t p2 = u_hi * v_lo;
110 const std::uint64_t p3 = u_hi * v_hi;
111
112 const std::uint64_t p0_hi = p0 >> 32u;
113 const std::uint64_t p1_lo = p1 & 0xFFFFFFFFu;
114 const std::uint64_t p1_hi = p1 >> 32u;
115 const std::uint64_t p2_lo = p2 & 0xFFFFFFFFu;
116 const std::uint64_t p2_hi = p2 >> 32u;
117
118 std::uint64_t Q = p0_hi + p1_lo + p2_lo;
119
120 // The full product might now be computed as
121 //
122 // p_hi = p3 + p2_hi + p1_hi + (Q >> 32)
123 // p_lo = p0_lo + (Q << 32)
124 //
125 // But in this particular case here, the full p_lo is not required.
126 // Effectively we only need to add the highest bit in p_lo to p_hi (and
127 // Q_hi + 1 does not overflow).
128
129 Q += std::uint64_t{1} << (64u - 32u - 1u); // round, ties up
130
131 const std::uint64_t h = p3 + p2_hi + p1_hi + (Q >> 32u);
132
133 return {h, x.e + y.e + 64};
134 }
135
136 /*!
137 @brief normalize x such that the significand is >= 2^(q-1)
138 @pre x.f != 0
139 */
140 static diyfp normalize(diyfp x) noexcept
141 {
142 assert(x.f != 0);
143
144 while ((x.f >> 63u) == 0)
145 {
146 x.f <<= 1u;
147 x.e--;
148 }
149
150 return x;
151 }
152
153 /*!
154 @brief normalize x such that the result has the exponent E
155 @pre e >= x.e and the upper e - x.e bits of x.f must be zero.
156 */
157 static diyfp normalize_to(const diyfp& x, const int target_exponent) noexcept
158 {
159 const int delta = x.e - target_exponent;
160
161 assert(delta >= 0);
162 assert(((x.f << delta) >> delta) == x.f);
163
164 return {x.f << delta, target_exponent};
165 }
166};
167
168struct boundaries
169{
170 diyfp w;
171 diyfp minus;
172 diyfp plus;
173};
174
175/*!
176Compute the (normalized) diyfp representing the input number 'value' and its
177boundaries.
178
179@pre value must be finite and positive
180*/
181template <typename FloatType>
182boundaries compute_boundaries(FloatType value)
183{
184 assert(std::isfinite(value));
185 assert(value > 0);
186
187 // Convert the IEEE representation into a diyfp.
188 //
189 // If v is denormal:
190 // value = 0.F * 2^(1 - bias) = ( F) * 2^(1 - bias - (p-1))
191 // If v is normalized:
192 // value = 1.F * 2^(E - bias) = (2^(p-1) + F) * 2^(E - bias - (p-1))
193
194 static_assert(std::numeric_limits<FloatType>::is_iec559,
195 "internal error: dtoa_short requires an IEEE-754 floating-point implementation");
196
197 constexpr int kPrecision = std::numeric_limits<FloatType>::digits; // = p (includes the hidden bit)
198 constexpr int kBias = std::numeric_limits<FloatType>::max_exponent - 1 + (kPrecision - 1);
199 constexpr int kMinExp = 1 - kBias;
200 constexpr std::uint64_t kHiddenBit = std::uint64_t{1} << (kPrecision - 1); // = 2^(p-1)
201
202 using bits_type = typename std::conditional<kPrecision == 24, std::uint32_t, std::uint64_t >::type;
203
204 const std::uint64_t bits = reinterpret_bits<bits_type>(value);
205 const std::uint64_t E = bits >> (kPrecision - 1);
206 const std::uint64_t F = bits & (kHiddenBit - 1);
207
208 const bool is_denormal = E == 0;
209 const diyfp v = is_denormal
210 ? diyfp(F, kMinExp)
211 : diyfp(F + kHiddenBit, static_cast<int>(E) - kBias);
212
213 // Compute the boundaries m- and m+ of the floating-point value
214 // v = f * 2^e.
215 //
216 // Determine v- and v+, the floating-point predecessor and successor if v,
217 // respectively.
218 //
219 // v- = v - 2^e if f != 2^(p-1) or e == e_min (A)
220 // = v - 2^(e-1) if f == 2^(p-1) and e > e_min (B)
221 //
222 // v+ = v + 2^e
223 //
224 // Let m- = (v- + v) / 2 and m+ = (v + v+) / 2. All real numbers _strictly_
225 // between m- and m+ round to v, regardless of how the input rounding
226 // algorithm breaks ties.
227 //
228 // ---+-------------+-------------+-------------+-------------+--- (A)
229 // v- m- v m+ v+
230 //
231 // -----------------+------+------+-------------+-------------+--- (B)
232 // v- m- v m+ v+
233
234 const bool lower_boundary_is_closer = F == 0 and E > 1;
235 const diyfp m_plus = diyfp(2 * v.f + 1, v.e - 1);
236 const diyfp m_minus = lower_boundary_is_closer
237 ? diyfp(4 * v.f - 1, v.e - 2) // (B)
238 : diyfp(2 * v.f - 1, v.e - 1); // (A)
239
240 // Determine the normalized w+ = m+.
241 const diyfp w_plus = diyfp::normalize(m_plus);
242
243 // Determine w- = m- such that e_(w-) = e_(w+).
244 const diyfp w_minus = diyfp::normalize_to(m_minus, w_plus.e);
245
246 return {diyfp::normalize(v), w_minus, w_plus};
247}
248
249// Given normalized diyfp w, Grisu needs to find a (normalized) cached
250// power-of-ten c, such that the exponent of the product c * w = f * 2^e lies
251// within a certain range [alpha, gamma] (Definition 3.2 from [1])
252//
253// alpha <= e = e_c + e_w + q <= gamma
254//
255// or
256//
257// f_c * f_w * 2^alpha <= f_c 2^(e_c) * f_w 2^(e_w) * 2^q
258// <= f_c * f_w * 2^gamma
259//
260// Since c and w are normalized, i.e. 2^(q-1) <= f < 2^q, this implies
261//
262// 2^(q-1) * 2^(q-1) * 2^alpha <= c * w * 2^q < 2^q * 2^q * 2^gamma
263//
264// or
265//
266// 2^(q - 2 + alpha) <= c * w < 2^(q + gamma)
267//
268// The choice of (alpha,gamma) determines the size of the table and the form of
269// the digit generation procedure. Using (alpha,gamma)=(-60,-32) works out well
270// in practice:
271//
272// The idea is to cut the number c * w = f * 2^e into two parts, which can be
273// processed independently: An integral part p1, and a fractional part p2:
274//
275// f * 2^e = ( (f div 2^-e) * 2^-e + (f mod 2^-e) ) * 2^e
276// = (f div 2^-e) + (f mod 2^-e) * 2^e
277// = p1 + p2 * 2^e
278//
279// The conversion of p1 into decimal form requires a series of divisions and
280// modulos by (a power of) 10. These operations are faster for 32-bit than for
281// 64-bit integers, so p1 should ideally fit into a 32-bit integer. This can be
282// achieved by choosing
283//
284// -e >= 32 or e <= -32 := gamma
285//
286// In order to convert the fractional part
287//
288// p2 * 2^e = p2 / 2^-e = d[-1] / 10^1 + d[-2] / 10^2 + ...
289//
290// into decimal form, the fraction is repeatedly multiplied by 10 and the digits
291// d[-i] are extracted in order:
292//
293// (10 * p2) div 2^-e = d[-1]
294// (10 * p2) mod 2^-e = d[-2] / 10^1 + ...
295//
296// The multiplication by 10 must not overflow. It is sufficient to choose
297//
298// 10 * p2 < 16 * p2 = 2^4 * p2 <= 2^64.
299//
300// Since p2 = f mod 2^-e < 2^-e,
301//
302// -e <= 60 or e >= -60 := alpha
303
304constexpr int kAlpha = -60;
305constexpr int kGamma = -32;
306
307struct cached_power // c = f * 2^e ~= 10^k
308{
309 std::uint64_t f;
310 int e;
311 int k;
312};
313
314/*!
315For a normalized diyfp w = f * 2^e, this function returns a (normalized) cached
316power-of-ten c = f_c * 2^e_c, such that the exponent of the product w * c
317satisfies (Definition 3.2 from [1])
318
319 alpha <= e_c + e + q <= gamma.
320*/
321inline cached_power get_cached_power_for_binary_exponent(int e)
322{
323 // Now
324 //
325 // alpha <= e_c + e + q <= gamma (1)
326 // ==> f_c * 2^alpha <= c * 2^e * 2^q
327 //
328 // and since the c's are normalized, 2^(q-1) <= f_c,
329 //
330 // ==> 2^(q - 1 + alpha) <= c * 2^(e + q)
331 // ==> 2^(alpha - e - 1) <= c
332 //
333 // If c were an exact power of ten, i.e. c = 10^k, one may determine k as
334 //
335 // k = ceil( log_10( 2^(alpha - e - 1) ) )
336 // = ceil( (alpha - e - 1) * log_10(2) )
337 //
338 // From the paper:
339 // "In theory the result of the procedure could be wrong since c is rounded,
340 // and the computation itself is approximated [...]. In practice, however,
341 // this simple function is sufficient."
342 //
343 // For IEEE double precision floating-point numbers converted into
344 // normalized diyfp's w = f * 2^e, with q = 64,
345 //
346 // e >= -1022 (min IEEE exponent)
347 // -52 (p - 1)
348 // -52 (p - 1, possibly normalize denormal IEEE numbers)
349 // -11 (normalize the diyfp)
350 // = -1137
351 //
352 // and
353 //
354 // e <= +1023 (max IEEE exponent)
355 // -52 (p - 1)
356 // -11 (normalize the diyfp)
357 // = 960
358 //
359 // This binary exponent range [-1137,960] results in a decimal exponent
360 // range [-307,324]. One does not need to store a cached power for each
361 // k in this range. For each such k it suffices to find a cached power
362 // such that the exponent of the product lies in [alpha,gamma].
363 // This implies that the difference of the decimal exponents of adjacent
364 // table entries must be less than or equal to
365 //
366 // floor( (gamma - alpha) * log_10(2) ) = 8.
367 //
368 // (A smaller distance gamma-alpha would require a larger table.)
369
370 // NB:
371 // Actually this function returns c, such that -60 <= e_c + e + 64 <= -34.
372
373 constexpr int kCachedPowersMinDecExp = -300;
374 constexpr int kCachedPowersDecStep = 8;
375
376 static constexpr std::array<cached_power, 79> kCachedPowers =
377 {
378 {
379 { 0xAB70FE17C79AC6CA, -1060, -300 },
380 { 0xFF77B1FCBEBCDC4F, -1034, -292 },
381 { 0xBE5691EF416BD60C, -1007, -284 },
382 { 0x8DD01FAD907FFC3C, -980, -276 },
383 { 0xD3515C2831559A83, -954, -268 },
384 { 0x9D71AC8FADA6C9B5, -927, -260 },
385 { 0xEA9C227723EE8BCB, -901, -252 },
386 { 0xAECC49914078536D, -874, -244 },
387 { 0x823C12795DB6CE57, -847, -236 },
388 { 0xC21094364DFB5637, -821, -228 },
389 { 0x9096EA6F3848984F, -794, -220 },
390 { 0xD77485CB25823AC7, -768, -212 },
391 { 0xA086CFCD97BF97F4, -741, -204 },
392 { 0xEF340A98172AACE5, -715, -196 },
393 { 0xB23867FB2A35B28E, -688, -188 },
394 { 0x84C8D4DFD2C63F3B, -661, -180 },
395 { 0xC5DD44271AD3CDBA, -635, -172 },
396 { 0x936B9FCEBB25C996, -608, -164 },
397 { 0xDBAC6C247D62A584, -582, -156 },
398 { 0xA3AB66580D5FDAF6, -555, -148 },
399 { 0xF3E2F893DEC3F126, -529, -140 },
400 { 0xB5B5ADA8AAFF80B8, -502, -132 },
401 { 0x87625F056C7C4A8B, -475, -124 },
402 { 0xC9BCFF6034C13053, -449, -116 },
403 { 0x964E858C91BA2655, -422, -108 },
404 { 0xDFF9772470297EBD, -396, -100 },
405 { 0xA6DFBD9FB8E5B88F, -369, -92 },
406 { 0xF8A95FCF88747D94, -343, -84 },
407 { 0xB94470938FA89BCF, -316, -76 },
408 { 0x8A08F0F8BF0F156B, -289, -68 },
409 { 0xCDB02555653131B6, -263, -60 },
410 { 0x993FE2C6D07B7FAC, -236, -52 },
411 { 0xE45C10C42A2B3B06, -210, -44 },
412 { 0xAA242499697392D3, -183, -36 },
413 { 0xFD87B5F28300CA0E, -157, -28 },
414 { 0xBCE5086492111AEB, -130, -20 },
415 { 0x8CBCCC096F5088CC, -103, -12 },
416 { 0xD1B71758E219652C, -77, -4 },
417 { 0x9C40000000000000, -50, 4 },
418 { 0xE8D4A51000000000, -24, 12 },
419 { 0xAD78EBC5AC620000, 3, 20 },
420 { 0x813F3978F8940984, 30, 28 },
421 { 0xC097CE7BC90715B3, 56, 36 },
422 { 0x8F7E32CE7BEA5C70, 83, 44 },
423 { 0xD5D238A4ABE98068, 109, 52 },
424 { 0x9F4F2726179A2245, 136, 60 },
425 { 0xED63A231D4C4FB27, 162, 68 },
426 { 0xB0DE65388CC8ADA8, 189, 76 },
427 { 0x83C7088E1AAB65DB, 216, 84 },
428 { 0xC45D1DF942711D9A, 242, 92 },
429 { 0x924D692CA61BE758, 269, 100 },
430 { 0xDA01EE641A708DEA, 295, 108 },
431 { 0xA26DA3999AEF774A, 322, 116 },
432 { 0xF209787BB47D6B85, 348, 124 },
433 { 0xB454E4A179DD1877, 375, 132 },
434 { 0x865B86925B9BC5C2, 402, 140 },
435 { 0xC83553C5C8965D3D, 428, 148 },
436 { 0x952AB45CFA97A0B3, 455, 156 },
437 { 0xDE469FBD99A05FE3, 481, 164 },
438 { 0xA59BC234DB398C25, 508, 172 },
439 { 0xF6C69A72A3989F5C, 534, 180 },
440 { 0xB7DCBF5354E9BECE, 561, 188 },
441 { 0x88FCF317F22241E2, 588, 196 },
442 { 0xCC20CE9BD35C78A5, 614, 204 },
443 { 0x98165AF37B2153DF, 641, 212 },
444 { 0xE2A0B5DC971F303A, 667, 220 },
445 { 0xA8D9D1535CE3B396, 694, 228 },
446 { 0xFB9B7CD9A4A7443C, 720, 236 },
447 { 0xBB764C4CA7A44410, 747, 244 },
448 { 0x8BAB8EEFB6409C1A, 774, 252 },
449 { 0xD01FEF10A657842C, 800, 260 },
450 { 0x9B10A4E5E9913129, 827, 268 },
451 { 0xE7109BFBA19C0C9D, 853, 276 },
452 { 0xAC2820D9623BF429, 880, 284 },
453 { 0x80444B5E7AA7CF85, 907, 292 },
454 { 0xBF21E44003ACDD2D, 933, 300 },
455 { 0x8E679C2F5E44FF8F, 960, 308 },
456 { 0xD433179D9C8CB841, 986, 316 },
457 { 0x9E19DB92B4E31BA9, 1013, 324 },
458 }
459 };
460
461 // This computation gives exactly the same results for k as
462 // k = ceil((kAlpha - e - 1) * 0.30102999566398114)
463 // for |e| <= 1500, but doesn't require floating-point operations.
464 // NB: log_10(2) ~= 78913 / 2^18
465 assert(e >= -1500);
466 assert(e <= 1500);
467 const int f = kAlpha - e - 1;
468 const int k = (f * 78913) / (1 << 18) + static_cast<int>(f > 0);
469
470 const int index = (-kCachedPowersMinDecExp + k + (kCachedPowersDecStep - 1)) / kCachedPowersDecStep;
471 assert(index >= 0);
472 assert(static_cast<std::size_t>(index) < kCachedPowers.size());
473
474 const cached_power cached = kCachedPowers[static_cast<std::size_t>(index)];
475 assert(kAlpha <= cached.e + e + 64);
476 assert(kGamma >= cached.e + e + 64);
477
478 return cached;
479}
480
481/*!
482For n != 0, returns k, such that pow10 := 10^(k-1) <= n < 10^k.
483For n == 0, returns 1 and sets pow10 := 1.
484*/
485inline int find_largest_pow10(const std::uint32_t n, std::uint32_t& pow10)
486{
487 // LCOV_EXCL_START
488 if (n >= 1000000000)
489 {
490 pow10 = 1000000000;
491 return 10;
492 }
493 // LCOV_EXCL_STOP
494 else if (n >= 100000000)
495 {
496 pow10 = 100000000;
497 return 9;
498 }
499 else if (n >= 10000000)
500 {
501 pow10 = 10000000;
502 return 8;
503 }
504 else if (n >= 1000000)
505 {
506 pow10 = 1000000;
507 return 7;
508 }
509 else if (n >= 100000)
510 {
511 pow10 = 100000;
512 return 6;
513 }
514 else if (n >= 10000)
515 {
516 pow10 = 10000;
517 return 5;
518 }
519 else if (n >= 1000)
520 {
521 pow10 = 1000;
522 return 4;
523 }
524 else if (n >= 100)
525 {
526 pow10 = 100;
527 return 3;
528 }
529 else if (n >= 10)
530 {
531 pow10 = 10;
532 return 2;
533 }
534 else
535 {
536 pow10 = 1;
537 return 1;
538 }
539}
540
541inline void grisu2_round(char* buf, int len, std::uint64_t dist, std::uint64_t delta,
542 std::uint64_t rest, std::uint64_t ten_k)
543{
544 assert(len >= 1);
545 assert(dist <= delta);
546 assert(rest <= delta);
547 assert(ten_k > 0);
548
549 // <--------------------------- delta ---->
550 // <---- dist --------->
551 // --------------[------------------+-------------------]--------------
552 // M- w M+
553 //
554 // ten_k
555 // <------>
556 // <---- rest ---->
557 // --------------[------------------+----+--------------]--------------
558 // w V
559 // = buf * 10^k
560 //
561 // ten_k represents a unit-in-the-last-place in the decimal representation
562 // stored in buf.
563 // Decrement buf by ten_k while this takes buf closer to w.
564
565 // The tests are written in this order to avoid overflow in unsigned
566 // integer arithmetic.
567
568 while (rest < dist
569 and delta - rest >= ten_k
570 and (rest + ten_k < dist or dist - rest > rest + ten_k - dist))
571 {
572 assert(buf[len - 1] != '0');
573 buf[len - 1]--;
574 rest += ten_k;
575 }
576}
577
578/*!
579Generates V = buffer * 10^decimal_exponent, such that M- <= V <= M+.
580M- and M+ must be normalized and share the same exponent -60 <= e <= -32.
581*/
582inline void grisu2_digit_gen(char* buffer, int& length, int& decimal_exponent,
583 diyfp M_minus, diyfp w, diyfp M_plus)
584{
585 static_assert(kAlpha >= -60, "internal error");
586 static_assert(kGamma <= -32, "internal error");
587
588 // Generates the digits (and the exponent) of a decimal floating-point
589 // number V = buffer * 10^decimal_exponent in the range [M-, M+]. The diyfp's
590 // w, M- and M+ share the same exponent e, which satisfies alpha <= e <= gamma.
591 //
592 // <--------------------------- delta ---->
593 // <---- dist --------->
594 // --------------[------------------+-------------------]--------------
595 // M- w M+
596 //
597 // Grisu2 generates the digits of M+ from left to right and stops as soon as
598 // V is in [M-,M+].
599
600 assert(M_plus.e >= kAlpha);
601 assert(M_plus.e <= kGamma);
602
603 std::uint64_t delta = diyfp::sub(M_plus, M_minus).f; // (significand of (M+ - M-), implicit exponent is e)
604 std::uint64_t dist = diyfp::sub(M_plus, w ).f; // (significand of (M+ - w ), implicit exponent is e)
605
606 // Split M+ = f * 2^e into two parts p1 and p2 (note: e < 0):
607 //
608 // M+ = f * 2^e
609 // = ((f div 2^-e) * 2^-e + (f mod 2^-e)) * 2^e
610 // = ((p1 ) * 2^-e + (p2 )) * 2^e
611 // = p1 + p2 * 2^e
612
613 const diyfp one(std::uint64_t{1} << -M_plus.e, M_plus.e);
614
615 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.)
616 std::uint64_t p2 = M_plus.f & (one.f - 1); // p2 = f mod 2^-e
617
618 // 1)
619 //
620 // Generate the digits of the integral part p1 = d[n-1]...d[1]d[0]
621
622 assert(p1 > 0);
623
624 std::uint32_t pow10;
625 const int k = find_largest_pow10(p1, pow10);
626
627 // 10^(k-1) <= p1 < 10^k, pow10 = 10^(k-1)
628 //
629 // p1 = (p1 div 10^(k-1)) * 10^(k-1) + (p1 mod 10^(k-1))
630 // = (d[k-1] ) * 10^(k-1) + (p1 mod 10^(k-1))
631 //
632 // M+ = p1 + p2 * 2^e
633 // = d[k-1] * 10^(k-1) + (p1 mod 10^(k-1)) + p2 * 2^e
634 // = d[k-1] * 10^(k-1) + ((p1 mod 10^(k-1)) * 2^-e + p2) * 2^e
635 // = d[k-1] * 10^(k-1) + ( rest) * 2^e
636 //
637 // Now generate the digits d[n] of p1 from left to right (n = k-1,...,0)
638 //
639 // p1 = d[k-1]...d[n] * 10^n + d[n-1]...d[0]
640 //
641 // but stop as soon as
642 //
643 // rest * 2^e = (d[n-1]...d[0] * 2^-e + p2) * 2^e <= delta * 2^e
644
645 int n = k;
646 while (n > 0)
647 {
648 // Invariants:
649 // M+ = buffer * 10^n + (p1 + p2 * 2^e) (buffer = 0 for n = k)
650 // pow10 = 10^(n-1) <= p1 < 10^n
651 //
652 const std::uint32_t d = p1 / pow10; // d = p1 div 10^(n-1)
653 const std::uint32_t r = p1 % pow10; // r = p1 mod 10^(n-1)
654 //
655 // M+ = buffer * 10^n + (d * 10^(n-1) + r) + p2 * 2^e
656 // = (buffer * 10 + d) * 10^(n-1) + (r + p2 * 2^e)
657 //
658 assert(d <= 9);
659 buffer[length++] = static_cast<char>('0' + d); // buffer := buffer * 10 + d
660 //
661 // M+ = buffer * 10^(n-1) + (r + p2 * 2^e)
662 //
663 p1 = r;
664 n--;
665 //
666 // M+ = buffer * 10^n + (p1 + p2 * 2^e)
667 // pow10 = 10^n
668 //
669
670 // Now check if enough digits have been generated.
671 // Compute
672 //
673 // p1 + p2 * 2^e = (p1 * 2^-e + p2) * 2^e = rest * 2^e
674 //
675 // Note:
676 // Since rest and delta share the same exponent e, it suffices to
677 // compare the significands.
678 const std::uint64_t rest = (std::uint64_t{p1} << -one.e) + p2;
679 if (rest <= delta)
680 {
681 // V = buffer * 10^n, with M- <= V <= M+.
682
683 decimal_exponent += n;
684
685 // We may now just stop. But instead look if the buffer could be
686 // decremented to bring V closer to w.
687 //
688 // pow10 = 10^n is now 1 ulp in the decimal representation V.
689 // The rounding procedure works with diyfp's with an implicit
690 // exponent of e.
691 //
692 // 10^n = (10^n * 2^-e) * 2^e = ulp * 2^e
693 //
694 const std::uint64_t ten_n = std::uint64_t{pow10} << -one.e;
695 grisu2_round(buffer, length, dist, delta, rest, ten_n);
696
697 return;
698 }
699
700 pow10 /= 10;
701 //
702 // pow10 = 10^(n-1) <= p1 < 10^n
703 // Invariants restored.
704 }
705
706 // 2)
707 //
708 // The digits of the integral part have been generated:
709 //
710 // M+ = d[k-1]...d[1]d[0] + p2 * 2^e
711 // = buffer + p2 * 2^e
712 //
713 // Now generate the digits of the fractional part p2 * 2^e.
714 //
715 // Note:
716 // No decimal point is generated: the exponent is adjusted instead.
717 //
718 // p2 actually represents the fraction
719 //
720 // p2 * 2^e
721 // = p2 / 2^-e
722 // = d[-1] / 10^1 + d[-2] / 10^2 + ...
723 //
724 // Now generate the digits d[-m] of p1 from left to right (m = 1,2,...)
725 //
726 // p2 * 2^e = d[-1]d[-2]...d[-m] * 10^-m
727 // + 10^-m * (d[-m-1] / 10^1 + d[-m-2] / 10^2 + ...)
728 //
729 // using
730 //
731 // 10^m * p2 = ((10^m * p2) div 2^-e) * 2^-e + ((10^m * p2) mod 2^-e)
732 // = ( d) * 2^-e + ( r)
733 //
734 // or
735 // 10^m * p2 * 2^e = d + r * 2^e
736 //
737 // i.e.
738 //
739 // M+ = buffer + p2 * 2^e
740 // = buffer + 10^-m * (d + r * 2^e)
741 // = (buffer * 10^m + d) * 10^-m + 10^-m * r * 2^e
742 //
743 // and stop as soon as 10^-m * r * 2^e <= delta * 2^e
744
745 assert(p2 > delta);
746
747 int m = 0;
748 for (;;)
749 {
750 // Invariant:
751 // M+ = buffer * 10^-m + 10^-m * (d[-m-1] / 10 + d[-m-2] / 10^2 + ...) * 2^e
752 // = buffer * 10^-m + 10^-m * (p2 ) * 2^e
753 // = buffer * 10^-m + 10^-m * (1/10 * (10 * p2) ) * 2^e
754 // = buffer * 10^-m + 10^-m * (1/10 * ((10*p2 div 2^-e) * 2^-e + (10*p2 mod 2^-e)) * 2^e
755 //
756 assert(p2 <= (std::numeric_limits<std::uint64_t>::max)() / 10);
757 p2 *= 10;
758 const std::uint64_t d = p2 >> -one.e; // d = (10 * p2) div 2^-e
759 const std::uint64_t r = p2 & (one.f - 1); // r = (10 * p2) mod 2^-e
760 //
761 // M+ = buffer * 10^-m + 10^-m * (1/10 * (d * 2^-e + r) * 2^e
762 // = buffer * 10^-m + 10^-m * (1/10 * (d + r * 2^e))
763 // = (buffer * 10 + d) * 10^(-m-1) + 10^(-m-1) * r * 2^e
764 //
765 assert(d <= 9);
766 buffer[length++] = static_cast<char>('0' + d); // buffer := buffer * 10 + d
767 //
768 // M+ = buffer * 10^(-m-1) + 10^(-m-1) * r * 2^e
769 //
770 p2 = r;
771 m++;
772 //
773 // M+ = buffer * 10^-m + 10^-m * p2 * 2^e
774 // Invariant restored.
775
776 // Check if enough digits have been generated.
777 //
778 // 10^-m * p2 * 2^e <= delta * 2^e
779 // p2 * 2^e <= 10^m * delta * 2^e
780 // p2 <= 10^m * delta
781 delta *= 10;
782 dist *= 10;
783 if (p2 <= delta)
784 {
785 break;
786 }
787 }
788
789 // V = buffer * 10^-m, with M- <= V <= M+.
790
791 decimal_exponent -= m;
792
793 // 1 ulp in the decimal representation is now 10^-m.
794 // Since delta and dist are now scaled by 10^m, we need to do the
795 // same with ulp in order to keep the units in sync.
796 //
797 // 10^m * 10^-m = 1 = 2^-e * 2^e = ten_m * 2^e
798 //
799 const std::uint64_t ten_m = one.f;
800 grisu2_round(buffer, length, dist, delta, p2, ten_m);
801
802 // By construction this algorithm generates the shortest possible decimal
803 // number (Loitsch, Theorem 6.2) which rounds back to w.
804 // For an input number of precision p, at least
805 //
806 // N = 1 + ceil(p * log_10(2))
807 //
808 // decimal digits are sufficient to identify all binary floating-point
809 // numbers (Matula, "In-and-Out conversions").
810 // This implies that the algorithm does not produce more than N decimal
811 // digits.
812 //
813 // N = 17 for p = 53 (IEEE double precision)
814 // N = 9 for p = 24 (IEEE single precision)
815}
816
817/*!
818v = buf * 10^decimal_exponent
819len is the length of the buffer (number of decimal digits)
820The buffer must be large enough, i.e. >= max_digits10.
821*/
822JSON_HEDLEY_NON_NULL(1)
823inline void grisu2(char* buf, int& len, int& decimal_exponent,
824 diyfp m_minus, diyfp v, diyfp m_plus)
825{
826 assert(m_plus.e == m_minus.e);
827 assert(m_plus.e == v.e);
828
829 // --------(-----------------------+-----------------------)-------- (A)
830 // m- v m+
831 //
832 // --------------------(-----------+-----------------------)-------- (B)
833 // m- v m+
834 //
835 // First scale v (and m- and m+) such that the exponent is in the range
836 // [alpha, gamma].
837
838 const cached_power cached = get_cached_power_for_binary_exponent(m_plus.e);
839
840 const diyfp c_minus_k(cached.f, cached.e); // = c ~= 10^-k
841
842 // The exponent of the products is = v.e + c_minus_k.e + q and is in the range [alpha,gamma]
843 const diyfp w = diyfp::mul(v, c_minus_k);
844 const diyfp w_minus = diyfp::mul(m_minus, c_minus_k);
845 const diyfp w_plus = diyfp::mul(m_plus, c_minus_k);
846
847 // ----(---+---)---------------(---+---)---------------(---+---)----
848 // w- w w+
849 // = c*m- = c*v = c*m+
850 //
851 // diyfp::mul rounds its result and c_minus_k is approximated too. w, w- and
852 // w+ are now off by a small amount.
853 // In fact:
854 //
855 // w - v * 10^k < 1 ulp
856 //
857 // To account for this inaccuracy, add resp. subtract 1 ulp.
858 //
859 // --------+---[---------------(---+---)---------------]---+--------
860 // w- M- w M+ w+
861 //
862 // Now any number in [M-, M+] (bounds included) will round to w when input,
863 // regardless of how the input rounding algorithm breaks ties.
864 //
865 // And digit_gen generates the shortest possible such number in [M-, M+].
866 // Note that this does not mean that Grisu2 always generates the shortest
867 // possible number in the interval (m-, m+).
868 const diyfp M_minus(w_minus.f + 1, w_minus.e);
869 const diyfp M_plus (w_plus.f - 1, w_plus.e );
870
871 decimal_exponent = -cached.k; // = -(-k) = k
872
873 grisu2_digit_gen(buf, len, decimal_exponent, M_minus, w, M_plus);
874}
875
876/*!
877v = buf * 10^decimal_exponent
878len is the length of the buffer (number of decimal digits)
879The buffer must be large enough, i.e. >= max_digits10.
880*/
881template <typename FloatType>
882JSON_HEDLEY_NON_NULL(1)
883void grisu2(char* buf, int& len, int& decimal_exponent, FloatType value)
884{
885 static_assert(diyfp::kPrecision >= std::numeric_limits<FloatType>::digits + 3,
886 "internal error: not enough precision");
887
888 assert(std::isfinite(value));
889 assert(value > 0);
890
891 // If the neighbors (and boundaries) of 'value' are always computed for double-precision
892 // numbers, all float's can be recovered using strtod (and strtof). However, the resulting
893 // decimal representations are not exactly "short".
894 //
895 // The documentation for 'std::to_chars' (https://en.cppreference.com/w/cpp/utility/to_chars)
896 // says "value is converted to a string as if by std::sprintf in the default ("C") locale"
897 // and since sprintf promotes float's to double's, I think this is exactly what 'std::to_chars'
898 // does.
899 // On the other hand, the documentation for 'std::to_chars' requires that "parsing the
900 // representation using the corresponding std::from_chars function recovers value exactly". That
901 // indicates that single precision floating-point numbers should be recovered using
902 // 'std::strtof'.
903 //
904 // NB: If the neighbors are computed for single-precision numbers, there is a single float
905 // (7.0385307e-26f) which can't be recovered using strtod. The resulting double precision
906 // value is off by 1 ulp.
907#if 0
908 const boundaries w = compute_boundaries(static_cast<double>(value));
909#else
910 const boundaries w = compute_boundaries(value);
911#endif
912
913 grisu2(buf, len, decimal_exponent, w.minus, w.w, w.plus);
914}
915
916/*!
917@brief appends a decimal representation of e to buf
918@return a pointer to the element following the exponent.
919@pre -1000 < e < 1000
920*/
921JSON_HEDLEY_NON_NULL(1)
922JSON_HEDLEY_RETURNS_NON_NULL
923inline char* append_exponent(char* buf, int e)
924{
925 assert(e > -1000);
926 assert(e < 1000);
927
928 if (e < 0)
929 {
930 e = -e;
931 *buf++ = '-';
932 }
933 else
934 {
935 *buf++ = '+';
936 }
937
938 auto k = static_cast<std::uint32_t>(e);
939 if (k < 10)
940 {
941 // Always print at least two digits in the exponent.
942 // This is for compatibility with printf("%g").
943 *buf++ = '0';
944 *buf++ = static_cast<char>('0' + k);
945 }
946 else if (k < 100)
947 {
948 *buf++ = static_cast<char>('0' + k / 10);
949 k %= 10;
950 *buf++ = static_cast<char>('0' + k);
951 }
952 else
953 {
954 *buf++ = static_cast<char>('0' + k / 100);
955 k %= 100;
956 *buf++ = static_cast<char>('0' + k / 10);
957 k %= 10;
958 *buf++ = static_cast<char>('0' + k);
959 }
960
961 return buf;
962}
963
964/*!
965@brief prettify v = buf * 10^decimal_exponent
966
967If v is in the range [10^min_exp, 10^max_exp) it will be printed in fixed-point
968notation. Otherwise it will be printed in exponential notation.
969
970@pre min_exp < 0
971@pre max_exp > 0
972*/
973JSON_HEDLEY_NON_NULL(1)
974JSON_HEDLEY_RETURNS_NON_NULL
975inline char* format_buffer(char* buf, int len, int decimal_exponent,
976 int min_exp, int max_exp)
977{
978 assert(min_exp < 0);
979 assert(max_exp > 0);
980
981 const int k = len;
982 const int n = len + decimal_exponent;
983
984 // v = buf * 10^(n-k)
985 // k is the length of the buffer (number of decimal digits)
986 // n is the position of the decimal point relative to the start of the buffer.
987
988 if (k <= n and n <= max_exp)
989 {
990 // digits[000]
991 // len <= max_exp + 2
992
993 std::memset(buf + k, '0', static_cast<size_t>(n - k));
994 // Make it look like a floating-point number (#362, #378)
995 buf[n + 0] = '.';
996 buf[n + 1] = '0';
997 return buf + (n + 2);
998 }
999
1000 if (0 < n and n <= max_exp)
1001 {
1002 // dig.its
1003 // len <= max_digits10 + 1
1004
1005 assert(k > n);
1006
1007 std::memmove(buf + (n + 1), buf + n, static_cast<size_t>(k - n));
1008 buf[n] = '.';
1009 return buf + (k + 1);
1010 }
1011
1012 if (min_exp < n and n <= 0)
1013 {
1014 // 0.[000]digits
1015 // len <= 2 + (-min_exp - 1) + max_digits10
1016
1017 std::memmove(buf + (2 + -n), buf, static_cast<size_t>(k));
1018 buf[0] = '0';
1019 buf[1] = '.';
1020 std::memset(buf + 2, '0', static_cast<size_t>(-n));
1021 return buf + (2 + (-n) + k);
1022 }
1023
1024 if (k == 1)
1025 {
1026 // dE+123
1027 // len <= 1 + 5
1028
1029 buf += 1;
1030 }
1031 else
1032 {
1033 // d.igitsE+123
1034 // len <= max_digits10 + 1 + 5
1035
1036 std::memmove(buf + 2, buf + 1, static_cast<size_t>(k - 1));
1037 buf[1] = '.';
1038 buf += 1 + k;
1039 }
1040
1041 *buf++ = 'e';
1042 return append_exponent(buf, n - 1);
1043}
1044
1045} // namespace dtoa_impl
1046
1047/*!
1048@brief generates a decimal representation of the floating-point number value in [first, last).
1049
1050The format of the resulting decimal representation is similar to printf's %g
1051format. Returns an iterator pointing past-the-end of the decimal representation.
1052
1053@note The input number must be finite, i.e. NaN's and Inf's are not supported.
1054@note The buffer must be large enough.
1055@note The result is NOT null-terminated.
1056*/
1057template <typename FloatType>
1058JSON_HEDLEY_NON_NULL(1, 2)
1059JSON_HEDLEY_RETURNS_NON_NULL
1060char* to_chars(char* first, const char* last, FloatType value)
1061{
1062 static_cast<void>(last); // maybe unused - fix warning
1063 assert(std::isfinite(value));
1064
1065 // Use signbit(value) instead of (value < 0) since signbit works for -0.
1066 if (std::signbit(value))
1067 {
1068 value = -value;
1069 *first++ = '-';
1070 }
1071
1072 if (value == 0) // +-0
1073 {
1074 *first++ = '0';
1075 // Make it look like a floating-point number (#362, #378)
1076 *first++ = '.';
1077 *first++ = '0';
1078 return first;
1079 }
1080
1081 assert(last - first >= std::numeric_limits<FloatType>::max_digits10);
1082
1083 // Compute v = buffer * 10^decimal_exponent.
1084 // The decimal digits are stored in the buffer, which needs to be interpreted
1085 // as an unsigned decimal integer.
1086 // len is the length of the buffer, i.e. the number of decimal digits.
1087 int len = 0;
1088 int decimal_exponent = 0;
1089 dtoa_impl::grisu2(first, len, decimal_exponent, value);
1090
1091 assert(len <= std::numeric_limits<FloatType>::max_digits10);
1092
1093 // Format the buffer like printf("%.*g", prec, value)
1094 constexpr int kMinExp = -4;
1095 // Use digits10 here to increase compatibility with version 2.
1096 constexpr int kMaxExp = std::numeric_limits<FloatType>::digits10;
1097
1098 assert(last - first >= kMaxExp + 2);
1099 assert(last - first >= 2 + (-kMinExp - 1) + std::numeric_limits<FloatType>::max_digits10);
1100 assert(last - first >= std::numeric_limits<FloatType>::max_digits10 + 6);
1101
1102 return dtoa_impl::format_buffer(first, len, decimal_exponent, kMinExp, kMaxExp);
1103}
1104
1105} // namespace detail
1106} // namespace nlohmann
1107