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