]>
git.proxmox.com Git - mirror_edk2.git/blob - StdLib/LibC/gdtoa/dtoa.c
1 /* $NetBSD: dtoa.c,v 1.3.4.1.4.1 2008/04/08 21:10:55 jdc 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 "."). */
33 #include <LibConfig.h>
37 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
39 * Inspired by "How to Print Floating-Point Numbers Accurately" by
40 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
43 * 1. Rather than iterating, we use a simple numeric overestimate
44 * to determine k = floor(log10(d)). We scale relevant
45 * quantities using O(log2(k)) rather than O(k) multiplications.
46 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
47 * try to generate digits strictly left to right. Instead, we
48 * compute with fewer bits and propagate the carry if necessary
49 * when rounding the final digit up. This is often faster.
50 * 3. Under the assumption that input will be rounded nearest,
51 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
52 * That is, we allow equality in stopping tests when the
53 * round-nearest rule will give the same floating-point value
54 * as would satisfaction of the stopping test with strict
56 * 4. We remove common factors of powers of 2 from relevant
58 * 5. When converting floating-point integers less than 1e16,
59 * we use floating-point arithmetic rather than resorting
60 * to multiple-precision integers.
61 * 6. When asked to produce fewer than 15 digits, we first try
62 * to get by with floating-point arithmetic; we resort to
63 * multiple-precision integer arithmetic only if we cannot
64 * guarantee that the floating-point calculation has given
65 * the correctly rounded result. For k requested digits and
66 * "uniformly" distributed input, the probability is
67 * something like 10^(k-15) that we must resort to the Long
71 #ifdef Honor_FLT_ROUNDS
72 #define Rounding rounding
73 #undef Check_FLT_ROUNDS
74 #define Check_FLT_ROUNDS
76 #define Rounding Flt_Rounds
79 #if defined(_MSC_VER) /* Handle Microsoft VC++ compiler specifics. */
80 // Disable: warning C4700: uninitialized local variable 'xx' used
81 #pragma warning ( disable : 4700 )
82 #endif /* defined(_MSC_VER) */
87 (d
, mode
, ndigits
, decpt
, sign
, rve
)
88 double d
; int mode
, ndigits
, *decpt
, *sign
; char **rve
;
90 (double d
, int mode
, int ndigits
, int *decpt
, int *sign
, char **rve
)
93 /* Arguments ndigits, decpt, sign are similar to those
94 of ecvt and fcvt; trailing zeros are suppressed from
95 the returned string. If not null, *rve is set to point
96 to the end of the return value. If d is +-Infinity or NaN,
97 then *decpt is set to 9999.
100 0 ==> shortest string that yields d when read in
101 and rounded to nearest.
102 1 ==> like 0, but with Steele & White stopping rule;
103 e.g. with IEEE P754 arithmetic , mode 0 gives
104 1e23 whereas mode 1 gives 9.999999999999999e22.
105 2 ==> max(1,ndigits) significant digits. This gives a
106 return value similar to that of ecvt, except
107 that trailing zeros are suppressed.
108 3 ==> through ndigits past the decimal point. This
109 gives a return value similar to that from fcvt,
110 except that trailing zeros are suppressed, and
111 ndigits can be negative.
112 4,5 ==> similar to 2 and 3, respectively, but (in
113 round-nearest mode) with the tests of mode 0 to
114 possibly return a shorter string that rounds to d.
115 With IEEE arithmetic and compilation with
116 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
117 as modes 2 and 3 when FLT_ROUNDS != 1.
118 6-9 ==> Debugging modes similar to mode - 4: don't try
119 fast floating-point estimate (if applicable).
121 Values of mode other than 0-9 are treated as mode 0.
123 Sufficient space is allocated to the return value
124 to hold the suppressed trailing zeros.
127 int bbits
, b2
, b5
, be
, dig
, i
, ieps
, ilim0
,
128 j
, jj1
, k
, k0
, k_check
, leftright
, m2
, m5
, s2
, s5
,
129 spec_case
, try_quick
;
130 int ilim
= 0, ilim1
= 0; /* pacify gcc */
132 #ifndef Sudden_Underflow
136 Bigint
*b
, *b1
, *delta
, *mhi
, *S
;
137 Bigint
*mlo
= NULL
; /* pacify gcc */
140 #ifdef Honor_FLT_ROUNDS
144 int inexact
, oldinexact
;
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
191 if ((rounding
= Flt_Rounds
) >= 2) {
193 rounding
= rounding
== 2 ? 0 : 2;
200 b
= d2b(dval(d
), &be
, &bbits
);
203 #ifdef Sudden_Underflow
204 i
= (int)(word0(d
) >> Exp_shift1
& (Exp_mask
>>Exp_shift1
));
206 if (( i
= (int)(word0(d
) >> Exp_shift1
& (Exp_mask
>>Exp_shift1
)) )!=0) {
209 word0(d2
) &= Frac_mask1
;
212 if (( j
= 11 - hi0bits(word0(d2
) & Frac_mask
) )!=0)
216 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
217 * log10(x) = log(x) / log(10)
218 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
219 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
221 * This suggests computing an approximation k to log10(d) by
223 * k = (i - Bias)*0.301029995663981
224 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
226 * We want k to be too large rather than too small.
227 * The error in the first-order Taylor series approximation
228 * is in our favor, so we just round up the constant enough
229 * to compensate for any error in the multiplication of
230 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
231 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
232 * adding 1e-13 to the constant term more than suffices.
233 * Hence we adjust the constant term to 0.1760912590558.
234 * (We could get a more accurate k by invoking log10,
235 * but this is probably not worthwhile.)
243 #ifndef Sudden_Underflow
247 /* d is denormalized */
249 i
= bbits
+ be
+ (Bias
+ (P
-1) - 1);
250 x
= i
> 32 ? word0(d
) << (64 - i
) | word1(d
) >> (i
- 32)
251 : word1(d
) << (32 - i
);
252 dval(d2
) = (double)x
;
253 word0(d2
) -= 31*Exp_msk1
; /* adjust exponent */
254 i
-= (Bias
+ (P
-1) - 1) + 1;
258 ds
= (dval(d2
)-1.5)*0.289529654602168 + 0.1760912590558 + i
*0.301029995663981;
260 if (ds
< 0. && ds
!= k
)
261 k
--; /* want k = floor(ds) */
263 if (k
>= 0 && k
<= Ten_pmax
) {
264 if (dval(d
) < tens
[k
])
287 if (mode
< 0 || mode
> 9)
291 #ifdef Check_FLT_ROUNDS
292 try_quick
= Rounding
== 1;
296 #endif /*SET_INEXACT*/
316 ilim
= ilim1
= i
= ndigits
;
328 s
= s0
= rv_alloc((size_t)i
);
332 #ifdef Honor_FLT_ROUNDS
333 if (mode
> 1 && rounding
!= 1)
337 if (ilim
>= 0 && ilim
<= Quick_max
&& try_quick
) {
339 /* Try to get by with floating-point arithmetic. */
345 ieps
= 2; /* conservative */
348 j
= (unsigned int)k
>> 4;
350 /* prevent overflows */
352 dval(d
) /= bigtens
[n_bigtens
-1];
355 for(; j
; j
= (unsigned int)j
>> 1, i
++)
362 else if (( jj1
= -k
)!=0) {
363 dval(d
) *= tens
[jj1
& 0xf];
364 for(j
= jj1
>> 4; j
; j
>>= 1, i
++)
367 dval(d
) *= bigtens
[i
];
370 if (k_check
&& dval(d
) < 1. && ilim
> 0) {
378 dval(eps
) = ieps
*dval(d
) + 7.;
379 word0(eps
) -= (P
-1)*Exp_msk1
;
383 if (dval(d
) > dval(eps
))
385 if (dval(d
) < -dval(eps
))
391 /* Use Steele & White method of only
392 * generating digits needed.
394 dval(eps
) = 0.5/tens
[ilim
-1] - dval(eps
);
398 *s
++ = (char)('0' + (int)L
);
399 if (dval(d
) < dval(eps
))
401 if (1. - dval(d
) < dval(eps
))
411 /* Generate ilim digits, then fix them up. */
412 dval(eps
) *= tens
[ilim
-1];
413 for(i
= 1;; i
++, dval(d
) *= 10.) {
417 *s
++ = (char)('0' + (int)L
);
419 if (dval(d
) > 0.5 + dval(eps
))
421 else if (dval(d
) < 0.5 - dval(eps
)) {
439 /* Do we have a "small" integer? */
441 if (be
>= 0 && k
<= Int_max
) {
444 if (ndigits
< 0 && ilim
<= 0) {
446 if (ilim
< 0 || dval(d
) <= 5*ds
)
450 for(i
= 1;; i
++, dval(d
) *= 10.) {
451 L
= (Long
)(dval(d
) / ds
);
453 #ifdef Check_FLT_ROUNDS
454 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
460 *s
++ = (char)('0' + (int)L
);
468 #ifdef Honor_FLT_ROUNDS
472 case 2: goto bump_up
;
476 if (dval(d
) > ds
|| (dval(d
) == ds
&& L
& 1)) {
497 #ifndef Sudden_Underflow
498 denorm
? be
+ (Bias
+ (P
-1) - 1 + 1) :
501 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
);
529 if (( j
= b5
- m5
)!=0)
548 /* Check for special case that d is a normalized power of 2. */
551 if ((mode
< 2 || leftright
)
552 #ifdef Honor_FLT_ROUNDS
556 if (!word1(d
) && !(word0(d
) & Bndry_mask
)
557 #ifndef Sudden_Underflow
558 && word0(d
) & (Exp_mask
& ~Exp_msk1
)
561 /* The special case */
568 /* Arrange for convenient computation of quotients:
569 * shift left if necessary so divisor has 4 leading 0 bits.
571 * Perhaps we should just compute leading 28 bits of S once
572 * and for all and pass them and a shift to quorem, so it
573 * can do shifts and ors to compute the numerator for q.
576 if (( i
= ((s5
? 32 - hi0bits(S
->x
[S
->wds
-1]) : 1) + s2
) & 0x1f )!=0)
579 if (( i
= ((s5
? 32 - hi0bits(S
->x
[S
->wds
-1]) : 1) + s2
) & 0xf )!=0)
607 b
= multadd(b
, 10, 0); /* we botched the k estimate */
611 mhi
= multadd(mhi
, 10, 0);
618 if (ilim
<= 0 && (mode
== 3 || mode
== 5)) {
619 if (ilim
< 0 || cmp(b
,S
= multadd(S
,5,0)) <= 0) {
620 /* no digits, fcvt style */
632 mhi
= lshift(mhi
, m2
);
637 /* Compute mlo -- check for special case
638 * that d is a normalized power of 2.
643 mhi
= Balloc(mhi
->k
);
647 mhi
= lshift(mhi
, Log2P
);
653 dig
= quorem(b
,S
) + '0';
654 /* Do we yet have the shortest decimal string
655 * that will round to d?
658 delta
= diff(S
, mhi
);
661 jj1
= delta
->sign
? 1 : cmp(b
, delta
);
664 if (jj1
== 0 && mode
!= 1 && !(word1(d
) & 1)
665 #ifdef Honor_FLT_ROUNDS
674 else if (!b
->x
[0] && b
->wds
<= 1)
681 if (j
< 0 || (j
== 0 && mode
!= 1
686 if (!b
->x
[0] && b
->wds
<= 1) {
692 #ifdef Honor_FLT_ROUNDS
695 case 0: goto accept_dig
;
696 case 2: goto keep_dig
;
698 #endif /*Honor_FLT_ROUNDS*/
704 if ((jj1
> 0 || (jj1
== 0 && dig
& 1))
713 #ifdef Honor_FLT_ROUNDS
717 if (dig
== '9') { /* possible if i == 1 */
722 *s
++ = (char)(dig
+ 1);
725 #ifdef Honor_FLT_ROUNDS
731 b
= multadd(b
, 10, 0);
735 mlo
= mhi
= multadd(mhi
, 10, 0);
740 mlo
= multadd(mlo
, 10, 0);
743 mhi
= multadd(mhi
, 10, 0);
751 *s
++ = (char)(dig
= (int)(quorem(b
,S
) + '0'));
752 if (!b
->x
[0] && b
->wds
<= 1) {
760 b
= multadd(b
, 10, 0);
765 /* Round off last digit */
767 #ifdef Honor_FLT_ROUNDS
769 case 0: goto trimzeros
;
770 case 2: goto roundoff
;
775 if (j
> 0 || (j
== 0 && dig
& 1)) {
786 #ifdef Honor_FLT_ROUNDS
795 if (mlo
&& mlo
!= mhi
)
803 word0(d
) = Exp_1
+ (70 << Exp_shift
);
808 else if (!oldinexact
)
812 if (s
== s0
) { /* don't return empty string */