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