fiss

strtod.c (9133B)

---

 /* Copyright (c) 2002-2006 Lucent Technologies; see LICENSE */
 #include "common.h"
 #include "fmt.h"
 #include "fmtdef.h"
 #include "plan9.h"
 
 #include <ctype.h>
 #include <errno.h>
 #include <math.h>
 #include <stdlib.h>
 #include <string.h>
 
 static ulong
 umuldiv(ulong a, ulong b, ulong c) {
 	double d;
 
 	d = ((double) a * (double) b) / (double) c;
 	if (d >= 4294967295.)
 		d = 4294967295.;
 	return (ulong) d;
 }
 
 /*
  * This routine will convert to arbitrary precision
  * floating point entirely in multi-precision fixed.
  * The answer is the closest floating point number to
  * the given decimal number. Exactly half way are
  * rounded ala ieee rules.
  * Method is to scale input decimal between .500 and .999...
  * with external power of 2, then binary search for the
  * closest mantissa to this decimal number.
  * Nmant is is the required precision. (53 for ieee dp)
  * Nbits is the max number of bits/word. (must be <= 28)
  * Prec is calculated - the number of words of fixed mantissa.
  */
 enum {
 	Nbits  = 28,                          /* bits safely represented in a ulong */
 	Nmant  = 53,                          /* bits of precision required */
 	Prec   = (Nmant + Nbits + 1) / Nbits, /* words of Nbits each to represent mantissa */
 	Sigbit = 1 << (Prec * Nbits - Nmant), /* first significant bit of Prec-th word */
 	Ndig   = 1500,
 	One    = (ulong) (1 << Nbits),
 	Half   = (ulong) (One >> 1),
 	Maxe   = 310,
 
 	Fsign   = 1 << 0, /* found - */
 	Fesign  = 1 << 1, /* found e- */
 	Fdpoint = 1 << 2, /* found . */
 
 	S0 = 0, /* _		_S0	+S1	#S2	.S3 */
 	S1,     /* _+		#S2	.S3 */
 	S2,     /* _+#		#S2	.S4	eS5 */
 	S3,     /* _+.		#S4 */
 	S4,     /* _+#.#	#S4	eS5 */
 	S5,     /* _+#.#e	+S6	#S7 */
 	S6,     /* _+#.#e+	#S7 */
 	S7      /* _+#.#e+#	#S7 */
 };
 
 static int  xcmp(char*, char*);
 static int  fpcmp(char*, ulong*);
 static void frnorm(ulong*);
 static void divascii(char*, int*, int*, int*);
 static void mulascii(char*, int*, int*, int*);
 
 typedef struct Tab Tab;
 struct Tab {
 	int   bp;
 	int   siz;
 	char* cmp;
 };
 
 double
 fmtstrtod(const char* as, char** aas) {
 	int    na, ex, dp, bp, c, i, flag, state;
 	ulong  low[Prec], hig[Prec], mid[Prec];
 	double d;
 	char * s, a[Ndig];
 
 	flag = 0; /* Fsign, Fesign, Fdpoint */
 	na   = 0; /* number of digits of a[] */
 	dp   = 0; /* na of decimal point */
 	ex   = 0; /* exonent */
 
 	state = S0;
 	for (s = (char*) as;; s++) {
 		c = *s;
 		if (c >= '0' && c <= '9') {
 			switch (state) {
 				case S0:
 				case S1:
 				case S2:
 					state = S2;
 					break;
 				case S3:
 				case S4:
 					state = S4;
 					break;
 
 				case S5:
 				case S6:
 				case S7:
 					state = S7;
 					ex    = ex * 10 + (c - '0');
 					continue;
 			}
 			if (na == 0 && c == '0') {
 				dp--;
 				continue;
 			}
 			if (na < Ndig - 50)
 				a[na++] = c;
 			continue;
 		}
 		switch (c) {
 			case '\t':
 			case '\n':
 			case '\v':
 			case '\f':
 			case '\r':
 			case ' ':
 				if (state == S0)
 					continue;
 				break;
 			case '-':
 				if (state == S0)
 					flag |= Fsign;
 				else
 					flag |= Fesign;
 				FALLTHROUGH;
 			case '+':
 				if (state == S0)
 					state = S1;
 				else if (state == S5)
 					state = S6;
 				else
 					break; /* syntax */
 				continue;
 			case '.':
 				flag |= Fdpoint;
 				dp = na;
 				if (state == S0 || state == S1) {
 					state = S3;
 					continue;
 				}
 				if (state == S2) {
 					state = S4;
 					continue;
 				}
 				break;
 			case 'e':
 			case 'E':
 				if (state == S2 || state == S4) {
 					state = S5;
 					continue;
 				}
 				break;
 		}
 		break;
 	}
 
 	/*
 	 * clean up return char-pointer
 	 */
 	switch (state) {
 		case S0:
 			if (xcmp(s, "nan") == 0) {
 				if (aas != nil)
 					*aas = s + 3;
 				goto retnan;
 			}
 			FALLTHROUGH;
 		case S1:
 			if (xcmp(s, "infinity") == 0) {
 				if (aas != nil)
 					*aas = s + 8;
 				goto retinf;
 			}
 			if (xcmp(s, "inf") == 0) {
 				if (aas != nil)
 					*aas = s + 3;
 				goto retinf;
 			}
 			FALLTHROUGH;
 		case S3:
 			if (aas != nil)
 				*aas = (char*) as;
 			goto ret0; /* no digits found */
 			FALLTHROUGH;
 		case S6:
 			s--; /* back over +- */
 			FALLTHROUGH;
 		case S5:
 			s--; /* back over e */
 			break;
 	}
 	if (aas != nil)
 		*aas = s;
 
 	if (flag & Fdpoint)
 		while (na > 0 && a[na - 1] == '0')
 			na--;
 	if (na == 0)
 		goto ret0; /* zero */
 	a[na] = 0;
 	if (!(flag & Fdpoint))
 		dp = na;
 	if (flag & Fesign)
 		ex = -ex;
 	dp += ex;
 	if (dp < -Maxe) {
 		errno = ERANGE;
 		goto ret0; /* underflow by exp */
 	} else if (dp > +Maxe)
 		goto retinf; /* overflow by exp */
 
 	/*
 	 * normalize the decimal ascii number
 	 * to range .[5-9][0-9]* e0
 	 */
 	bp = 0; /* binary exponent */
 	while (dp > 0)
 		divascii(a, &na, &dp, &bp);
 	while (dp < 0 || a[0] < '5')
 		mulascii(a, &na, &dp, &bp);
 
 	/* close approx by naive conversion */
 	mid[0] = 0;
 	mid[1] = 1;
 	for (i = 0; (c = a[i]) != '\0'; i++) {
 		mid[0] = mid[0] * 10 + (c - '0');
 		mid[1] = mid[1] * 10;
 		if (i >= 8)
 			break;
 	}
 	low[0] = umuldiv(mid[0], One, mid[1]);
 	hig[0] = umuldiv(mid[0] + 1, One, mid[1]);
 	for (i = 1; i < Prec; i++) {
 		low[i] = 0;
 		hig[i] = One - 1;
 	}
 
 	/* binary search for closest mantissa */
 	for (;;) {
 		/* mid = (hig + low) / 2 */
 		c = 0;
 		for (i = 0; i < Prec; i++) {
 			mid[i] = hig[i] + low[i];
 			if (c)
 				mid[i] += One;
 			c = mid[i] & 1;
 			mid[i] >>= 1;
 		}
 		frnorm(mid);
 
 		/* compare */
 		c = fpcmp(a, mid);
 		if (c > 0) {
 			c = 1;
 			for (i = 0; i < Prec; i++)
 				if (low[i] != mid[i]) {
 					c      = 0;
 					low[i] = mid[i];
 				}
 			if (c)
 				break; /* between mid and hig */
 			continue;
 		}
 		if (c < 0) {
 			for (i = 0; i < Prec; i++)
 				hig[i] = mid[i];
 			continue;
 		}
 
 		/* only hard part is if even/odd roundings wants to go up */
 		c = mid[Prec - 1] & (Sigbit - 1);
 		if (c == Sigbit / 2 && (mid[Prec - 1] & Sigbit) == 0)
 			mid[Prec - 1] -= c;
 		break; /* exactly mid */
 	}
 
 	/* normal rounding applies */
 	c = mid[Prec - 1] & (Sigbit - 1);
 	mid[Prec - 1] -= c;
 	if (c >= Sigbit / 2) {
 		mid[Prec - 1] += Sigbit;
 		frnorm(mid);
 	}
 	goto out;
 
 ret0:
 	return 0;
 
 retnan:
 	return __NaN();
 
 retinf:
 	/*
 	 * Unix strtod requires these.  Plan 9 would return Inf(0) or Inf(-1). */
 	errno = ERANGE;
 	if (flag & Fsign)
 		return -HUGE_VAL;
 	return HUGE_VAL;
 
 out:
 	d = 0;
 	for (i = 0; i < Prec; i++)
 		d = d * One + mid[i];
 	if (flag & Fsign)
 		d = -d;
 	d = ldexp(d, bp - Prec * Nbits);
 	if (d == 0) { /* underflow */
 		errno = ERANGE;
 	}
 	return d;
 }
 
 static void
 frnorm(ulong* f) {
 	int i, c;
 
 	c = 0;
 	for (i = Prec - 1; i > 0; i--) {
 		f[i] += c;
 		c = f[i] >> Nbits;
 		f[i] &= One - 1;
 	}
 	f[0] += c;
 }
 
 static int
 fpcmp(char* a, ulong* f) {
 	ulong tf[Prec];
 	int   i, d, c;
 
 	for (i = 0; i < Prec; i++)
 		tf[i] = f[i];
 
 	for (;;) {
 		/* tf *= 10 */
 		for (i = 0; i < Prec; i++)
 			tf[i] = tf[i] * 10;
 		frnorm(tf);
 		d = (tf[0] >> Nbits) + '0';
 		tf[0] &= One - 1;
 
 		/* compare next digit */
 		c = *a;
 		if (c == 0) {
 			if ('0' < d)
 				return -1;
 			if (tf[0] != 0)
 				goto cont;
 			for (i = 1; i < Prec; i++)
 				if (tf[i] != 0)
 					goto cont;
 			return 0;
 		}
 		if (c > d)
 			return +1;
 		if (c < d)
 			return -1;
 		a++;
 	cont:;
 	}
 }
 
 static void
 divby(char* a, int* na, int b) {
 	int   n, c;
 	char* p;
 
 	p = a;
 	n = 0;
 	while (n >> b == 0) {
 		c = *a++;
 		if (c == 0) {
 			while (n) {
 				c = n * 10;
 				if (c >> b)
 					break;
 				n = c;
 			}
 			goto xx;
 		}
 		n = n * 10 + c - '0';
 		(*na)--;
 	}
 	for (;;) {
 		c = n >> b;
 		n -= c << b;
 		*p++ = c + '0';
 		c    = *a++;
 		if (c == 0)
 			break;
 		n = n * 10 + c - '0';
 	}
 	(*na)++;
 xx:
 	while (n) {
 		n = n * 10;
 		c = n >> b;
 		n -= c << b;
 		*p++ = c + '0';
 		(*na)++;
 	}
 	*p = 0;
 }
 
 static Tab tab1[] = {
 	{ 1, 0, "" },
 	{ 3, 1, "7" },
 	{ 6, 2, "63" },
 	{ 9, 3, "511" },
 	{ 13, 4, "8191" },
 	{ 16, 5, "65535" },
 	{ 19, 6, "524287" },
 	{ 23, 7, "8388607" },
 	{ 26, 8, "67108863" },
 	{ 27, 9, "134217727" }
 };
 
 static void
 divascii(char* a, int* na, int* dp, int* bp) {
 	int  b, d;
 	Tab* t;
 
 	d = *dp;
 	if (d >= (int) (nelem(tab1)))
 		d = (int) (nelem(tab1)) - 1;
 	t = tab1 + d;
 	b = t->bp;
 	if (memcmp(a, t->cmp, t->siz) > 0)
 		d--;
 	*dp -= d;
 	*bp += b;
 	divby(a, na, b);
 }
 
 static void
 mulby(char* a, char* p, char* q, int b) {
 	int n, c;
 
 	n  = 0;
 	*p = 0;
 	for (;;) {
 		q--;
 		if (q < a)
 			break;
 		c = *q - '0';
 		c = (c << b) + n;
 		n = c / 10;
 		c -= n * 10;
 		p--;
 		*p = c + '0';
 	}
 	while (n) {
 		c = n;
 		n = c / 10;
 		c -= n * 10;
 		p--;
 		*p = c + '0';
 	}
 }
 
 static Tab tab2[] = {
 	{ 1, 1, "" }, /* dp = 0-0 */
 	{ 3, 3, "125" },
 	{ 6, 5, "15625" },
 	{ 9, 7, "1953125" },
 	{ 13, 10, "1220703125" },
 	{ 16, 12, "152587890625" },
 	{ 19, 14, "19073486328125" },
 	{ 23, 17, "11920928955078125" },
 	{ 26, 19, "1490116119384765625" },
 	{ 27, 19, "7450580596923828125" }, /* dp 8-9 */
 };
 
 static void
 mulascii(char* a, int* na, int* dp, int* bp) {
 	char* p;
 	int   d, b;
 	Tab*  t;
 
 	d = -*dp;
 	if (d >= (int) (nelem(tab2)))
 		d = (int) (nelem(tab2)) - 1;
 	t = tab2 + d;
 	b = t->bp;
 	if (memcmp(a, t->cmp, t->siz) < 0)
 		d--;
 	p = a + *na;
 	*bp -= b;
 	*dp += d;
 	*na += d;
 	mulby(a, p + d, p, b);
 }
 
 static int
 xcmp(char* a, char* b) {
 	int c1, c2;
 
 	while ((c1 = *b++) != '\0') {
 		c2 = *a++;
 		if (isupper(c2))
 			c2 = tolower(c2);
 		if (c1 != c2)
 			return 1;
 	}
 	return 0;
 }