1 /* derived from /netlib/fp/dtoa.c assuming IEEE, Standard C */
2 /* kudos to dmg@bell-labs.com, gripes to ehg@bell-labs.com */
3
4 /* Let x be the exact mathematical number defined by a decimal
5 * string s. Then atof(s) is the round-nearest-even IEEE
6 * floating point value.
7 * Let y be an IEEE floating point value and let s be the string
8 * printed as %.17g. Then atof(s) is exactly y.
9 */
10 #include <u.h>
11 #include <libc.h>
12
13 static Lock _dtoalk[2];
14 #define ACQUIRE_DTOA_LOCK(n) lock(&_dtoalk[n])
15 #define FREE_DTOA_LOCK(n) unlock(&_dtoalk[n])
16
17 #define PRIVATE_mem ((2000+sizeof(double)-1)/sizeof(double))
18 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
19
20 #define FLT_ROUNDS 1
21 #define DBL_DIG 15
22 #define DBL_MAX_10_EXP 308
23 #define DBL_MAX_EXP 1024
24 #define FLT_RADIX 2
25 #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
26
27 /* Ten_pmax = floor(P*log(2)/log(5)) */
28 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
29 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
30 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
31
32 #define Exp_shift 20
33 #define Exp_shift1 20
34 #define Exp_msk1 0x100000
35 #define Exp_msk11 0x100000
36 #define Exp_mask 0x7ff00000
37 #define P 53
38 #define Bias 1023
39 #define Emin (-1022)
40 #define Exp_1 0x3ff00000
41 #define Exp_11 0x3ff00000
42 #define Ebits 11
43 #define Frac_mask 0xfffff
44 #define Frac_mask1 0xfffff
45 #define Ten_pmax 22
46 #define Bletch 0x10
47 #define Bndry_mask 0xfffff
48 #define Bndry_mask1 0xfffff
49 #define LSB 1
50 #define Sign_bit 0x80000000
51 #define Log2P 1
52 #define Tiny0 0
53 #define Tiny1 1
54 #define Quick_max 14
55 #define Int_max 14
56 #define Avoid_Underflow
57
58 #define rounded_product(a,b) a *= b
59 #define rounded_quotient(a,b) a /= b
60
61 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
62 #define Big1 0xffffffff
63
64 #define FFFFFFFF 0xffffffffUL
65
66 #define Kmax 15
67
68 typedef struct Bigint Bigint;
69 typedef struct Ulongs Ulongs;
70
71 struct Bigint {
72 Bigint *next;
73 int k, maxwds, sign, wds;
74 unsigned x[1];
75 };
76
77 struct Ulongs {
78 ulong hi;
79 ulong lo;
80 };
81
82 static Bigint *freelist[Kmax+1];
83
84 Ulongs
double2ulongs(double d)85 double2ulongs(double d)
86 {
87 FPdbleword dw;
88 Ulongs uls;
89
90 dw.x = d;
91 uls.hi = dw.hi;
92 uls.lo = dw.lo;
93 return uls;
94 }
95
96 double
ulongs2double(Ulongs uls)97 ulongs2double(Ulongs uls)
98 {
99 FPdbleword dw;
100
101 dw.hi = uls.hi;
102 dw.lo = uls.lo;
103 return dw.x;
104 }
105
106 static Bigint *
Balloc(int k)107 Balloc(int k)
108 {
109 int x;
110 Bigint * rv;
111 unsigned int len;
112
113 ACQUIRE_DTOA_LOCK(0);
114 if (rv = freelist[k]) {
115 freelist[k] = rv->next;
116 } else {
117 x = 1 << k;
118 len = (sizeof(Bigint) + (x - 1) * sizeof(unsigned int) + sizeof(double) -1)
119 / sizeof(double);
120 if (pmem_next - private_mem + len <= PRIVATE_mem) {
121 rv = (Bigint * )pmem_next;
122 pmem_next += len;
123 } else
124 rv = (Bigint * )malloc(len * sizeof(double));
125 rv->k = k;
126 rv->maxwds = x;
127 }
128 FREE_DTOA_LOCK(0);
129 rv->sign = rv->wds = 0;
130 return rv;
131 }
132
133 static void
Bfree(Bigint * v)134 Bfree(Bigint *v)
135 {
136 if (v) {
137 ACQUIRE_DTOA_LOCK(0);
138 v->next = freelist[v->k];
139 freelist[v->k] = v;
140 FREE_DTOA_LOCK(0);
141 }
142 }
143
144 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
145 y->wds*sizeof(int) + 2*sizeof(int))
146
147 static Bigint *
multadd(Bigint * b,int m,int a)148 multadd(Bigint *b, int m, int a) /* multiply by m and add a */
149 {
150 int i, wds;
151 unsigned int carry, *x, y;
152 unsigned int xi, z;
153 Bigint * b1;
154
155 wds = b->wds;
156 x = b->x;
157 i = 0;
158 carry = a;
159 do {
160 xi = *x;
161 y = (xi & 0xffff) * m + carry;
162 z = (xi >> 16) * m + (y >> 16);
163 carry = z >> 16;
164 *x++ = (z << 16) + (y & 0xffff);
165 } while (++i < wds);
166 if (carry) {
167 if (wds >= b->maxwds) {
168 b1 = Balloc(b->k + 1);
169 Bcopy(b1, b);
170 Bfree(b);
171 b = b1;
172 }
173 b->x[wds++] = carry;
174 b->wds = wds;
175 }
176 return b;
177 }
178
179 static Bigint *
s2b(const char * s,int nd0,int nd,unsigned int y9)180 s2b(const char *s, int nd0, int nd, unsigned int y9)
181 {
182 Bigint * b;
183 int i, k;
184 int x, y;
185
186 x = (nd + 8) / 9;
187 for (k = 0, y = 1; x > y; y <<= 1, k++)
188 ;
189 b = Balloc(k);
190 b->x[0] = y9;
191 b->wds = 1;
192
193 i = 9;
194 if (9 < nd0) {
195 s += 9;
196 do
197 b = multadd(b, 10, *s++ - '0');
198 while (++i < nd0);
199 s++;
200 } else
201 s += 10;
202 for (; i < nd; i++)
203 b = multadd(b, 10, *s++ - '0');
204 return b;
205 }
206
207 static int
hi0bits(register unsigned int x)208 hi0bits(register unsigned int x)
209 {
210 register int k = 0;
211
212 if (!(x & 0xffff0000)) {
213 k = 16;
214 x <<= 16;
215 }
216 if (!(x & 0xff000000)) {
217 k += 8;
218 x <<= 8;
219 }
220 if (!(x & 0xf0000000)) {
221 k += 4;
222 x <<= 4;
223 }
224 if (!(x & 0xc0000000)) {
225 k += 2;
226 x <<= 2;
227 }
228 if (!(x & 0x80000000)) {
229 k++;
230 if (!(x & 0x40000000))
231 return 32;
232 }
233 return k;
234 }
235
236 static int
lo0bits(unsigned int * y)237 lo0bits(unsigned int *y)
238 {
239 register int k;
240 register unsigned int x = *y;
241
242 if (x & 7) {
243 if (x & 1)
244 return 0;
245 if (x & 2) {
246 *y = x >> 1;
247 return 1;
248 }
249 *y = x >> 2;
250 return 2;
251 }
252 k = 0;
253 if (!(x & 0xffff)) {
254 k = 16;
255 x >>= 16;
256 }
257 if (!(x & 0xff)) {
258 k += 8;
259 x >>= 8;
260 }
261 if (!(x & 0xf)) {
262 k += 4;
263 x >>= 4;
264 }
265 if (!(x & 0x3)) {
266 k += 2;
267 x >>= 2;
268 }
269 if (!(x & 1)) {
270 k++;
271 x >>= 1;
272 if (!x & 1)
273 return 32;
274 }
275 *y = x;
276 return k;
277 }
278
279 static Bigint *
i2b(int i)280 i2b(int i)
281 {
282 Bigint * b;
283
284 b = Balloc(1);
285 b->x[0] = i;
286 b->wds = 1;
287 return b;
288 }
289
290 static Bigint *
mult(Bigint * a,Bigint * b)291 mult(Bigint *a, Bigint *b)
292 {
293 Bigint * c;
294 int k, wa, wb, wc;
295 unsigned int * x, *xa, *xae, *xb, *xbe, *xc, *xc0;
296 unsigned int y;
297 unsigned int carry, z;
298 unsigned int z2;
299
300 if (a->wds < b->wds) {
301 c = a;
302 a = b;
303 b = c;
304 }
305 k = a->k;
306 wa = a->wds;
307 wb = b->wds;
308 wc = wa + wb;
309 if (wc > a->maxwds)
310 k++;
311 c = Balloc(k);
312 for (x = c->x, xa = x + wc; x < xa; x++)
313 *x = 0;
314 xa = a->x;
315 xae = xa + wa;
316 xb = b->x;
317 xbe = xb + wb;
318 xc0 = c->x;
319 for (; xb < xbe; xb++, xc0++) {
320 if (y = *xb & 0xffff) {
321 x = xa;
322 xc = xc0;
323 carry = 0;
324 do {
325 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
326 carry = z >> 16;
327 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
328 carry = z2 >> 16;
329 Storeinc(xc, z2, z);
330 } while (x < xae);
331 *xc = carry;
332 }
333 if (y = *xb >> 16) {
334 x = xa;
335 xc = xc0;
336 carry = 0;
337 z2 = *xc;
338 do {
339 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
340 carry = z >> 16;
341 Storeinc(xc, z, z2);
342 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
343 carry = z2 >> 16;
344 } while (x < xae);
345 *xc = z2;
346 }
347 }
348 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc)
349 ;
350 c->wds = wc;
351 return c;
352 }
353
354 static Bigint *p5s;
355
356 static Bigint *
pow5mult(Bigint * b,int k)357 pow5mult(Bigint *b, int k)
358 {
359 Bigint * b1, *p5, *p51;
360 int i;
361 static int p05[3] = {
362 5, 25, 125 };
363
364 if (i = k & 3)
365 b = multadd(b, p05[i-1], 0);
366
367 if (!(k >>= 2))
368 return b;
369 if (!(p5 = p5s)) {
370 /* first time */
371 ACQUIRE_DTOA_LOCK(1);
372 if (!(p5 = p5s)) {
373 p5 = p5s = i2b(625);
374 p5->next = 0;
375 }
376 FREE_DTOA_LOCK(1);
377 }
378 for (; ; ) {
379 if (k & 1) {
380 b1 = mult(b, p5);
381 Bfree(b);
382 b = b1;
383 }
384 if (!(k >>= 1))
385 break;
386 if (!(p51 = p5->next)) {
387 ACQUIRE_DTOA_LOCK(1);
388 if (!(p51 = p5->next)) {
389 p51 = p5->next = mult(p5, p5);
390 p51->next = 0;
391 }
392 FREE_DTOA_LOCK(1);
393 }
394 p5 = p51;
395 }
396 return b;
397 }
398
399 static Bigint *
lshift(Bigint * b,int k)400 lshift(Bigint *b, int k)
401 {
402 int i, k1, n, n1;
403 Bigint * b1;
404 unsigned int * x, *x1, *xe, z;
405
406 n = k >> 5;
407 k1 = b->k;
408 n1 = n + b->wds + 1;
409 for (i = b->maxwds; n1 > i; i <<= 1)
410 k1++;
411 b1 = Balloc(k1);
412 x1 = b1->x;
413 for (i = 0; i < n; i++)
414 *x1++ = 0;
415 x = b->x;
416 xe = x + b->wds;
417 if (k &= 0x1f) {
418 k1 = 32 - k;
419 z = 0;
420 do {
421 *x1++ = *x << k | z;
422 z = *x++ >> k1;
423 } while (x < xe);
424 if (*x1 = z)
425 ++n1;
426 } else
427 do
428 *x1++ = *x++;
429 while (x < xe);
430 b1->wds = n1 - 1;
431 Bfree(b);
432 return b1;
433 }
434
435 static int
cmp(Bigint * a,Bigint * b)436 cmp(Bigint *a, Bigint *b)
437 {
438 unsigned int * xa, *xa0, *xb, *xb0;
439 int i, j;
440
441 i = a->wds;
442 j = b->wds;
443 if (i -= j)
444 return i;
445 xa0 = a->x;
446 xa = xa0 + j;
447 xb0 = b->x;
448 xb = xb0 + j;
449 for (; ; ) {
450 if (*--xa != *--xb)
451 return * xa < *xb ? -1 : 1;
452 if (xa <= xa0)
453 break;
454 }
455 return 0;
456 }
457
458 static Bigint *
diff(Bigint * a,Bigint * b)459 diff(Bigint *a, Bigint *b)
460 {
461 Bigint * c;
462 int i, wa, wb;
463 unsigned int * xa, *xae, *xb, *xbe, *xc;
464 unsigned int borrow, y;
465 unsigned int z;
466
467 i = cmp(a, b);
468 if (!i) {
469 c = Balloc(0);
470 c->wds = 1;
471 c->x[0] = 0;
472 return c;
473 }
474 if (i < 0) {
475 c = a;
476 a = b;
477 b = c;
478 i = 1;
479 } else
480 i = 0;
481 c = Balloc(a->k);
482 c->sign = i;
483 wa = a->wds;
484 xa = a->x;
485 xae = xa + wa;
486 wb = b->wds;
487 xb = b->x;
488 xbe = xb + wb;
489 xc = c->x;
490 borrow = 0;
491 do {
492 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
493 borrow = (y & 0x10000) >> 16;
494 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
495 borrow = (z & 0x10000) >> 16;
496 Storeinc(xc, z, y);
497 } while (xb < xbe);
498 while (xa < xae) {
499 y = (*xa & 0xffff) - borrow;
500 borrow = (y & 0x10000) >> 16;
501 z = (*xa++ >> 16) - borrow;
502 borrow = (z & 0x10000) >> 16;
503 Storeinc(xc, z, y);
504 }
505 while (!*--xc)
506 wa--;
507 c->wds = wa;
508 return c;
509 }
510
511 static double
ulp(double x)512 ulp(double x)
513 {
514 ulong L;
515 Ulongs uls;
516
517 uls = double2ulongs(x);
518 L = (uls.hi & Exp_mask) - (P - 1) * Exp_msk1;
519 return ulongs2double((Ulongs){L, 0});
520 }
521
522 static double
b2d(Bigint * a,int * e)523 b2d(Bigint *a, int *e)
524 {
525 unsigned *xa, *xa0, w, y, z;
526 int k;
527 ulong d0, d1;
528
529 xa0 = a->x;
530 xa = xa0 + a->wds;
531 y = *--xa;
532 k = hi0bits(y);
533 *e = 32 - k;
534 if (k < Ebits) {
535 w = xa > xa0 ? *--xa : 0;
536 d1 = y << (32 - Ebits) + k | w >> Ebits - k;
537 return ulongs2double((Ulongs){Exp_1 | y >> Ebits - k, d1});
538 }
539 z = xa > xa0 ? *--xa : 0;
540 if (k -= Ebits) {
541 d0 = Exp_1 | y << k | z >> 32 - k;
542 y = xa > xa0 ? *--xa : 0;
543 d1 = z << k | y >> 32 - k;
544 } else {
545 d0 = Exp_1 | y;
546 d1 = z;
547 }
548 return ulongs2double((Ulongs){d0, d1});
549 }
550
551 static Bigint *
d2b(double d,int * e,int * bits)552 d2b(double d, int *e, int *bits)
553 {
554 Bigint * b;
555 int de, i, k;
556 unsigned *x, y, z;
557 Ulongs uls;
558
559 b = Balloc(1);
560 x = b->x;
561
562 uls = double2ulongs(d);
563 z = uls.hi & Frac_mask;
564 uls.hi &= 0x7fffffff; /* clear sign bit, which we ignore */
565 de = (int)(uls.hi >> Exp_shift);
566 z |= Exp_msk11;
567 if (y = uls.lo) { /* assignment = */
568 if (k = lo0bits(&y)) { /* assignment = */
569 x[0] = y | z << 32 - k;
570 z >>= k;
571 } else
572 x[0] = y;
573 i = b->wds = (x[1] = z) ? 2 : 1;
574 } else {
575 k = lo0bits(&z);
576 x[0] = z;
577 i = b->wds = 1;
578 k += 32;
579 }
580 USED(i);
581 *e = de - Bias - (P - 1) + k;
582 *bits = P - k;
583 return b;
584 }
585
586 static double
ratio(Bigint * a,Bigint * b)587 ratio(Bigint *a, Bigint *b)
588 {
589 double da, db;
590 int k, ka, kb;
591 Ulongs uls;
592
593 da = b2d(a, &ka);
594 db = b2d(b, &kb);
595 k = ka - kb + 32 * (a->wds - b->wds);
596 if (k > 0) {
597 uls = double2ulongs(da);
598 uls.hi += k * Exp_msk1;
599 da = ulongs2double(uls);
600 } else {
601 k = -k;
602 uls = double2ulongs(db);
603 uls.hi += k * Exp_msk1;
604 db = ulongs2double(uls);
605 }
606 return da / db;
607 }
608
609 static const double
610 tens[] = {
611 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
612 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
613 1e20, 1e21, 1e22
614 };
615
616 static const double
617 bigtens[] = {
618 1e16, 1e32, 1e64, 1e128, 1e256 };
619
620 static const double tinytens[] = {
621 1e-16, 1e-32, 1e-64, 1e-128,
622 9007199254740992.e-256
623 };
624
625 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
626 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
627 #define Scale_Bit 0x10
628 #define n_bigtens 5
629
630 #define NAN_WORD0 0x7ff80000
631
632 #define NAN_WORD1 0
633
634 static int
match(const char ** sp,char * t)635 match(const char **sp, char *t)
636 {
637 int c, d;
638 const char * s = *sp;
639
640 while (d = *t++) {
641 if ((c = *++s) >= 'A' && c <= 'Z')
642 c += 'a' - 'A';
643 if (c != d)
644 return 0;
645 }
646 *sp = s + 1;
647 return 1;
648 }
649
650 static void
gethex(double * rvp,const char ** sp)651 gethex(double *rvp, const char **sp)
652 {
653 unsigned int c, x[2];
654 const char * s;
655 int havedig, udx0, xshift;
656
657 x[0] = x[1] = 0;
658 havedig = xshift = 0;
659 udx0 = 1;
660 s = *sp;
661 while (c = *(const unsigned char * )++s) {
662 if (c >= '0' && c <= '9')
663 c -= '0';
664 else if (c >= 'a' && c <= 'f')
665 c += 10 - 'a';
666 else if (c >= 'A' && c <= 'F')
667 c += 10 - 'A';
668 else if (c <= ' ') {
669 if (udx0 && havedig) {
670 udx0 = 0;
671 xshift = 1;
672 }
673 continue;
674 } else if (/*(*/ c == ')') {
675 *sp = s + 1;
676 break;
677 } else
678 return; /* invalid form: don't change *sp */
679 havedig = 1;
680 if (xshift) {
681 xshift = 0;
682 x[0] = x[1];
683 x[1] = 0;
684 }
685 if (udx0)
686 x[0] = (x[0] << 4) | (x[1] >> 28);
687 x[1] = (x[1] << 4) | c;
688 }
689 if ((x[0] &= 0xfffff) || x[1])
690 *rvp = ulongs2double((Ulongs){Exp_mask | x[0], x[1]});
691 }
692
693 static int
quorem(Bigint * b,Bigint * S)694 quorem(Bigint *b, Bigint *S)
695 {
696 int n;
697 unsigned int * bx, *bxe, q, *sx, *sxe;
698 unsigned int borrow, carry, y, ys;
699 unsigned int si, z, zs;
700
701 n = S->wds;
702 if (b->wds < n)
703 return 0;
704 sx = S->x;
705 sxe = sx + --n;
706 bx = b->x;
707 bxe = bx + n;
708 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
709 if (q) {
710 borrow = 0;
711 carry = 0;
712 do {
713 si = *sx++;
714 ys = (si & 0xffff) * q + carry;
715 zs = (si >> 16) * q + (ys >> 16);
716 carry = zs >> 16;
717 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
718 borrow = (y & 0x10000) >> 16;
719 z = (*bx >> 16) - (zs & 0xffff) - borrow;
720 borrow = (z & 0x10000) >> 16;
721 Storeinc(bx, z, y);
722 } while (sx <= sxe);
723 if (!*bxe) {
724 bx = b->x;
725 while (--bxe > bx && !*bxe)
726 --n;
727 b->wds = n;
728 }
729 }
730 if (cmp(b, S) >= 0) {
731 q++;
732 borrow = 0;
733 carry = 0;
734 bx = b->x;
735 sx = S->x;
736 do {
737 si = *sx++;
738 ys = (si & 0xffff) + carry;
739 zs = (si >> 16) + (ys >> 16);
740 carry = zs >> 16;
741 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
742 borrow = (y & 0x10000) >> 16;
743 z = (*bx >> 16) - (zs & 0xffff) - borrow;
744 borrow = (z & 0x10000) >> 16;
745 Storeinc(bx, z, y);
746 } while (sx <= sxe);
747 bx = b->x;
748 bxe = bx + n;
749 if (!*bxe) {
750 while (--bxe > bx && !*bxe)
751 --n;
752 b->wds = n;
753 }
754 }
755 return q;
756 }
757
758 static char *
rv_alloc(int i)759 rv_alloc(int i)
760 {
761 int j, k, *r;
762
763 j = sizeof(unsigned int);
764 for (k = 0;
765 sizeof(Bigint) - sizeof(unsigned int) - sizeof(int) + j <= i;
766 j <<= 1)
767 k++;
768 r = (int * )Balloc(k);
769 *r = k;
770 return
771 (char *)(r + 1);
772 }
773
774 static char *
nrv_alloc(char * s,char ** rve,int n)775 nrv_alloc(char *s, char **rve, int n)
776 {
777 char *rv, *t;
778
779 t = rv = rv_alloc(n);
780 while (*t = *s++)
781 t++;
782 if (rve)
783 *rve = t;
784 return rv;
785 }
786
787 /* freedtoa(s) must be used to free values s returned by dtoa
788 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
789 * but for consistency with earlier versions of dtoa, it is optional
790 * when MULTIPLE_THREADS is not defined.
791 */
792
793 void
freedtoa(char * s)794 freedtoa(char *s)
795 {
796 Bigint * b = (Bigint * )((int *)s - 1);
797 b->maxwds = 1 << (b->k = *(int * )b);
798 Bfree(b);
799 }
800
801 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
802 *
803 * Inspired by "How to Print Floating-Point Numbers Accurately" by
804 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
805 *
806 * Modifications:
807 * 1. Rather than iterating, we use a simple numeric overestimate
808 * to determine k = floor(log10(d)). We scale relevant
809 * quantities using O(log2(k)) rather than O(k) multiplications.
810 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
811 * try to generate digits strictly left to right. Instead, we
812 * compute with fewer bits and propagate the carry if necessary
813 * when rounding the final digit up. This is often faster.
814 * 3. Under the assumption that input will be rounded nearest,
815 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
816 * That is, we allow equality in stopping tests when the
817 * round-nearest rule will give the same floating-point value
818 * as would satisfaction of the stopping test with strict
819 * inequality.
820 * 4. We remove common factors of powers of 2 from relevant
821 * quantities.
822 * 5. When converting floating-point integers less than 1e16,
823 * we use floating-point arithmetic rather than resorting
824 * to multiple-precision integers.
825 * 6. When asked to produce fewer than 15 digits, we first try
826 * to get by with floating-point arithmetic; we resort to
827 * multiple-precision integer arithmetic only if we cannot
828 * guarantee that the floating-point calculation has given
829 * the correctly rounded result. For k requested digits and
830 * "uniformly" distributed input, the probability is
831 * something like 10^(k-15) that we must resort to the int
832 * calculation.
833 */
834
835 char *
dtoa(double d,int mode,int ndigits,int * decpt,int * sign,char ** rve)836 dtoa(double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
837 {
838 /* Arguments ndigits, decpt, sign are similar to those
839 of ecvt and fcvt; trailing zeros are suppressed from
840 the returned string. If not null, *rve is set to point
841 to the end of the return value. If d is +-Infinity or NaN,
842 then *decpt is set to 9999.
843
844 mode:
845 0 ==> shortest string that yields d when read in
846 and rounded to nearest.
847 1 ==> like 0, but with Steele & White stopping rule;
848 e.g. with IEEE P754 arithmetic , mode 0 gives
849 1e23 whereas mode 1 gives 9.999999999999999e22.
850 2 ==> max(1,ndigits) significant digits. This gives a
851 return value similar to that of ecvt, except
852 that trailing zeros are suppressed.
853 3 ==> through ndigits past the decimal point. This
854 gives a return value similar to that from fcvt,
855 except that trailing zeros are suppressed, and
856 ndigits can be negative.
857 4-9 should give the same return values as 2-3, i.e.,
858 4 <= mode <= 9 ==> same return as mode
859 2 + (mode & 1). These modes are mainly for
860 debugging; often they run slower but sometimes
861 faster than modes 2-3.
862 4,5,8,9 ==> left-to-right digit generation.
863 6-9 ==> don't try fast floating-point estimate
864 (if applicable).
865
866 Values of mode other than 0-9 are treated as mode 0.
867
868 Sufficient space is allocated to the return value
869 to hold the suppressed trailing zeros.
870 */
871
872 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
873 j, j1, k, k0, k_check, L, leftright, m2, m5, s2, s5,
874 spec_case, try_quick;
875 Bigint * b, *b1, *delta, *mlo=nil, *mhi, *S;
876 double d2, ds, eps;
877 char *s, *s0;
878 Ulongs ulsd, ulsd2;
879
880 ulsd = double2ulongs(d);
881 if (ulsd.hi & Sign_bit) {
882 /* set sign for everything, including 0's and NaNs */
883 *sign = 1;
884 ulsd.hi &= ~Sign_bit; /* clear sign bit */
885 } else
886 *sign = 0;
887
888 if ((ulsd.hi & Exp_mask) == Exp_mask) {
889 /* Infinity or NaN */
890 *decpt = 9999;
891 if (!ulsd.lo && !(ulsd.hi & 0xfffff))
892 return nrv_alloc("Infinity", rve, 8);
893 return nrv_alloc("NaN", rve, 3);
894 }
895 d = ulongs2double(ulsd);
896
897 if (!d) {
898 *decpt = 1;
899 return nrv_alloc("0", rve, 1);
900 }
901
902 b = d2b(d, &be, &bbits);
903 i = (int)(ulsd.hi >> Exp_shift1 & (Exp_mask >> Exp_shift1));
904
905 ulsd2 = ulsd;
906 ulsd2.hi &= Frac_mask1;
907 ulsd2.hi |= Exp_11;
908 d2 = ulongs2double(ulsd2);
909
910 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
911 * log10(x) = log(x) / log(10)
912 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
913 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
914 *
915 * This suggests computing an approximation k to log10(d) by
916 *
917 * k = (i - Bias)*0.301029995663981
918 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
919 *
920 * We want k to be too large rather than too small.
921 * The error in the first-order Taylor series approximation
922 * is in our favor, so we just round up the constant enough
923 * to compensate for any error in the multiplication of
924 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
925 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
926 * adding 1e-13 to the constant term more than suffices.
927 * Hence we adjust the constant term to 0.1760912590558.
928 * (We could get a more accurate k by invoking log10,
929 * but this is probably not worthwhile.)
930 */
931
932 i -= Bias;
933 ds = (d2 - 1.5) * 0.289529654602168 + 0.1760912590558 + i * 0.301029995663981;
934 k = (int)ds;
935 if (ds < 0. && ds != k)
936 k--; /* want k = floor(ds) */
937 k_check = 1;
938 if (k >= 0 && k <= Ten_pmax) {
939 if (d < tens[k])
940 k--;
941 k_check = 0;
942 }
943 j = bbits - i - 1;
944 if (j >= 0) {
945 b2 = 0;
946 s2 = j;
947 } else {
948 b2 = -j;
949 s2 = 0;
950 }
951 if (k >= 0) {
952 b5 = 0;
953 s5 = k;
954 s2 += k;
955 } else {
956 b2 -= k;
957 b5 = -k;
958 s5 = 0;
959 }
960 if (mode < 0 || mode > 9)
961 mode = 0;
962 try_quick = 1;
963 if (mode > 5) {
964 mode -= 4;
965 try_quick = 0;
966 }
967 leftright = 1;
968 switch (mode) {
969 case 0:
970 case 1:
971 default:
972 ilim = ilim1 = -1;
973 i = 18;
974 ndigits = 0;
975 break;
976 case 2:
977 leftright = 0;
978 /* no break */
979 case 4:
980 if (ndigits <= 0)
981 ndigits = 1;
982 ilim = ilim1 = i = ndigits;
983 break;
984 case 3:
985 leftright = 0;
986 /* no break */
987 case 5:
988 i = ndigits + k + 1;
989 ilim = i;
990 ilim1 = i - 1;
991 if (i <= 0)
992 i = 1;
993 }
994 s = s0 = rv_alloc(i);
995
996 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
997
998 /* Try to get by with floating-point arithmetic. */
999
1000 i = 0;
1001 d2 = d;
1002 k0 = k;
1003 ilim0 = ilim;
1004 ieps = 2; /* conservative */
1005 if (k > 0) {
1006 ds = tens[k&0xf];
1007 j = k >> 4;
1008 if (j & Bletch) {
1009 /* prevent overflows */
1010 j &= Bletch - 1;
1011 d /= bigtens[n_bigtens-1];
1012 ieps++;
1013 }
1014 for (; j; j >>= 1, i++)
1015 if (j & 1) {
1016 ieps++;
1017 ds *= bigtens[i];
1018 }
1019 d /= ds;
1020 } else if (j1 = -k) {
1021 d *= tens[j1 & 0xf];
1022 for (j = j1 >> 4; j; j >>= 1, i++)
1023 if (j & 1) {
1024 ieps++;
1025 d *= bigtens[i];
1026 }
1027 }
1028 if (k_check && d < 1. && ilim > 0) {
1029 if (ilim1 <= 0)
1030 goto fast_failed;
1031 ilim = ilim1;
1032 k--;
1033 d *= 10.;
1034 ieps++;
1035 }
1036 eps = ieps * d + 7.;
1037
1038 ulsd = double2ulongs(eps);
1039 ulsd.hi -= (P - 1) * Exp_msk1;
1040 eps = ulongs2double(ulsd);
1041
1042 if (ilim == 0) {
1043 S = mhi = 0;
1044 d -= 5.;
1045 if (d > eps)
1046 goto one_digit;
1047 if (d < -eps)
1048 goto no_digits;
1049 goto fast_failed;
1050 }
1051 /* Generate ilim digits, then fix them up. */
1052 eps *= tens[ilim-1];
1053 for (i = 1; ; i++, d *= 10.) {
1054 L = d;
1055 // assert(L < 10);
1056 d -= L;
1057 *s++ = '0' + (int)L;
1058 if (i == ilim) {
1059 if (d > 0.5 + eps)
1060 goto bump_up;
1061 else if (d < 0.5 - eps) {
1062 while (*--s == '0')
1063 ;
1064 s++;
1065 goto ret1;
1066 }
1067 break;
1068 }
1069 }
1070 fast_failed:
1071 s = s0;
1072 d = d2;
1073 k = k0;
1074 ilim = ilim0;
1075 }
1076
1077 /* Do we have a "small" integer? */
1078
1079 if (be >= 0 && k <= Int_max) {
1080 /* Yes. */
1081 ds = tens[k];
1082 if (ndigits < 0 && ilim <= 0) {
1083 S = mhi = 0;
1084 if (ilim < 0 || d <= 5 * ds)
1085 goto no_digits;
1086 goto one_digit;
1087 }
1088 for (i = 1; ; i++) {
1089 L = d / ds;
1090 d -= L * ds;
1091 *s++ = '0' + (int)L;
1092 if (i == ilim) {
1093 d += d;
1094 if (d > ds || d == ds && L & 1) {
1095 bump_up:
1096 while (*--s == '9')
1097 if (s == s0) {
1098 k++;
1099 *s = '0';
1100 break;
1101 }
1102 ++ * s++;
1103 }
1104 break;
1105 }
1106 if (!(d *= 10.))
1107 break;
1108 }
1109 goto ret1;
1110 }
1111
1112 m2 = b2;
1113 m5 = b5;
1114 mhi = mlo = 0;
1115 if (leftright) {
1116 if (mode < 2) {
1117 i =
1118 1 + P - bbits;
1119 } else {
1120 j = ilim - 1;
1121 if (m5 >= j)
1122 m5 -= j;
1123 else {
1124 s5 += j -= m5;
1125 b5 += j;
1126 m5 = 0;
1127 }
1128 if ((i = ilim) < 0) {
1129 m2 -= i;
1130 i = 0;
1131 }
1132 }
1133 b2 += i;
1134 s2 += i;
1135 mhi = i2b(1);
1136 }
1137 if (m2 > 0 && s2 > 0) {
1138 i = m2 < s2 ? m2 : s2;
1139 b2 -= i;
1140 m2 -= i;
1141 s2 -= i;
1142 }
1143 if (b5 > 0) {
1144 if (leftright) {
1145 if (m5 > 0) {
1146 mhi = pow5mult(mhi, m5);
1147 b1 = mult(mhi, b);
1148 Bfree(b);
1149 b = b1;
1150 }
1151 if (j = b5 - m5)
1152 b = pow5mult(b, j);
1153 } else
1154 b = pow5mult(b, b5);
1155 }
1156 S = i2b(1);
1157 if (s5 > 0)
1158 S = pow5mult(S, s5);
1159
1160 /* Check for special case that d is a normalized power of 2. */
1161
1162 spec_case = 0;
1163 if (mode < 2) {
1164 ulsd = double2ulongs(d);
1165 if (!ulsd.lo && !(ulsd.hi & Bndry_mask)) {
1166 /* The special case */
1167 b2 += Log2P;
1168 s2 += Log2P;
1169 spec_case = 1;
1170 }
1171 }
1172
1173 /* Arrange for convenient computation of quotients:
1174 * shift left if necessary so divisor has 4 leading 0 bits.
1175 *
1176 * Perhaps we should just compute leading 28 bits of S once
1177 * and for all and pass them and a shift to quorem, so it
1178 * can do shifts and ors to compute the numerator for q.
1179 */
1180 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
1181 i = 32 - i;
1182 if (i > 4) {
1183 i -= 4;
1184 b2 += i;
1185 m2 += i;
1186 s2 += i;
1187 } else if (i < 4) {
1188 i += 28;
1189 b2 += i;
1190 m2 += i;
1191 s2 += i;
1192 }
1193 if (b2 > 0)
1194 b = lshift(b, b2);
1195 if (s2 > 0)
1196 S = lshift(S, s2);
1197 if (k_check) {
1198 if (cmp(b, S) < 0) {
1199 k--;
1200 b = multadd(b, 10, 0); /* we botched the k estimate */
1201 if (leftright)
1202 mhi = multadd(mhi, 10, 0);
1203 ilim = ilim1;
1204 }
1205 }
1206 if (ilim <= 0 && mode > 2) {
1207 if (ilim < 0 || cmp(b, S = multadd(S, 5, 0)) <= 0) {
1208 /* no digits, fcvt style */
1209 no_digits:
1210 k = -1 - ndigits;
1211 goto ret;
1212 }
1213 one_digit:
1214 *s++ = '1';
1215 k++;
1216 goto ret;
1217 }
1218 if (leftright) {
1219 if (m2 > 0)
1220 mhi = lshift(mhi, m2);
1221
1222 /* Compute mlo -- check for special case
1223 * that d is a normalized power of 2.
1224 */
1225
1226 mlo = mhi;
1227 if (spec_case) {
1228 mhi = Balloc(mhi->k);
1229 Bcopy(mhi, mlo);
1230 mhi = lshift(mhi, Log2P);
1231 }
1232
1233 for (i = 1; ; i++) {
1234 dig = quorem(b, S) + '0';
1235 /* Do we yet have the shortest decimal string
1236 * that will round to d?
1237 */
1238 j = cmp(b, mlo);
1239 delta = diff(S, mhi);
1240 j1 = delta->sign ? 1 : cmp(b, delta);
1241 Bfree(delta);
1242 ulsd = double2ulongs(d);
1243 if (j1 == 0 && !mode && !(ulsd.lo & 1)) {
1244 if (dig == '9')
1245 goto round_9_up;
1246 if (j > 0)
1247 dig++;
1248 *s++ = dig;
1249 goto ret;
1250 }
1251 if (j < 0 || j == 0 && !mode && !(ulsd.lo & 1)) {
1252 if (j1 > 0) {
1253 b = lshift(b, 1);
1254 j1 = cmp(b, S);
1255 if ((j1 > 0 || j1 == 0 && dig & 1)
1256 && dig++ == '9')
1257 goto round_9_up;
1258 }
1259 *s++ = dig;
1260 goto ret;
1261 }
1262 if (j1 > 0) {
1263 if (dig == '9') { /* possible if i == 1 */
1264 round_9_up:
1265 *s++ = '9';
1266 goto roundoff;
1267 }
1268 *s++ = dig + 1;
1269 goto ret;
1270 }
1271 *s++ = dig;
1272 if (i == ilim)
1273 break;
1274 b = multadd(b, 10, 0);
1275 if (mlo == mhi)
1276 mlo = mhi = multadd(mhi, 10, 0);
1277 else {
1278 mlo = multadd(mlo, 10, 0);
1279 mhi = multadd(mhi, 10, 0);
1280 }
1281 }
1282 } else
1283 for (i = 1; ; i++) {
1284 *s++ = dig = quorem(b, S) + '0';
1285 if (i >= ilim)
1286 break;
1287 b = multadd(b, 10, 0);
1288 }
1289
1290 /* Round off last digit */
1291
1292 b = lshift(b, 1);
1293 j = cmp(b, S);
1294 if (j > 0 || j == 0 && dig & 1) {
1295 roundoff:
1296 while (*--s == '9')
1297 if (s == s0) {
1298 k++;
1299 *s++ = '1';
1300 goto ret;
1301 }
1302 ++ * s++;
1303 } else {
1304 while (*--s == '0')
1305 ;
1306 s++;
1307 }
1308 ret:
1309 Bfree(S);
1310 if (mhi) {
1311 if (mlo && mlo != mhi)
1312 Bfree(mlo);
1313 Bfree(mhi);
1314 }
1315 ret1:
1316 Bfree(b);
1317 *s = 0;
1318 *decpt = k + 1;
1319 if (rve)
1320 *rve = s;
1321 return s0;
1322 }
1323