![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to _dtoa.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10no / libc / stdio|
Return to _dtoa.c CVS log | Up to [Research Unix] / researchv10no / libc / stdio |
1.1 root 1: #include "fconv.h"
2:
3: static int quorem(Bigint *, Bigint *);
4:
5: /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
6: *
7: * Inspired by "How to Print Floating-Point Numbers Accurately" by
8: * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
9: *
10: * Modifications:
11: * 1. Rather than iterating, we use a simple numeric overestimate
12: * to determine k = floor(log10(d)). We scale relevant
13: * quantities using O(log2(k)) rather than O(k) multiplications.
14: * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
15: * try to generate digits strictly left to right. Instead, we
16: * compute with fewer bits and propagate the carry if necessary
17: * when rounding the final digit up. This is often faster.
18: * 3. Under the assumption that input will be rounded nearest,
19: * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
20: * That is, we allow equality in stopping tests when the
21: * round-nearest rule will give the same floating-point value
22: * as would satisfaction of the stopping test with strict
23: * inequality.
24: * 4. We remove common factors of powers of 2 from relevant
25: * quantities.
26: * 5. When converting floating-point integers less than 1e16,
27: * we use floating-point arithmetic rather than resorting
28: * to multiple-precision integers.
29: * 6. When asked to produce fewer than 15 digits, we first try
30: * to get by with floating-point arithmetic; we resort to
31: * multiple-precision integer arithmetic only if we cannot
32: * guarantee that the floating-point calculation has given
33: * the correctly rounded result. For k requested digits and
34: * "uniformly" distributed input, the probability is
35: * something like 10^(k-15) that we must resort to the long
36: * calculation.
37: */
38:
39: char *
40: _dtoa(double darg, int mode, int ndigits, int *decpt, int *sign, char **rve)
41: {
42: /* Arguments ndigits, decpt, sign are similar to those
43: of ecvt and fcvt; trailing zeros are suppressed from
44: the returned string. If not null, *rve is set to point
45: to the end of the return value. If d is +-Infinity or NaN,
46: then *decpt is set to 9999.
47:
48: mode:
49: 0 ==> shortest string that yields d when read in
50: and rounded to nearest.
51: 1 ==> like 0, but with Steele & White stopping rule;
52: e.g. with IEEE P754 arithmetic , mode 0 gives
53: 1e23 whereas mode 1 gives 9.999999999999999e22.
54: 2 ==> max(1,ndigits) significant digits. This gives a
55: return value similar to that of ecvt, except
56: that trailing zeros are suppressed.
57: 3 ==> through ndigits past the decimal point. This
58: gives a return value similar to that from fcvt,
59: except that trailing zeros are suppressed, and
60: ndigits can be negative.
61: 4-9 should give the same return values as 2-3, i.e.,
62: 4 <= mode <= 9 ==> same return as mode
63: 2 + (mode & 1). These modes are mainly for
64: debugging; often they run slower but sometimes
65: faster than modes 2-3.
66: 4,5,8,9 ==> left-to-right digit generation.
67: 6-9 ==> don't try fast floating-point estimate
68: (if applicable).
69:
70: Values of mode other than 0-9 are treated as mode 0.
71:
72: Sufficient space is allocated to the return value
73: to hold the suppressed trailing zeros.
74: */
75:
76: int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
77: j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
78: spec_case, try_quick;
79: long L;
80: #ifndef Sudden_Underflow
81: int denorm;
82: unsigned long x;
83: #endif
84: Bigint *b, *b1, *delta, *mlo, *mhi, *S;
85: double ds;
86: Dul d2, eps;
87: char *s, *s0;
88: static Bigint *result;
89: static int result_k;
90: Dul d;
91:
92: d.d = darg;
93: if (result) {
94: result->k = result_k;
95: result->maxwds = 1 << result_k;
96: Bfree(result);
97: result = 0;
98: }
99:
100: if (word0(d) & Sign_bit) {
101: /* set sign for everything, including 0's and NaNs */
102: *sign = 1;
103: word0(d) &= ~Sign_bit; /* clear sign bit */
104: }
105: else
106: *sign = 0;
107:
108: #if defined(IEEE_Arith) + defined(VAX)
109: #ifdef IEEE_Arith
110: if ((word0(d) & Exp_mask) == Exp_mask)
111: #else
112: if (word0(d) == 0x8000)
113: #endif
114: {
115: /* Infinity or NaN */
116: *decpt = 9999;
117: s =
118: #ifdef IEEE_Arith
119: !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" :
120: #endif
121: "NaN";
122: if (rve)
123: *rve =
124: #ifdef IEEE_Arith
125: s[3] ? s + 8 :
126: #endif
127: s + 3;
128: return s;
129: }
130: #endif
131: #ifdef IBM
132: d.d += 0; /* normalize */
133: #endif
134: if (!d.d) {
135: *decpt = 1;
136: s = "0";
137: if (rve)
138: *rve = s + 1;
139: return s;
140: }
141:
142: b = d2b(d.d, &be, &bbits);
143: #ifdef Sudden_Underflow
144: i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
145: #else
146: if (i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) {
147: #endif
148: d2.d = d.d;
149: word0(d2) &= Frac_mask1;
150: word0(d2) |= Exp_11;
151: #ifdef IBM
152: if (j = 11 - hi0bits(word0(d2) & Frac_mask))
153: d2.d /= 1 << j;
154: #endif
155:
156: /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
157: * log10(x) = log(x) / log(10)
158: * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
159: * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
160: *
161: * This suggests computing an approximation k to log10(d) by
162: *
163: * k = (i - Bias)*0.301029995663981
164: * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
165: *
166: * We want k to be too large rather than too small.
167: * The error in the first-order Taylor series approximation
168: * is in our favor, so we just round up the constant enough
169: * to compensate for any error in the multiplication of
170: * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
171: * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
172: * adding 1e-13 to the constant term more than suffices.
173: * Hence we adjust the constant term to 0.1760912590558.
174: * (We could get a more accurate k by invoking log10,
175: * but this is probably not worthwhile.)
176: */
177:
178: i -= Bias;
179: #ifdef IBM
180: i <<= 2;
181: i += j;
182: #endif
183: #ifndef Sudden_Underflow
184: denorm = 0;
185: }
186: else {
187: /* d is denormalized */
188:
189: i = bbits + be + (Bias + (P-1) - 1);
190: x = i > 32 ? word0(d) << 64 - i | word1(d) >> i - 32
191: : word1(d) << 32 - i;
192: d2.d = x;
193: word0(d2) -= 31*Exp_msk1; /* adjust exponent */
194: i -= (Bias + (P-1) - 1) + 1;
195: denorm = 1;
196: }
197: #endif
198: ds = (d2.d-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
199: k = (int)ds;
200: if (ds < 0. && ds != k)
201: k--; /* want k = floor(ds) */
202: k_check = 1;
203: if (k >= 0 && k <= Ten_pmax) {
204: if (d.d < tens[k])
205: k--;
206: k_check = 0;
207: }
208: j = bbits - i - 1;
209: if (j >= 0) {
210: b2 = 0;
211: s2 = j;
212: }
213: else {
214: b2 = -j;
215: s2 = 0;
216: }
217: if (k >= 0) {
218: b5 = 0;
219: s5 = k;
220: s2 += k;
221: }
222: else {
223: b2 -= k;
224: b5 = -k;
225: s5 = 0;
226: }
227: if (mode < 0 || mode > 9)
228: mode = 0;
229: try_quick = 1;
230: if (mode > 5) {
231: mode -= 4;
232: try_quick = 0;
233: }
234: leftright = 1;
235: switch(mode) {
236: case 0:
237: case 1:
238: ilim = ilim1 = -1;
239: i = 18;
240: ndigits = 0;
241: break;
242: case 2:
243: leftright = 0;
244: /* no break */
245: case 4:
246: if (ndigits <= 0)
247: ndigits = 1;
248: ilim = ilim1 = i = ndigits;
249: break;
250: case 3:
251: leftright = 0;
252: /* no break */
253: case 5:
254: i = ndigits + k + 1;
255: ilim = i;
256: ilim1 = i - 1;
257: if (i <= 0)
258: i = 1;
259: }
260: j = sizeof(unsigned long);
261: for(result_k = 0; sizeof(Bigint) - sizeof(unsigned long) + j < i;
262: j <<= 1) result_k++;
263: result = Balloc(result_k);
264: s = s0 = (char *)result;
265:
266: if (ilim >= 0 && ilim <= Quick_max && try_quick) {
267:
268: /* Try to get by with floating-point arithmetic. */
269:
270: i = 0;
271: d2.d = d.d;
272: k0 = k;
273: ilim0 = ilim;
274: ieps = 2; /* conservative */
275: if (k > 0) {
276: ds = tens[k&0xf];
277: j = k >> 4;
278: if (j & Bletch) {
279: /* prevent overflows */
280: j &= Bletch - 1;
281: d.d /= bigtens[n_bigtens-1];
282: ieps++;
283: }
284: for(; j; j >>= 1, i++)
285: if (j & 1) {
286: ieps++;
287: ds *= bigtens[i];
288: }
289: d.d /= ds;
290: }
291: else if (j1 = -k) {
292: d.d *= tens[j1 & 0xf];
293: for(j = j1 >> 4; j; j >>= 1, i++)
294: if (j & 1) {
295: ieps++;
296: d.d *= bigtens[i];
297: }
298: }
299: if (k_check && d.d < 1. && ilim > 0) {
300: if (ilim1 <= 0)
301: goto fast_failed;
302: ilim = ilim1;
303: k--;
304: d.d *= 10.;
305: ieps++;
306: }
307: eps.d = ieps*d.d + 7.;
308: word0(eps) -= (P-1)*Exp_msk1;
309: if (ilim == 0) {
310: S = mhi = 0;
311: d.d -= 5.;
312: if (d.d > eps.d)
313: goto one_digit;
314: if (d.d < -eps.d)
315: goto no_digits;
316: goto fast_failed;
317: }
318: #ifndef No_leftright
319: if (leftright) {
320: /* Use Steele & White method of only
321: * generating digits needed.
322: */
323: eps.d = 0.5/tens[ilim-1] - eps.d;
324: for(i = 0;;) {
325: L = d.d;
326: d.d -= L;
327: *s++ = '0' + (int)L;
328: if (d.d < eps.d)
329: goto ret1;
330: if (1. - d.d < eps.d)
331: goto bump_up;
332: if (++i >= ilim)
333: break;
334: eps.d *= 10.;
335: d.d *= 10.;
336: }
337: }
338: else {
339: #endif
340: /* Generate ilim digits, then fix them up. */
341: eps.d *= tens[ilim-1];
342: for(i = 1;; i++, d.d *= 10.) {
343: L = d.d;
344: d.d -= L;
345: *s++ = '0' + (int)L;
346: if (i == ilim) {
347: if (d.d > 0.5 + eps.d)
348: goto bump_up;
349: else if (d.d < 0.5 - eps.d) {
350: while(*--s == '0');
351: s++;
352: goto ret1;
353: }
354: break;
355: }
356: }
357: #ifndef No_leftright
358: }
359: #endif
360: fast_failed:
361: s = s0;
362: d.d = d2.d;
363: k = k0;
364: ilim = ilim0;
365: }
366:
367: /* Do we have a "small" integer? */
368:
369: if (be >= 0 && k <= Int_max) {
370: /* Yes. */
371: ds = tens[k];
372: if (ndigits < 0 && ilim <= 0) {
373: S = mhi = 0;
374: if (ilim < 0 || d.d <= 5*ds)
375: goto no_digits;
376: goto one_digit;
377: }
378: for(i = 1;; i++) {
379: L = d.d / ds;
380: d.d -= L*ds;
381: #ifdef Check_FLT_ROUNDS
382: /* If FLT_ROUNDS == 2, L will usually be high by 1 */
383: if (d.d < 0) {
384: L--;
385: d.d += ds;
386: }
387: #endif
388: *s++ = '0' + (int)L;
389: if (i == ilim) {
390: d.d += d.d;
391: if (d.d > ds || d.d == ds && L & 1) {
392: bump_up:
393: while(*--s == '9')
394: if (s == s0) {
395: k++;
396: *s = '0';
397: break;
398: }
399: ++*s++;
400: }
401: break;
402: }
403: if (!(d.d *= 10.))
404: break;
405: }
406: goto ret1;
407: }
408:
409: m2 = b2;
410: m5 = b5;
411: mhi = mlo = 0;
412: if (leftright) {
413: if (mode < 2) {
414: i =
415: #ifndef Sudden_Underflow
416: denorm ? be + (Bias + (P-1) - 1 + 1) :
417: #endif
418: #ifdef IBM
419: 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
420: #else
421: 1 + P - bbits;
422: #endif
423: }
424: else {
425: j = ilim - 1;
426: if (m5 >= j)
427: m5 -= j;
428: else {
429: s5 += j -= m5;
430: b5 += j;
431: m5 = 0;
432: }
433: if ((i = ilim) < 0) {
434: m2 -= i;
435: i = 0;
436: }
437: }
438: b2 += i;
439: s2 += i;
440: mhi = i2b(1);
441: }
442: if (m2 > 0 && s2 > 0) {
443: i = m2 < s2 ? m2 : s2;
444: b2 -= i;
445: m2 -= i;
446: s2 -= i;
447: }
448: if (b5 > 0) {
449: if (leftright) {
450: if (m5 > 0) {
451: mhi = pow5mult(mhi, m5);
452: b1 = mult(mhi, b);
453: Bfree(b);
454: b = b1;
455: }
456: if (j = b5 - m5)
457: b = pow5mult(b, j);
458: }
459: else
460: b = pow5mult(b, b5);
461: }
462: S = i2b(1);
463: if (s5 > 0)
464: S = pow5mult(S, s5);
465:
466: /* Check for special case that d is a normalized power of 2. */
467:
468: if (mode < 2) {
469: if (!word1(d) && !(word0(d) & Bndry_mask)
470: #ifndef Sudden_Underflow
471: && word0(d) & Exp_mask
472: #endif
473: ) {
474: /* The special case */
475: b2 += Log2P;
476: s2 += Log2P;
477: spec_case = 1;
478: }
479: else
480: spec_case = 0;
481: }
482:
483: /* Arrange for convenient computation of quotients:
484: * shift left if necessary so divisor has 4 leading 0 bits.
485: *
486: * Perhaps we should just compute leading 28 bits of S once
487: * and for all and pass them and a shift to quorem, so it
488: * can do shifts and ors to compute the numerator for q.
489: */
490: #ifdef Pack_32
491: if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
492: i = 32 - i;
493: #else
494: if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
495: i = 16 - i;
496: #endif
497: if (i > 4) {
498: i -= 4;
499: b2 += i;
500: m2 += i;
501: s2 += i;
502: }
503: else if (i < 4) {
504: i += 28;
505: b2 += i;
506: m2 += i;
507: s2 += i;
508: }
509: if (b2 > 0)
510: b = lshift(b, b2);
511: if (s2 > 0)
512: S = lshift(S, s2);
513: if (k_check) {
514: if (cmp(b,S) < 0) {
515: k--;
516: b = multadd(b, 10, 0); /* we botched the k estimate */
517: if (leftright)
518: mhi = multadd(mhi, 10, 0);
519: ilim = ilim1;
520: }
521: }
522: if (ilim <= 0 && mode > 2) {
523: if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
524: /* no digits, fcvt style */
525: no_digits:
526: k = -1 - ndigits;
527: goto ret;
528: }
529: one_digit:
530: *s++ = '1';
531: k++;
532: goto ret;
533: }
534: if (leftright) {
535: if (m2 > 0)
536: mhi = lshift(mhi, m2);
537:
538: /* Compute mlo -- check for special case
539: * that d is a normalized power of 2.
540: */
541:
542: mlo = mhi;
543: if (spec_case) {
544: mhi = Balloc(mhi->k);
545: Bcopy(mhi, mlo);
546: mhi = lshift(mhi, Log2P);
547: }
548:
549: for(i = 1;;i++) {
550: dig = quorem(b,S) + '0';
551: /* Do we yet have the shortest decimal string
552: * that will round to d?
553: */
554: j = cmp(b, mlo);
555: delta = diff(S, mhi);
556: j1 = delta->sign ? 1 : cmp(b, delta);
557: Bfree(delta);
558: #ifndef ROUND_BIASED
559: if (j1 == 0 && !mode && !(word1(d) & 1)) {
560: if (dig == '9')
561: goto round_9_up;
562: if (j > 0)
563: dig++;
564: *s++ = dig;
565: goto ret;
566: }
567: #endif
568: if (j < 0 || j == 0 && !mode
569: #ifndef ROUND_BIASED
570: && !(word1(d) & 1)
571: #endif
572: ) {
573: if (j1 > 0) {
574: b = lshift(b, 1);
575: j1 = cmp(b, S);
576: if ((j1 > 0 || j1 == 0 && dig & 1)
577: && dig++ == '9')
578: goto round_9_up;
579: }
580: *s++ = dig;
581: goto ret;
582: }
583: if (j1 > 0) {
584: if (dig == '9') { /* possible if i == 1 */
585: round_9_up:
586: *s++ = '9';
587: goto roundoff;
588: }
589: *s++ = dig + 1;
590: goto ret;
591: }
592: *s++ = dig;
593: if (i == ilim)
594: break;
595: b = multadd(b, 10, 0);
596: if (mlo == mhi)
597: mlo = mhi = multadd(mhi, 10, 0);
598: else {
599: mlo = multadd(mlo, 10, 0);
600: mhi = multadd(mhi, 10, 0);
601: }
602: }
603: }
604: else
605: for(i = 1;; i++) {
606: *s++ = dig = quorem(b,S) + '0';
607: if (i >= ilim)
608: break;
609: b = multadd(b, 10, 0);
610: }
611:
612: /* Round off last digit */
613:
614: b = lshift(b, 1);
615: j = cmp(b, S);
616: if (j > 0 || j == 0 && dig & 1) {
617: roundoff:
618: while(*--s == '9')
619: if (s == s0) {
620: k++;
621: *s++ = '1';
622: goto ret;
623: }
624: ++*s++;
625: }
626: else {
627: while(*--s == '0');
628: s++;
629: }
630: ret:
631: Bfree(S);
632: if (mhi) {
633: if (mlo && mlo != mhi)
634: Bfree(mlo);
635: Bfree(mhi);
636: }
637: ret1:
638: Bfree(b);
639: *s = 0;
640: *decpt = k + 1;
641: if (rve)
642: *rve = s;
643: return s0;
644: }
645:
646: static int
647: quorem(Bigint *b, Bigint *S)
648: {
649: int n;
650: long borrow, y;
651: unsigned long carry, q, ys;
652: unsigned long *bx, *bxe, *sx, *sxe;
653: #ifdef Pack_32
654: long z;
655: unsigned long si, zs;
656: #endif
657:
658: n = S->wds;
659: #ifdef DEBUG
660: /*debug*/ if (b->wds > n)
661: /*debug*/ Bug("oversize b in quorem");
662: #endif
663: if (b->wds < n)
664: return 0;
665: sx = S->x;
666: sxe = sx + --n;
667: bx = b->x;
668: bxe = bx + n;
669: q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
670: #ifdef DEBUG
671: /*debug*/ if (q > 9)
672: /*debug*/ Bug("oversized quotient in quorem");
673: #endif
674: if (q) {
675: borrow = 0;
676: carry = 0;
677: do {
678: #ifdef Pack_32
679: si = *sx++;
680: ys = (si & 0xffff) * q + carry;
681: zs = (si >> 16) * q + (ys >> 16);
682: carry = zs >> 16;
683: y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
684: borrow = y >> 16;
685: Sign_Extend(borrow, y);
686: z = (*bx >> 16) - (zs & 0xffff) + borrow;
687: borrow = z >> 16;
688: Sign_Extend(borrow, z);
689: Storeinc(bx, z, y);
690: #else
691: ys = *sx++ * q + carry;
692: carry = ys >> 16;
693: y = *bx - (ys & 0xffff) + borrow;
694: borrow = y >> 16;
695: Sign_Extend(borrow, y);
696: *bx++ = y & 0xffff;
697: #endif
698: }
699: while(sx <= sxe);
700: if (!*bxe) {
701: bx = b->x;
702: while(--bxe > bx && !*bxe)
703: --n;
704: b->wds = n;
705: }
706: }
707: if (cmp(b, S) >= 0) {
708: q++;
709: borrow = 0;
710: carry = 0;
711: bx = b->x;
712: sx = S->x;
713: do {
714: #ifdef Pack_32
715: si = *sx++;
716: ys = (si & 0xffff) + carry;
717: zs = (si >> 16) + (ys >> 16);
718: carry = zs >> 16;
719: y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
720: borrow = y >> 16;
721: Sign_Extend(borrow, y);
722: z = (*bx >> 16) - (zs & 0xffff) + borrow;
723: borrow = z >> 16;
724: Sign_Extend(borrow, z);
725: Storeinc(bx, z, y);
726: #else
727: ys = *sx++ + carry;
728: carry = ys >> 16;
729: y = *bx - (ys & 0xffff) + borrow;
730: borrow = y >> 16;
731: Sign_Extend(borrow, y);
732: *bx++ = y & 0xffff;
733: #endif
734: }
735: while(sx <= sxe);
736: bx = b->x;
737: bxe = bx + n;
738: if (!*bxe) {
739: while(--bxe > bx && !*bxe)
740: --n;
741: b->wds = n;
742: }
743: }
744: return q;
745: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.