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