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 #include "lib9.h"
4
5 #ifdef __APPLE__
6 #pragma clang diagnostic ignored "-Wlogical-op-parentheses"
7 #pragma clang diagnostic ignored "-Wparentheses"
8 #endif
9 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/
10 #define FREE_DTOA_LOCK(n) /*nothing*/
11
12 /* let's provide reasonable defaults for usual implementation of IEEE f.p. */
13 #ifndef DBL_DIG
14 #define DBL_DIG 15
15 #endif
16 #ifndef DBL_MAX_10_EXP
17 #define DBL_MAX_10_EXP 308
18 #endif
19 #ifndef DBL_MAX_EXP
20 #define DBL_MAX_EXP 1024
21 #endif
22 #ifndef FLT_RADIX
23 #define FLT_RADIX 2
24 #endif
25 #ifndef FLT_ROUNDS
26 #define FLT_ROUNDS 1
27 #endif
28 #ifndef Storeinc
29 #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
30 #endif
31
32 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
33
34 #ifdef USE_FPdbleword
35 #define word0(x) ((FPdbleword*)&x)->hi
36 #define word1(x) ((FPdbleword*)&x)->lo
37 #else
38 #ifdef __LITTLE_ENDIAN
39 #define word0(x) ((unsigned long *)&x)[1]
40 #define word1(x) ((unsigned long *)&x)[0]
41 #else
42 #define word0(x) ((unsigned long *)&x)[0]
43 #define word1(x) ((unsigned long *)&x)[1]
44 #endif
45 #endif
46
47 /* #define P DBL_MANT_DIG */
48 /* Ten_pmax = floor(P*log(2)/log(5)) */
49 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
50 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
51 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
52
53 #define Exp_shift 20
54 #define Exp_shift1 20
55 #define Exp_msk1 0x100000
56 #define Exp_msk11 0x100000
57 #define Exp_mask 0x7ff00000
58 #define P 53
59 #define Bias 1023
60 #define Emin (-1022)
61 #define Exp_1 0x3ff00000
62 #define Exp_11 0x3ff00000
63 #define Ebits 11
64 #define Frac_mask 0xfffff
65 #define Frac_mask1 0xfffff
66 #define Ten_pmax 22
67 #define Bletch 0x10
68 #define Bndry_mask 0xfffff
69 #define Bndry_mask1 0xfffff
70 #define LSB 1
71 #define Sign_bit 0x80000000
72 #define Log2P 1
73 #define Tiny0 0
74 #define Tiny1 1
75 #define Quick_max 14
76 #define Int_max 14
77 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
78 #define Avoid_Underflow
79
80 #define rounded_product(a,b) a *= b
81 #define rounded_quotient(a,b) a /= b
82
83 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
84 #define Big1 0xffffffff
85
86 #define Kmax 15
87
88 struct
89 Bigint {
90 struct Bigint *next;
91 int k, maxwds, sign, wds;
92 unsigned long x[1];
93 };
94
95 typedef struct Bigint Bigint;
96
97 static Bigint *freelist[Kmax+1];
98
99 static Bigint *
Balloc(int k)100 Balloc(int k)
101 {
102 int x;
103 Bigint * rv;
104
105 ACQUIRE_DTOA_LOCK(0);
106 if (rv = freelist[k]) {
107 freelist[k] = rv->next;
108 } else {
109 x = 1 << k;
110 rv = (Bigint * )malloc(sizeof(Bigint) + (x - 1) * sizeof(unsigned long));
111 if(rv == nil)
112 return nil;
113 rv->k = k;
114 rv->maxwds = x;
115 }
116 FREE_DTOA_LOCK(0);
117 rv->sign = rv->wds = 0;
118 return rv;
119 }
120
121 static void
Bfree(Bigint * v)122 Bfree(Bigint *v)
123 {
124 if (v) {
125 ACQUIRE_DTOA_LOCK(0);
126 v->next = freelist[v->k];
127 freelist[v->k] = v;
128 FREE_DTOA_LOCK(0);
129 }
130 }
131
132 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
133 y->wds*sizeof(long) + 2*sizeof(int))
134
135 static Bigint *
multadd(Bigint * b,int m,int a)136 multadd(Bigint *b, int m, int a) /* multiply by m and add a */
137 {
138 int i, wds;
139 unsigned long * x, y;
140 unsigned long xi, z;
141 Bigint * b1;
142
143 wds = b->wds;
144 x = b->x;
145 i = 0;
146 do {
147 xi = *x;
148 y = (xi & 0xffff) * m + a;
149 z = (xi >> 16) * m + (y >> 16);
150 a = (int)(z >> 16);
151 *x++ = (z << 16) + (y & 0xffff);
152 } while (++i < wds);
153 if (a) {
154 if (wds >= b->maxwds) {
155 b1 = Balloc(b->k + 1);
156 Bcopy(b1, b);
157 Bfree(b);
158 b = b1;
159 }
160 b->x[wds++] = a;
161 b->wds = wds;
162 }
163 return b;
164 }
165
166 static Bigint *
s2b(const char * s,int nd0,int nd,unsigned long y9)167 s2b(const char *s, int nd0, int nd, unsigned long y9)
168 {
169 Bigint * b;
170 int i, k;
171 long x, y;
172
173 x = (nd + 8) / 9;
174 for (k = 0, y = 1; x > y; y <<= 1, k++)
175 ;
176 b = Balloc(k);
177 b->x[0] = y9;
178 b->wds = 1;
179
180 i = 9;
181 if (9 < nd0) {
182 s += 9;
183 do
184 b = multadd(b, 10, *s++ - '0');
185 while (++i < nd0);
186 s++;
187 } else
188 s += 10;
189 for (; i < nd; i++)
190 b = multadd(b, 10, *s++ - '0');
191 return b;
192 }
193
194 static int
hi0bits(register unsigned long x)195 hi0bits(register unsigned long x)
196 {
197 register int k = 0;
198
199 if (!(x & 0xffff0000)) {
200 k = 16;
201 x <<= 16;
202 }
203 if (!(x & 0xff000000)) {
204 k += 8;
205 x <<= 8;
206 }
207 if (!(x & 0xf0000000)) {
208 k += 4;
209 x <<= 4;
210 }
211 if (!(x & 0xc0000000)) {
212 k += 2;
213 x <<= 2;
214 }
215 if (!(x & 0x80000000)) {
216 k++;
217 if (!(x & 0x40000000))
218 return 32;
219 }
220 return k;
221 }
222
223 static int
lo0bits(unsigned long * y)224 lo0bits(unsigned long *y)
225 {
226 register int k;
227 register unsigned long x = *y;
228
229 if (x & 7) {
230 if (x & 1)
231 return 0;
232 if (x & 2) {
233 *y = x >> 1;
234 return 1;
235 }
236 *y = x >> 2;
237 return 2;
238 }
239 k = 0;
240 if (!(x & 0xffff)) {
241 k = 16;
242 x >>= 16;
243 }
244 if (!(x & 0xff)) {
245 k += 8;
246 x >>= 8;
247 }
248 if (!(x & 0xf)) {
249 k += 4;
250 x >>= 4;
251 }
252 if (!(x & 0x3)) {
253 k += 2;
254 x >>= 2;
255 }
256 if (!(x & 1)) {
257 k++;
258 x >>= 1;
259 if (!x & 1)
260 return 32;
261 }
262 *y = x;
263 return k;
264 }
265
266 static Bigint *
i2b(int i)267 i2b(int i)
268 {
269 Bigint * b;
270
271 b = Balloc(1);
272 b->x[0] = i;
273 b->wds = 1;
274 return b;
275 }
276
277 static Bigint *
mult(Bigint * a,Bigint * b)278 mult(Bigint *a, Bigint *b)
279 {
280 Bigint * c;
281 int k, wa, wb, wc;
282 unsigned long carry, y, z;
283 unsigned long * x, *xa, *xae, *xb, *xbe, *xc, *xc0;
284 unsigned long z2;
285
286 if (a->wds < b->wds) {
287 c = a;
288 a = b;
289 b = c;
290 }
291 k = a->k;
292 wa = a->wds;
293 wb = b->wds;
294 wc = wa + wb;
295 if (wc > a->maxwds)
296 k++;
297 c = Balloc(k);
298 for (x = c->x, xa = x + wc; x < xa; x++)
299 *x = 0;
300 xa = a->x;
301 xae = xa + wa;
302 xb = b->x;
303 xbe = xb + wb;
304 xc0 = c->x;
305 for (; xb < xbe; xb++, xc0++) {
306 if (y = *xb & 0xffff) {
307 x = xa;
308 xc = xc0;
309 carry = 0;
310 do {
311 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
312 carry = z >> 16;
313 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
314 carry = z2 >> 16;
315 Storeinc(xc, z2, z);
316 } while (x < xae);
317 *xc = carry;
318 }
319 if (y = *xb >> 16) {
320 x = xa;
321 xc = xc0;
322 carry = 0;
323 z2 = *xc;
324 do {
325 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
326 carry = z >> 16;
327 Storeinc(xc, z, z2);
328 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
329 carry = z2 >> 16;
330 } while (x < xae);
331 *xc = z2;
332 }
333 }
334 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc)
335 ;
336 c->wds = wc;
337 return c;
338 }
339
340 static Bigint *p5s;
341
342 static Bigint *
pow5mult(Bigint * b,int k)343 pow5mult(Bigint *b, int k)
344 {
345 Bigint * b1, *p5, *p51;
346 int i;
347 static int p05[3] = {
348 5, 25, 125 };
349
350 if (i = k & 3)
351 b = multadd(b, p05[i-1], 0);
352
353 if (!(k >>= 2))
354 return b;
355 if (!(p5 = p5s)) {
356 /* first time */
357 ACQUIRE_DTOA_LOCK(1);
358 if (!(p5 = p5s)) {
359 p5 = p5s = i2b(625);
360 p5->next = 0;
361 }
362 FREE_DTOA_LOCK(1);
363 }
364 for (; ; ) {
365 if (k & 1) {
366 b1 = mult(b, p5);
367 Bfree(b);
368 b = b1;
369 }
370 if (!(k >>= 1))
371 break;
372 if (!(p51 = p5->next)) {
373 ACQUIRE_DTOA_LOCK(1);
374 if (!(p51 = p5->next)) {
375 p51 = p5->next = mult(p5, p5);
376 p51->next = 0;
377 }
378 FREE_DTOA_LOCK(1);
379 }
380 p5 = p51;
381 }
382 return b;
383 }
384
385 static Bigint *
lshift(Bigint * b,int k)386 lshift(Bigint *b, int k)
387 {
388 int i, k1, n, n1;
389 Bigint * b1;
390 unsigned long * x, *x1, *xe, z;
391
392 n = k >> 5;
393 k1 = b->k;
394 n1 = n + b->wds + 1;
395 for (i = b->maxwds; n1 > i; i <<= 1)
396 k1++;
397 b1 = Balloc(k1);
398 x1 = b1->x;
399 for (i = 0; i < n; i++)
400 *x1++ = 0;
401 x = b->x;
402 xe = x + b->wds;
403 if (k &= 0x1f) {
404 k1 = 32 - k;
405 z = 0;
406 do {
407 *x1++ = *x << k | z;
408 z = *x++ >> k1;
409 } while (x < xe);
410 if (*x1 = z)
411 ++n1;
412 } else
413 do
414 *x1++ = *x++;
415 while (x < xe);
416 b1->wds = n1 - 1;
417 Bfree(b);
418 return b1;
419 }
420
421 static int
cmp(Bigint * a,Bigint * b)422 cmp(Bigint *a, Bigint *b)
423 {
424 unsigned long * xa, *xa0, *xb, *xb0;
425 int i, j;
426
427 i = a->wds;
428 j = b->wds;
429 if (i -= j)
430 return i;
431 xa0 = a->x;
432 xa = xa0 + j;
433 xb0 = b->x;
434 xb = xb0 + j;
435 for (; ; ) {
436 if (*--xa != *--xb)
437 return * xa < *xb ? -1 : 1;
438 if (xa <= xa0)
439 break;
440 }
441 return 0;
442 }
443
444 static Bigint *
diff(Bigint * a,Bigint * b)445 diff(Bigint *a, Bigint *b)
446 {
447 Bigint * c;
448 int i, wa, wb;
449 long borrow, y; /* We need signed shifts here. */
450 unsigned long * xa, *xae, *xb, *xbe, *xc;
451 long z;
452
453 i = cmp(a, b);
454 if (!i) {
455 c = Balloc(0);
456 c->wds = 1;
457 c->x[0] = 0;
458 return c;
459 }
460 if (i < 0) {
461 c = a;
462 a = b;
463 b = c;
464 i = 1;
465 } else
466 i = 0;
467 c = Balloc(a->k);
468 c->sign = i;
469 wa = a->wds;
470 xa = a->x;
471 xae = xa + wa;
472 wb = b->wds;
473 xb = b->x;
474 xbe = xb + wb;
475 xc = c->x;
476 borrow = 0;
477 do {
478 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
479 borrow = y >> 16;
480 Sign_Extend(borrow, y);
481 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
482 borrow = z >> 16;
483 Sign_Extend(borrow, z);
484 Storeinc(xc, z, y);
485 } while (xb < xbe);
486 while (xa < xae) {
487 y = (*xa & 0xffff) + borrow;
488 borrow = y >> 16;
489 Sign_Extend(borrow, y);
490 z = (*xa++ >> 16) + borrow;
491 borrow = z >> 16;
492 Sign_Extend(borrow, z);
493 Storeinc(xc, z, y);
494 }
495 while (!*--xc)
496 wa--;
497 c->wds = wa;
498 return c;
499 }
500
501 static double
ulp(double x)502 ulp(double x)
503 {
504 register long L;
505 double a;
506
507 L = (word0(x) & Exp_mask) - (P - 1) * Exp_msk1;
508 #ifndef Sudden_Underflow
509 if (L > 0) {
510 #endif
511 word0(a) = L;
512 word1(a) = 0;
513 #ifndef Sudden_Underflow
514 } else {
515 L = -L >> Exp_shift;
516 if (L < Exp_shift) {
517 word0(a) = 0x80000 >> L;
518 word1(a) = 0;
519 } else {
520 word0(a) = 0;
521 L -= Exp_shift;
522 word1(a) = L >= 31 ? 1 : 1 << 31 - L;
523 }
524 }
525 #endif
526 return a;
527 }
528
529 static double
b2d(Bigint * a,int * e)530 b2d(Bigint *a, int *e)
531 {
532 unsigned long * xa, *xa0, w, y, z;
533 int k;
534 double d;
535 #define d0 word0(d)
536 #define d1 word1(d)
537
538 xa0 = a->x;
539 xa = xa0 + a->wds;
540 y = *--xa;
541 k = hi0bits(y);
542 *e = 32 - k;
543 if (k < Ebits) {
544 d0 = Exp_1 | y >> Ebits - k;
545 w = xa > xa0 ? *--xa : 0;
546 d1 = y << (32 - Ebits) + k | w >> Ebits - k;
547 goto ret_d;
548 }
549 z = xa > xa0 ? *--xa : 0;
550 if (k -= Ebits) {
551 d0 = Exp_1 | y << k | z >> 32 - k;
552 y = xa > xa0 ? *--xa : 0;
553 d1 = z << k | y >> 32 - k;
554 } else {
555 d0 = Exp_1 | y;
556 d1 = z;
557 }
558 ret_d:
559 #undef d0
560 #undef d1
561 return d;
562 }
563
564 static Bigint *
d2b(double d,int * e,int * bits)565 d2b(double d, int *e, int *bits)
566 {
567 Bigint * b;
568 int de, i, k;
569 unsigned long * x, y, z;
570 #define d0 word0(d)
571 #define d1 word1(d)
572
573 b = Balloc(1);
574 x = b->x;
575
576 z = d0 & Frac_mask;
577 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
578 #ifdef Sudden_Underflow
579 de = (int)(d0 >> Exp_shift);
580 z |= Exp_msk11;
581 #else
582 if (de = (int)(d0 >> Exp_shift))
583 z |= Exp_msk1;
584 #endif
585 if (y = d1) {
586 if (k = lo0bits(&y)) {
587 x[0] = y | z << 32 - k;
588 z >>= k;
589 } else
590 x[0] = y;
591 i = b->wds = (x[1] = z) ? 2 : 1;
592 } else {
593 k = lo0bits(&z);
594 x[0] = z;
595 i = b->wds = 1;
596 k += 32;
597 }
598 #ifndef Sudden_Underflow
599 if (de) {
600 #endif
601 *e = de - Bias - (P - 1) + k;
602 *bits = P - k;
603 #ifndef Sudden_Underflow
604 } else {
605 *e = de - Bias - (P - 1) + 1 + k;
606 *bits = 32 * i - hi0bits(x[i-1]);
607 }
608 #endif
609 return b;
610 }
611
612 #undef d0
613 #undef d1
614
615 static double
ratio(Bigint * a,Bigint * b)616 ratio(Bigint *a, Bigint *b)
617 {
618 double da, db;
619 int k, ka, kb;
620
621 da = b2d(a, &ka);
622 db = b2d(b, &kb);
623 k = ka - kb + 32 * (a->wds - b->wds);
624 if (k > 0)
625 word0(da) += k * Exp_msk1;
626 else {
627 k = -k;
628 word0(db) += k * Exp_msk1;
629 }
630 return da / db;
631 }
632
633 static const double
634 tens[] = {
635 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
636 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
637 1e20, 1e21, 1e22
638 };
639
640 static const double
641 bigtens[] = {
642 1e16, 1e32, 1e64, 1e128, 1e256 };
643
644 static const double tinytens[] = {
645 1e-16, 1e-32, 1e-64, 1e-128,
646 9007199254740992.e-256
647 };
648
649 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
650 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
651 #define Scale_Bit 0x10
652 #define n_bigtens 5
653
654 #define NAN_WORD0 0x7ff80000
655
656 #define NAN_WORD1 0
657
658 static int
match(const char ** sp,char * t)659 match(const char **sp, char *t)
660 {
661 int c, d;
662 const char * s = *sp;
663
664 while (d = *t++) {
665 if ((c = *++s) >= 'A' && c <= 'Z')
666 c += 'a' - 'A';
667 if (c != d)
668 return 0;
669 }
670 *sp = s + 1;
671 return 1;
672 }
673
674 double
strtod(const char * s00,char ** se)675 strtod(const char *s00, char **se)
676 {
677 int scale;
678 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
679 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
680 const char * s, *s0, *s1;
681 double aadj, aadj1, adj, rv, rv0;
682 long L;
683 unsigned long y, z;
684 Bigint * bb, *bb1, *bd, *bd0, *bs, *delta;
685 sign = nz0 = nz = 0;
686 rv = 0.;
687 for (s = s00; ; s++)
688 switch (*s) {
689 case '-':
690 sign = 1;
691 /* no break */
692 case '+':
693 if (*++s)
694 goto break2;
695 /* no break */
696 case 0:
697 s = s00;
698 goto ret;
699 case '\t':
700 case '\n':
701 case '\v':
702 case '\f':
703 case '\r':
704 case ' ':
705 continue;
706 default:
707 goto break2;
708 }
709 break2:
710 if (*s == '0') {
711 nz0 = 1;
712 while (*++s == '0')
713 ;
714 if (!*s)
715 goto ret;
716 }
717 s0 = s;
718 y = z = 0;
719 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
720 if (nd < 9)
721 y = 10 * y + c - '0';
722 else if (nd < 16)
723 z = 10 * z + c - '0';
724 nd0 = nd;
725 if (c == '.') {
726 c = *++s;
727 if (!nd) {
728 for (; c == '0'; c = *++s)
729 nz++;
730 if (c > '0' && c <= '9') {
731 s0 = s;
732 nf += nz;
733 nz = 0;
734 goto have_dig;
735 }
736 goto dig_done;
737 }
738 for (; c >= '0' && c <= '9'; c = *++s) {
739 have_dig:
740 nz++;
741 if (c -= '0') {
742 nf += nz;
743 for (i = 1; i < nz; i++)
744 if (nd++ < 9)
745 y *= 10;
746 else if (nd <= DBL_DIG + 1)
747 z *= 10;
748 if (nd++ < 9)
749 y = 10 * y + c;
750 else if (nd <= DBL_DIG + 1)
751 z = 10 * z + c;
752 nz = 0;
753 }
754 }
755 }
756 dig_done:
757 e = 0;
758 if (c == 'e' || c == 'E') {
759 if (!nd && !nz && !nz0) {
760 s = s00;
761 goto ret;
762 }
763 s00 = s;
764 esign = 0;
765 switch (c = *++s) {
766 case '-':
767 esign = 1;
768 case '+':
769 c = *++s;
770 }
771 if (c >= '0' && c <= '9') {
772 while (c == '0')
773 c = *++s;
774 if (c > '0' && c <= '9') {
775 L = c - '0';
776 s1 = s;
777 while ((c = *++s) >= '0' && c <= '9')
778 L = 10 * L + c - '0';
779 if (s - s1 > 8 || L > 19999)
780 /* Avoid confusion from exponents
781 * so large that e might overflow.
782 */
783 e = 19999; /* safe for 16 bit ints */
784 else
785 e = (int)L;
786 if (esign)
787 e = -e;
788 } else
789 e = 0;
790 } else
791 s = s00;
792 }
793 if (!nd) {
794 if (!nz && !nz0) {
795 /* Check for Nan and Infinity */
796 switch (c) {
797 case 'i':
798 case 'I':
799 if (match(&s, "nfinity")) {
800 word0(rv) = 0x7ff00000;
801 word1(rv) = 0;
802 goto ret;
803 }
804 break;
805 case 'n':
806 case 'N':
807 if (match(&s, "an")) {
808 word0(rv) = NAN_WORD0;
809 word1(rv) = NAN_WORD1;
810 goto ret;
811 }
812 }
813 s = s00;
814 }
815 goto ret;
816 }
817 e1 = e -= nf;
818
819 /* Now we have nd0 digits, starting at s0, followed by a
820 * decimal point, followed by nd-nd0 digits. The number we're
821 * after is the integer represented by those digits times
822 * 10**e */
823
824 if (!nd0)
825 nd0 = nd;
826 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
827 rv = y;
828 if (k > 9)
829 rv = tens[k - 9] * rv + z;
830 bd0 = 0;
831 if (nd <= DBL_DIG
832 && FLT_ROUNDS == 1
833 ) {
834 if (!e)
835 goto ret;
836 if (e > 0) {
837 if (e <= Ten_pmax) {
838 /* rv = */ rounded_product(rv, tens[e]);
839 goto ret;
840 }
841 i = DBL_DIG - nd;
842 if (e <= Ten_pmax + i) {
843 /* A fancier test would sometimes let us do
844 * this for larger i values.
845 */
846 e -= i;
847 rv *= tens[i];
848 /* rv = */ rounded_product(rv, tens[e]);
849 goto ret;
850 }
851 } else if (e >= -Ten_pmax) {
852 /* rv = */ rounded_quotient(rv, tens[-e]);
853 goto ret;
854 }
855 }
856 e1 += nd - k;
857
858 scale = 0;
859
860 /* Get starting approximation = rv * 10**e1 */
861
862 if (e1 > 0) {
863 if (i = e1 & 15)
864 rv *= tens[i];
865 if (e1 &= ~15) {
866 if (e1 > DBL_MAX_10_EXP) {
867 ovfl:
868 /* Can't trust HUGE_VAL */
869 word0(rv) = Exp_mask;
870 word1(rv) = 0;
871 if (bd0)
872 goto retfree;
873 goto ret;
874 }
875 if (e1 >>= 4) {
876 for (j = 0; e1 > 1; j++, e1 >>= 1)
877 if (e1 & 1)
878 rv *= bigtens[j];
879 /* The last multiplication could overflow. */
880 word0(rv) -= P * Exp_msk1;
881 rv *= bigtens[j];
882 if ((z = word0(rv) & Exp_mask)
883 > Exp_msk1 * (DBL_MAX_EXP + Bias - P))
884 goto ovfl;
885 if (z > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P)) {
886 /* set to largest number */
887 /* (Can't trust DBL_MAX) */
888 word0(rv) = Big0;
889 word1(rv) = Big1;
890 } else
891 word0(rv) += P * Exp_msk1;
892 }
893
894 }
895 } else if (e1 < 0) {
896 e1 = -e1;
897 if (i = e1 & 15)
898 rv /= tens[i];
899 if (e1 &= ~15) {
900 e1 >>= 4;
901 if (e1 >= 1 << n_bigtens)
902 goto undfl;
903 if (e1 & Scale_Bit)
904 scale = P;
905 for (j = 0; e1 > 0; j++, e1 >>= 1)
906 if (e1 & 1)
907 rv *= tinytens[j];
908 if (!rv) {
909 undfl:
910 rv = 0.;
911 if (bd0)
912 goto retfree;
913 goto ret;
914 }
915 }
916 }
917
918 /* Now the hard part -- adjusting rv to the correct value.*/
919
920 /* Put digits into bd: true value = bd * 10^e */
921
922 bd0 = s2b(s0, nd0, nd, y);
923
924 for (; ; ) {
925 bd = Balloc(bd0->k);
926 Bcopy(bd, bd0);
927 bb = d2b(rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
928 bs = i2b(1);
929
930 if (e >= 0) {
931 bb2 = bb5 = 0;
932 bd2 = bd5 = e;
933 } else {
934 bb2 = bb5 = -e;
935 bd2 = bd5 = 0;
936 }
937 if (bbe >= 0)
938 bb2 += bbe;
939 else
940 bd2 -= bbe;
941 bs2 = bb2;
942 #ifdef Sudden_Underflow
943 j = P + 1 - bbbits;
944 #else
945 i = bbe + bbbits - 1; /* logb(rv) */
946 if (i < Emin) /* denormal */
947 j = bbe + (P - Emin);
948 else
949 j = P + 1 - bbbits;
950 #endif
951 bb2 += j;
952 bd2 += j;
953 bd2 += scale;
954 i = bb2 < bd2 ? bb2 : bd2;
955 if (i > bs2)
956 i = bs2;
957 if (i > 0) {
958 bb2 -= i;
959 bd2 -= i;
960 bs2 -= i;
961 }
962 if (bb5 > 0) {
963 bs = pow5mult(bs, bb5);
964 bb1 = mult(bs, bb);
965 Bfree(bb);
966 bb = bb1;
967 }
968 if (bb2 > 0)
969 bb = lshift(bb, bb2);
970 if (bd5 > 0)
971 bd = pow5mult(bd, bd5);
972 if (bd2 > 0)
973 bd = lshift(bd, bd2);
974 if (bs2 > 0)
975 bs = lshift(bs, bs2);
976 delta = diff(bb, bd);
977 dsign = delta->sign;
978 delta->sign = 0;
979 i = cmp(delta, bs);
980 if (i < 0) {
981 /* Error is less than half an ulp -- check for
982 * special case of mantissa a power of two.
983 */
984 if (dsign || word1(rv) || word0(rv) & Bndry_mask
985 || (word0(rv) & Exp_mask) <= Exp_msk1
986 ) {
987 if (!delta->x[0] && delta->wds == 1)
988 dsign = 2;
989 break;
990 }
991 delta = lshift(delta, Log2P);
992 if (cmp(delta, bs) > 0)
993 goto drop_down;
994 break;
995 }
996 if (i == 0) {
997 /* exactly half-way between */
998 if (dsign) {
999 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
1000 && word1(rv) == 0xffffffff) {
1001 /*boundary case -- increment exponent*/
1002 word0(rv) = (word0(rv) & Exp_mask)
1003 + Exp_msk1
1004 ;
1005 word1(rv) = 0;
1006 dsign = 0;
1007 break;
1008 }
1009 } else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
1010 dsign = 2;
1011 drop_down:
1012 /* boundary case -- decrement exponent */
1013 #ifdef Sudden_Underflow
1014 L = word0(rv) & Exp_mask;
1015 if (L <= Exp_msk1)
1016 goto undfl;
1017 L -= Exp_msk1;
1018 #else
1019 L = (word0(rv) & Exp_mask) - Exp_msk1;
1020 #endif
1021 word0(rv) = L | Bndry_mask1;
1022 word1(rv) = 0xffffffff;
1023 break;
1024 }
1025 if (!(word1(rv) & LSB))
1026 break;
1027 if (dsign)
1028 rv += ulp(rv);
1029 else {
1030 rv -= ulp(rv);
1031 #ifndef Sudden_Underflow
1032 if (!rv)
1033 goto undfl;
1034 #endif
1035 }
1036 dsign = 1 - dsign;
1037 break;
1038 }
1039 if ((aadj = ratio(delta, bs)) <= 2.) {
1040 if (dsign)
1041 aadj = aadj1 = 1.;
1042 else if (word1(rv) || word0(rv) & Bndry_mask) {
1043 #ifndef Sudden_Underflow
1044 if (word1(rv) == Tiny1 && !word0(rv))
1045 goto undfl;
1046 #endif
1047 aadj = 1.;
1048 aadj1 = -1.;
1049 } else {
1050 /* special case -- power of FLT_RADIX to be */
1051 /* rounded down... */
1052
1053 if (aadj < 2. / FLT_RADIX)
1054 aadj = 1. / FLT_RADIX;
1055 else
1056 aadj *= 0.5;
1057 aadj1 = -aadj;
1058 }
1059 } else {
1060 aadj *= 0.5;
1061 aadj1 = dsign ? aadj : -aadj;
1062 if (FLT_ROUNDS == 0)
1063 aadj1 += 0.5;
1064 }
1065 y = word0(rv) & Exp_mask;
1066
1067 /* Check for overflow */
1068
1069 if (y == Exp_msk1 * (DBL_MAX_EXP + Bias - 1)) {
1070 rv0 = rv;
1071 word0(rv) -= P * Exp_msk1;
1072 adj = aadj1 * ulp(rv);
1073 rv += adj;
1074 if ((word0(rv) & Exp_mask) >=
1075 Exp_msk1 * (DBL_MAX_EXP + Bias - P)) {
1076 if (word0(rv0) == Big0 && word1(rv0) == Big1)
1077 goto ovfl;
1078 word0(rv) = Big0;
1079 word1(rv) = Big1;
1080 goto cont;
1081 } else
1082 word0(rv) += P * Exp_msk1;
1083 } else {
1084 #ifdef Sudden_Underflow
1085 if ((word0(rv) & Exp_mask) <= P * Exp_msk1) {
1086 rv0 = rv;
1087 word0(rv) += P * Exp_msk1;
1088 adj = aadj1 * ulp(rv);
1089 rv += adj;
1090 if ((word0(rv) & Exp_mask) <= P * Exp_msk1) {
1091 if (word0(rv0) == Tiny0
1092 && word1(rv0) == Tiny1)
1093 goto undfl;
1094 word0(rv) = Tiny0;
1095 word1(rv) = Tiny1;
1096 goto cont;
1097 } else
1098 word0(rv) -= P * Exp_msk1;
1099 } else {
1100 adj = aadj1 * ulp(rv);
1101 rv += adj;
1102 }
1103 #else
1104 /* Compute adj so that the IEEE rounding rules will
1105 * correctly round rv + adj in some half-way cases.
1106 * If rv * ulp(rv) is denormalized (i.e.,
1107 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1108 * trouble from bits lost to denormalization;
1109 * example: 1.2e-307 .
1110 */
1111 if (y <= (P - 1) * Exp_msk1 && aadj >= 1.) {
1112 aadj1 = (double)(int)(aadj + 0.5);
1113 if (!dsign)
1114 aadj1 = -aadj1;
1115 }
1116 adj = aadj1 * ulp(rv);
1117 rv += adj;
1118 #endif
1119 }
1120 z = word0(rv) & Exp_mask;
1121 if (!scale)
1122 if (y == z) {
1123 /* Can we stop now? */
1124 L = aadj;
1125 aadj -= L;
1126 /* The tolerances below are conservative. */
1127 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
1128 if (aadj < .4999999 || aadj > .5000001)
1129 break;
1130 } else if (aadj < .4999999 / FLT_RADIX)
1131 break;
1132 }
1133 cont:
1134 Bfree(bb);
1135 Bfree(bd);
1136 Bfree(bs);
1137 Bfree(delta);
1138 }
1139 if (scale) {
1140 if ((word0(rv) & Exp_mask) <= P * Exp_msk1
1141 && word1(rv) & 1
1142 && dsign != 2)
1143 if (dsign)
1144 rv += ulp(rv);
1145 else
1146 word1(rv) &= ~1;
1147 word0(rv0) = Exp_1 - P * Exp_msk1;
1148 word1(rv0) = 0;
1149 rv *= rv0;
1150 }
1151 retfree:
1152 Bfree(bb);
1153 Bfree(bd);
1154 Bfree(bs);
1155 Bfree(bd0);
1156 Bfree(delta);
1157 ret:
1158 if (se)
1159 *se = (char *)s;
1160 return sign ? -rv : rv;
1161 }
1162
1163 static int
quorem(Bigint * b,Bigint * S)1164 quorem(Bigint *b, Bigint *S)
1165 {
1166 int n;
1167 long borrow, y;
1168 unsigned long carry, q, ys;
1169 unsigned long * bx, *bxe, *sx, *sxe;
1170 long z;
1171 unsigned long si, zs;
1172
1173 n = S->wds;
1174 if (b->wds < n)
1175 return 0;
1176 sx = S->x;
1177 sxe = sx + --n;
1178 bx = b->x;
1179 bxe = bx + n;
1180 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1181 if (q) {
1182 borrow = 0;
1183 carry = 0;
1184 do {
1185 si = *sx++;
1186 ys = (si & 0xffff) * q + carry;
1187 zs = (si >> 16) * q + (ys >> 16);
1188 carry = zs >> 16;
1189 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1190 borrow = y >> 16;
1191 Sign_Extend(borrow, y);
1192 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1193 borrow = z >> 16;
1194 Sign_Extend(borrow, z);
1195 Storeinc(bx, z, y);
1196 } while (sx <= sxe);
1197 if (!*bxe) {
1198 bx = b->x;
1199 while (--bxe > bx && !*bxe)
1200 --n;
1201 b->wds = n;
1202 }
1203 }
1204 if (cmp(b, S) >= 0) {
1205 q++;
1206 borrow = 0;
1207 carry = 0;
1208 bx = b->x;
1209 sx = S->x;
1210 do {
1211 si = *sx++;
1212 ys = (si & 0xffff) + carry;
1213 zs = (si >> 16) + (ys >> 16);
1214 carry = zs >> 16;
1215 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1216 borrow = y >> 16;
1217 Sign_Extend(borrow, y);
1218 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1219 borrow = z >> 16;
1220 Sign_Extend(borrow, z);
1221 Storeinc(bx, z, y);
1222 } while (sx <= sxe);
1223 bx = b->x;
1224 bxe = bx + n;
1225 if (!*bxe) {
1226 while (--bxe > bx && !*bxe)
1227 --n;
1228 b->wds = n;
1229 }
1230 }
1231 return q;
1232 }
1233
1234 static char *
rv_alloc(int i)1235 rv_alloc(int i)
1236 {
1237 int j, k, *r;
1238
1239 j = sizeof(unsigned long);
1240 for (k = 0;
1241 sizeof(Bigint) - sizeof(unsigned long) - sizeof(int) + j <= i;
1242 j <<= 1)
1243 k++;
1244 r = (int * )Balloc(k);
1245 *r = k;
1246 return
1247 (char *)(r + 1);
1248 }
1249
1250 static char *
nrv_alloc(char * s,char ** rve,int n)1251 nrv_alloc(char *s, char **rve, int n)
1252 {
1253 char *rv, *t;
1254
1255 t = rv = rv_alloc(n);
1256 while (*t = *s++)
1257 t++;
1258 if (rve)
1259 *rve = t;
1260 return rv;
1261 }
1262
1263 /* freedtoa(s) must be used to free values s returned by dtoa
1264 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
1265 * but for consistency with earlier versions of dtoa, it is optional
1266 * when MULTIPLE_THREADS is not defined.
1267 */
1268
1269 void
freedtoa(char * s)1270 freedtoa(char *s)
1271 {
1272 Bigint * b = (Bigint * )((int *)s - 1);
1273 b->maxwds = 1 << (b->k = *(int * )b);
1274 Bfree(b);
1275 }
1276
1277 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1278 *
1279 * Inspired by "How to Print Floating-Point Numbers Accurately" by
1280 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1281 *
1282 * Modifications:
1283 * 1. Rather than iterating, we use a simple numeric overestimate
1284 * to determine k = floor(log10(d)). We scale relevant
1285 * quantities using O(log2(k)) rather than O(k) multiplications.
1286 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1287 * try to generate digits strictly left to right. Instead, we
1288 * compute with fewer bits and propagate the carry if necessary
1289 * when rounding the final digit up. This is often faster.
1290 * 3. Under the assumption that input will be rounded nearest,
1291 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1292 * That is, we allow equality in stopping tests when the
1293 * round-nearest rule will give the same floating-point value
1294 * as would satisfaction of the stopping test with strict
1295 * inequality.
1296 * 4. We remove common factors of powers of 2 from relevant
1297 * quantities.
1298 * 5. When converting floating-point integers less than 1e16,
1299 * we use floating-point arithmetic rather than resorting
1300 * to multiple-precision integers.
1301 * 6. When asked to produce fewer than 15 digits, we first try
1302 * to get by with floating-point arithmetic; we resort to
1303 * multiple-precision integer arithmetic only if we cannot
1304 * guarantee that the floating-point calculation has given
1305 * the correctly rounded result. For k requested digits and
1306 * "uniformly" distributed input, the probability is
1307 * something like 10^(k-15) that we must resort to the long
1308 * calculation.
1309 */
1310
1311 char *
dtoa(double d,int mode,int ndigits,int * decpt,int * sign,char ** rve)1312 dtoa(double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
1313 {
1314 /* Arguments ndigits, decpt, sign are similar to those
1315 of ecvt and fcvt; trailing zeros are suppressed from
1316 the returned string. If not null, *rve is set to point
1317 to the end of the return value. If d is +-Infinity or NaN,
1318 then *decpt is set to 9999.
1319
1320 mode:
1321 0 ==> shortest string that yields d when read in
1322 and rounded to nearest.
1323 1 ==> like 0, but with Steele & White stopping rule;
1324 e.g. with IEEE P754 arithmetic , mode 0 gives
1325 1e23 whereas mode 1 gives 9.999999999999999e22.
1326 2 ==> max(1,ndigits) significant digits. This gives a
1327 return value similar to that of ecvt, except
1328 that trailing zeros are suppressed.
1329 3 ==> through ndigits past the decimal point. This
1330 gives a return value similar to that from fcvt,
1331 except that trailing zeros are suppressed, and
1332 ndigits can be negative.
1333 4-9 should give the same return values as 2-3, i.e.,
1334 4 <= mode <= 9 ==> same return as mode
1335 2 + (mode & 1). These modes are mainly for
1336 debugging; often they run slower but sometimes
1337 faster than modes 2-3.
1338 4,5,8,9 ==> left-to-right digit generation.
1339 6-9 ==> don't try fast floating-point estimate
1340 (if applicable).
1341
1342 Values of mode other than 0-9 are treated as mode 0.
1343
1344 Sufficient space is allocated to the return value
1345 to hold the suppressed trailing zeros.
1346 */
1347
1348 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
1349 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1350 spec_case, try_quick;
1351 long L;
1352 #ifndef Sudden_Underflow
1353 int denorm;
1354 unsigned long x;
1355 #endif
1356 Bigint * b, *b1, *delta, *mlo, *mhi, *S;
1357 double d2, ds, eps;
1358 char *s, *s0;
1359
1360 if (word0(d) & Sign_bit) {
1361 /* set sign for everything, including 0's and NaNs */
1362 *sign = 1;
1363 word0(d) &= ~Sign_bit; /* clear sign bit */
1364 } else
1365 *sign = 0;
1366
1367 if ((word0(d) & Exp_mask) == Exp_mask) {
1368 /* Infinity or NaN */
1369 *decpt = 9999;
1370 if (!word1(d) && !(word0(d) & 0xfffff))
1371 return nrv_alloc("Infinity", rve, 8);
1372 return nrv_alloc("NaN", rve, 3);
1373 }
1374 if (!d) {
1375 *decpt = 1;
1376 return nrv_alloc("0", rve, 1);
1377 }
1378
1379 b = d2b(d, &be, &bbits);
1380 #ifdef Sudden_Underflow
1381 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask >> Exp_shift1));
1382 #else
1383 if (i = (int)(word0(d) >> Exp_shift1 & (Exp_mask >> Exp_shift1))) {
1384 #endif
1385 word0(d2) = (word0(d) & Frac_mask1) | Exp_11;
1386 word1(d2) = word1(d);
1387
1388 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
1389 * log10(x) = log(x) / log(10)
1390 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1391 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
1392 *
1393 * This suggests computing an approximation k to log10(d) by
1394 *
1395 * k = (i - Bias)*0.301029995663981
1396 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1397 *
1398 * We want k to be too large rather than too small.
1399 * The error in the first-order Taylor series approximation
1400 * is in our favor, so we just round up the constant enough
1401 * to compensate for any error in the multiplication of
1402 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
1403 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
1404 * adding 1e-13 to the constant term more than suffices.
1405 * Hence we adjust the constant term to 0.1760912590558.
1406 * (We could get a more accurate k by invoking log10,
1407 * but this is probably not worthwhile.)
1408 */
1409
1410 i -= Bias;
1411 #ifndef Sudden_Underflow
1412 denorm = 0;
1413 } else {
1414 /* d is denormalized */
1415
1416 i = bbits + be + (Bias + (P - 1) - 1);
1417 x = i > 32 ? word0(d) << 64 - i | word1(d) >> i - 32
1418 : word1(d) << 32 - i;
1419 d2 = x;
1420 word0(d2) -= 31 * Exp_msk1; /* adjust exponent */
1421 i -= (Bias + (P - 1) - 1) + 1;
1422 denorm = 1;
1423 }
1424 #endif
1425 ds = (d2 - 1.5) * 0.289529654602168 + 0.1760912590558 + i * 0.301029995663981;
1426 k = (int)ds;
1427 if (ds < 0. && ds != k)
1428 k--; /* want k = floor(ds) */
1429 k_check = 1;
1430 if (k >= 0 && k <= Ten_pmax) {
1431 if (d < tens[k])
1432 k--;
1433 k_check = 0;
1434 }
1435 j = bbits - i - 1;
1436 if (j >= 0) {
1437 b2 = 0;
1438 s2 = j;
1439 } else {
1440 b2 = -j;
1441 s2 = 0;
1442 }
1443 if (k >= 0) {
1444 b5 = 0;
1445 s5 = k;
1446 s2 += k;
1447 } else {
1448 b2 -= k;
1449 b5 = -k;
1450 s5 = 0;
1451 }
1452 if (mode < 0 || mode > 9)
1453 mode = 0;
1454 try_quick = 1;
1455 if (mode > 5) {
1456 mode -= 4;
1457 try_quick = 0;
1458 }
1459 leftright = 1;
1460 switch (mode) {
1461 case 0:
1462 case 1:
1463 ilim = ilim1 = -1;
1464 i = 18;
1465 ndigits = 0;
1466 break;
1467 case 2:
1468 leftright = 0;
1469 /* no break */
1470 case 4:
1471 if (ndigits <= 0)
1472 ndigits = 1;
1473 ilim = ilim1 = i = ndigits;
1474 break;
1475 case 3:
1476 leftright = 0;
1477 /* no break */
1478 case 5:
1479 i = ndigits + k + 1;
1480 ilim = i;
1481 ilim1 = i - 1;
1482 if (i <= 0)
1483 i = 1;
1484 }
1485 s = s0 = rv_alloc(i);
1486
1487 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
1488
1489 /* Try to get by with floating-point arithmetic. */
1490
1491 i = 0;
1492 d2 = d;
1493 k0 = k;
1494 ilim0 = ilim;
1495 ieps = 2; /* conservative */
1496 if (k > 0) {
1497 ds = tens[k&0xf];
1498 j = k >> 4;
1499 if (j & Bletch) {
1500 /* prevent overflows */
1501 j &= Bletch - 1;
1502 d /= bigtens[n_bigtens-1];
1503 ieps++;
1504 }
1505 for (; j; j >>= 1, i++)
1506 if (j & 1) {
1507 ieps++;
1508 ds *= bigtens[i];
1509 }
1510 d /= ds;
1511 } else if (j1 = -k) {
1512 d *= tens[j1 & 0xf];
1513 for (j = j1 >> 4; j; j >>= 1, i++)
1514 if (j & 1) {
1515 ieps++;
1516 d *= bigtens[i];
1517 }
1518 }
1519 if (k_check && d < 1. && ilim > 0) {
1520 if (ilim1 <= 0)
1521 goto fast_failed;
1522 ilim = ilim1;
1523 k--;
1524 d *= 10.;
1525 ieps++;
1526 }
1527 eps = ieps * d + 7.;
1528 word0(eps) -= (P - 1) * Exp_msk1;
1529 if (ilim == 0) {
1530 S = mhi = 0;
1531 d -= 5.;
1532 if (d > eps)
1533 goto one_digit;
1534 if (d < -eps)
1535 goto no_digits;
1536 goto fast_failed;
1537 }
1538 /* Generate ilim digits, then fix them up. */
1539 eps *= tens[ilim-1];
1540 for (i = 1; ; i++, d *= 10.) {
1541 L = d;
1542 d -= L;
1543 *s++ = '0' + (int)L;
1544 if (i == ilim) {
1545 if (d > 0.5 + eps)
1546 goto bump_up;
1547 else if (d < 0.5 - eps) {
1548 while (*--s == '0')
1549 ;
1550 s++;
1551 goto ret1;
1552 }
1553 break;
1554 }
1555 }
1556 fast_failed:
1557 s = s0;
1558 d = d2;
1559 k = k0;
1560 ilim = ilim0;
1561 }
1562
1563 /* Do we have a "small" integer? */
1564
1565 if (be >= 0 && k <= Int_max) {
1566 /* Yes. */
1567 ds = tens[k];
1568 if (ndigits < 0 && ilim <= 0) {
1569 S = mhi = 0;
1570 if (ilim < 0 || d <= 5 * ds)
1571 goto no_digits;
1572 goto one_digit;
1573 }
1574 for (i = 1; ; i++) {
1575 L = d / ds;
1576 d -= L * ds;
1577 *s++ = '0' + (int)L;
1578 if (i == ilim) {
1579 d += d;
1580 if (d > ds || d == ds && L & 1) {
1581 bump_up:
1582 while (*--s == '9')
1583 if (s == s0) {
1584 k++;
1585 *s = '0';
1586 break;
1587 }
1588 ++ * s++;
1589 }
1590 break;
1591 }
1592 if (!(d *= 10.))
1593 break;
1594 }
1595 goto ret1;
1596 }
1597
1598 m2 = b2;
1599 m5 = b5;
1600 mhi = mlo = 0;
1601 if (leftright) {
1602 if (mode < 2) {
1603 i =
1604 #ifndef Sudden_Underflow
1605 denorm ? be + (Bias + (P - 1) - 1 + 1) :
1606 #endif
1607 1 + P - bbits;
1608 } else {
1609 j = ilim - 1;
1610 if (m5 >= j)
1611 m5 -= j;
1612 else {
1613 s5 += j -= m5;
1614 b5 += j;
1615 m5 = 0;
1616 }
1617 if ((i = ilim) < 0) {
1618 m2 -= i;
1619 i = 0;
1620 }
1621 }
1622 b2 += i;
1623 s2 += i;
1624 mhi = i2b(1);
1625 }
1626 if (m2 > 0 && s2 > 0) {
1627 i = m2 < s2 ? m2 : s2;
1628 b2 -= i;
1629 m2 -= i;
1630 s2 -= i;
1631 }
1632 if (b5 > 0) {
1633 if (leftright) {
1634 if (m5 > 0) {
1635 mhi = pow5mult(mhi, m5);
1636 b1 = mult(mhi, b);
1637 Bfree(b);
1638 b = b1;
1639 }
1640 if (j = b5 - m5)
1641 b = pow5mult(b, j);
1642 } else
1643 b = pow5mult(b, b5);
1644 }
1645 S = i2b(1);
1646 if (s5 > 0)
1647 S = pow5mult(S, s5);
1648
1649 /* Check for special case that d is a normalized power of 2. */
1650
1651 spec_case = 0;
1652 if (mode < 2) {
1653 if (!word1(d) && !(word0(d) & Bndry_mask)
1654 #ifndef Sudden_Underflow
1655 && word0(d) & Exp_mask
1656 #endif
1657 ) {
1658 /* The special case */
1659 b2 += Log2P;
1660 s2 += Log2P;
1661 spec_case = 1;
1662 }
1663 }
1664
1665 /* Arrange for convenient computation of quotients:
1666 * shift left if necessary so divisor has 4 leading 0 bits.
1667 *
1668 * Perhaps we should just compute leading 28 bits of S once
1669 * and for all and pass them and a shift to quorem, so it
1670 * can do shifts and ors to compute the numerator for q.
1671 */
1672 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
1673 i = 32 - i;
1674 if (i > 4) {
1675 i -= 4;
1676 b2 += i;
1677 m2 += i;
1678 s2 += i;
1679 } else if (i < 4) {
1680 i += 28;
1681 b2 += i;
1682 m2 += i;
1683 s2 += i;
1684 }
1685 if (b2 > 0)
1686 b = lshift(b, b2);
1687 if (s2 > 0)
1688 S = lshift(S, s2);
1689 if (k_check) {
1690 if (cmp(b, S) < 0) {
1691 k--;
1692 b = multadd(b, 10, 0); /* we botched the k estimate */
1693 if (leftright)
1694 mhi = multadd(mhi, 10, 0);
1695 ilim = ilim1;
1696 }
1697 }
1698 if (ilim <= 0 && mode > 2) {
1699 if (ilim < 0 || cmp(b, S = multadd(S, 5, 0)) <= 0) {
1700 /* no digits, fcvt style */
1701 no_digits:
1702 k = -1 - ndigits;
1703 goto ret;
1704 }
1705 one_digit:
1706 *s++ = '1';
1707 k++;
1708 goto ret;
1709 }
1710 if (leftright) {
1711 if (m2 > 0)
1712 mhi = lshift(mhi, m2);
1713
1714 /* Compute mlo -- check for special case
1715 * that d is a normalized power of 2.
1716 */
1717
1718 mlo = mhi;
1719 if (spec_case) {
1720 mhi = Balloc(mhi->k);
1721 Bcopy(mhi, mlo);
1722 mhi = lshift(mhi, Log2P);
1723 }
1724
1725 for (i = 1; ; i++) {
1726 dig = quorem(b, S) + '0';
1727 /* Do we yet have the shortest decimal string
1728 * that will round to d?
1729 */
1730 j = cmp(b, mlo);
1731 delta = diff(S, mhi);
1732 j1 = delta->sign ? 1 : cmp(b, delta);
1733 Bfree(delta);
1734 if (j1 == 0 && !mode && !(word1(d) & 1)) {
1735 if (dig == '9')
1736 goto round_9_up;
1737 if (j > 0)
1738 dig++;
1739 *s++ = dig;
1740 goto ret;
1741 }
1742 if (j < 0 || j == 0 && !mode
1743 && !(word1(d) & 1)
1744 ) {
1745 if (j1 > 0) {
1746 b = lshift(b, 1);
1747 j1 = cmp(b, S);
1748 if ((j1 > 0 || j1 == 0 && dig & 1)
1749 && dig++ == '9')
1750 goto round_9_up;
1751 }
1752 *s++ = dig;
1753 goto ret;
1754 }
1755 if (j1 > 0) {
1756 if (dig == '9') { /* possible if i == 1 */
1757 round_9_up:
1758 *s++ = '9';
1759 goto roundoff;
1760 }
1761 *s++ = dig + 1;
1762 goto ret;
1763 }
1764 *s++ = dig;
1765 if (i == ilim)
1766 break;
1767 b = multadd(b, 10, 0);
1768 if (mlo == mhi)
1769 mlo = mhi = multadd(mhi, 10, 0);
1770 else {
1771 mlo = multadd(mlo, 10, 0);
1772 mhi = multadd(mhi, 10, 0);
1773 }
1774 }
1775 } else
1776 for (i = 1; ; i++) {
1777 *s++ = dig = quorem(b, S) + '0';
1778 if (i >= ilim)
1779 break;
1780 b = multadd(b, 10, 0);
1781 }
1782
1783 /* Round off last digit */
1784
1785 b = lshift(b, 1);
1786 j = cmp(b, S);
1787 if (j > 0 || j == 0 && dig & 1) {
1788 roundoff:
1789 while (*--s == '9')
1790 if (s == s0) {
1791 k++;
1792 *s++ = '1';
1793 goto ret;
1794 }
1795 ++ * s++;
1796 } else {
1797 while (*--s == '0')
1798 ;
1799 s++;
1800 }
1801 ret:
1802 Bfree(S);
1803 if (mhi) {
1804 if (mlo && mlo != mhi)
1805 Bfree(mlo);
1806 Bfree(mhi);
1807 }
1808 ret1:
1809 Bfree(b);
1810 *s = 0;
1811 *decpt = k + 1;
1812 if (rve)
1813 *rve = s;
1814 return s0;
1815 }
1816
1817