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