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