Browse Source

JSON: Add musl functions for float parsing without standard libc

Lukas M 5 years ago
parent
commit
af9c7c2522
4 changed files with 1320 additions and 0 deletions
  1. 719 0
      deps/musl/floatscan.c
  2. 18 0
      deps/musl/floatscan.h
  3. 514 0
      deps/musl/vfprintf.c
  4. 69 0
      deps/musl/vfprintf.h

+ 719 - 0
deps/musl/floatscan.c

@@ -0,0 +1,719 @@
+/* Originally released by the musl project (http://www.musl-libc.org/) under the
+ * MIT license. Taken from the file src/internal/floatscan.c*/
+
+#include <stdint.h>
+#include <math.h>
+#include <limits.h>
+#include <errno.h>
+#include <ctype.h>
+
+#include "floatscan.h"
+#include "vfprintf.h"
+int shgetc(char* input, int *index);
+void shunget(int *index);
+int shlim(int a, int b);
+
+int shgetc(char* input, int *index){
+    int res = input[*index];
+    (*index)++;
+    return res;
+}
+
+void shunget(int *index){
+    (*index)--;
+}
+int shlim(int a, int b){
+    return '0';
+}
+
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double fmodl(long double x, long double y)
+{
+	return fmod(x, y);
+}
+#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
+long double fmodl(long double x, long double y)
+{
+	union ldshape ux = {x}, uy = {y};
+	int ex = ux.i.se & 0x7fff;
+	int ey = uy.i.se & 0x7fff;
+	int sx = ux.i.se & 0x8000;
+
+	if (y == 0 || isnan(y) || ex == 0x7fff)
+		return (x*y)/(x*y);
+	ux.i.se = (uint16_t)ex;
+	uy.i.se = (uint16_t)ey;
+	if (ux.f <= uy.f) {
+		if (ux.f == uy.f)
+			return 0*x;
+		return x;
+	}
+
+	/* normalize x and y */
+	if (!ex) {
+		ux.f *= 0x1p120f;
+		ex = ux.i.se - 120;
+	}
+	if (!ey) {
+		uy.f *= 0x1p120f;
+		ey = uy.i.se - 120;
+	}
+
+	/* x mod y */
+#if LDBL_MANT_DIG == 64
+	uint64_t i, mx, my;
+	mx = ux.i.m;
+	my = uy.i.m;
+	for (; ex > ey; ex--) {
+		i = mx - my;
+		if (mx >= my) {
+			if (i == 0)
+				return 0*x;
+			mx = 2*i;
+		} else if (2*mx < mx) {
+			mx = 2*mx - my;
+		} else {
+			mx = 2*mx;
+		}
+	}
+	i = mx - my;
+	if (mx >= my) {
+		if (i == 0)
+			return 0*x;
+		mx = i;
+	}
+	for (; mx >> 63 == 0; mx *= 2, ex--);
+	ux.i.m = mx;
+#elif LDBL_MANT_DIG == 113
+	uint64_t hi, lo, xhi, xlo, yhi, ylo;
+	xhi = (ux.i2.hi & -1ULL>>16) | 1ULL<<48;
+	yhi = (uy.i2.hi & -1ULL>>16) | 1ULL<<48;
+	xlo = ux.i2.lo;
+	ylo = uy.i2.lo;
+	for (; ex > ey; ex--) {
+		hi = xhi - yhi;
+		lo = xlo - ylo;
+		if (xlo < ylo)
+			hi -= 1;
+		if (hi >> 63 == 0) {
+			if ((hi|lo) == 0)
+				return 0*x;
+			xhi = 2*hi + (lo>>63);
+			xlo = 2*lo;
+		} else {
+			xhi = 2*xhi + (xlo>>63);
+			xlo = 2*xlo;
+		}
+	}
+	hi = xhi - yhi;
+	lo = xlo - ylo;
+	if (xlo < ylo)
+		hi -= 1;
+	if (hi >> 63 == 0) {
+		if ((hi|lo) == 0)
+			return 0*x;
+		xhi = hi;
+		xlo = lo;
+	}
+	for (; xhi >> 48 == 0; xhi = 2*xhi + (xlo>>63), xlo = 2*xlo, ex--);
+	ux.i2.hi = xhi;
+	ux.i2.lo = xlo;
+#endif
+
+	/* scale result */
+	if (ex <= 0) {
+		ux.i.se = (uint16_t)((ex+120)|sx);
+		ux.f *= (uint16_t)(0x1p-120f);
+	} else
+		ux.i.se = (uint16_t)(ex|sx);
+	return ux.f;
+}
+#endif
+
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double copysignl(long double x, long double y)
+{
+	return copysign(x, y);
+}
+#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
+long double copysignl(long double x, long double y)
+{
+	union ldshape ux = {x}, uy = {y};
+	ux.i.se &= 0x7fff;
+	ux.i.se = (uint16_t)(ux.i.se | (uy.i.se & 0x8000));
+	return ux.f;
+}
+#endif
+
+double scalbn(double x, int n)
+{
+	union {double f; uint64_t i;} u;
+	double_t y = x;
+
+	if (n > 1023) {
+		y *= 0x1p1023;
+		n -= 1023;
+		if (n > 1023) {
+			y *= 0x1p1023;
+			n -= 1023;
+			if (n > 1023)
+				n = 1023;
+		}
+	} else if (n < -1022) {
+		/* make sure final n < -53 to avoid double
+		   rounding in the subnormal range */
+		y *= 0x1p-1022 * 0x1p53;
+		n += 1022 - 53;
+		if (n < -1022) {
+			y *= 0x1p-1022 * 0x1p53;
+			n += 1022 - 53;
+			if (n < -1022)
+				n = -1022;
+		}
+	}
+	u.i = (uint64_t)(0x3ff+n)<<52;
+	x = y * u.f;
+	return x;
+}
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+
+#define LD_B1B_DIG 2
+#define LD_B1B_MAX 9007199, 254740991
+#define KMAX 128
+
+#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
+
+#define LD_B1B_DIG 3
+#define LD_B1B_MAX 18, 446744073, 709551615
+#define KMAX 2048
+
+#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
+
+#define LD_B1B_DIG 4
+#define LD_B1B_MAX 10384593, 717069655, 257060992, 658440191
+#define KMAX 2048
+
+#else
+#error Unsupported long double representation
+#endif
+
+#define MASK (KMAX-1)
+
+#define CONCAT2(x,y) x ## y
+#define CONCAT(x,y) CONCAT2(x,y)
+
+
+#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 * 0x1p113L;
+		n += 16382 - 113;
+		if (n < -16382) {
+			x *= 0x1p-16382L * 0x1p113L;
+			n += 16382 - 113;
+			if (n < -16382)
+				n = -16382;
+		}
+	}
+	u.f = 1.0;
+	u.i.se = (uint16_t)(0x3fff + n);
+	return x * u.f;
+}
+#endif
+
+static long long scanexp(char* input, int *index, int pok)
+{
+	int c;
+	int x;
+	long long y;
+	int neg = 0;
+	
+	c = shgetc(input, index);
+	if (c=='+' || c=='-') {
+		neg = (c=='-');
+		c = shgetc(input, index);
+		if ((unsigned)(c-'0')>=10U && pok) shunget(index);
+	}
+	if ((unsigned)(c-'0')>=10U) {
+		shunget(index);
+		return LLONG_MIN;
+	}
+	for (x=0; (unsigned)(c-'0')<10U && x<INT_MAX/10; c = shgetc(input, index))
+		x = 10*x + c-'0';
+	for (y=x; (unsigned)(c-'0')<10U && y<LLONG_MAX/100; c = shgetc(input, index))
+		y = 10*y + c-'0';
+	for (; (unsigned)(c-'0')<10U; c = shgetc(input, index));
+	shunget(index);
+	return neg ? -y : y;
+}
+
+
+static long double decfloat(char *input, int *index, 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(input, index)) gotdig=1;
+	if (c=='.') {
+		gotrad = 1;
+		for (c = shgetc(input, index); c=='0'; c = shgetc(input, index)) gotdig=1, lrp--;
+	}
+
+	x[0] = 0;
+	for (; (unsigned)(c-'0')<10U || c=='.'; c = shgetc(input, index)) {
+		if (c == '.') {
+			if (gotrad) break;
+			gotrad = 1;
+			lrp = dc;
+		} else if (k < KMAX-3) {
+			dc++;
+			if (c!='0') lnz = (int)dc;
+			if (j) x[k] = (x[k]*10 + (uint32_t)(c-'0'));
+			else x[k] = (uint32_t)(c-'0');
+			if (++j==9) {
+				k++;
+				j=0;
+			}
+			gotdig=1;
+		} else {
+			dc++;
+			if (c!='0') {
+				lnz = (KMAX-4)*9;
+				x[KMAX-4] |= 1;
+			}
+		}
+	}
+	if (!gotrad) lrp=dc;
+
+	if (gotdig && (c|32)=='e') {
+		e10 = scanexp(input, index, pok);
+		if (e10 == LLONG_MIN) {
+			if (pok) {
+				shunget(index);
+			} else {
+				//shlim(f, 0);
+				return 0;
+			}
+			e10 = 0;
+		}
+		lrp += e10;
+	} else if (c>=0) {
+		shunget(index);
+	}
+	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 = (int)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];
+	}
+
+	/* Drop trailing zeros */
+	for (; !x[z-1]; z--);
+
+	/* 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] % (uint32_t)p10);
+			x[k] = x[k]/(uint32_t)p10 + carry;
+			carry = 1000000000/(uint32_t)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 = (uint32_t)(tmp / 1000000000);
+				x[k] = (uint32_t)(tmp % 1000000000);
+			} else {
+				carry = 0;
+				x[k] = (uint32_t)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] & (uint32_t)((1<<sh)-1));
+			x[k] = (x[k]>>sh) + carry;
+			carry = ((uint32_t)(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)) implicit conversion turns floating-point number into integer:
+                if (LDBL_MANT_DIG-bits >= 2 && !((_Bool)fmodl(frac, 1)))
+			frac++;
+	}
+
+	y += frac;
+	y -= bias;
+
+	if (((e2+LDBL_MANT_DIG) & INT_MAX) > emax-5) {
+		if (fabs((double)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)) implicit conversion turns floating-point number into integer:
+		if (e2+LDBL_MANT_DIG>emax || ((_Bool)denormal && (_Bool)frac))
+			errno = ERANGE;
+	}
+
+	return scalbnl(y, e2);
+}
+
+static long double hexfloat(char *input, int *index, 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(input, index);
+
+	/* Skip leading zeros */
+	for (; c=='0'; c = shgetc(input, index)) gotdig = 1;
+
+	if (c=='.') {
+		gotrad = 1;
+		c = shgetc(input, index);
+		/* Count zeros after the radix point before significand */
+		for (rp=0; c=='0'; c = shgetc(input, index), rp--) gotdig = 1;
+	}
+
+	for (; (unsigned)(c-'0')<10U || (unsigned)((c|32)-'a')<6U || c=='.'; c = shgetc(input, index)) {
+		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 + (uint32_t)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(index);
+		if (pok) {
+			shunget(index);
+			if (gotrad) shunget(index);
+		} else {
+			//shlim(f, 0);
+		}
+		return sign * 0.0;
+	}
+	if (!gotrad) rp = dc;
+	while (dc<8) x *= 16, dc++;
+	if ((c|32)=='p') {
+		e2 = scanexp(input, index, pok);
+		if (e2 == LLONG_MIN) {
+			if (pok) {
+				shunget(index);
+			} else {
+				//shlim(f, 0);
+				return 0;
+			}
+			e2 = 0;
+		}
+	} else {
+		shunget(index);
+	}
+	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 =(int)(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 && (_Bool)y && !(x&1)) x++, y=0;
+
+	y = bias + sign*(long double)x + sign*y;
+	y -= bias;
+
+	if (!((_Bool)y)) errno = ERANGE;
+
+	return scalbnl(y, (int)e2);
+}
+
+long double __floatscan(char* input, int prec, int pok)
+{
+    int index = 0;
+	int sign = 1;
+	//size_t i;
+        int 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(input, &index))));
+
+	if (c=='+' || c=='-') {
+		sign -= 2*(c=='-');
+		c = shgetc(input, &index);
+	}
+
+	for (i=0; i<8 && (c|32)=="infinity"[i]; i++)
+		if (i<7) c = shgetc(input, &index);
+	if (i==3 || i==8 || (i>3 && pok)) {
+		if (i!=8) {
+			shunget(&index);
+			if (pok) for (; i>3; i--) shunget(&index);
+		}
+		return ((float)sign * INFINITY);
+	}
+	if (!i) for (i=0; i<3 && (c|32)=="nan"[i]; i++)
+		if (i<2) c = shgetc(input, &index);
+	if (i==3) {
+		if (shgetc(input, &index) != '(') {
+			shunget(&index);
+			return NAN;
+		}
+		for (i=1; ; i++) {
+			c = shgetc(input, &index);
+			if ((unsigned)(c-'0')<10U || (unsigned)(c-'A')<26U || (unsigned)(c-'a')<26U || c=='_')
+				continue;
+			if (c==')') return NAN;
+			shunget(&index);
+			if (!pok) {
+				errno = EINVAL;
+				//shlim(0, 0);
+				return 0;
+			}
+			while (i--) shunget(&index);
+			return NAN;
+		}
+		return NAN;
+	}
+
+	if (i) {
+		shunget(&index);
+		errno = EINVAL;
+		//shlim(0, 0);
+		return 0;
+	}
+
+	if (c=='0') {
+		c = shgetc(input, &index);
+		if ((c|32) == 'x')
+			return hexfloat(input, &index, bits, emin, sign, pok);
+		shunget(&index);
+		c = '0';
+	}
+
+	return decfloat(input, &index, c, bits, emin, sign, pok);
+}

+ 18 - 0
deps/musl/floatscan.h

@@ -0,0 +1,18 @@
+#ifndef FLOATSCAN_H
+#define FLOATSCAN_H
+
+//#include <stdio.h>
+
+#define FLT_MANT_DIG 24
+#define FLT_MIN_EXP (-125)
+#define FLT_MAX_EXP 128
+#define FLT_HAS_SUBNORM 1
+
+#define DBL_MANT_DIG 53
+#define DBL_MIN_EXP (-1021)
+#define DBL_MAX_EXP 1024
+#define DBL_HAS_SUBNORM 1
+
+long double __floatscan(char*, int, int);
+
+#endif

+ 514 - 0
deps/musl/vfprintf.c

@@ -0,0 +1,514 @@
+/* Originally released by the musl project (http://www.musl-libc.org/) under the
+ * MIT license. Taken from the file src/stdio/vfprintf.c */
+
+#include <ctype.h>
+#include "floatscan.h"
+#include "vfprintf.h"
+//int isdigit(int);
+//#define isdigit(a) (0 ? isdigit(a) : ((unsigned)(a)-'0') < 10)
+
+long double frexpl(long double x, int *e);
+//#include <math.h>
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double frexpl(long double x, int *e)
+{
+	return frexp(x, e);
+}
+#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
+long double frexpl(long double x, int *e)
+{
+	union ldshape u = {x};
+	int ee = u.i.se & 0x7fff;
+
+	if (!ee) {
+		if ((_Bool)x) {
+			x = frexpl(x*0x1p120, e);
+			*e -= 120;
+		} else *e = 0;
+		return x;
+	} else if (ee == 0x7fff) {
+		return x;
+	}
+
+	*e = ee - 0x3ffe;
+	u.i.se &= 0x8000;
+	u.i.se |= 0x3ffe;
+	return u.f;
+}
+#endif
+
+int __signbitl(long double x);
+#if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
+int __signbitl(long double x)
+{
+	union ldshape u = {x};
+	return u.i.se >> 15;
+}
+#elif LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+int __signbitl(long double x)
+{
+	return __signbit(x);
+}
+#endif
+
+static __inline unsigned __FLOAT_BITS(float __f)
+{
+	union {float __f; unsigned __i;} __u;
+	__u.__f = __f;
+	return __u.__i;
+}
+static __inline unsigned long long __DOUBLE_BITS(double __f)
+{
+	union {double __f; unsigned long long __i;} __u;
+	__u.__f = __f;
+	return __u.__i;
+}
+
+#define signbit(x) ( \
+	sizeof(x) == sizeof(float) ? (int)(__FLOAT_BITS((float)x)>>31) : \
+	sizeof(x) == sizeof(double) ? (int)(__DOUBLE_BITS((double)x)>>63) : \
+	__signbitl(x) )
+
+#define FP_NAN       0
+#define FP_INFINITE  1
+#define FP_ZERO      2
+#define FP_SUBNORMAL 3
+#define FP_NORMAL    4
+
+
+int __fpclassifyl(long double x);
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+int __fpclassifyl(long double x)
+{
+	return __fpclassify(x);
+}
+#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
+int __fpclassifyl(long double x)
+{
+	union ldshape u = {x};
+	int e = u.i.se & 0x7fff;
+	int msb = (int)(u.i.m>>63);
+	if (!e && !msb)
+		return u.i.m ? FP_SUBNORMAL : FP_ZERO;
+	if (!msb)
+		return FP_NAN;
+	if (e == 0x7fff)
+		return u.i.m << 1 ? FP_NAN : FP_INFINITE;
+	return FP_NORMAL;
+}
+#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
+int __fpclassifyl(long double x)
+{
+	union ldshape u = {x};
+	int e = u.i.se & 0x7fff;
+	u.i.se = 0;
+	if (!e)
+		return u.i2.lo | u.i2.hi ? FP_SUBNORMAL : FP_ZERO;
+	if (e == 0x7fff)
+		return u.i2.lo | u.i2.hi ? FP_NAN : FP_INFINITE;
+	return FP_NORMAL;
+}
+#endif
+
+#define isfinite(x) ( \
+	sizeof(x) == sizeof(float) ? (__FLOAT_BITS((float)x) & 0x7fffffff) < 0x7f800000 : \
+	sizeof(x) == sizeof(double) ? (__DOUBLE_BITS((double)x) & -1ULL>>1) < 0x7ffULL<<52 : \
+	__fpclassifyl(x) > FP_INFINITE)
+
+#include "vfprintf.h"
+/* Some useful macros */
+
+#define MAX(a,b) ((a)>(b) ? (a) : (b))
+#define MIN(a,b) ((a)<(b) ? (a) : (b))
+
+/* Convenient bit representation for modifier flags, which all fall
+ * within 31 codepoints of the space character. */
+
+#define ALT_FORM   (1U<<('#'-' '))
+#define ZERO_PAD   (1U<<('0'-' '))
+#define LEFT_ADJ   (1U<<('-'-' '))
+#define PAD_POS    (1U<<(' '-' '))
+#define MARK_POS   (1U<<('+'-' '))
+#define GROUPED    (1U<<('\''-' '))
+
+#define FLAGMASK (ALT_FORM|ZERO_PAD|LEFT_ADJ|PAD_POS|MARK_POS|GROUPED)
+
+/* State machine to accept length modifiers + conversion specifiers.
+ * Result is 0 on failure, or an argument type to pop on success. */
+
+enum {
+	BARE, LPRE, LLPRE, HPRE, HHPRE, BIGLPRE,
+	ZTPRE, JPRE,
+	STOP,
+	PTR, INT, UINT, ULLONG,
+	LONG, ULONG,
+	SHORT, USHORT, CHAR, UCHAR,
+	LLONG, SIZET, IMAX, UMAX, PDIFF, UIPTR,
+	DBL, LDBL,
+	NOARG,
+	MAXSTATE
+};
+
+#define S(x) [(x)-'A']
+/*
+static const unsigned char states[]['z'-'A'+1] = {
+	{ // 0: bare types 
+		S('d') = INT, S('i') = INT,
+		S('o') = UINT, S('u') = UINT, S('x') = UINT, S('X') = UINT,
+		S('e') = DBL, S('f') = DBL, S('g') = DBL, S('a') = DBL,
+		S('E') = DBL, S('F') = DBL, S('G') = DBL, S('A') = DBL,
+		S('c') = CHAR, S('C') = INT,
+		S('s') = PTR, S('S') = PTR, S('p') = UIPTR, S('n') = PTR,
+		S('m') = NOARG,
+		S('l') = LPRE, S('h') = HPRE, S('L') = BIGLPRE,
+		S('z') = ZTPRE, S('j') = JPRE, S('t') = ZTPRE,
+	}, { // 1: l-prefixed 
+		S('d') = LONG, S('i') = LONG,
+		S('o') = ULONG, S('u') = ULONG, S('x') = ULONG, S('X') = ULONG,
+		S('e') = DBL, S('f') = DBL, S('g') = DBL, S('a') = DBL,
+		S('E') = DBL, S('F') = DBL, S('G') = DBL, S('A') = DBL,
+		S('c') = INT, S('s') = PTR, S('n') = PTR,
+		S('l') = LLPRE,
+	}, { // 2: ll-prefixed 
+		S('d') = LLONG, S('i') = LLONG,
+		S('o') = ULLONG, S('u') = ULLONG,
+		S('x') = ULLONG, S('X') = ULLONG,
+		S('n') = PTR,
+	}, { //3: h-prefixed 
+		S('d') = SHORT, S('i') = SHORT,
+		S('o') = USHORT, S('u') = USHORT,
+		S('x') = USHORT, S('X') = USHORT,
+		S('n') = PTR,
+		S('h') = HHPRE,
+	}, { // 4: hh-prefixed 
+		S('d') = CHAR, S('i') = CHAR,
+		S('o') = UCHAR, S('u') = UCHAR,
+		S('x') = UCHAR, S('X') = UCHAR,
+		S('n') = PTR,
+	}, { // 5: L-prefixed 
+		S('e') = LDBL, S('f') = LDBL, S('g') = LDBL, S('a') = LDBL,
+		S('E') = LDBL, S('F') = LDBL, S('G') = LDBL, S('A') = LDBL,
+		S('n') = PTR,
+	}, { // 6: z- or t-prefixed (assumed to be same size) 
+		S('d') = PDIFF, S('i') = PDIFF,
+		S('o') = SIZET, S('u') = SIZET,
+		S('x') = SIZET, S('X') = SIZET,
+		S('n') = PTR,
+	}, { //7: j-prefixed 
+		S('d') = IMAX, S('i') = IMAX,
+		S('o') = UMAX, S('u') = UMAX,
+		S('x') = UMAX, S('X') = UMAX,
+		S('n') = PTR,
+	}
+};
+*/
+#define OOB(x) ((unsigned)(x)-'A' > 'z'-'A')
+
+
+
+static void out(char **sp, const char *s, size_t l)
+{
+	//if (!(f->flags & F_ERR)) __fwritex((void *)s, l, f);
+        while (l--) {
+                **sp = *s;
+                (*sp)++;
+                s++;
+        }
+}
+
+static void pad(char **sp, char c, int w, int l, int fl)
+{
+	char pad[256];
+	if ((unsigned int)fl & (LEFT_ADJ | ZERO_PAD) || l >= w) return;
+	l = w - l;
+	memset(pad, c, (long unsigned int)l>sizeof pad ? sizeof pad : (long unsigned int)l);
+	for (; (long unsigned int)l >= sizeof pad; l = l - (int)(sizeof pad))
+		out(sp, pad, sizeof pad);
+	out(sp, pad, (size_t)l);
+}
+
+static const char xdigits[17] = {"0123456789ABCDEF"};
+
+/*
+static char *fmt_x(uintmax_t x, char *s, int lower)
+{
+	for (; x; x>>=4) *--s = xdigits[(x&15)]|lower;
+	return s;
+}
+
+static char *fmt_o(uintmax_t x, char *s)
+{
+	for (; x; x>>=3) *--s = '0' + (x&7);
+	return s;
+}
+*/
+static char *fmt_u(uintmax_t x, char *s)
+{
+	unsigned long y;
+	for (   ; x>ULONG_MAX; x/=10) *--s = (char)('0' + x%10);
+	for (y=x;           y; y/=10) *--s = (char)('0' + y%10);
+	return s;
+}
+
+/* Do not override this check. The floating point printing code below
+ * depends on the float.h constants being right. If they are wrong, it
+ * may overflow the stack. */
+#if LDBL_MANT_DIG == 53
+typedef char compiler_defines_long_double_incorrectly[9-(int)sizeof(long double)];
+#endif
+
+
+/*
+ * w = string width, p = precision, fl = flags, t = type. "%20.5g gives w=20, p=5, fl=0, t='g'
+ */
+int fmt_fp(char *output, long double y, int w, int p, int fl, int t)
+{
+    char* sp = output;
+	uint32_t big[(LDBL_MANT_DIG+28)/29 + 1          // mantissa expansion
+		+ (LDBL_MAX_EXP+LDBL_MANT_DIG+28+8)/9]; // exponent expansion
+	uint32_t *a, *d, *r, *z;
+	int e2=0, e, i, j, l;
+	char buf[9+LDBL_MANT_DIG/4], *s;
+	const char *prefix="-0X+0X 0X-0x+0x 0x";
+	int pl;
+	char ebuf0[3*sizeof(int)], *ebuf=&ebuf0[3*sizeof(int)], *estr = NULL;
+
+	pl=1;
+	if (signbit(y)) {
+		y=-y;
+	} else if ((unsigned int)fl & MARK_POS) {
+		prefix+=3;
+	} else if ((unsigned int)fl & PAD_POS) {
+		prefix+=6;
+	} else prefix++, pl=0;
+
+	if (!isfinite(y)) {
+		s = (t&32)?"inf":"INF";
+		if (y!=y) s=(t&32)?"nan":"NAN";
+		pad(&sp, ' ', w, 3+pl, (int)((unsigned int)fl &~ ZERO_PAD));
+		out(&sp, prefix, (size_t)pl);
+		out(&sp, s, 3);
+		pad(&sp, ' ', w, 3+pl, (int)((unsigned int)fl^LEFT_ADJ));
+		return MAX(w, 3+pl);
+	}
+
+	y = frexpl(y, &e2) * 2;
+	if ((_Bool)y) e2--;
+
+	if ((t|32)=='a') {
+		long double round = 8.0;
+		int re;
+
+		if (t&32) prefix += 9;
+		pl += 2;
+
+		if (p<0 || p>=LDBL_MANT_DIG/4-1) re=0;
+		else re=LDBL_MANT_DIG/4-1-p;
+
+		if (re) {
+			while (re--) round*=16;
+			if (*prefix=='-') {
+				y=-y;
+				y-=round;
+				y+=round;
+				y=-y;
+			} else {
+				y+=round;
+				y-=round;
+			}
+		}
+
+		estr=fmt_u((uintmax_t )(e2<0 ? -e2 : e2), ebuf);
+		if (estr==ebuf) *--estr='0';
+		*--estr = (e2<0 ? '-' : '+');
+		*--estr = (char)(t+('p'-'a'));
+
+		s=buf;
+		do {
+			int x=(int)y;
+			*s++=(char)(xdigits[x]|(t&32));
+			y=16*(y-x);
+			if (s-buf==1 && ((_Bool)y||p>0||((unsigned int)fl&ALT_FORM))) *s++='.';
+		} while ((_Bool)y);
+
+		if (p > INT_MAX-2-(ebuf-estr)-pl)
+			return -1;
+		if (p && s-buf-2 < p)
+			l = (int)((p+2) + (ebuf-estr));
+		else
+			l = (int)((s-buf) + (ebuf-estr));
+
+		pad(&sp, ' ', w, pl+l, fl);
+		out(&sp, prefix, (size_t)pl);
+		pad(&sp, '0', w, pl+l, (int)((unsigned int)fl^ZERO_PAD));
+		out(&sp, buf, (size_t)(s-buf));
+		pad(&sp, '0', (int)(l-(ebuf-estr)-(s-buf)), 0, 0);
+		out(&sp, estr, (size_t)(ebuf-estr));
+		pad(&sp, ' ', w, pl+l, (int)((unsigned int)fl^LEFT_ADJ));
+		return MAX(w, pl+l);
+	}
+	if (p<0) p=6;
+
+	if ((_Bool)y) y *= 0x1p28, e2-=28;
+
+	if (e2<0) a=r=z=big;
+	else a=r=z=big+sizeof(big)/sizeof(*big) - LDBL_MANT_DIG - 1;
+
+	do {
+		*z = (uint32_t)y;
+		y = 1000000000*(y-*z++);
+	} while ((_Bool)y);
+
+	while (e2>0) {
+		uint32_t carry=0;
+		int sh=MIN(29,e2);
+		for (d=z-1; d>=a; d--) {
+			uint64_t x = ((uint64_t)*d<<sh)+carry;
+			*d = (uint32_t)(x % 1000000000);
+			carry = (uint32_t)(x / 1000000000);
+		}
+		if (carry) *--a = carry;
+		while (z>a && !z[-1]) z--;
+		e2-=sh;
+	}
+	while (e2<0) {
+		uint32_t carry=0, *b;
+		int sh=MIN(9,-e2), need=(int)(1+((unsigned int)p+LDBL_MANT_DIG/3U+8)/9);
+		for (d=a; d<z; d++) {
+			uint32_t rm = (*d & (uint32_t)((1<<sh)-1));
+			*d = (*d>>sh) + carry;
+			carry = ((uint32_t)(1000000000>>sh) * rm);
+		}
+		if (!*a) a++;
+		if (carry) *z++ = carry;
+		/* Avoid (slow!) computation past requested precision */
+		b = (t|32)=='f' ? r : a;
+		if (z-b > need) z = b+need;
+		e2+=sh;
+	}
+
+	if (a<z) for (i=10, e=(int)(9*(r-a)); *a>=(unsigned)i; i*=10, e++);
+	else e=0;
+
+	/* Perform rounding: j is precision after the radix (possibly neg) */
+	j = p - ((t|32)!='f')*e - ((t|32)=='g' && p);
+	if (j < 9*(z-r-1)) {
+		uint32_t x;
+		/* We avoid C's broken division of negative numbers */
+		d = r + 1 + ((j+9*LDBL_MAX_EXP)/9 - LDBL_MAX_EXP);
+		j += 9*LDBL_MAX_EXP;
+		j %= 9;
+		for (i=10, j++; j<9; i*=10, j++);
+		x = (*d % (uint32_t)i);
+		/* Are there any significant digits past j? */
+		if (x || d+1!=z) {
+			long double round = 2/LDBL_EPSILON;
+			long double small;
+			if ((*d/(uint32_t)(i) & 1) || (i==1000000000 && d>a && (d[-1]&1)))
+				round += 2;
+			if (x<(unsigned)i/2) small=0x0.8p0;
+			else if (x==(unsigned)i/2 && d+1==z) small=0x1.0p0;
+			else small=0x1.8p0;
+			if (pl && *prefix=='-') round*=-1, small*=-1;
+			*d -= x;
+			/* Decide whether to round by probing round+small */
+			if (round+small != round) {
+				*d = *d + (uint32_t)i;
+				while (*d > 999999999) {
+					*d--=0;
+					if (d<a) *--a=0;
+					(*d)++;
+				}
+				for (i=10, e=(int)(9*(r-a)); *a>=(unsigned)i; i*=10, e++);
+			}
+		}
+		if (z>d+1) z=d+1;
+	}
+	for (; z>a && !z[-1]; z--);
+	
+	if ((t|32)=='g') {
+		if (!p) p++;
+		if (p>e && e>=-4) {
+			t--;
+			p-=e+1;
+		} else {
+			t-=2;
+			p--;
+		}
+		if (!((uint32_t)fl&ALT_FORM)) {
+			/* Count trailing zeros in last place */
+			if (z>a && z[-1]) for (i=10, j=0; (z[-1]%(uint32_t)i)==0; i*=10, j++);
+			else j=9;
+			if ((t|32)=='f')
+				p = (int)MIN(p,MAX(0,9*(z-r-1)-j));
+			else
+				p = (int)MIN(p,MAX(0,9*(z-r-1)+e-j));
+		}
+	}
+	if (p > INT_MAX-1-(p || ((unsigned int)fl&ALT_FORM)))
+		return -1;
+	l = 1 + p + (p || ((unsigned int)fl&ALT_FORM));
+	if ((t|32)=='f') {
+		if (e > INT_MAX-l) return -1;
+		if (e>0) l+=e;
+	} else {
+		estr=fmt_u((uintmax_t)(e<0 ? -e : e), ebuf);
+		while(ebuf-estr<2) *--estr='0';
+		*--estr = (e<0 ? '-' : '+');
+		*--estr = (char)t;
+		if (ebuf-estr > INT_MAX-l) return -1;
+		l += (int)(ebuf-estr);
+	}
+
+	if (l > INT_MAX-pl) return -1;
+	pad(&sp, ' ', w, pl+l, fl);
+	out(&sp, prefix, (size_t)pl);
+	pad(&sp, '0', w, pl+l, (int)((unsigned int)fl^ZERO_PAD));
+
+	if ((t|32)=='f') {
+		if (a>r) a=r;
+		for (d=a; d<=r; d++) {
+			s = fmt_u(*d, buf+9); //@@@
+			if (d!=a) while (s>buf) *--s='0';
+			else if (s==buf+9) *--s='0';
+			out(&sp, s, (size_t)(buf+9-s));
+		}
+		if (p || ((unsigned int)fl&ALT_FORM)) out(&sp, ".", 1);
+		for (; d<z && p>0; d++, p-=9) {
+			s = fmt_u(*d, buf+9); //@@@
+			while (s>buf) *--s='0';
+			out(&sp, s, (size_t)(MIN(9,p)));
+		}
+		pad(&sp, '0', p+9, 9, 0);
+	} else {
+		if (z<=a) z=a+1;
+		for (d=a; d<z && p>=0; d++) {
+			s = fmt_u(*d, buf+9); //@@@
+			if (s==buf+9) *--s='0';
+			if (d!=a) while (s>buf) *--s='0';
+			else {
+				out(&sp, s++, 1);
+				if (p>0||((unsigned int)fl&ALT_FORM)) out(&sp, ".", 1);
+			}
+			out(&sp, s, (size_t)(MIN(buf+9-s, p)));
+			p -= (int)(buf+9-s);
+		}
+		pad(&sp, '0', p+18, 18, 0);
+		out(&sp, estr, (size_t)(ebuf-estr));
+	}
+
+	pad(&sp, ' ', w, pl+l, (int)((unsigned int)fl^LEFT_ADJ));
+
+	return MAX(w, pl+l);
+}
+
+/*
+static int getint(char **s) {
+	int i;
+	for (i=0; isdigit(**s); (*s)++) {
+		if (i > INT_MAX/10U || **s-'0' > INT_MAX-10*i) i = -1;
+		else i = 10*i + (**s-'0');
+	}
+	return i;
+}
+*/

+ 69 - 0
deps/musl/vfprintf.h

@@ -0,0 +1,69 @@
+#ifndef VFPRINTF_H
+#define VFPRINTF_H
+
+#include <limits.h>
+#include <string.h>
+//#include <stdarg.h>
+//#include <stddef.h>
+//#include <wchar.h>
+#include <inttypes.h>
+#include <endian.h>
+
+
+//#include <stdio.h>
+
+//#include <float.h>
+#define LDBL_TRUE_MIN 3.6451995318824746025e-4951L
+#define LDBL_MIN     3.3621031431120935063e-4932L
+#define LDBL_MAX     1.1897314953572317650e+4932L
+#define LDBL_EPSILON 1.0842021724855044340e-19L
+
+#define LDBL_MANT_DIG 64
+#define LDBL_MIN_EXP (-16381)
+#define LDBL_MAX_EXP 16384
+//libm.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;
+};
+#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __BIG_ENDIAN
+union ldshape {
+	long double f;
+	struct {
+		uint16_t se;
+		uint16_t top;
+		uint32_t mid;
+		uint64_t lo;
+	} i;
+	struct {
+		uint64_t hi;
+		uint64_t lo;
+	} i2;
+};
+#else
+#error Unsupported long double representation
+#endif
+
+
+int fmt_fp(char*, long double, int, int, int, int);
+
+#endif