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