diff options
Diffstat (limited to 'src/core/fops.h')
| -rw-r--r-- | src/core/fops.h | 1048 |
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 |