path: root/src/core/fops.h

/* 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