1 /* This file is part of an implementation of the "grisu3" double to string
2         conversion algorithm described in the research paper
3
4         "Printing Floating-Point Numbers Quickly And Accurately with Integers"
5         by Florian Loitsch, available at
6         http://www.cs.tufts.edu/~nr/cs257/archive/florian-loitsch/printf.pdf */
7
8 #include <stdint.h> // uint64_t etc.
9 #include <assert.h> // assert
10 #include <math.h> // ceil
11 #include <stdio.h> // sprintf
12
13 #ifdef _MSC_VER
14 #pragma warning(disable : 4204) // nonstandard extension used : non-constant aggregate initializer
15 #endif
16
17 #define D64_SIGN         0x8000000000000000ULL
18 #define D64_EXP_MASK     0x7FF0000000000000ULL
19 #define D64_FRACT_MASK   0x000FFFFFFFFFFFFFULL
20 #define D64_IMPLICIT_ONE 0x0010000000000000ULL
21 #define D64_EXP_POS      52
22 #define D64_EXP_BIAS     1075
23 #define DIYFP_FRACT_SIZE 64
24 #define D_1_LOG2_10      0.30102999566398114 // 1 / lg(10)
25 #define MIN_TARGET_EXP   -60
26 #define MASK32           0xFFFFFFFFULL
27
28 #define CAST_U64(d) (*(uint64_t*)&d)
29 #define MIN(x,y) ((x) <= (y) ? (x) : (y))
30 #define MAX(x,y) ((x) >= (y) ? (x) : (y))
31
32 #define MIN_CACHED_EXP -348
33 #define CACHED_EXP_STEP 8
34
35 typedef struct diy_fp
36 {
37         uint64_t f;
38         int e;
39 } diy_fp;
40
41 typedef struct power
42 {
43         uint64_t fract;
44         int16_t b_exp, d_exp;
45 } power;
46
47 static const power pow_cache[] =
48 {
49         { 0xfa8fd5a0081c0288ULL, -1220, -348 },
50         { 0xbaaee17fa23ebf76ULL, -1193, -340 },
51         { 0x8b16fb203055ac76ULL, -1166, -332 },
52         { 0xcf42894a5dce35eaULL, -1140, -324 },
53         { 0x9a6bb0aa55653b2dULL, -1113, -316 },
54         { 0xe61acf033d1a45dfULL, -1087, -308 },
55         { 0xab70fe17c79ac6caULL, -1060, -300 },
56         { 0xff77b1fcbebcdc4fULL, -1034, -292 },
57         { 0xbe5691ef416bd60cULL, -1007, -284 },
58         { 0x8dd01fad907ffc3cULL,  -980, -276 },
59         { 0xd3515c2831559a83ULL,  -954, -268 },
60         { 0x9d71ac8fada6c9b5ULL,  -927, -260 },
61         { 0xea9c227723ee8bcbULL,  -901, -252 },
62         { 0xaecc49914078536dULL,  -874, -244 },
63         { 0x823c12795db6ce57ULL,  -847, -236 },
64         { 0xc21094364dfb5637ULL,  -821, -228 },
65         { 0x9096ea6f3848984fULL,  -794, -220 },
66         { 0xd77485cb25823ac7ULL,  -768, -212 },
67         { 0xa086cfcd97bf97f4ULL,  -741, -204 },
68         { 0xef340a98172aace5ULL,  -715, -196 },
69         { 0xb23867fb2a35b28eULL,  -688, -188 },
70         { 0x84c8d4dfd2c63f3bULL,  -661, -180 },
71         { 0xc5dd44271ad3cdbaULL,  -635, -172 },
72         { 0x936b9fcebb25c996ULL,  -608, -164 },
73         { 0xdbac6c247d62a584ULL,  -582, -156 },
74         { 0xa3ab66580d5fdaf6ULL,  -555, -148 },
75         { 0xf3e2f893dec3f126ULL,  -529, -140 },
76         { 0xb5b5ada8aaff80b8ULL,  -502, -132 },
77         { 0x87625f056c7c4a8bULL,  -475, -124 },
78         { 0xc9bcff6034c13053ULL,  -449, -116 },
79         { 0x964e858c91ba2655ULL,  -422, -108 },
80         { 0xdff9772470297ebdULL,  -396, -100 },
81         { 0xa6dfbd9fb8e5b88fULL,  -369,  -92 },
82         { 0xf8a95fcf88747d94ULL,  -343,  -84 },
83         { 0xb94470938fa89bcfULL,  -316,  -76 },
84         { 0x8a08f0f8bf0f156bULL,  -289,  -68 },
85         { 0xcdb02555653131b6ULL,  -263,  -60 },
86         { 0x993fe2c6d07b7facULL,  -236,  -52 },
87         { 0xe45c10c42a2b3b06ULL,  -210,  -44 },
88         { 0xaa242499697392d3ULL,  -183,  -36 },
89         { 0xfd87b5f28300ca0eULL,  -157,  -28 },
90         { 0xbce5086492111aebULL,  -130,  -20 },
91         { 0x8cbccc096f5088ccULL,  -103,  -12 },
92         { 0xd1b71758e219652cULL,   -77,   -4 },
93         { 0x9c40000000000000ULL,   -50,    4 },
94         { 0xe8d4a51000000000ULL,   -24,   12 },
95         { 0xad78ebc5ac620000ULL,     3,   20 },
96         { 0x813f3978f8940984ULL,    30,   28 },
97         { 0xc097ce7bc90715b3ULL,    56,   36 },
98         { 0x8f7e32ce7bea5c70ULL,    83,   44 },
99         { 0xd5d238a4abe98068ULL,   109,   52 },
100         { 0x9f4f2726179a2245ULL,   136,   60 },
101         { 0xed63a231d4c4fb27ULL,   162,   68 },
102         { 0xb0de65388cc8ada8ULL,   189,   76 },
103         { 0x83c7088e1aab65dbULL,   216,   84 },
104         { 0xc45d1df942711d9aULL,   242,   92 },
105         { 0x924d692ca61be758ULL,   269,  100 },
106         { 0xda01ee641a708deaULL,   295,  108 },
107         { 0xa26da3999aef774aULL,   322,  116 },
108         { 0xf209787bb47d6b85ULL,   348,  124 },
109         { 0xb454e4a179dd1877ULL,   375,  132 },
110         { 0x865b86925b9bc5c2ULL,   402,  140 },
111         { 0xc83553c5c8965d3dULL,   428,  148 },
112         { 0x952ab45cfa97a0b3ULL,   455,  156 },
113         { 0xde469fbd99a05fe3ULL,   481,  164 },
114         { 0xa59bc234db398c25ULL,   508,  172 },
115         { 0xf6c69a72a3989f5cULL,   534,  180 },
116         { 0xb7dcbf5354e9beceULL,   561,  188 },
117         { 0x88fcf317f22241e2ULL,   588,  196 },
118         { 0xcc20ce9bd35c78a5ULL,   614,  204 },
119         { 0x98165af37b2153dfULL,   641,  212 },
120         { 0xe2a0b5dc971f303aULL,   667,  220 },
121         { 0xa8d9d1535ce3b396ULL,   694,  228 },
122         { 0xfb9b7cd9a4a7443cULL,   720,  236 },
123         { 0xbb764c4ca7a44410ULL,   747,  244 },
124         { 0x8bab8eefb6409c1aULL,   774,  252 },
125         { 0xd01fef10a657842cULL,   800,  260 },
126         { 0x9b10a4e5e9913129ULL,   827,  268 },
127         { 0xe7109bfba19c0c9dULL,   853,  276 },
128         { 0xac2820d9623bf429ULL,   880,  284 },
129         { 0x80444b5e7aa7cf85ULL,   907,  292 },
130         { 0xbf21e44003acdd2dULL,   933,  300 },
131         { 0x8e679c2f5e44ff8fULL,   960,  308 },
132         { 0xd433179d9c8cb841ULL,   986,  316 },
133         { 0x9e19db92b4e31ba9ULL,  1013,  324 },
134         { 0xeb96bf6ebadf77d9ULL,  1039,  332 },
135         { 0xaf87023b9bf0ee6bULL,  1066,  340 }
136 };
137
138 static int cached_pow(int exp, diy_fp *p)
139 {
140         int k = (int)ceil((exp+DIYFP_FRACT_SIZE-1) * D_1_LOG2_10);
141         int i = (k-MIN_CACHED_EXP-1) / CACHED_EXP_STEP + 1;
142         p->f = pow_cache[i].fract;
143         p->e = pow_cache[i].b_exp;
144         return pow_cache[i].d_exp;
145 }
146
147 static diy_fp minus(diy_fp x, diy_fp y)
148 {
149         diy_fp d; d.f = x.f - y.f; d.e = x.e;
150         assert(x.e == y.e && x.f >= y.f);
151         return d;
152 }
153
154 static diy_fp multiply(diy_fp x, diy_fp y)
155 {
156         uint64_t a, b, c, d, ac, bc, ad, bd, tmp;
157         diy_fp r;
158         a = x.f >> 32; b = x.f & MASK32;
159         c = y.f >> 32; d = y.f & MASK32;
160         ac = a*c; bc = b*c;
161         ad = a*d; bd = b*d;
162         tmp = (bd >> 32) + (ad & MASK32) + (bc & MASK32);
163         tmp += 1U << 31; // round
164         r.f = ac + (ad >> 32) + (bc >> 32) + (tmp >> 32);
165         r.e = x.e + y.e + 64;
166         return r;
167 }
168
169 static diy_fp normalize_diy_fp(diy_fp n)
170 {
171         assert(n.f != 0);
172         while(!(n.f & 0xFFC0000000000000ULL)) { n.f <<= 10; n.e -= 10; }
173         while(!(n.f & D64_SIGN)) { n.f <<= 1; --n.e; }
174         return n;
175 }
176
177 static diy_fp double2diy_fp(double d)
178 {
179         diy_fp fp;
180         uint64_t u64 = CAST_U64(d);
181         if (!(u64 & D64_EXP_MASK)) { fp.f = u64 & D64_FRACT_MASK; fp.e = 1 - D64_EXP_BIAS; }
182         else { fp.f = (u64 & D64_FRACT_MASK) + D64_IMPLICIT_ONE; fp.e = (int)((u64 & D64_EXP_MASK) >> D64_EXP_POS) - D64_EXP_BIAS; }
183         return fp;
184 }
185
186 // pow10_cache[i] = 10^(i-1)
187 static const unsigned int pow10_cache[] = { 0, 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000 };
188
189 static int largest_pow10(uint32_t n, int n_bits, uint32_t *power)
190 {
191         int guess = ((n_bits + 1) * 1233 >> 12) + 1/*skip first entry*/;
192         if (n < pow10_cache[guess]) --guess; // We don't have any guarantees that 2^n_bits <= n.
193         *power = pow10_cache[guess];
194         return guess;
195 }
196
197 static int round_weed(char *buffer, int len, uint64_t wp_W, uint64_t delta, uint64_t rest, uint64_t ten_kappa, uint64_t ulp)
198 {
199         uint64_t wp_Wup = wp_W - ulp;
200         uint64_t wp_Wdown = wp_W + ulp;
201         while(rest < wp_Wup && delta - rest >= ten_kappa
202                 && (rest + ten_kappa < wp_Wup || wp_Wup - rest >= rest + ten_kappa - wp_Wup))
203         {
204                 --buffer[len-1];
205                 rest += ten_kappa;
206         }
207         if (rest < wp_Wdown && delta - rest >= ten_kappa
208                 && (rest + ten_kappa < wp_Wdown || wp_Wdown - rest > rest + ten_kappa - wp_Wdown))
209                 return 0;
210
211         return 2*ulp <= rest && rest <= delta - 4*ulp;
212 }
213
214 static int digit_gen(diy_fp low, diy_fp w, diy_fp high, char *buffer, int *length, int *kappa)
215 {
216         uint64_t unit = 1;
217         diy_fp too_low = { low.f - unit, low.e };
218         diy_fp too_high = { high.f + unit, high.e };
219         diy_fp unsafe_interval = minus(too_high, too_low);
220         diy_fp one = { 1ULL << -w.e, w.e };
221         uint32_t p1 = (uint32_t)(too_high.f >> -one.e);
222         uint64_t p2 = too_high.f & (one.f - 1);
223         uint32_t div;
224         *kappa = largest_pow10(p1, DIYFP_FRACT_SIZE + one.e, &div);
225         *length = 0;
226
227         while(*kappa > 0)
228         {
229                 uint64_t rest;
230                 int digit = p1 / div;
231                 buffer[*length] = (char)('0' + digit);
232                 ++*length;
233                 p1 %= div;
234                 --*kappa;
235                 rest = ((uint64_t)p1 << -one.e) + p2;
236                 if (rest < unsafe_interval.f) return round_weed(buffer, *length, minus(too_high, w).f, unsafe_interval.f, rest, (uint64_t)div << -one.e, unit);
237                 div /= 10;
238         }
239
240         for(;;)
241         {
242                 int digit;
243                 p2 *= 10;
244                 unit *= 10;
245                 unsafe_interval.f *= 10;
246                 // Integer division by one.
247                 digit = (int)(p2 >> -one.e);
248                 buffer[*length] = (char)('0' + digit);
249                 ++*length;
250                 p2 &= one.f - 1;  // Modulo by one.
251                 --*kappa;
252                 if (p2 < unsafe_interval.f) return round_weed(buffer, *length, minus(too_high, w).f * unit, unsafe_interval.f, p2, one.f, unit);
253         }
254 }
255
256 static int grisu3(double v, char *buffer, int *length, int *d_exp)
257 {
258         int mk, kappa, success;
259         diy_fp dfp = double2diy_fp(v);
260         diy_fp w = normalize_diy_fp(dfp);
261
262         // normalize boundaries
263         diy_fp t = { (dfp.f << 1) + 1, dfp.e - 1 };
264         diy_fp b_plus = normalize_diy_fp(t);
265         diy_fp b_minus;
266         diy_fp c_mk; // Cached power of ten: 10^-k
267         uint64_t u64 = CAST_U64(v);
268         assert(v > 0 && v <= 1.7976931348623157e308); // Grisu only handles strictly positive finite numbers.
269         if (!(u64 & D64_FRACT_MASK) && (u64 & D64_EXP_MASK) != 0) { b_minus.f = (dfp.f << 2) - 1; b_minus.e =  dfp.e - 2;} // lower boundary is closer?
270         else { b_minus.f = (dfp.f << 1) - 1; b_minus.e = dfp.e - 1; }
271         b_minus.f = b_minus.f << (b_minus.e - b_plus.e);
272         b_minus.e = b_plus.e;
273
274         mk = cached_pow(MIN_TARGET_EXP - DIYFP_FRACT_SIZE - w.e, &c_mk);
275
276         w = multiply(w, c_mk);
277         b_minus = multiply(b_minus, c_mk);
278         b_plus  = multiply(b_plus,  c_mk);
279
280         success = digit_gen(b_minus, w, b_plus, buffer, length, &kappa);
281         *d_exp = kappa - mk;
282         return success;
283 }
284
285 // Returns length of u when written out as a string, u must be in the range of [-9999, 9999].
286 static int exp_len(int u)
287 {
288         if (u > 0) return u >= 1000 ? 4 : (u >= 100 ? 3 : (u >= 10 ? 2 : 1));
289         else if (u < 0) return u <= -1000 ? 5 : (u <= -100 ? 4 : (u <= -10 ? 3 : 2));
290         else return 1;
291 }
292
293 static int i_to_str(int val, char *str)
294 {
295         int len, i;
296         char *s;
297         char *begin = str;
298         if (val < 0) { *str++ = '-'; val = -val; }
299         s = str;
300
301         for(;;)
302         {
303                 int ni = val / 10;
304                 int digit = val - ni*10;
305                 *s++ = (char)('0' + digit);
306                 if (ni == 0)
307                         break;
308                 val = ni;
309         }
310         *s = '\0';
311         len = (int)(s - str);
312         for(i = 0; i < len/2; ++i)
313         {
314                 char ch = str[i];
315                 str[i] = str[len-1-i];
316                 str[len-1-i] = ch;
317         }
318
319         return (int)(s - begin);
320 }
321
322 int dtoa_grisu3(double v, char *dst)
323 {
324         int d_exp, len, success, decimals, i;
325         uint64_t u64 = CAST_U64(v);
326         char *s2 = dst;
327         assert(dst);
328
329         // Prehandle NaNs
330         if ((u64 << 1) > 0xFFE0000000000000ULL) return sprintf(dst, "NaN(%08X%08X)", (uint32_t)(u64 >> 32), (uint32_t)u64);
331         // Prehandle negative values.
332         if ((u64 & D64_SIGN) != 0) { *s2++ = '-'; v = -v; u64 ^= D64_SIGN; }
333         // Prehandle zero.
334         if (!u64) { *s2++ = '0'; *s2 = '\0'; return (int)(s2 - dst); }
335         // Prehandle infinity.
336         if (u64 == D64_EXP_MASK) { *s2++ = 'i'; *s2++ = 'n'; *s2++ = 'f'; *s2 = '\0'; return (int)(s2 - dst); }
337
338         success = grisu3(v, s2, &len, &d_exp);
339         // If grisu3 was not able to convert the number to a string, then use old sprintf (suboptimal).
340         if (!success) return sprintf(s2, "%.17g", v) + (int)(s2 - dst);
341
342         decimals = MIN(-d_exp, MAX(1, len-1));
343         if (d_exp < 0 && (len >= -d_exp || exp_len(d_exp+decimals)+1 <= exp_len(d_exp))) // Add decimal point?
344         {
345                 for(i = 0; i < decimals; ++i) s2[len-i] = s2[len-i-1];
346                 s2[len++ - decimals] = '.';
347                 d_exp += decimals;
348                 // Need scientific notation as well?
349                 if (d_exp != 0) { s2[len++] = 'e'; len += i_to_str(d_exp, s2+len); }
350         }
351         else if (d_exp < 0 && d_exp >= -3) // Add decimal point for numbers of form 0.000x where it's shorter?
352         {
353                 for(i = 0; i < len; ++i) s2[len-d_exp-1-i] = s2[len-i-1];
354                 s2[0] = '.';
355                 for(i = 1; i < -d_exp; ++i) s2[i] = '0';
356                 len += -d_exp;
357         }
358         // Add scientific notation?
359         else if (d_exp < 0 || d_exp > 2) { s2[len++] = 'e'; len += i_to_str(d_exp, s2+len); }
360         // Add zeroes instead of scientific notation?
361         else if (d_exp > 0) { while(d_exp-- > 0) s2[len++] = '0'; }
362         s2[len] = '\0'; // grisu3 doesn't null terminate, so ensure termination.
363         return (int)(s2+len-dst);
364 }