7 #define fracmask 0x000FFFFFFFFFFFFFU
8 #define expmask 0x7FF0000000000000U
9 #define hiddenbit 0x0010000000000000U
10 #define signmask 0x8000000000000000U
11 #define expbias (1023 + 52)
13 #define absv(n) ((n) < 0 ? -(n) : (n))
14 #define minv(a, b) ((a) < (b) ? (a) : (b))
16 static uint64_t tens[] = {
17 10000000000000000000U, 1000000000000000000U, 100000000000000000U,
18 10000000000000000U, 1000000000000000U, 100000000000000U,
19 10000000000000U, 1000000000000U, 100000000000U,
20 10000000000U, 1000000000U, 100000000U,
21 10000000U, 1000000U, 100000U,
26 static inline uint64_t get_dbits(double d)
36 static Fp build_fp(double d)
38 uint64_t bits = get_dbits(d);
41 fp.frac = bits & fracmask;
42 fp.exp = (bits & expmask) >> 52;
49 fp.exp = -expbias + 1;
55 static void normalize(Fp* fp)
57 while ((fp->frac & hiddenbit) == 0) {
62 int shift = 64 - 52 - 1;
67 static void get_normalized_boundaries(Fp* fp, Fp* lower, Fp* upper)
69 upper->frac = (fp->frac << 1) + 1;
70 upper->exp = fp->exp - 1;
72 while ((upper->frac & (hiddenbit << 1)) == 0) {
77 int u_shift = 64 - 52 - 2;
79 upper->frac <<= u_shift;
80 upper->exp = upper->exp - u_shift;
83 int l_shift = fp->frac == hiddenbit ? 2 : 1;
85 lower->frac = (fp->frac << l_shift) - 1;
86 lower->exp = fp->exp - l_shift;
89 lower->frac <<= lower->exp - upper->exp;
90 lower->exp = upper->exp;
93 static Fp multiply(Fp* a, Fp* b)
95 const uint64_t lomask = 0x00000000FFFFFFFF;
97 uint64_t ah_bl = (a->frac >> 32) * (b->frac & lomask);
98 uint64_t al_bh = (a->frac & lomask) * (b->frac >> 32);
99 uint64_t al_bl = (a->frac & lomask) * (b->frac & lomask);
100 uint64_t ah_bh = (a->frac >> 32) * (b->frac >> 32);
102 uint64_t tmp = (ah_bl & lomask) + (al_bh & lomask) + (al_bl >> 32);
107 ah_bh + (ah_bl >> 32) + (al_bh >> 32) + (tmp >> 32),
114 static void round_digit(char* digits, int ndigits, uint64_t delta, uint64_t rem, uint64_t kappa, uint64_t frac)
116 while (rem < frac && delta - rem >= kappa &&
117 (rem + kappa < frac || frac - rem > rem + kappa - frac)) {
119 digits[ndigits - 1]--;
124 static int generate_digits(Fp* fp, Fp* upper, Fp* lower, char* digits, int* K)
126 uint64_t wfrac = upper->frac - fp->frac;
127 uint64_t delta = upper->frac - lower->frac;
130 one.frac = 1ULL << -upper->exp;
131 one.exp = upper->exp;
133 uint64_t part1 = upper->frac >> -one.exp;
134 uint64_t part2 = upper->frac & (one.frac - 1);
136 int idx = 0, kappa = 10;
139 for(divp = tens + 10; kappa > 0; divp++) {
141 uint64_t div = *divp;
142 unsigned digit = part1 / div;
145 digits[idx++] = digit + '0';
148 part1 -= digit * div;
151 uint64_t tmp = (part1 <<-one.exp) + part2;
154 round_digit(digits, idx, delta, tmp, div << -one.exp, wfrac);
161 uint64_t* unit = tens + 18;
168 unsigned digit = part2 >> -one.exp;
170 digits[idx++] = digit + '0';
173 part2 &= one.frac - 1;
176 round_digit(digits, idx, delta, part2, one.frac, wfrac * *unit);
185 static int grisu2(double d, char* digits, int* K)
190 get_normalized_boundaries(&w, &lower, &upper);
195 Fp cp = find_cachedpow10(upper.exp, &k);
197 w = multiply(&w, &cp);
198 upper = multiply(&upper, &cp);
199 lower = multiply(&lower, &cp);
206 return generate_digits(&w, &upper, &lower, digits, K);
209 static int emit_digits(char* digits, int ndigits, char* dest, int K, bool neg)
211 int exp = absv(K + ndigits - 1);
213 /* write plain integer */
214 if(K >= 0 && (exp < (ndigits + 7))) {
215 memcpy(dest, digits, ndigits);
216 memset(dest + ndigits, '0', K);
221 /* write decimal w/o scientific notation */
222 if(K < 0 && (K > -7 || exp < 4)) {
223 int offset = ndigits - absv(K);
224 /* fp < 1.0 -> write leading zero */
229 memset(dest + 2, '0', offset);
230 memcpy(dest + offset + 2, digits, ndigits);
232 return ndigits + 2 + offset;
236 memcpy(dest, digits, offset);
238 memcpy(dest + offset + 1, digits + offset, ndigits - offset);
244 /* write decimal w/ scientific notation */
245 ndigits = minv(ndigits, 18 - neg);
248 dest[idx++] = digits[0];
252 memcpy(dest + idx, digits + 1, ndigits - 1);
258 char sign = K + ndigits - 1 < 0 ? '-' : '+';
265 dest[idx++] = cent + '0';
270 dest[idx++] = dec + '0';
277 dest[idx++] = exp % 10 + '0';
282 static int filter_special(double fp, char* dest)
289 uint64_t bits = get_dbits(fp);
291 bool nan = (bits & expmask) == expmask;
297 if(bits & fracmask) {
298 dest[0] = 'n'; dest[1] = 'a'; dest[2] = 'n';
301 dest[0] = 'i'; dest[1] = 'n'; dest[2] = 'f';
307 int fpconv_dtoa(double d, char dest[24])
314 if(get_dbits(d) & signmask) {
320 int spec = filter_special(d, dest + str_len);
323 return str_len + spec;
327 int ndigits = grisu2(d, digits, &K);
329 str_len += emit_digits(digits, ndigits, dest + str_len, K, neg);