Subversion Repositories planix.SVN

#include "fconv.h"

/* strtod for IEEE-, VAX-, and IBM-arithmetic machines (dmg).

 *

 * This strtod returns a nearest machine number to the input decimal

 * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are

 * broken by the IEEE round-even rule.  Otherwise ties are broken by

 * biased rounding (add half and chop).

 *

 * Inspired loosely by William D. Clinger's paper "How to Read Floating

 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].

 *

 * Modifications:

 *

 *	1. We only require IEEE, IBM, or VAX double-precision

 *		arithmetic (not IEEE double-extended).

 *	2. We get by with floating-point arithmetic in a case that

 *		Clinger missed -- when we're computing d * 10^n

 *		for a small integer d and the integer n is not too

 *		much larger than 22 (the maximum integer k for which

 *		we can represent 10^k exactly), we may be able to

 *		compute (d*10^k) * 10^(e-k) with just one roundoff.

 *	3. Rather than a bit-at-a-time adjustment of the binary

 *		result in the hard case, we use floating-point

 *		arithmetic to determine the adjustment to within

 *		one bit; only in really hard cases do we need to

 *		compute a second residual.

 *	4. Because of 3., we don't need a large table of powers of 10

 *		for ten-to-e (just some small tables, e.g. of 10^k

 *		for 0 <= k <= 22).

 */

#ifdef RND_PRODQUOT

#define rounded_product(a,b) a = rnd_prod(a, b)

#define rounded_quotient(a,b) a = rnd_quot(a, b)

extern double rnd_prod(double, double), rnd_quot(double, double);

#else

#define rounded_product(a,b) a *= b

#define rounded_quotient(a,b) a /= b

#endif

 static double

ulp(double xarg)

{

	register long L;

	Dul a;

	Dul x;

	x.d = xarg;

	L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;

#ifndef Sudden_Underflow

	if (L > 0) {

#endif

#ifdef IBM

		L |= Exp_msk1 >> 4;

#endif

		word0(a) = L;

		word1(a) = 0;

#ifndef Sudden_Underflow

		}

	else {

		L = -L >> Exp_shift;

		if (L < Exp_shift) {

			word0(a) = 0x80000 >> L;

			word1(a) = 0;

			}

		else {

			word0(a) = 0;

			L -= Exp_shift;

			word1(a) = L >= 31 ? 1 : 1 << 31 - L;

			}

		}

#endif

	return a.d;

	}

 static Bigint *

s2b(CONST char *s, int nd0, int nd, unsigned long y9)

{

	Bigint *b;

	int i, k;

	long x, y;

	x = (nd + 8) / 9;

	for(k = 0, y = 1; x > y; y <<= 1, k++) ;

#ifdef Pack_32

	b = Balloc(k);

	b->x[0] = y9;

	b->wds = 1;

#else

	b = Balloc(k+1);

	b->x[0] = y9 & 0xffff;

	b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;

#endif

	i = 9;

	if (9 < nd0) {

		s += 9;

		do b = multadd(b, 10, *s++ - '0');

			while(++i < nd0);

		s++;

		}

	else

		s += 10;

	for(; i < nd; i++)

		b = multadd(b, 10, *s++ - '0');

	return b;

	}

 static double

b2d(Bigint *a, int *e)

{

	unsigned long *xa, *xa0, w, y, z;

	int k;

	Dul d;

#ifdef VAX

	unsigned long d0, d1;

#else

#define d0 word0(d)

#define d1 word1(d)

#endif

	xa0 = a->x;

	xa = xa0 + a->wds;

	y = *--xa;

#ifdef DEBUG

	if (!y) Bug("zero y in b2d");

#endif

	k = hi0bits(y);

	*e = 32 - k;

#ifdef Pack_32

	if (k < Ebits) {

		d0 = Exp_1 | y >> Ebits - k;

		w = xa > xa0 ? *--xa : 0;

		d1 = y << (32-Ebits) + k | w >> Ebits - k;

		goto ret_d;

		}

	z = xa > xa0 ? *--xa : 0;

	if (k -= Ebits) {

		d0 = Exp_1 | y << k | z >> 32 - k;

		y = xa > xa0 ? *--xa : 0;

		d1 = z << k | y >> 32 - k;

		}

	else {

		d0 = Exp_1 | y;

		d1 = z;

		}

#else

	if (k < Ebits + 16) {

		z = xa > xa0 ? *--xa : 0;

		d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;

		w = xa > xa0 ? *--xa : 0;

		y = xa > xa0 ? *--xa : 0;

		d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;

		goto ret_d;

		}

	z = xa > xa0 ? *--xa : 0;

	w = xa > xa0 ? *--xa : 0;

	k -= Ebits + 16;

	d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;

	y = xa > xa0 ? *--xa : 0;

	d1 = w << k + 16 | y << k;

#endif

 ret_d:

#ifdef VAX

	word0(d) = d0 >> 16 | d0 << 16;

	word1(d) = d1 >> 16 | d1 << 16;

#else

#undef d0

#undef d1

#endif

	return d.d;

	}

 static double

ratio(Bigint *a, Bigint *b)

{

	Dul da, db;

	int k, ka, kb;

	da.d = b2d(a, &ka);

	db.d = b2d(b, &kb);

#ifdef Pack_32

	k = ka - kb + 32*(a->wds - b->wds);

#else

	k = ka - kb + 16*(a->wds - b->wds);

#endif

#ifdef IBM

	if (k > 0) {

		word0(da) += (k >> 2)*Exp_msk1;

		if (k &= 3)

			da *= 1 << k;

		}

	else {

		k = -k;

		word0(db) += (k >> 2)*Exp_msk1;

		if (k &= 3)

			db *= 1 << k;

		}

#else

	if (k > 0)

		word0(da) += k*Exp_msk1;

	else {

		k = -k;

		word0(db) += k*Exp_msk1;

		}

#endif

	return da.d / db.d;

	}

 double

strtod(CONST char *s00, char **se)

{

	int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,

		 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;

	CONST char *s, *s0, *s1;

	double aadj, aadj1, adj;

	Dul rv, rv0;

	long L;

	unsigned long y, z;

	Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;

	sign = nz0 = nz = 0;

	rv.d = 0.;

	for(s = s00;;s++) switch(*s) {

		case '-':

			sign = 1;

			/* no break */

		case '+':

			if (*++s)

				goto break2;

			/* no break */

		case 0:

			s = s00;

			goto ret;

		case '\t':

		case '\n':

		case '\v':

		case '\f':

		case '\r':

		case ' ':

			continue;

		default:

			goto break2;

		}

 break2:

	if (*s == '0') {

		nz0 = 1;

		while(*++s == '0') ;

		if (!*s)

			goto ret;

		}

	s0 = s;

	y = z = 0;

	for(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;

	if (c == '.') {

		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) {

			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') {

				e = c - '0';

				s1 = s;

				while((c = *++s) >= '0' && c <= '9')

					e = 10*e + c - '0';

				if (s - s1 > 8)

					/* Avoid confusion from exponents

					 * so large that e might overflow.

					 */

					e = 9999999;

				if (esign)

					e = -e;

				}

			else

				e = 0;

			}

		else

			s = s00;

		}

	if (!nd) {

		if (!nz && !nz0)

			s = s00;

		goto ret;

		}

	e1 = e -= nf;

	/* 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;

	rv.d = y;

	if (k > 9)

		rv.d = tens[k - 9] * rv.d + z;

	bd0 = 0;

	if (nd <= DBL_DIG

#ifndef RND_PRODQUOT

		&& FLT_ROUNDS == 1

#endif

			) {

		if (!e)

			goto ret;

		if (e > 0) {

			if (e <= Ten_pmax) {

#ifdef VAX

				goto vax_ovfl_check;

#else

				/* rv = */ rounded_product(rv.d, tens[e]);

				goto ret;

#endif

				}

			i = DBL_DIG - nd;

			if (e <= Ten_pmax + i) {

				/* A fancier test would sometimes let us do

				 * this for larger i values.

				 */

				e -= i;

				rv.d *= tens[i];

#ifdef VAX

				/* VAX exponent range is so narrow we must

				 * worry about overflow here...

				 */

 vax_ovfl_check:

				word0(rv) -= P*Exp_msk1;

				/* rv = */ rounded_product(rv.d, tens[e]);

				if ((word0(rv) & Exp_mask)

				 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))

					goto ovfl;

				word0(rv) += P*Exp_msk1;

#else

				/* rv = */ rounded_product(rv.d, tens[e]);

#endif

				goto ret;

				}

			}

		else if (e >= -Ten_pmax) {

			/* rv = */ rounded_quotient(rv.d, tens[-e]);

			goto ret;

			}

		}

	e1 += nd - k;

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

	if (e1 > 0) {

		if (nd0 + e1 - 1 > DBL_MAX_10_EXP)

			goto ovfl;

		if (i = e1 & 15)

			rv.d *= tens[i];

		if (e1 &= ~15) {

			if (e1 > DBL_MAX_10_EXP) {

 ovfl:

				errno = ERANGE;

				rv.d = HUGE_VAL;

				if (bd0)

					goto retfree;

				goto ret;

				}

			if (e1 >>= 4) {

				for(j = 0; e1 > 1; j++, e1 >>= 1)

					if (e1 & 1)

						rv.d *= bigtens[j];

			/* The last multiplication could overflow. */

				word0(rv) -= P*Exp_msk1;

				rv.d *= bigtens[j];

				if ((z = word0(rv) & Exp_mask)

				 > Exp_msk1*(DBL_MAX_EXP+Bias-P))

					goto ovfl;

				if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {

					/* set to largest number */

					/* (Can't trust DBL_MAX) */

					word0(rv) = Big0;

					word1(rv) = Big1;

					}

				else

					word0(rv) += P*Exp_msk1;

				}

			}

		}

	else if (e1 < 0) {

		e1 = -e1;

		if (i = e1 & 15)

			rv.d /= tens[i];

		if (e1 &= ~15) {

			e1 >>= 4;

			if (e1 >= 1 << n_bigtens)

				goto undfl;

			for(j = 0; e1 > 1; j++, e1 >>= 1)

				if (e1 & 1)

					rv.d *= tinytens[j];

			/* The last multiplication could underflow. */

			rv0.d = rv.d;

			rv.d *= tinytens[j];

			if (rv.d == 0) {

				rv.d = 2.*rv0.d;

				rv.d *= tinytens[j];

				if (rv.d == 0) {

 undfl:

					rv.d = 0.;

					errno = ERANGE;

					if (bd0)

						goto retfree;

					goto ret;

					}

				word0(rv) = Tiny0;

				word1(rv) = Tiny1;

				/* The refinement below will clean

				 * this approximation up.

				 */

				}

			}

		}

	/* 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 = d2b(rv.d, &bbe, &bbbits);	/* rv = bb * 2^bbe */

		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;

#ifdef Sudden_Underflow

#ifdef IBM

		j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);

#else

		j = P + 1 - bbbits;

#endif

#else

		i = bbe + bbbits - 1;	/* logb(rv) */

		if (i < Emin)	/* denormal */

			j = bbe + (P-Emin);

		else

			j = P + 1 - bbbits;

#endif

		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;

			}

		if (bb2 > 0)

			bb = lshift(bb, bb2);

		if (bd5 > 0)

			bd = pow5mult(bd, bd5);

		if (bd2 > 0)

			bd = lshift(bd, bd2);

		if (bs2 > 0)

			bs = lshift(bs, bs2);

		delta = diff(bb, bd);

		dsign = delta->sign;

		delta->sign = 0;

		i = cmp(delta, bs);

		if (i < 0) {

			/* Error is less than half an ulp -- check for

			 * special case of mantissa a power of two.

			 */

			if (dsign || word1(rv) || word0(rv) & Bndry_mask)

				break;

			delta = lshift(delta,Log2P);

			if (cmp(delta, bs) > 0)

				goto drop_down;

			break;

			}

		if (i == 0) {

			/* exactly half-way between */

			if (dsign) {

				if ((word0(rv) & Bndry_mask1) == Bndry_mask1

				 &&  word1(rv) == 0xffffffff) {

					/*boundary case -- increment exponent*/

					word0(rv) = (word0(rv) & Exp_mask)

						+ Exp_msk1

#ifdef IBM

						| Exp_msk1 >> 4

#endif

						;

					word1(rv) = 0;

					break;

					}

				}

			else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {

 drop_down:

				/* boundary case -- decrement exponent */

#ifdef Sudden_Underflow

				L = word0(rv) & Exp_mask;

#ifdef IBM

				if (L <  Exp_msk1)

#else

				if (L <= Exp_msk1)

#endif

					goto undfl;

				L -= Exp_msk1;

#else

				L = (word0(rv) & Exp_mask) - Exp_msk1;

#endif

				word0(rv) = L | Bndry_mask1;

				word1(rv) = 0xffffffff;

#ifdef IBM

				goto cont;

#else

				break;

#endif

				}

#ifndef ROUND_BIASED

			if (!(word1(rv) & LSB))

				break;

#endif

			if (dsign)

				rv.d += ulp(rv.d);

#ifndef ROUND_BIASED

			else {

				rv.d -= ulp(rv.d);

#ifndef Sudden_Underflow

				if (rv.d == 0)

					goto undfl;

#endif

				}

#endif

			break;

			}

		if ((aadj = ratio(delta, bs)) <= 2.) {

			if (dsign)

				aadj = aadj1 = 1.;

			else if (word1(rv) || word0(rv) & Bndry_mask) {

#ifndef Sudden_Underflow

				if (word1(rv) == Tiny1 && !word0(rv))

					goto undfl;

#endif

				aadj = 1.;

				aadj1 = -1.;

				}

			else {

				/* special case -- power of FLT_RADIX to be */

				/* rounded down... */

				if (aadj < 2./FLT_RADIX)

					aadj = 1./FLT_RADIX;

				else

					aadj *= 0.5;

				aadj1 = -aadj;

				}

			}

		else {

			aadj *= 0.5;

			aadj1 = dsign ? aadj : -aadj;

#ifdef Check_FLT_ROUNDS

			switch(FLT_ROUNDS) {

				case 2: /* towards +infinity */

					aadj1 -= 0.5;

					break;

				case 0: /* towards 0 */

				case 3: /* towards -infinity */

					aadj1 += 0.5;

				}

#else

			if (FLT_ROUNDS == 0)

				aadj1 += 0.5;

#endif

			}

		y = word0(rv) & Exp_mask;

		/* Check for overflow */

		if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {

			rv0.d = rv.d;

			word0(rv) -= P*Exp_msk1;

			adj = aadj1 * ulp(rv.d);

			rv.d += adj;

			if ((word0(rv) & Exp_mask) >=

					Exp_msk1*(DBL_MAX_EXP+Bias-P)) {

				if (word0(rv0) == Big0 && word1(rv0) == Big1)

					goto ovfl;

				word0(rv) = Big0;

				word1(rv) = Big1;

				goto cont;

				}

			else

				word0(rv) += P*Exp_msk1;

			}

		else {

#ifdef Sudden_Underflow

			if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {

				rv0.d = rv.d;

				word0(rv) += P*Exp_msk1;

				adj = aadj1 * ulp(rv.d);

				rv.d += adj;

#ifdef IBM

				if ((word0(rv) & Exp_mask) <  P*Exp_msk1)

#else

				if ((word0(rv) & Exp_mask) <= P*Exp_msk1)

#endif

					{

					if (word0(rv0) == Tiny0

					 && word1(rv0) == Tiny1)

						goto undfl;

					word0(rv) = Tiny0;

					word1(rv) = Tiny1;

					goto cont;

					}

				else

					word0(rv) -= P*Exp_msk1;

				}

			else {

				adj = aadj1 * ulp(rv.d);

				rv.d += adj;

				}

#else

			/* Compute adj so that the IEEE rounding rules will

			 * correctly round rv + adj in some half-way cases.

			 * If rv * ulp(rv) is denormalized (i.e.,

			 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid

			 * trouble from bits lost to denormalization;

			 * example: 1.2e-307 .

			 */

			if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {

				aadj1 = (double)(int)(aadj + 0.5);

				if (!dsign)

					aadj1 = -aadj1;

				}

			adj = aadj1 * ulp(rv.d);

			rv.d += adj;

#endif

			}

		z = word0(rv) & Exp_mask;

		if (y == z) {

			/* Can we stop now? */

			L = aadj;

			aadj -= L;

			/* The tolerances below are conservative. */

			if (dsign || word1(rv) || word0(rv) & Bndry_mask) {

				if (aadj < .4999999 || aadj > .5000001)

					break;

				}

			else if (aadj < .4999999/FLT_RADIX)

				break;

			}

 cont:

		Bfree(bb);

		Bfree(bd);

		Bfree(bs);

		Bfree(delta);

		}

 retfree:

	Bfree(bb);

	Bfree(bd);

	Bfree(bs);

	Bfree(bd0);

	Bfree(delta);

 ret:

	if (se)

		*se = (char *)s;

	return sign ? -rv.d : rv.d;

	}