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