a

1 files changed, 1048 insertions, 0 deletions

diff --git a/src/core/fops.h b/src/core/fops.h
new file mode 100644
index 0000000..9213603
--- /dev/null
+++ b/src/core/fops.h
@@ -0,0 +1,1048 @@
+/* IEEE-754 binary64, packed in one 64-bit word.
+   a float on the stack is the raw 64 bits. the odd/even tag is only for the
+   instruction stream, so a float literal is two words: opcode 600 (the even
+   word 1200) and then the raw bits, which the decoder must consume itself
+   because bit 0 is often set.
+
+   two arithmetic paths, chosen when the image is built:
+     - hardware (the default). fop_* are volatile host doubles, so gcc emits
+       the FPU (x86_64 SSE2, AArch64 FP, RV64 D). the arm64 jit and the llvm
+       backend do not call these for 601..609; they emit fadd/fmul/fcmp/scvtf
+       and the matching llvm double ops. bare metal turns the FPU on in the
+       boot stub and compiles this file with a hard-float ABI.
+     - soft (compile with -DZY_SOFTFLOAT). integer IEEE, round ties to even,
+       no host double. a port with no FPU selects this. print and parse stay
+       on the integer side either way, so the text does not depend on libc.
+
+   print and parse are integer, on both paths, so the text does not depend
+   on libc or the host libm. opcode 610 prints six digits after the dot and
+   rounds a tie in the next digit away from zero. opcode 618 (the word 1236)
+   takes (bits prec mode --): mode 0 is the shortest spelling that parses
+   back to the same bits, mode 1 is fixed (prec digits after the dot) and
+   mode 2 is scientific (prec significant digits). prec is capped at 17.
+   the parser keeps the digits it is given, rounds ties to even, and accepts
+   an exponent, inf and nan. */
+#ifndef FOPS_H
+#define FOPS_H
+#include <stdint.h>
+
+#ifdef ZY_SOFTFLOAT
+/* ---- soft-float: binary64 in integer registers, selected with -DZY_SOFTFLOAT ---- */
+typedef struct { uint64_t sig; int exp; int sign; int cls; } FopSF;	/* cls 0 finite, 1 inf, 2 nan */
+
+static int fop_clz(uint64_t x)
+{
+	int n = 0;
+	if (!x) return 64;
+	if (!(x & 0xffffffff00000000ull)) { n += 32; x <<= 32; }
+	if (!(x & 0xffff000000000000ull)) { n += 16; x <<= 16; }
+	if (!(x & 0xff00000000000000ull)) { n += 8; x <<= 8; }
+	if (!(x & 0xf000000000000000ull)) { n += 4; x <<= 4; }
+	if (!(x & 0xc000000000000000ull)) { n += 2; x <<= 2; }
+	if (!(x & 0x8000000000000000ull)) n++;
+	return n;
+}
+static int fop_clz128(unsigned __int128 x)
+{
+	uint64_t hi = (uint64_t)(x >> 64);
+	if (hi) return fop_clz(hi);
+	return 64 + fop_clz((uint64_t)x);
+}
+static void fop_dec(uint64_t u, FopSF *o)
+{
+	int eb = (int)((u >> 52) & 0x7ff);
+	uint64_t frac = u & 0x000fffffffffffffull;
+	o->sign = (int)(u >> 63);
+	o->sig = frac;
+	o->exp = 0;
+	if (eb == 0x7ff) { o->cls = frac ? 2 : 1; return; }
+	o->cls = 0;
+	if (!eb) { o->exp = -1074; return; }
+	o->sig = frac | (1ull << 52);
+	o->exp = eb - 1023 - 52;
+}
+/* value = sig * 2^exp, ties to even. sticky is bits below sig.
+   a magnitude below 2^-1075 is a zero; that exact power is the tie between
+   zero and the smallest subnormal, and it rounds back to zero. */
+static uint64_t fop_sf_round(int sign, int exp, uint64_t sig, int sticky)
+{
+	int round = 0, bit, topw;
+	uint64_t sgn = (uint64_t)sign << 63;
+	if (!sig) return sgn;
+	topw = exp + 63 - fop_clz(sig);
+	if (topw < -1075) return sgn;
+	if (topw == -1075) return (sticky || (sig & (sig - 1))) ? sgn | 1 : sgn;
+	while (sig < (1ull << 52) && exp > -1074) { sig <<= 1; exp--; }
+	while (sig >= (1ull << 53) || (exp < -1074 && sig)) {
+		bit = (int)(sig & 1);
+		sig >>= 1;
+		exp++;
+		sticky = sticky || round;
+		round = bit;
+	}
+	if (round && (sticky || (sig & 1))) {
+		sig++;
+		if (sig >= (1ull << 53)) { sig >>= 1; exp++; }
+	}
+	if (!sig) return sgn;
+	if (exp > 971) return sgn | 0x7ff0000000000000ull;	/* exp + 52 > 1023 */
+	if (sig < (1ull << 52)) return sgn | sig;
+	return sgn | ((uint64_t)(exp + 52 + 1023) << 52) | (sig & 0x000fffffffffffffull);
+}
+static uint64_t fop_add(uint64_t x, uint64_t y)
+{
+	FopSF a, b;
+	unsigned __int128 as, bs, sum;
+	int sh, sticky = 0, sign, lz;
+	uint64_t top;
+	fop_dec(x, &a);
+	fop_dec(y, &b);
+	if (a.cls == 2 || b.cls == 2) return 0x7ff8000000000000ull;
+	if (a.cls == 1 && b.cls == 1) {
+		if (a.sign != b.sign) return 0x7ff8000000000000ull;
+		return ((uint64_t)a.sign << 63) | 0x7ff0000000000000ull;
+	}
+	if (a.cls == 1) return ((uint64_t)a.sign << 63) | 0x7ff0000000000000ull;
+	if (b.cls == 1) return ((uint64_t)b.sign << 63) | 0x7ff0000000000000ull;
+	if (!a.sig && !b.sig) return (a.sign && b.sign) ? 0x8000000000000000ull : 0;
+	if (!a.sig) return y;
+	if (!b.sig) return x;
+	as = (unsigned __int128)a.sig << 60;
+	bs = (unsigned __int128)b.sig << 60;
+	if (a.exp < b.exp) { FopSF t = a; a = b; b = t; as = bs; bs = (unsigned __int128)b.sig << 60; }
+	sh = a.exp - b.exp;
+	if (sh >= 128) { sticky = b.sig != 0; bs = 0; }
+	else if (sh) {
+		sticky = (bs & (((unsigned __int128)1 << sh) - 1)) != 0;
+		bs >>= sh;
+	}
+	if (a.sign == b.sign) { sum = as + bs; sign = a.sign; }
+	else if (as == bs) return 0;
+	else if (as > bs) { sum = as - bs; sign = a.sign; }
+	else { sum = bs - as; sign = b.sign; }
+	lz = fop_clz128(sum);
+	sum <<= lz;
+	top = (uint64_t)(sum >> 64);
+	if ((uint64_t)sum) sticky = 1;
+	return fop_sf_round(sign, a.exp + 4 - lz, top, sticky);
+}
+static uint64_t fop_sub(uint64_t x, uint64_t y) { return fop_add(x, y ^ 0x8000000000000000ull); }
+static uint64_t fop_mul(uint64_t x, uint64_t y)
+{
+	FopSF a, b;
+	unsigned __int128 prod;
+	int sign, lz, sticky = 0;
+	uint64_t top;
+	fop_dec(x, &a);
+	fop_dec(y, &b);
+	sign = a.sign ^ b.sign;
+	if (a.cls == 2 || b.cls == 2) return 0x7ff8000000000000ull;
+	if ((a.cls == 1 && !b.sig) || (b.cls == 1 && !a.sig)) return 0x7ff8000000000000ull;
+	if (a.cls == 1 || b.cls == 1) return ((uint64_t)sign << 63) | 0x7ff0000000000000ull;
+	if (!a.sig || !b.sig) return (uint64_t)sign << 63;
+	prod = (unsigned __int128)a.sig * b.sig;
+	lz = fop_clz128(prod);
+	prod <<= lz;
+	top = (uint64_t)(prod >> 64);
+	if ((uint64_t)prod) sticky = 1;
+	return fop_sf_round(sign, a.exp + b.exp + 64 - lz, top, sticky);
+}
+static uint64_t fop_div(uint64_t x, uint64_t y)
+{
+	FopSF a, b;
+	unsigned __int128 num, q128;
+	int sign, lz, sticky;
+	fop_dec(x, &a);
+	fop_dec(y, &b);
+	sign = a.sign ^ b.sign;
+	if (a.cls == 2 || b.cls == 2) return 0x7ff8000000000000ull;
+	if (a.cls == 1 && b.cls == 1) return 0x7ff8000000000000ull;
+	if (a.cls == 1) return ((uint64_t)sign << 63) | 0x7ff0000000000000ull;
+	if (!b.sig || b.cls == 1) {
+		if (!b.sig && !a.sig) return 0x7ff8000000000000ull;
+		if (!b.sig) return ((uint64_t)sign << 63) | 0x7ff0000000000000ull;
+		return (uint64_t)sign << 63;	/* finite / inf */
+	}
+	if (!a.sig) return (uint64_t)sign << 63;
+	/* (sig << 64) / sig fits in 128 bits even for a subnormal divisor of 1. */
+	num = (unsigned __int128)a.sig << 64;
+	q128 = num / b.sig;
+	sticky = (num % b.sig) != 0;
+	lz = fop_clz128(q128);
+	q128 <<= lz;
+	if ((uint64_t)q128) sticky = 1;
+	return fop_sf_round(sign, a.exp - b.exp - lz, (uint64_t)(q128 >> 64), sticky);
+}
+static int fop_eq(uint64_t x, uint64_t y)
+{
+	FopSF a, b;
+	fop_dec(x, &a);
+	fop_dec(y, &b);
+	if (a.cls == 2 || b.cls == 2) return 0;
+	if (!a.cls && !b.cls && !a.sig && !b.sig) return 1;
+	return x == y;
+}
+static int fop_lt(uint64_t x, uint64_t y)
+{
+	FopSF a, b;
+	int mag;
+	fop_dec(x, &a);
+	fop_dec(y, &b);
+	if (a.cls == 2 || b.cls == 2) return 0;
+	if (a.cls == 1 && b.cls == 1) return a.sign > b.sign;
+	if (a.cls == 1) return a.sign;
+	if (b.cls == 1) return !b.sign;
+	if (!a.sig && !b.sig) return 0;
+	if (a.sign != b.sign) return a.sign;
+	if (a.exp != b.exp) mag = a.exp < b.exp;
+	else mag = a.sig < b.sig;
+	if (a.exp == b.exp && a.sig == b.sig) return 0;
+	return a.sign ? !mag : mag;
+}
+static int fop_le(uint64_t x, uint64_t y) { return fop_eq(x, y) || fop_lt(x, y); }
+static uint64_t fop_i2f(int64_t i)
+{
+	uint64_t mag;
+	if (!i) return 0;
+	if (i < 0) {
+		mag = (i == (int64_t)((uint64_t)1 << 63)) ? ((uint64_t)1 << 63) : (uint64_t)(-i);
+		return fop_sf_round(1, 0, mag, 0);
+	}
+	return fop_sf_round(0, 0, (uint64_t)i, 0);
+}
+static int64_t fop_f2i(uint64_t u)
+{
+	FopSF a;
+	uint64_t mag;
+	fop_dec(u, &a);
+	if (a.cls == 2 || !a.sig) return 0;
+	if (a.cls == 1 || a.exp >= 63) return a.sign ? (int64_t)((uint64_t)1 << 63) : 9223372036854775807LL;
+	if (a.exp >= 0) {
+		if (a.exp > 0 && a.sig > (~0ull >> a.exp)) return a.sign ? (int64_t)((uint64_t)1 << 63) : 9223372036854775807LL;
+		mag = a.sig << a.exp;
+	} else {
+		if (-a.exp >= 64) return 0;
+		mag = a.sig >> -a.exp;
+	}
+	if (!a.sign) {
+		if (mag > 9223372036854775807ull) return 9223372036854775807LL;
+		return (int64_t)mag;
+	}
+	if (mag > 9223372036854775808ull) return (int64_t)((uint64_t)1 << 63);
+	if (mag == 9223372036854775808ull) return (int64_t)((uint64_t)1 << 63);
+	return -(int64_t)mag;
+}
+
+#else
+static double fop_d(uint64_t u)
+{
+	union { uint64_t u; double d; } x;
+	x.u = u;
+	return x.d;
+}
+static uint64_t fop_u(double d)
+{
+	union { uint64_t u; double d; } x;
+	x.d = d;
+	return x.u;
+}
+static uint64_t fop_add(uint64_t x, uint64_t y)
+{
+	volatile double a = fop_d(x), b = fop_d(y), c = a + b;
+	return fop_u(c);
+}
+static uint64_t fop_sub(uint64_t x, uint64_t y)
+{
+	volatile double a = fop_d(x), b = fop_d(y), c = a - b;
+	return fop_u(c);
+}
+static uint64_t fop_mul(uint64_t x, uint64_t y)
+{
+	volatile double a = fop_d(x), b = fop_d(y), c = a * b;
+	return fop_u(c);
+}
+static uint64_t fop_div(uint64_t x, uint64_t y)
+{
+	volatile double a = fop_d(x), b = fop_d(y), c = a / b;
+	return fop_u(c);
+}
+static int fop_eq(uint64_t x, uint64_t y)
+{
+	volatile double a = fop_d(x), b = fop_d(y);
+	return a == b;
+}
+static int fop_lt(uint64_t x, uint64_t y)
+{
+	volatile double a = fop_d(x), b = fop_d(y);
+	return a < b;
+}
+static int fop_le(uint64_t x, uint64_t y)
+{
+	volatile double a = fop_d(x), b = fop_d(y);
+	return a <= b;
+}
+static uint64_t fop_i2f(int64_t i)
+{
+	volatile double d = (double)i;
+	return fop_u(d);
+}
+/* trunc toward 0. values outside int64 saturate, and a nan becomes 0.
+   tests stay in range, where this matches a hardware fcvtzs. */
+static int64_t fop_f2i(uint64_t x)
+{
+	volatile double d = fop_d(x);
+	if (!(d == d)) return 0;
+	if (d >= 9223372036854775808.0) return 9223372036854775807LL;
+	if (d <= -9223372036854775808.0) return (int64_t)((uint64_t)1 << 63);
+	return (int64_t)d;
+}
+#endif
+
+/* binary64 <-> binary32, in integer arithmetic so every build agrees (no libm,
+   no compiler runtime). f2s rounds to nearest even and takes a double's bits to
+   a single's bits (in the low 32); s2f is exact. a nan keeps its sign and the top
+   of its payload, and comes out quiet. */
+static uint64_t fop_f2s(uint64_t x)
+{
+	uint64_t sign = (x >> 63) << 31, frac = x & 0xFFFFFFFFFFFFFull, sig, q, rem, half;
+	int e = (int)((x >> 52) & 0x7FF), E, shift;
+	if (e == 0x7FF) {
+		if (!frac) return sign | 0x7F800000u;
+		return sign | 0x7FC00000u | (frac >> 29);
+	}
+	if (e == 0) return sign;	/* zero, or a subnormal double: far below the smallest single */
+	E = e - 1023;
+	sig = frac | (1ull << 52);
+	if (E > 127) return sign | 0x7F800000u;
+	shift = 29;
+	if (E < -126) { shift = 29 + (-126 - E); if (shift > 54) return sign; }
+	q = sig >> shift;
+	rem = sig & ((1ull << shift) - 1);
+	half = 1ull << (shift - 1);
+	if (rem > half || (rem == half && (q & 1))) q++;
+	if (shift == 29) {
+		if (q == (1ull << 24)) { q >>= 1; E++; if (E > 127) return sign | 0x7F800000u; }
+		return sign | ((uint64_t)(E + 127) << 23) | (q & 0x7FFFFF);
+	}
+	return sign | q;	/* a subnormal; if it rounded up to the smallest normal these are that normal's bits */
+}
+static uint64_t fop_s2f(uint64_t x)
+{
+	uint64_t sign, frac;
+	int e, p;
+	x &= 0xFFFFFFFFull;
+	sign = (x >> 31) << 63;
+	frac = x & 0x7FFFFF;
+	e = (int)((x >> 23) & 0xFF);
+	if (e == 0xFF) {
+		if (!frac) return sign | 0x7FF0000000000000ull;
+		return sign | 0x7FF8000000000000ull | (frac << 29);
+	}
+	if (e == 0) {
+		if (!frac) return sign;
+		p = 22;
+		while (!(frac >> p)) p--;
+		frac = (frac << (52 - p)) & 0xFFFFFFFFFFFFFull;
+		return sign | ((uint64_t)(p - 149 + 1023) << 52) | frac;
+	}
+	return sign | ((uint64_t)(e - 127 + 1023) << 52) | (frac << 29);
+}
+
+/* decimal conversion in base 2^32. 96 limbs is room for a 120-digit
+   significand times 5^445, which covers every double and the digits a
+   127-character literal can hold. copies go through a volatile store so a
+   freestanding build does not emit memcpy. */
+#define FOP_BI 96
+typedef struct { uint32_t d[FOP_BI]; int n; } FopBi;
+
+static void bi_norm(FopBi *a)
+{
+	while (a->n > 0 && !a->d[a->n - 1]) a->n--;
+}
+static void bi_clear(FopBi *a) { a->n = 0; }
+static void bi_copy(FopBi *dst, const FopBi *src)
+{
+	volatile uint32_t *w = dst->d;
+	int i;
+	dst->n = src->n;
+	for (i = 0; i < src->n; i++) w[i] = src->d[i];
+}
+static void bi_set(FopBi *a, uint64_t v)
+{
+	a->n = 0;
+	if (!v) return;
+	a->d[0] = (uint32_t)v;
+	a->n = 1;
+	if (v >> 32) { a->d[1] = (uint32_t)(v >> 32); a->n = 2; }
+}
+static void bi_add_u32(FopBi *a, uint32_t v)
+{
+	uint64_t c = v;
+	int i = 0;
+	if (!v) return;
+	while (c) {
+		if (i >= FOP_BI) return;
+		if (i == a->n) a->d[a->n++] = 0;
+		c += a->d[i];
+		a->d[i++] = (uint32_t)c;
+		c >>= 32;
+	}
+}
+static void bi_mul_u32(FopBi *a, uint32_t m)
+{
+	uint64_t c = 0;
+	int i;
+	if (!m || !a->n) { a->n = 0; return; }
+	for (i = 0; i < a->n; i++) {
+		c += (uint64_t)a->d[i] * m;
+		a->d[i] = (uint32_t)c;
+		c >>= 32;
+	}
+	if (c && a->n < FOP_BI) a->d[a->n++] = (uint32_t)c;
+}
+static void bi_mul(FopBi *a, const FopBi *b)
+{
+	FopBi r;
+	int i, j;
+	uint64_t c;
+	if (!a->n || !b->n) { a->n = 0; return; }
+	bi_clear(&r);
+	r.n = a->n + b->n;
+	if (r.n > FOP_BI) r.n = FOP_BI;
+	for (i = 0; i < r.n; i++) r.d[i] = 0;
+	for (i = 0; i < a->n; i++) {
+		c = 0;
+		for (j = 0; j < b->n && i + j < FOP_BI; j++) {
+			c += (uint64_t)r.d[i + j] + (uint64_t)a->d[i] * b->d[j];
+			r.d[i + j] = (uint32_t)c;
+			c >>= 32;
+		}
+		if (c && i + b->n < r.n) r.d[i + b->n] = (uint32_t)c;
+	}
+	bi_norm(&r);
+	bi_copy(a, &r);
+}
+static void bi_pow5(FopBi *out, int e)
+{
+	FopBi base, acc;
+	bi_set(&acc, 1);
+	bi_set(&base, 5);
+	if (e < 0) e = 0;
+	while (e) {
+		if (e & 1) bi_mul(&acc, &base);
+		e >>= 1;
+		if (e) bi_mul(&base, &base);
+	}
+	bi_copy(out, &acc);
+}
+static int bi_bitlen(const FopBi *a)
+{
+	uint32_t x;
+	int n;
+	if (!a->n) return 0;
+	x = a->d[a->n - 1];
+	n = (a->n - 1) * 32;
+	while (x) { n++; x >>= 1; }
+	return n;
+}
+static int bi_bit(const FopBi *a, int i)
+{
+	if (i < 0 || (i >> 5) >= a->n) return 0;
+	return (int)((a->d[i >> 5] >> (i & 31)) & 1);
+}
+static int bi_low_any(const FopBi *a, int nbits)
+{
+	int i, w, b;
+	if (nbits <= 0) return 0;
+	w = nbits >> 5;
+	b = nbits & 31;
+	for (i = 0; i < w && i < a->n; i++) if (a->d[i]) return 1;
+	if (b && w < a->n && (a->d[w] & ((1u << b) - 1))) return 1;
+	return 0;
+}
+static uint64_t bi_extract(const FopBi *a, int lo, int nbits)
+{
+	uint64_t r = 0;
+	int i;
+	if (nbits > 63) nbits = 63;
+	for (i = 0; i < nbits; i++) if (bi_bit(a, lo + i)) r |= (uint64_t)1 << i;
+	return r;
+}
+static void bi_shl(FopBi *a, int s)
+{
+	int w, b, i;
+	uint32_t prev;
+	if (s <= 0 || !a->n) return;
+	w = s >> 5;
+	b = s & 31;
+	if (a->n + w + (b ? 1 : 0) > FOP_BI) {
+		for (i = 0; i < FOP_BI; i++) a->d[i] = 0xffffffffu;
+		a->n = FOP_BI;
+		return;
+	}
+	if (b) {
+		prev = 0;
+		for (i = 0; i < a->n; i++) {
+			uint32_t cur = a->d[i];
+			a->d[i] = (cur << b) | prev;
+			prev = cur >> (32 - b);
+		}
+		if (prev) a->d[a->n++] = prev;
+	}
+	if (w) {
+		for (i = a->n - 1; i >= 0; i--) a->d[i + w] = a->d[i];
+		for (i = 0; i < w; i++) a->d[i] = 0;
+		a->n += w;
+	}
+}
+static int bi_cmp(const FopBi *a, const FopBi *b)
+{
+	int i;
+	if (a->n != b->n) return a->n < b->n ? -1 : 1;
+	for (i = a->n - 1; i >= 0; i--) if (a->d[i] != b->d[i]) return a->d[i] < b->d[i] ? -1 : 1;
+	return 0;
+}
+static void bi_sub(FopBi *a, const FopBi *b)
+{
+	uint64_t br = 0;
+	int i;
+	for (i = 0; i < a->n; i++) {
+		uint64_t av = a->d[i];
+		uint64_t bv = (uint64_t)(i < b->n ? b->d[i] : 0) + br;
+		if (av < bv) { a->d[i] = (uint32_t)(av + (1ull << 32) - bv); br = 1; }
+		else { a->d[i] = (uint32_t)(av - bv); br = 0; }
+	}
+	bi_norm(a);
+}
+/* q = floor(num/den), rem = num % den. den must be nonzero. */
+static void bi_divmod(const FopBi *num, const FopBi *den, FopBi *q, FopBi *rem)
+{
+	int i, nb = bi_bitlen(num);
+	bi_clear(q);
+	bi_clear(rem);
+	for (i = nb - 1; i >= 0; i--) {
+		bi_shl(rem, 1);
+		if (bi_bit(num, i)) bi_add_u32(rem, 1);
+		if (bi_cmp(rem, den) >= 0) {
+			bi_sub(rem, den);
+			{
+				int w = i >> 5;
+				while (q->n <= w) q->d[q->n++] = 0;
+				q->d[w] |= 1u << (i & 31);
+			}
+		}
+	}
+	bi_norm(q);
+}
+/* half up: remainder >= half the divisor increments the quotient. */
+static void bi_div_halfup(const FopBi *num, const FopBi *den, FopBi *q)
+{
+	FopBi rem, twice;
+	bi_divmod(num, den, q, &rem);
+	bi_copy(&twice, &rem);
+	bi_shl(&twice, 1);
+	if (bi_cmp(&twice, den) >= 0) bi_add_u32(q, 1);
+}
+static uint32_t bi_div10(FopBi *a)
+{
+	uint64_t rem = 0;
+	int i;
+	for (i = a->n - 1; i >= 0; i--) {
+		rem = (rem << 32) | a->d[i];
+		a->d[i] = (uint32_t)(rem / 10);
+		rem %= 10;
+	}
+	bi_norm(a);
+	return (uint32_t)rem;
+}
+static int bi_to_digits(const FopBi *a, char *dig)
+{
+	FopBi x;
+	char rev[400];
+	int n = 0, i;
+	bi_copy(&x, a);
+	if (!x.n) { dig[0] = '0'; dig[1] = 0; return 1; }
+	while (x.n && n < 390) rev[n++] = (char)('0' + (int)bi_div10(&x));
+	for (i = 0; i < n; i++) dig[i] = rev[n - 1 - i];
+	dig[n] = 0;
+	return n;
+}
+/* right shift by s, then +1 if the highest discarded bit was set. */
+static void bi_shr_halfup(const FopBi *a, int s, FopBi *dst)
+{
+	int up, w, b, i, len;
+	if (s <= 0) { bi_copy(dst, a); return; }
+	len = bi_bitlen(a);
+	up = bi_bit(a, s - 1);
+	if (s >= len) {
+		bi_clear(dst);
+		if (up) bi_set(dst, 1);
+		return;
+	}
+	w = s >> 5;
+	b = s & 31;
+	dst->n = 0;
+	if (!b) {
+		for (i = w; i < a->n; i++) dst->d[dst->n++] = a->d[i];
+	} else {
+		for (i = w; i < a->n; i++) {
+			uint32_t cur = a->d[i] >> b;
+			if (i + 1 < a->n) cur |= a->d[i + 1] << (32 - b);
+			dst->d[dst->n++] = cur;
+		}
+	}
+	bi_norm(dst);
+	if (up) bi_add_u32(dst, 1);
+}
+
+/* q * 2^scale, plus a nonzero fraction below bit 0 when extra is set.
+   round ties to even. a pure power of two at 2^-1075 is the tie between
+   zero and the smallest subnormal, and it rounds back to zero. */
+static uint64_t fop_qbits(const FopBi *q, int scale, int extra, int neg)
+{
+	uint64_t sgn = (uint64_t)neg << 63;
+	int top, unbiased, low, roundb, sticky;
+	uint64_t sig;
+	if (!q->n) return sgn;
+	top = bi_bitlen(q) - 1;
+	unbiased = top + scale;
+	if (unbiased > 1023) return sgn | 0x7ff0000000000000ull;
+	if (unbiased < -1075) return sgn;
+	if (unbiased == -1075) {
+		if (!extra && !bi_low_any(q, top)) return sgn;
+		return sgn | 1;
+	}
+	if (unbiased < -1022) {
+		int rsh = -(scale + 1074);
+		if (rsh >= 0) {
+			sig = bi_extract(q, rsh, top - rsh + 1);
+			roundb = rsh > 0 ? bi_bit(q, rsh - 1) : 0;
+			sticky = extra || (rsh > 1 && bi_low_any(q, rsh - 1));
+		} else {
+			sig = bi_extract(q, 0, top + 1) << (-rsh);
+			roundb = 0;
+			sticky = extra;
+		}
+		if (roundb && (sticky || (sig & 1))) sig++;
+		if (sig >= (1ull << 52)) return sgn | (1ull << 52);
+		return sgn | sig;
+	}
+	low = top - 52;
+	if (low <= 0) {
+		sig = bi_extract(q, 0, top + 1);
+		if (low < 0) sig <<= -low;
+		roundb = 0;
+		sticky = extra;
+	} else {
+		sig = bi_extract(q, low, 53);
+		roundb = bi_bit(q, low - 1);
+		sticky = extra || bi_low_any(q, low - 1);
+	}
+	if (roundb && (sticky || (sig & 1))) {
+		sig++;
+		if (sig == (1ull << 53)) {
+			sig >>= 1;
+			unbiased++;
+			if (unbiased > 1023) return sgn | 0x7ff0000000000000ull;
+		}
+	}
+	return sgn | ((uint64_t)(unbiased + 1023) << 52) | (sig & 0xfffffffffffffull);
+}
+
+/* sig * 10^exp10, with sticky meaning the true value is strictly larger. */
+static uint64_t fop_dec_bits(const FopBi *sig, int exp10, int neg, int sticky)
+{
+	FopBi num, den, q, rem;
+	int bexp, approx, target, sh, guard = 0;
+	int sgn = neg ? 1 : 0;
+	if (!sig->n) return sgn ? 0x8000000000000000ull : 0;
+	if (exp10 >= 0) {
+		bi_copy(&num, sig);
+		bi_pow5(&den, exp10);
+		bi_mul(&num, &den);
+		bi_set(&den, 1);
+		bexp = exp10;
+	} else {
+		bi_copy(&num, sig);
+		bi_pow5(&den, -exp10);
+		bexp = exp10;
+	}
+	approx = bi_bitlen(&num) - bi_bitlen(&den) + bexp;
+	if (approx > 1080) return sgn ? 0xfff0000000000000ull : 0x7ff0000000000000ull;
+	if (approx < -1140) return sgn ? 0x8000000000000000ull : 0;
+again:
+	if (approx < -902) target = -1088;
+	else target = approx - 80;
+	target -= guard * 32;
+	sh = bexp - target;
+	{
+		FopBi n2, d2;
+		bi_copy(&n2, &num);
+		bi_copy(&d2, &den);
+		if (sh >= 0) bi_shl(&n2, sh);
+		else bi_shl(&d2, -sh);
+		bi_divmod(&n2, &d2, &q, &rem);
+		{
+			int top = q.n ? bi_bitlen(&q) - 1 : 0;
+			int low = top - 52;
+			int unbiased = q.n ? top + target : -9999;
+			if (q.n && unbiased >= -1022 && low < 4 && guard < 4) { guard++; goto again; }
+		}
+		return fop_qbits(&q, target, sticky || rem.n != 0, sgn);
+	}
+}
+
+/* parse directly from byte text or vm words, without truncating a lexer buffer. */
+typedef struct { const char *bytes; const int64_t *words; uint64_t length; } FopText;
+static char fop_char(FopText *s, uint64_t at)
+{
+	if (s->words) return at < s->length ? (char)(s->words[at] & 255) : 0;
+	return s->bytes[at];
+}
+static int fop_word(FopText *s, uint64_t p, const char *w)
+{
+	while (*w) {
+		char c = fop_char(s, p);
+		if (c >= 'A' && c <= 'Z') c = (char)(c - 'A' + 'a');
+		if (c != *w) return 0;
+		p++;
+		w++;
+	}
+	return fop_char(s, p) == 0;
+}
+
+static uint64_t fop_parse_text(FopText *s)
+{
+	uint64_t p = 0;
+	int neg = 0, exp10 = 0, nexp = 0, eneg = 0, frac = 0, saw = 0, nd = 0, sticky = 0;
+	int dec;
+	FopBi sig;
+	if (fop_char(s, p) == '+') p++;
+	else if (fop_char(s, p) == '-') { neg = 1; p++; }
+	if (fop_word(s, p, "inf")) return neg ? 0xfff0000000000000ull : 0x7ff0000000000000ull;
+	if (fop_word(s, p, "nan")) return neg ? 0xfff8000000000000ull : 0x7ff8000000000000ull;
+	bi_clear(&sig);
+	while (fop_char(s, p) == '.' || (fop_char(s, p) >= '0' && fop_char(s, p) <= '9')) {
+		if (fop_char(s, p) == '.') {
+			if (frac) break;
+			frac = 1;
+			p++;
+			continue;
+		}
+		saw = 1;
+		if (fop_char(s, p) == '0' && !sig.n) {
+			if (frac) exp10--;
+		} else if (nd < 80) {
+			bi_mul_u32(&sig, 10);
+			bi_add_u32(&sig, (uint32_t)(fop_char(s, p) - '0'));
+			nd++;
+			if (frac) exp10--;
+		} else {
+			if (fop_char(s, p) != '0') sticky = 1;
+			if (!frac) exp10++;
+		}
+		p++;
+	}
+	if (fop_char(s, p) == 'e' || fop_char(s, p) == 'E') {
+		p++;
+		if (fop_char(s, p) == '+') p++;
+		else if (fop_char(s, p) == '-') { eneg = 1; p++; }
+		while (fop_char(s, p) >= '0' && fop_char(s, p) <= '9') {
+			if (nexp < 1000000) nexp = nexp * 10 + (fop_char(s, p) - '0');
+			p++;
+		}
+		if (eneg) exp10 -= nexp;
+		else exp10 += nexp;
+	}
+	if (!saw || !sig.n) return neg ? 0x8000000000000000ull : 0;
+	dec = exp10 + nd - 1;
+	if (dec > 309) return neg ? 0xfff0000000000000ull : 0x7ff0000000000000ull;
+	if (dec < -326) return neg ? 0x8000000000000000ull : 0;
+	return fop_dec_bits(&sig, exp10, neg, sticky);
+}
+
+static uint64_t fop_parse(const char *bytes)
+{
+	FopText text = { bytes, 0, 0 };
+	return fop_parse_text(&text);
+}
+static uint64_t fop_parse_words(const int64_t *words, uint64_t length)
+{
+	FopText text = { 0, words, length };
+	return length ? fop_parse_text(&text) : 0;
+}
+
+static void fop_decode(uint64_t bits, uint64_t *sig, int *e2, int *neg, int *cls)
+{
+	uint64_t exp = (bits >> 52) & 0x7ff;
+	uint64_t frac = bits & 0x000fffffffffffffull;
+	*neg = (int)(bits >> 63);
+	*cls = 0;
+	*sig = 0;
+	*e2 = 0;
+	if (exp == 0x7ff) { *cls = frac ? 2 : 1; return; }
+	if (!exp) {
+		*sig = frac;
+		*e2 = -1074;
+		return;
+	}
+	*sig = frac | (1ull << 52);
+	*e2 = (int)exp - 1023 - 52;
+}
+
+/* x ? 10^k. x = sig * 2^e2, sig > 0. */
+static int fop_cmp10(uint64_t sig, int e2, int k)
+{
+	FopBi a, b;
+	int q;
+	if (k >= 0) {
+		bi_set(&a, sig);
+		bi_pow5(&b, k);
+		if (e2 >= k) bi_shl(&a, e2 - k);
+		else bi_shl(&b, k - e2);
+		return bi_cmp(&a, &b);
+	}
+	q = -k;
+	bi_set(&a, sig);
+	bi_pow5(&b, q);
+	bi_mul(&a, &b);
+	if (e2 + q >= 0) {
+		bi_shl(&a, e2 + q);
+		if (!a.n) return -1;
+		if (a.n == 1 && a.d[0] == 1) return 0;
+		return 1;
+	}
+	bi_set(&b, 1);
+	bi_shl(&b, -(e2 + q));
+	return bi_cmp(&a, &b);
+}
+
+/* no __builtin_clzll: a freestanding riscv link has no __clzdi2. */
+static int fop_nlz(uint64_t x)
+{
+	int n = 0;
+	if (!x) return 64;
+	if (!(x & 0xffffffff00000000ull)) { n += 32; x <<= 32; }
+	if (!(x & 0xffff000000000000ull)) { n += 16; x <<= 16; }
+	if (!(x & 0xff00000000000000ull)) { n += 8; x <<= 8; }
+	if (!(x & 0xf000000000000000ull)) { n += 4; x <<= 4; }
+	if (!(x & 0xc000000000000000ull)) { n += 2; x <<= 2; }
+	if (!(x & 0x8000000000000000ull)) n++;
+	return n;
+}
+static int fop_decexp(uint64_t sig, int e2)
+{
+	int log2 = 63 - fop_nlz(sig) + e2;
+	int est, guard = 0;
+	if (log2 >= 0) est = (int)(((unsigned __int128)log2 * 1292913986ull) >> 32);
+	else {
+		unsigned __int128 m = (unsigned __int128)(-log2) * 1292913986ull;
+		est = -(int)(m >> 32);
+		if (m & 0xffffffffull) est--;
+	}
+	while (guard++ < 12 && fop_cmp10(sig, e2, est + 1) >= 0) est++;
+	guard = 0;
+	while (guard++ < 12 && fop_cmp10(sig, e2, est) < 0) est--;
+	return est;
+}
+
+/* round_half_up(sig * 2^e2 * 10^scale) as an integer. */
+static void fop_scale(uint64_t sig, int e2, int scale, FopBi *mag)
+{
+	FopBi num, den;
+	int i, q10;
+	if (scale >= 0) {
+		bi_set(&num, sig);
+		for (i = 0; i < scale; i++) bi_mul_u32(&num, 5);
+		if (e2 + scale >= 0) bi_shl(&num, e2 + scale);
+		else { bi_shr_halfup(&num, -(e2 + scale), mag); return; }
+		bi_copy(mag, &num);
+		return;
+	}
+	q10 = -scale;
+	bi_set(&num, sig);
+	bi_pow5(&den, q10);
+	if (e2 >= q10) bi_shl(&num, e2 - q10);
+	else bi_shl(&den, q10 - e2);
+	bi_div_halfup(&num, &den, mag);
+}
+
+static int fop_putu(char *out, int n, int v)
+{
+	char b[16];
+	int i = 0, j;
+	if (v < 0) v = -v;
+	if (!v) { out[n++] = '0'; return n; }
+	while (v) { b[i++] = (char)('0' + v % 10); v /= 10; }
+	for (j = i - 1; j >= 0; j--) out[n++] = b[j];
+	return n;
+}
+
+/* plain when the exponent is in [-4, 15], otherwise scientific. */
+static int fop_spell(const char *dig, int nd, int dec, int neg, int scientific, char *out)
+{
+	int n = 0, i;
+	if (neg) out[n++] = '-';
+	if (!scientific) {
+		if (dec >= 0) {
+			int whole = dec + 1;
+			if (nd <= whole) {
+				for (i = 0; i < nd; i++) out[n++] = dig[i];
+				for (; i < whole; i++) out[n++] = '0';
+			} else {
+				for (i = 0; i < whole; i++) out[n++] = dig[i];
+				out[n++] = '.';
+				for (; i < nd; i++) out[n++] = dig[i];
+			}
+		} else {
+			out[n++] = '0';
+			out[n++] = '.';
+			for (i = 0; i < -dec - 1; i++) out[n++] = '0';
+			for (i = 0; i < nd; i++) out[n++] = dig[i];
+		}
+	} else {
+		out[n++] = dig[0];
+		if (nd > 1) {
+			out[n++] = '.';
+			for (i = 1; i < nd; i++) out[n++] = dig[i];
+		}
+		out[n++] = 'e';
+		out[n++] = dec < 0 ? '-' : '+';
+		n = fop_putu(out, n, dec < 0 ? -dec : dec);
+	}
+	out[n] = 0;
+	return n;
+}
+
+static void fop_sigdigs(uint64_t sig, int e2, int prec, char *dig, int *dec)
+{
+	int d = fop_decexp(sig, e2);
+	FopBi mag;
+	char tmp[48];
+	int nd = 0, i, guard = 0;
+	while (guard++ < 8) {
+		int zeros = 1;
+		fop_scale(sig, e2, prec - 1 - d, &mag);
+		nd = bi_to_digits(&mag, tmp);
+		if (nd == prec) break;
+		if (nd == prec + 1) {
+			for (i = 1; i < nd; i++) if (tmp[i] != '0') zeros = 0;
+			if (tmp[0] == '1' && zeros) {
+				d++;
+				dig[0] = '1';
+				for (i = 1; i < prec; i++) dig[i] = '0';
+				dig[prec] = 0;
+				*dec = d;
+				return;
+			}
+		}
+		if (nd > prec) d += nd - prec;
+		else d -= prec - nd;
+	}
+	if (nd > prec) nd = prec;
+	for (i = 0; i < prec - nd; i++) dig[i] = '0';
+	for (i = 0; i < nd && i < prec; i++) dig[prec - nd + i] = tmp[i];
+	dig[prec] = 0;
+	*dec = d;
+}
+
+static int fop_format_fixed(uint64_t sig, int e2, int neg, int prec, char *out)
+{
+	FopBi t, mag;
+	char dig[400];
+	int n = 0, nd, i, whole;
+	bi_set(&t, sig ? sig : 0);
+	if (!sig) bi_clear(&mag);
+	else {
+		int i5;
+		for (i5 = 0; i5 < prec; i5++) bi_mul_u32(&t, 5);
+		if (e2 + prec >= 0) { bi_shl(&t, e2 + prec); bi_copy(&mag, &t); }
+		else bi_shr_halfup(&t, -(e2 + prec), &mag);
+	}
+	if (neg) out[n++] = '-';
+	nd = bi_to_digits(&mag, dig);
+	if (nd <= prec) {
+		out[n++] = '0';
+		if (prec) {
+			out[n++] = '.';
+			for (i = 0; i < prec - nd; i++) out[n++] = '0';
+			for (i = 0; i < nd; i++) out[n++] = dig[i];
+		}
+	} else {
+		whole = nd - prec;
+		for (i = 0; i < whole && n < 380; i++) out[n++] = dig[i];
+		if (prec && n < 390) {
+			out[n++] = '.';
+			for (i = 0; i < prec; i++) out[n++] = dig[whole + i];
+		}
+	}
+	out[n] = 0;
+	return n;
+}
+
+static int fop_special_text(int cls, int neg, char *out)
+{
+	int n = 0;
+	if (cls == 2) { out[0] = 'n'; out[1] = 'a'; out[2] = 'n'; out[3] = 0; return 3; }
+	if (neg) out[n++] = '-';
+	out[n++] = 'i'; out[n++] = 'n'; out[n++] = 'f';
+	out[n] = 0;
+	return n;
+}
+
+/* mode 0 shortest, 1 fixed, 2 scientific. out must hold 400 bytes. */
+static int fop_format_mode(uint64_t bits, int64_t mode, int64_t prec, char *out)
+{
+	uint64_t sig;
+	int e2, neg, cls, n, k, dec, scientific;
+	char dig[24], buf[80];
+	fop_decode(bits, &sig, &e2, &neg, &cls);
+	if (cls) return fop_special_text(cls, neg, out);
+	if (mode == 1) {
+		if (prec < 0) prec = 0;
+		if (prec > 17) prec = 17;
+		return fop_format_fixed(sig, e2, neg, (int)prec, out);
+	}
+	if (!sig) {
+		if (mode == 2) {
+			int p = (int)prec;
+			int i;
+			if (p < 1) p = 1;
+			if (p > 17) p = 17;
+			n = 0;
+			if (neg) out[n++] = '-';
+			out[n++] = '0';
+			if (p > 1) {
+				out[n++] = '.';
+				for (i = 1; i < p; i++) out[n++] = '0';
+			}
+			out[n++] = 'e'; out[n++] = '+'; out[n++] = '0';
+			out[n] = 0;
+			return n;
+		}
+		out[0] = '0';
+		if (neg) { out[0] = '-'; out[1] = '0'; out[2] = 0; return 2; }
+		out[1] = 0;
+		return 1;
+	}
+	if (mode == 2) {
+		int p = (int)prec;
+		if (p < 1) p = 1;
+		if (p > 17) p = 17;
+		fop_sigdigs(sig, e2, p, dig, &dec);
+		return fop_spell(dig, p, dec, neg, 1, out);
+	}
+	for (k = 1; k <= 18; k++) {
+		fop_sigdigs(sig, e2, k, dig, &dec);
+		scientific = !(dec >= -4 && dec <= 15);
+		fop_spell(dig, k, dec, neg, scientific, buf);
+		if (k == 18 || fop_parse(buf) == bits) return fop_spell(dig, k, dec, neg, scientific, out);
+	}
+	return fop_spell(dig, 18, dec, neg, 1, out);
+}
+
+static int fop_format(uint64_t bits, char *out)
+{
+	return fop_format_mode(bits, 1, 6, out);
+}
+#endif