1 /* $NetBSD: dtoa.c,v 1.5 2008/03/21 23:13:48 christos Exp $ */
3 /****************************************************************
5 The author of this software is David M. Gay.
7 Copyright (C) 1998, 1999 by Lucent Technologies
10 Permission to use, copy, modify, and distribute this software and
11 its documentation for any purpose and without fee is hereby
12 granted, provided that the above copyright notice appear in all
13 copies and that both that the copyright notice and this
14 permission notice and warranty disclaimer appear in supporting
15 documentation, and that the name of Lucent or any of its entities
16 not be used in advertising or publicity pertaining to
17 distribution of the software without specific, written prior
20 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
21 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
22 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
23 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
24 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
25 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
26 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
29 ****************************************************************/
31 /* Please send bug reports to David M. Gay (dmg at acm dot org,
32 * with " at " changed at "@" and " dot " changed to "."). */
36 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
38 * Inspired by "How to Print Floating-Point Numbers Accurately" by
39 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
42 * 1. Rather than iterating, we use a simple numeric overestimate
43 * to determine k = floor(log10(d)). We scale relevant
44 * quantities using O(log2(k)) rather than O(k) multiplications.
45 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
46 * try to generate digits strictly left to right. Instead, we
47 * compute with fewer bits and propagate the carry if necessary
48 * when rounding the final digit up. This is often faster.
49 * 3. Under the assumption that input will be rounded nearest,
50 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
51 * That is, we allow equality in stopping tests when the
52 * round-nearest rule will give the same floating-point value
53 * as would satisfaction of the stopping test with strict
55 * 4. We remove common factors of powers of 2 from relevant
57 * 5. When converting floating-point integers less than 1e16,
58 * we use floating-point arithmetic rather than resorting
59 * to multiple-precision integers.
60 * 6. When asked to produce fewer than 15 digits, we first try
61 * to get by with floating-point arithmetic; we resort to
62 * multiple-precision integer arithmetic only if we cannot
63 * guarantee that the floating-point calculation has given
64 * the correctly rounded result. For k requested digits and
65 * "uniformly" distributed input, the probability is
66 * something like 10^(k-15) that we must resort to the Long
70 #ifdef Honor_FLT_ROUNDS
71 #define Rounding rounding
72 #undef Check_FLT_ROUNDS
73 #define Check_FLT_ROUNDS
75 #define Rounding Flt_Rounds
81 (d
, mode
, ndigits
, decpt
, sign
, rve
)
82 double d
; int mode
, ndigits
, *decpt
, *sign
; char **rve
;
84 (double d
, int mode
, int ndigits
, int *decpt
, int *sign
, char **rve
)
87 /* Arguments ndigits, decpt, sign are similar to those
88 of ecvt and fcvt; trailing zeros are suppressed from
89 the returned string. If not null, *rve is set to point
90 to the end of the return value. If d is +-Infinity or NaN,
91 then *decpt is set to 9999.
94 0 ==> shortest string that yields d when read in
95 and rounded to nearest.
96 1 ==> like 0, but with Steele & White stopping rule;
97 e.g. with IEEE P754 arithmetic , mode 0 gives
98 1e23 whereas mode 1 gives 9.999999999999999e22.
99 2 ==> max(1,ndigits) significant digits. This gives a
100 return value similar to that of ecvt, except
101 that trailing zeros are suppressed.
102 3 ==> through ndigits past the decimal point. This
103 gives a return value similar to that from fcvt,
104 except that trailing zeros are suppressed, and
105 ndigits can be negative.
106 4,5 ==> similar to 2 and 3, respectively, but (in
107 round-nearest mode) with the tests of mode 0 to
108 possibly return a shorter string that rounds to d.
109 With IEEE arithmetic and compilation with
110 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
111 as modes 2 and 3 when FLT_ROUNDS != 1.
112 6-9 ==> Debugging modes similar to mode - 4: don't try
113 fast floating-point estimate (if applicable).
115 Values of mode other than 0-9 are treated as mode 0.
117 Sufficient space is allocated to the return value
118 to hold the suppressed trailing zeros.
121 int bbits
, b2
, b5
, be
, dig
, i
, ieps
, ilim0
,
122 j
, jj1
, k
, k0
, k_check
, leftright
, m2
, m5
, s2
, s5
,
123 spec_case
, try_quick
;
124 int ilim
= 0, ilim1
= 0; /* pacify gcc */
126 #ifndef Sudden_Underflow
130 Bigint
*b
, *b1
, *delta
, *mhi
, *S
;
131 Bigint
*mlo
= NULL
; /* pacify gcc */
134 #ifdef Honor_FLT_ROUNDS
138 int inexact
, oldinexact
;
141 #ifndef MULTIPLE_THREADS
143 freedtoa(dtoa_result
);
148 if (word0(d
) & Sign_bit
) {
149 /* set sign for everything, including 0's and NaNs */
151 word0(d
) &= ~Sign_bit
; /* clear sign bit */
156 #if defined(IEEE_Arith) + defined(VAX)
158 if ((word0(d
) & Exp_mask
) == Exp_mask
)
160 if (word0(d
) == 0x8000)
163 /* Infinity or NaN */
166 if (!word1(d
) && !(word0(d
) & 0xfffff))
167 return nrv_alloc("Infinity", rve
, 8);
169 return nrv_alloc("NaN", rve
, 3);
173 dval(d
) += 0; /* normalize */
177 return nrv_alloc("0", rve
, 1);
181 try_quick
= oldinexact
= get_inexact();
184 #ifdef Honor_FLT_ROUNDS
185 if ((rounding
= Flt_Rounds
) >= 2) {
187 rounding
= rounding
== 2 ? 0 : 2;
194 b
= d2b(dval(d
), &be
, &bbits
);
197 #ifdef Sudden_Underflow
198 i
= (int)(word0(d
) >> Exp_shift1
& (Exp_mask
>>Exp_shift1
));
200 if (( i
= (int)(word0(d
) >> Exp_shift1
& (Exp_mask
>>Exp_shift1
)) )!=0) {
203 word0(d2
) &= Frac_mask1
;
206 if (( j
= 11 - hi0bits(word0(d2
) & Frac_mask
) )!=0)
210 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
211 * log10(x) = log(x) / log(10)
212 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
213 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
215 * This suggests computing an approximation k to log10(d) by
217 * k = (i - Bias)*0.301029995663981
218 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
220 * We want k to be too large rather than too small.
221 * The error in the first-order Taylor series approximation
222 * is in our favor, so we just round up the constant enough
223 * to compensate for any error in the multiplication of
224 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
225 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
226 * adding 1e-13 to the constant term more than suffices.
227 * Hence we adjust the constant term to 0.1760912590558.
228 * (We could get a more accurate k by invoking log10,
229 * but this is probably not worthwhile.)
237 #ifndef Sudden_Underflow
241 /* d is denormalized */
243 i
= bbits
+ be
+ (Bias
+ (P
-1) - 1);
244 x
= i
> 32 ? word0(d
) << (64 - i
) | word1(d
) >> (i
- 32)
245 : word1(d
) << (32 - i
);
247 word0(d2
) -= 31*Exp_msk1
; /* adjust exponent */
248 i
-= (Bias
+ (P
-1) - 1) + 1;
252 ds
= (dval(d2
)-1.5)*0.289529654602168 + 0.1760912590558 + i
*0.301029995663981;
254 if (ds
< 0. && ds
!= k
)
255 k
--; /* want k = floor(ds) */
257 if (k
>= 0 && k
<= Ten_pmax
) {
258 if (dval(d
) < tens
[k
])
281 if (mode
< 0 || mode
> 9)
285 #ifdef Check_FLT_ROUNDS
286 try_quick
= Rounding
== 1;
290 #endif /*SET_INEXACT*/
310 ilim
= ilim1
= i
= ndigits
;
322 s
= s0
= rv_alloc((size_t)i
);
326 #ifdef Honor_FLT_ROUNDS
327 if (mode
> 1 && rounding
!= 1)
331 if (ilim
>= 0 && ilim
<= Quick_max
&& try_quick
) {
333 /* Try to get by with floating-point arithmetic. */
339 ieps
= 2; /* conservative */
342 j
= (unsigned int)k
>> 4;
344 /* prevent overflows */
346 dval(d
) /= bigtens
[n_bigtens
-1];
349 for(; j
; j
= (unsigned int)j
>> 1, i
++)
356 else if (( jj1
= -k
)!=0) {
357 dval(d
) *= tens
[jj1
& 0xf];
358 for(j
= jj1
>> 4; j
; j
>>= 1, i
++)
361 dval(d
) *= bigtens
[i
];
364 if (k_check
&& dval(d
) < 1. && ilim
> 0) {
372 dval(eps
) = ieps
*dval(d
) + 7.;
373 word0(eps
) -= (P
-1)*Exp_msk1
;
377 if (dval(d
) > dval(eps
))
379 if (dval(d
) < -dval(eps
))
385 /* Use Steele & White method of only
386 * generating digits needed.
388 dval(eps
) = 0.5/tens
[ilim
-1] - dval(eps
);
393 if (dval(d
) < dval(eps
))
395 if (1. - dval(d
) < dval(eps
))
405 /* Generate ilim digits, then fix them up. */
406 dval(eps
) *= tens
[ilim
-1];
407 for(i
= 1;; i
++, dval(d
) *= 10.) {
413 if (dval(d
) > 0.5 + dval(eps
))
415 else if (dval(d
) < 0.5 - dval(eps
)) {
433 /* Do we have a "small" integer? */
435 if (be
>= 0 && k
<= Int_max
) {
438 if (ndigits
< 0 && ilim
<= 0) {
440 if (ilim
< 0 || dval(d
) <= 5*ds
)
444 for(i
= 1;; i
++, dval(d
) *= 10.) {
445 L
= (Long
)(dval(d
) / ds
);
447 #ifdef Check_FLT_ROUNDS
448 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
462 #ifdef Honor_FLT_ROUNDS
466 case 2: goto bump_up
;
470 if (dval(d
) > ds
|| (dval(d
) == ds
&& L
& 1)) {
491 #ifndef Sudden_Underflow
492 denorm
? be
+ (Bias
+ (P
-1) - 1 + 1) :
495 1 + 4*P
- 3 - bbits
+ ((bbits
+ be
- 1) & 3);
505 if (m2
> 0 && s2
> 0) {
506 i
= m2
< s2
? m2
: s2
;
514 mhi
= pow5mult(mhi
, m5
);
523 if (( j
= b5
- m5
)!=0)
542 /* Check for special case that d is a normalized power of 2. */
545 if ((mode
< 2 || leftright
)
546 #ifdef Honor_FLT_ROUNDS
550 if (!word1(d
) && !(word0(d
) & Bndry_mask
)
551 #ifndef Sudden_Underflow
552 && word0(d
) & (Exp_mask
& ~Exp_msk1
)
555 /* The special case */
562 /* Arrange for convenient computation of quotients:
563 * shift left if necessary so divisor has 4 leading 0 bits.
565 * Perhaps we should just compute leading 28 bits of S once
566 * and for all and pass them and a shift to quorem, so it
567 * can do shifts and ors to compute the numerator for q.
570 if (( i
= ((s5
? 32 - hi0bits(S
->x
[S
->wds
-1]) : 1) + s2
) & 0x1f )!=0)
573 if (( i
= ((s5
? 32 - hi0bits(S
->x
[S
->wds
-1]) : 1) + s2
) & 0xf )!=0)
601 b
= multadd(b
, 10, 0); /* we botched the k estimate */
605 mhi
= multadd(mhi
, 10, 0);
612 if (ilim
<= 0 && (mode
== 3 || mode
== 5)) {
613 if (ilim
< 0 || cmp(b
,S
= multadd(S
,5,0)) <= 0) {
614 /* no digits, fcvt style */
626 mhi
= lshift(mhi
, m2
);
631 /* Compute mlo -- check for special case
632 * that d is a normalized power of 2.
637 mhi
= Balloc(mhi
->k
);
641 mhi
= lshift(mhi
, Log2P
);
647 dig
= quorem(b
,S
) + '0';
648 /* Do we yet have the shortest decimal string
649 * that will round to d?
652 delta
= diff(S
, mhi
);
655 jj1
= delta
->sign
? 1 : cmp(b
, delta
);
658 if (jj1
== 0 && mode
!= 1 && !(word1(d
) & 1)
659 #ifdef Honor_FLT_ROUNDS
668 else if (!b
->x
[0] && b
->wds
<= 1)
675 if (j
< 0 || (j
== 0 && mode
!= 1
680 if (!b
->x
[0] && b
->wds
<= 1) {
686 #ifdef Honor_FLT_ROUNDS
689 case 0: goto accept_dig
;
690 case 2: goto keep_dig
;
692 #endif /*Honor_FLT_ROUNDS*/
698 if ((jj1
> 0 || (jj1
== 0 && dig
& 1))
707 #ifdef Honor_FLT_ROUNDS
711 if (dig
== '9') { /* possible if i == 1 */
719 #ifdef Honor_FLT_ROUNDS
725 b
= multadd(b
, 10, 0);
729 mlo
= mhi
= multadd(mhi
, 10, 0);
734 mlo
= multadd(mlo
, 10, 0);
737 mhi
= multadd(mhi
, 10, 0);
745 *s
++ = dig
= quorem(b
,S
) + '0';
746 if (!b
->x
[0] && b
->wds
<= 1) {
754 b
= multadd(b
, 10, 0);
759 /* Round off last digit */
761 #ifdef Honor_FLT_ROUNDS
763 case 0: goto trimzeros
;
764 case 2: goto roundoff
;
769 if (j
> 0 || (j
== 0 && dig
& 1)) {
780 #ifdef Honor_FLT_ROUNDS
789 if (mlo
&& mlo
!= mhi
)
797 word0(d
) = Exp_1
+ (70 << Exp_shift
);
802 else if (!oldinexact
)
806 if (s
== s0
) { /* don't return empty string */