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