ref: 93cc0584dfd56f0d65ff4f1185ba5f85cc4656da
parent: 1f137b8e48ba8345cc7a45494009b1cf69e6a47b
author: Tor Andersson <[email protected]>
date: Mon Sep 1 12:56:38 EDT 2014
Add portable strtod implementation. Ignore locale problems with strtod.
--- a/jsbuiltin.c
+++ b/jsbuiltin.c
@@ -74,7 +74,7 @@
else if (!strncmp(s, "-Infinity", 9))
js_pushnumber(J, -INFINITY);
else {- n = js_strtod(s, &e);
+ n = js_stringtofloat(s, &e);
if (e == s)
js_pushnumber(J, NAN);
else
--- a/jsdtoa.c
+++ b/jsdtoa.c
@@ -18,6 +18,10 @@
#include <stdlib.h>
#include <errno.h>
+#include "jsi.h"
+
+typedef unsigned long ulong;
+
enum { NSIGNIF = 17 };/*
@@ -44,7 +48,7 @@
1e150, 1e151, 1e152, 1e153, 1e154, 1e155, 1e156, 1e157, 1e158, 1e159,
};
#define npows10 ((int)(sizeof(pows10)/sizeof(pows10[0])))
-#define pow10(x) fmtpow10(x)
+#define pow10(x) fmtpow10(x)
static double
pow10(int n)
@@ -101,7 +105,7 @@
* need to overflow adding digit.
* shift number down and insert 1 at beginning.
* decimal is known to be 0s or we wouldn't
- * have gotten this far. (e.g., 99999+1 => 00000)
+ * have gotten this far. (e.g., 99999+1 => 00000)
*/
a[0] = '1';
return 1;
@@ -111,7 +115,7 @@
* subtract 1 from the decimal integer string a.
* if 10000 underflows into 09999, make it 99999
* and return 1 to tell caller to move the virtual
- * decimal point. this way, xsub1 is inverse of xadd1.
+ * decimal point. this way, xsub1 is inverse of xadd1.
*/
static int
xsub1(char *a, int n)
@@ -128,7 +132,7 @@
/*
* just zeroed the top digit; shift everyone up.
* decimal is known to be 9s or we wouldn't
- * have gotten this far. (e.g., 10000-1 => 09999)
+ * have gotten this far. (e.g., 10000-1 => 09999)
*/
*b = '9';
return 1;
@@ -139,7 +143,7 @@
*b = '9';
}
/*
- * can't get here. the number a is always normalized
+ * can't get here. the number a is always normalized
* so that it has a nonzero first digit.
*/
return 0;
@@ -149,7 +153,7 @@
* format exponent like sprintf(p, "e%+d", e)
*/
void
-jsV_fmtexp(char *p, int e)
+js_fmtexp(char *p, int e)
{char se[9];
int i;
@@ -179,7 +183,7 @@
* assumes special cases (NaN, +Inf, -Inf) have been handled.
*/
void
-jsV_dtoa(double f, char *s, int *exp, int *neg, int *ns)
+js_dtoa(double f, char *s, int *exp, int *neg, int *ns)
{int c, d, e2, e, ee, i, ndigit, oerrno;
char tmp[NSIGNIF+10];
@@ -237,18 +241,18 @@
* adjust e because s is 314159... not 3.14159...
*/
e -= NSIGNIF-1;
- jsV_fmtexp(s+NSIGNIF, e);
+ js_fmtexp(s+NSIGNIF, e);
/*
* adjust conversion until strtod(s) == f exactly.
*/
for(i=0; i<10; i++) {- g = strtod(s, NULL);
+ g = js_strtod(s, NULL);
if(f > g) { if(xadd1(s, NSIGNIF)) {/* gained a digit */
e--;
- jsV_fmtexp(s+NSIGNIF, e);
+ js_fmtexp(s+NSIGNIF, e);
}
continue;
}
@@ -256,7 +260,7 @@
if(xsub1(s, NSIGNIF)) {/* lost a digit */
e++;
- jsV_fmtexp(s+NSIGNIF, e);
+ js_fmtexp(s+NSIGNIF, e);
}
continue;
}
@@ -274,7 +278,7 @@
c = s[i];
if(c != '9') {s[i] = '9';
- g = strtod(s, NULL);
+ g = js_strtod(s, NULL);
if(g != f) {s[i] = c;
break;
@@ -290,9 +294,9 @@
ee = e;
if(xadd1(tmp, NSIGNIF)) {ee--;
- jsV_fmtexp(tmp+NSIGNIF, ee);
+ js_fmtexp(tmp+NSIGNIF, ee);
}
- g = strtod(tmp, NULL);
+ g = js_strtod(tmp, NULL);
if(g == f) {strcpy(s, tmp);
e = ee;
@@ -306,7 +310,7 @@
c = s[i];
if(c != '0') {s[i] = '0';
- g = strtod(s, NULL);
+ g = js_strtod(s, NULL);
if(g != f) {s[i] = c;
break;
@@ -326,4 +330,514 @@
*exp = e;
*ns = ndigit;
errno = oerrno;
+}
+
+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
+js_strtod(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;
+ 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 != NULL)
+ *aas = s+3;
+ goto retnan;
+ }
+ case S1:
+ if(xcmp(s, "infinity") == 0) {+ if(aas != NULL)
+ *aas = s+8;
+ goto retinf;
+ }
+ if(xcmp(s, "inf") == 0) {+ if(aas != NULL)
+ *aas = s+3;
+ goto retinf;
+ }
+ case S3:
+ if(aas != NULL)
+ *aas = (char*)as;
+ goto ret0; /* no digits found */
+ case S6:
+ s--; /* back over +- */
+ case S5:
+ s--; /* back over e */
+ break;
+ }
+ if(aas != NULL)
+ *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(c2 >= 'A' && c2 <= 'Z')
+ c2 = c2 - 'A' + 'a';
+ if(c1 != c2)
+ return 1;
+ }
+ return 0;
}
--- a/jsi.h
+++ b/jsi.h
@@ -58,6 +58,12 @@
void jsS_dumpstrings(js_State *J);
void jsS_freestrings(js_State *J);
+/* Portable strtod and printf float formatting */
+
+void js_fmtexp(char *p, int e);
+void js_dtoa(double f, char *digits, int *exp, int *neg, int *ndigits);
+double js_strtod(const char *as, char **aas);
+
/* Private stack functions */
void js_newfunction(js_State *J, js_Function *function, js_Environment *scope);
--- a/jslex.c
+++ b/jslex.c
@@ -310,7 +310,7 @@
if (jsY_isidentifierstart(PEEK))
jsY_error(J, "number with letter suffix");
- J->number = strtod(s, NULL);
+ J->number = js_strtod(s, NULL);
return TK_NUMBER;
}
--- a/jsvalue.c
+++ b/jsvalue.c
@@ -122,7 +122,7 @@
}
}
-double js_strtod(const char *s, char **ep)
+double js_stringtofloat(const char *s, char **ep)
{char *end;
double n;
@@ -136,7 +136,7 @@
if (*e == '+' || *e == '-') ++e;
while (*e >= '0' && *e <= '9') ++e;
}
- n = strtod(s, &end);
+ n = js_strtod(s, &end);
if (end == e) {*ep = (char*)e;
return n;
@@ -160,7 +160,7 @@
else if (!strncmp(s, "-Infinity", 9))
n = -INFINITY, e = (char*)s + 9;
else
- n = js_strtod(s, &e);
+ n = js_stringtofloat(s, &e);
while (jsY_iswhite(*e) || jsY_isnewline(*e)) ++e;
if (*e) return NAN;
return n;
@@ -199,7 +199,7 @@
if (isinf(f)) return f < 0 ? "-Infinity" : "Infinity";
if (f == 0) return "0";
- jsV_dtoa(f, digits, &exp, &neg, &ndigits);
+ js_dtoa(f, digits, &exp, &neg, &ndigits);
point = ndigits + exp;
if (neg)
@@ -213,7 +213,7 @@
while (n--)
*p++ = *s++;
}
- jsV_fmtexp(p, point - 1);
+ js_fmtexp(p, point - 1);
}
else if (point <= 0) {--- a/jsvalue.h
+++ b/jsvalue.h
@@ -130,7 +130,7 @@
js_Object *jsV_toobject(js_State *J, const js_Value *v);
js_Value jsV_toprimitive(js_State *J, const js_Value *v, int preferred);
-double js_strtod(const char *s, char **ep);
+double js_stringtofloat(const char *s, char **ep);
double jsV_numbertointeger(double n);
int jsV_numbertoint32(double n);
unsigned int jsV_numbertouint32(double n);
@@ -138,10 +138,6 @@
unsigned short jsV_numbertouint16(double n);
const char *jsV_numbertostring(js_State *J, double number);
double jsV_stringtonumber(js_State *J, const char *string);
-
-/* jsdtoa.c */
-void jsV_fmtexp(char *p, int e);
-void jsV_dtoa(double f, char *digits, int *exp, int *neg, int *ndigits);
/* jsproperty.c */
js_Object *jsV_newobject(js_State *J, enum js_Class type, js_Object *prototype);