1 /****************************************************************
2
3 The author of this software is David M. Gay.
4
5 Copyright (C) 1998, 1999 by Lucent Technologies
6 All Rights Reserved
7
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
16 permission.
17
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
25 THIS 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
bitstob(ULong * bits,int nbits,int * bbits)34 static 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
__gdtoa(FPI * fpi,int be,ULong * bits,int * kindp,int mode,int ndigits,int * decpt,char ** rve)106 char *__gdtoa (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 /* no break */
274 case 4:
275 if (ndigits <= 0)
276 ndigits = 1;
277 ilim = ilim1 = i = ndigits;
278 break;
279 case 3:
280 leftright = 0;
281 /* no break */
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 ( (rdir = fpi->rounding - 1) !=0) {
292 if (rdir < 0)
293 rdir = 2;
294 if (kind & STRTOG_Neg)
295 rdir = 3 - rdir;
296 }
297
298 /* Now rdir = 0 ==> round near, 1 ==> round up, 2 ==> round down. */
299
300 if (ilim >= 0 && ilim <= Quick_max && try_quick && !rdir
301 #ifndef IMPRECISE_INEXACT
302 && k == 0
303 #endif
304 ) {
305
306 /* Try to get by with floating-point arithmetic. */
307
308 i = 0;
309 d2 = dval(&d);
310 k0 = k;
311 ilim0 = ilim;
312 ieps = 2; /* conservative */
313 if (k > 0) {
314 ds = tens[k&0xf];
315 j = k >> 4;
316 if (j & Bletch) {
317 /* prevent overflows */
318 j &= Bletch - 1;
319 dval(&d) /= bigtens[n_bigtens-1];
320 ieps++;
321 }
322 for(; j; j >>= 1, i++)
323 if (j & 1) {
324 ieps++;
325 ds *= bigtens[i];
326 }
327 }
328 else {
329 ds = 1.;
330 if ( (j2 = -k) !=0) {
331 dval(&d) *= tens[j2 & 0xf];
332 for(j = j2 >> 4; j; j >>= 1, i++)
333 if (j & 1) {
334 ieps++;
335 dval(&d) *= bigtens[i];
336 }
337 }
338 }
339 if (k_check && dval(&d) < 1. && ilim > 0) {
340 if (ilim1 <= 0)
341 goto fast_failed;
342 ilim = ilim1;
343 k--;
344 dval(&d) *= 10.;
345 ieps++;
346 }
347 dval(&eps) = ieps*dval(&d) + 7.;
348 word0(&eps) -= (P-1)*Exp_msk1;
349 if (ilim == 0) {
350 S = mhi = 0;
351 dval(&d) -= 5.;
352 if (dval(&d) > dval(&eps))
353 goto one_digit;
354 if (dval(&d) < -dval(&eps))
355 goto no_digits;
356 goto fast_failed;
357 }
358 #ifndef No_leftright
359 if (leftright) {
360 /* Use Steele & White method of only
361 * generating digits needed.
362 */
363 dval(&eps) = ds*0.5/tens[ilim-1] - dval(&eps);
364 for(i = 0;;) {
365 L = (Long)(dval(&d)/ds);
366 dval(&d) -= L*ds;
367 *s++ = '0' + (int)L;
368 if (dval(&d) < dval(&eps)) {
369 if (dval(&d))
370 inex = STRTOG_Inexlo;
371 goto ret1;
372 }
373 if (ds - dval(&d) < dval(&eps))
374 goto bump_up;
375 if (++i >= ilim)
376 break;
377 dval(&eps) *= 10.;
378 dval(&d) *= 10.;
379 }
380 }
381 else {
382 #endif
383 /* Generate ilim digits, then fix them up. */
384 dval(&eps) *= tens[ilim-1];
385 for(i = 1;; i++, dval(&d) *= 10.) {
386 if ( (L = (Long)(dval(&d)/ds)) !=0)
387 dval(&d) -= L*ds;
388 *s++ = '0' + (int)L;
389 if (i == ilim) {
390 ds *= 0.5;
391 if (dval(&d) > ds + dval(&eps))
392 goto bump_up;
393 else if (dval(&d) < ds - dval(&eps)) {
394 if (dval(&d))
395 inex = STRTOG_Inexlo;
396 goto clear_trailing0;
397 }
398 break;
399 }
400 }
401 #ifndef No_leftright
402 }
403 #endif
404 fast_failed:
405 s = s0;
406 dval(&d) = d2;
407 k = k0;
408 ilim = ilim0;
409 }
410
411 /* Do we have a "small" integer? */
412
413 if (be >= 0 && k <= fpi->int_max) {
414 /* Yes. */
415 ds = tens[k];
416 if (ndigits < 0 && ilim <= 0) {
417 S = mhi = 0;
418 if (ilim < 0 || dval(&d) <= 5*ds)
419 goto no_digits;
420 goto one_digit;
421 }
422 for(i = 1;; i++, dval(&d) *= 10.) {
423 L = dval(&d) / ds;
424 dval(&d) -= L*ds;
425 #ifdef Check_FLT_ROUNDS
426 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
427 if (dval(&d) < 0) {
428 L--;
429 dval(&d) += ds;
430 }
431 #endif
432 *s++ = '0' + (int)L;
433 if (dval(&d) == 0.)
434 break;
435 if (i == ilim) {
436 if (rdir) {
437 if (rdir == 1)
438 goto bump_up;
439 inex = STRTOG_Inexlo;
440 goto ret1;
441 }
442 dval(&d) += dval(&d);
443 #ifdef ROUND_BIASED
444 if (dval(&d) >= ds)
445 #else
446 if (dval(&d) > ds || (dval(&d) == ds && L & 1))
447 #endif
448 {
449 bump_up:
450 inex = STRTOG_Inexhi;
451 while(*--s == '9')
452 if (s == s0) {
453 k++;
454 *s = '0';
455 break;
456 }
457 ++*s++;
458 }
459 else {
460 inex = STRTOG_Inexlo;
461 clear_trailing0:
462 while(*--s == '0'){}
463 ++s;
464 }
465 break;
466 }
467 }
468 goto ret1;
469 }
470
471 m2 = b2;
472 m5 = b5;
473 mhi = mlo = 0;
474 if (leftright) {
475 i = nbits - bbits;
476 if (be - i++ < fpi->emin && mode != 3 && mode != 5) {
477 /* denormal */
478 i = be - fpi->emin + 1;
479 if (mode >= 2 && ilim > 0 && ilim < i)
480 goto small_ilim;
481 }
482 else if (mode >= 2) {
483 small_ilim:
484 j = ilim - 1;
485 if (m5 >= j)
486 m5 -= j;
487 else {
488 s5 += j -= m5;
489 b5 += j;
490 m5 = 0;
491 }
492 if ((i = ilim) < 0) {
493 m2 -= i;
494 i = 0;
495 }
496 }
497 b2 += i;
498 s2 += i;
499 mhi = i2b(1);
500 }
501 if (m2 > 0 && s2 > 0) {
502 i = m2 < s2 ? m2 : s2;
503 b2 -= i;
504 m2 -= i;
505 s2 -= i;
506 }
507 if (b5 > 0) {
508 if (leftright) {
509 if (m5 > 0) {
510 mhi = pow5mult(mhi, m5);
511 b1 = mult(mhi, b);
512 Bfree(b);
513 b = b1;
514 }
515 if ( (j = b5 - m5) !=0)
516 b = pow5mult(b, j);
517 }
518 else
519 b = pow5mult(b, b5);
520 }
521 S = i2b(1);
522 if (s5 > 0)
523 S = pow5mult(S, s5);
524
525 /* Check for special case that d is a normalized power of 2. */
526
527 spec_case = 0;
528 if (mode < 2) {
529 if (bbits == 1 && be0 > fpi->emin + 1) {
530 /* The special case */
531 b2++;
532 s2++;
533 spec_case = 1;
534 }
535 }
536
537 /* Arrange for convenient computation of quotients:
538 * shift left if necessary so divisor has 4 leading 0 bits.
539 *
540 * Perhaps we should just compute leading 28 bits of S once
541 * and for all and pass them and a shift to quorem, so it
542 * can do shifts and ors to compute the numerator for q.
543 */
544 i = ((s5 ? hi0bits(S->x[S->wds-1]) : ULbits - 1) - s2 - 4) & kmask;
545 m2 += i;
546 if ((b2 += i) > 0)
547 b = lshift(b, b2);
548 if ((s2 += i) > 0)
549 S = lshift(S, s2);
550 if (k_check) {
551 if (cmp(b,S) < 0) {
552 k--;
553 b = multadd(b, 10, 0); /* we botched the k estimate */
554 if (leftright)
555 mhi = multadd(mhi, 10, 0);
556 ilim = ilim1;
557 }
558 }
559 if (ilim <= 0 && mode > 2) {
560 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
561 /* no digits, fcvt style */
562 no_digits:
563 k = -1 - ndigits;
564 inex = STRTOG_Inexlo;
565 goto ret;
566 }
567 one_digit:
568 inex = STRTOG_Inexhi;
569 *s++ = '1';
570 k++;
571 goto ret;
572 }
573 if (leftright) {
574 if (m2 > 0)
575 mhi = lshift(mhi, m2);
576
577 /* Compute mlo -- check for special case
578 * that d is a normalized power of 2.
579 */
580
581 mlo = mhi;
582 if (spec_case) {
583 mhi = Balloc(mhi->k);
584 Bcopy(mhi, mlo);
585 mhi = lshift(mhi, 1);
586 }
587
588 for(i = 1;;i++) {
589 dig = quorem(b,S) + '0';
590 /* Do we yet have the shortest decimal string
591 * that will round to d?
592 */
593 j = cmp(b, mlo);
594 delta = diff(S, mhi);
595 j2 = delta->sign ? 1 : cmp(b, delta);
596 Bfree(delta);
597 #ifndef ROUND_BIASED
598 if (j2 == 0 && !mode && !(bits[0] & 1) && !rdir) {
599 if (dig == '9')
600 goto round_9_up;
601 if (j <= 0) {
602 if (b->wds > 1 || b->x[0])
603 inex = STRTOG_Inexlo;
604 }
605 else {
606 dig++;
607 inex = STRTOG_Inexhi;
608 }
609 *s++ = dig;
610 goto ret;
611 }
612 #endif
613 if (j < 0 || (j == 0 && !mode
614 #ifndef ROUND_BIASED
615 && !(bits[0] & 1)
616 #endif
617 )) {
618 if (rdir && (b->wds > 1 || b->x[0])) {
619 if (rdir == 2) {
620 inex = STRTOG_Inexlo;
621 goto accept;
622 }
623 while (cmp(S,mhi) > 0) {
624 *s++ = dig;
625 mhi1 = multadd(mhi, 10, 0);
626 if (mlo == mhi)
627 mlo = mhi1;
628 mhi = mhi1;
629 b = multadd(b, 10, 0);
630 dig = quorem(b,S) + '0';
631 }
632 if (dig++ == '9')
633 goto round_9_up;
634 inex = STRTOG_Inexhi;
635 goto accept;
636 }
637 if (j2 > 0) {
638 b = lshift(b, 1);
639 j2 = cmp(b, S);
640 #ifdef ROUND_BIASED
641 if (j2 >= 0 /*)*/
642 #else
643 if ((j2 > 0 || (j2 == 0 && dig & 1))
644 #endif
645 && dig++ == '9')
646 goto round_9_up;
647 inex = STRTOG_Inexhi;
648 }
649 if (b->wds > 1 || b->x[0])
650 inex = STRTOG_Inexlo;
651 accept:
652 *s++ = dig;
653 goto ret;
654 }
655 if (j2 > 0 && rdir != 2) {
656 if (dig == '9') { /* possible if i == 1 */
657 round_9_up:
658 *s++ = '9';
659 inex = STRTOG_Inexhi;
660 goto roundoff;
661 }
662 inex = STRTOG_Inexhi;
663 *s++ = dig + 1;
664 goto ret;
665 }
666 *s++ = dig;
667 if (i == ilim)
668 break;
669 b = multadd(b, 10, 0);
670 if (mlo == mhi)
671 mlo = mhi = multadd(mhi, 10, 0);
672 else {
673 mlo = multadd(mlo, 10, 0);
674 mhi = multadd(mhi, 10, 0);
675 }
676 }
677 }
678 else
679 for(i = 1;; i++) {
680 *s++ = dig = quorem(b,S) + '0';
681 if (i >= ilim)
682 break;
683 b = multadd(b, 10, 0);
684 }
685
686 /* Round off last digit */
687
688 if (rdir) {
689 if (rdir == 2 || (b->wds <= 1 && !b->x[0]))
690 goto chopzeros;
691 goto roundoff;
692 }
693 b = lshift(b, 1);
694 j = cmp(b, S);
695 #ifdef ROUND_BIASED
696 if (j >= 0)
697 #else
698 if (j > 0 || (j == 0 && dig & 1))
699 #endif
700 {
701 roundoff:
702 inex = STRTOG_Inexhi;
703 while(*--s == '9')
704 if (s == s0) {
705 k++;
706 *s++ = '1';
707 goto ret;
708 }
709 ++*s++;
710 }
711 else {
712 chopzeros:
713 if (b->wds > 1 || b->x[0])
714 inex = STRTOG_Inexlo;
715 while(*--s == '0'){}
716 ++s;
717 }
718 ret:
719 Bfree(S);
720 if (mhi) {
721 if (mlo && mlo != mhi)
722 Bfree(mlo);
723 Bfree(mhi);
724 }
725 ret1:
726 Bfree(b);
727 *s = 0;
728 *decpt = k + 1;
729 if (rve)
730 *rve = s;
731 *kindp |= inex;
732 return s0;
733 }
734