1 #include "fconv.h"
2
3 static int quorem(Bigint *, Bigint *);
4
5 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
6 *
7 * Inspired by "How to Print Floating-Point Numbers Accurately" by
8 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
9 *
10 * Modifications:
11 * 1. Rather than iterating, we use a simple numeric overestimate
12 * to determine k = floor(log10(d)). We scale relevant
13 * quantities using O(log2(k)) rather than O(k) multiplications.
14 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
15 * try to generate digits strictly left to right. Instead, we
16 * compute with fewer bits and propagate the carry if necessary
17 * when rounding the final digit up. This is often faster.
18 * 3. Under the assumption that input will be rounded nearest,
19 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
20 * That is, we allow equality in stopping tests when the
21 * round-nearest rule will give the same floating-point value
22 * as would satisfaction of the stopping test with strict
23 * inequality.
24 * 4. We remove common factors of powers of 2 from relevant
25 * quantities.
26 * 5. When converting floating-point integers less than 1e16,
27 * we use floating-point arithmetic rather than resorting
28 * to multiple-precision integers.
29 * 6. When asked to produce fewer than 15 digits, we first try
30 * to get by with floating-point arithmetic; we resort to
31 * multiple-precision integer arithmetic only if we cannot
32 * guarantee that the floating-point calculation has given
33 * the correctly rounded result. For k requested digits and
34 * "uniformly" distributed input, the probability is
35 * something like 10^(k-15) that we must resort to the long
36 * calculation.
37 */
38
39 char *
_dtoa(double darg,int mode,int ndigits,int * decpt,int * sign,char ** rve)40 _dtoa(double darg, int mode, int ndigits, int *decpt, int *sign, char **rve)
41 {
42 /* Arguments ndigits, decpt, sign are similar to those
43 of ecvt and fcvt; trailing zeros are suppressed from
44 the returned string. If not null, *rve is set to point
45 to the end of the return value. If d is +-Infinity or NaN,
46 then *decpt is set to 9999.
47
48 mode:
49 0 ==> shortest string that yields d when read in
50 and rounded to nearest.
51 1 ==> like 0, but with Steele & White stopping rule;
52 e.g. with IEEE P754 arithmetic , mode 0 gives
53 1e23 whereas mode 1 gives 9.999999999999999e22.
54 2 ==> max(1,ndigits) significant digits. This gives a
55 return value similar to that of ecvt, except
56 that trailing zeros are suppressed.
57 3 ==> through ndigits past the decimal point. This
58 gives a return value similar to that from fcvt,
59 except that trailing zeros are suppressed, and
60 ndigits can be negative.
61 4-9 should give the same return values as 2-3, i.e.,
62 4 <= mode <= 9 ==> same return as mode
63 2 + (mode & 1). These modes are mainly for
64 debugging; often they run slower but sometimes
65 faster than modes 2-3.
66 4,5,8,9 ==> left-to-right digit generation.
67 6-9 ==> don't try fast floating-point estimate
68 (if applicable).
69
70 Values of mode other than 0-9 are treated as mode 0.
71
72 Sufficient space is allocated to the return value
73 to hold the suppressed trailing zeros.
74 */
75
76 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
77 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
78 spec_case, try_quick;
79 long L;
80 #ifndef Sudden_Underflow
81 int denorm;
82 unsigned long x;
83 #endif
84 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
85 double ds;
86 Dul d2, eps;
87 char *s, *s0;
88 static Bigint *result;
89 static int result_k;
90 Dul d;
91
92 d.d = darg;
93 if (result) {
94 result->k = result_k;
95 result->maxwds = 1 << result_k;
96 Bfree(result);
97 result = 0;
98 }
99
100 if (word0(d) & Sign_bit) {
101 /* set sign for everything, including 0's and NaNs */
102 *sign = 1;
103 word0(d) &= ~Sign_bit; /* clear sign bit */
104 }
105 else
106 *sign = 0;
107
108 #if defined(IEEE_Arith) + defined(VAX)
109 #ifdef IEEE_Arith
110 if ((word0(d) & Exp_mask) == Exp_mask)
111 #else
112 if (word0(d) == 0x8000)
113 #endif
114 {
115 /* Infinity or NaN */
116 *decpt = 9999;
117 s =
118 #ifdef IEEE_Arith
119 !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" :
120 #endif
121 "NaN";
122 if (rve)
123 *rve =
124 #ifdef IEEE_Arith
125 s[3] ? s + 8 :
126 #endif
127 s + 3;
128 return s;
129 }
130 #endif
131 #ifdef IBM
132 d.d += 0; /* normalize */
133 #endif
134 if (!d.d) {
135 *decpt = 1;
136 s = "0";
137 if (rve)
138 *rve = s + 1;
139 return s;
140 }
141
142 b = d2b(d.d, &be, &bbits);
143 #ifdef Sudden_Underflow
144 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
145 #else
146 if (i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) {
147 #endif
148 d2.d = d.d;
149 word0(d2) &= Frac_mask1;
150 word0(d2) |= Exp_11;
151 #ifdef IBM
152 if (j = 11 - hi0bits(word0(d2) & Frac_mask))
153 d2.d /= 1 << j;
154 #endif
155
156 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
157 * log10(x) = log(x) / log(10)
158 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
159 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
160 *
161 * This suggests computing an approximation k to log10(d) by
162 *
163 * k = (i - Bias)*0.301029995663981
164 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
165 *
166 * We want k to be too large rather than too small.
167 * The error in the first-order Taylor series approximation
168 * is in our favor, so we just round up the constant enough
169 * to compensate for any error in the multiplication of
170 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
171 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
172 * adding 1e-13 to the constant term more than suffices.
173 * Hence we adjust the constant term to 0.1760912590558.
174 * (We could get a more accurate k by invoking log10,
175 * but this is probably not worthwhile.)
176 */
177
178 i -= Bias;
179 #ifdef IBM
180 i <<= 2;
181 i += j;
182 #endif
183 #ifndef Sudden_Underflow
184 denorm = 0;
185 }
186 else {
187 /* d is denormalized */
188
189 i = bbits + be + (Bias + (P-1) - 1);
190 x = i > 32 ? word0(d) << 64 - i | word1(d) >> i - 32
191 : word1(d) << 32 - i;
192 d2.d = x;
193 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
194 i -= (Bias + (P-1) - 1) + 1;
195 denorm = 1;
196 }
197 #endif
198 ds = (d2.d-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
199 k = floor(ds);
200 k_check = 1;
201 if (k >= 0 && k <= Ten_pmax) {
202 if (d.d < tens[k])
203 k--;
204 k_check = 0;
205 }
206 j = bbits - i - 1;
207 if (j >= 0) {
208 b2 = 0;
209 s2 = j;
210 }
211 else {
212 b2 = -j;
213 s2 = 0;
214 }
215 if (k >= 0) {
216 b5 = 0;
217 s5 = k;
218 s2 += k;
219 }
220 else {
221 b2 -= k;
222 b5 = -k;
223 s5 = 0;
224 }
225 if (mode < 0 || mode > 9)
226 mode = 0;
227 try_quick = 1;
228 if (mode > 5) {
229 mode -= 4;
230 try_quick = 0;
231 }
232 leftright = 1;
233 switch(mode) {
234 case 0:
235 case 1:
236 ilim = ilim1 = -1;
237 i = 18;
238 ndigits = 0;
239 break;
240 case 2:
241 leftright = 0;
242 /* no break */
243 case 4:
244 if (ndigits <= 0)
245 ndigits = 1;
246 ilim = ilim1 = i = ndigits;
247 break;
248 case 3:
249 leftright = 0;
250 /* no break */
251 case 5:
252 i = ndigits + k + 1;
253 ilim = i;
254 ilim1 = i - 1;
255 if (i <= 0)
256 i = 1;
257 }
258 j = sizeof(unsigned long);
259 for(result_k = 0; sizeof(Bigint) - sizeof(unsigned long) + j <= i;
260 j <<= 1) result_k++;
261 result = Balloc(result_k);
262 s = s0 = (char *)result;
263
264 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
265
266 /* Try to get by with floating-point arithmetic. */
267
268 i = 0;
269 d2.d = d.d;
270 k0 = k;
271 ilim0 = ilim;
272 ieps = 2; /* conservative */
273 if (k > 0) {
274 ds = tens[k&0xf];
275 j = k >> 4;
276 if (j & Bletch) {
277 /* prevent overflows */
278 j &= Bletch - 1;
279 d.d /= bigtens[n_bigtens-1];
280 ieps++;
281 }
282 for(; j; j >>= 1, i++)
283 if (j & 1) {
284 ieps++;
285 ds *= bigtens[i];
286 }
287 d.d /= ds;
288 }
289 else if (j1 = -k) {
290 d.d *= tens[j1 & 0xf];
291 for(j = j1 >> 4; j; j >>= 1, i++)
292 if (j & 1) {
293 ieps++;
294 d.d *= bigtens[i];
295 }
296 }
297 if (k_check && d.d < 1. && ilim > 0) {
298 if (ilim1 <= 0)
299 goto fast_failed;
300 ilim = ilim1;
301 k--;
302 d.d *= 10.;
303 ieps++;
304 }
305 eps.d = ieps*d.d + 7.;
306 word0(eps) -= (P-1)*Exp_msk1;
307 if (ilim == 0) {
308 S = mhi = 0;
309 d.d -= 5.;
310 if (d.d > eps.d)
311 goto one_digit;
312 if (d.d < -eps.d)
313 goto no_digits;
314 goto fast_failed;
315 }
316 #ifndef No_leftright
317 if (leftright) {
318 /* Use Steele & White method of only
319 * generating digits needed.
320 */
321 eps.d = 0.5/tens[ilim-1] - eps.d;
322 for(i = 0;;) {
323 L = floor(d.d);
324 d.d -= L;
325 *s++ = '0' + (int)L;
326 if (d.d < eps.d)
327 goto ret1;
328 if (1. - d.d < eps.d)
329 goto bump_up;
330 if (++i >= ilim)
331 break;
332 eps.d *= 10.;
333 d.d *= 10.;
334 }
335 }
336 else {
337 #endif
338 /* Generate ilim digits, then fix them up. */
339 eps.d *= tens[ilim-1];
340 for(i = 1;; i++, d.d *= 10.) {
341 L = floor(d.d);
342 d.d -= L;
343 *s++ = '0' + (int)L;
344 if (i == ilim) {
345 if (d.d > 0.5 + eps.d)
346 goto bump_up;
347 else if (d.d < 0.5 - eps.d) {
348 while(*--s == '0');
349 s++;
350 goto ret1;
351 }
352 break;
353 }
354 }
355 #ifndef No_leftright
356 }
357 #endif
358 fast_failed:
359 s = s0;
360 d.d = d2.d;
361 k = k0;
362 ilim = ilim0;
363 }
364
365 /* Do we have a "small" integer? */
366
367 if (be >= 0 && k <= Int_max) {
368 /* Yes. */
369 ds = tens[k];
370 if (ndigits < 0 && ilim <= 0) {
371 S = mhi = 0;
372 if (ilim < 0 || d.d <= 5*ds)
373 goto no_digits;
374 goto one_digit;
375 }
376 for(i = 1;; i++) {
377 L = floor(d.d / ds);
378 d.d -= L*ds;
379 #ifdef Check_FLT_ROUNDS
380 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
381 if (d.d < 0) {
382 L--;
383 d.d += ds;
384 }
385 #endif
386 *s++ = '0' + (int)L;
387 if (i == ilim) {
388 d.d += d.d;
389 if (d.d > ds || d.d == ds && L & 1) {
390 bump_up:
391 while(*--s == '9')
392 if (s == s0) {
393 k++;
394 *s = '0';
395 break;
396 }
397 ++*s++;
398 }
399 break;
400 }
401 d.d *= 10.;
402 if (d.d == 0.)
403 break;
404 }
405 goto ret1;
406 }
407
408 m2 = b2;
409 m5 = b5;
410 mhi = mlo = 0;
411 if (leftright) {
412 if (mode < 2) {
413 i =
414 #ifndef Sudden_Underflow
415 denorm ? be + (Bias + (P-1) - 1 + 1) :
416 #endif
417 #ifdef IBM
418 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
419 #else
420 1 + P - bbits;
421 #endif
422 }
423 else {
424 j = ilim - 1;
425 if (m5 >= j)
426 m5 -= j;
427 else {
428 s5 += j -= m5;
429 b5 += j;
430 m5 = 0;
431 }
432 if ((i = ilim) < 0) {
433 m2 -= i;
434 i = 0;
435 }
436 }
437 b2 += i;
438 s2 += i;
439 mhi = i2b(1);
440 }
441 if (m2 > 0 && s2 > 0) {
442 i = m2 < s2 ? m2 : s2;
443 b2 -= i;
444 m2 -= i;
445 s2 -= i;
446 }
447 if (b5 > 0) {
448 if (leftright) {
449 if (m5 > 0) {
450 mhi = pow5mult(mhi, m5);
451 b1 = mult(mhi, b);
452 Bfree(b);
453 b = b1;
454 }
455 if (j = b5 - m5)
456 b = pow5mult(b, j);
457 }
458 else
459 b = pow5mult(b, b5);
460 }
461 S = i2b(1);
462 if (s5 > 0)
463 S = pow5mult(S, s5);
464
465 /* Check for special case that d is a normalized power of 2. */
466
467 if (mode < 2) {
468 if (!word1(d) && !(word0(d) & Bndry_mask)
469 #ifndef Sudden_Underflow
470 && word0(d) & Exp_mask
471 #endif
472 ) {
473 /* The special case */
474 b2 += Log2P;
475 s2 += Log2P;
476 spec_case = 1;
477 }
478 else
479 spec_case = 0;
480 }
481
482 /* Arrange for convenient computation of quotients:
483 * shift left if necessary so divisor has 4 leading 0 bits.
484 *
485 * Perhaps we should just compute leading 28 bits of S once
486 * and for all and pass them and a shift to quorem, so it
487 * can do shifts and ors to compute the numerator for q.
488 */
489 #ifdef Pack_32
490 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
491 i = 32 - i;
492 #else
493 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
494 i = 16 - i;
495 #endif
496 if (i > 4) {
497 i -= 4;
498 b2 += i;
499 m2 += i;
500 s2 += i;
501 }
502 else if (i < 4) {
503 i += 28;
504 b2 += i;
505 m2 += i;
506 s2 += i;
507 }
508 if (b2 > 0)
509 b = lshift(b, b2);
510 if (s2 > 0)
511 S = lshift(S, s2);
512 if (k_check) {
513 if (cmp(b,S) < 0) {
514 k--;
515 b = multadd(b, 10, 0); /* we botched the k estimate */
516 if (leftright)
517 mhi = multadd(mhi, 10, 0);
518 ilim = ilim1;
519 }
520 }
521 if (ilim <= 0 && mode > 2) {
522 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
523 /* no digits, fcvt style */
524 no_digits:
525 k = -1 - ndigits;
526 goto ret;
527 }
528 one_digit:
529 *s++ = '1';
530 k++;
531 goto ret;
532 }
533 if (leftright) {
534 if (m2 > 0)
535 mhi = lshift(mhi, m2);
536
537 /* Compute mlo -- check for special case
538 * that d is a normalized power of 2.
539 */
540
541 mlo = mhi;
542 if (spec_case) {
543 mhi = Balloc(mhi->k);
544 Bcopy(mhi, mlo);
545 mhi = lshift(mhi, Log2P);
546 }
547
548 for(i = 1;;i++) {
549 dig = quorem(b,S) + '0';
550 /* Do we yet have the shortest decimal string
551 * that will round to d?
552 */
553 j = cmp(b, mlo);
554 delta = diff(S, mhi);
555 j1 = delta->sign ? 1 : cmp(b, delta);
556 Bfree(delta);
557 #ifndef ROUND_BIASED
558 if (j1 == 0 && !mode && !(word1(d) & 1)) {
559 if (dig == '9')
560 goto round_9_up;
561 if (j > 0)
562 dig++;
563 *s++ = dig;
564 goto ret;
565 }
566 #endif
567 if (j < 0 || j == 0 && !mode
568 #ifndef ROUND_BIASED
569 && !(word1(d) & 1)
570 #endif
571 ) {
572 if (j1 > 0) {
573 b = lshift(b, 1);
574 j1 = cmp(b, S);
575 if ((j1 > 0 || j1 == 0 && dig & 1)
576 && dig++ == '9')
577 goto round_9_up;
578 }
579 *s++ = dig;
580 goto ret;
581 }
582 if (j1 > 0) {
583 if (dig == '9') { /* possible if i == 1 */
584 round_9_up:
585 *s++ = '9';
586 goto roundoff;
587 }
588 *s++ = dig + 1;
589 goto ret;
590 }
591 *s++ = dig;
592 if (i == ilim)
593 break;
594 b = multadd(b, 10, 0);
595 if (mlo == mhi)
596 mlo = mhi = multadd(mhi, 10, 0);
597 else {
598 mlo = multadd(mlo, 10, 0);
599 mhi = multadd(mhi, 10, 0);
600 }
601 }
602 }
603 else
604 for(i = 1;; i++) {
605 *s++ = dig = quorem(b,S) + '0';
606 if (i >= ilim)
607 break;
608 b = multadd(b, 10, 0);
609 }
610
611 /* Round off last digit */
612
613 b = lshift(b, 1);
614 j = cmp(b, S);
615 if (j > 0 || j == 0 && dig & 1) {
616 roundoff:
617 while(*--s == '9')
618 if (s == s0) {
619 k++;
620 *s++ = '1';
621 goto ret;
622 }
623 ++*s++;
624 }
625 else {
626 while(*--s == '0');
627 s++;
628 }
629 ret:
630 Bfree(S);
631 if (mhi) {
632 if (mlo && mlo != mhi)
633 Bfree(mlo);
634 Bfree(mhi);
635 }
636 ret1:
637 Bfree(b);
638 *s = 0;
639 *decpt = k + 1;
640 if (rve)
641 *rve = s;
642 return s0;
643 }
644
645 static int
quorem(Bigint * b,Bigint * S)646 quorem(Bigint *b, Bigint *S)
647 {
648 int n;
649 long borrow, y;
650 unsigned long carry, q, ys;
651 unsigned long *bx, *bxe, *sx, *sxe;
652 #ifdef Pack_32
653 long z;
654 unsigned long si, zs;
655 #endif
656
657 n = S->wds;
658 #ifdef DEBUG
659 /*debug*/ if (b->wds > n)
660 /*debug*/ Bug("oversize b in quorem");
661 #endif
662 if (b->wds < n)
663 return 0;
664 sx = S->x;
665 sxe = sx + --n;
666 bx = b->x;
667 bxe = bx + n;
668 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
669 #ifdef DEBUG
670 /*debug*/ if (q > 9)
671 /*debug*/ Bug("oversized quotient in quorem");
672 #endif
673 if (q) {
674 borrow = 0;
675 carry = 0;
676 do {
677 #ifdef Pack_32
678 si = *sx++;
679 ys = (si & 0xffff) * q + carry;
680 zs = (si >> 16) * q + (ys >> 16);
681 carry = zs >> 16;
682 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
683 borrow = y >> 16;
684 Sign_Extend(borrow, y);
685 z = (*bx >> 16) - (zs & 0xffff) + borrow;
686 borrow = z >> 16;
687 Sign_Extend(borrow, z);
688 Storeinc(bx, z, y);
689 #else
690 ys = *sx++ * q + carry;
691 carry = ys >> 16;
692 y = *bx - (ys & 0xffff) + borrow;
693 borrow = y >> 16;
694 Sign_Extend(borrow, y);
695 *bx++ = y & 0xffff;
696 #endif
697 }
698 while(sx <= sxe);
699 if (!*bxe) {
700 bx = b->x;
701 while(--bxe > bx && !*bxe)
702 --n;
703 b->wds = n;
704 }
705 }
706 if (cmp(b, S) >= 0) {
707 q++;
708 borrow = 0;
709 carry = 0;
710 bx = b->x;
711 sx = S->x;
712 do {
713 #ifdef Pack_32
714 si = *sx++;
715 ys = (si & 0xffff) + carry;
716 zs = (si >> 16) + (ys >> 16);
717 carry = zs >> 16;
718 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
719 borrow = y >> 16;
720 Sign_Extend(borrow, y);
721 z = (*bx >> 16) - (zs & 0xffff) + borrow;
722 borrow = z >> 16;
723 Sign_Extend(borrow, z);
724 Storeinc(bx, z, y);
725 #else
726 ys = *sx++ + carry;
727 carry = ys >> 16;
728 y = *bx - (ys & 0xffff) + borrow;
729 borrow = y >> 16;
730 Sign_Extend(borrow, y);
731 *bx++ = y & 0xffff;
732 #endif
733 }
734 while(sx <= sxe);
735 bx = b->x;
736 bxe = bx + n;
737 if (!*bxe) {
738 while(--bxe > bx && !*bxe)
739 --n;
740 b->wds = n;
741 }
742 }
743 return q;
744 }
745