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