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