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