1 /*
2 * go-dtoa.c: double-to-string conversion.
3 *
4 * Copyright 2014 Morten Welinder <terra@gnome.org>
5 *
6 *
7 *
8 * Large portions of this code was copied from http://www.musl-libc.org/
9 * under the following license:
10 *
11 * Copyright © 2005-2014 Rich Felker, et al.
12 *
13 * Permission is hereby granted, free of charge, to any person obtaining
14 * a copy of this software and associated documentation files (the
15 * "Software"), to deal in the Software without restriction, including
16 * without limitation the rights to use, copy, modify, merge, publish,
17 * distribute, sublicense, and/or sell copies of the Software, and to
18 * permit persons to whom the Software is furnished to do so, subject to
19 * the following conditions:
20 *
21 * The above copyright notice and this permission notice shall be
22 * included in all copies or substantial portions of the Software.
23 *
24 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
25 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
26 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
27 * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
28 * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
29 * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
30 * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
31 */
32
33 #include <goffice/goffice-config.h>
34 #include <goffice/math/go-math.h>
35 #include "go-dtoa.h"
36 #include <string.h>
37 #include <inttypes.h>
38 #include <stdarg.h>
39
40 #if (defined(i386) || defined(__i386__) || defined(__i386) || defined(__x86_64__) || defined(__x86_64)) && HAVE_FPU_CONTROL_H
41 #define ENSURE_FPU_STATE
42 #include <fpu_control.h>
43 #endif
44
45 typedef GString FAKE_FILE;
46
47 #ifndef GOFFICE_WITH_LONG_DOUBLE
48 #define go_finitel isfinite
49 #endif
50
51
52 /* musl code starts here */
53
54 /* Some useful macros */
55
56 #ifdef MUSL_ORIGINAL
57 #define MAX(a,b) ((a)>(b) ? (a) : (b))
58 #define MIN(a,b) ((a)<(b) ? (a) : (b))
59 #define CONCAT2(x,y) x ## y
60 #define CONCAT(x,y) CONCAT2(x,y)
61 #endif
62
63 /* Convenient bit representation for modifier flags, which all fall
64 * within 31 codepoints of the space character. */
65
66 #define ALT_FORM (1U<<('#'-' '))
67 #define ZERO_PAD (1U<<('0'-' '))
68 #define LEFT_ADJ (1U<<('-'-' '))
69 #define PAD_POS (1U<<(' '-' '))
70 #define MARK_POS (1U<<('+'-' '))
71 #define GROUPED (1U<<('\''-' '))
72 #ifdef MUSL_ORIGINAL
73 #define FLAGMASK (ALT_FORM|ZERO_PAD|LEFT_ADJ|PAD_POS|MARK_POS|GROUPED)
74 #else
75 #define FLAG_TIES_AWAY_0 (1U<<('"'-' '))
76 #define FLAG_SHORTEST (1U<<('!'-' '))
77 #define FLAG_TRUNCATE (1U<<('='-' '))
78 #define FLAG_ASCII (1U<<(','-' '))
79 #endif
80
81
out(FAKE_FILE * f,const char * s,size_t l)82 static void out(FAKE_FILE *f, const char *s, size_t l)
83 {
84 #ifdef MUSL_ORIGINAL
85 __fwritex((void *)s, l, f);
86 #else
87 g_string_append_len (f, s, l);
88 #endif
89 }
90
91 #ifndef MUSL_ORIGINAL
92 static void
outdecimal(FAKE_FILE * f,int fl)93 outdecimal (FAKE_FILE *f, int fl)
94 {
95 if (fl & FLAG_ASCII)
96 out(f, ".", 1);
97 else {
98 GString const *dec = go_locale_get_decimal();
99 out(f, dec->str, dec->len);
100 }
101 }
102 #endif
103
104
pad(FAKE_FILE * f,char c,int w,int l,int fl)105 static void pad(FAKE_FILE *f, char c, int w, int l, int fl)
106 {
107 char pad[256];
108 if (fl & (LEFT_ADJ | ZERO_PAD) || l >= w) return;
109 l = w - l;
110 #ifdef MUSL_ORIGINAL
111 memset(pad, c, l>sizeof pad ? sizeof pad : l);
112 for (; l >= sizeof pad; l -= sizeof pad)
113 out(f, pad, sizeof pad);
114 #else
115 /* Just kill some warnings */
116 memset(pad, c, (size_t)l>sizeof pad ? sizeof pad : (size_t)l);
117 for (; (size_t)l >= sizeof pad; l -= sizeof pad)
118 out(f, pad, sizeof pad);
119 #endif
120 out(f, pad, l);
121 }
122
123 static const char xdigits[16] = {
124 "0123456789ABCDEF"
125 };
126
fmt_u(uintmax_t x,char * s)127 static char *fmt_u(uintmax_t x, char *s)
128 {
129 unsigned long y;
130 for ( ; x>ULONG_MAX; x/=10) *--s = '0' + x%10;
131 for (y=x; y; y/=10) *--s = '0' + y%10;
132 return s;
133 }
134
135 /* Do not override this check. The floating point printing code below
136 * depends on the float.h constants being right. If they are wrong, it
137 * may overflow the stack. */
138 #if LDBL_MANT_DIG == 53
139 typedef char compiler_defines_long_double_incorrectly[9-(int)sizeof(long double)];
140 #endif
141
fmt_fp(FAKE_FILE * f,long double y,int w,int p,int fl,int t)142 static int fmt_fp(FAKE_FILE *f, long double y, int w, int p, int fl, int t)
143 {
144 uint32_t big[(LDBL_MANT_DIG+28)/29 + 1 // mantissa expansion
145 + (LDBL_MAX_EXP+LDBL_MANT_DIG+28+8)/9]; // exponent expansion
146 uint32_t *a, *d, *r, *z;
147 #ifdef MUSL_ORIGINAL
148 int e2=0, e, i, j, l;
149 #else
150 /* Make i unsigned kills some warnings */
151 int e2=0, e, j, l;
152 unsigned i;
153 #endif
154 char buf[9+LDBL_MANT_DIG/4], *s;
155 const char *prefix="-0X+0X 0X-0x+0x 0x";
156 int pl;
157 char ebuf0[3*sizeof(int)], *ebuf=&ebuf0[3*sizeof(int)], *estr;
158
159 #ifndef MUSL_ORIGINAL
160 if (fl & FLAG_TRUNCATE) g_string_truncate (f, 0);
161 #endif
162
163 pl=1;
164 if (signbit(y)) {
165 y=-y;
166 } else if (fl & MARK_POS) {
167 prefix+=3;
168 } else if (fl & PAD_POS) {
169 prefix+=6;
170 } else prefix++, pl=0;
171
172 if (!go_finitel(y)) {
173 const char *s = (t&32)?"inf":"INF";
174 if (y!=y) s=(t&32)?"nan":"NAN", pl=0;
175 pad(f, ' ', w, 3+pl, fl&~ZERO_PAD);
176 out(f, prefix, pl);
177 out(f, s, 3);
178 pad(f, ' ', w, 3+pl, fl^LEFT_ADJ);
179 return MAX(w, 3+pl);
180 }
181
182 y = frexpl(y, &e2) * 2;
183 if (y) e2--;
184
185 if ((t|32)=='a') {
186 long double round = 8.0;
187 int re;
188
189 if (t&32) prefix += 9;
190 pl += 2;
191
192 if (p<0 || p>=LDBL_MANT_DIG/4-1) re=0;
193 else re=LDBL_MANT_DIG/4-1-p;
194
195 if (re) {
196 while (re--) round*=16;
197 if (*prefix=='-') {
198 y=-y;
199 y-=round;
200 y+=round;
201 y=-y;
202 } else {
203 y+=round;
204 y-=round;
205 }
206 }
207
208 estr=fmt_u(e2<0 ? -e2 : e2, ebuf);
209 if (estr==ebuf) *--estr='0';
210 *--estr = (e2<0 ? '-' : '+');
211 *--estr = t+('p'-'a');
212
213 s=buf;
214 do {
215 int x=y;
216 *s++=xdigits[x]|(t&32);
217 y=16*(y-x);
218 if (s-buf==1 && (y||p>0||(fl&ALT_FORM))) *s++='.';
219 } while (y);
220
221 if (p && s-buf-2 < p)
222 l = (p+2) + (ebuf-estr);
223 else
224 l = (s-buf) + (ebuf-estr);
225
226 pad(f, ' ', w, pl+l, fl);
227 out(f, prefix, pl);
228 pad(f, '0', w, pl+l, fl^ZERO_PAD);
229 out(f, buf, s-buf);
230 pad(f, '0', l-(ebuf-estr)-(s-buf), 0, 0);
231 out(f, estr, ebuf-estr);
232 pad(f, ' ', w, pl+l, fl^LEFT_ADJ);
233 return MAX(w, pl+l);
234 }
235 if (p<0) p=6;
236
237 if (y) y *= 0x1p28, e2-=28;
238
239 if (e2<0) a=r=z=big;
240 else a=r=z=big+sizeof(big)/sizeof(*big) - LDBL_MANT_DIG - 1;
241
242 do {
243 *z = y;
244 y = 1000000000*(y-*z++);
245 } while (y);
246
247 while (e2>0) {
248 uint32_t carry=0;
249 int sh=MIN(29,e2);
250 for (d=z-1; d>=a; d--) {
251 uint64_t x = ((uint64_t)*d<<sh)+carry;
252 *d = x % 1000000000;
253 carry = x / 1000000000;
254 }
255 if (carry) *--a = carry;
256 while (z>a && !z[-1]) z--;
257 e2-=sh;
258 }
259 while (e2<0) {
260 uint32_t carry=0, *b;
261 int sh=MIN(9,-e2), need=1+(p+LDBL_MANT_DIG/3+8)/9;
262 for (d=a; d<z; d++) {
263 uint32_t rm = *d & ((1<<sh)-1);
264 *d = (*d>>sh) + carry;
265 carry = (1000000000>>sh) * rm;
266 }
267 if (!*a) a++;
268 if (carry) *z++ = carry;
269 /* Avoid (slow!) computation past requested precision */
270 b = (t|32)=='f' ? r : a;
271 if (z-b > need) z = b+need;
272 e2+=sh;
273 }
274
275 if (a<z) for (i=10, e=9*(r-a); *a>=i; i*=10, e++);
276 else e=0;
277
278 /* Perform rounding: j is precision after the radix (possibly neg) */
279 j = p - ((t|32)!='f')*e - ((t|32)=='g' && p);
280 if (j < 9*(z-r-1)) {
281 uint32_t x;
282 /* We avoid C's broken division of negative numbers */
283 d = r + 1 + ((j+9*LDBL_MAX_EXP)/9 - LDBL_MAX_EXP);
284 j += 9*LDBL_MAX_EXP;
285 j %= 9;
286 for (i=10, j++; j<9; i*=10, j++);
287 x = *d % i;
288 /* Are there any significant digits past j? */
289 if (x || d+1!=z) {
290 #ifdef MUSL_ORIGINAL
291 long double round = CONCAT(0x1p,LDBL_MANT_DIG);
292 #else
293 long double round = 2 / LDBL_EPSILON;
294 #endif
295 long double small;
296 #ifdef MUSL_ORIGINAL
297 if (*d/i & 1) round += 2;
298 #else
299 /*
300 * For ties-away-from-zero, just pretend we have
301 * an odd digit.
302 */
303 if ((fl & FLAG_TIES_AWAY_0) || (*d/i & 1)) round += 2;
304 #endif
305 if (x<i/2) small=0x0.8p0;
306 else if (x==i/2 && d+1==z) small=0x1.0p0;
307 else small=0x1.8p0;
308 if (pl && *prefix=='-') round*=-1, small*=-1;
309 *d -= x;
310 /* Decide whether to round by probing round+small */
311 if (round+small != round) {
312 *d = *d + i;
313 while (*d > 999999999) {
314 *d--=0;
315 if (d<a) *--a=0;
316 (*d)++;
317 }
318 for (i=10, e=9*(r-a); *a>=i; i*=10, e++);
319 }
320 }
321 if (z>d+1) z=d+1;
322 }
323 for (; z>a && !z[-1]; z--);
324
325 if ((t|32)=='g') {
326 if (!p) p++;
327 if (p>e && e>=-4) {
328 t--;
329 p-=e+1;
330 } else {
331 t-=2;
332 p--;
333 }
334 if (!(fl&ALT_FORM)) {
335 /* Count trailing zeros in last place */
336 if (z>a && z[-1]) for (i=10, j=0; z[-1]%i==0; i*=10, j++);
337 else j=9;
338 if ((t|32)=='f')
339 p = MIN(p,MAX(0,9*(z-r-1)-j));
340 else
341 p = MIN(p,MAX(0,9*(z-r-1)+e-j));
342 }
343 }
344 l = 1 + p + (p || (fl&ALT_FORM));
345 if ((t|32)=='f') {
346 if (e>0) l+=e;
347 } else {
348 estr=fmt_u(e<0 ? -e : e, ebuf);
349 while(ebuf-estr<2) *--estr='0';
350 *--estr = (e<0 ? '-' : '+');
351 *--estr = t;
352 l += ebuf-estr;
353 }
354
355 pad(f, ' ', w, pl+l, fl);
356 out(f, prefix, pl);
357 pad(f, '0', w, pl+l, fl^ZERO_PAD);
358
359 if ((t|32)=='f') {
360 if (a>r) a=r;
361 for (d=a; d<=r; d++) {
362 char *s = fmt_u(*d, buf+9);
363 if (d!=a) while (s>buf) *--s='0';
364 else if (s==buf+9) *--s='0';
365 out(f, s, buf+9-s);
366 }
367 #ifdef MUSL_ORIGINAL
368 if (p || (fl&ALT_FORM)) out(f, ".", 1);
369 #else
370 if (p || (fl&ALT_FORM)) outdecimal(f, fl);
371 #endif
372 for (; d<z && p>0; d++, p-=9) {
373 char *s = fmt_u(*d, buf+9);
374 while (s>buf) *--s='0';
375 out(f, s, MIN(9,p));
376 }
377 pad(f, '0', p+9, 9, 0);
378 } else {
379 if (z<=a) z=a+1;
380 for (d=a; d<z && p>=0; d++) {
381 char *s = fmt_u(*d, buf+9);
382 if (s==buf+9) *--s='0';
383 if (d!=a) while (s>buf) *--s='0';
384 else {
385 out(f, s++, 1);
386 #ifdef MUSL_ORIGINAL
387 if (p>0||(fl&ALT_FORM)) out(f, ".", 1);
388 #else
389 if (p>0||(fl&ALT_FORM)) outdecimal(f, fl);
390 #endif
391 }
392 out(f, s, MIN(buf+9-s, p));
393 p -= buf+9-s;
394 }
395 pad(f, '0', p+18, 18, 0);
396 out(f, estr, ebuf-estr);
397 }
398
399 pad(f, ' ', w, pl+l, fl^LEFT_ADJ);
400
401 return MAX(w, pl+l);
402 }
403
404 /* musl code ends here */
405
406 static void
parse_fmt(const char * fmt,va_list args,gboolean * is_long,int * w,int * p,int * fl,int * t,long double * d)407 parse_fmt (const char *fmt, va_list args, gboolean *is_long,
408 int *w, int *p, int *fl, int *t, long double *d)
409 {
410 *is_long = FALSE;
411 *w = 1;
412 *p = -1;
413 *fl = 0;
414 *d = 0;
415 *t = 'g';
416
417 while (1) {
418 switch (*fmt) {
419 case '0': *fl |= ZERO_PAD; fmt++; continue;
420 case '^': *fl |= FLAG_TIES_AWAY_0; fmt++; continue;
421 case '+': *fl |= MARK_POS; fmt++; continue;
422 case '-': *fl |= LEFT_ADJ; fmt++; continue;
423 case '!': *fl |= FLAG_SHORTEST; fmt++; continue;
424 case '=': *fl |= FLAG_TRUNCATE; fmt++; continue;
425 case ',': *fl |= FLAG_ASCII; fmt++; continue;
426 }
427 break;
428 }
429
430 while (g_ascii_isdigit (*fmt))
431 *w = *w * 10 + (*fmt++ - '0');
432
433 if (*fmt == '.') {
434 if (fmt[1] == '*') {
435 fmt += 2;
436 *p = va_arg (args, int);
437 } else {
438 *p = 0;
439 fmt++;
440 while (g_ascii_isdigit (*fmt))
441 *p = *p * 10 + (*fmt++ - '0');
442 }
443 }
444
445 if (*fmt == 'L') {
446 *is_long = TRUE;
447 fmt++;
448 }
449
450 if (!strchr ("efgaEFGA", *fmt))
451 return;
452 *t = *fmt;
453
454 if (*is_long)
455 *d = va_arg (args, long double);
456 else
457 *d = va_arg (args, double);
458 }
459
460 static long double
strto(const char * str,gboolean is_long,gboolean is_ascii)461 strto (const char *str, gboolean is_long, gboolean is_ascii)
462 {
463 if (is_long) {
464 #ifdef GOFFICE_WITH_LONG_DOUBLE
465 /*
466 * FIXME: ignore ascii flag for now. Given the limited
467 * and checked usage of this function this should still
468 * be safe, if suboptimal.
469 */
470 return go_strtold (str, NULL);
471 #else
472 return go_nan;
473 #endif
474 } else {
475 return is_ascii
476 ? g_ascii_strtod (str, NULL)
477 : go_strtod (str, NULL);
478 }
479 }
480
481
482 void
go_dtoa(GString * dst,const char * fmt,...)483 go_dtoa (GString *dst, const char *fmt, ...)
484 {
485 int w, p, fl, t;
486 va_list args;
487 long double d;
488 gboolean is_long;
489 gboolean debug = FALSE;
490 size_t oldlen;
491 #ifdef ENSURE_FPU_STATE
492 fpu_control_t oldstate;
493 const fpu_control_t mask = _FPU_EXTENDED | _FPU_DOUBLE | _FPU_SINGLE;
494 #endif
495
496 va_start (args, fmt);
497 parse_fmt (fmt, args, &is_long, &w, &p, &fl, &t, &d);
498 va_end (args);
499
500 if (fl & FLAG_SHORTEST) p = is_long ? 20 : 17;
501 oldlen = (fl & FLAG_TRUNCATE) ? 0 : dst->len;
502
503 #ifdef ENSURE_FPU_STATE
504 // fmt_fp depends on "long double" behaving right. That means that the
505 // fpu must not be in round-to-double mode.
506 // This code ought to do nothing on Linux, but Windows and FreeBSD seem
507 // to have round-to-double as default.
508 _FPU_GETCW (oldstate);
509 if ((oldstate & mask) != _FPU_EXTENDED) {
510 fpu_control_t newstate = (oldstate & ~mask) | _FPU_EXTENDED;
511 _FPU_SETCW (newstate);
512 }
513 #endif
514
515 if (debug) g_printerr ("%Lg [%s] t=%c p=%d\n", d, fmt, t, p);
516 fmt_fp (dst, d, w, p, fl, t);
517 if (debug) g_printerr (" --> %s\n", dst->str);
518
519 if (fl & FLAG_SHORTEST) {
520 GString *alt = g_string_new (NULL);
521 const char *dec = (fl & FLAG_ASCII)
522 ? "."
523 : go_locale_get_decimal()->str;
524 while (p > 0) {
525 ssize_t len_dif;
526 const char *dot = strstr (dst->str + oldlen, dec);
527 long double dalt;
528 if (!dot)
529 break;
530
531 /*
532 * This is crude. We have a dot, so try to render
533 * the same number with less precision and check
534 * that the result round-trips.
535 */
536
537 fmt_fp (alt, d, w, p - 1, fl | FLAG_TRUNCATE, t);
538 len_dif = dst->len - oldlen - alt->len;
539
540 if (len_dif == 0) {
541 /* Same string so we've already removed 0s. */
542 break;
543 }
544
545 dalt = strto (alt->str, is_long, (fl & FLAG_ASCII));
546 if (dalt != d)
547 break;
548
549 g_string_truncate (dst, oldlen);
550 go_string_append_gstring (dst, alt);
551 if (debug) g_printerr (" --> %s\n", dst->str);
552
553 p--;
554
555 if (len_dif > 1)
556 break;
557 }
558 g_string_free (alt, TRUE);
559 }
560
561 #ifdef ENSURE_FPU_STATE
562 if ((oldstate & mask) != _FPU_EXTENDED) {
563 _FPU_SETCW (oldstate);
564 }
565 #endif
566 }
567