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