1 #include <stdbool.h>
2 #include <string.h>
3 #include <sys/param.h>
4
5 #include "fpconv.h"
6 #include "powers.h"
7
8 #define fracmask 0x000FFFFFFFFFFFFFU
9 #define expmask 0x7FF0000000000000U
10 #define hiddenbit 0x0010000000000000U
11 #define signmask 0x8000000000000000U
12 #define expbias (1023 + 52)
13
14 #define absv(n) ((n) < 0 ? -(n) : (n))
15 #define minv(a, b) ((a) < (b) ? (a) : (b))
16
17 static uint64_t tens[] = {
18 10000000000000000000U, 1000000000000000000U, 100000000000000000U,
19 10000000000000000U, 1000000000000000U, 100000000000000U,
20 10000000000000U, 1000000000000U, 100000000000U,
21 10000000000U, 1000000000U, 100000000U,
22 10000000U, 1000000U, 100000U,
23 10000U, 1000U, 100U,
24 10U, 1U
25 };
26
get_dbits(double d)27 static inline uint64_t get_dbits (double d) {
28 union {
29 double dbl;
30 uint64_t i;
31 } dbl_bits = {d};
32
33 return dbl_bits.i;
34 }
35
build_fp(double d)36 static Fp build_fp (double d) {
37 uint64_t bits = get_dbits (d);
38
39 Fp fp;
40 fp.frac = bits & fracmask;
41 fp.exp = (bits & expmask) >> 52u;
42
43 if (fp.exp) {
44 fp.frac += hiddenbit;
45 fp.exp -= expbias;
46
47 }
48 else {
49 fp.exp = -expbias + 1;
50 }
51
52 return fp;
53 }
54
normalize(Fp * fp)55 static void normalize (Fp *fp) {
56 while ((fp->frac & hiddenbit) == 0) {
57 fp->frac <<= 1u;
58 fp->exp--;
59 }
60
61 const unsigned int shift = 64 - 52 - 1;
62 fp->frac <<= shift;
63 fp->exp -= shift;
64 }
65
get_normalized_boundaries(Fp * fp,Fp * lower,Fp * upper)66 static void get_normalized_boundaries (Fp *fp, Fp *lower, Fp *upper) {
67 upper->frac = (fp->frac << 1u) + 1u;
68 upper->exp = fp->exp - 1u;
69
70 while ((upper->frac & (hiddenbit << 1u)) == 0) {
71 upper->frac <<= 1u;
72 upper->exp--;
73 }
74
75 const unsigned int u_shift = 64 - 52 - 2;
76
77 upper->frac <<= u_shift;
78 upper->exp = upper->exp - u_shift;
79
80
81 unsigned int l_shift = fp->frac == hiddenbit ? 2u : 1u;
82
83 lower->frac = (fp->frac << l_shift) - 1;
84 lower->exp = fp->exp - l_shift;
85
86
87 lower->frac <<= lower->exp - upper->exp;
88 lower->exp = upper->exp;
89 }
90
multiply(Fp * a,Fp * b)91 static Fp multiply (Fp *a, Fp *b) {
92 const uint64_t lomask = 0x00000000FFFFFFFFu;
93
94 uint64_t ah_bl = (a->frac >> 32u) * (b->frac & lomask);
95 uint64_t al_bh = (a->frac & lomask) * (b->frac >> 32u);
96 uint64_t al_bl = (a->frac & lomask) * (b->frac & lomask);
97 uint64_t ah_bh = (a->frac >> 32u) * (b->frac >> 32u);
98
99 uint64_t tmp = (ah_bl & lomask) + (al_bh & lomask) + (al_bl >> 32u);
100 /* round up */
101 tmp += 1U << 31u;
102
103 Fp fp = {
104 ah_bh + (ah_bl >> 32u) + (al_bh >> 32u) + (tmp >> 32u),
105 a->exp + b->exp + 64u
106 };
107
108 return fp;
109 }
110
round_digit(char * digits,int ndigits,uint64_t delta,uint64_t rem,uint64_t kappa,uint64_t frac)111 static void round_digit (char *digits, int ndigits, uint64_t delta, uint64_t rem, uint64_t kappa, uint64_t frac) {
112 while (rem < frac && delta - rem >= kappa &&
113 (rem + kappa < frac || frac - rem > rem + kappa - frac)) {
114
115 digits[ndigits - 1]--;
116 rem += kappa;
117 }
118 }
119
generate_digits(Fp * fp,Fp * upper,Fp * lower,char * digits,int * K)120 static int generate_digits (Fp *fp, Fp *upper, Fp *lower, char *digits, int *K) {
121 uint64_t wfrac = upper->frac - fp->frac;
122 uint64_t delta = upper->frac - lower->frac;
123
124 Fp one;
125 one.frac = 1ULL << -upper->exp;
126 one.exp = upper->exp;
127
128 uint64_t part1 = upper->frac >> -one.exp;
129 uint64_t part2 = upper->frac & (one.frac - 1);
130
131 int idx = 0, kappa = 10;
132 uint64_t *divp;
133 /* 1000000000 */
134 for (divp = tens + 10; kappa > 0; divp++) {
135
136 uint64_t div = *divp;
137 unsigned digit = part1 / div;
138
139 if (digit || idx) {
140 digits[idx++] = digit + '0';
141 }
142
143 part1 -= digit * div;
144 kappa--;
145
146 uint64_t tmp = (part1 << -one.exp) + part2;
147 if (tmp <= delta) {
148 *K += kappa;
149 round_digit (digits, idx, delta, tmp, div << -one.exp, wfrac);
150
151 return idx;
152 }
153 }
154
155 /* 10 */
156 uint64_t *unit = tens + 18;
157
158 while (true) {
159 part2 *= 10;
160 delta *= 10;
161 kappa--;
162
163 unsigned digit = part2 >> -one.exp;
164 if (digit || idx) {
165 digits[idx++] = digit + '0';
166 }
167
168 part2 &= one.frac - 1;
169 if (part2 < delta) {
170 *K += kappa;
171 round_digit (digits, idx, delta, part2, one.frac, wfrac * *unit);
172
173 return idx;
174 }
175
176 unit--;
177 }
178 }
179
grisu2(double d,char * digits,int * K)180 static int grisu2 (double d, char *digits, int *K) {
181 Fp w = build_fp (d);
182
183 Fp lower, upper;
184 get_normalized_boundaries (&w, &lower, &upper);
185
186 normalize (&w);
187
188 int k;
189 Fp cp = find_cachedpow10 (upper.exp, &k);
190
191 w = multiply (&w, &cp);
192 upper = multiply (&upper, &cp);
193 lower = multiply (&lower, &cp);
194
195 lower.frac++;
196 upper.frac--;
197
198 *K = -k;
199
200 return generate_digits (&w, &upper, &lower, digits, K);
201 }
202
emit_integer(char * digits,int ndigits,char * dest,int K,bool neg,unsigned precision)203 static inline int emit_integer (char *digits, int ndigits,
204 char *dest, int K, bool neg,
205 unsigned precision)
206 {
207 char *d = dest;
208
209 memcpy (d, digits, ndigits);
210 d += ndigits;
211 memset (d, '0', K);
212 d += K;
213
214 precision = MIN(precision, FPCONV_BUFLEN - (ndigits + K + 1));
215
216 if (precision) {
217 *d++ = '.';
218 memset (d, '0', precision);
219 d += precision;
220 }
221
222 return d - dest;
223 }
224
emit_scientific_digits(char * digits,int ndigits,char * dest,int K,bool neg,unsigned precision,int exp)225 static inline int emit_scientific_digits (char *digits, int ndigits,
226 char *dest, int K, bool neg,
227 unsigned precision, int exp)
228 {
229 /* write decimal w/ scientific notation */
230 ndigits = minv(ndigits, 18 - neg);
231
232 int idx = 0;
233 dest[idx++] = digits[0];
234
235 if (ndigits > 1) {
236 dest[idx++] = '.';
237 memcpy(dest + idx, digits + 1, ndigits - 1);
238 idx += ndigits - 1;
239 }
240
241 dest[idx++] = 'e';
242
243 char sign = K + ndigits - 1 < 0 ? '-' : '+';
244 dest[idx++] = sign;
245
246 int cent = 0;
247
248 if (exp > 99) {
249 cent = exp / 100;
250 dest[idx++] = cent + '0';
251 exp -= cent * 100;
252 }
253 if (exp > 9) {
254 int dec = exp / 10;
255 dest[idx++] = dec + '0';
256 exp -= dec * 10;
257
258 }
259 else if (cent) {
260 dest[idx++] = '0';
261 }
262
263 dest[idx++] = exp % 10 + '0';
264
265 return idx;
266 }
267
emit_fixed_digits(char * digits,int ndigits,char * dest,int K,bool neg,unsigned precision,int exp)268 static inline int emit_fixed_digits (char *digits, int ndigits,
269 char *dest, int K, bool neg,
270 unsigned precision, int exp)
271 {
272 int offset = ndigits - absv(K), to_print;
273 /* fp < 1.0 -> write leading zero */
274 if (K < 0) {
275 if (offset <= 0) {
276 if (precision) {
277 if (-offset >= precision) {
278 /* Just print 0.[0]{precision} */
279 dest[0] = '0';
280 dest[1] = '.';
281 memset(dest + 2, '0', precision);
282
283 return precision + 2;
284 }
285
286 to_print = MAX(ndigits - offset, precision);
287 }
288 else {
289 to_print = ndigits - offset;
290 }
291
292 if (to_print <= FPCONV_BUFLEN - 3) {
293 offset = -offset;
294 dest[0] = '0';
295 dest[1] = '.';
296 memset(dest + 2, '0', offset);
297
298 if (precision) {
299 /* The case where offset > precision is covered previously */
300 precision -= offset;
301
302 if (precision <= ndigits) {
303 /* Truncate or leave as is */
304 memcpy(dest + offset + 2, digits, precision);
305
306 return precision + 2 + offset;
307 }
308 else {
309 /* Expand */
310 memcpy(dest + offset + 2, digits, ndigits);
311 precision -= ndigits;
312 memset(dest + offset + 2 + ndigits, '0', precision);
313
314 return ndigits + 2 + offset + precision;
315 }
316 }
317 else {
318 memcpy(dest + offset + 2, digits, ndigits);
319 }
320
321 return ndigits + 2 + offset;
322 }
323 else {
324 return emit_scientific_digits (digits, ndigits, dest, K, neg, precision, exp);
325 }
326 }
327 else {
328 /*
329 * fp > 1.0, if offset > 0 then we have less digits than
330 * fp exponent, so we need to switch to scientific notation to
331 * display number at least more or less precisely
332 */
333 if (offset > 0 && ndigits <= FPCONV_BUFLEN - 3) {
334 char *d = dest;
335 memcpy(d, digits, offset);
336 d += offset;
337 *d++ = '.';
338
339 ndigits -= offset;
340
341 if (precision) {
342 if (ndigits >= precision) {
343 /* Truncate or leave as is */
344 memcpy(d, digits + offset, precision);
345 d += precision;
346 }
347 else {
348 /* Expand */
349 memcpy(d, digits + offset, ndigits);
350 precision -= ndigits;
351 d += ndigits;
352
353 /* Check if we have enough bufspace */
354 if ((d - dest) + precision <= FPCONV_BUFLEN) {
355 memset (d, '0', precision);
356 d += precision;
357 }
358 else {
359 memset (d, '0', FPCONV_BUFLEN - (d - dest));
360 d += FPCONV_BUFLEN - (d - dest);
361 }
362 }
363 }
364 else {
365 memcpy(d, digits + offset, ndigits);
366 d += ndigits;
367 }
368
369 return d - dest;
370 }
371 }
372 }
373
374 return emit_scientific_digits (digits, ndigits, dest, K, neg, precision, exp);
375 }
376
emit_digits(char * digits,int ndigits,char * dest,int K,bool neg,unsigned precision,bool scientific)377 static int emit_digits (char *digits, int ndigits, char *dest, int K, bool neg,
378 unsigned precision, bool scientific)
379 {
380 int exp = absv(K + ndigits - 1);
381
382 /* write plain integer */
383 if (K >= 0 && (exp < (ndigits + 7))) {
384 return emit_integer (digits, ndigits, dest, K, neg, precision);
385 }
386
387 /* write decimal w/o scientific notation */
388 if (!scientific || (K < 0 && (K > -7 || exp < 4))) {
389 return emit_fixed_digits (digits, ndigits, dest, K, neg, precision, exp);
390 }
391
392 return emit_scientific_digits (digits, ndigits, dest, K, neg, precision, exp);
393 }
394
filter_special(double fp,char * dest,unsigned precision)395 static int filter_special (double fp, char *dest, unsigned precision)
396 {
397 int nchars = 3;
398 char *d = dest;
399
400 if (fp == 0.0) {
401 if (get_dbits (fp) & signmask) {
402 *d++ = '-';
403 *d++ = '0';
404 }
405 else {
406 *d++ = '0';
407 }
408
409 if (precision) {
410 *d ++ = '.';
411 memset (d, '0', precision);
412 }
413
414 return d - dest + precision;
415 }
416
417 uint64_t bits = get_dbits (fp);
418
419 bool nan = (bits & expmask) == expmask;
420
421 if (!nan) {
422 return 0;
423 }
424
425 if (bits & fracmask) {
426 dest[0] = 'n';
427 dest[1] = 'a';
428 dest[2] = 'n';
429 }
430 else {
431 if (get_dbits (fp) & signmask) {
432 dest[0] = '-';
433 dest[1] = 'i';
434 dest[2] = 'n';
435 dest[3] = 'f';
436 nchars = 4;
437 }
438 else {
439 dest[0] = 'i';
440 dest[1] = 'n';
441 dest[2] = 'f';
442 }
443 }
444
445 return nchars;
446 }
447
448 int
fpconv_dtoa(double d,char dest[FPCONV_BUFLEN],unsigned precision,bool scientific)449 fpconv_dtoa (double d, char dest[FPCONV_BUFLEN],
450 unsigned precision, bool scientific)
451 {
452 char digits[18];
453
454 int str_len = 0;
455 bool neg = false;
456
457 if (precision > FPCONV_BUFLEN - 5) {
458 precision = FPCONV_BUFLEN - 5;
459 }
460
461 int spec = filter_special (d, dest, precision);
462
463 if (spec) {
464 return spec;
465 }
466
467 if (get_dbits (d) & signmask) {
468 dest[0] = '-';
469 str_len++;
470 neg = true;
471 }
472
473 int K = 0;
474 int ndigits = grisu2 (d, digits, &K);
475
476 str_len += emit_digits (digits, ndigits, dest + str_len, K, neg, precision,
477 scientific);
478
479 return str_len;
480 }
481