Import sources from 2011-03-30 iso image

1 files changed, 731 insertions, 0 deletions

diff --git a/sys/src/ape/lib/ap/stdio/strtod.c b/sys/src/ape/lib/ap/stdio/strtod.c
new file mode 100755
index 000000000..ff820ef62
--- /dev/null
+++ b/sys/src/ape/lib/ap/stdio/strtod.c
@@ -0,0 +1,731 @@
+#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;
+	}