Add strtod(), wcstod(), wcstol() and friends.

This implementation of strtod() replaces the old as it implements
support for parsing hex constants which is needed by various tests.

17 files changed, 1114 insertions, 2 deletions

diff --git a/system/lib/libc.symbols b/system/lib/libc.symbols
index 931e9e5f..d68b62cb 100644
--- a/system/lib/libc.symbols
+++ b/system/lib/libc.symbols
@@ -39,6 +39,17 @@
          W _ZnajRKSt9nothrow_t
          W _Znwj
          W _ZnwjRKSt9nothrow_t
+         T __floatscan
+         T __overflow
+         T __seek_on_exit
+         T __shgetc
+         T __shlim
+         W __strtod_l
+         W __strtof_l
+         W __strtold_l
+         T __toread
+         W __towrite_used
+         T __uflow
          T _err
          T _errx
          T _verr
@@ -78,11 +89,13 @@
          W pvalloc
          W realloc
          W realloc_in_place
+         T scalbnl
          T strtod
-         T strtod_l
+         W strtod_l
          T strtof
+         W strtof_l
          T strtold
-         T strtold_l
+         W strtold_l
          W valloc
          W verr
          W verrx
diff --git a/system/lib/libc/musl/src/internal/floatscan.c b/system/lib/libc/musl/src/internal/floatscan.c
new file mode 100644
index 00000000..f6e331d4
--- /dev/null
+++ b/system/lib/libc/musl/src/internal/floatscan.c
@@ -0,0 +1,496 @@
+#include <stdint.h>
+#include <stdio.h>
+#include <math.h>
+#include <float.h>
+#include <limits.h>
+#include <errno.h>
+#include <ctype.h>
+
+#include "shgetc.h"
+#include "floatscan.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+
+#define LD_B1B_DIG 2
+#define LD_B1B_MAX 9007199, 254740991
+#define KMAX 128
+
+#else /* LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 */
+
+#define LD_B1B_DIG 3
+#define LD_B1B_MAX 18, 446744073, 709551615
+#define KMAX 2048
+
+#endif
+
+#define MASK (KMAX-1)
+
+#define CONCAT2(x,y) x ## y
+#define CONCAT(x,y) CONCAT2(x,y)
+
+static long long scanexp(FILE *f, int pok)
+{
+	int c;
+	int x;
+	long long y;
+	int neg = 0;
+	
+	c = shgetc(f);
+	if (c=='+' || c=='-') {
+		neg = (c=='-');
+		c = shgetc(f);
+		if (c-'0'>=10U && pok) shunget(f);
+	}
+	if (c-'0'>=10U) {
+		shunget(f);
+		return LLONG_MIN;
+	}
+	for (x=0; c-'0'<10U && x<INT_MAX/10; c = shgetc(f))
+		x = 10*x + c-'0';
+	for (y=x; c-'0'<10U && y<LLONG_MAX/100; c = shgetc(f))
+		y = 10*y + c-'0';
+	for (; c-'0'<10U; c = shgetc(f));
+	shunget(f);
+	return neg ? -y : y;
+}
+
+
+static long double decfloat(FILE *f, int c, int bits, int emin, int sign, int pok)
+{
+	uint32_t x[KMAX];
+	static const uint32_t th[] = { LD_B1B_MAX };
+	int i, j, k, a, z;
+	long long lrp=0, dc=0;
+	long long e10=0;
+	int lnz = 0;
+	int gotdig = 0, gotrad = 0;
+	int rp;
+	int e2;
+	int emax = -emin-bits+3;
+	int denormal = 0;
+	long double y;
+	long double frac=0;
+	long double bias=0;
+	static const int p10s[] = { 10, 100, 1000, 10000,
+		100000, 1000000, 10000000, 100000000 };
+
+	j=0;
+	k=0;
+
+	/* Don't let leading zeros consume buffer space */
+	for (; c=='0'; c = shgetc(f)) gotdig=1;
+	if (c=='.') {
+		gotrad = 1;
+		for (c = shgetc(f); c=='0'; c = shgetc(f)) gotdig=1, lrp--;
+	}
+
+	x[0] = 0;
+	for (; c-'0'<10U || c=='.'; c = shgetc(f)) {
+		if (c == '.') {
+			if (gotrad) break;
+			gotrad = 1;
+			lrp = dc;
+		} else if (k < KMAX-3) {
+			dc++;
+			if (c!='0') lnz = dc;
+			if (j) x[k] = x[k]*10 + c-'0';
+			else x[k] = c-'0';
+			if (++j==9) {
+				k++;
+				j=0;
+			}
+			gotdig=1;
+		} else {
+			dc++;
+			if (c!='0') x[KMAX-4] |= 1;
+		}
+	}
+	if (!gotrad) lrp=dc;
+
+	if (gotdig && (c|32)=='e') {
+		e10 = scanexp(f, pok);
+		if (e10 == LLONG_MIN) {
+			if (pok) {
+				shunget(f);
+			} else {
+				shlim(f, 0);
+				return 0;
+			}
+			e10 = 0;
+		}
+		lrp += e10;
+	} else if (c>=0) {
+		shunget(f);
+	}
+	if (!gotdig) {
+		errno = EINVAL;
+		shlim(f, 0);
+		return 0;
+	}
+
+	/* Handle zero specially to avoid nasty special cases later */
+	if (!x[0]) return sign * 0.0;
+
+	/* Optimize small integers (w/no exponent) and over/under-flow */
+	if (lrp==dc && dc<10 && (bits>30 || x[0]>>bits==0))
+		return sign * (long double)x[0];
+	if (lrp > -emin/2) {
+		errno = ERANGE;
+		return sign * LDBL_MAX * LDBL_MAX;
+	}
+	if (lrp < emin-2*LDBL_MANT_DIG) {
+		errno = ERANGE;
+		return sign * LDBL_MIN * LDBL_MIN;
+	}
+
+	/* Align incomplete final B1B digit */
+	if (j) {
+		for (; j<9; j++) x[k]*=10;
+		k++;
+		j=0;
+	}
+
+	a = 0;
+	z = k;
+	e2 = 0;
+	rp = lrp;
+
+	/* Optimize small to mid-size integers (even in exp. notation) */
+	if (lnz<9 && lnz<=rp && rp < 18) {
+		if (rp == 9) return sign * (long double)x[0];
+		if (rp < 9) return sign * (long double)x[0] / p10s[8-rp];
+		int bitlim = bits-3*(int)(rp-9);
+		if (bitlim>30 || x[0]>>bitlim==0)
+			return sign * (long double)x[0] * p10s[rp-10];
+	}
+
+	/* Align radix point to B1B digit boundary */
+	if (rp % 9) {
+		int rpm9 = rp>=0 ? rp%9 : rp%9+9;
+		int p10 = p10s[8-rpm9];
+		uint32_t carry = 0;
+		for (k=a; k!=z; k++) {
+			uint32_t tmp = x[k] % p10;
+			x[k] = x[k]/p10 + carry;
+			carry = 1000000000/p10 * tmp;
+			if (k==a && !x[k]) {
+				a = (a+1 & MASK);
+				rp -= 9;
+			}
+		}
+		if (carry) x[z++] = carry;
+		rp += 9-rpm9;
+	}
+
+	/* Upscale until desired number of bits are left of radix point */
+	while (rp < 9*LD_B1B_DIG || (rp == 9*LD_B1B_DIG && x[a]<th[0])) {
+		uint32_t carry = 0;
+		e2 -= 29;
+		for (k=(z-1 & MASK); ; k=(k-1 & MASK)) {
+			uint64_t tmp = ((uint64_t)x[k] << 29) + carry;
+			if (tmp > 1000000000) {
+				carry = tmp / 1000000000;
+				x[k] = tmp % 1000000000;
+			} else {
+				carry = 0;
+				x[k] = tmp;
+			}
+			if (k==(z-1 & MASK) && k!=a && !x[k]) z = k;
+			if (k==a) break;
+		}
+		if (carry) {
+			rp += 9;
+			a = (a-1 & MASK);
+			if (a == z) {
+				z = (z-1 & MASK);
+				x[z-1 & MASK] |= x[z];
+			}
+			x[a] = carry;
+		}
+	}
+
+	/* Downscale until exactly number of bits are left of radix point */
+	for (;;) {
+		uint32_t carry = 0;
+		int sh = 1;
+		for (i=0; i<LD_B1B_DIG; i++) {
+			k = (a+i & MASK);
+			if (k == z || x[k] < th[i]) {
+				i=LD_B1B_DIG;
+				break;
+			}
+			if (x[a+i & MASK] > th[i]) break;
+		}
+		if (i==LD_B1B_DIG && rp==9*LD_B1B_DIG) break;
+		/* FIXME: find a way to compute optimal sh */
+		if (rp > 9+9*LD_B1B_DIG) sh = 9;
+		e2 += sh;
+		for (k=a; k!=z; k=(k+1 & MASK)) {
+			uint32_t tmp = x[k] & (1<<sh)-1;
+			x[k] = (x[k]>>sh) + carry;
+			carry = (1000000000>>sh) * tmp;
+			if (k==a && !x[k]) {
+				a = (a+1 & MASK);
+				i--;
+				rp -= 9;
+			}
+		}
+		if (carry) {
+			if ((z+1 & MASK) != a) {
+				x[z] = carry;
+				z = (z+1 & MASK);
+			} else x[z-1 & MASK] |= 1;
+		}
+	}
+
+	/* Assemble desired bits into floating point variable */
+	for (y=i=0; i<LD_B1B_DIG; i++) {
+		if ((a+i & MASK)==z) x[(z=(z+1 & MASK))-1] = 0;
+		y = 1000000000.0L * y + x[a+i & MASK];
+	}
+
+	y *= sign;
+
+	/* Limit precision for denormal results */
+	if (bits > LDBL_MANT_DIG+e2-emin) {
+		bits = LDBL_MANT_DIG+e2-emin;
+		if (bits<0) bits=0;
+		denormal = 1;
+	}
+
+	/* Calculate bias term to force rounding, move out lower bits */
+	if (bits < LDBL_MANT_DIG) {
+		bias = copysignl(scalbn(1, 2*LDBL_MANT_DIG-bits-1), y);
+		frac = fmodl(y, scalbn(1, LDBL_MANT_DIG-bits));
+		y -= frac;
+		y += bias;
+	}
+
+	/* Process tail of decimal input so it can affect rounding */
+	if ((a+i & MASK) != z) {
+		uint32_t t = x[a+i & MASK];
+		if (t < 500000000 && (t || (a+i+1 & MASK) != z))
+			frac += 0.25*sign;
+		else if (t > 500000000)
+			frac += 0.75*sign;
+		else if (t == 500000000) {
+			if ((a+i+1 & MASK) == z)
+				frac += 0.5*sign;
+			else
+				frac += 0.75*sign;
+		}
+		if (LDBL_MANT_DIG-bits >= 2 && !fmodl(frac, 1))
+			frac++;
+	}
+
+	y += frac;
+	y -= bias;
+
+	if ((e2+LDBL_MANT_DIG & INT_MAX) > emax-5) {
+		if (fabs(y) >= CONCAT(0x1p, LDBL_MANT_DIG)) {
+			if (denormal && bits==LDBL_MANT_DIG+e2-emin)
+				denormal = 0;
+			y *= 0.5;
+			e2++;
+		}
+		if (e2+LDBL_MANT_DIG>emax || (denormal && frac))
+			errno = ERANGE;
+	}
+
+	return scalbnl(y, e2);
+}
+
+static long double hexfloat(FILE *f, int bits, int emin, int sign, int pok)
+{
+	uint32_t x = 0;
+	long double y = 0;
+	long double scale = 1;
+	long double bias = 0;
+	int gottail = 0, gotrad = 0, gotdig = 0;
+	long long rp = 0;
+	long long dc = 0;
+	long long e2 = 0;
+	int d;
+	int c;
+
+	c = shgetc(f);
+
+	/* Skip leading zeros */
+	for (; c=='0'; c = shgetc(f)) gotdig = 1;
+
+	if (c=='.') {
+		gotrad = 1;
+		c = shgetc(f);
+		/* Count zeros after the radix point before significand */
+		for (rp=0; c=='0'; c = shgetc(f), rp--) gotdig = 1;
+	}
+
+	for (; c-'0'<10U || (c|32)-'a'<6U || c=='.'; c = shgetc(f)) {
+		if (c=='.') {
+			if (gotrad) break;
+			rp = dc;
+			gotrad = 1;
+		} else {
+			gotdig = 1;
+			if (c > '9') d = (c|32)+10-'a';
+			else d = c-'0';
+			if (dc<8) {
+				x = x*16 + d;
+			} else if (dc < LDBL_MANT_DIG/4+1) {
+				y += d*(scale/=16);
+			} else if (d && !gottail) {
+				y += 0.5*scale;
+				gottail = 1;
+			}
+			dc++;
+		}
+	}
+	if (!gotdig) {
+		shunget(f);
+		if (pok) {
+			shunget(f);
+			if (gotrad) shunget(f);
+		} else {
+			shlim(f, 0);
+		}
+		return sign * 0.0;
+	}
+	if (!gotrad) rp = dc;
+	while (dc<8) x *= 16, dc++;
+	if ((c|32)=='p') {
+		e2 = scanexp(f, pok);
+		if (e2 == LLONG_MIN) {
+			if (pok) {
+				shunget(f);
+			} else {
+				shlim(f, 0);
+				return 0;
+			}
+			e2 = 0;
+		}
+	} else {
+		shunget(f);
+	}
+	e2 += 4*rp - 32;
+
+	if (!x) return sign * 0.0;
+	if (e2 > -emin) {
+		errno = ERANGE;
+		return sign * LDBL_MAX * LDBL_MAX;
+	}
+	if (e2 < emin-2*LDBL_MANT_DIG) {
+		errno = ERANGE;
+		return sign * LDBL_MIN * LDBL_MIN;
+	}
+
+	while (x < 0x80000000) {
+		if (y>=0.5) {
+			x += x + 1;
+			y += y - 1;
+		} else {
+			x += x;
+			y += y;
+		}
+		e2--;
+	}
+
+	if (bits > 32+e2-emin) {
+		bits = 32+e2-emin;
+		if (bits<0) bits=0;
+	}
+
+	if (bits < LDBL_MANT_DIG)
+		bias = copysignl(scalbn(1, 32+LDBL_MANT_DIG-bits-1), sign);
+
+	if (bits<32 && y && !(x&1)) x++, y=0;
+
+	y = bias + sign*(long double)x + sign*y;
+	y -= bias;
+
+	if (!y) errno = ERANGE;
+
+	return scalbnl(y, e2);
+}
+
+long double __floatscan(FILE *f, int prec, int pok)
+{
+	int sign = 1;
+	size_t i;
+	int bits;
+	int emin;
+	int c;
+
+	switch (prec) {
+	case 0:
+		bits = FLT_MANT_DIG;
+		emin = FLT_MIN_EXP-bits;
+		break;
+	case 1:
+		bits = DBL_MANT_DIG;
+		emin = DBL_MIN_EXP-bits;
+		break;
+	case 2:
+		bits = LDBL_MANT_DIG;
+		emin = LDBL_MIN_EXP-bits;
+		break;
+	default:
+		return 0;
+	}
+
+	while (isspace((c=shgetc(f))));
+
+	if (c=='+' || c=='-') {
+		sign -= 2*(c=='-');
+		c = shgetc(f);
+	}
+
+	for (i=0; i<8 && (c|32)=="infinity"[i]; i++)
+		if (i<7) c = shgetc(f);
+	if (i==3 || i==8 || (i>3 && pok)) {
+		if (i!=8) {
+			shunget(f);
+			if (pok) for (; i>3; i--) shunget(f);
+		}
+		return sign * INFINITY;
+	}
+	if (!i) for (i=0; i<3 && (c|32)=="nan"[i]; i++)
+		if (i<2) c = shgetc(f);
+	if (i==3) {
+		if (shgetc(f) != '(') {
+			shunget(f);
+			return NAN;
+		}
+		for (i=1; ; i++) {
+			c = shgetc(f);
+			if (c-'0'<10U || c-'A'<26U || c-'a'<26U || c=='_')
+				continue;
+			if (c==')') return NAN;
+			shunget(f);
+			if (!pok) {
+				errno = EINVAL;
+				shlim(f, 0);
+				return 0;
+			}
+			while (i--) shunget(f);
+			return NAN;
+		}
+		return NAN;
+	}
+
+	if (i) {
+		shunget(f);
+		errno = EINVAL;
+		shlim(f, 0);
+		return 0;
+	}
+
+	if (c=='0') {
+		c = shgetc(f);
+		if ((c|32) == 'x')
+			return hexfloat(f, bits, emin, sign, pok);
+		shunget(f);
+		c = '0';
+	}
+
+	return decfloat(f, c, bits, emin, sign, pok);
+}
diff --git a/system/lib/libc/musl/src/internal/floatscan.h b/system/lib/libc/musl/src/internal/floatscan.h
new file mode 100644
index 00000000..e027fa08
--- /dev/null
+++ b/system/lib/libc/musl/src/internal/floatscan.h
@@ -0,0 +1,8 @@
+#ifndef FLOATSCAN_H
+#define FLOATSCAN_H
+
+#include <stdio.h>
+
+long double __floatscan(FILE *, int, int);
+
+#endif
diff --git a/system/lib/libc/musl/src/internal/intscan.c b/system/lib/libc/musl/src/internal/intscan.c
new file mode 100644
index 00000000..69350efa
--- /dev/null
+++ b/system/lib/libc/musl/src/internal/intscan.c
@@ -0,0 +1,99 @@
+#include <limits.h>
+#include <errno.h>
+#include <ctype.h>
+#include "shgetc.h"
+
+/* Lookup table for digit values. -1==255>=36 -> invalid */
+static const unsigned char table[] = { -1,
+-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
+-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
+-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
+ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,-1,-1,-1,-1,-1,-1,
+-1,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,
+25,26,27,28,29,30,31,32,33,34,35,-1,-1,-1,-1,-1,
+-1,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,
+25,26,27,28,29,30,31,32,33,34,35,-1,-1,-1,-1,-1,
+-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
+-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
+-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
+-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
+-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
+-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
+-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
+-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
+};
+
+unsigned long long __intscan(FILE *f, unsigned base, int pok, unsigned long long lim)
+{
+	const unsigned char *val = table+1;
+	int c, neg=0;
+	unsigned x;
+	unsigned long long y;
+	if (base > 36) {
+		errno = EINVAL;
+		return 0;
+	}
+	while (isspace((c=shgetc(f))));
+	if (c=='+' || c=='-') {
+		neg = -(c=='-');
+		c = shgetc(f);
+	}
+	if ((base == 0 || base == 16) && c=='0') {
+		c = shgetc(f);
+		if ((c|32)=='x') {
+			c = shgetc(f);
+			if (val[c]>=16) {
+				shunget(f);
+				if (pok) shunget(f);
+				else shlim(f, 0);
+				return 0;
+			}
+			base = 16;
+		} else if (base == 0) {
+			base = 8;
+		}
+	} else {
+		if (base == 0) base = 10;
+		if (val[c] >= base) {
+			shunget(f);
+			shlim(f, 0);
+			errno = EINVAL;
+			return 0;
+		}
+	}
+	if (base == 10) {
+		for (x=0; c-'0'<10U && x<=UINT_MAX/10-1; c=shgetc(f))
+			x = x*10 + (c-'0');
+		for (y=x; c-'0'<10U && y<=ULLONG_MAX/10 && 10*y<=ULLONG_MAX-(c-'0'); c=shgetc(f))
+			y = y*10 + (c-'0');
+		if (c-'0'>=10U) goto done;
+	} else if (!(base & base-1)) {
+		int bs = "\0\1\2\4\7\3\6\5"[(0x17*base)>>5&7];
+		for (x=0; val[c]<base && x<=UINT_MAX/32; c=shgetc(f))
+			x = x<<bs | val[c];
+		for (y=x; val[c]<base && y<=ULLONG_MAX>>bs; c=shgetc(f))
+			y = y<<bs | val[c];
+	} else {
+		for (x=0; val[c]<base && x<=UINT_MAX/36-1; c=shgetc(f))
+			x = x*base + val[c];
+		for (y=x; val[c]<base && y<=ULLONG_MAX/base && base*y<=ULLONG_MAX-val[c]; c=shgetc(f))
+			y = y*base + val[c];
+	}
+	if (val[c]<base) {
+		for (; val[c]<base; c=shgetc(f));
+		errno = ERANGE;
+		y = lim;
+	}
+done:
+	shunget(f);
+	if (y>=lim) {
+		if (!(lim&1) && !neg) {
+			errno = ERANGE;
+			return lim-1;
+		} else if (y>lim) {
+			errno = ERANGE;
+			return lim;
+		}
+	}
+	return (y^neg)-neg;
+}
diff --git a/system/lib/libc/musl/src/internal/intscan.h b/system/lib/libc/musl/src/internal/intscan.h
new file mode 100644
index 00000000..994c5e7d
--- /dev/null
+++ b/system/lib/libc/musl/src/internal/intscan.h
@@ -0,0 +1,8 @@
+#ifndef INTSCAN_H
+#define INTSCAN_H
+
+#include <stdio.h>
+
+unsigned long long __intscan(FILE *, unsigned, int, unsigned long long);
+
+#endif
diff --git a/system/lib/libc/musl/src/internal/libm.h b/system/lib/libc/musl/src/internal/libm.h
new file mode 100644
index 00000000..9f0d3bc8
--- /dev/null
+++ b/system/lib/libc/musl/src/internal/libm.h
@@ -0,0 +1,169 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/math_private.h */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#ifndef _LIBM_H
+#define _LIBM_H
+
+#include <stdint.h>
+#include <float.h>
+#include <math.h>
+#include <complex.h>
+#include <endian.h>
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN
+union ldshape {
+	long double f;
+	struct {
+		uint64_t m;
+		uint16_t se;
+	} i;
+};
+#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN
+union ldshape {
+	long double f;
+	struct {
+		uint64_t lo;
+		uint32_t mid;
+		uint16_t top;
+		uint16_t se;
+	} i;
+	struct {
+		uint64_t lo;
+		uint64_t hi;
+	} i2;
+};
+#else
+#error Unsupported long double representation
+#endif
+
+#define FORCE_EVAL(x) do {                        \
+	if (sizeof(x) == sizeof(float)) {         \
+		volatile float __x;               \
+		__x = (x);                        \
+	} else if (sizeof(x) == sizeof(double)) { \
+		volatile double __x;              \
+		__x = (x);                        \
+	} else {                                  \
+		volatile long double __x;         \
+		__x = (x);                        \
+	}                                         \
+} while(0)
+
+/* Get two 32 bit ints from a double.  */
+#define EXTRACT_WORDS(hi,lo,d)                    \
+do {                                              \
+  union {double f; uint64_t i;} __u;              \
+  __u.f = (d);                                    \
+  (hi) = __u.i >> 32;                             \
+  (lo) = (uint32_t)__u.i;                         \
+} while (0)
+
+/* Get the more significant 32 bit int from a double.  */
+#define GET_HIGH_WORD(hi,d)                       \
+do {                                              \
+  union {double f; uint64_t i;} __u;              \
+  __u.f = (d);                                    \
+  (hi) = __u.i >> 32;                             \
+} while (0)
+
+/* Get the less significant 32 bit int from a double.  */
+#define GET_LOW_WORD(lo,d)                        \
+do {                                              \
+  union {double f; uint64_t i;} __u;              \
+  __u.f = (d);                                    \
+  (lo) = (uint32_t)__u.i;                         \
+} while (0)
+
+/* Set a double from two 32 bit ints.  */
+#define INSERT_WORDS(d,hi,lo)                     \
+do {                                              \
+  union {double f; uint64_t i;} __u;              \
+  __u.i = ((uint64_t)(hi)<<32) | (uint32_t)(lo);  \
+  (d) = __u.f;                                    \
+} while (0)
+
+/* Set the more significant 32 bits of a double from an int.  */
+#define SET_HIGH_WORD(d,hi)                       \
+do {                                              \
+  union {double f; uint64_t i;} __u;              \
+  __u.f = (d);                                    \
+  __u.i &= 0xffffffff;                            \
+  __u.i |= (uint64_t)(hi) << 32;                  \
+  (d) = __u.f;                                    \
+} while (0)
+
+/* Set the less significant 32 bits of a double from an int.  */
+#define SET_LOW_WORD(d,lo)                        \
+do {                                              \
+  union {double f; uint64_t i;} __u;              \
+  __u.f = (d);                                    \
+  __u.i &= 0xffffffff00000000ull;                 \
+  __u.i |= (uint32_t)(lo);                        \
+  (d) = __u.f;                                    \
+} while (0)
+
+/* Get a 32 bit int from a float.  */
+#define GET_FLOAT_WORD(w,d)                       \
+do {                                              \
+  union {float f; uint32_t i;} __u;               \
+  __u.f = (d);                                    \
+  (w) = __u.i;                                    \
+} while (0)
+
+/* Set a float from a 32 bit int.  */
+#define SET_FLOAT_WORD(d,w)                       \
+do {                                              \
+  union {float f; uint32_t i;} __u;               \
+  __u.i = (w);                                    \
+  (d) = __u.f;                                    \
+} while (0)
+
+/* fdlibm kernel functions */
+
+int    __rem_pio2_large(double*,double*,int,int,int);
+
+int    __rem_pio2(double,double*);
+double __sin(double,double,int);
+double __cos(double,double);
+double __tan(double,double,int);
+double __expo2(double);
+double complex __ldexp_cexp(double complex,int);
+
+int    __rem_pio2f(float,double*);
+float  __sindf(double);
+float  __cosdf(double);
+float  __tandf(double,int);
+float  __expo2f(float);
+float complex __ldexp_cexpf(float complex,int);
+
+int __rem_pio2l(long double, long double *);
+long double __sinl(long double, long double, int);
+long double __cosl(long double, long double);
+long double __tanl(long double, long double, int);
+
+/* polynomial evaluation */
+long double __polevll(long double, const long double *, int);
+long double __p1evll(long double, const long double *, int);
+
+#if 0
+/* Attempt to get strict C99 semantics for assignment with non-C99 compilers. */
+#define STRICT_ASSIGN(type, lval, rval) do {    \
+        volatile type __v = (rval);             \
+        (lval) = __v;                           \
+} while (0)
+#else
+/* Should work with -fexcess-precision=standard (>=gcc-4.5) or -ffloat-store */
+#define STRICT_ASSIGN(type, lval, rval) ((lval) = (type)(rval))
+#endif
+
+#endif
diff --git a/system/lib/libc/musl/src/internal/shgetc.c b/system/lib/libc/musl/src/internal/shgetc.c
new file mode 100644
index 00000000..e878b00a
--- /dev/null
+++ b/system/lib/libc/musl/src/internal/shgetc.c
@@ -0,0 +1,27 @@
+#include "shgetc.h"
+
+void __shlim(FILE *f, off_t lim)
+{
+	f->shlim = lim;
+	f->shcnt = f->rend - f->rpos;
+	if (lim && f->shcnt > lim)
+		f->shend = f->rpos + lim;
+	else
+		f->shend = f->rend;
+}
+
+int __shgetc(FILE *f)
+{
+	int c;
+	if (f->shlim && f->shcnt >= f->shlim || (c=__uflow(f)) < 0) {
+		f->shend = 0;
+		return EOF;
+	}
+	if (f->shlim && f->rend - f->rpos > f->shlim - f->shcnt - 1)
+		f->shend = f->rpos + (f->shlim - f->shcnt - 1);
+	else
+		f->shend = f->rend;
+	if (f->rend) f->shcnt += f->rend - f->rpos + 1;
+	if (f->rpos[-1] != c) f->rpos[-1] = c;
+	return c;
+}
diff --git a/system/lib/libc/musl/src/internal/shgetc.h b/system/lib/libc/musl/src/internal/shgetc.h
new file mode 100644
index 00000000..7beb8ce6
--- /dev/null
+++ b/system/lib/libc/musl/src/internal/shgetc.h
@@ -0,0 +1,9 @@
+#include "stdio_impl.h"
+
+void __shlim(FILE *, off_t);
+int __shgetc(FILE *);
+
+#define shcnt(f) ((f)->shcnt + ((f)->rpos - (f)->rend))
+#define shlim(f, lim) __shlim((f), (lim))
+#define shgetc(f) (((f)->rpos < (f)->shend) ? *(f)->rpos++ : __shgetc(f))
+#define shunget(f) ((f)->shend ? (void)(f)->rpos-- : (void)0)
diff --git a/system/lib/libc/musl/src/math/scalbnl.c b/system/lib/libc/musl/src/math/scalbnl.c
new file mode 100644
index 00000000..08a4c587
--- /dev/null
+++ b/system/lib/libc/musl/src/math/scalbnl.c
@@ -0,0 +1,36 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double scalbnl(long double x, int n)
+{
+	return scalbn(x, n);
+}
+#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
+long double scalbnl(long double x, int n)
+{
+	union ldshape u;
+
+	if (n > 16383) {
+		x *= 0x1p16383L;
+		n -= 16383;
+		if (n > 16383) {
+			x *= 0x1p16383L;
+			n -= 16383;
+			if (n > 16383)
+				n = 16383;
+		}
+	} else if (n < -16382) {
+		x *= 0x1p-16382L;
+		n += 16382;
+		if (n < -16382) {
+			x *= 0x1p-16382L;
+			n += 16382;
+			if (n < -16382)
+				n = -16382;
+		}
+	}
+	u.f = 1.0;
+	u.i.se = 0x3fff + n;
+	return x * u.f;
+}
+#endif
diff --git a/system/lib/libc/musl/src/stdio/__overflow.c b/system/lib/libc/musl/src/stdio/__overflow.c
new file mode 100644
index 00000000..3bb37923
--- /dev/null
+++ b/system/lib/libc/musl/src/stdio/__overflow.c
@@ -0,0 +1,10 @@
+#include "stdio_impl.h"
+
+int __overflow(FILE *f, int _c)
+{
+	unsigned char c = _c;
+	if (!f->wend && __towrite(f)) return EOF;
+	if (f->wpos < f->wend && c != f->lbf) return *f->wpos++ = c;
+	if (f->write(f, &c, 1)!=1) return EOF;
+	return c;
+}
diff --git a/system/lib/libc/musl/src/stdio/__toread.c b/system/lib/libc/musl/src/stdio/__toread.c
new file mode 100644
index 00000000..2e804f64
--- /dev/null
+++ b/system/lib/libc/musl/src/stdio/__toread.c
@@ -0,0 +1,24 @@
+#include <stdio_impl.h>
+
+int __toread(FILE *f)
+{
+	f->mode |= f->mode-1;
+	if (f->wpos > f->buf) f->write(f, 0, 0);
+	f->wpos = f->wbase = f->wend = 0;
+	if (f->flags & (F_EOF|F_NORD)) {
+		if (f->flags & F_NORD) f->flags |= F_ERR;
+		return EOF;
+	}
+	f->rpos = f->rend = f->buf;
+	return 0;
+}
+
+static const int dummy = 0;
+weak_alias(dummy, __towrite_used);
+
+void __stdio_exit(void);
+
+void __seek_on_exit()
+{
+	if (!__towrite_used) __stdio_exit();
+}
diff --git a/system/lib/libc/musl/src/stdio/__uflow.c b/system/lib/libc/musl/src/stdio/__uflow.c
new file mode 100644
index 00000000..e28922c2
--- /dev/null
+++ b/system/lib/libc/musl/src/stdio/__uflow.c
@@ -0,0 +1,11 @@
+#include "stdio_impl.h"
+
+/* This function will never be called if there is already data
+ * buffered for reading. Thus we can get by with very few branches. */
+
+int __uflow(FILE *f)
+{
+	unsigned char c;
+	if ((f->rend || !__toread(f)) && f->read(f, &c, 1)==1) return c;
+	return EOF;
+}
diff --git a/system/lib/libc/musl/src/stdlib/atof.c b/system/lib/libc/musl/src/stdlib/atof.c
new file mode 100644
index 00000000..f7fcd826
--- /dev/null
+++ b/system/lib/libc/musl/src/stdlib/atof.c
@@ -0,0 +1,6 @@
+#include <stdlib.h>
+
+double atof(const char *s)
+{
+	return strtod(s, 0);
+}
diff --git a/system/lib/libc/musl/src/stdlib/strtod.c b/system/lib/libc/musl/src/stdlib/strtod.c
new file mode 100644
index 00000000..461dcf85
--- /dev/null
+++ b/system/lib/libc/musl/src/stdlib/strtod.c
@@ -0,0 +1,40 @@
+#include <stdlib.h>
+#include "shgetc.h"
+#include "floatscan.h"
+#include "stdio_impl.h"
+#include "libc.h"
+
+static long double strtox(const char *s, char **p, int prec)
+{
+	FILE f = {
+		.buf = (void *)s, .rpos = (void *)s,
+		.rend = (void *)-1, .lock = -1
+	};
+	shlim(&f, 0);
+	long double y = __floatscan(&f, prec, 1);
+	off_t cnt = shcnt(&f);
+	if (p) *p = cnt ? (char *)s + cnt : (char *)s;
+	return y;
+}
+
+float strtof(const char *restrict s, char **restrict p)
+{
+	return strtox(s, p, 0);
+}
+
+double strtod(const char *restrict s, char **restrict p)
+{
+	return strtox(s, p, 1);
+}
+
+long double strtold(const char *restrict s, char **restrict p)
+{
+	return strtox(s, p, 2);
+}
+
+weak_alias(strtof, strtof_l);
+weak_alias(strtod, strtod_l);
+weak_alias(strtold, strtold_l);
+weak_alias(strtof, __strtof_l);
+weak_alias(strtod, __strtod_l);
+weak_alias(strtold, __strtold_l);
diff --git a/system/lib/libc/musl/src/stdlib/wcstod.c b/system/lib/libc/musl/src/stdlib/wcstod.c
new file mode 100644
index 00000000..83f308d3
--- /dev/null
+++ b/system/lib/libc/musl/src/stdlib/wcstod.c
@@ -0,0 +1,64 @@
+#include "shgetc.h"
+#include "floatscan.h"
+#include "stdio_impl.h"
+#include <wctype.h>
+
+/* This read function heavily cheats. It knows:
+ *  (1) len will always be 1
+ *  (2) non-ascii characters don't matter */
+
+static size_t do_read(FILE *f, unsigned char *buf, size_t len)
+{
+	size_t i;
+	const wchar_t *wcs = f->cookie;
+
+	if (!wcs[0]) wcs=L"@";
+	for (i=0; i<f->buf_size && wcs[i]; i++)
+		f->buf[i] = wcs[i] < 128 ? wcs[i] : '@';
+	f->rpos = f->buf;
+	f->rend = f->buf + i;
+	f->cookie = (void *)(wcs+i);
+
+	if (i && len) {
+		*buf = *f->rpos++;
+		return 1;
+	}
+	return 0;
+}
+
+static long double wcstox(const wchar_t *s, wchar_t **p, int prec)
+{
+	wchar_t *t = (wchar_t *)s;
+	unsigned char buf[64];
+	FILE f = {0};
+	f.flags = 0;
+	f.rpos = f.rend = 0;
+	f.buf = buf + 4;
+	f.buf_size = sizeof buf - 4;
+	f.lock = -1;
+	f.read = do_read;
+	while (iswspace(*t)) t++;
+	f.cookie = (void *)t;
+	shlim(&f, 0);
+	long double y = __floatscan(&f, prec, 1);
+	if (p) {
+		size_t cnt = shcnt(&f);
+		*p = cnt ? t + cnt : (wchar_t *)s;
+	}
+	return y;
+}
+
+float wcstof(const wchar_t *restrict s, wchar_t **restrict p)
+{
+	return wcstox(s, p, 0);
+}
+
+double wcstod(const wchar_t *restrict s, wchar_t **restrict p)
+{
+	return wcstox(s, p, 1);
+}
+
+long double wcstold(const wchar_t *restrict s, wchar_t **restrict p)
+{
+	return wcstox(s, p, 2);
+}
diff --git a/system/lib/libc/musl/src/stdlib/wcstol.c b/system/lib/libc/musl/src/stdlib/wcstol.c
new file mode 100644
index 00000000..4443f577
--- /dev/null
+++ b/system/lib/libc/musl/src/stdlib/wcstol.c
@@ -0,0 +1,82 @@
+#include "stdio_impl.h"
+#include "intscan.h"
+#include "shgetc.h"
+#include <inttypes.h>
+#include <limits.h>
+#include <wctype.h>
+#include <wchar.h>
+
+/* This read function heavily cheats. It knows:
+ *  (1) len will always be 1
+ *  (2) non-ascii characters don't matter */
+
+static size_t do_read(FILE *f, unsigned char *buf, size_t len)
+{
+	size_t i;
+	const wchar_t *wcs = f->cookie;
+
+	if (!wcs[0]) wcs=L"@";
+	for (i=0; i<f->buf_size && wcs[i]; i++)
+		f->buf[i] = wcs[i] < 128 ? wcs[i] : '@';
+	f->rpos = f->buf;
+	f->rend = f->buf + i;
+	f->cookie = (void *)(wcs+i);
+
+	if (i && len) {
+		*buf = *f->rpos++;
+		return 1;
+	}
+	return 0;
+}
+
+static unsigned long long wcstox(const wchar_t *s, wchar_t **p, int base, unsigned long long lim)
+{
+	wchar_t *t = (wchar_t *)s;
+	unsigned char buf[64];
+	FILE f = {0};
+	f.flags = 0;
+	f.rpos = f.rend = 0;
+	f.buf = buf + 4;
+	f.buf_size = sizeof buf - 4;
+	f.lock = -1;
+	f.read = do_read;
+	while (iswspace(*t)) t++;
+	f.cookie = (void *)t;
+	shlim(&f, 0);
+	unsigned long long y = __intscan(&f, base, 1, lim);
+	if (p) {
+		size_t cnt = shcnt(&f);
+		*p = cnt ? t + cnt : (wchar_t *)s;
+	}
+	return y;
+}
+
+unsigned long long wcstoull(const wchar_t *restrict s, wchar_t **restrict p, int base)
+{
+	return wcstox(s, p, base, ULLONG_MAX);
+}
+
+long long wcstoll(const wchar_t *restrict s, wchar_t **restrict p, int base)
+{
+	return wcstox(s, p, base, LLONG_MIN);
+}
+
+unsigned long wcstoul(const wchar_t *restrict s, wchar_t **restrict p, int base)
+{
+	return wcstox(s, p, base, ULONG_MAX);
+}
+
+long wcstol(const wchar_t *restrict s, wchar_t **restrict p, int base)
+{
+	return wcstox(s, p, base, 0UL+LONG_MIN);
+}
+
+intmax_t wcstoimax(const wchar_t *restrict s, wchar_t **restrict p, int base)
+{
+	return wcstoll(s, p, base);
+}
+
+uintmax_t wcstoumax(const wchar_t *restrict s, wchar_t **restrict p, int base)
+{
+	return wcstoull(s, p, base);
+}
diff --git a/system/lib/libcextra.symbols b/system/lib/libcextra.symbols
index eb8526de..a32bdbdb 100644
--- a/system/lib/libcextra.symbols
+++ b/system/lib/libcextra.symbols
@@ -1,3 +1,4 @@
+         T __intscan
          W __iswctype_l
          T __memrchr
          T __strchrnul
@@ -108,8 +109,17 @@
          T wcsrtombs
          T wcsspn
          T wcsstr
+         T wcstod
+         T wcstof
+         T wcstoimax
          T wcstok
+         T wcstol
+         T wcstold
+         T wcstoll
          T wcstombs
+         T wcstoul
+         T wcstoull
+         T wcstoumax
          T wcswcs
          T wcswidth
          T wcsxfrm