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