34 #define copysign _copysign
37 #define UINT64_C2(h, l) ((static_cast<uint64_t>(h) << 32) | static_cast<uint64_t>(l))
42 DiyFp(uint64_t f,
int e) : f(f), e(e) {}
50 int biased_e = (u.u64 & kDpExponentMask) >> kDpSignificandSize;
51 uint64_t significand = (u.u64 & kDpSignificandMask);
53 f = significand + kDpHiddenBit;
54 e = biased_e - kDpExponentBias;
58 e = kDpMinExponent + 1;
65 return DiyFp(f - rhs.f, e);
69 #if defined(_MSC_VER) && defined(_M_AMD64)
71 uint64_t l = _umul128(f, rhs.f, &h);
72 if (l & (uint64_t(1) << 63))
74 return DiyFp(h, e + rhs.e + 64);
75 #elif (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6)) && defined(__x86_64__)
76 unsigned __int128 p =
static_cast<unsigned __int128
>(f) * static_cast<unsigned __int128>(rhs.f);
78 uint64_t l =
static_cast<uint64_t
>(p);
79 if (l & (uint64_t(1) << 63))
81 return DiyFp(h, e + rhs.e + 64);
83 const uint64_t M32 = 0xFFFFFFFF;
84 const uint64_t a = f >> 32;
85 const uint64_t b = f & M32;
86 const uint64_t c = rhs.f >> 32;
87 const uint64_t d = rhs.f & M32;
88 const uint64_t ac = a * c;
89 const uint64_t bc = b * c;
90 const uint64_t ad = a * d;
91 const uint64_t bd = b * d;
92 uint64_t tmp = (bd >> 32) + (ad & M32) + (bc & M32);
94 return DiyFp(ac + (ad >> 32) + (bc >> 32) + (tmp >> 32), e + rhs.e + 64);
98 DiyFp Normalize()
const {
99 #if defined(_MSC_VER) && defined(_M_AMD64)
101 _BitScanReverse64(&index, f);
102 return DiyFp(f << (63 - index), e - (63 - index));
103 #elif defined(__GNUC__)
104 int s = __builtin_clzll(f);
105 return DiyFp(f << s, e - s);
108 while (!(res.f & kDpHiddenBit)) {
112 res.f <<= (kDiySignificandSize - kDpSignificandSize - 1);
113 res.e = res.e - (kDiySignificandSize - kDpSignificandSize - 1);
118 DiyFp NormalizeBoundary()
const {
119 #if defined(_MSC_VER) && defined(_M_AMD64)
121 _BitScanReverse64(&index, f);
122 return DiyFp (f << (63 - index), e - (63 - index));
125 while (!(res.f & (kDpHiddenBit << 1))) {
129 res.f <<= (kDiySignificandSize - kDpSignificandSize - 2);
130 res.e = res.e - (kDiySignificandSize - kDpSignificandSize - 2);
135 void NormalizedBoundaries(
DiyFp* minus,
DiyFp* plus)
const {
136 DiyFp pl =
DiyFp((f << 1) + 1, e - 1).NormalizeBoundary();
137 DiyFp mi = (f == kDpHiddenBit) ?
DiyFp((f << 2) - 1, e - 2) :
DiyFp((f << 1) - 1, e - 1);
138 mi.f <<= mi.e - pl.e;
144 static const int kDiySignificandSize = 64;
145 static const int kDpSignificandSize = 52;
146 static const int kDpExponentBias = 0x3FF + kDpSignificandSize;
147 static const int kDpMinExponent = -kDpExponentBias;
148 static const uint64_t kDpExponentMask = UINT64_C2(0x7FF00000, 0x00000000);
149 static const uint64_t kDpSignificandMask = UINT64_C2(0x000FFFFF, 0xFFFFFFFF);
150 static const uint64_t kDpHiddenBit = UINT64_C2(0x00100000, 0x00000000);
156 inline static DiyFp GetCachedPower(
int e,
int* K) {
158 static const uint64_t kCachedPowers_F[] = {
159 UINT64_C2(0xfa8fd5a0, 0x081c0288), UINT64_C2(0xbaaee17f, 0xa23ebf76),
160 UINT64_C2(0x8b16fb20, 0x3055ac76), UINT64_C2(0xcf42894a, 0x5dce35ea),
161 UINT64_C2(0x9a6bb0aa, 0x55653b2d), UINT64_C2(0xe61acf03, 0x3d1a45df),
162 UINT64_C2(0xab70fe17, 0xc79ac6ca), UINT64_C2(0xff77b1fc, 0xbebcdc4f),
163 UINT64_C2(0xbe5691ef, 0x416bd60c), UINT64_C2(0x8dd01fad, 0x907ffc3c),
164 UINT64_C2(0xd3515c28, 0x31559a83), UINT64_C2(0x9d71ac8f, 0xada6c9b5),
165 UINT64_C2(0xea9c2277, 0x23ee8bcb), UINT64_C2(0xaecc4991, 0x4078536d),
166 UINT64_C2(0x823c1279, 0x5db6ce57), UINT64_C2(0xc2109436, 0x4dfb5637),
167 UINT64_C2(0x9096ea6f, 0x3848984f), UINT64_C2(0xd77485cb, 0x25823ac7),
168 UINT64_C2(0xa086cfcd, 0x97bf97f4), UINT64_C2(0xef340a98, 0x172aace5),
169 UINT64_C2(0xb23867fb, 0x2a35b28e), UINT64_C2(0x84c8d4df, 0xd2c63f3b),
170 UINT64_C2(0xc5dd4427, 0x1ad3cdba), UINT64_C2(0x936b9fce, 0xbb25c996),
171 UINT64_C2(0xdbac6c24, 0x7d62a584), UINT64_C2(0xa3ab6658, 0x0d5fdaf6),
172 UINT64_C2(0xf3e2f893, 0xdec3f126), UINT64_C2(0xb5b5ada8, 0xaaff80b8),
173 UINT64_C2(0x87625f05, 0x6c7c4a8b), UINT64_C2(0xc9bcff60, 0x34c13053),
174 UINT64_C2(0x964e858c, 0x91ba2655), UINT64_C2(0xdff97724, 0x70297ebd),
175 UINT64_C2(0xa6dfbd9f, 0xb8e5b88f), UINT64_C2(0xf8a95fcf, 0x88747d94),
176 UINT64_C2(0xb9447093, 0x8fa89bcf), UINT64_C2(0x8a08f0f8, 0xbf0f156b),
177 UINT64_C2(0xcdb02555, 0x653131b6), UINT64_C2(0x993fe2c6, 0xd07b7fac),
178 UINT64_C2(0xe45c10c4, 0x2a2b3b06), UINT64_C2(0xaa242499, 0x697392d3),
179 UINT64_C2(0xfd87b5f2, 0x8300ca0e), UINT64_C2(0xbce50864, 0x92111aeb),
180 UINT64_C2(0x8cbccc09, 0x6f5088cc), UINT64_C2(0xd1b71758, 0xe219652c),
181 UINT64_C2(0x9c400000, 0x00000000), UINT64_C2(0xe8d4a510, 0x00000000),
182 UINT64_C2(0xad78ebc5, 0xac620000), UINT64_C2(0x813f3978, 0xf8940984),
183 UINT64_C2(0xc097ce7b, 0xc90715b3), UINT64_C2(0x8f7e32ce, 0x7bea5c70),
184 UINT64_C2(0xd5d238a4, 0xabe98068), UINT64_C2(0x9f4f2726, 0x179a2245),
185 UINT64_C2(0xed63a231, 0xd4c4fb27), UINT64_C2(0xb0de6538, 0x8cc8ada8),
186 UINT64_C2(0x83c7088e, 0x1aab65db), UINT64_C2(0xc45d1df9, 0x42711d9a),
187 UINT64_C2(0x924d692c, 0xa61be758), UINT64_C2(0xda01ee64, 0x1a708dea),
188 UINT64_C2(0xa26da399, 0x9aef774a), UINT64_C2(0xf209787b, 0xb47d6b85),
189 UINT64_C2(0xb454e4a1, 0x79dd1877), UINT64_C2(0x865b8692, 0x5b9bc5c2),
190 UINT64_C2(0xc83553c5, 0xc8965d3d), UINT64_C2(0x952ab45c, 0xfa97a0b3),
191 UINT64_C2(0xde469fbd, 0x99a05fe3), UINT64_C2(0xa59bc234, 0xdb398c25),
192 UINT64_C2(0xf6c69a72, 0xa3989f5c), UINT64_C2(0xb7dcbf53, 0x54e9bece),
193 UINT64_C2(0x88fcf317, 0xf22241e2), UINT64_C2(0xcc20ce9b, 0xd35c78a5),
194 UINT64_C2(0x98165af3, 0x7b2153df), UINT64_C2(0xe2a0b5dc, 0x971f303a),
195 UINT64_C2(0xa8d9d153, 0x5ce3b396), UINT64_C2(0xfb9b7cd9, 0xa4a7443c),
196 UINT64_C2(0xbb764c4c, 0xa7a44410), UINT64_C2(0x8bab8eef, 0xb6409c1a),
197 UINT64_C2(0xd01fef10, 0xa657842c), UINT64_C2(0x9b10a4e5, 0xe9913129),
198 UINT64_C2(0xe7109bfb, 0xa19c0c9d), UINT64_C2(0xac2820d9, 0x623bf429),
199 UINT64_C2(0x80444b5e, 0x7aa7cf85), UINT64_C2(0xbf21e440, 0x03acdd2d),
200 UINT64_C2(0x8e679c2f, 0x5e44ff8f), UINT64_C2(0xd433179d, 0x9c8cb841),
201 UINT64_C2(0x9e19db92, 0xb4e31ba9), UINT64_C2(0xeb96bf6e, 0xbadf77d9),
202 UINT64_C2(0xaf87023b, 0x9bf0ee6b)
204 static const int16_t kCachedPowers_E[] = {
205 -1220, -1193, -1166, -1140, -1113, -1087, -1060, -1034, -1007, -980,
206 -954, -927, -901, -874, -847, -821, -794, -768, -741, -715,
207 -688, -661, -635, -608, -582, -555, -529, -502, -475, -449,
208 -422, -396, -369, -343, -316, -289, -263, -236, -210, -183,
209 -157, -130, -103, -77, -50, -24, 3, 30, 56, 83,
210 109, 136, 162, 189, 216, 242, 269, 295, 322, 348,
211 375, 402, 428, 455, 481, 508, 534, 561, 588, 614,
212 641, 667, 694, 720, 747, 774, 800, 827, 853, 880,
213 907, 933, 960, 986, 1013, 1039, 1066
217 double dk = (-61 - e) * 0.30102999566398114 + 347;
218 int k =
static_cast<int>(dk);
222 unsigned index =
static_cast<unsigned>((k >> 3) + 1);
223 *K = -(-348 +
static_cast<int>(index << 3));
225 assert(index <
sizeof(kCachedPowers_F) /
sizeof(kCachedPowers_F[0]));
226 return DiyFp(kCachedPowers_F[index], kCachedPowers_E[index]);
229 inline static void GrisuRound(
char* buffer,
int len, uint64_t delta, uint64_t rest, uint64_t ten_kappa, uint64_t wp_w) {
230 while (rest < wp_w && delta - rest >= ten_kappa &&
231 (rest + ten_kappa < wp_w ||
232 wp_w - rest > rest + ten_kappa - wp_w)) {
238 inline static unsigned CountDecimalDigit32(uint32_t n) {
240 if (n < 10)
return 1;
241 if (n < 100)
return 2;
242 if (n < 1000)
return 3;
243 if (n < 10000)
return 4;
244 if (n < 100000)
return 5;
245 if (n < 1000000)
return 6;
246 if (n < 10000000)
return 7;
247 if (n < 100000000)
return 8;
248 if (n < 1000000000)
return 9;
252 inline static void DigitGen(
const DiyFp& W,
const DiyFp& Mp, uint64_t delta,
char* buffer,
int* len,
int* K) {
253 static const uint32_t kPow10[] = { 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000 };
254 const DiyFp one(uint64_t(1) << -Mp.e, Mp.e);
255 const DiyFp wp_w = Mp - W;
256 uint32_t p1 =
static_cast<uint32_t
>(Mp.f >> -one.e);
257 uint64_t p2 = Mp.f & (one.f - 1);
258 int kappa =
static_cast<int>(CountDecimalDigit32(p1));
264 case 10: d = p1 / 1000000000; p1 %= 1000000000;
break;
265 case 9: d = p1 / 100000000; p1 %= 100000000;
break;
266 case 8: d = p1 / 10000000; p1 %= 10000000;
break;
267 case 7: d = p1 / 1000000; p1 %= 1000000;
break;
268 case 6: d = p1 / 100000; p1 %= 100000;
break;
269 case 5: d = p1 / 10000; p1 %= 10000;
break;
270 case 4: d = p1 / 1000; p1 %= 1000;
break;
271 case 3: d = p1 / 100; p1 %= 100;
break;
272 case 2: d = p1 / 10; p1 %= 10;
break;
273 case 1: d = p1; p1 = 0;
break;
275 #if defined(_MSC_VER)
277 #elif __GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 5)
278 __builtin_unreachable();
284 buffer[(*len)++] =
'0' +
static_cast<char>(d);
286 uint64_t tmp = (
static_cast<uint64_t
>(p1) << -one.e) + p2;
289 GrisuRound(buffer, *len, delta, tmp, static_cast<uint64_t>(kPow10[kappa]) << -one.e, wp_w.f);
298 char d =
static_cast<char>(p2 >> -one.e);
300 buffer[(*len)++] =
'0' + d;
305 GrisuRound(buffer, *len, delta, p2, one.f, wp_w.f * kPow10[-kappa]);
311 inline static void Grisu2(
double value,
char* buffer,
int* length,
int* K) {
312 const DiyFp v(value);
314 v.NormalizedBoundaries(&w_m, &w_p);
316 const DiyFp c_mk = GetCachedPower(w_p.e, K);
317 const DiyFp W = v.Normalize() * c_mk;
318 DiyFp Wp = w_p * c_mk;
319 DiyFp Wm = w_m * c_mk;
322 DigitGen(W, Wp, Wp.f - Wm.f, buffer, length, K);
325 static const char cDigitsLut[200] = {
326 '0',
'0',
'0',
'1',
'0',
'2',
'0',
'3',
'0',
'4',
'0',
'5',
'0',
'6',
'0',
'7',
'0',
'8',
'0',
'9',
327 '1',
'0',
'1',
'1',
'1',
'2',
'1',
'3',
'1',
'4',
'1',
'5',
'1',
'6',
'1',
'7',
'1',
'8',
'1',
'9',
328 '2',
'0',
'2',
'1',
'2',
'2',
'2',
'3',
'2',
'4',
'2',
'5',
'2',
'6',
'2',
'7',
'2',
'8',
'2',
'9',
329 '3',
'0',
'3',
'1',
'3',
'2',
'3',
'3',
'3',
'4',
'3',
'5',
'3',
'6',
'3',
'7',
'3',
'8',
'3',
'9',
330 '4',
'0',
'4',
'1',
'4',
'2',
'4',
'3',
'4',
'4',
'4',
'5',
'4',
'6',
'4',
'7',
'4',
'8',
'4',
'9',
331 '5',
'0',
'5',
'1',
'5',
'2',
'5',
'3',
'5',
'4',
'5',
'5',
'5',
'6',
'5',
'7',
'5',
'8',
'5',
'9',
332 '6',
'0',
'6',
'1',
'6',
'2',
'6',
'3',
'6',
'4',
'6',
'5',
'6',
'6',
'6',
'7',
'6',
'8',
'6',
'9',
333 '7',
'0',
'7',
'1',
'7',
'2',
'7',
'3',
'7',
'4',
'7',
'5',
'7',
'6',
'7',
'7',
'7',
'8',
'7',
'9',
334 '8',
'0',
'8',
'1',
'8',
'2',
'8',
'3',
'8',
'4',
'8',
'5',
'8',
'6',
'8',
'7',
'8',
'8',
'8',
'9',
335 '9',
'0',
'9',
'1',
'9',
'2',
'9',
'3',
'9',
'4',
'9',
'5',
'9',
'6',
'9',
'7',
'9',
'8',
'9',
'9'
338 inline void WriteExponent(
int K,
char* buffer) {
345 *buffer++ =
'0' +
static_cast<char>(K / 100);
347 const char* d = cDigitsLut + K * 2;
352 const char* d = cDigitsLut + K * 2;
357 *buffer++ =
'0' +
static_cast<char>(K);
362 inline static void Prettify(
char* buffer,
int length,
int k) {
363 const int kk = length + k;
365 if (length <= kk && kk <= 21) {
367 for (
int i = length; i < kk; i++)
370 buffer[kk + 1] =
'0';
371 buffer[kk + 2] =
'\0';
373 else if (0 < kk && kk <= 21) {
375 memmove(&buffer[kk + 1], &buffer[kk], length - kk);
377 buffer[length + 1] =
'\0';
379 else if (-6 < kk && kk <= 0) {
381 const int offset = 2 - kk;
382 memmove(&buffer[offset], &buffer[0], length);
385 for (
int i = 2; i < offset; i++)
387 buffer[length + offset] =
'\0';
389 else if (length == 1) {
392 WriteExponent(kk - 1, &buffer[2]);
396 memmove(&buffer[2], &buffer[1], length - 1);
398 buffer[length + 1] =
'e';
399 WriteExponent(kk - 1, &buffer[0 + length + 2]);
403 void pdtoa(
double value,
char *buffer) {
405 if (copysign(1.0, value) < 0) {
407 if (signbit(value)) {
417 }
else if (cnan(value)) {
422 }
else if (value == 0.0) {
427 }
else if (value == 1.0) {
434 Grisu2(value, buffer, &length, &K);
435 Prettify(buffer, length, K);