[apple/libc.git] / gdtoa / FreeBSD / gdtoa-strtodg.c

/****************************************************************

The author of this software is David M. Gay.

Copyright (C) 1998-2001 by Lucent Technologies
All Rights Reserved

Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.

LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.

****************************************************************/

/* Please send bug reports to David M. Gay (dmg at acm dot org,
 * with " at " changed at "@" and " dot " changed to ".").	*/

#include "gdtoaimp.h"

#ifdef USE_LOCALE
#include "locale.h"
#endif

 static CONST int
fivesbits[] = {	 0,  3,  5,  7, 10, 12, 14, 17, 19, 21,
		24, 26, 28, 31, 33, 35, 38, 40, 42, 45,
		47, 49, 52
#ifdef VAX
		, 54, 56
#endif
		};

 Bigint *
#ifdef KR_headers
increment(b) Bigint *b;
#else
increment(Bigint *b)
#endif
{
	ULong *x, *xe;
	Bigint *b1;
#ifdef Pack_16
	ULong carry = 1, y;
#endif

	x = b->x;
	xe = x + b->wds;
#ifdef Pack_32
	do {
		if (*x < (ULong)0xffffffffL) {
			++*x;
			return b;
			}
		*x++ = 0;
		} while(x < xe);
#else
	do {
		y = *x + carry;
		carry = y >> 16;
		*x++ = y & 0xffff;
		if (!carry)
			return b;
		} while(x < xe);
	if (carry)
#endif
	{
		if (b->wds >= b->maxwds) {
			b1 = Balloc(b->k+1);
			Bcopy(b1,b);
			Bfree(b);
			b = b1;
			}
		b->x[b->wds++] = 1;
		}
	return b;
	}

 int
#ifdef KR_headers
decrement(b) Bigint *b;
#else
decrement(Bigint *b)
#endif
{
	ULong *x, *xe;
#ifdef Pack_16
	ULong borrow = 1, y;
#endif

	x = b->x;
	xe = x + b->wds;
#ifdef Pack_32
	do {
		if (*x) {
			--*x;
			break;
			}
		*x++ = 0xffffffffL;
		}
		while(x < xe);
#else
	do {
		y = *x - borrow;
		borrow = (y & 0x10000) >> 16;
		*x++ = y & 0xffff;
		} while(borrow && x < xe);
#endif
	return STRTOG_Inexlo;
	}

 static int
#ifdef KR_headers
all_on(b, n) Bigint *b; int n;
#else
all_on(Bigint *b, int n)
#endif
{
	ULong *x, *xe;

	x = b->x;
	xe = x + (n >> kshift);
	while(x < xe)
		if ((*x++ & ALL_ON) != ALL_ON)
			return 0;
	if (n &= kmask)
		return ((*x | (ALL_ON << n)) & ALL_ON) == ALL_ON;
	return 1;
	}

 Bigint *
#ifdef KR_headers
set_ones(b, n) Bigint *b; int n;
#else
set_ones(Bigint *b, int n)
#endif
{
	int k;
	ULong *x, *xe;

	k = (n + ((1 << kshift) - 1)) >> kshift;
	if (b->k < k) {
		Bfree(b);
		b = Balloc(k);
		}
	k = n >> kshift;
	if (n &= kmask)
		k++;
	b->wds = k;
	x = b->x;
	xe = x + k;
	while(x < xe)
		*x++ = ALL_ON;
	if (n)
		x[-1] >>= ULbits - n;
	return b;
	}

 static int
rvOK
#ifdef KR_headers
 (d, fpi, exp, bits, exact, rd, irv)
 double d; FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv;
#else
 (double d, FPI *fpi, Long *exp, ULong *bits, int exact, int rd, int *irv)
#endif
{
	Bigint *b;
	ULong carry, inex, lostbits;
	int bdif, e, j, k, k1, nb, rv;

	carry = rv = 0;
	b = d2b(d, &e, &bdif);
	bdif -= nb = fpi->nbits;
	e += bdif;
	if (bdif <= 0) {
		if (exact)
			goto trunc;
		goto ret;
		}
	if (P == nb) {
		if (
#ifndef IMPRECISE_INEXACT
			exact &&
#endif
			fpi->rounding ==
#ifdef RND_PRODQUOT
					FPI_Round_near
#else
					Flt_Rounds
#endif
			) goto trunc;
		goto ret;
		}
	switch(rd) {
	  case 1:
		goto trunc;
	  case 2:
		break;
	  default: /* round near */
		k = bdif - 1;
		if (k < 0)
			goto trunc;
		if (!k) {
			if (!exact)
				goto ret;
			if (b->x[0] & 2)
				break;
			goto trunc;
			}
		if (b->x[k>>kshift] & ((ULong)1 << (k & kmask)))
			break;
		goto trunc;
	  }
	/* "break" cases: round up 1 bit, then truncate; bdif > 0 */
	carry = 1;
 trunc:
	inex = lostbits = 0;
	if (bdif > 0) {
		if ( (lostbits = any_on(b, bdif)) !=0)
			inex = STRTOG_Inexlo;
		rshift(b, bdif);
		if (carry) {
			inex = STRTOG_Inexhi;
			b = increment(b);
			if ( (j = nb & kmask) !=0)
				j = ULbits - j;
			if (hi0bits(b->x[b->wds - 1]) != j) {
				if (!lostbits)
					lostbits = b->x[0] & 1;
				rshift(b, 1);
				e++;
				}
			}
		}
	else if (bdif < 0)
		b = lshift(b, -bdif);
	if (e < fpi->emin) {
		k = fpi->emin - e;
		e = fpi->emin;
		if (k > nb || fpi->sudden_underflow) {
			b->wds = inex = 0;
			*irv = STRTOG_Underflow | STRTOG_Inexlo;
			}
		else {
			k1 = k - 1;
			if (k1 > 0 && !lostbits)
				lostbits = any_on(b, k1);
			if (!lostbits && !exact)
				goto ret;
			lostbits |=
			  carry = b->x[k1>>kshift] & (1 << (k1 & kmask));
			rshift(b, k);
			*irv = STRTOG_Denormal;
			if (carry) {
				b = increment(b);
				inex = STRTOG_Inexhi | STRTOG_Underflow;
				}
			else if (lostbits)
				inex = STRTOG_Inexlo | STRTOG_Underflow;
			}
		}
	else if (e > fpi->emax) {
		e = fpi->emax + 1;
		*irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
#ifndef NO_ERRNO
		errno = ERANGE;
#endif
		b->wds = inex = 0;
		}
	*exp = e;
	copybits(bits, nb, b);
	*irv |= inex;
	rv = 1;
 ret:
	Bfree(b);
	return rv;
	}

 static int
#ifdef KR_headers
mantbits(d) double d;
#else
mantbits(double d)
#endif
{
	ULong L;
#ifdef VAX
	L = word1(d) << 16 | word1(d) >> 16;
	if (L)
#else
	if ( (L = word1(d)) !=0)
#endif
		return P - lo0bits(&L);
#ifdef VAX
	L = word0(d) << 16 | word0(d) >> 16 | Exp_msk11;
#else
	L = word0(d) | Exp_msk1;
#endif
	return P - 32 - lo0bits(&L);
	}

 int
strtodg
#ifdef KR_headers
	(s00, se, fpi, exp, bits)
	CONST char *s00; char **se; FPI *fpi; Long *exp; ULong *bits;
#else
	(CONST char *s00, char **se, FPI *fpi, Long *exp, ULong *bits)
#endif
{
	int abe, abits, asub;
	int bb0, bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, denorm;
	int dsign, e, e1, e2, emin, esign, finished, i, inex, irv;
	int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign;
	int sudden_underflow;
	CONST char *s, *s0, *s1;
	double adj, adj0, rv, tol;
	Long L;
	ULong y, z;
	Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0;

	irv = STRTOG_Zero;
	denorm = sign = nz0 = nz = 0;
	dval(rv) = 0.;
	rvb = 0;
	nbits = fpi->nbits;
	for(s = s00;;s++) switch(*s) {
		case '-':
			sign = 1;
			/* no break */
		case '+':
			if (*++s)
				goto break2;
			/* no break */
		case 0:
			sign = 0;
			irv = STRTOG_NoNumber;
			s = s00;
			goto ret;
		case '\t':
		case '\n':
		case '\v':
		case '\f':
		case '\r':
		case ' ':
			continue;
		default:
			goto break2;
		}
 break2:
	if (*s == '0') {
#ifndef NO_HEX_FP
		switch(s[1]) {
		  case 'x':
		  case 'X':
			irv = gethex(&s, fpi, exp, &rvb, sign);
			if (irv == STRTOG_NoNumber) {
				s = s00;
				sign = 0;
				}
			goto ret;
		  }
#endif
		nz0 = 1;
		while(*++s == '0') ;
		if (!*s)
			goto ret;
		}
	sudden_underflow = fpi->sudden_underflow;
	s0 = s;
	y = z = 0;
	for(decpt = nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
		if (nd < 9)
			y = 10*y + c - '0';
		else if (nd < 16)
			z = 10*z + c - '0';
	nd0 = nd;
#ifdef USE_LOCALE
	if (c == *localeconv()->decimal_point)
#else
	if (c == '.')
#endif
		{
		decpt = 1;
		c = *++s;
		if (!nd) {
			for(; c == '0'; c = *++s)
				nz++;
			if (c > '0' && c <= '9') {
				s0 = s;
				nf += nz;
				nz = 0;
				goto have_dig;
				}
			goto dig_done;
			}
		for(; c >= '0' && c <= '9'; c = *++s) {
 have_dig:
			nz++;
			if (c -= '0') {
				nf += nz;
				for(i = 1; i < nz; i++)
					if (nd++ < 9)
						y *= 10;
					else if (nd <= DBL_DIG + 1)
						z *= 10;
				if (nd++ < 9)
					y = 10*y + c;
				else if (nd <= DBL_DIG + 1)
					z = 10*z + c;
				nz = 0;
				}
			}
		}
 dig_done:
	e = 0;
	if (c == 'e' || c == 'E') {
		if (!nd && !nz && !nz0) {
			irv = STRTOG_NoNumber;
			s = s00;
			goto ret;
			}
		s00 = s;
		esign = 0;
		switch(c = *++s) {
			case '-':
				esign = 1;
			case '+':
				c = *++s;
			}
		if (c >= '0' && c <= '9') {
			while(c == '0')
				c = *++s;
			if (c > '0' && c <= '9') {
				L = c - '0';
				s1 = s;
				while((c = *++s) >= '0' && c <= '9')
					L = 10*L + c - '0';
				if (s - s1 > 8 || L > 19999)
					/* Avoid confusion from exponents
					 * so large that e might overflow.
					 */
					e = 19999; /* safe for 16 bit ints */
				else
					e = (int)L;
				if (esign)
					e = -e;
				}
			else
				e = 0;
			}
		else
			s = s00;
		}
	if (!nd) {
		if (!nz && !nz0) {
#ifdef INFNAN_CHECK
			/* Check for Nan and Infinity */
			if (!decpt)
			 switch(c) {
			  case 'i':
			  case 'I':
				if (match(&s,"nf")) {
					--s;
					if (!match(&s,"inity"))
						++s;
					irv = STRTOG_Infinite;
					goto infnanexp;
					}
				break;
			  case 'n':
			  case 'N':
				if (match(&s, "an")) {
					irv = STRTOG_NaN;
					*exp = fpi->emax + 1;
#ifndef No_Hex_NaN
					if (*s == '(') /*)*/
						irv = hexnan(&s, fpi, bits);
#endif
					goto infnanexp;
					}
			  }
#endif /* INFNAN_CHECK */
			irv = STRTOG_NoNumber;
			s = s00;
			}
		goto ret;
		}

	irv = STRTOG_Normal;
	e1 = e -= nf;
	rd = 0;
	switch(fpi->rounding & 3) {
	  case FPI_Round_up:
		rd = 2 - sign;
		break;
	  case FPI_Round_zero:
		rd = 1;
		break;
	  case FPI_Round_down:
		rd = 1 + sign;
	  }

	/* Now we have nd0 digits, starting at s0, followed by a
	 * decimal point, followed by nd-nd0 digits.  The number we're
	 * after is the integer represented by those digits times
	 * 10**e */

	if (!nd0)
		nd0 = nd;
	k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
	dval(rv) = y;
	if (k > 9)
		dval(rv) = tens[k - 9] * dval(rv) + z;
	bd0 = 0;
	if (nbits <= P && nd <= DBL_DIG) {
		if (!e) {
			if (rvOK(dval(rv), fpi, exp, bits, 1, rd, &irv))
				goto ret;
			}
		else if (e > 0) {
			if (e <= Ten_pmax) {
#ifdef VAX
				goto vax_ovfl_check;
#else
				i = fivesbits[e] + mantbits(dval(rv)) <= P;
				/* rv = */ rounded_product(dval(rv), tens[e]);
				if (rvOK(dval(rv), fpi, exp, bits, i, rd, &irv))
					goto ret;
				e1 -= e;
				goto rv_notOK;
#endif
				}
			i = DBL_DIG - nd;
			if (e <= Ten_pmax + i) {
				/* A fancier test would sometimes let us do
				 * this for larger i values.
				 */
				e2 = e - i;
				e1 -= i;
				dval(rv) *= tens[i];
#ifdef VAX
				/* VAX exponent range is so narrow we must
				 * worry about overflow here...
				 */
 vax_ovfl_check:
				dval(adj) = dval(rv);
				word0(adj) -= P*Exp_msk1;
				/* adj = */ rounded_product(dval(adj), tens[e2]);
				if ((word0(adj) & Exp_mask)
				 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
					goto rv_notOK;
				word0(adj) += P*Exp_msk1;
				dval(rv) = dval(adj);
#else
				/* rv = */ rounded_product(dval(rv), tens[e2]);
#endif
				if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv))
					goto ret;
				e1 -= e2;
				}
			}
#ifndef Inaccurate_Divide
		else if (e >= -Ten_pmax) {
			/* rv = */ rounded_quotient(dval(rv), tens[-e]);
			if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv))
				goto ret;
			e1 -= e;
			}
#endif
		}
 rv_notOK:
	e1 += nd - k;

	/* Get starting approximation = rv * 10**e1 */

	e2 = 0;
	if (e1 > 0) {
		if ( (i = e1 & 15) !=0)
			dval(rv) *= tens[i];
		if (e1 &= ~15) {
			e1 >>= 4;
			while(e1 >= (1 << n_bigtens-1)) {
				e2 += ((word0(rv) & Exp_mask)
					>> Exp_shift1) - Bias;
				word0(rv) &= ~Exp_mask;
				word0(rv) |= Bias << Exp_shift1;
				dval(rv) *= bigtens[n_bigtens-1];
				e1 -= 1 << n_bigtens-1;
				}
			e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
			word0(rv) &= ~Exp_mask;
			word0(rv) |= Bias << Exp_shift1;
			for(j = 0; e1 > 0; j++, e1 >>= 1)
				if (e1 & 1)
					dval(rv) *= bigtens[j];
			}
		}
	else if (e1 < 0) {
		e1 = -e1;
		if ( (i = e1 & 15) !=0)
			dval(rv) /= tens[i];
		if (e1 &= ~15) {
			e1 >>= 4;
			while(e1 >= (1 << n_bigtens-1)) {
				e2 += ((word0(rv) & Exp_mask)
					>> Exp_shift1) - Bias;
				word0(rv) &= ~Exp_mask;
				word0(rv) |= Bias << Exp_shift1;
				dval(rv) *= tinytens[n_bigtens-1];
				e1 -= 1 << n_bigtens-1;
				}
			e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
			word0(rv) &= ~Exp_mask;
			word0(rv) |= Bias << Exp_shift1;
			for(j = 0; e1 > 0; j++, e1 >>= 1)
				if (e1 & 1)
					dval(rv) *= tinytens[j];
			}
		}
#ifdef IBM
	/* e2 is a correction to the (base 2) exponent of the return
	 * value, reflecting adjustments above to avoid overflow in the
	 * native arithmetic.  For native IBM (base 16) arithmetic, we
	 * must multiply e2 by 4 to change from base 16 to 2.
	 */
	e2 <<= 2;
#endif
	rvb = d2b(dval(rv), &rve, &rvbits);	/* rv = rvb * 2^rve */
	rve += e2;
	if ((j = rvbits - nbits) > 0) {
		rshift(rvb, j);
		rvbits = nbits;
		rve += j;
		}
	bb0 = 0;	/* trailing zero bits in rvb */
	e2 = rve + rvbits - nbits;
	if (e2 > fpi->emax + 1)
		goto huge;
	rve1 = rve + rvbits - nbits;
	if (e2 < (emin = fpi->emin)) {
		denorm = 1;
		j = rve - emin;
		if (j > 0) {
			rvb = lshift(rvb, j);
			rvbits += j;
			}
		else if (j < 0) {
			rvbits += j;
			if (rvbits <= 0) {
				if (rvbits < -1) {
 ufl:
					rvb->wds = 0;
					rvb->x[0] = 0;
					*exp = emin;
					irv = STRTOG_Underflow | STRTOG_Inexlo;
					goto ret;
					}
				rvb->x[0] = rvb->wds = rvbits = 1;
				}
			else
				rshift(rvb, -j);
			}
		rve = rve1 = emin;
		if (sudden_underflow && e2 + 1 < emin)
			goto ufl;
		}

	/* Now the hard part -- adjusting rv to the correct value.*/

	/* Put digits into bd: true value = bd * 10^e */

	bd0 = s2b(s0, nd0, nd, y);

	for(;;) {
		bd = Balloc(bd0->k);
		Bcopy(bd, bd0);
		bb = Balloc(rvb->k);
		Bcopy(bb, rvb);
		bbbits = rvbits - bb0;
		bbe = rve + bb0;
		bs = i2b(1);

		if (e >= 0) {
			bb2 = bb5 = 0;
			bd2 = bd5 = e;
			}
		else {
			bb2 = bb5 = -e;
			bd2 = bd5 = 0;
			}
		if (bbe >= 0)
			bb2 += bbe;
		else
			bd2 -= bbe;
		bs2 = bb2;
		j = nbits + 1 - bbbits;
		i = bbe + bbbits - nbits;
		if (i < emin)	/* denormal */
			j += i - emin;
		bb2 += j;
		bd2 += j;
		i = bb2 < bd2 ? bb2 : bd2;
		if (i > bs2)
			i = bs2;
		if (i > 0) {
			bb2 -= i;
			bd2 -= i;
			bs2 -= i;
			}
		if (bb5 > 0) {
			bs = pow5mult(bs, bb5);
			bb1 = mult(bs, bb);
			Bfree(bb);
			bb = bb1;
			}
		bb2 -= bb0;
		if (bb2 > 0)
			bb = lshift(bb, bb2);
		else if (bb2 < 0)
			rshift(bb, -bb2);
		if (bd5 > 0)
			bd = pow5mult(bd, bd5);
		if (bd2 > 0)
			bd = lshift(bd, bd2);
		if (bs2 > 0)
			bs = lshift(bs, bs2);
		asub = 1;
		inex = STRTOG_Inexhi;
		delta = diff(bb, bd);
		if (delta->wds <= 1 && !delta->x[0])
			break;
		dsign = delta->sign;
		delta->sign = finished = 0;
		L = 0;
		i = cmp(delta, bs);
		if (rd && i <= 0) {
			irv = STRTOG_Normal;
			if ( (finished = dsign ^ (rd&1)) !=0) {
				if (dsign != 0) {
					irv |= STRTOG_Inexhi;
					goto adj1;
					}
				irv |= STRTOG_Inexlo;
				if (rve1 == emin)
					goto adj1;
				for(i = 0, j = nbits; j >= ULbits;
						i++, j -= ULbits) {
					if (rvb->x[i] & ALL_ON)
						goto adj1;
					}
				if (j > 1 && lo0bits(rvb->x + i) < j - 1)
					goto adj1;
				rve = rve1 - 1;
				rvb = set_ones(rvb, rvbits = nbits);
				break;
				}
			irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi;
			break;
			}
		if (i < 0) {
			/* Error is less than half an ulp -- check for
			 * special case of mantissa a power of two.
			 */
			irv = dsign
				? STRTOG_Normal | STRTOG_Inexlo
				: STRTOG_Normal | STRTOG_Inexhi;
			if (dsign || bbbits > 1 || denorm || rve1 == emin)
				break;
			delta = lshift(delta,1);
			if (cmp(delta, bs) > 0) {
				irv = STRTOG_Normal | STRTOG_Inexlo;
				goto drop_down;
				}
			break;
			}
		if (i == 0) {
			/* exactly half-way between */
			if (dsign) {
				if (denorm && all_on(rvb, rvbits)) {
					/*boundary case -- increment exponent*/
					rvb->wds = 1;
					rvb->x[0] = 1;
					rve = emin + nbits - (rvbits = 1);
					irv = STRTOG_Normal | STRTOG_Inexhi;
					denorm = 0;
					break;
					}
				irv = STRTOG_Normal | STRTOG_Inexlo;
				}
			else if (bbbits == 1) {
				irv = STRTOG_Normal;
 drop_down:
				/* boundary case -- decrement exponent */
				if (rve1 == emin) {
					irv = STRTOG_Normal | STRTOG_Inexhi;
					if (rvb->wds == 1 && rvb->x[0] == 1)
						sudden_underflow = 1;
					break;
					}
				rve -= nbits;
				rvb = set_ones(rvb, rvbits = nbits);
				break;
				}
			else
				irv = STRTOG_Normal | STRTOG_Inexhi;
			if (bbbits < nbits && !denorm || !(rvb->x[0] & 1))
				break;
			if (dsign) {
				rvb = increment(rvb);
				if ( (j = rvbits & kmask) !=0)
					j = ULbits - j;
				if (hi0bits(rvb->x[(rvb->wds - 1) >> kshift])
						!= j)
					rvbits++;
				irv = STRTOG_Normal | STRTOG_Inexhi;
				}
			else {
				if (bbbits == 1)
					goto undfl;
				decrement(rvb);
				irv = STRTOG_Normal | STRTOG_Inexlo;
				}
			break;
			}
		if ((dval(adj) = ratio(delta, bs)) <= 2.) {
 adj1:
			inex = STRTOG_Inexlo;
			if (dsign) {
				asub = 0;
				inex = STRTOG_Inexhi;
				}
			else if (denorm && bbbits <= 1) {
 undfl:
				rvb->wds = 0;
				rve = emin;
				irv = STRTOG_Underflow | STRTOG_Inexlo;
				break;
				}
			adj0 = dval(adj) = 1.;
			}
		else {
			adj0 = dval(adj) *= 0.5;
			if (dsign) {
				asub = 0;
				inex = STRTOG_Inexlo;
				}
			if (dval(adj) < 2147483647.) {
				L = adj0;
				adj0 -= L;
				switch(rd) {
				  case 0:
					if (adj0 >= .5)
						goto inc_L;
					break;
				  case 1:
					if (asub && adj0 > 0.)
						goto inc_L;
					break;
				  case 2:
					if (!asub && adj0 > 0.) {
 inc_L:
						L++;
						inex = STRTOG_Inexact - inex;
						}
				  }
				dval(adj) = L;
				}
			}
		y = rve + rvbits;

		/* adj *= ulp(dval(rv)); */
		/* if (asub) rv -= adj; else rv += adj; */

		if (!denorm && rvbits < nbits) {
			rvb = lshift(rvb, j = nbits - rvbits);
			rve -= j;
			rvbits = nbits;
			}
		ab = d2b(dval(adj), &abe, &abits);
		if (abe < 0)
			rshift(ab, -abe);
		else if (abe > 0)
			ab = lshift(ab, abe);
		rvb0 = rvb;
		if (asub) {
			/* rv -= adj; */
			j = hi0bits(rvb->x[rvb->wds-1]);
			rvb = diff(rvb, ab);
			k = rvb0->wds - 1;
			if (denorm)
				/* do nothing */;
			else if (rvb->wds <= k
				|| hi0bits( rvb->x[k]) >
				   hi0bits(rvb0->x[k])) {
				/* unlikely; can only have lost 1 high bit */
				if (rve1 == emin) {
					--rvbits;
					denorm = 1;
					}
				else {
					rvb = lshift(rvb, 1);
					--rve;
					--rve1;
					L = finished = 0;
					}
				}
			}
		else {
			rvb = sum(rvb, ab);
			k = rvb->wds - 1;
			if (k >= rvb0->wds
			 || hi0bits(rvb->x[k]) < hi0bits(rvb0->x[k])) {
				if (denorm) {
					if (++rvbits == nbits)
						denorm = 0;
					}
				else {
					rshift(rvb, 1);
					rve++;
					rve1++;
					L = 0;
					}
				}
			}
		Bfree(ab);
		Bfree(rvb0);
		if (finished)
			break;

		z = rve + rvbits;
		if (y == z && L) {
			/* Can we stop now? */
			tol = dval(adj) * 5e-16; /* > max rel error */
			dval(adj) = adj0 - .5;
			if (dval(adj) < -tol) {
				if (adj0 > tol) {
					irv |= inex;
					break;
					}
				}
			else if (dval(adj) > tol && adj0 < 1. - tol) {
				irv |= inex;
				break;
				}
			}
		bb0 = denorm ? 0 : trailz(rvb);
		Bfree(bb);
		Bfree(bd);
		Bfree(bs);
		Bfree(delta);
		}
	if (!denorm && (j = nbits - rvbits)) {
		if (j > 0)
			rvb = lshift(rvb, j);
		else
			rshift(rvb, -j);
		rve -= j;
		}
	*exp = rve;
	Bfree(bb);
	Bfree(bd);
	Bfree(bs);
	Bfree(bd0);
	Bfree(delta);
	if (rve > fpi->emax) {
 huge:
		rvb->wds = 0;
		irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
#ifndef NO_ERRNO
		errno = ERANGE;
#endif
 infnanexp:
		*exp = fpi->emax + 1;
		}
 ret:
	if (denorm) {
		if (sudden_underflow) {
			rvb->wds = 0;
			irv = STRTOG_Underflow | STRTOG_Inexlo;
			}
		else  {
			irv = (irv & ~STRTOG_Retmask) |
				(rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero);
			if (irv & STRTOG_Inexact)
				irv |= STRTOG_Underflow;
			}
		}
	if (se)
		*se = (char *)s;
	if (sign)
		irv |= STRTOG_Neg;
	if (rvb) {
		copybits(bits, nbits, rvb);
		Bfree(rvb);
		}
	return irv;
	}
Commit	Line	Data
9385eb3d A	1	/****************************************************************
	2
	3	The author of this software is David M. Gay.
	4
	5	Copyright (C) 1998-2001 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
3d9156a7 A	29	/* Please send bug reports to David M. Gay (dmg at acm dot org,
3d9156a7 A	30	* with " at " changed at "@" and " dot " changed to "."). */
9385eb3d A	31
	32	#include "gdtoaimp.h"
	33
	34	#ifdef USE_LOCALE
	35	#include "locale.h"
	36	#endif
	37
	38	static CONST int
	39	fivesbits[] = { 0, 3, 5, 7, 10, 12, 14, 17, 19, 21,
	40	24, 26, 28, 31, 33, 35, 38, 40, 42, 45,
	41	47, 49, 52
	42	#ifdef VAX
	43	, 54, 56
	44	#endif
	45	};
	46
	47	Bigint *
	48	#ifdef KR_headers
	49	increment(b) Bigint *b;
	50	#else
	51	increment(Bigint *b)
	52	#endif
	53	{
	54	ULong x, xe;
	55	Bigint *b1;
	56	#ifdef Pack_16
	57	ULong carry = 1, y;
	58	#endif
	59
	60	x = b->x;
	61	xe = x + b->wds;
	62	#ifdef Pack_32
	63	do {
	64	if (*x < (ULong)0xffffffffL) {
	65	++*x;
	66	return b;
	67	}
	68	*x++ = 0;
	69	} while(x < xe);
	70	#else
	71	do {
	72	y = *x + carry;
	73	carry = y >> 16;
	74	*x++ = y & 0xffff;
	75	if (!carry)
	76	return b;
	77	} while(x < xe);
	78	if (carry)
	79	#endif
	80	{
	81	if (b->wds >= b->maxwds) {
	82	b1 = Balloc(b->k+1);
	83	Bcopy(b1,b);
	84	Bfree(b);
	85	b = b1;
	86	}
	87	b->x[b->wds++] = 1;
	88	}
	89	return b;
	90	}
	91
	92	int
	93	#ifdef KR_headers
	94	decrement(b) Bigint *b;
95	#else
96	decrement(Bigint *b)
97	#endif
98	{
99	ULong x, xe;
100	#ifdef Pack_16
101	ULong borrow = 1, y;
102	#endif
103
104	x = b->x;
105	xe = x + b->wds;
106	#ifdef Pack_32
107	do {
108	if (*x) {
109	--*x;
110	break;
111	}
112	*x++ = 0xffffffffL;
113	}
114	while(x < xe);
115	#else
116	do {
117	y = *x - borrow;
118	borrow = (y & 0x10000) >> 16;
119	*x++ = y & 0xffff;
120	} while(borrow && x < xe);
121	#endif
122	return STRTOG_Inexlo;
123	}
124
125	static int
126	#ifdef KR_headers
127	all_on(b, n) Bigint *b; int n;
128	#else
129	all_on(Bigint *b, int n)
130	#endif
131	{
132	ULong x, xe;
133
134	x = b->x;
135	xe = x + (n >> kshift);
136	while(x < xe)
137	if ((*x++ & ALL_ON) != ALL_ON)
138	return 0;
139	if (n &= kmask)
140	return ((*x \| (ALL_ON << n)) & ALL_ON) == ALL_ON;
141	return 1;
142	}
143
144	Bigint *
145	#ifdef KR_headers
146	set_ones(b, n) Bigint *b; int n;
147	#else
148	set_ones(Bigint *b, int n)
149	#endif
150	{
151	int k;
152	ULong x, xe;
153
154	k = (n + ((1 << kshift) - 1)) >> kshift;
155	if (b->k < k) {
156	Bfree(b);
157	b = Balloc(k);
158	}
159	k = n >> kshift;
160	if (n &= kmask)
161	k++;
162	b->wds = k;
163	x = b->x;
164	xe = x + k;
165	while(x < xe)
166	*x++ = ALL_ON;
167	if (n)
168	x[-1] >>= ULbits - n;
169	return b;
170	}
171
172	static int
173	rvOK
174	#ifdef KR_headers
175	(d, fpi, exp, bits, exact, rd, irv)
176	double d; FPI fpi; Long exp; ULong bits; int exact, rd, irv;
177	#else
178	(double d, FPI fpi, Long exp, ULong bits, int exact, int rd, int irv)
179	#endif
180	{
181	Bigint *b;
182	ULong carry, inex, lostbits;
183	int bdif, e, j, k, k1, nb, rv;
184
185	carry = rv = 0;
186	b = d2b(d, &e, &bdif);
187	bdif -= nb = fpi->nbits;
188	e += bdif;
189	if (bdif <= 0) {
190	if (exact)
191	goto trunc;
192	goto ret;
193	}
194	if (P == nb) {
195	if (
196	#ifndef IMPRECISE_INEXACT
197	exact &&
198	#endif
199	fpi->rounding ==
200	#ifdef RND_PRODQUOT
201	FPI_Round_near
202	#else
203	Flt_Rounds
204	#endif
205	) goto trunc;
206	goto ret;
207	}
208	switch(rd) {
209	case 1:
210	goto trunc;
211	case 2:
212	break;
213	default: /* round near */
214	k = bdif - 1;
215	if (k < 0)
216	goto trunc;
217	if (!k) {
218	if (!exact)
219	goto ret;
220	if (b->x[0] & 2)
221	break;
222	goto trunc;
223	}
224	if (b->x[k>>kshift] & ((ULong)1 << (k & kmask)))
225	break;
226	goto trunc;
227	}
228	/* "break" cases: round up 1 bit, then truncate; bdif > 0 */
229	carry = 1;
230	trunc:
231	inex = lostbits = 0;
232	if (bdif > 0) {
233	if ( (lostbits = any_on(b, bdif)) !=0)
234	inex = STRTOG_Inexlo;
235	rshift(b, bdif);
236	if (carry) {
237	inex = STRTOG_Inexhi;
238	b = increment(b);
239	if ( (j = nb & kmask) !=0)
3d9156a7	240	j = ULbits - j;
9385eb3d A	241	if (hi0bits(b->x[b->wds - 1]) != j) {
	242	if (!lostbits)
	243	lostbits = b->x[0] & 1;
	244	rshift(b, 1);
	245	e++;
	246	}
	247	}
	248	}
	249	else if (bdif < 0)
	250	b = lshift(b, -bdif);
	251	if (e < fpi->emin) {
	252	k = fpi->emin - e;
	253	e = fpi->emin;
	254	if (k > nb \|\| fpi->sudden_underflow) {
	255	b->wds = inex = 0;
	256	*irv = STRTOG_Underflow \| STRTOG_Inexlo;
	257	}
	258	else {
	259	k1 = k - 1;
	260	if (k1 > 0 && !lostbits)
	261	lostbits = any_on(b, k1);
	262	if (!lostbits && !exact)
	263	goto ret;
	264	lostbits \|=
	265	carry = b->x[k1>>kshift] & (1 << (k1 & kmask));
	266	rshift(b, k);
	267	*irv = STRTOG_Denormal;
	268	if (carry) {
	269	b = increment(b);
	270	inex = STRTOG_Inexhi \| STRTOG_Underflow;
	271	}
	272	else if (lostbits)
	273	inex = STRTOG_Inexlo \| STRTOG_Underflow;
	274	}
	275	}
	276	else if (e > fpi->emax) {
	277	e = fpi->emax + 1;
	278	*irv = STRTOG_Infinite \| STRTOG_Overflow \| STRTOG_Inexhi;
	279	#ifndef NO_ERRNO
	280	errno = ERANGE;
	281	#endif
	282	b->wds = inex = 0;
	283	}
	284	*exp = e;
	285	copybits(bits, nb, b);
	286	*irv \|= inex;
	287	rv = 1;
	288	ret:
	289	Bfree(b);
	290	return rv;
	291	}
	292
	293	static int
	294	#ifdef KR_headers
	295	mantbits(d) double d;
	296	#else
	297	mantbits(double d)
	298	#endif
	299	{
	300	ULong L;
	301	#ifdef VAX
	302	L = word1(d) << 16 \| word1(d) >> 16;
	303	if (L)
	304	#else
305	if ( (L = word1(d)) !=0)
306	#endif
307	return P - lo0bits(&L);
308	#ifdef VAX
309	L = word0(d) << 16 \| word0(d) >> 16 \| Exp_msk11;
310	#else
311	L = word0(d) \| Exp_msk1;
312	#endif
313	return P - 32 - lo0bits(&L);
314	}
315
316	int
317	strtodg
318	#ifdef KR_headers
319	(s00, se, fpi, exp, bits)
320	CONST char s00; char se; FPI fpi; Long exp; ULong bits;
321	#else
322	(CONST char s00, char se, FPI fpi, Long exp, ULong bits)
323	#endif
324	{
325	int abe, abits, asub;
3d9156a7 A	326	int bb0, bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, denorm;
3d9156a7 A	327	int dsign, e, e1, e2, emin, esign, finished, i, inex, irv;
9385eb3d A	328	int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign;
	329	int sudden_underflow;
	330	CONST char s, s0, *s1;
	331	double adj, adj0, rv, tol;
	332	Long L;
	333	ULong y, z;
	334	Bigint ab, bb, bb1, bd, bd0, bs, delta, rvb, *rvb0;
	335
	336	irv = STRTOG_Zero;
	337	denorm = sign = nz0 = nz = 0;
	338	dval(rv) = 0.;
	339	rvb = 0;
	340	nbits = fpi->nbits;
	341	for(s = s00;;s++) switch(*s) {
	342	case '-':
	343	sign = 1;
	344	/* no break */
	345	case '+':
	346	if (*++s)
	347	goto break2;
	348	/* no break */
	349	case 0:
	350	sign = 0;
	351	irv = STRTOG_NoNumber;
	352	s = s00;
	353	goto ret;
	354	case '\t':
	355	case '\n':
	356	case '\v':
	357	case '\f':
	358	case '\r':
	359	case ' ':
	360	continue;
	361	default:
	362	goto break2;
	363	}
	364	break2:
	365	if (*s == '0') {
	366	#ifndef NO_HEX_FP
	367	switch(s[1]) {
	368	case 'x':
	369	case 'X':
	370	irv = gethex(&s, fpi, exp, &rvb, sign);
	371	if (irv == STRTOG_NoNumber) {
	372	s = s00;
	373	sign = 0;
	374	}
	375	goto ret;
	376	}
	377	#endif
	378	nz0 = 1;
	379	while(*++s == '0') ;
	380	if (!*s)
	381	goto ret;
	382	}
	383	sudden_underflow = fpi->sudden_underflow;
	384	s0 = s;
	385	y = z = 0;
3d9156a7	386	for(decpt = nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
9385eb3d A	387	if (nd < 9)
	388	y = 10*y + c - '0';
	389	else if (nd < 16)
	390	z = 10*z + c - '0';
	391	nd0 = nd;
	392	#ifdef USE_LOCALE
	393	if (c == *localeconv()->decimal_point)
	394	#else
	395	if (c == '.')
	396	#endif
	397	{
3d9156a7	398	decpt = 1;
9385eb3d A	399	c = *++s;
	400	if (!nd) {
	401	for(; c == '0'; c = *++s)
	402	nz++;
	403	if (c > '0' && c <= '9') {
	404	s0 = s;
	405	nf += nz;
	406	nz = 0;
	407	goto have_dig;
	408	}
	409	goto dig_done;
	410	}
	411	for(; c >= '0' && c <= '9'; c = *++s) {
	412	have_dig:
	413	nz++;
	414	if (c -= '0') {
	415	nf += nz;
	416	for(i = 1; i < nz; i++)
	417	if (nd++ < 9)
	418	y *= 10;
	419	else if (nd <= DBL_DIG + 1)
	420	z *= 10;
	421	if (nd++ < 9)
	422	y = 10*y + c;
	423	else if (nd <= DBL_DIG + 1)
	424	z = 10*z + c;
	425	nz = 0;
	426	}
	427	}
	428	}
	429	dig_done:
	430	e = 0;
	431	if (c == 'e' \|\| c == 'E') {
	432	if (!nd && !nz && !nz0) {
	433	irv = STRTOG_NoNumber;
	434	s = s00;
	435	goto ret;
	436	}
	437	s00 = s;
	438	esign = 0;
	439	switch(c = *++s) {
	440	case '-':
	441	esign = 1;
	442	case '+':
	443	c = *++s;
	444	}
	445	if (c >= '0' && c <= '9') {
	446	while(c == '0')
	447	c = *++s;
	448	if (c > '0' && c <= '9') {
	449	L = c - '0';
	450	s1 = s;
	451	while((c = *++s) >= '0' && c <= '9')
	452	L = 10*L + c - '0';
	453	if (s - s1 > 8 \|\| L > 19999)
	454	/* Avoid confusion from exponents
	455	* so large that e might overflow.
	456	*/
	457	e = 19999; /* safe for 16 bit ints */
	458	else
	459	e = (int)L;
	460	if (esign)
	461	e = -e;
	462	}
463	else
464	e = 0;
465	}
466	else
467	s = s00;
468	}
469	if (!nd) {
470	if (!nz && !nz0) {
471	#ifdef INFNAN_CHECK
472	/* Check for Nan and Infinity */
3d9156a7 A	473	if (!decpt)
3d9156a7 A	474	switch(c) {
9385eb3d A	475	case 'i':
	476	case 'I':
	477	if (match(&s,"nf")) {
	478	--s;
	479	if (!match(&s,"inity"))
	480	++s;
	481	irv = STRTOG_Infinite;
	482	goto infnanexp;
	483	}
	484	break;
	485	case 'n':
	486	case 'N':
	487	if (match(&s, "an")) {
	488	irv = STRTOG_NaN;
	489	*exp = fpi->emax + 1;
	490	#ifndef No_Hex_NaN
	491	if (s == '(') /)*/
	492	irv = hexnan(&s, fpi, bits);
	493	#endif
	494	goto infnanexp;
	495	}
	496	}
	497	#endif /* INFNAN_CHECK */
	498	irv = STRTOG_NoNumber;
	499	s = s00;
	500	}
	501	goto ret;
	502	}
	503
	504	irv = STRTOG_Normal;
	505	e1 = e -= nf;
	506	rd = 0;
	507	switch(fpi->rounding & 3) {
	508	case FPI_Round_up:
	509	rd = 2 - sign;
	510	break;
	511	case FPI_Round_zero:
	512	rd = 1;
	513	break;
	514	case FPI_Round_down:
	515	rd = 1 + sign;
	516	}
	517
	518	/* Now we have nd0 digits, starting at s0, followed by a
	519	* decimal point, followed by nd-nd0 digits. The number we're
	520	* after is the integer represented by those digits times
	521	* 10*e /
	522
	523	if (!nd0)
	524	nd0 = nd;
	525	k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
	526	dval(rv) = y;
	527	if (k > 9)
	528	dval(rv) = tens[k - 9] * dval(rv) + z;
	529	bd0 = 0;
	530	if (nbits <= P && nd <= DBL_DIG) {
	531	if (!e) {
	532	if (rvOK(dval(rv), fpi, exp, bits, 1, rd, &irv))
	533	goto ret;
	534	}
	535	else if (e > 0) {
	536	if (e <= Ten_pmax) {
	537	#ifdef VAX
	538	goto vax_ovfl_check;
539	#else
540	i = fivesbits[e] + mantbits(dval(rv)) <= P;
541	/* rv = */ rounded_product(dval(rv), tens[e]);
542	if (rvOK(dval(rv), fpi, exp, bits, i, rd, &irv))
543	goto ret;
544	e1 -= e;
545	goto rv_notOK;
546	#endif
547	}
548	i = DBL_DIG - nd;
549	if (e <= Ten_pmax + i) {
550	/* A fancier test would sometimes let us do
551	* this for larger i values.
552	*/
553	e2 = e - i;
554	e1 -= i;
555	dval(rv) *= tens[i];
556	#ifdef VAX
557	/* VAX exponent range is so narrow we must
558	* worry about overflow here...
559	*/
560	vax_ovfl_check:
561	dval(adj) = dval(rv);
562	word0(adj) -= P*Exp_msk1;
563	/* adj = */ rounded_product(dval(adj), tens[e2]);
564	if ((word0(adj) & Exp_mask)
565	> Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
566	goto rv_notOK;
567	word0(adj) += P*Exp_msk1;
568	dval(rv) = dval(adj);
569	#else
570	/* rv = */ rounded_product(dval(rv), tens[e2]);
571	#endif
572	if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv))
573	goto ret;
574	e1 -= e2;
575	}
576	}
577	#ifndef Inaccurate_Divide
578	else if (e >= -Ten_pmax) {
579	/* rv = */ rounded_quotient(dval(rv), tens[-e]);
580	if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv))
581	goto ret;
582	e1 -= e;
583	}
584	#endif
585	}
586	rv_notOK:
587	e1 += nd - k;
588
589	/* Get starting approximation = rv * 10*e1 /
590
591	e2 = 0;
592	if (e1 > 0) {
593	if ( (i = e1 & 15) !=0)
594	dval(rv) *= tens[i];
595	if (e1 &= ~15) {
596	e1 >>= 4;
597	while(e1 >= (1 << n_bigtens-1)) {
598	e2 += ((word0(rv) & Exp_mask)
599	>> Exp_shift1) - Bias;
600	word0(rv) &= ~Exp_mask;
601	word0(rv) \|= Bias << Exp_shift1;
602	dval(rv) *= bigtens[n_bigtens-1];
603	e1 -= 1 << n_bigtens-1;
604	}
605	e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
606	word0(rv) &= ~Exp_mask;
607	word0(rv) \|= Bias << Exp_shift1;
608	for(j = 0; e1 > 0; j++, e1 >>= 1)
609	if (e1 & 1)
610	dval(rv) *= bigtens[j];
611	}
612	}
613	else if (e1 < 0) {
614	e1 = -e1;
615	if ( (i = e1 & 15) !=0)
616	dval(rv) /= tens[i];
617	if (e1 &= ~15) {
618	e1 >>= 4;
619	while(e1 >= (1 << n_bigtens-1)) {
620	e2 += ((word0(rv) & Exp_mask)
621	>> Exp_shift1) - Bias;
622	word0(rv) &= ~Exp_mask;
623	word0(rv) \|= Bias << Exp_shift1;
624	dval(rv) *= tinytens[n_bigtens-1];
625	e1 -= 1 << n_bigtens-1;
626	}
627	e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
628	word0(rv) &= ~Exp_mask;
629	word0(rv) \|= Bias << Exp_shift1;
630	for(j = 0; e1 > 0; j++, e1 >>= 1)
631	if (e1 & 1)
632	dval(rv) *= tinytens[j];
633	}
634	}
3d9156a7 A	635	#ifdef IBM
	636	/* e2 is a correction to the (base 2) exponent of the return
	637	* value, reflecting adjustments above to avoid overflow in the
	638	* native arithmetic. For native IBM (base 16) arithmetic, we
	639	* must multiply e2 by 4 to change from base 16 to 2.
	640	*/
	641	e2 <<= 2;
	642	#endif
9385eb3d A	643	rvb = d2b(dval(rv), &rve, &rvbits); /* rv = rvb * 2^rve */
	644	rve += e2;
	645	if ((j = rvbits - nbits) > 0) {
	646	rshift(rvb, j);
	647	rvbits = nbits;
	648	rve += j;
	649	}
	650	bb0 = 0; /* trailing zero bits in rvb */
	651	e2 = rve + rvbits - nbits;
eb1cde05 A	652	if (e2 > fpi->emax + 1)
eb1cde05 A	653	goto huge;
9385eb3d A	654	rve1 = rve + rvbits - nbits;
	655	if (e2 < (emin = fpi->emin)) {
	656	denorm = 1;
	657	j = rve - emin;
	658	if (j > 0) {
	659	rvb = lshift(rvb, j);
	660	rvbits += j;
	661	}
	662	else if (j < 0) {
	663	rvbits += j;
	664	if (rvbits <= 0) {
	665	if (rvbits < -1) {
	666	ufl:
	667	rvb->wds = 0;
	668	rvb->x[0] = 0;
	669	*exp = emin;
	670	irv = STRTOG_Underflow \| STRTOG_Inexlo;
	671	goto ret;
	672	}
	673	rvb->x[0] = rvb->wds = rvbits = 1;
	674	}
	675	else
	676	rshift(rvb, -j);
	677	}
	678	rve = rve1 = emin;
	679	if (sudden_underflow && e2 + 1 < emin)
	680	goto ufl;
	681	}
	682
	683	/* Now the hard part -- adjusting rv to the correct value.*/
	684
	685	/* Put digits into bd: true value = bd * 10^e */
	686
	687	bd0 = s2b(s0, nd0, nd, y);
	688
	689	for(;;) {
	690	bd = Balloc(bd0->k);
	691	Bcopy(bd, bd0);
	692	bb = Balloc(rvb->k);
	693	Bcopy(bb, rvb);
	694	bbbits = rvbits - bb0;
	695	bbe = rve + bb0;
	696	bs = i2b(1);
	697
	698	if (e >= 0) {
	699	bb2 = bb5 = 0;
	700	bd2 = bd5 = e;
	701	}
	702	else {
	703	bb2 = bb5 = -e;
	704	bd2 = bd5 = 0;
	705	}
	706	if (bbe >= 0)
	707	bb2 += bbe;
	708	else
	709	bd2 -= bbe;
	710	bs2 = bb2;
	711	j = nbits + 1 - bbbits;
	712	i = bbe + bbbits - nbits;
	713	if (i < emin) /* denormal */
	714	j += i - emin;
	715	bb2 += j;
	716	bd2 += j;
	717	i = bb2 < bd2 ? bb2 : bd2;
718	if (i > bs2)
719	i = bs2;
720	if (i > 0) {
721	bb2 -= i;
722	bd2 -= i;
723	bs2 -= i;
724	}
725	if (bb5 > 0) {
726	bs = pow5mult(bs, bb5);
727	bb1 = mult(bs, bb);
728	Bfree(bb);
729	bb = bb1;
730	}
731	bb2 -= bb0;
732	if (bb2 > 0)
733	bb = lshift(bb, bb2);
734	else if (bb2 < 0)
735	rshift(bb, -bb2);
736	if (bd5 > 0)
737	bd = pow5mult(bd, bd5);
738	if (bd2 > 0)
739	bd = lshift(bd, bd2);
740	if (bs2 > 0)
741	bs = lshift(bs, bs2);
742	asub = 1;
743	inex = STRTOG_Inexhi;
744	delta = diff(bb, bd);
745	if (delta->wds <= 1 && !delta->x[0])
746	break;
747	dsign = delta->sign;
748	delta->sign = finished = 0;
749	L = 0;
750	i = cmp(delta, bs);
751	if (rd && i <= 0) {
752	irv = STRTOG_Normal;
753	if ( (finished = dsign ^ (rd&1)) !=0) {
754	if (dsign != 0) {
755	irv \|= STRTOG_Inexhi;
756	goto adj1;
757	}
758	irv \|= STRTOG_Inexlo;
759	if (rve1 == emin)
760	goto adj1;
761	for(i = 0, j = nbits; j >= ULbits;
762	i++, j -= ULbits) {
763	if (rvb->x[i] & ALL_ON)
764	goto adj1;
765	}
766	if (j > 1 && lo0bits(rvb->x + i) < j - 1)
767	goto adj1;
768	rve = rve1 - 1;
769	rvb = set_ones(rvb, rvbits = nbits);
770	break;
771	}
772	irv \|= dsign ? STRTOG_Inexlo : STRTOG_Inexhi;
773	break;
774	}
775	if (i < 0) {
776	/* Error is less than half an ulp -- check for
777	* special case of mantissa a power of two.
778	*/
779	irv = dsign
780	? STRTOG_Normal \| STRTOG_Inexlo
781	: STRTOG_Normal \| STRTOG_Inexhi;
782	if (dsign \|\| bbbits > 1 \|\| denorm \|\| rve1 == emin)
783	break;
784	delta = lshift(delta,1);
785	if (cmp(delta, bs) > 0) {
786	irv = STRTOG_Normal \| STRTOG_Inexlo;
787	goto drop_down;
788	}
789	break;
790	}
791	if (i == 0) {
792	/* exactly half-way between */
793	if (dsign) {
794	if (denorm && all_on(rvb, rvbits)) {
795	/boundary case -- increment exponent/
796	rvb->wds = 1;
797	rvb->x[0] = 1;
798	rve = emin + nbits - (rvbits = 1);
799	irv = STRTOG_Normal \| STRTOG_Inexhi;
800	denorm = 0;
801	break;
802	}
803	irv = STRTOG_Normal \| STRTOG_Inexlo;
804	}
805	else if (bbbits == 1) {
806	irv = STRTOG_Normal;
807	drop_down:
808	/* boundary case -- decrement exponent */
809	if (rve1 == emin) {
810	irv = STRTOG_Normal \| STRTOG_Inexhi;
811	if (rvb->wds == 1 && rvb->x[0] == 1)
812	sudden_underflow = 1;
813	break;
814	}
815	rve -= nbits;
816	rvb = set_ones(rvb, rvbits = nbits);
817	break;
818	}
819	else
820	irv = STRTOG_Normal \| STRTOG_Inexhi;
821	if (bbbits < nbits && !denorm \|\| !(rvb->x[0] & 1))
822	break;
823	if (dsign) {
824	rvb = increment(rvb);
3d9156a7 A	825	if ( (j = rvbits & kmask) !=0)
3d9156a7 A	826	j = ULbits - j;
9385eb3d A	827	if (hi0bits(rvb->x[(rvb->wds - 1) >> kshift])
	828	!= j)
	829	rvbits++;
	830	irv = STRTOG_Normal \| STRTOG_Inexhi;
	831	}
	832	else {
	833	if (bbbits == 1)
	834	goto undfl;
	835	decrement(rvb);
	836	irv = STRTOG_Normal \| STRTOG_Inexlo;
	837	}
	838	break;
	839	}
	840	if ((dval(adj) = ratio(delta, bs)) <= 2.) {
	841	adj1:
	842	inex = STRTOG_Inexlo;
	843	if (dsign) {
	844	asub = 0;
	845	inex = STRTOG_Inexhi;
	846	}
	847	else if (denorm && bbbits <= 1) {
	848	undfl:
	849	rvb->wds = 0;
	850	rve = emin;
	851	irv = STRTOG_Underflow \| STRTOG_Inexlo;
	852	break;
	853	}
	854	adj0 = dval(adj) = 1.;
	855	}
	856	else {
	857	adj0 = dval(adj) *= 0.5;
	858	if (dsign) {
	859	asub = 0;
	860	inex = STRTOG_Inexlo;
	861	}
	862	if (dval(adj) < 2147483647.) {
	863	L = adj0;
	864	adj0 -= L;
	865	switch(rd) {
	866	case 0:
	867	if (adj0 >= .5)
	868	goto inc_L;
	869	break;
	870	case 1:
	871	if (asub && adj0 > 0.)
	872	goto inc_L;
	873	break;
	874	case 2:
	875	if (!asub && adj0 > 0.) {
	876	inc_L:
	877	L++;
	878	inex = STRTOG_Inexact - inex;
	879	}
	880	}
	881	dval(adj) = L;
	882	}
	883	}
	884	y = rve + rvbits;
	885
	886	/* adj = ulp(dval(rv)); /
	887	/* if (asub) rv -= adj; else rv += adj; */
	888
	889	if (!denorm && rvbits < nbits) {
	890	rvb = lshift(rvb, j = nbits - rvbits);
891	rve -= j;
892	rvbits = nbits;
893	}
894	ab = d2b(dval(adj), &abe, &abits);
895	if (abe < 0)
896	rshift(ab, -abe);
897	else if (abe > 0)
898	ab = lshift(ab, abe);
899	rvb0 = rvb;
900	if (asub) {
901	/* rv -= adj; */
902	j = hi0bits(rvb->x[rvb->wds-1]);
903	rvb = diff(rvb, ab);
904	k = rvb0->wds - 1;
905	if (denorm)
906	/* do nothing */;
907	else if (rvb->wds <= k
908	\|\| hi0bits( rvb->x[k]) >
909	hi0bits(rvb0->x[k])) {
910	/* unlikely; can only have lost 1 high bit */
911	if (rve1 == emin) {
912	--rvbits;
913	denorm = 1;
914	}
915	else {
916	rvb = lshift(rvb, 1);
917	--rve;
918	--rve1;
919	L = finished = 0;
920	}
921	}
922	}
923	else {
924	rvb = sum(rvb, ab);
925	k = rvb->wds - 1;
926	if (k >= rvb0->wds
927	\|\| hi0bits(rvb->x[k]) < hi0bits(rvb0->x[k])) {
928	if (denorm) {
929	if (++rvbits == nbits)
930	denorm = 0;
931	}
932	else {
933	rshift(rvb, 1);
934	rve++;
935	rve1++;
936	L = 0;
937	}
938	}
939	}
940	Bfree(ab);
941	Bfree(rvb0);
942	if (finished)
943	break;
944
945	z = rve + rvbits;
946	if (y == z && L) {
947	/* Can we stop now? */
948	tol = dval(adj) * 5e-16; /* > max rel error */
949	dval(adj) = adj0 - .5;
950	if (dval(adj) < -tol) {
951	if (adj0 > tol) {
952	irv \|= inex;
953	break;
954	}
955	}
956	else if (dval(adj) > tol && adj0 < 1. - tol) {
957	irv \|= inex;
958	break;
959	}
960	}
961	bb0 = denorm ? 0 : trailz(rvb);
962	Bfree(bb);
963	Bfree(bd);
964	Bfree(bs);
965	Bfree(delta);
966	}
3d9156a7 A	967	if (!denorm && (j = nbits - rvbits)) {
	968	if (j > 0)
	969	rvb = lshift(rvb, j);
	970	else
	971	rshift(rvb, -j);
9385eb3d A	972	rve -= j;
	973	}
	974	*exp = rve;
	975	Bfree(bb);
	976	Bfree(bd);
	977	Bfree(bs);
	978	Bfree(bd0);
	979	Bfree(delta);
eb1cde05 A	980	if (rve > fpi->emax) {
	981	huge:
	982	rvb->wds = 0;
	983	irv = STRTOG_Infinite \| STRTOG_Overflow \| STRTOG_Inexhi;
	984	#ifndef NO_ERRNO
	985	errno = ERANGE;
	986	#endif
	987	infnanexp:
	988	*exp = fpi->emax + 1;
	989	}
9385eb3d A	990	ret:
	991	if (denorm) {
	992	if (sudden_underflow) {
	993	rvb->wds = 0;
	994	irv = STRTOG_Underflow \| STRTOG_Inexlo;
	995	}
	996	else {
	997	irv = (irv & ~STRTOG_Retmask) \|
	998	(rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero);
	999	if (irv & STRTOG_Inexact)
	1000	irv \|= STRTOG_Underflow;
	1001	}
	1002	}
	1003	if (se)
	1004	se = (char )s;
	1005	if (sign)
	1006	irv \|= STRTOG_Neg;
	1007	if (rvb) {
	1008	copybits(bits, nbits, rvb);
	1009	Bfree(rvb);
	1010	}
	1011	return irv;
	1012	}