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