xref: /netbsd-src/lib/libc/gdtoa/gdtoa.c (revision 8b0f9554ff8762542c4defc4f70e1eb76fb508fa)
1 /* $NetBSD: gdtoa.c,v 1.3 2007/02/03 18:09:20 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
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 	be = bits + (((unsigned int)nbits - 1) >> kshift);
59 	x = x0 = b->x;
60 	do {
61 		*x++ = *bits & ALL_ON;
62 #ifdef Pack_16
63 		*x++ = (*bits >> 16) & ALL_ON;
64 #endif
65 		} while(++bits <= be);
66 	i = x - x0;
67 	while(!x0[--i])
68 		if (!i) {
69 			b->wds = 0;
70 			*bbits = 0;
71 			goto ret;
72 			}
73 	b->wds = i + 1;
74 	*bbits = i*ULbits + 32 - hi0bits(b->x[i]);
75  ret:
76 	return b;
77 	}
78 
79 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
80  *
81  * Inspired by "How to Print Floating-Point Numbers Accurately" by
82  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
83  *
84  * Modifications:
85  *	1. Rather than iterating, we use a simple numeric overestimate
86  *	   to determine k = floor(log10(d)).  We scale relevant
87  *	   quantities using O(log2(k)) rather than O(k) multiplications.
88  *	2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
89  *	   try to generate digits strictly left to right.  Instead, we
90  *	   compute with fewer bits and propagate the carry if necessary
91  *	   when rounding the final digit up.  This is often faster.
92  *	3. Under the assumption that input will be rounded nearest,
93  *	   mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
94  *	   That is, we allow equality in stopping tests when the
95  *	   round-nearest rule will give the same floating-point value
96  *	   as would satisfaction of the stopping test with strict
97  *	   inequality.
98  *	4. We remove common factors of powers of 2 from relevant
99  *	   quantities.
100  *	5. When converting floating-point integers less than 1e16,
101  *	   we use floating-point arithmetic rather than resorting
102  *	   to multiple-precision integers.
103  *	6. When asked to produce fewer than 15 digits, we first try
104  *	   to get by with floating-point arithmetic; we resort to
105  *	   multiple-precision integer arithmetic only if we cannot
106  *	   guarantee that the floating-point calculation has given
107  *	   the correctly rounded result.  For k requested digits and
108  *	   "uniformly" distributed input, the probability is
109  *	   something like 10^(k-15) that we must resort to the Long
110  *	   calculation.
111  */
112 
113  char *
114 gdtoa
115 #ifdef KR_headers
116 	(fpi, be, bits, kindp, mode, ndigits, decpt, rve)
117 	FPI *fpi; int be; ULong *bits;
118 	int *kindp, mode, ndigits, *decpt; char **rve;
119 #else
120 	(FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, int *decpt, char **rve)
121 #endif
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 	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 ( (jj1 = j & 3) !=0)
251 		dval(d) *= 1 << jj1;
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 	leftright = 1;
288 	switch(mode) {
289 		case 0:
290 		case 1:
291 			ilim = ilim1 = -1;
292 			i = (int)(nbits * .30103) + 3;
293 			ndigits = 0;
294 			break;
295 		case 2:
296 			leftright = 0;
297 			/*FALLTHROUGH*/
298 		case 4:
299 			if (ndigits <= 0)
300 				ndigits = 1;
301 			ilim = ilim1 = i = ndigits;
302 			break;
303 		case 3:
304 			leftright = 0;
305 			/*FALLTHROUGH*/
306 		case 5:
307 			i = ndigits + k + 1;
308 			ilim = i;
309 			ilim1 = i - 1;
310 			if (i <= 0)
311 				i = 1;
312 		}
313 	s = s0 = rv_alloc((size_t)i);
314 
315 	if ( (rdir = fpi->rounding - 1) !=0) {
316 		if (rdir < 0)
317 			rdir = 2;
318 		if (kind & STRTOG_Neg)
319 			rdir = 3 - rdir;
320 		}
321 
322 	/* Now rdir = 0 ==> round near, 1 ==> round up, 2 ==> round down. */
323 
324 	if (ilim >= 0 && ilim <= Quick_max && try_quick && !rdir
325 #ifndef IMPRECISE_INEXACT
326 		&& k == 0
327 #endif
328 								) {
329 
330 		/* Try to get by with floating-point arithmetic. */
331 
332 		i = 0;
333 		d2 = dval(d);
334 #ifdef IBM
335 		if ( (j = 11 - hi0bits(word0(d) & Frac_mask)) !=0)
336 			dval(d) /= 1 << j;
337 #endif
338 		k0 = k;
339 		ilim0 = ilim;
340 		ieps = 2; /* conservative */
341 		if (k > 0) {
342 			ds = tens[k&0xf];
343 			j = (unsigned int)k >> 4;
344 			if (j & Bletch) {
345 				/* prevent overflows */
346 				j &= Bletch - 1;
347 				dval(d) /= bigtens[n_bigtens-1];
348 				ieps++;
349 				}
350 			for(; j; j /= 2, i++)
351 				if (j & 1) {
352 					ieps++;
353 					ds *= bigtens[i];
354 					}
355 			}
356 		else  {
357 			ds = 1.;
358 			if ( (jj1 = -k) !=0) {
359 				dval(d) *= tens[jj1 & 0xf];
360 				for(j = jj1 >> 4; j; j >>= 1, i++)
361 					if (j & 1) {
362 						ieps++;
363 						dval(d) *= bigtens[i];
364 						}
365 				}
366 			}
367 		if (k_check && dval(d) < 1. && ilim > 0) {
368 			if (ilim1 <= 0)
369 				goto fast_failed;
370 			ilim = ilim1;
371 			k--;
372 			dval(d) *= 10.;
373 			ieps++;
374 			}
375 		dval(eps) = ieps*dval(d) + 7.;
376 		word0(eps) -= (P-1)*Exp_msk1;
377 		if (ilim == 0) {
378 			S = mhi = 0;
379 			dval(d) -= 5.;
380 			if (dval(d) > dval(eps))
381 				goto one_digit;
382 			if (dval(d) < -dval(eps))
383 				goto no_digits;
384 			goto fast_failed;
385 			}
386 #ifndef No_leftright
387 		if (leftright) {
388 			/* Use Steele & White method of only
389 			 * generating digits needed.
390 			 */
391 			dval(eps) = ds*0.5/tens[ilim-1] - dval(eps);
392 			for(i = 0;;) {
393 				L = (Long)(dval(d)/ds);
394 				dval(d) -= L*ds;
395 				*s++ = '0' + (int)L;
396 				if (dval(d) < dval(eps)) {
397 					if (dval(d))
398 						inex = STRTOG_Inexlo;
399 					goto ret1;
400 					}
401 				if (ds - dval(d) < dval(eps))
402 					goto bump_up;
403 				if (++i >= ilim)
404 					break;
405 				dval(eps) *= 10.;
406 				dval(d) *= 10.;
407 				}
408 			}
409 		else {
410 #endif
411 			/* Generate ilim digits, then fix them up. */
412 			dval(eps) *= tens[ilim-1];
413 			for(i = 1;; i++, dval(d) *= 10.) {
414 				if ( (L = (Long)(dval(d)/ds)) !=0)
415 					dval(d) -= L*ds;
416 				*s++ = '0' + (int)L;
417 				if (i == ilim) {
418 					ds *= 0.5;
419 					if (dval(d) > ds + dval(eps))
420 						goto bump_up;
421 					else if (dval(d) < ds - dval(eps)) {
422 						while(*--s == '0'){}
423 						s++;
424 						if (dval(d))
425 							inex = STRTOG_Inexlo;
426 						goto ret1;
427 						}
428 					break;
429 					}
430 				}
431 #ifndef No_leftright
432 			}
433 #endif
434  fast_failed:
435 		s = s0;
436 		dval(d) = d2;
437 		k = k0;
438 		ilim = ilim0;
439 		}
440 
441 	/* Do we have a "small" integer? */
442 
443 	if (be >= 0 && k <= Int_max) {
444 		/* Yes. */
445 		ds = tens[k];
446 		if (ndigits < 0 && ilim <= 0) {
447 			S = mhi = 0;
448 			if (ilim < 0 || dval(d) <= 5*ds)
449 				goto no_digits;
450 			goto one_digit;
451 			}
452 		for(i = 1;; i++, dval(d) *= 10.) {
453 			L = dval(d) / ds;
454 			dval(d) -= L*ds;
455 #ifdef Check_FLT_ROUNDS
456 			/* If FLT_ROUNDS == 2, L will usually be high by 1 */
457 			if (dval(d) < 0) {
458 				L--;
459 				dval(d) += ds;
460 				}
461 #endif
462 			*s++ = '0' + (int)L;
463 			if (dval(d) == 0.)
464 				break;
465 			if (i == ilim) {
466 				if (rdir) {
467 					if (rdir == 1)
468 						goto bump_up;
469 					inex = STRTOG_Inexlo;
470 					goto ret1;
471 					}
472 				dval(d) += dval(d);
473 				if (dval(d) > ds || (dval(d) == ds && L & 1)) {
474  bump_up:
475 					inex = STRTOG_Inexhi;
476 					while(*--s == '9')
477 						if (s == s0) {
478 							k++;
479 							*s = '0';
480 							break;
481 							}
482 					++*s++;
483 					}
484 				else
485 					inex = STRTOG_Inexlo;
486 				break;
487 				}
488 			}
489 		goto ret1;
490 		}
491 
492 	m2 = b2;
493 	m5 = b5;
494 	mhi = mlo = 0;
495 	if (leftright) {
496 		if (mode < 2) {
497 			i = nbits - bbits;
498 			if (be - i++ < fpi->emin)
499 				/* denormal */
500 				i = be - fpi->emin + 1;
501 			}
502 		else {
503 			j = ilim - 1;
504 			if (m5 >= j)
505 				m5 -= j;
506 			else {
507 				s5 += j -= m5;
508 				b5 += j;
509 				m5 = 0;
510 				}
511 			if ((i = ilim) < 0) {
512 				m2 -= i;
513 				i = 0;
514 				}
515 			}
516 		b2 += i;
517 		s2 += i;
518 		mhi = i2b(1);
519 		}
520 	if (m2 > 0 && s2 > 0) {
521 		i = m2 < s2 ? m2 : s2;
522 		b2 -= i;
523 		m2 -= i;
524 		s2 -= i;
525 		}
526 	if (b5 > 0) {
527 		if (leftright) {
528 			if (m5 > 0) {
529 				mhi = pow5mult(mhi, m5);
530 				b1 = mult(mhi, b);
531 				Bfree(b);
532 				b = b1;
533 				}
534 			if ( (j = b5 - m5) !=0)
535 				b = pow5mult(b, j);
536 			}
537 		else
538 			b = pow5mult(b, b5);
539 		}
540 	S = i2b(1);
541 	if (s5 > 0)
542 		S = pow5mult(S, s5);
543 
544 	/* Check for special case that d is a normalized power of 2. */
545 
546 	spec_case = 0;
547 	if (mode < 2) {
548 		if (bbits == 1 && be0 > fpi->emin + 1) {
549 			/* The special case */
550 			b2++;
551 			s2++;
552 			spec_case = 1;
553 			}
554 		}
555 
556 	/* Arrange for convenient computation of quotients:
557 	 * shift left if necessary so divisor has 4 leading 0 bits.
558 	 *
559 	 * Perhaps we should just compute leading 28 bits of S once
560 	 * and for all and pass them and a shift to quorem, so it
561 	 * can do shifts and ors to compute the numerator for q.
562 	 */
563 #ifdef Pack_32
564 	if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) !=0)
565 		i = 32 - i;
566 #else
567 	if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) !=0)
568 		i = 16 - i;
569 #endif
570 	if (i > 4) {
571 		i -= 4;
572 		b2 += i;
573 		m2 += i;
574 		s2 += i;
575 		}
576 	else if (i < 4) {
577 		i += 28;
578 		b2 += i;
579 		m2 += i;
580 		s2 += i;
581 		}
582 	if (b2 > 0)
583 		b = lshift(b, b2);
584 	if (s2 > 0)
585 		S = lshift(S, s2);
586 	if (k_check) {
587 		if (cmp(b,S) < 0) {
588 			k--;
589 			b = multadd(b, 10, 0);	/* we botched the k estimate */
590 			if (leftright)
591 				mhi = multadd(mhi, 10, 0);
592 			ilim = ilim1;
593 			}
594 		}
595 	if (ilim <= 0 && mode > 2) {
596 		if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
597 			/* no digits, fcvt style */
598  no_digits:
599 			k = -1 - ndigits;
600 			inex = STRTOG_Inexlo;
601 			goto ret;
602 			}
603  one_digit:
604 		inex = STRTOG_Inexhi;
605 		*s++ = '1';
606 		k++;
607 		goto ret;
608 		}
609 	if (leftright) {
610 		if (m2 > 0)
611 			mhi = lshift(mhi, m2);
612 
613 		/* Compute mlo -- check for special case
614 		 * that d is a normalized power of 2.
615 		 */
616 
617 		mlo = mhi;
618 		if (spec_case) {
619 			mhi = Balloc(mhi->k);
620 			Bcopy(mhi, mlo);
621 			mhi = lshift(mhi, 1);
622 			}
623 
624 		for(i = 1;;i++) {
625 			dig = quorem(b,S) + '0';
626 			/* Do we yet have the shortest decimal string
627 			 * that will round to d?
628 			 */
629 			j = cmp(b, mlo);
630 			delta = diff(S, mhi);
631 			jj1 = delta->sign ? 1 : cmp(b, delta);
632 			Bfree(delta);
633 #ifndef ROUND_BIASED
634 			if (jj1 == 0 && !mode && !(bits[0] & 1) && !rdir) {
635 				if (dig == '9')
636 					goto round_9_up;
637 				if (j <= 0) {
638 					if (b->wds > 1 || b->x[0])
639 						inex = STRTOG_Inexlo;
640 					}
641 				else {
642 					dig++;
643 					inex = STRTOG_Inexhi;
644 					}
645 				*s++ = dig;
646 				goto ret;
647 				}
648 #endif
649 			if (j < 0 || (j == 0 && !mode
650 #ifndef ROUND_BIASED
651 							&& !(bits[0] & 1)
652 #endif
653 					)) {
654 				if (rdir && (b->wds > 1 || b->x[0])) {
655 					if (rdir == 2) {
656 						inex = STRTOG_Inexlo;
657 						goto accept;
658 						}
659 					while (cmp(S,mhi) > 0) {
660 						*s++ = dig;
661 						mhi1 = multadd(mhi, 10, 0);
662 						if (mlo == mhi)
663 							mlo = mhi1;
664 						mhi = mhi1;
665 						b = multadd(b, 10, 0);
666 						dig = quorem(b,S) + '0';
667 						}
668 					if (dig++ == '9')
669 						goto round_9_up;
670 					inex = STRTOG_Inexhi;
671 					goto accept;
672 					}
673 				if (jj1 > 0) {
674 					b = lshift(b, 1);
675 					jj1 = cmp(b, S);
676 					if ((jj1 > 0 || (jj1 == 0 && dig & 1))
677 					&& dig++ == '9')
678 						goto round_9_up;
679 					inex = STRTOG_Inexhi;
680 					}
681 				if (b->wds > 1 || b->x[0])
682 					inex = STRTOG_Inexlo;
683  accept:
684 				*s++ = dig;
685 				goto ret;
686 				}
687 			if (jj1 > 0 && rdir != 2) {
688 				if (dig == '9') { /* possible if i == 1 */
689  round_9_up:
690 					*s++ = '9';
691 					inex = STRTOG_Inexhi;
692 					goto roundoff;
693 					}
694 				inex = STRTOG_Inexhi;
695 				*s++ = dig + 1;
696 				goto ret;
697 				}
698 			*s++ = dig;
699 			if (i == ilim)
700 				break;
701 			b = multadd(b, 10, 0);
702 			if (mlo == mhi)
703 				mlo = mhi = multadd(mhi, 10, 0);
704 			else {
705 				mlo = multadd(mlo, 10, 0);
706 				mhi = multadd(mhi, 10, 0);
707 				}
708 			}
709 		}
710 	else
711 		for(i = 1;; i++) {
712 			*s++ = dig = quorem(b,S) + '0';
713 			if (i >= ilim)
714 				break;
715 			b = multadd(b, 10, 0);
716 			}
717 
718 	/* Round off last digit */
719 
720 	if (rdir) {
721 		if (rdir == 2 || (b->wds <= 1 && !b->x[0]))
722 			goto chopzeros;
723 		goto roundoff;
724 		}
725 	b = lshift(b, 1);
726 	j = cmp(b, S);
727 	if (j > 0 || (j == 0 && dig & 1)) {
728  roundoff:
729 		inex = STRTOG_Inexhi;
730 		while(*--s == '9')
731 			if (s == s0) {
732 				k++;
733 				*s++ = '1';
734 				goto ret;
735 				}
736 		++*s++;
737 		}
738 	else {
739  chopzeros:
740 		if (b->wds > 1 || b->x[0])
741 			inex = STRTOG_Inexlo;
742 		while(*--s == '0'){}
743 		s++;
744 		}
745  ret:
746 	Bfree(S);
747 	if (mhi) {
748 		if (mlo && mlo != mhi)
749 			Bfree(mlo);
750 		Bfree(mhi);
751 		}
752  ret1:
753 	Bfree(b);
754 	*s = 0;
755 	*decpt = k + 1;
756 	if (rve)
757 		*rve = s;
758 	*kindp |= inex;
759 	return s0;
760 	}
761