xref: /netbsd-src/lib/libc/gdtoa/dtoa.c (revision 274254cdae52594c1aa480a736aef78313d15c9c)
1 /* $NetBSD: dtoa.c,v 1.5 2008/03/21 23:13:48 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 /* 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 	if (b == NULL)
196 		return NULL;
197 #ifdef Sudden_Underflow
198 	i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
199 #else
200 	if (( i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)) )!=0) {
201 #endif
202 		dval(d2) = dval(d);
203 		word0(d2) &= Frac_mask1;
204 		word0(d2) |= Exp_11;
205 #ifdef IBM
206 		if (( j = 11 - hi0bits(word0(d2) & Frac_mask) )!=0)
207 			dval(d2) /= 1 << j;
208 #endif
209 
210 		/* log(x)	~=~ log(1.5) + (x-1.5)/1.5
211 		 * log10(x)	 =  log(x) / log(10)
212 		 *		~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
213 		 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
214 		 *
215 		 * This suggests computing an approximation k to log10(d) by
216 		 *
217 		 * k = (i - Bias)*0.301029995663981
218 		 *	+ ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
219 		 *
220 		 * We want k to be too large rather than too small.
221 		 * The error in the first-order Taylor series approximation
222 		 * is in our favor, so we just round up the constant enough
223 		 * to compensate for any error in the multiplication of
224 		 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
225 		 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
226 		 * adding 1e-13 to the constant term more than suffices.
227 		 * Hence we adjust the constant term to 0.1760912590558.
228 		 * (We could get a more accurate k by invoking log10,
229 		 *  but this is probably not worthwhile.)
230 		 */
231 
232 		i -= Bias;
233 #ifdef IBM
234 		i <<= 2;
235 		i += j;
236 #endif
237 #ifndef Sudden_Underflow
238 		denorm = 0;
239 		}
240 	else {
241 		/* d is denormalized */
242 
243 		i = bbits + be + (Bias + (P-1) - 1);
244 		x = i > 32  ? word0(d) << (64 - i) | word1(d) >> (i - 32)
245 			    : word1(d) << (32 - i);
246 		dval(d2) = x;
247 		word0(d2) -= 31*Exp_msk1; /* adjust exponent */
248 		i -= (Bias + (P-1) - 1) + 1;
249 		denorm = 1;
250 		}
251 #endif
252 	ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
253 	k = (int)ds;
254 	if (ds < 0. && ds != k)
255 		k--;	/* want k = floor(ds) */
256 	k_check = 1;
257 	if (k >= 0 && k <= Ten_pmax) {
258 		if (dval(d) < tens[k])
259 			k--;
260 		k_check = 0;
261 		}
262 	j = bbits - i - 1;
263 	if (j >= 0) {
264 		b2 = 0;
265 		s2 = j;
266 		}
267 	else {
268 		b2 = -j;
269 		s2 = 0;
270 		}
271 	if (k >= 0) {
272 		b5 = 0;
273 		s5 = k;
274 		s2 += k;
275 		}
276 	else {
277 		b2 -= k;
278 		b5 = -k;
279 		s5 = 0;
280 		}
281 	if (mode < 0 || mode > 9)
282 		mode = 0;
283 
284 #ifndef SET_INEXACT
285 #ifdef Check_FLT_ROUNDS
286 	try_quick = Rounding == 1;
287 #else
288 	try_quick = 1;
289 #endif
290 #endif /*SET_INEXACT*/
291 
292 	if (mode > 5) {
293 		mode -= 4;
294 		try_quick = 0;
295 		}
296 	leftright = 1;
297 	switch(mode) {
298 		case 0:
299 		case 1:
300 			ilim = ilim1 = -1;
301 			i = 18;
302 			ndigits = 0;
303 			break;
304 		case 2:
305 			leftright = 0;
306 			/* FALLTHROUGH */
307 		case 4:
308 			if (ndigits <= 0)
309 				ndigits = 1;
310 			ilim = ilim1 = i = ndigits;
311 			break;
312 		case 3:
313 			leftright = 0;
314 			/* FALLTHROUGH */
315 		case 5:
316 			i = ndigits + k + 1;
317 			ilim = i;
318 			ilim1 = i - 1;
319 			if (i <= 0)
320 				i = 1;
321 		}
322 	s = s0 = rv_alloc((size_t)i);
323 	if (s == NULL)
324 		return NULL;
325 
326 #ifdef Honor_FLT_ROUNDS
327 	if (mode > 1 && rounding != 1)
328 		leftright = 0;
329 #endif
330 
331 	if (ilim >= 0 && ilim <= Quick_max && try_quick) {
332 
333 		/* Try to get by with floating-point arithmetic. */
334 
335 		i = 0;
336 		dval(d2) = dval(d);
337 		k0 = k;
338 		ilim0 = ilim;
339 		ieps = 2; /* conservative */
340 		if (k > 0) {
341 			ds = tens[k&0xf];
342 			j = (unsigned int)k >> 4;
343 			if (j & Bletch) {
344 				/* prevent overflows */
345 				j &= Bletch - 1;
346 				dval(d) /= bigtens[n_bigtens-1];
347 				ieps++;
348 				}
349 			for(; j; j = (unsigned int)j >> 1, i++)
350 				if (j & 1) {
351 					ieps++;
352 					ds *= bigtens[i];
353 					}
354 			dval(d) /= ds;
355 			}
356 		else if (( jj1 = -k )!=0) {
357 			dval(d) *= tens[jj1 & 0xf];
358 			for(j = jj1 >> 4; j; j >>= 1, i++)
359 				if (j & 1) {
360 					ieps++;
361 					dval(d) *= bigtens[i];
362 					}
363 			}
364 		if (k_check && dval(d) < 1. && ilim > 0) {
365 			if (ilim1 <= 0)
366 				goto fast_failed;
367 			ilim = ilim1;
368 			k--;
369 			dval(d) *= 10.;
370 			ieps++;
371 			}
372 		dval(eps) = ieps*dval(d) + 7.;
373 		word0(eps) -= (P-1)*Exp_msk1;
374 		if (ilim == 0) {
375 			S = mhi = 0;
376 			dval(d) -= 5.;
377 			if (dval(d) > dval(eps))
378 				goto one_digit;
379 			if (dval(d) < -dval(eps))
380 				goto no_digits;
381 			goto fast_failed;
382 			}
383 #ifndef No_leftright
384 		if (leftright) {
385 			/* Use Steele & White method of only
386 			 * generating digits needed.
387 			 */
388 			dval(eps) = 0.5/tens[ilim-1] - dval(eps);
389 			for(i = 0;;) {
390 				L = dval(d);
391 				dval(d) -= L;
392 				*s++ = '0' + (int)L;
393 				if (dval(d) < dval(eps))
394 					goto ret1;
395 				if (1. - dval(d) < dval(eps))
396 					goto bump_up;
397 				if (++i >= ilim)
398 					break;
399 				dval(eps) *= 10.;
400 				dval(d) *= 10.;
401 				}
402 			}
403 		else {
404 #endif
405 			/* Generate ilim digits, then fix them up. */
406 			dval(eps) *= tens[ilim-1];
407 			for(i = 1;; i++, dval(d) *= 10.) {
408 				L = (Long)(dval(d));
409 				if (!(dval(d) -= L))
410 					ilim = i;
411 				*s++ = '0' + (int)L;
412 				if (i == ilim) {
413 					if (dval(d) > 0.5 + dval(eps))
414 						goto bump_up;
415 					else if (dval(d) < 0.5 - dval(eps)) {
416 						while(*--s == '0');
417 						s++;
418 						goto ret1;
419 						}
420 					break;
421 					}
422 				}
423 #ifndef No_leftright
424 			}
425 #endif
426  fast_failed:
427 		s = s0;
428 		dval(d) = dval(d2);
429 		k = k0;
430 		ilim = ilim0;
431 		}
432 
433 	/* Do we have a "small" integer? */
434 
435 	if (be >= 0 && k <= Int_max) {
436 		/* Yes. */
437 		ds = tens[k];
438 		if (ndigits < 0 && ilim <= 0) {
439 			S = mhi = 0;
440 			if (ilim < 0 || dval(d) <= 5*ds)
441 				goto no_digits;
442 			goto one_digit;
443 			}
444 		for(i = 1;; i++, dval(d) *= 10.) {
445 			L = (Long)(dval(d) / ds);
446 			dval(d) -= L*ds;
447 #ifdef Check_FLT_ROUNDS
448 			/* If FLT_ROUNDS == 2, L will usually be high by 1 */
449 			if (dval(d) < 0) {
450 				L--;
451 				dval(d) += ds;
452 				}
453 #endif
454 			*s++ = '0' + (int)L;
455 			if (!dval(d)) {
456 #ifdef SET_INEXACT
457 				inexact = 0;
458 #endif
459 				break;
460 				}
461 			if (i == ilim) {
462 #ifdef Honor_FLT_ROUNDS
463 				if (mode > 1)
464 				switch(rounding) {
465 				  case 0: goto ret1;
466 				  case 2: goto bump_up;
467 				  }
468 #endif
469 				dval(d) += dval(d);
470 				if (dval(d) > ds || (dval(d) == ds && L & 1)) {
471  bump_up:
472 					while(*--s == '9')
473 						if (s == s0) {
474 							k++;
475 							*s = '0';
476 							break;
477 							}
478 					++*s++;
479 					}
480 				break;
481 				}
482 			}
483 		goto ret1;
484 		}
485 
486 	m2 = b2;
487 	m5 = b5;
488 	mhi = mlo = 0;
489 	if (leftright) {
490 		i =
491 #ifndef Sudden_Underflow
492 			denorm ? be + (Bias + (P-1) - 1 + 1) :
493 #endif
494 #ifdef IBM
495 			1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
496 #else
497 			1 + P - bbits;
498 #endif
499 		b2 += i;
500 		s2 += i;
501 		mhi = i2b(1);
502 		if (mhi == NULL)
503 			return NULL;
504 		}
505 	if (m2 > 0 && s2 > 0) {
506 		i = m2 < s2 ? m2 : s2;
507 		b2 -= i;
508 		m2 -= i;
509 		s2 -= i;
510 		}
511 	if (b5 > 0) {
512 		if (leftright) {
513 			if (m5 > 0) {
514 				mhi = pow5mult(mhi, m5);
515 				if (mhi == NULL)
516 					return NULL;
517 				b1 = mult(mhi, b);
518 				if (b1 == NULL)
519 					return NULL;
520 				Bfree(b);
521 				b = b1;
522 				}
523 			if (( j = b5 - m5 )!=0)
524 				b = pow5mult(b, j);
525 				if (b == NULL)
526 					return NULL;
527 			}
528 		else
529 			b = pow5mult(b, b5);
530 			if (b == NULL)
531 				return NULL;
532 		}
533 	S = i2b(1);
534 	if (S == NULL)
535 		return NULL;
536 	if (s5 > 0) {
537 		S = pow5mult(S, s5);
538 		if (S == NULL)
539 			return NULL;
540 	}
541 
542 	/* Check for special case that d is a normalized power of 2. */
543 
544 	spec_case = 0;
545 	if ((mode < 2 || leftright)
546 #ifdef Honor_FLT_ROUNDS
547 			&& rounding == 1
548 #endif
549 				) {
550 		if (!word1(d) && !(word0(d) & Bndry_mask)
551 #ifndef Sudden_Underflow
552 		 && word0(d) & (Exp_mask & ~Exp_msk1)
553 #endif
554 				) {
555 			/* The special case */
556 			b2 += Log2P;
557 			s2 += Log2P;
558 			spec_case = 1;
559 			}
560 		}
561 
562 	/* Arrange for convenient computation of quotients:
563 	 * shift left if necessary so divisor has 4 leading 0 bits.
564 	 *
565 	 * Perhaps we should just compute leading 28 bits of S once
566 	 * and for all and pass them and a shift to quorem, so it
567 	 * can do shifts and ors to compute the numerator for q.
568 	 */
569 #ifdef Pack_32
570 	if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f )!=0)
571 		i = 32 - i;
572 #else
573 	if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf )!=0)
574 		i = 16 - i;
575 #endif
576 	if (i > 4) {
577 		i -= 4;
578 		b2 += i;
579 		m2 += i;
580 		s2 += i;
581 		}
582 	else if (i < 4) {
583 		i += 28;
584 		b2 += i;
585 		m2 += i;
586 		s2 += i;
587 		}
588 	if (b2 > 0) {
589 		b = lshift(b, b2);
590 		if (b == NULL)
591 			return NULL;
592 	}
593 	if (s2 > 0) {
594 		S = lshift(S, s2);
595 		if (S == NULL)
596 			return NULL;
597 	}
598 	if (k_check) {
599 		if (cmp(b,S) < 0) {
600 			k--;
601 			b = multadd(b, 10, 0);	/* we botched the k estimate */
602 			if (b == NULL)
603 				return NULL;
604 			if (leftright) {
605 				mhi = multadd(mhi, 10, 0);
606 				if (mhi == NULL)
607 					return NULL;
608 			}
609 			ilim = ilim1;
610 			}
611 		}
612 	if (ilim <= 0 && (mode == 3 || mode == 5)) {
613 		if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
614 			/* no digits, fcvt style */
615  no_digits:
616 			k = -1 - ndigits;
617 			goto ret;
618 			}
619  one_digit:
620 		*s++ = '1';
621 		k++;
622 		goto ret;
623 		}
624 	if (leftright) {
625 		if (m2 > 0) {
626 			mhi = lshift(mhi, m2);
627 			if (mhi == NULL)
628 				return NULL;
629 		}
630 
631 		/* Compute mlo -- check for special case
632 		 * that d is a normalized power of 2.
633 		 */
634 
635 		mlo = mhi;
636 		if (spec_case) {
637 			mhi = Balloc(mhi->k);
638 			if (mhi == NULL)
639 				return NULL;
640 			Bcopy(mhi, mlo);
641 			mhi = lshift(mhi, Log2P);
642 			if (mhi == NULL)
643 				return NULL;
644 			}
645 
646 		for(i = 1;;i++) {
647 			dig = quorem(b,S) + '0';
648 			/* Do we yet have the shortest decimal string
649 			 * that will round to d?
650 			 */
651 			j = cmp(b, mlo);
652 			delta = diff(S, mhi);
653 			if (delta == NULL)
654 				return NULL;
655 			jj1 = delta->sign ? 1 : cmp(b, delta);
656 			Bfree(delta);
657 #ifndef ROUND_BIASED
658 			if (jj1 == 0 && mode != 1 && !(word1(d) & 1)
659 #ifdef Honor_FLT_ROUNDS
660 				&& rounding >= 1
661 #endif
662 								   ) {
663 				if (dig == '9')
664 					goto round_9_up;
665 				if (j > 0)
666 					dig++;
667 #ifdef SET_INEXACT
668 				else if (!b->x[0] && b->wds <= 1)
669 					inexact = 0;
670 #endif
671 				*s++ = dig;
672 				goto ret;
673 				}
674 #endif
675 			if (j < 0 || (j == 0 && mode != 1
676 #ifndef ROUND_BIASED
677 							&& !(word1(d) & 1)
678 #endif
679 					)) {
680 				if (!b->x[0] && b->wds <= 1) {
681 #ifdef SET_INEXACT
682 					inexact = 0;
683 #endif
684 					goto accept_dig;
685 					}
686 #ifdef Honor_FLT_ROUNDS
687 				if (mode > 1)
688 				 switch(rounding) {
689 				  case 0: goto accept_dig;
690 				  case 2: goto keep_dig;
691 				  }
692 #endif /*Honor_FLT_ROUNDS*/
693 				if (jj1 > 0) {
694 					b = lshift(b, 1);
695 					if (b == NULL)
696 						return NULL;
697 					jj1 = cmp(b, S);
698 					if ((jj1 > 0 || (jj1 == 0 && dig & 1))
699 					&& dig++ == '9')
700 						goto round_9_up;
701 					}
702  accept_dig:
703 				*s++ = dig;
704 				goto ret;
705 				}
706 			if (jj1 > 0) {
707 #ifdef Honor_FLT_ROUNDS
708 				if (!rounding)
709 					goto accept_dig;
710 #endif
711 				if (dig == '9') { /* possible if i == 1 */
712  round_9_up:
713 					*s++ = '9';
714 					goto roundoff;
715 					}
716 				*s++ = dig + 1;
717 				goto ret;
718 				}
719 #ifdef Honor_FLT_ROUNDS
720  keep_dig:
721 #endif
722 			*s++ = dig;
723 			if (i == ilim)
724 				break;
725 			b = multadd(b, 10, 0);
726 			if (b == NULL)
727 				return NULL;
728 			if (mlo == mhi) {
729 				mlo = mhi = multadd(mhi, 10, 0);
730 				if (mlo == NULL)
731 					return NULL;
732 				}
733 			else {
734 				mlo = multadd(mlo, 10, 0);
735 				if (mlo == NULL)
736 					return NULL;
737 				mhi = multadd(mhi, 10, 0);
738 				if (mhi == NULL)
739 					return NULL;
740 				}
741 			}
742 		}
743 	else
744 		for(i = 1;; i++) {
745 			*s++ = dig = quorem(b,S) + '0';
746 			if (!b->x[0] && b->wds <= 1) {
747 #ifdef SET_INEXACT
748 				inexact = 0;
749 #endif
750 				goto ret;
751 				}
752 			if (i >= ilim)
753 				break;
754 			b = multadd(b, 10, 0);
755 			if (b == NULL)
756 				return NULL;
757 			}
758 
759 	/* Round off last digit */
760 
761 #ifdef Honor_FLT_ROUNDS
762 	switch(rounding) {
763 	  case 0: goto trimzeros;
764 	  case 2: goto roundoff;
765 	  }
766 #endif
767 	b = lshift(b, 1);
768 	j = cmp(b, S);
769 	if (j > 0 || (j == 0 && dig & 1)) {
770  roundoff:
771 		while(*--s == '9')
772 			if (s == s0) {
773 				k++;
774 				*s++ = '1';
775 				goto ret;
776 				}
777 		++*s++;
778 		}
779 	else {
780 #ifdef Honor_FLT_ROUNDS
781  trimzeros:
782 #endif
783 		while(*--s == '0');
784 		s++;
785 		}
786  ret:
787 	Bfree(S);
788 	if (mhi) {
789 		if (mlo && mlo != mhi)
790 			Bfree(mlo);
791 		Bfree(mhi);
792 		}
793  ret1:
794 #ifdef SET_INEXACT
795 	if (inexact) {
796 		if (!oldinexact) {
797 			word0(d) = Exp_1 + (70 << Exp_shift);
798 			word1(d) = 0;
799 			dval(d) += 1.;
800 			}
801 		}
802 	else if (!oldinexact)
803 		clear_inexact();
804 #endif
805 	Bfree(b);
806 	if (s == s0) {			/* don't return empty string */
807 		*s++ = '0';
808 		k = 0;
809 	}
810 	*s = 0;
811 	*decpt = k + 1;
812 	if (rve)
813 		*rve = s;
814 	return s0;
815 	}
816