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