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 "."). */
36 bitstob(bits
, nbits
, bbits
) ULong
*bits
; int nbits
; int *bbits
;
38 bitstob(ULong
*bits
, int nbits
, int *bbits
)
56 be
= bits
+ ((nbits
- 1) >> kshift
);
59 *x
++ = *bits
& ALL_ON
;
61 *x
++ = (*bits
>> 16) & ALL_ON
;
63 } while(++bits
<= be
);
72 *bbits
= i
*ULbits
+ 32 - hi0bits(b
->x
[i
]);
77 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
79 * Inspired by "How to Print Floating-Point Numbers Accurately" by
80 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
83 * 1. Rather than iterating, we use a simple numeric overestimate
84 * to determine k = floor(log10(d)). We scale relevant
85 * quantities using O(log2(k)) rather than O(k) multiplications.
86 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
87 * try to generate digits strictly left to right. Instead, we
88 * compute with fewer bits and propagate the carry if necessary
89 * when rounding the final digit up. This is often faster.
90 * 3. Under the assumption that input will be rounded nearest,
91 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
92 * That is, we allow equality in stopping tests when the
93 * round-nearest rule will give the same floating-point value
94 * as would satisfaction of the stopping test with strict
96 * 4. We remove common factors of powers of 2 from relevant
98 * 5. When converting floating-point integers less than 1e16,
99 * we use floating-point arithmetic rather than resorting
100 * to multiple-precision integers.
101 * 6. When asked to produce fewer than 15 digits, we first try
102 * to get by with floating-point arithmetic; we resort to
103 * multiple-precision integer arithmetic only if we cannot
104 * guarantee that the floating-point calculation has given
105 * the correctly rounded result. For k requested digits and
106 * "uniformly" distributed input, the probability is
107 * something like 10^(k-15) that we must resort to the Long
114 (fpi
, be
, bits
, kindp
, mode
, ndigits
, decpt
, rve
)
115 FPI
*fpi
; int be
; ULong
*bits
;
116 int *kindp
, mode
, ndigits
, *decpt
; char **rve
;
118 (FPI
*fpi
, int be
, ULong
*bits
, int *kindp
, int mode
, int ndigits
, int *decpt
, char **rve
)
121 /* Arguments ndigits and decpt are similar to the second and third
122 arguments of ecvt and fcvt; trailing zeros are suppressed from
123 the returned string. If not null, *rve is set to point
124 to the end of the return value. If d is +-Infinity or NaN,
125 then *decpt is set to 9999.
126 be = exponent: value = (integer represented by bits) * (2 to the power of be).
129 0 ==> shortest string that yields d when read in
130 and rounded to nearest.
131 1 ==> like 0, but with Steele & White stopping rule;
132 e.g. with IEEE P754 arithmetic , mode 0 gives
133 1e23 whereas mode 1 gives 9.999999999999999e22.
134 2 ==> max(1,ndigits) significant digits. This gives a
135 return value similar to that of ecvt, except
136 that trailing zeros are suppressed.
137 3 ==> through ndigits past the decimal point. This
138 gives a return value similar to that from fcvt,
139 except that trailing zeros are suppressed, and
140 ndigits can be negative.
141 4-9 should give the same return values as 2-3, i.e.,
142 4 <= mode <= 9 ==> same return as mode
143 2 + (mode & 1). These modes are mainly for
144 debugging; often they run slower but sometimes
145 faster than modes 2-3.
146 4,5,8,9 ==> left-to-right digit generation.
147 6-9 ==> don't try fast floating-point estimate
150 Values of mode other than 0-9 are treated as mode 0.
152 Sufficient space is allocated to the return value
153 to hold the suppressed trailing zeros.
156 int bbits
, b2
, b5
, be0
, dig
, i
, ieps
, ilim
, ilim0
, ilim1
, inex
;
157 int j
, j1
, k
, k0
, k_check
, kind
, leftright
, m2
, m5
, nbits
;
158 int rdir
, s2
, s5
, spec_case
, try_quick
;
160 Bigint
*b
, *b1
, *delta
, *mlo
, *mhi
, *mhi1
, *S
;
165 #ifndef MULTIPLE_THREADS
167 freedtoa(dtoa_result
);
172 kind
= *kindp
&= ~STRTOG_Inexact
;
173 switch(kind
& STRTOG_Retmask
) {
177 case STRTOG_Denormal
:
179 case STRTOG_Infinite
:
181 return nrv_alloc("Infinity", rve
, 8);
184 return nrv_alloc("NaN", rve
, 3);
188 b
= bitstob(bits
, nbits
= fpi
->nbits
, &bbits
);
190 if ( (i
= trailz(b
)) !=0) {
199 return nrv_alloc("0", rve
, 1);
202 dval(&d
) = b2d(b
, &i
);
204 word0(&d
) &= Frac_mask1
;
207 if ( (j
= 11 - hi0bits(word0(&d
) & Frac_mask
)) !=0)
211 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
212 * log10(x) = log(x) / log(10)
213 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
214 * log10(&d) = (i-Bias)*log(2)/log(10) + log10(d2)
216 * This suggests computing an approximation k to log10(&d) by
218 * k = (i - Bias)*0.301029995663981
219 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
221 * We want k to be too large rather than too small.
222 * The error in the first-order Taylor series approximation
223 * is in our favor, so we just round up the constant enough
224 * to compensate for any error in the multiplication of
225 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
226 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
227 * adding 1e-13 to the constant term more than suffices.
228 * Hence we adjust the constant term to 0.1760912590558.
229 * (We could get a more accurate k by invoking log10,
230 * but this is probably not worthwhile.)
236 ds
= (dval(&d
)-1.5)*0.289529654602168 + 0.1760912590558 + i
*0.301029995663981;
238 /* correct assumption about exponent range */
245 if (ds
< 0. && ds
!= k
)
246 k
--; /* want k = floor(ds) */
250 if ( (j1
= j
& 3) !=0)
252 word0(&d
) += j
<< Exp_shift
- 2 & Exp_mask
;
254 word0(&d
) += (be
+ bbits
- 1) << Exp_shift
;
256 if (k
>= 0 && k
<= Ten_pmax
) {
257 if (dval(&d
) < tens
[k
])
280 if (mode
< 0 || mode
> 9)
287 else if (i
>= -4 - Emin
|| i
< Emin
)
290 ilim
= ilim1
= -1; /* Values for cases 0 and 1; done here to */
291 /* silence erroneous "gcc -Wall" warning. */
295 i
= (int)(nbits
* .30103) + 3;
304 ilim
= ilim1
= i
= ndigits
;
316 s
= s0
= rv_alloc(i
);
318 if ( (rdir
= fpi
->rounding
- 1) !=0) {
321 if (kind
& STRTOG_Neg
)
325 /* Now rdir = 0 ==> round near, 1 ==> round up, 2 ==> round down. */
327 if (ilim
>= 0 && ilim
<= Quick_max
&& try_quick
&& !rdir
328 #ifndef IMPRECISE_INEXACT
333 /* Try to get by with floating-point arithmetic. */
338 if ( (j
= 11 - hi0bits(word0(&d
) & Frac_mask
)) !=0)
343 ieps
= 2; /* conservative */
348 /* prevent overflows */
350 dval(&d
) /= bigtens
[n_bigtens
-1];
353 for(; j
; j
>>= 1, i
++)
361 if ( (j1
= -k
) !=0) {
362 dval(&d
) *= tens
[j1
& 0xf];
363 for(j
= j1
>> 4; j
; j
>>= 1, i
++)
366 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
) = ds
*0.5/tens
[ilim
-1] - dval(&eps
);
396 L
= (Long
)(dval(&d
)/ds
);
399 if (dval(&d
) < dval(&eps
)) {
401 inex
= STRTOG_Inexlo
;
404 if (ds
- dval(&d
) < dval(&eps
))
414 /* Generate ilim digits, then fix them up. */
415 dval(&eps
) *= tens
[ilim
-1];
416 for(i
= 1;; i
++, dval(&d
) *= 10.) {
417 if ( (L
= (Long
)(dval(&d
)/ds
)) !=0)
422 if (dval(&d
) > ds
+ dval(&eps
))
424 else if (dval(&d
) < ds
- dval(&eps
)) {
426 inex
= STRTOG_Inexlo
;
427 goto clear_trailing0
;
442 /* Do we have a "small" integer? */
444 if (be
>= 0 && k
<= Int_max
) {
447 if (ndigits
< 0 && ilim
<= 0) {
449 if (ilim
< 0 || dval(&d
) <= 5*ds
)
453 for(i
= 1;; i
++, dval(&d
) *= 10.) {
456 #ifdef Check_FLT_ROUNDS
457 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
470 inex
= STRTOG_Inexlo
;
473 dval(&d
) += dval(&d
);
477 if (dval(&d
) > ds
|| (dval(&d
) == ds
&& L
& 1))
481 inex
= STRTOG_Inexhi
;
491 inex
= STRTOG_Inexlo
;
507 if (be
- i
++ < fpi
->emin
&& mode
!= 3 && mode
!= 5) {
509 i
= be
- fpi
->emin
+ 1;
510 if (mode
>= 2 && ilim
> 0 && ilim
< i
)
513 else if (mode
>= 2) {
523 if ((i
= ilim
) < 0) {
532 if (m2
> 0 && s2
> 0) {
533 i
= m2
< s2
? m2
: s2
;
541 mhi
= pow5mult(mhi
, m5
);
546 if ( (j
= b5
- m5
) !=0)
556 /* Check for special case that d is a normalized power of 2. */
560 if (bbits
== 1 && be0
> fpi
->emin
+ 1) {
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.
575 i
= ((s5
? hi0bits(S
->x
[S
->wds
-1]) : ULbits
- 1) - s2
- 4) & kmask
;
584 b
= multadd(b
, 10, 0); /* we botched the k estimate */
586 mhi
= multadd(mhi
, 10, 0);
590 if (ilim
<= 0 && mode
> 2) {
591 if (ilim
< 0 || cmp(b
,S
= multadd(S
,5,0)) <= 0) {
592 /* no digits, fcvt style */
595 inex
= STRTOG_Inexlo
;
599 inex
= STRTOG_Inexhi
;
606 mhi
= lshift(mhi
, m2
);
608 /* Compute mlo -- check for special case
609 * that d is a normalized power of 2.
614 mhi
= Balloc(mhi
->k
);
616 mhi
= lshift(mhi
, 1);
620 dig
= quorem(b
,S
) + '0';
621 /* Do we yet have the shortest decimal string
622 * that will round to d?
625 delta
= diff(S
, mhi
);
626 j1
= delta
->sign
? 1 : cmp(b
, delta
);
629 if (j1
== 0 && !mode
&& !(bits
[0] & 1) && !rdir
) {
633 if (b
->wds
> 1 || b
->x
[0])
634 inex
= STRTOG_Inexlo
;
638 inex
= STRTOG_Inexhi
;
644 if (j
< 0 || (j
== 0 && !mode
649 if (rdir
&& (b
->wds
> 1 || b
->x
[0])) {
651 inex
= STRTOG_Inexlo
;
654 while (cmp(S
,mhi
) > 0) {
656 mhi1
= multadd(mhi
, 10, 0);
660 b
= multadd(b
, 10, 0);
661 dig
= quorem(b
,S
) + '0';
665 inex
= STRTOG_Inexhi
;
674 if ((j1
> 0 || (j1
== 0 && dig
& 1))
678 inex
= STRTOG_Inexhi
;
680 if (b
->wds
> 1 || b
->x
[0])
681 inex
= STRTOG_Inexlo
;
686 if (j1
> 0 && rdir
!= 2) {
687 if (dig
== '9') { /* possible if i == 1 */
690 inex
= STRTOG_Inexhi
;
693 inex
= STRTOG_Inexhi
;
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';
714 b
= multadd(b
, 10, 0);
717 /* Round off last digit */
720 if (rdir
== 2 || (b
->wds
<= 1 && !b
->x
[0]))
729 if (j
> 0 || (j
== 0 && dig
& 1))
733 inex
= STRTOG_Inexhi
;
744 if (b
->wds
> 1 || b
->x
[0])
745 inex
= STRTOG_Inexlo
;
752 if (mlo
&& mlo
!= mhi
)