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.
128 0 ==> shortest string that yields d when read in
129 and rounded to nearest.
130 1 ==> like 0, but with Steele & White stopping rule;
131 e.g. with IEEE P754 arithmetic , mode 0 gives
132 1e23 whereas mode 1 gives 9.999999999999999e22.
133 2 ==> max(1,ndigits) significant digits. This gives a
134 return value similar to that of ecvt, except
135 that trailing zeros are suppressed.
136 3 ==> through ndigits past the decimal point. This
137 gives a return value similar to that from fcvt,
138 except that trailing zeros are suppressed, and
139 ndigits can be negative.
140 4-9 should give the same return values as 2-3, i.e.,
141 4 <= mode <= 9 ==> same return as mode
142 2 + (mode & 1). These modes are mainly for
143 debugging; often they run slower but sometimes
144 faster than modes 2-3.
145 4,5,8,9 ==> left-to-right digit generation.
146 6-9 ==> don't try fast floating-point estimate
149 Values of mode other than 0-9 are treated as mode 0.
151 Sufficient space is allocated to the return value
152 to hold the suppressed trailing zeros.
155 int bbits, b2, b5, be0, dig, i, ieps, ilim, ilim0, ilim1, inex;
156 int j, j1, k, k0, k_check, kind, leftright, m2, m5, nbits;
157 int rdir, s2, s5, spec_case, try_quick;
159 Bigint *b, *b1, *delta, *mlo, *mhi, *mhi1, *S;
164 #ifndef MULTIPLE_THREADS
166 freedtoa(dtoa_result);
171 kind = *kindp &= ~STRTOG_Inexact;
172 switch(kind & STRTOG_Retmask) {
176 case STRTOG_Denormal:
178 case STRTOG_Infinite:
180 return nrv_alloc("Infinity", rve, 8);
183 return nrv_alloc("NaN", rve, 3);
187 b = bitstob(bits, nbits = fpi->nbits, &bbits);
189 if ( (i = trailz(b)) !=0) {
198 return nrv_alloc("0", rve, 1);
201 dval(&d) = b2d(b, &i);
203 word0(&d) &= Frac_mask1;
206 if ( (j = 11 - hi0bits(word0(&d) & 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.)
235 ds = (dval(&d)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
237 /* correct assumption about exponent range */
244 if (ds < 0. && ds != k)
245 k--; /* want k = floor(ds) */
249 if ( (j1 = j & 3) !=0)
251 word0(&d) += j << Exp_shift - 2 & Exp_mask;
253 word0(&d) += (be + bbits - 1) << Exp_shift;
255 if (k >= 0 && k <= Ten_pmax) {
256 if (dval(&d) < tens[k])
279 if (mode < 0 || mode > 9)
287 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
288 /* silence erroneous "gcc -Wall" warning. */
292 i = (int)(nbits * .30103) + 3;
301 ilim = ilim1 = i = ndigits;
313 s = s0 = rv_alloc(i);
315 if ( (rdir = fpi->rounding - 1) !=0) {
318 if (kind & STRTOG_Neg)
322 /* Now rdir = 0 ==> round near, 1 ==> round up, 2 ==> round down. */
324 if (ilim >= 0 && ilim <= Quick_max && try_quick && !rdir
325 #ifndef IMPRECISE_INEXACT
330 /* Try to get by with floating-point arithmetic. */
335 if ( (j = 11 - hi0bits(word0(&d) & Frac_mask)) !=0)
340 ieps = 2; /* conservative */
345 /* prevent overflows */
347 dval(&d) /= bigtens[n_bigtens-1];
350 for(; j; j >>= 1, i++)
358 if ( (j1 = -k) !=0) {
359 dval(&d) *= tens[j1 & 0xf];
360 for(j = j1 >> 4; j; j >>= 1, i++)
363 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) = ds*0.5/tens[ilim-1] - dval(&eps);
393 L = (Long)(dval(&d)/ds);
396 if (dval(&d) < dval(&eps)) {
398 inex = STRTOG_Inexlo;
401 if (ds - 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.) {
414 if ( (L = (Long)(dval(&d)/ds)) !=0)
419 if (dval(&d) > ds + dval(&eps))
421 else if (dval(&d) < ds - dval(&eps)) {
423 inex = STRTOG_Inexlo;
424 goto clear_trailing0;
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.) {
453 #ifdef Check_FLT_ROUNDS
454 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
467 inex = STRTOG_Inexlo;
470 dval(&d) += dval(&d);
474 if (dval(&d) > ds || (dval(&d) == ds && L & 1))
478 inex = STRTOG_Inexhi;
488 inex = STRTOG_Inexlo;
505 if (be - i++ < fpi->emin)
507 i = be - fpi->emin + 1;
518 if ((i = ilim) < 0) {
527 if (m2 > 0 && s2 > 0) {
528 i = m2 < s2 ? m2 : s2;
536 mhi = pow5mult(mhi, m5);
541 if ( (j = b5 - m5) !=0)
551 /* Check for special case that d is a normalized power of 2. */
555 if (bbits == 1 && be0 > fpi->emin + 1) {
556 /* The special case */
563 /* Arrange for convenient computation of quotients:
564 * shift left if necessary so divisor has 4 leading 0 bits.
566 * Perhaps we should just compute leading 28 bits of S once
567 * and for all and pass them and a shift to quorem, so it
568 * can do shifts and ors to compute the numerator for q.
570 i = ((s5 ? hi0bits(S->x[S->wds-1]) : ULbits - 1) - s2 - 4) & kmask;
579 b = multadd(b, 10, 0); /* we botched the k estimate */
581 mhi = multadd(mhi, 10, 0);
585 if (ilim <= 0 && mode > 2) {
586 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
587 /* no digits, fcvt style */
590 inex = STRTOG_Inexlo;
594 inex = STRTOG_Inexhi;
601 mhi = lshift(mhi, m2);
603 /* Compute mlo -- check for special case
604 * that d is a normalized power of 2.
609 mhi = Balloc(mhi->k);
611 mhi = lshift(mhi, 1);
615 dig = quorem(b,S) + '0';
616 /* Do we yet have the shortest decimal string
617 * that will round to d?
620 delta = diff(S, mhi);
621 j1 = delta->sign ? 1 : cmp(b, delta);
624 if (j1 == 0 && !mode && !(bits[0] & 1) && !rdir) {
628 if (b->wds > 1 || b->x[0])
629 inex = STRTOG_Inexlo;
633 inex = STRTOG_Inexhi;
639 if (j < 0 || (j == 0 && !mode
644 if (rdir && (b->wds > 1 || b->x[0])) {
646 inex = STRTOG_Inexlo;
649 while (cmp(S,mhi) > 0) {
651 mhi1 = multadd(mhi, 10, 0);
655 b = multadd(b, 10, 0);
656 dig = quorem(b,S) + '0';
660 inex = STRTOG_Inexhi;
669 if ((j1 > 0 || (j1 == 0 && dig & 1))
673 inex = STRTOG_Inexhi;
675 if (b->wds > 1 || b->x[0])
676 inex = STRTOG_Inexlo;
681 if (j1 > 0 && rdir != 2) {
682 if (dig == '9') { /* possible if i == 1 */
685 inex = STRTOG_Inexhi;
688 inex = STRTOG_Inexhi;
695 b = multadd(b, 10, 0);
697 mlo = mhi = multadd(mhi, 10, 0);
699 mlo = multadd(mlo, 10, 0);
700 mhi = multadd(mhi, 10, 0);
706 *s++ = dig = quorem(b,S) + '0';
709 b = multadd(b, 10, 0);
712 /* Round off last digit */
715 if (rdir == 2 || (b->wds <= 1 && !b->x[0]))
724 if (j > 0 || (j == 0 && dig & 1))
728 inex = STRTOG_Inexhi;
739 if (b->wds > 1 || b->x[0])
740 inex = STRTOG_Inexlo;
747 if (mlo && mlo != mhi)