1 /****************************************************************
3 The author of this software is David M. Gay.
5 Copyright (C) 1998, 1999 by Lucent Technologies
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
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
27 ****************************************************************/
29 /* Please send bug reports to David M. Gay (dmg at acm dot org,
30 * with " at " changed at "@" and " dot " changed to "."). */
34 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
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].
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
53 * 4. We remove common factors of powers of 2 from relevant
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
68 #ifdef Honor_FLT_ROUNDS
69 #undef Check_FLT_ROUNDS
70 #define Check_FLT_ROUNDS
72 #define Rounding Flt_Rounds
78 (d0
, mode
, ndigits
, decpt
, sign
, rve
)
79 double d0
; int mode
, ndigits
, *decpt
, *sign
; char **rve
;
81 (double d0
, int mode
, int ndigits
, int *decpt
, int *sign
, char **rve
)
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.
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).
112 Values of mode other than 0-9 are treated as mode 0.
114 Sufficient space is allocated to the return value
115 to hold the suppressed trailing zeros.
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
;
122 #ifndef Sudden_Underflow
126 Bigint
*b
, *b1
, *delta
, *mlo
, *mhi
, *S
;
131 int inexact
, oldinexact
;
133 #ifdef Honor_FLT_ROUNDS /*{*/
135 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
136 Rounding
= Flt_Rounds
;
139 switch(fegetround()) {
140 case FE_TOWARDZERO
: Rounding
= 0; break;
141 case FE_UPWARD
: Rounding
= 2; break;
142 case FE_DOWNWARD
: Rounding
= 3;
147 #ifndef MULTIPLE_THREADS
149 freedtoa(dtoa_result
);
154 if (word0(&d
) & Sign_bit
) {
155 /* set sign for everything, including 0's and NaNs */
157 word0(&d
) &= ~Sign_bit
; /* clear sign bit */
162 #if defined(IEEE_Arith) + defined(VAX)
164 if ((word0(&d
) & Exp_mask
) == Exp_mask
)
166 if (word0(&d
) == 0x8000)
169 /* Infinity or NaN */
172 if (!word1(&d
) && !(word0(&d
) & 0xfffff))
173 return nrv_alloc("Infinity", rve
, 8);
175 return nrv_alloc("NaN", rve
, 3);
179 dval(&d
) += 0; /* normalize */
183 return nrv_alloc("0", rve
, 1);
187 try_quick
= oldinexact
= get_inexact();
190 #ifdef Honor_FLT_ROUNDS
193 Rounding
= Rounding
== 2 ? 0 : 2;
200 b
= d2b(dval(&d
), &be
, &bbits
);
201 #ifdef Sudden_Underflow
202 i
= (int)(word0(&d
) >> Exp_shift1
& (Exp_mask
>>Exp_shift1
));
204 if (( i
= (int)(word0(&d
) >> Exp_shift1
& (Exp_mask
>>Exp_shift1
)) )!=0) {
206 dval(&d2
) = dval(&d
);
207 word0(&d2
) &= Frac_mask1
;
208 word0(&d2
) |= Exp_11
;
210 if (( j
= 11 - hi0bits(word0(&d2
) & Frac_mask
) )!=0)
214 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
215 * log10(x) = log(x) / log(10)
216 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
217 * log10(&d) = (i-Bias)*log(2)/log(10) + log10(&d2)
219 * This suggests computing an approximation k to log10(&d) by
221 * k = (i - Bias)*0.301029995663981
222 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
224 * We want k to be too large rather than too small.
225 * The error in the first-order Taylor series approximation
226 * is in our favor, so we just round up the constant enough
227 * to compensate for any error in the multiplication of
228 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
229 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
230 * adding 1e-13 to the constant term more than suffices.
231 * Hence we adjust the constant term to 0.1760912590558.
232 * (We could get a more accurate k by invoking log10,
233 * but this is probably not worthwhile.)
241 #ifndef Sudden_Underflow
245 /* d is denormalized */
247 i
= bbits
+ be
+ (Bias
+ (P
-1) - 1);
248 x
= i
> 32 ? word0(&d
) << (64 - i
) | word1(&d
) >> (i
- 32)
249 : word1(&d
) << (32 - i
);
251 word0(&d2
) -= 31*Exp_msk1
; /* adjust exponent */
252 i
-= (Bias
+ (P
-1) - 1) + 1;
256 ds
= (dval(&d2
)-1.5)*0.289529654602168 + 0.1760912590558 + i
*0.301029995663981;
258 if (ds
< 0. && ds
!= k
)
259 k
--; /* want k = floor(ds) */
261 if (k
>= 0 && k
<= Ten_pmax
) {
262 if (dval(&d
) < tens
[k
])
285 if (mode
< 0 || mode
> 9)
289 #ifdef Check_FLT_ROUNDS
290 try_quick
= Rounding
== 1;
294 #endif /*SET_INEXACT*/
301 ilim
= ilim1
= -1; /* Values for cases 0 and 1; done here to */
302 /* silence erroneous "gcc -Wall" warning. */
315 ilim
= ilim1
= i
= ndigits
;
327 s
= s0
= rv_alloc(i
);
329 #ifdef Honor_FLT_ROUNDS
330 if (mode
> 1 && Rounding
!= 1)
334 if (ilim
>= 0 && ilim
<= Quick_max
&& try_quick
) {
336 /* Try to get by with floating-point arithmetic. */
339 dval(&d2
) = dval(&d
);
342 ieps
= 2; /* conservative */
347 /* prevent overflows */
349 dval(&d
) /= bigtens
[n_bigtens
-1];
352 for(; j
; j
>>= 1, i
++)
359 else if (( j1
= -k
)!=0) {
360 dval(&d
) *= tens
[j1
& 0xf];
361 for(j
= j1
>> 4; j
; j
>>= 1, i
++)
364 dval(&d
) *= bigtens
[i
];
367 if (k_check
&& dval(&d
) < 1. && ilim
> 0) {
375 dval(&eps
) = ieps
*dval(&d
) + 7.;
376 word0(&eps
) -= (P
-1)*Exp_msk1
;
380 if (dval(&d
) > dval(&eps
))
382 if (dval(&d
) < -dval(&eps
))
388 /* Use Steele & White method of only
389 * generating digits needed.
391 dval(&eps
) = 0.5/tens
[ilim
-1] - dval(&eps
);
396 if (dval(&d
) < dval(&eps
))
398 if (1. - dval(&d
) < dval(&eps
))
408 /* Generate ilim digits, then fix them up. */
409 dval(&eps
) *= tens
[ilim
-1];
410 for(i
= 1;; i
++, dval(&d
) *= 10.) {
411 L
= (Long
)(dval(&d
));
412 if (!(dval(&d
) -= L
))
416 if (dval(&d
) > 0.5 + dval(&eps
))
418 else if (dval(&d
) < 0.5 - dval(&eps
)) {
431 dval(&d
) = dval(&d2
);
436 /* Do we have a "small" integer? */
438 if (be
>= 0 && k
<= Int_max
) {
441 if (ndigits
< 0 && ilim
<= 0) {
443 if (ilim
< 0 || dval(&d
) <= 5*ds
)
447 for(i
= 1;; i
++, dval(&d
) *= 10.) {
448 L
= (Long
)(dval(&d
) / ds
);
450 #ifdef Check_FLT_ROUNDS
451 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
465 #ifdef Honor_FLT_ROUNDS
469 case 2: goto bump_up
;
472 dval(&d
) += dval(&d
);
476 if (dval(&d
) > ds
|| (dval(&d
) == ds
&& L
& 1))
499 #ifndef Sudden_Underflow
500 denorm
? be
+ (Bias
+ (P
-1) - 1 + 1) :
503 1 + 4*P
- 3 - bbits
+ ((bbits
+ be
- 1) & 3);
511 if (m2
> 0 && s2
> 0) {
512 i
= m2
< s2
? m2
: s2
;
520 mhi
= pow5mult(mhi
, m5
);
525 if (( j
= b5
- m5
)!=0)
535 /* Check for special case that d is a normalized power of 2. */
538 if ((mode
< 2 || leftright
)
539 #ifdef Honor_FLT_ROUNDS
543 if (!word1(&d
) && !(word0(&d
) & Bndry_mask
)
544 #ifndef Sudden_Underflow
545 && word0(&d
) & (Exp_mask
& ~Exp_msk1
)
548 /* The special case */
555 /* Arrange for convenient computation of quotients:
556 * shift left if necessary so divisor has 4 leading 0 bits.
558 * Perhaps we should just compute leading 28 bits of S once
559 * and for all and pass them and a shift to quorem, so it
560 * can do shifts and ors to compute the numerator for q.
563 if (( i
= ((s5
? 32 - hi0bits(S
->x
[S
->wds
-1]) : 1) + s2
) & 0x1f )!=0)
566 if (( i
= ((s5
? 32 - hi0bits(S
->x
[S
->wds
-1]) : 1) + s2
) & 0xf )!=0)
588 b
= multadd(b
, 10, 0); /* we botched the k estimate */
590 mhi
= multadd(mhi
, 10, 0);
594 if (ilim
<= 0 && (mode
== 3 || mode
== 5)) {
595 if (ilim
< 0 || cmp(b
,S
= multadd(S
,5,0)) <= 0) {
596 /* no digits, fcvt style */
608 mhi
= lshift(mhi
, m2
);
610 /* Compute mlo -- check for special case
611 * that d is a normalized power of 2.
616 mhi
= Balloc(mhi
->k
);
618 mhi
= lshift(mhi
, Log2P
);
622 dig
= quorem(b
,S
) + '0';
623 /* Do we yet have the shortest decimal string
624 * that will round to d?
627 delta
= diff(S
, mhi
);
628 j1
= delta
->sign
? 1 : cmp(b
, delta
);
631 if (j1
== 0 && mode
!= 1 && !(word1(&d
) & 1)
632 #ifdef Honor_FLT_ROUNDS
641 else if (!b
->x
[0] && b
->wds
<= 1)
648 if (j
< 0 || (j
== 0 && mode
!= 1
653 if (!b
->x
[0] && b
->wds
<= 1) {
659 #ifdef Honor_FLT_ROUNDS
662 case 0: goto accept_dig
;
663 case 2: goto keep_dig
;
665 #endif /*Honor_FLT_ROUNDS*/
672 if ((j1
> 0 || (j1
== 0 && dig
& 1))
682 #ifdef Honor_FLT_ROUNDS
686 if (dig
== '9') { /* possible if i == 1 */
694 #ifdef Honor_FLT_ROUNDS
700 b
= multadd(b
, 10, 0);
702 mlo
= mhi
= multadd(mhi
, 10, 0);
704 mlo
= multadd(mlo
, 10, 0);
705 mhi
= multadd(mhi
, 10, 0);
711 *s
++ = dig
= quorem(b
,S
) + '0';
712 if (!b
->x
[0] && b
->wds
<= 1) {
720 b
= multadd(b
, 10, 0);
723 /* Round off last digit */
725 #ifdef Honor_FLT_ROUNDS
727 case 0: goto trimzeros
;
728 case 2: goto roundoff
;
736 if (j
> 0 || (j
== 0 && dig
& 1))
749 #ifdef Honor_FLT_ROUNDS
758 if (mlo
&& mlo
!= mhi
)
766 word0(&d
) = Exp_1
+ (70 << Exp_shift
);
771 else if (!oldinexact
)