[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) {
		rvb->wds = 0;
		irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
#ifndef NO_ERRNO
		errno = ERANGE;
#endif
 infnanexp:
		*exp = fpi->emax + 1;
		goto ret;
		}
	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);
 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;
f74c7596 A	652	if (e2 > fpi->emax) {
	653	rvb->wds = 0;
	654	irv = STRTOG_Infinite \| STRTOG_Overflow \| STRTOG_Inexhi;
	655	#ifndef NO_ERRNO
	656	errno = ERANGE;
	657	#endif
	658	infnanexp:
	659	*exp = fpi->emax + 1;
	660	goto ret;
	661	}
9385eb3d A	662	rve1 = rve + rvbits - nbits;
	663	if (e2 < (emin = fpi->emin)) {
	664	denorm = 1;
	665	j = rve - emin;
	666	if (j > 0) {
	667	rvb = lshift(rvb, j);
	668	rvbits += j;
	669	}
	670	else if (j < 0) {
	671	rvbits += j;
	672	if (rvbits <= 0) {
	673	if (rvbits < -1) {
	674	ufl:
	675	rvb->wds = 0;
	676	rvb->x[0] = 0;
	677	*exp = emin;
	678	irv = STRTOG_Underflow \| STRTOG_Inexlo;
	679	goto ret;
	680	}
	681	rvb->x[0] = rvb->wds = rvbits = 1;
	682	}
	683	else
	684	rshift(rvb, -j);
	685	}
	686	rve = rve1 = emin;
	687	if (sudden_underflow && e2 + 1 < emin)
	688	goto ufl;
	689	}
	690
	691	/* Now the hard part -- adjusting rv to the correct value.*/
	692
	693	/* Put digits into bd: true value = bd * 10^e */
	694
	695	bd0 = s2b(s0, nd0, nd, y);
	696
	697	for(;;) {
	698	bd = Balloc(bd0->k);
	699	Bcopy(bd, bd0);
	700	bb = Balloc(rvb->k);
	701	Bcopy(bb, rvb);
	702	bbbits = rvbits - bb0;
	703	bbe = rve + bb0;
	704	bs = i2b(1);
	705
	706	if (e >= 0) {
	707	bb2 = bb5 = 0;
	708	bd2 = bd5 = e;
	709	}
	710	else {
	711	bb2 = bb5 = -e;
	712	bd2 = bd5 = 0;
	713	}
	714	if (bbe >= 0)
	715	bb2 += bbe;
	716	else
	717	bd2 -= bbe;
	718	bs2 = bb2;
	719	j = nbits + 1 - bbbits;
	720	i = bbe + bbbits - nbits;
	721	if (i < emin) /* denormal */
	722	j += i - emin;
	723	bb2 += j;
	724	bd2 += j;
	725	i = bb2 < bd2 ? bb2 : bd2;
726	if (i > bs2)
727	i = bs2;
728	if (i > 0) {
729	bb2 -= i;
730	bd2 -= i;
731	bs2 -= i;
732	}
733	if (bb5 > 0) {
734	bs = pow5mult(bs, bb5);
735	bb1 = mult(bs, bb);
736	Bfree(bb);
737	bb = bb1;
738	}
739	bb2 -= bb0;
740	if (bb2 > 0)
741	bb = lshift(bb, bb2);
742	else if (bb2 < 0)
743	rshift(bb, -bb2);
744	if (bd5 > 0)
745	bd = pow5mult(bd, bd5);
746	if (bd2 > 0)
747	bd = lshift(bd, bd2);
748	if (bs2 > 0)
749	bs = lshift(bs, bs2);
750	asub = 1;
751	inex = STRTOG_Inexhi;
752	delta = diff(bb, bd);
753	if (delta->wds <= 1 && !delta->x[0])
754	break;
755	dsign = delta->sign;
756	delta->sign = finished = 0;
757	L = 0;
758	i = cmp(delta, bs);
759	if (rd && i <= 0) {
760	irv = STRTOG_Normal;
761	if ( (finished = dsign ^ (rd&1)) !=0) {
762	if (dsign != 0) {
763	irv \|= STRTOG_Inexhi;
764	goto adj1;
765	}
766	irv \|= STRTOG_Inexlo;
767	if (rve1 == emin)
768	goto adj1;
769	for(i = 0, j = nbits; j >= ULbits;
770	i++, j -= ULbits) {
771	if (rvb->x[i] & ALL_ON)
772	goto adj1;
773	}
774	if (j > 1 && lo0bits(rvb->x + i) < j - 1)
775	goto adj1;
776	rve = rve1 - 1;
777	rvb = set_ones(rvb, rvbits = nbits);
778	break;
779	}
780	irv \|= dsign ? STRTOG_Inexlo : STRTOG_Inexhi;
781	break;
782	}
783	if (i < 0) {
784	/* Error is less than half an ulp -- check for
785	* special case of mantissa a power of two.
786	*/
787	irv = dsign
788	? STRTOG_Normal \| STRTOG_Inexlo
789	: STRTOG_Normal \| STRTOG_Inexhi;
790	if (dsign \|\| bbbits > 1 \|\| denorm \|\| rve1 == emin)
791	break;
792	delta = lshift(delta,1);
793	if (cmp(delta, bs) > 0) {
794	irv = STRTOG_Normal \| STRTOG_Inexlo;
795	goto drop_down;
796	}
797	break;
798	}
799	if (i == 0) {
800	/* exactly half-way between */
801	if (dsign) {
802	if (denorm && all_on(rvb, rvbits)) {
803	/boundary case -- increment exponent/
804	rvb->wds = 1;
805	rvb->x[0] = 1;
806	rve = emin + nbits - (rvbits = 1);
807	irv = STRTOG_Normal \| STRTOG_Inexhi;
808	denorm = 0;
809	break;
810	}
811	irv = STRTOG_Normal \| STRTOG_Inexlo;
812	}
813	else if (bbbits == 1) {
814	irv = STRTOG_Normal;
815	drop_down:
816	/* boundary case -- decrement exponent */
817	if (rve1 == emin) {
818	irv = STRTOG_Normal \| STRTOG_Inexhi;
819	if (rvb->wds == 1 && rvb->x[0] == 1)
820	sudden_underflow = 1;
821	break;
822	}
823	rve -= nbits;
824	rvb = set_ones(rvb, rvbits = nbits);
825	break;
826	}
827	else
828	irv = STRTOG_Normal \| STRTOG_Inexhi;
829	if (bbbits < nbits && !denorm \|\| !(rvb->x[0] & 1))
830	break;
831	if (dsign) {
832	rvb = increment(rvb);
3d9156a7 A	833	if ( (j = rvbits & kmask) !=0)
3d9156a7 A	834	j = ULbits - j;
9385eb3d A	835	if (hi0bits(rvb->x[(rvb->wds - 1) >> kshift])
	836	!= j)
	837	rvbits++;
	838	irv = STRTOG_Normal \| STRTOG_Inexhi;
	839	}
	840	else {
	841	if (bbbits == 1)
	842	goto undfl;
	843	decrement(rvb);
	844	irv = STRTOG_Normal \| STRTOG_Inexlo;
	845	}
	846	break;
	847	}
	848	if ((dval(adj) = ratio(delta, bs)) <= 2.) {
	849	adj1:
	850	inex = STRTOG_Inexlo;
	851	if (dsign) {
	852	asub = 0;
	853	inex = STRTOG_Inexhi;
	854	}
	855	else if (denorm && bbbits <= 1) {
	856	undfl:
	857	rvb->wds = 0;
	858	rve = emin;
	859	irv = STRTOG_Underflow \| STRTOG_Inexlo;
	860	break;
	861	}
	862	adj0 = dval(adj) = 1.;
	863	}
	864	else {
	865	adj0 = dval(adj) *= 0.5;
	866	if (dsign) {
	867	asub = 0;
	868	inex = STRTOG_Inexlo;
	869	}
	870	if (dval(adj) < 2147483647.) {
	871	L = adj0;
	872	adj0 -= L;
	873	switch(rd) {
	874	case 0:
	875	if (adj0 >= .5)
	876	goto inc_L;
	877	break;
	878	case 1:
	879	if (asub && adj0 > 0.)
	880	goto inc_L;
	881	break;
	882	case 2:
	883	if (!asub && adj0 > 0.) {
	884	inc_L:
	885	L++;
	886	inex = STRTOG_Inexact - inex;
	887	}
	888	}
	889	dval(adj) = L;
	890	}
	891	}
	892	y = rve + rvbits;
	893
	894	/* adj = ulp(dval(rv)); /
	895	/* if (asub) rv -= adj; else rv += adj; */
	896
	897	if (!denorm && rvbits < nbits) {
	898	rvb = lshift(rvb, j = nbits - rvbits);
899	rve -= j;
900	rvbits = nbits;
901	}
902	ab = d2b(dval(adj), &abe, &abits);
903	if (abe < 0)
904	rshift(ab, -abe);
905	else if (abe > 0)
906	ab = lshift(ab, abe);
907	rvb0 = rvb;
908	if (asub) {
909	/* rv -= adj; */
910	j = hi0bits(rvb->x[rvb->wds-1]);
911	rvb = diff(rvb, ab);
912	k = rvb0->wds - 1;
913	if (denorm)
914	/* do nothing */;
915	else if (rvb->wds <= k
916	\|\| hi0bits( rvb->x[k]) >
917	hi0bits(rvb0->x[k])) {
918	/* unlikely; can only have lost 1 high bit */
919	if (rve1 == emin) {
920	--rvbits;
921	denorm = 1;
922	}
923	else {
924	rvb = lshift(rvb, 1);
925	--rve;
926	--rve1;
927	L = finished = 0;
928	}
929	}
930	}
931	else {
932	rvb = sum(rvb, ab);
933	k = rvb->wds - 1;
934	if (k >= rvb0->wds
935	\|\| hi0bits(rvb->x[k]) < hi0bits(rvb0->x[k])) {
936	if (denorm) {
937	if (++rvbits == nbits)
938	denorm = 0;
939	}
940	else {
941	rshift(rvb, 1);
942	rve++;
943	rve1++;
944	L = 0;
945	}
946	}
947	}
948	Bfree(ab);
949	Bfree(rvb0);
950	if (finished)
951	break;
952
953	z = rve + rvbits;
954	if (y == z && L) {
955	/* Can we stop now? */
956	tol = dval(adj) * 5e-16; /* > max rel error */
957	dval(adj) = adj0 - .5;
958	if (dval(adj) < -tol) {
959	if (adj0 > tol) {
960	irv \|= inex;
961	break;
962	}
963	}
964	else if (dval(adj) > tol && adj0 < 1. - tol) {
965	irv \|= inex;
966	break;
967	}
968	}
969	bb0 = denorm ? 0 : trailz(rvb);
970	Bfree(bb);
971	Bfree(bd);
972	Bfree(bs);
973	Bfree(delta);
974	}
3d9156a7 A	975	if (!denorm && (j = nbits - rvbits)) {
	976	if (j > 0)
	977	rvb = lshift(rvb, j);
	978	else
	979	rshift(rvb, -j);
9385eb3d A	980	rve -= j;
	981	}
	982	*exp = rve;
	983	Bfree(bb);
	984	Bfree(bd);
	985	Bfree(bs);
	986	Bfree(bd0);
	987	Bfree(delta);
	988	ret:
	989	if (denorm) {
	990	if (sudden_underflow) {
	991	rvb->wds = 0;
	992	irv = STRTOG_Underflow \| STRTOG_Inexlo;
	993	}
	994	else {
	995	irv = (irv & ~STRTOG_Retmask) \|
	996	(rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero);
	997	if (irv & STRTOG_Inexact)
	998	irv \|= STRTOG_Underflow;
	999	}
	1000	}
	1001	if (se)
	1002	se = (char )s;
	1003	if (sign)
	1004	irv \|= STRTOG_Neg;
	1005	if (rvb) {
	1006	copybits(bits, nbits, rvb);
	1007	Bfree(rvb);
	1008	}
	1009	return irv;
	1010	}