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