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