Home | History | Annotate | Line # | Download | only in gdtoa
dtoa.c revision 1.4
      1 /* $NetBSD: dtoa.c,v 1.4 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 /* 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((size_t)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