Blame for firmware/graphcalc.c

// Carbon In-Weights Graphing & Symbolic/Numeric Calculator Firmware (RV32IM + Zba + Bignum).
// Executes directly inside compiled Qwen3.5 Gated-DeltaNet + SwiGLU layers (layers[0..10]).
// Strictly integer & fixed-point (Q32.32 + 128-digit hardware bignum register file) — zero FPU.

#include <stdint.h>
#include <stddef.h>
#include "rt.h"
#include "graphcalc_tables.h"

#ifdef QEMU
// Software emulation of chipformer's custom-0/custom-1 128-digit bignum unit for qemu-riscv32.
#define BN_OPS(n, off, m, d, s, ovw)                                                        \
    ((uint32_t)(d) | (uint32_t)(s) << 3 | (uint32_t)(ovw) << 6 | (uint32_t)(n) << 7 |       \
     (uint32_t)(off) << 13 | (uint32_t)(m) << 20)
#define BN_OFF(off) ((uint32_t)(off) << 13)
#define BN_M(m) ((uint32_t)(m) << 20)

static uint8_t sw_bn[8][128];

static uint32_t bn_mac(uint32_t ops) {
    uint32_t dd = ops & 7, ds = (ops >> 3) & 7, ovw = (ops >> 6) & 1;
    uint32_t n = (ops >> 7) & 63, off = (ops >> 13) & 127, m = (ops >> 20) & 15;
    if (m > 9) m = 9;
    uint32_t c = 0;
    for (uint32_t i = 0; i < n; i++) {
        uint32_t r = i + off;
        uint32_t d = r < 128 ? sw_bn[dd][r] : 0;
        uint32_t v = (ovw ? 0 : d) + sw_bn[ds][i] * m + c;
        c = v / 10;
        if (r < 128) sw_bn[dd][r] = (uint8_t)(v % 10);
    }
    return c;
}

static uint32_t bn_msb(uint32_t ops) {
    uint32_t dd = ops & 7, ds = (ops >> 3) & 7;
    uint32_t n = (ops >> 7) & 63, off = (ops >> 13) & 127, m = (ops >> 20) & 15;
    if (m > 9) m = 9;
    int32_t c = 0;
    for (uint32_t i = 0; i < n; i++) {
        uint32_t r = i + off;
        int32_t d = r < 128 ? (int32_t)sw_bn[dd][r] : 0;
        int32_t v = d - (int32_t)(sw_bn[ds][i] * m) - c;
        int32_t rem = ((v % 10) + 10) % 10;
        c = (rem - v) / 10;
        if (r < 128) sw_bn[dd][r] = (uint8_t)rem;
    }
    return (uint32_t)c;
}

static void bn_ld(const char *end, uint32_t ops) {
    uint32_t dd = ops & 7, n = (ops >> 7) & 63, off = (ops >> 13) & 127;
    for (uint32_t i = 0; i < n; i++) {
        uint32_t r = i + off;
        uint32_t v = (uint32_t)((unsigned char)end[-1 - (int)i]) & 15u;
        if (v > 9) v = 9;
        if (r < 128) sw_bn[dd][r] = (uint8_t)v;
    }
}

static uint32_t bn_get(uint32_t ops) {
    uint32_t dd = ops & 7, off = (ops >> 13) & 127;
    return sw_bn[dd][off];
}

static uint32_t bn_scan(const char *s) {
    uint32_t n = 0;
    while (n < 63 && (unsigned char)s[n] >= '0' && (unsigned char)s[n] <= '9') n++;
    return n;
}

static void bn_print(int r, int len) {
    for (int i = len - 1; i >= 0; i--) rv_putchar((char)('0' + sw_bn[r][i]));
}
#else
#include "bn.h"
#endif

// ============================================================================
// Basic Output & String Helpers (no libc printf bloat)
// ============================================================================

static void puts_raw(const char *s) {
    while (*s) rv_putchar(*s++);
}

static void putc_rep(char c, int n) {
    for (int i = 0; i < n; i++) rv_putchar(c);
}

static int slen(const char *s) {
    int n = 0;
    while (s[n]) n++;
    return n;
}

static int streq(const char *a, const char *b) {
    while (*a && *b && *a == *b) { a++; b++; }
    return *a == 0 && *b == 0;
}

static int streq_ci(const char *a, const char *b) {
    while (*a && *b) {
        char ca = *a++, cb = *b++;
        if (ca >= 'A' && ca <= 'Z') ca = (char)(ca + 32);
        if (cb >= 'A' && cb <= 'Z') cb = (char)(cb + 32);
        if (ca != cb) return 0;
    }
    return *a == 0 && *b == 0;
}

static int starts_with_ci(const char *s, const char *p) {
    while (*p) {
        char cs = *s++, cp = *p++;
        if (cs >= 'A' && cs <= 'Z') cs = (char)(cs + 32);
        if (cp >= 'A' && cp <= 'Z') cp = (char)(cp + 32);
        if (cs != cp) return 0;
    }
    return 1;
}

static int is_space(char c) {
    return c == ' ' || c == '\t' || c == '\r' || c == '\n';
}

static int is_alpha(char c) {
    return (c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z') || c == '_';
}

static int is_digit(char c) {
    return c >= '0' && c <= '9';
}

static void copy_trim(char *dst, int cap, const char *src, int len) {
    int a = 0, b = len;
    while (a < b && is_space(src[a])) a++;
    while (b > a && is_space(src[b - 1])) b--;
    if (b - a >= 2 && ((src[a] == '"' && src[b - 1] == '"') ||
                       (src[a] == '\'' && src[b - 1] == '\'') ||
                       (src[a] == '$' && src[b - 1] == '$'))) {
        a++;
        b--;
    }
    int o = 0;
    for (int i = a; i < b && o < cap - 1; i++) {
        dst[o++] = src[i] == '\n' || src[i] == '\r' ? ' ' : src[i];
    }
    dst[o] = 0;
}

static void append_str(char *buf, int cap, int *pos, const char *s) {
    while (*s && *pos < cap - 1) buf[(*pos)++] = *s++;
    buf[*pos] = 0;
}

static void append_char(char *buf, int cap, int *pos, char c) {
    if (*pos < cap - 1) {
        buf[(*pos)++] = c;
        buf[*pos] = 0;
    }
}

// ============================================================================
// Engine A: 62-Digit Hardware Bignum + Exact Reduced Rationals (bn.h)
// ============================================================================

#define MAXD 62
#define ONE_REG 7

typedef uint32_t num;
#define REG(a) ((a) & 7u)
#define NEG(a) (((a) >> 8) & 1u)
#define LEN(a) ((int)((a) >> 16))
#define NUM(r, neg, len) ((uint32_t)(r) | ((uint32_t)(neg) << 8) | ((uint32_t)(len) << 16))

static uint32_t bn_dirty[8];
static uint32_t bn_used;
static int bn_err;
static char one_ch = '1';  // MUST stay in RAM (.data), NOT const (.rodata)!

static char bn_steps[512];
static int bn_step_pos;
static int bn_step_cnt;

static int bn_dig(num a, int i) { return (int)bn_get(BN_OPS(0, i, 0, REG(a), 0, 0)); }

static void bn_wrote(num a, int n) {
    if ((uint32_t)n > bn_dirty[REG(a)]) bn_dirty[REG(a)] = (uint32_t)n;
}

static num bn_alloc(int len) {
    for (int r = 0; r < ONE_REG; r++) {
        if (!(bn_used >> r & 1u)) {
            bn_used |= 1u << r;
            for (uint32_t o = 0; o < bn_dirty[r]; o += 63) {
                uint32_t chunk = bn_dirty[r] - o < 63 ? bn_dirty[r] - o : 63;
                bn_mac(BN_OPS(chunk, o, 0, r, r, 1));
            }
            bn_dirty[r] = 0;
            return NUM(r, 0, len);
        }
    }
    bn_err = 2;
    return NUM(0, 0, len);
}

static void bn_release(num a) { bn_used &= ~(1u << REG(a)); }

static int bn_is_zero(num a) { return LEN(a) == 1 && !bn_dig(a, 0); }

static num bn_trim(num a) {
    int n = LEN(a);
    while (n > 1 && !bn_dig(a, n - 1)) n--;
    a = (a & 0xFFFFu) | ((uint32_t)n << 16);
    return (n == 1 && !bn_dig(a, 0)) ? (a & ~0x100u) : a;
}

static int bn_cmp_mag(num a, num b) {
    if (LEN(a) != LEN(b)) return LEN(a) < LEN(b) ? -1 : 1;
    for (int i = LEN(a) - 1; i >= 0; i--) {
        int x = bn_dig(a, i), y = bn_dig(b, i);
        if (x != y) return x < y ? -1 : 1;
    }
    return 0;
}

static num bn_add_mag(num a, num b) {
    int n = (LEN(a) > LEN(b) ? LEN(a) : LEN(b)) + 1;
    if (n > MAXD) bn_err = 1;
    bn_mac(BN_OPS(n, 0, 1, REG(a), REG(b), 0));
    bn_wrote(a, n);
    if (!bn_dig(a, n - 1)) n--;
    return (a & 0xFFFFu) | ((uint32_t)n << 16);
}

static num bn_sub_mag(num a, num b) {
    bn_msb(BN_OPS(LEN(a), 0, 1, REG(a), REG(b), 0));
    return bn_trim(a);
}

__attribute__((noinline)) static num bn_add(num a, num b) {
    num s;
    if (NEG(a) == NEG(b)) {
        s = bn_add_mag(a, b);
        bn_release(b);
    } else if (bn_cmp_mag(a, b) >= 0) {
        s = bn_sub_mag(a, b);
        bn_release(b);
    } else {
        s = bn_sub_mag(b, a);
        bn_release(a);
    }
    return s;
}

__attribute__((noinline)) static num bn_mul(num a, num b) {
    if (LEN(a) < LEN(b)) { num t = a; a = b; b = t; }
    int la = LEN(a), lb = LEN(b);
    num c = bn_alloc(la + lb);
    if (la + lb > MAXD) {
        bn_err = 1;
        return c;
    }
    uint32_t base = BN_OPS(la + 1, 0, 0, REG(c), REG(a), 0);
    for (int j = 0; j < lb; j++) {
        uint32_t d = (uint32_t)bn_dig(b, j);
        if (d) bn_mac(base | BN_OFF(j) | BN_M(d));
    }
    bn_wrote(c, la + lb);
    return bn_trim(c | ((NEG(a) ^ NEG(b)) << 8));
}

typedef struct { num q, r; } bn_qr;

__attribute__((noinline)) static bn_qr bn_divmod(num a, num b) {
    if (bn_is_zero(b)) bn_err = 1;
    if (bn_err || bn_cmp_mag(a, b) < 0) return (bn_qr){bn_alloc(1), a & ~0x100u};
    int la = LEN(a), lb = LEN(b), ra = (int)REG(a);
    num q = bn_alloc(la - lb + 1);
    int b1 = bn_dig(b, lb - 1), den = lb >= 2 ? b1 * 10 + bn_dig(b, lb - 2) : b1;
    uint32_t msb = BN_OPS(lb + 1, 0, 0, ra, REG(b), 0), mac = BN_OPS(lb + 1, 0, 1, ra, REG(b), 0);
    uint32_t setq = BN_OPS(1, 0, 0, REG(q), ONE_REG, 1);
    for (int t = la - lb; t >= 0; t--) {
        int top = t + lb;
        int w = bn_dig(a, top) * 10 + bn_dig(a, top - 1);
        if (lb >= 2) w = w * 10 + bn_dig(a, top - 2);
        int qd = (lb >= 3 ? w + 1 : w) / den;
        if (qd > 9) qd = 9;
        if (qd) {
            uint32_t o = BN_OFF(t);
            if (bn_msb(msb | o | BN_M(qd)))
                do qd--;
                while (!bn_mac(mac | o) && qd > 0);
            if (qd) bn_mac(setq | o | BN_M(qd));
        }
    }
    bn_wrote(q, la - lb + 1);
    return (bn_qr){bn_trim(q), bn_trim(NUM(ra, 0, lb))};
}

__attribute__((noinline)) static num bn_divide(num a, num b, int want_mod) {
    uint32_t sa = NEG(a), sb = NEG(b);
    bn_qr d = bn_divmod(a, b);
    if (sa != sb && !bn_is_zero(d.r)) {
        d.q = bn_add_mag(d.q, NUM(ONE_REG, 0, 1));
        num r2 = bn_sub_mag(b & ~0x100u, d.r);
        bn_release(d.r);
        d.r = r2;
    } else {
        bn_release(b);
    }
    bn_release(want_mod ? d.q : d.r);
    num out = want_mod ? d.r : d.q;
    out = (out & ~0x100u) | ((want_mod ? sb : (sa ^ sb)) << 8);
    return bn_is_zero(out) ? (out & ~0x100u) : out;
}

__attribute__((noinline)) static num bn_gcd(num a, num b) {
    a &= ~0x100u;
    b &= ~0x100u;
    while (!bn_err && !bn_is_zero(b)) {
        bn_qr d = bn_divmod(a, b);
        bn_release(d.q);
        a = b;
        b = d.r;
    }
    bn_release(b);
    return a;
}

__attribute__((noinline)) static num bn_power(num a, num e) {
    int x = 0;
    if (NEG(e) || LEN(e) > 3) bn_err = 1;
    for (int i = LEN(e) - 1; i >= 0 && !bn_err; i--) x = x * 10 + bn_dig(e, i);
    bn_release(e);
    num res = bn_alloc(1);
    bn_mac(BN_OPS(1, 0, 1, REG(res), ONE_REG, 1));
    bn_wrote(res, 1);
    while (x && !bn_err) {
        if (x & 1) {
            num t = bn_mul(res, a);
            bn_release(res);
            res = t;
        }
        x >>= 1;
        if (x) {
            num t = bn_mul(a, a);
            bn_release(a);
            a = t;
        }
    }
    bn_release(a);
    return res;
}

static num bn_vals[16];
static num bn_vden[16];
static int bn_hasden[16];
static int bn_vsc[16];
static char bn_ops[16];

static int bn_prec(int op) {
    return op == '^' ? 4 : op == 'n' ? 3 : (op == '*' || op == 'f' || op == 'd' || op == '%') ? 2 : 1;
}

static num bn_copy_of(num a) {
    num c = bn_alloc(LEN(a));
    bn_mac(BN_OPS(LEN(a), 0, 1, REG(c), REG(a), 1));
    bn_wrote(c, LEN(a));
    return NUM(REG(c), NEG(a), LEN(a));
}

__attribute__((noinline)) static void bn_reduce(int i) {
    num n = bn_vals[i], d = bn_vden[i];
    bn_hasden[i] = 0;
    if (bn_err) return;
    if (bn_is_zero(n)) {
        bn_release(d);
        return;
    }
    uint32_t neg = NEG(n);
    num g = bn_gcd(bn_copy_of(n), bn_copy_of(d));
    bn_qr x = bn_divmod(n & ~0x100u, g);
    bn_release(x.r);
    bn_qr y = bn_divmod(d, g);
    bn_release(y.r);
    bn_release(g);
    bn_vals[i] = x.q | (neg << 8);
    if (LEN(y.q) == 1 && bn_dig(y.q, 0) == 1) {
        bn_release(y.q);
    } else {
        bn_vden[i] = y.q;
        bn_hasden[i] = 1;
    }
}

__attribute__((noinline)) static void bn_frac_op(int op, int i, num a, num da, int fa, num b, num db, int fb) {
    num n, d;
    if (op == '*') {
        n = bn_mul(a, b);
        bn_release(a);
        bn_release(b);
        if (fa && fb) {
            d = bn_mul(da, db);
            bn_release(da);
            bn_release(db);
        } else {
            d = fa ? da : db;
        }
    } else if (op == 'd') {
        if (bn_is_zero(b)) bn_err = 1;
        if (fb) { n = bn_mul(a, db); bn_release(a); bn_release(db); }
        else n = a;
        if (fa) { d = bn_mul(da, b); bn_release(da); bn_release(b); }
        else d = b;
        if (NEG(d)) {
            d &= ~0x100u;
            if (!bn_is_zero(n)) n ^= 0x100u;
        }
    } else {
        num t1 = a, t2 = b;
        if (fb) { t1 = bn_mul(a, db); bn_release(a); }
        if (fa) { t2 = bn_mul(b, da); bn_release(b); }
        if (fa && fb) { d = bn_mul(da, db); bn_release(da); bn_release(db); }
        else d = fa ? da : db;
        n = bn_add(t1, (op == '-' && !bn_is_zero(t2)) ? (t2 ^ 0x100u) : t2);
    }
    bn_vals[i] = n;
    bn_vden[i] = d;
    bn_reduce(i);
}

static num bn_shift(num a, int k) {
    int n = LEN(a) + k;
    if (n > MAXD) { bn_err = 1; return a; }
    num c = bn_alloc(n);
    bn_mac(BN_OPS(LEN(a), k, 1, REG(c), REG(a), 1));
    bn_wrote(c, n);
    bn_release(a);
    return NUM(REG(c), NEG(a), n);
}

static void bn_fmt_slot(int idx, char *buf, int cap) {
    int p = 0;
    num v = bn_vals[idx];
    int sc = bn_vsc[idx], n = LEN(v);
    if (NEG(v)) append_char(buf, cap, &p, '-');
    if (!sc) {
        for (int i = n - 1; i >= 0; i--) append_char(buf, cap, &p, (char)('0' + bn_dig(v, i)));
    } else {
        int z = 0;
        while (z < sc && !bn_dig(v, z)) z++;
        if (n <= sc) append_char(buf, cap, &p, '0');
        for (int i = n > sc ? n - 1 : sc - 1; i >= z; i--) {
            if (i == sc - 1) append_char(buf, cap, &p, '.');
            append_char(buf, cap, &p, (char)('0' + (i < n ? bn_dig(v, i) : 0)));
        }
    }
    if (bn_hasden[idx]) {
        num dv = bn_vden[idx];
        append_char(buf, cap, &p, '/');
        for (int i = LEN(dv) - 1; i >= 0; i--) append_char(buf, cap, &p, (char)('0' + bn_dig(dv, i)));
    }
}

static int bn_apply(int op, int nv) {
    if (op == 'n') {
        num v = bn_vals[nv - 1];
        bn_vals[nv - 1] = bn_is_zero(v) ? v : (v ^ 0x100u);
        return nv;
    }
    if (nv < 2) { bn_err = 1; return nv; }

    char sa_str[72], sb_str[72], sr_str[72];
    int can_log = (!bn_err && bn_step_cnt < 8 && LEN(bn_vals[nv - 2]) <= 30 && LEN(bn_vals[nv - 1]) <= 30);
    if (can_log) {
        bn_fmt_slot(nv - 2, sa_str, sizeof sa_str);
        bn_fmt_slot(nv - 1, sb_str, sizeof sb_str);
    }

    int sb = bn_vsc[nv - 1], sa = bn_vsc[nv - 2];
    int fb = bn_hasden[nv - 1], fa = bn_hasden[nv - 2];
    num b = bn_vals[--nv], a = bn_vals[nv - 1], r;
    int rational = op == 'd' || ((fa || fb) && (op == '+' || op == '-' || op == '*'));
    if ((fa || fb) && !rational) bn_err = 1;
    if (rational && (sa || sb)) bn_err = 1;
    if (bn_err) {
        r = a;
        bn_hasden[nv - 1] = 0;
        bn_vals[nv - 1] = r;
        bn_vsc[nv - 1] = sa;
    } else if (rational) {
        bn_frac_op(op, nv - 1, a, bn_vden[nv - 1], fa, b, bn_vden[nv], fb);
        bn_vsc[nv - 1] = 0;
    } else {
        if (op == '+' || op == '-') {
            if (sa < sb) { a = bn_shift(a, sb - sa); sa = sb; }
            if (sb < sa) b = bn_shift(b, sa - sb);
            r = bn_add(a, (op == '-' && !bn_is_zero(b)) ? (b ^ 0x100u) : b);
        } else if (op == '*') {
            r = bn_mul(a, b);
            bn_release(a);
            bn_release(b);
            sa += sb;
        } else if (sa || sb) {
            bn_err = 1;
            r = a;
        } else if (op == '^') {
            r = bn_power(a, b);
        } else if (op == 'g') {
            r = bn_gcd(a, b);
        } else {
            r = bn_divide(a, b, op == '%');
        }
        bn_hasden[nv - 1] = 0;
        bn_vals[nv - 1] = r;
        bn_vsc[nv - 1] = sa;
    }

    if (can_log && !bn_err && LEN(bn_vals[nv - 1]) <= 62) {
        bn_fmt_slot(nv - 1, sr_str, sizeof sr_str);
        bn_step_cnt++;
        append_str(bn_steps, sizeof bn_steps, &bn_step_pos, "step ");
        append_char(bn_steps, sizeof bn_steps, &bn_step_pos, (char)('0' + bn_step_cnt));
        append_str(bn_steps, sizeof bn_steps, &bn_step_pos, ": ");
        if (op == 'g') {
            append_str(bn_steps, sizeof bn_steps, &bn_step_pos, "gcd(");
            append_str(bn_steps, sizeof bn_steps, &bn_step_pos, sa_str);
            append_str(bn_steps, sizeof bn_steps, &bn_step_pos, ", ");
            append_str(bn_steps, sizeof bn_steps, &bn_step_pos, sb_str);
            append_str(bn_steps, sizeof bn_steps, &bn_step_pos, ") = ");
        } else {
            const char *op_s = op == 'd' ? "/" : op == 'f' ? "//" : op == '^' ? "^" :
                               op == '*' ? "*" : op == '+' ? "+" : op == '-' ? "-" : "%";
            append_str(bn_steps, sizeof bn_steps, &bn_step_pos, sa_str);
            append_char(bn_steps, sizeof bn_steps, &bn_step_pos, ' ');
            append_str(bn_steps, sizeof bn_steps, &bn_step_pos, op_s);
            append_char(bn_steps, sizeof bn_steps, &bn_step_pos, ' ');
            append_str(bn_steps, sizeof bn_steps, &bn_step_pos, sb_str);
            append_str(bn_steps, sizeof bn_steps, &bn_step_pos, " = ");
        }
        append_str(bn_steps, sizeof bn_steps, &bn_step_pos, sr_str);
        append_char(bn_steps, sizeof bn_steps, &bn_step_pos, '\n');
    }
    return nv;
}

// Try evaluating a pure arithmetic/rational/bignum expression in Engine A.
// Returns 0 on success (having printed the structured trace + exact answer), 1 if unsupported/failed.
static int run_bignum_engine(const char *expr_ram) {
    bn_used = 0;
    bn_err = 0;
    bn_step_pos = 0;
    bn_step_cnt = 0;
    bn_steps[0] = 0;
    for (int r = 0; r < ONE_REG; r++) {
        bn_mac(BN_OPS(63, 0, 0, r, r, 1));
        bn_mac(BN_OPS(63, 63, 0, r, r, 1));
        bn_dirty[r] = 0;
    }
    bn_ld(&one_ch + 1, BN_OPS(1, 0, 0, ONE_REG, 0, 0));

    const char *s = expr_ram;
    int operand = 1, nv = 0, no = 0;
    for (;;) {
        char c = *s;
        if (c == ' ' || c == '\t') { s++; continue; }
        if (nv >= 15 || no >= 15) bn_err = 1;
        if (bn_err) break;
        if (operand) {
            if (is_digit(c)) {
                int n = (int)bn_scan(s);
                const char *end = s + n;
                int frac = 0;
                if (*end == '.' && is_digit(end[1])) {
                    static char dbuf[MAXD + 4];
                    frac = (int)bn_scan(end + 1);
                    if (n + frac > MAXD) { bn_err = 1; break; }
                    for (int i = 0; i < n; i++) dbuf[i] = s[i];
                    for (int i = 0; i < frac; i++) dbuf[n + i] = end[1 + i];
                    s = end + 1 + frac;
                    n += frac;
                    end = dbuf + n;
                } else {
                    s = end;
                }
                if (n > MAXD) { bn_err = 1; break; }
                num a = bn_alloc(n);
                bn_ld(end, BN_OPS(n, 0, 0, REG(a), 0, 0));
                bn_wrote(a, n);
                bn_vsc[nv] = frac;
                bn_hasden[nv] = 0;
                bn_vals[nv++] = bn_dig(a, n - 1) ? a : bn_trim(a);
                operand = 0;
            } else if (c == '(') {
                bn_ops[no++] = '(';
                s++;
            } else if (c == '-' || c == '+') {
                if (c == '-') bn_ops[no++] = 'n';
                s++;
            } else if (starts_with_ci(s, "gcd(")) {
                bn_ops[no++] = 'g';
                s += 4;
            } else {
                bn_err = 1;
            }
            continue;
        }
        // Implicit multiplication before '(' when preceded by number or ')'
        if (c == '(') {
            int pr = bn_prec('*');
            while (no && bn_ops[no - 1] != '(' && bn_ops[no - 1] != 'g' && bn_prec(bn_ops[no - 1]) >= pr)
                nv = bn_apply(bn_ops[--no], nv);
            bn_ops[no++] = '*';
            bn_ops[no++] = '(';
            s++;
            operand = 1;
            continue;
        }
        if (c == ')' || c == ',') {
            while (no && bn_ops[no - 1] != '(' && bn_ops[no - 1] != 'g') nv = bn_apply(bn_ops[--no], nv);
            if (!no) {
                bn_err = 1;
            } else if (c == ')') {
                if (bn_ops[--no] == 'g') nv = bn_apply('g', nv);
            } else if (bn_ops[no - 1] != 'g') {
                bn_err = 1;
            } else {
                operand = 1;
            }
            s++;
            continue;
        }
        if (!c) break;
        int op = (unsigned char)c, len = 1;
        if (c == '*' && s[1] == '*') {
            op = '^';
            len = 2;
        } else if (c == '/') {
            op = s[1] == '/' ? 'f' : 'd';
            len = s[1] == '/' ? 2 : 1;
        } else if (c != '+' && c != '-' && c != '*' && c != '%' && c != '^') {
            bn_err = 1;
            continue;
        }
        int pr = bn_prec(op);
        while (no && bn_ops[no - 1] != '(' && bn_ops[no - 1] != 'g' &&
               (bn_prec(bn_ops[no - 1]) > pr || (bn_prec(bn_ops[no - 1]) == pr && op != '^')))
            nv = bn_apply(bn_ops[--no], nv);
        bn_ops[no++] = (char)op;
        s += len;
        operand = 1;
    }
    while (!bn_err && no) {
        int op = bn_ops[--no];
        if (op == '(' || op == 'g') bn_err = 1;
        else nv = bn_apply(op, nv);
    }
    if (bn_err || nv != 1 || operand) return 1;

    puts_raw("expression: ");
    puts_raw(expr_ram);
    rv_putchar('\n');
    if (bn_step_cnt >= 2) puts_raw(bn_steps);

    puts_raw("exact value: ");
    num v = bn_vals[0];
    int sc = bn_vsc[0], n = LEN(v);
    if (NEG(v)) rv_putchar('-');
    if (!sc) {
        bn_print((int)REG(v), n);
    } else {
        int z = 0;
        while (z < sc && !bn_dig(v, z)) z++;
        if (n <= sc) rv_putchar('0');
        for (int i = n > sc ? n - 1 : sc - 1; i >= z; i--) {
            if (i == sc - 1) rv_putchar('.');
            rv_putchar((char)('0' + (i < n ? bn_dig(v, i) : 0)));
        }
    }
    if (bn_hasden[0]) {
        num den = bn_vden[0];
        rv_putchar('/');
        bn_print((int)REG(den), LEN(den));
        // Also print 6-decimal expansion when numerator*10^6 fits in MAXD
        if (n + 6 <= MAXD) {
            uint32_t neg = NEG(v);
            num scaled = bn_shift(bn_copy_of(v & ~0x100u), 6);
            bn_qr dec = bn_divmod(scaled, den);
            int exact_dec = bn_is_zero(dec.r);
            puts_raw(exact_dec ? " (= " : " (~ ");
            if (neg) rv_putchar('-');
            int dn = LEN(dec.q);
            int z = 0;
            while (z < 6 && exact_dec && (z >= dn || !bn_dig(dec.q, z))) z++;
            if (dn <= 6) rv_putchar('0');
            for (int i = dn > 6 ? dn - 1 : 5; i >= z; i--) {
                if (i == 5) rv_putchar('.');
                rv_putchar((char)('0' + (i < dn ? bn_dig(dec.q, i) : 0)));
            }
            rv_putchar(')');
            bn_release(dec.q);
            bn_release(dec.r);
        }
    }
    rv_putchar('\n');
    return 0;
}

// ============================================================================
// Engine B: Q32.32 Fixed-Point Scientific, Equation Solver & Graphing Engine
// ============================================================================

typedef int64_t q32_t;
#define Q32_MAX 0x7FFFFFFFFFFFFFFFLL
#define Q32_MIN (-0x7FFFFFFFFFFFFFFFLL - 1LL)

static int q_err;

static q32_t q32_from_int(int64_t x) {
    if (x > 2147483647LL || x < -2147483648LL) { q_err = 1; return x > 0 ? Q32_MAX : Q32_MIN; }
    return x << 32;
}

static q32_t q32_add(q32_t a, q32_t b) {
    uint64_t r = (uint64_t)a + (uint64_t)b;
    if (((a ^ (q32_t)r) & (b ^ (q32_t)r)) < 0) {
        q_err = 1;
        return a > 0 ? Q32_MAX : Q32_MIN;
    }
    return (q32_t)r;
}

static q32_t q32_sub(q32_t a, q32_t b) {
    uint64_t r = (uint64_t)a - (uint64_t)b;
    if (((a ^ b) & (a ^ (q32_t)r)) < 0) {
        q_err = 1;
        return a > 0 ? Q32_MAX : Q32_MIN;
    }
    return (q32_t)r;
}

static q32_t q32_mul(q32_t a, q32_t b) {
    int neg = (a < 0) ^ (b < 0);
    uint64_t ua = a < 0 ? -(uint64_t)a : (uint64_t)a;
    uint64_t ub = b < 0 ? -(uint64_t)b : (uint64_t)b;
    uint32_t al = (uint32_t)ua, ah = (uint32_t)(ua >> 32);
    uint32_t bl = (uint32_t)ub, bh = (uint32_t)(ub >> 32);
    if (ah && bh) {
        uint64_t hh = (uint64_t)ah * bh;
        if (hh >= 0x80000000ULL) { q_err = 1; return neg ? Q32_MIN : Q32_MAX; }
        uint64_t p_ll = (uint64_t)al * bl;
        uint64_t p_hl = (uint64_t)ah * bl;
        uint64_t p_lh = (uint64_t)al * bh;
        uint64_t mid = (p_ll >> 32) + (p_hl & 0xFFFFFFFFULL) + (p_lh & 0xFFFFFFFFULL);
        uint64_t hi = hh + (p_hl >> 32) + (p_lh >> 32) + (mid >> 32);
        if (hi >= 0x80000000ULL) { q_err = 1; return neg ? Q32_MIN : Q32_MAX; }
        uint64_t res = (hi << 32) | (mid & 0xFFFFFFFFULL);
        res += (p_ll >> 31) & 1ULL;
        return neg ? -(q32_t)res : (q32_t)res;
    }
    uint64_t p_ll = (uint64_t)al * bl;
    uint64_t p_hl = (uint64_t)ah * bl;
    uint64_t p_lh = (uint64_t)al * bh;
    uint64_t mid = (p_ll >> 32) + (p_hl & 0xFFFFFFFFULL) + (p_lh & 0xFFFFFFFFULL);
    uint64_t hi = (p_hl >> 32) + (p_lh >> 32) + (mid >> 32);
    if (hi >= 0x80000000ULL) { q_err = 1; return neg ? Q32_MIN : Q32_MAX; }
    uint64_t res = (hi << 32) | (mid & 0xFFFFFFFFULL);
    res += (p_ll >> 31) & 1ULL;
    return neg ? -(q32_t)res : (q32_t)res;
}

static q32_t q32_div(q32_t a, q32_t b) {
    if (b == 0) { q_err = 1; return 0; }
    int neg = (a < 0) ^ (b < 0);
    uint64_t ua = a < 0 ? -(uint64_t)a : (uint64_t)a;
    uint64_t ub = b < 0 ? -(uint64_t)b : (uint64_t)b;
    uint64_t q_int = ua / ub;
    if (q_int >= 0x80000000ULL) { q_err = 1; return neg ? Q32_MIN : Q32_MAX; }
    uint64_t rem = ua % ub;
    uint32_t frac = 0;
    if (ub <= 0xFFFFFFFFULL) {
        frac = (uint32_t)(((rem << 32) + (ub >> 1)) / ub);
        uint64_t res = (q_int << 32) + frac;
        return neg ? -(q32_t)res : (q32_t)res;
    }
    for (int i = 0; i < 32; i++) {
        uint64_t top = rem >> 63;
        rem <<= 1;
        frac <<= 1;
        if (top || rem >= ub) {
            rem -= ub;
            frac |= 1u;
        }
    }
    uint64_t top = rem >> 63;
    rem <<= 1;
    uint64_t res = (q_int << 32) | frac;
    if (top || rem >= ub) res++;
    return neg ? -(q32_t)res : (q32_t)res;
}

static q32_t q32_sqrt(q32_t x) {
    if (x < 0) { q_err = 1; return 0; }
    if (x == 0) return 0;
    uint64_t ux = (uint64_t)x;
    uint64_t rem = 0, root = 0;
    for (int i = 47; i >= 0; i--) {
        uint64_t two = i >= 16 ? ((ux >> (2 * (i - 16))) & 3ULL) : 0ULL;
        rem = (rem << 2) | two;
        uint64_t cand = (root << 2) | 1ULL;
        if (rem >= cand) {
            rem -= cand;
            root = (root << 1) | 1ULL;
        } else {
            root <<= 1;
        }
    }
    return (q32_t)root;
}

static void q32_sincos(q32_t theta, q32_t *out_s, q32_t *out_c) {
    q32_t r = theta;
    if (r > Q32_PI || r < -Q32_PI) {
        int save = q_err;
        int64_t turns = q32_mul(r, Q32_INV_TWO_PI) >> 32;
        q_err = save;
        r -= turns * Q32_TWO_PI;
        while (r > Q32_PI) r -= Q32_TWO_PI;
        while (r < -Q32_PI) r += Q32_TWO_PI;
    }
    int flip = 0;
    if (r > Q32_HALF_PI) { r = Q32_PI - r; flip = 1; }
    else if (r < -Q32_HALF_PI) { r = -Q32_PI - r; flip = 1; }

    q32_t x = CORDIC_K, y = 0, z = r;
    for (int i = 0; i < 30; i++) {
        q32_t dx = y >> i, dy = x >> i;
        if (z >= 0) { x -= dx; y += dy; z -= CORDIC_ATAN[i]; }
        else        { x += dx; y -= dy; z += CORDIC_ATAN[i]; }
    }
    if (flip) x = -x;
    if (x > -16 && x < 16) x = 0;
    if (y > -16 && y < 16) y = 0;
    *out_s = y;
    *out_c = x;
}

static q32_t q32_ln(q32_t x) {
    if (x <= 0) { q_err = 1; return 0; }
    uint64_t m = (uint64_t)x;
    int k = 0;
    while (m >= (2ULL << 32)) { m >>= 1; k++; }
    while (m < (1ULL << 32)) { m <<= 1; k--; }
    uint64_t prod = 1ULL << 32;
    q32_t acc = 0;
    for (int i = 1; i <= 30; i++) {
        uint64_t nxt = prod + (prod >> i);
        if (nxt <= m) {
            prod = nxt;
            acc += LN1P_POW2[i];
        }
    }
    if (m > prod) acc += q32_div((q32_t)(m - prod), (q32_t)prod);
    return (int64_t)k * Q32_LN2 + acc;
}

static q32_t q32_exp(q32_t x) {
    if (x > 92288000000LL) { q_err = 1; return Q32_MAX; }
    if (x < -96000000000LL) return 0;
    q32_t kq = q32_mul(x, Q32_INV_LN2);
    int32_t k = (int32_t)(kq >> 32);
    q32_t r = x - (int64_t)k * Q32_LN2;
    while (r < 0) { r += Q32_LN2; k--; }
    while (r >= Q32_LN2) { r -= Q32_LN2; k++; }
    uint64_t prod = 1ULL << 32;
    for (int i = 1; i <= 30; i++) {
        if (r >= LN1P_POW2[i]) {
            r -= LN1P_POW2[i];
            prod += prod >> i;
        }
    }
    if (r > 0) prod += (prod * (uint64_t)r) >> 32;
    if (k >= 31) { q_err = 1; return Q32_MAX; }
    if (k >= 0) {
        uint64_t res = prod << k;
        if (res >= 0x8000000000000000ULL) { q_err = 1; return Q32_MAX; }
        return (q32_t)res;
    }
    return (-k >= 63) ? 0 : (q32_t)(prod >> (-k));
}

static q32_t q32_ipow(q32_t a, int32_t e) {
    int inv = 0;
    if (e < 0) { inv = 1; e = -e; }
    q32_t res = Q32_ONE;
    while (e > 0 && !q_err) {
        if (e & 1) res = q32_mul(res, a);
        e >>= 1;
        if (e) a = q32_mul(a, a);
    }
    return inv ? q32_div(Q32_ONE, res) : res;
}

static q32_t q32_pow(q32_t a, q32_t b) {
    if ((b & 0xFFFFFFFFLL) == 0) {
        int64_t e = b >> 32;
        if (e >= -64 && e <= 64) return q32_ipow(a, (int32_t)e);
    }
    if (a == 0) {
        if (b <= 0) q_err = 1;
        return 0;
    }
    if (a < 0) { q_err = 1; return 0; }
    return q32_exp(q32_mul(b, q32_ln(a)));
}

static int parse_q32_num(const char *s, const char **endp, q32_t *out) {
    while (*s == ' ') s++;
    int neg = 0;
    if (*s == '-' || *s == '+') { neg = (*s == '-'); s++; }
    while (*s == ' ') s++;
    if (starts_with_ci(s, "pi") && !is_alpha(s[2])) {
        *out = neg ? -Q32_PI : Q32_PI;
        if (endp) *endp = s + 2;
        return 1;
    }
    if (!is_digit(*s) && !(*s == '.' && is_digit(s[1]))) return 0;
    uint64_t ip = 0;
    while (is_digit(*s)) {
        ip = ip * 10 + (uint64_t)(*s - '0');
        if (ip > 2147483647ULL) return 0;
        s++;
    }
    uint64_t fp_num = 0, fp_den = 1;
    if (*s == '.') {
        s++;
        while (is_digit(*s)) {
            if (fp_den < 100000000000ULL) {
                fp_num = fp_num * 10 + (uint64_t)(*s - '0');
                fp_den *= 10;
            }
            s++;
        }
    }
    uint64_t frac = 0;
    if (fp_num > 0) {
        // Compute (fp_num << 32) / fp_den cleanly via 32-bit binary division
        uint64_t rem = fp_num;
        for (int i = 0; i < 32; i++) {
            rem <<= 1;
            frac <<= 1;
            if (rem >= fp_den) { rem -= fp_den; frac |= 1ULL; }
        }
        if ((rem << 1) >= fp_den) frac++;
    }
    q32_t v = (q32_t)((ip << 32) + frac);
    // Optional /denom for inline bounds like "-pi/2" or "3/2"
    if (*s == '/' && is_digit(s[1])) {
        s++;
        int64_t d = 0;
        while (is_digit(*s) && d < 100000000LL) { d = d * 10 + (*s - '0'); s++; }
        if (d > 0) v = q32_div(v, q32_from_int(d));
    }
    *out = neg ? -v : v;
    if (endp) *endp = s;
    return 1;
}

static void fmt_q32(q32_t v, char *buf, int cap, int decimals) {
    int p = 0;
    if (v > -64 && v < 64) v = 0;  // snap sub-1.5e-8 noise to 0
    int neg = v < 0;
    uint64_t uv = neg ? -(uint64_t)v : (uint64_t)v;
    uint32_t scale = 1;
    for (int i = 0; i < decimals; i++) scale *= 10u;
    uint32_t ip = (uint32_t)(uv >> 32);
    uint64_t fl = uv & 0xFFFFFFFFULL;
    // Rounded decimal fraction: ((fl * scale) + 2^31) >> 32
    uint64_t prod = fl * (uint64_t)scale + (1ULL << 31);
    uint32_t dec = (uint32_t)(prod >> 32);
    if (dec >= scale) { ip++; dec -= scale; }
    if (ip == 0 && dec == 0) neg = 0;
    if (neg) append_char(buf, cap, &p, '-');

    char itmp[16];
    int in = 0;
    if (ip == 0) itmp[in++] = '0';
    while (ip > 0) { itmp[in++] = (char)('0' + (ip % 10u)); ip /= 10u; }
    while (in > 0) append_char(buf, cap, &p, itmp[--in]);

    if (decimals > 0 && dec > 0) {
        char dtmp[12];
        for (int i = decimals - 1; i >= 0; i--) {
            dtmp[i] = (char)('0' + (dec % 10u));
            dec /= 10u;
        }
        int len = decimals;
        while (len > 0 && dtmp[len - 1] == '0') len--;
        if (len > 0) {
            append_char(buf, cap, &p, '.');
            for (int i = 0; i < len; i++) append_char(buf, cap, &p, dtmp[i]);
        }
    }
}

// Continued-fraction recovery for exact small rationals p/q (1 < q <= 240, error <= 6 ULPs).
static void fmt_q32_rich(q32_t v, char *buf, int cap) {
    fmt_q32(v, buf, cap, 6);
    uint64_t uv = v < 0 ? -(uint64_t)v : (uint64_t)v;
    if ((uv & 0xFFFFFFFFULL) < 64 || (uv & 0xFFFFFFFFULL) > 0xFFFFFFC0ULL) return;
    int64_t p0 = 0, q0 = 1, p1 = 1, q1 = 0;
    uint64_t num = uv, den = 1ULL << 32;
    for (int step = 0; step < 12 && den != 0; step++) {
        int64_t a = (int64_t)(num / den);
        uint64_t r = num % den;
        int64_t p2 = a * p1 + p0;
        int64_t q2 = a * q1 + q0;
        if (q2 > 240) break;
        p0 = p1; q0 = q1; p1 = p2; q1 = q2;
        if (q1 > 1) {
            int save = q_err;
            q_err = 0;
            q32_t cand = q32_div(q32_from_int(v < 0 ? -p1 : p1), q32_from_int(q1));
            q_err = save;
            q32_t diff = v > cand ? v - cand : cand - v;
            if (diff <= 6) {
                char dec_buf[32];
                fmt_q32(v, dec_buf, sizeof dec_buf, 6);
                int pos = 0;
                if (v < 0) append_char(buf, cap, &pos, '-');
                char tmp[24];
                int tn = 0;
                uint32_t up = (uint32_t)p1;
                while (up > 0) { tmp[tn++] = (char)('0' + (up % 10u)); up /= 10u; }
                while (tn > 0) append_char(buf, cap, &pos, tmp[--tn]);
                append_char(buf, cap, &pos, '/');
                uint32_t uq = (uint32_t)q1;
                while (uq > 0) { tmp[tn++] = (char)('0' + (uq % 10u)); uq /= 10u; }
                while (tn > 0) append_char(buf, cap, &pos, tmp[--tn]);
                append_str(buf, cap, &pos, " (~");
                append_str(buf, cap, &pos, dec_buf);
                append_char(buf, cap, &pos, ')');
                return;
            }
        }
        num = den;
        den = r;
    }
}

// --- RPN Bytecode Compiler & Evaluator ---
enum {
    OP_NUM = 1, OP_VAR_X,
    OP_ADD, OP_SUB, OP_MUL, OP_DIV, OP_POW, OP_NEG,
    OP_SIN, OP_COS, OP_TAN, OP_SQRT, OP_ABS, OP_EXP, OP_LN, OP_LOG10
};

typedef struct {
    uint8_t op;
    q32_t val;
} RpnInst;

typedef struct {
    RpnInst code[96];
    int n;
    int uses_x;
} RpnProg;

static int rpn_prec(int op) {
    if (op >= OP_SIN) return 5;
    if (op == OP_POW) return 4;
    if (op == OP_NEG) return 3;
    if (op == OP_MUL || op == OP_DIV) return 2;
    return 1;
}

static int compile_rpn(const char *s, RpnProg *prog) {
    prog->n = 0;
    prog->uses_x = 0;
    uint8_t op_stk[64];
    int no = 0, operand = 1;

    while (*s) {
        if (is_space(*s)) { s++; continue; }
        if (prog->n >= 92 || no >= 60) return 0;

        // Implicit multiplication when an atom/group just finished (`operand == 0`)
        // and the next token starts an operand (`x`, `pi`, `e`, `(`, or function name, or digit after `)`)
        if (!operand) {
            if (*s == '(' || is_alpha(*s)) {
                int pr = rpn_prec(OP_MUL);
                while (no && op_stk[no - 1] != '(' && rpn_prec(op_stk[no - 1]) >= pr)
                    prog->code[prog->n++] = (RpnInst){op_stk[--no], 0};
                op_stk[no++] = OP_MUL;
                operand = 1;
            }
        }

        if (operand) {
            if (*s == '(') {
                op_stk[no++] = '(';
                s++;
            } else if (*s == '-' || *s == '+') {
                if (*s == '-') op_stk[no++] = OP_NEG;
                s++;
            } else if (is_digit(*s) || (*s == '.' && is_digit(s[1]))) {
                const char *end = s;
                // Do not let parse_q32_num swallow `/denom` inside general expressions!
                uint64_t ip = 0;
                while (is_digit(*end)) {
                    ip = ip * 10 + (uint64_t)(*end - '0');
                    if (ip > 2147483647ULL) return 0;
                    end++;
                }
                uint64_t fp_num = 0, fp_den = 1;
                if (*end == '.') {
                    end++;
                    while (is_digit(*end)) {
                        if (fp_den < 100000000000ULL) {
                            fp_num = fp_num * 10 + (uint64_t)(*end - '0');
                            fp_den *= 10;
                        }
                        end++;
                    }
                }
                uint64_t frac = 0;
                if (fp_num > 0) {
                    uint64_t rem = fp_num;
                    for (int i = 0; i < 32; i++) {
                        rem <<= 1; frac <<= 1;
                        if (rem >= fp_den) { rem -= fp_den; frac |= 1ULL; }
                    }
                    if ((rem << 1) >= fp_den) frac++;
                }
                prog->code[prog->n++] = (RpnInst){OP_NUM, (q32_t)((ip << 32) + frac)};
                s = end;
                operand = 0;
            } else if (is_alpha(*s)) {
                char id[16];
                int il = 0;
                while (is_alpha(*s) || is_digit(*s)) {
                    if (il < 15) id[il++] = (*s >= 'A' && *s <= 'Z') ? (char)(*s + 32) : *s;
                    s++;
                }
                id[il] = 0;
                while (*s == ' ') s++;
                if (streq(id, "x") || streq(id, "t")) {
                    prog->code[prog->n++] = (RpnInst){OP_VAR_X, 0};
                    prog->uses_x = 1;
                    operand = 0;
                } else if (streq(id, "pi")) {
                    prog->code[prog->n++] = (RpnInst){OP_NUM, Q32_PI};
                    operand = 0;
                } else if (streq(id, "e")) {
                    prog->code[prog->n++] = (RpnInst){OP_NUM, Q32_E};
                    operand = 0;
                } else if (*s == '(') {
                    uint8_t fn = 0;
                    if (streq(id, "sin")) fn = OP_SIN;
                    else if (streq(id, "cos")) fn = OP_COS;
                    else if (streq(id, "tan")) fn = OP_TAN;
                    else if (streq(id, "sqrt")) fn = OP_SQRT;
                    else if (streq(id, "abs")) fn = OP_ABS;
                    else if (streq(id, "exp")) fn = OP_EXP;
                    else if (streq(id, "ln") || streq(id, "log")) fn = OP_LN;
                    else if (streq(id, "log10")) fn = OP_LOG10;
                    else return 0;
                    op_stk[no++] = fn;
                    op_stk[no++] = '(';
                    s++;  // consume '('
                } else {
                    return 0;
                }
            } else {
                return 0;
            }
            continue;
        }

        if (*s == ')') {
            while (no && op_stk[no - 1] != '(')
                prog->code[prog->n++] = (RpnInst){op_stk[--no], 0};
            if (!no) return 0;
            no--;  // pop '('
            if (no && op_stk[no - 1] >= OP_SIN)
                prog->code[prog->n++] = (RpnInst){op_stk[--no], 0};
            s++;
            continue;
        }

        uint8_t op = 0;
        int adv = 1;
        if (*s == '+') op = OP_ADD;
        else if (*s == '-') op = OP_SUB;
        else if (*s == '*' && s[1] == '*') { op = OP_POW; adv = 2; }
        else if (*s == '*') op = OP_MUL;
        else if (*s == '/') op = OP_DIV;
        else if (*s == '^') op = OP_POW;
        else return 0;

        int pr = rpn_prec(op);
        while (no && op_stk[no - 1] != '(' &&
               (rpn_prec(op_stk[no - 1]) > pr || (rpn_prec(op_stk[no - 1]) == pr && op != OP_POW)))
            prog->code[prog->n++] = (RpnInst){op_stk[--no], 0};
        op_stk[no++] = op;
        s += adv;
        operand = 1;
    }

    if (operand) return 0;
    while (no) {
        if (op_stk[no - 1] == '(') return 0;
        prog->code[prog->n++] = (RpnInst){op_stk[--no], 0};
    }
    return prog->n > 0;
}

static q32_t eval_rpn(const RpnProg *prog, q32_t x_val, int *ok, char *step_buf, int step_cap) {
    q32_t stk[48];
    int sp = 0, step_pos = 0, step_cnt = 0;
    q_err = 0;
    if (step_buf && step_cap > 0) step_buf[0] = 0;

    for (int i = 0; i < prog->n && !q_err; i++) {
        uint8_t op = prog->code[i].op;
        if (op == OP_NUM) {
            stk[sp++] = prog->code[i].val;
        } else if (op == OP_VAR_X) {
            stk[sp++] = x_val;
        } else if (op == OP_NEG || op >= OP_SIN) {
            if (sp < 1) { q_err = 1; break; }
            q32_t a = stk[sp - 1], r = 0, s_tmp, c_tmp;
            const char *fn_name = NULL;
            if (op == OP_NEG) r = -a;
            else if (op == OP_SIN) { q32_sincos(a, &s_tmp, &c_tmp); r = s_tmp; fn_name = "sin"; }
            else if (op == OP_COS) { q32_sincos(a, &s_tmp, &c_tmp); r = c_tmp; fn_name = "cos"; }
            else if (op == OP_TAN) {
                q32_sincos(a, &s_tmp, &c_tmp);
                if (c_tmp > -4096 && c_tmp < 4096) q_err = 1;
                else r = q32_div(s_tmp, c_tmp);
                fn_name = "tan";
            }
            else if (op == OP_SQRT) { r = q32_sqrt(a); fn_name = "sqrt"; }
            else if (op == OP_ABS) { r = a < 0 ? -a : a; fn_name = "abs"; }
            else if (op == OP_EXP) { r = q32_exp(a); fn_name = "exp"; }
            else if (op == OP_LN) { r = q32_ln(a); fn_name = "ln"; }
            else if (op == OP_LOG10) { r = q32_div(q32_ln(a), Q32_LN10); fn_name = "log10"; }
            stk[sp - 1] = r;
            if (step_buf && fn_name && !q_err && step_cnt < 8) {
                char ba[32], br[32];
                fmt_q32(a, ba, sizeof ba, 6);
                fmt_q32(r, br, sizeof br, 6);
                step_cnt++;
                append_str(step_buf, step_cap, &step_pos, "step ");
                append_char(step_buf, step_cap, &step_pos, (char)('0' + step_cnt));
                append_str(step_buf, step_cap, &step_pos, ": ");
                append_str(step_buf, step_cap, &step_pos, fn_name);
                append_char(step_buf, step_cap, &step_pos, '(');
                append_str(step_buf, step_cap, &step_pos, ba);
                append_str(step_buf, step_cap, &step_pos, ") = ");
                append_str(step_buf, step_cap, &step_pos, br);
                append_char(step_buf, step_cap, &step_pos, '\n');
            }
        } else {
            if (sp < 2) { q_err = 1; break; }
            q32_t b = stk[--sp], a = stk[sp - 1], r = 0;
            const char *op_s = "+";
            if (op == OP_ADD) { r = q32_add(a, b); op_s = "+"; }
            else if (op == OP_SUB) { r = q32_sub(a, b); op_s = "-"; }
            else if (op == OP_MUL) { r = q32_mul(a, b); op_s = "*"; }
            else if (op == OP_DIV) { r = q32_div(a, b); op_s = "/"; }
            else if (op == OP_POW) { r = q32_pow(a, b); op_s = "^"; }
            stk[sp - 1] = r;
            if (step_buf && !q_err && step_cnt < 8) {
                char ba[32], bb[32], br[32];
                fmt_q32(a, ba, sizeof ba, 6);
                fmt_q32(b, bb, sizeof bb, 6);
                fmt_q32(r, br, sizeof br, 6);
                step_cnt++;
                append_str(step_buf, step_cap, &step_pos, "step ");
                append_char(step_buf, step_cap, &step_pos, (char)('0' + step_cnt));
                append_str(step_buf, step_cap, &step_pos, ": ");
                append_str(step_buf, step_cap, &step_pos, ba);
                append_char(step_buf, step_cap, &step_pos, ' ');
                append_str(step_buf, step_cap, &step_pos, op_s);
                append_char(step_buf, step_cap, &step_pos, ' ');
                append_str(step_buf, step_cap, &step_pos, bb);
                append_str(step_buf, step_cap, &step_pos, " = ");
                append_str(step_buf, step_cap, &step_pos, br);
                append_char(step_buf, step_cap, &step_pos, '\n');
            }
        }
    }
    if (q_err || sp != 1) {
        if (ok) *ok = 0;
        return 0;
    }
    if (ok) *ok = 1;
    return stk[0];
}

// Detect exact degree <= 2 polynomial f(x) = a*x^2 + b*x + c using 7-point verification.
// Returns degree (0, 1, 2) or -1 if non-quadratic / nonlinear.
static int detect_poly2(const RpnProg *prog, q32_t *out_a, q32_t *out_b, q32_t *out_c) {
    int ok0, ok1, okm1, ok2, okm2, okt1, okt2;
    q32_t y0 = eval_rpn(prog, 0, &ok0, NULL, 0);
    q32_t y1 = eval_rpn(prog, Q32_ONE, &ok1, NULL, 0);
    q32_t ym1 = eval_rpn(prog, -Q32_ONE, &okm1, NULL, 0);
    if (!ok0 || !ok1 || !okm1) return -1;

    q_err = 0;
    q32_t c = y0;
    q32_t b = q32_sub(y1, ym1) >> 1;
    q32_t a = q32_sub(q32_add(y1, ym1), q32_add(y0, y0)) >> 1;
    if (q_err) return -1;

    q32_t xs[4] = {2 * Q32_ONE, -2 * Q32_ONE, (11LL << 29), -(13LL << 30)};  // 2, -2, 1.375, -3.25
    int *oks[4] = {&ok2, &okm2, &okt1, &okt2};
    for (int i = 0; i < 4; i++) {
        q32_t y = eval_rpn(prog, xs[i], oks[i], NULL, 0);
        if (!*oks[i]) return -1;
        q_err = 0;
        q32_t pred = q32_add(q32_mul(q32_add(q32_mul(a, xs[i]), b), xs[i]), c);
        if (q_err) return -1;
        q32_t diff = y > pred ? y - pred : pred - y;
        if (diff > 256) return -1;  // > 6e-8 mismatch -> higher degree or transcendental
    }
    if (a > -64 && a < 64) a = 0;
    if (b > -64 && b < 64) b = 0;
    if (c > -64 && c < 64) c = 0;
    *out_a = a;
    *out_b = b;
    *out_c = c;
    return a != 0 ? 2 : (b != 0 ? 1 : 0);
}

static q32_t bisect_diff(const RpnProg *lhs, const RpnProg *rhs, q32_t rhs_const,
                         q32_t xa, q32_t xb, q32_t ga) {
    for (int it = 0; it < 28; it++) {
        q32_t xm = (xa >> 1) + (xb >> 1) + ((xa & xb) & 1LL);
        int ok1 = 1, ok2 = 1;
        q32_t lv = eval_rpn(lhs, xm, &ok1, NULL, 0);
        q32_t rv = (rhs && rhs->uses_x) ? eval_rpn(rhs, xm, &ok2, NULL, 0) : rhs_const;
        if (!ok1 || !ok2) break;
        q32_t gm = lv - rv;
        if (gm == 0) return xm;
        if ((ga < 0 && gm < 0) || (ga > 0 && gm > 0)) {
            xa = xm;
            ga = gm;
        } else {
            xb = xm;
        }
    }
    return (xa >> 1) + (xb >> 1);
}

static void add_unique_root(q32_t *roots, int *nroots, int cap, q32_t r) {
    for (int i = 0; i < *nroots; i++) {
        q32_t d = roots[i] > r ? roots[i] - r : r - roots[i];
        if (d < 400000LL) return;  // < 1e-4 duplicate
    }
    if (*nroots < cap) {
        int pos = *nroots;
        while (pos > 0 && roots[pos - 1] > r) {
            roots[pos] = roots[pos - 1];
            pos--;
        }
        roots[pos] = r;
        (*nroots)++;
    }
}

// --- Equation Solver Mode (`LHS = RHS`) ---
static int run_solver_mode(const char *full_eq, const char *lhs_s, const char *rhs_s,
                           int has_bounds, q32_t xmin, q32_t xmax) {
    RpnProg lhs, rhs;
    if (!compile_rpn(lhs_s, &lhs) || !compile_rpn(rhs_s, &rhs)) return 1;

    puts_raw("equation: ");
    puts_raw(full_eq);
    rv_putchar('\n');

    // Build combined difference program g(x) = LHS(x) - RHS(x) if small enough for poly detection
    RpnProg diff = lhs;
    int can_poly = (lhs.n + rhs.n + 1 <= 92);
    if (can_poly) {
        for (int i = 0; i < rhs.n; i++) diff.code[diff.n++] = rhs.code[i];
        diff.code[diff.n++] = (RpnInst){OP_SUB, 0};
        diff.uses_x = lhs.uses_x || rhs.uses_x;
    }

    q32_t pa = 0, pb = 0, pc = 0;
    int deg = (!has_bounds && can_poly) ? detect_poly2(&diff, &pa, &pb, &pc) : -1;

    char buf1[48], buf2[48];
    if (deg == 0) {
        puts_raw(pc == 0 ? "identity: true for all real x\n" : "no solution (contradiction)\n");
        return 0;
    }
    if (deg == 1) {
        q_err = 0;
        q32_t root = q32_div(-pc, pb);
        fmt_q32(pb, buf1, sizeof buf1, 6);
        fmt_q32(pc, buf2, sizeof buf2, 6);
        puts_raw("linear form: ");
        puts_raw(buf1);
        puts_raw("*x + (");
        puts_raw(buf2);
        puts_raw(") = 0\n");
        fmt_q32_rich(root, buf1, sizeof buf1);
        puts_raw("exact solution: x = ");
        puts_raw(buf1);
        rv_putchar('\n');
        int ok1, ok2;
        q32_t lv = eval_rpn(&lhs, root, &ok1, NULL, 0);
        q32_t rv = eval_rpn(&rhs, root, &ok2, NULL, 0);
        if (ok1 && ok2) {
            fmt_q32(lv, buf1, sizeof buf1, 6);
            fmt_q32(rv, buf2, sizeof buf2, 6);
            puts_raw("check: LHS = ");
            puts_raw(buf1);
            puts_raw(", RHS = ");
            puts_raw(buf2);
            rv_putchar('\n');
        }
        return 0;
    }
    if (deg == 2) {
        q_err = 0;
        q32_t disc = q32_sub(q32_mul(pb, pb), 4 * q32_mul(pa, pc));
        char ba[32], bb[32], bc[32], bd[32];
        fmt_q32(pa, ba, sizeof ba, 6);
        fmt_q32(pb, bb, sizeof bb, 6);
        fmt_q32(pc, bc, sizeof bc, 6);
        fmt_q32(disc, bd, sizeof bd, 6);
        puts_raw("quadratic form: (");
        puts_raw(ba);
        puts_raw(")*x^2 + (");
        puts_raw(bb);
        puts_raw(")*x + (");
        puts_raw(bc);
        puts_raw(") = 0\ndiscriminant: D = ");
        puts_raw(bd);
        rv_putchar('\n');
        q32_t two_a = 2 * pa;
        if (disc < -128) {
            q32_t re = q32_div(-pb, two_a);
            q32_t im = q32_div(q32_sqrt(-disc), two_a < 0 ? -two_a : two_a);
            fmt_q32(re, buf1, sizeof buf1, 6);
            fmt_q32(im, buf2, sizeof buf2, 6);
            puts_raw("real solutions: none (D < 0)\ncomplex solutions: x = ");
            puts_raw(buf1);
            puts_raw(" - ");
            puts_raw(buf2);
            puts_raw("*i, x = ");
            puts_raw(buf1);
            puts_raw(" + ");
            puts_raw(buf2);
            puts_raw("*i\n");
            return 0;
        }
        if (disc <= 128) {
            q32_t r = q32_div(-pb, two_a);
            fmt_q32_rich(r, buf1, sizeof buf1);
            puts_raw("double root: x = ");
            puts_raw(buf1);
            rv_putchar('\n');
            return 0;
        }
        q32_t sd = q32_sqrt(disc);
        q32_t r1 = q32_div(-pb - sd, two_a);
        q32_t r2 = q32_div(-pb + sd, two_a);
        if (r1 > r2) { q32_t t = r1; r1 = r2; r2 = t; }
        fmt_q32_rich(r1, buf1, sizeof buf1);
        fmt_q32_rich(r2, buf2, sizeof buf2);
        puts_raw("solutions: x = ");
        puts_raw(buf1);
        puts_raw(", x = ");
        puts_raw(buf2);
        rv_putchar('\n');
        int ok1, ok2;
        q32_t lv = eval_rpn(&lhs, r1, &ok1, NULL, 0);
        q32_t rv = eval_rpn(&rhs, r1, &ok2, NULL, 0);
        if (ok1 && ok2) {
            fmt_q32(lv, buf1, sizeof buf1, 6);
            fmt_q32(rv, buf2, sizeof buf2, 6);
            puts_raw("check (at x1): LHS = ");
            puts_raw(buf1);
            puts_raw(", RHS = ");
            puts_raw(buf2);
            rv_putchar('\n');
        }
        return 0;
    }

    // General nonlinear equation scan + 28-step bisection
    q32_t lo = has_bounds ? xmin : -(20LL << 32);
    q32_t hi = has_bounds ? xmax : (20LL << 32);
    int rhs_const_ok = 1;
    q32_t rhs_const = rhs.uses_x ? 0 : eval_rpn(&rhs, 0, &rhs_const_ok, NULL, 0);
    q32_t roots[8];
    int nroots = 0;
    for (int pass = 0; pass < (has_bounds ? 1 : 2) && nroots == 0; pass++) {
        if (pass == 1) { lo = -(100LL << 32); hi = 100LL << 32; }
        const int N = 64;
        q32_t step = (hi - lo) >> 6;
        q32_t px = lo, pg = 0;
        int pok = 0;
        for (int i = 0; i <= N; i++) {
            q32_t x = (i == N) ? hi : (lo + i * step);
            int ok1, ok2 = rhs_const_ok;
            q32_t lv = eval_rpn(&lhs, x, &ok1, NULL, 0);
            q32_t rv = rhs.uses_x ? eval_rpn(&rhs, x, &ok2, NULL, 0) : rhs_const;
            int ok = ok1 && ok2;
            q32_t g = ok ? (lv - rv) : 0;
            if (ok && g > -256 && g < 256) {
                add_unique_root(roots, &nroots, 8, x);
            } else if (ok && pok && ((pg < 0 && g > 0) || (pg > 0 && g < 0))) {
                q32_t r = bisect_diff(&lhs, &rhs, rhs_const, px, x, pg);
                // Reject vertical asymptotes (e.g. 1/x = 0)
                int okc1, okc2 = rhs_const_ok;
                q32_t cl = eval_rpn(&lhs, r, &okc1, NULL, 0);
                q32_t cr = rhs.uses_x ? eval_rpn(&rhs, r, &okc2, NULL, 0) : rhs_const;
                q32_t rem = (okc1 && okc2) ? (cl > cr ? cl - cr : cr - cl) : Q32_ONE;
                if (okc1 && okc2 && rem < (Q32_ONE >> 10))
                    add_unique_root(roots, &nroots, 8, r);
            }
            px = x; pg = g; pok = ok;
        }
    }
    if (nroots == 0) {
        puts_raw("solutions on search interval: none found\n");
    } else {
        puts_raw("solutions: ");
        for (int i = 0; i < nroots; i++) {
            if (i > 0) puts_raw(", ");
            puts_raw("x = ");
            fmt_q32_rich(roots[i], buf1, sizeof buf1);
            puts_raw(buf1);
        }
        rv_putchar('\n');
        int ok1, ok2;
        q32_t lv = eval_rpn(&lhs, roots[0], &ok1, NULL, 0);
        q32_t rv = eval_rpn(&rhs, roots[0], &ok2, NULL, 0);
        if (ok1 && ok2) {
            fmt_q32(lv, buf1, sizeof buf1, 6);
            fmt_q32(rv, buf2, sizeof buf2, 6);
            puts_raw("check: LHS = ");
            puts_raw(buf1);
            puts_raw(", RHS = ");
            puts_raw(buf2);
            rv_putchar('\n');
        }
    }
    return 0;
}

// Golden-section polish for local extremum on [a, b] (`find_max` = 1 for max, 0 for min)
static q32_t polish_extremum(const RpnProg *prog, q32_t a, q32_t b, int find_max, q32_t *out_y) {
    // inv_phi = (sqrt(5) - 1) / 2 = 0.61803398875 in Q32.32 = 2654435769LL
    const q32_t INV_PHI = 2654435769LL;
    q32_t c = b - q32_mul(b - a, INV_PHI);
    q32_t d = a + q32_mul(b - a, INV_PHI);
    int okc, okd;
    q32_t yc = eval_rpn(prog, c, &okc, NULL, 0);
    q32_t yd = eval_rpn(prog, d, &okd, NULL, 0);
    for (int it = 0; it < 15 && okc && okd; it++) {
        int pick_left = find_max ? (yc > yd) : (yc < yd);
        if (pick_left) {
            b = d; d = c; yd = yc;
            c = b - q32_mul(b - a, INV_PHI);
            yc = eval_rpn(prog, c, &okc, NULL, 0);
        } else {
            a = c; c = d; yc = yd;
            d = a + q32_mul(b - a, INV_PHI);
            yd = eval_rpn(prog, d, &okd, NULL, 0);
        }
    }
    q32_t xm = (a >> 1) + (b >> 1);
    int okm, okp, okn;
    q32_t ym = eval_rpn(prog, xm, &okm, NULL, 0);
    // Final 3-point parabolic vertex polish (`dh = 2^-11`) for exact 6-decimal extrema
    q32_t dh = (1LL << 21);
    q32_t yp = eval_rpn(prog, xm + dh, &okp, NULL, 0);
    q32_t yn = eval_rpn(prog, xm - dh, &okn, NULL, 0);
    if (okm && okp && okn) {
        q32_t denom = 2 * ((yp - 2 * ym) + yn);
        if (denom > 1024 || denom < -1024) {
            int save = q_err;
            q_err = 0;
            q32_t delta = q32_mul(dh, q32_div(yn - yp, denom));
            if (!q_err && delta > -dh && delta < dh) {
                xm += delta;
                ym = eval_rpn(prog, xm, &okm, NULL, 0);
            }
            q_err = save;
        }
    }
    *out_y = ym;
    return xm;
}

// --- Graphing Calculator & Table Mode (`y = f(x)`) ---
static int run_graph_mode(const char *expr_s, const RpnProg *prog, q32_t xmin, q32_t xmax, int table_only) {
    if (xmax <= xmin) xmax = xmin + (10LL << 32);

    // 48 subintervals -> 49 grid points:
    //   48 = 2 * 24 (exact 25 plot columns at i = 2*col)
    //   48 = 6 * 8  (exact 7 table points at i = 8*k)
    //   48 is even  (exact 48-panel Simpson's rule definite integral!)
    static q32_t gx[49], gy[49];
    static uint8_t gok[49];
    q32_t span = xmax - xmin;
    q32_t h = q32_div(span, 48 * Q32_ONE);

    q32_t ymin = Q32_MAX, ymax = Q32_MIN;
    int valid_cnt = 0, all_valid = 1;
    for (int i = 0; i <= 48; i++) {
        gx[i] = (i == 48) ? xmax : (xmin + i * h);
        int ok = 0;
        gy[i] = eval_rpn(prog, gx[i], &ok, NULL, 0);
        gok[i] = (uint8_t)ok;
        if (ok) {
            if (gy[i] < ymin) ymin = gy[i];
            if (gy[i] > ymax) ymax = gy[i];
            valid_cnt++;
        } else {
            all_valid = 0;
        }
    }

    char b1[40], b2[40], b3[40];
    fmt_q32(xmin, b1, sizeof b1, 4);
    fmt_q32(xmax, b2, sizeof b2, 4);
    puts_raw("function: f(x) = ");
    puts_raw(expr_s);
    puts_raw("  on [");
    puts_raw(b1);
    puts_raw(", ");
    puts_raw(b2);
    puts_raw("]\n");

    if (valid_cnt == 0) {
        puts_raw("error: function undefined on entire interval\n");
        return 0;
    }

    // 1. Roots on [xmin, xmax]
    q32_t roots[8];
    int nroots = 0;
    q32_t pa = 0, pb = 0, pc = 0;
    int deg = detect_poly2(prog, &pa, &pb, &pc);
    if (deg == 1) {
        q32_t r = q32_div(-pc, pb);
        if (r >= xmin - 1024 && r <= xmax + 1024) add_unique_root(roots, &nroots, 8, r);
    } else if (deg == 2) {
        q32_t disc = q32_sub(q32_mul(pb, pb), 4 * q32_mul(pa, pc));
        if (disc >= -128) {
            q32_t sd = disc > 0 ? q32_sqrt(disc) : 0;
            q32_t r1 = q32_div(-pb - sd, 2 * pa);
            q32_t r2 = q32_div(-pb + sd, 2 * pa);
            if (r1 >= xmin - 1024 && r1 <= xmax + 1024) add_unique_root(roots, &nroots, 8, r1);
            if (r2 >= xmin - 1024 && r2 <= xmax + 1024) add_unique_root(roots, &nroots, 8, r2);
        }
    } else {
        for (int i = 0; i <= 48; i++) {
            if (gok[i] && gy[i] > -512 && gy[i] < 512) {
                add_unique_root(roots, &nroots, 8, gx[i]);
            } else if (i < 48 && gok[i] && gok[i + 1] &&
                       ((gy[i] < 0 && gy[i + 1] > 0) || (gy[i] > 0 && gy[i + 1] < 0))) {
                q32_t r = bisect_diff(prog, NULL, 0, gx[i], gx[i + 1], gy[i]);
                int okr;
                q32_t yr = eval_rpn(prog, r, &okr, NULL, 0);
                q32_t ay = yr < 0 ? -yr : yr;
                if (okr && ay < (Q32_ONE >> 10))
                    add_unique_root(roots, &nroots, 8, r);
            }
        }
    }
    puts_raw("roots: ");
    if (nroots == 0) {
        puts_raw("none on interval\n");
    } else {
        for (int i = 0; i < nroots; i++) {
            if (i > 0) puts_raw(", ");
            puts_raw("x = ");
            fmt_q32_rich(roots[i], b1, sizeof b1);
            puts_raw(b1);
        }
        rv_putchar('\n');
    }

    // 2. Vertex / Local Extrema
    if (deg == 2) {
        q32_t vx = q32_div(-pb, 2 * pa);
        int okv;
        q32_t vy = eval_rpn(prog, vx, &okv, NULL, 0);
        if (okv && vy < ymin && vx >= xmin && vx <= xmax) ymin = vy;
        if (okv && vy > ymax && vx >= xmin && vx <= xmax) ymax = vy;
        fmt_q32(vx, b1, sizeof b1, 6);
        fmt_q32(vy, b2, sizeof b2, 6);
        puts_raw(pa > 0 ? "vertex (min): (" : "vertex (max): (");
        puts_raw(b1);
        puts_raw(", ");
        puts_raw(b2);
        puts_raw(")\n");
    } else {
        int printed_ext = 0;
        for (int i = 1; i < 48 && printed_ext < 3; i++) {
            if (!gok[i - 1] || !gok[i] || !gok[i + 1]) continue;
            int is_max = (gy[i] > gy[i - 1] && gy[i] >= gy[i + 1] && (gy[i] > gy[i - 1] || gy[i] > gy[i + 1]));
            int is_min = (gy[i] < gy[i - 1] && gy[i] <= gy[i + 1] && (gy[i] < gy[i - 1] || gy[i] < gy[i + 1]));
            if (is_max || is_min) {
                q32_t ey = 0;
                q32_t ex = polish_extremum(prog, gx[i - 1], gx[i + 1], is_max, &ey);
                if (ey < ymin) ymin = ey;
                if (ey > ymax) ymax = ey;
                if (printed_ext > 0) puts_raw(" | ");
                puts_raw(is_max ? "local max: (" : "local min: (");
                fmt_q32(ex, b1, sizeof b1, 6);
                fmt_q32(ey, b2, sizeof b2, 6);
                puts_raw(b1);
                puts_raw(", ");
                puts_raw(b2);
                rv_putchar(')');
                printed_ext++;
            }
        }
        if (printed_ext > 0) rv_putchar('\n');
    }

    // 3. y-intercept & 96-panel Simpson's definite integral
    int has_yint = 0;
    if (xmin <= 0 && xmax >= 0) {
        int ok0;
        q32_t y0 = eval_rpn(prog, 0, &ok0, NULL, 0);
        if (ok0) {
            fmt_q32(y0, b1, sizeof b1, 6);
            puts_raw("y-intercept: (0, ");
            puts_raw(b1);
            rv_putchar(')');
            has_yint = 1;
        }
    }
    if (all_valid) {
        q_err = 0;
        q32_t sum = q32_add(gy[0], gy[48]);
        for (int i = 1; i < 48 && !q_err; i++) {
            q32_t term = (i & 1) ? q32_mul(4 * Q32_ONE, gy[i]) : q32_mul(2 * Q32_ONE, gy[i]);
            sum = q32_add(sum, term);
        }
        q32_t integ = q32_mul(q32_div(h, 3 * Q32_ONE), sum);
        if (!q_err) {
            if (has_yint) puts_raw(" | ");
            puts_raw("integral: ");
            fmt_q32(integ, b1, sizeof b1, 6);
            puts_raw(b1);
            has_yint = 1;
        }
    }
    if (has_yint) rv_putchar('\n');

    // 4. Coordinate Table (7 clean points; with derivative f'(x) when in table mode)
    if (table_only) {
        puts_raw("x | f(x) | f'(x)\n");
        q32_t dh = (1LL << 20);  // 2^-12 (~2.44e-4) central difference step
        for (int k = 0; k <= 6; k++) {
            int idx = 8 * k;
            fmt_q32(gx[idx], b1, sizeof b1, 4);
            puts_raw(b1);
            puts_raw(" | ");
            if (!gok[idx]) {
                puts_raw("undef | undef\n");
            } else {
                fmt_q32(gy[idx], b2, sizeof b2, 6);
                puts_raw(b2);
                puts_raw(" | ");
                int okp, okm;
                q32_t yp = eval_rpn(prog, gx[idx] + dh, &okp, NULL, 0);
                q32_t ym = eval_rpn(prog, gx[idx] - dh, &okm, NULL, 0);
                if (okp && okm) {
                    q32_t deriv = (yp - ym) << 11;  // / (2 * 2^-12) = * 2^11
                    fmt_q32(deriv, b3, sizeof b3, 6);
                    puts_raw(b3);
                } else {
                    puts_raw("undef");
                }
                rv_putchar('\n');
            }
        }
        return 0;
    }

    puts_raw("table: ");
    for (int k = 0; k <= 6; k++) {
        int idx = 8 * k;
        if (k > 0) puts_raw(", ");
        rv_putchar('(');
        fmt_q32(gx[idx], b1, sizeof b1, 4);
        puts_raw(b1);
        puts_raw(", ");
        if (gok[idx]) {
            fmt_q32(gy[idx], b2, sizeof b2, 4);
            puts_raw(b2);
        } else {
            puts_raw("undef");
        }
        rv_putchar(')');
    }
    rv_putchar('\n');

    // 5. 9x25 ASCII Coordinate Plot
    q32_t pmin = ymin, pmax = ymax;
    // Clamp extreme asymptote spikes so e.g. 1/x or tan(x) remains readable
    if (pmax - pmin > (200LL << 32) && !all_valid) {
        pmin = -(20LL << 32);
        pmax = 20LL << 32;
    }
    if (pmax - pmin < (1LL << 28)) {
        pmin -= Q32_ONE;
        pmax += Q32_ONE;
    }
    q32_t yspan = pmax - pmin;

    char grid[9][25];
    for (int r = 0; r < 9; r++)
        for (int c = 0; c < 25; c++) grid[r][c] = ' ';

    int rx = -1, cy = -1;
    q32_t inv_yspan8 = q32_div(8 * Q32_ONE, yspan);
    if (pmin <= 0 && pmax >= 0) {
        rx = (int)((q32_mul(pmax, inv_yspan8) + (1LL << 31)) >> 32);
        if (rx >= 0 && rx < 9)
            for (int c = 0; c < 25; c++) grid[rx][c] = '-';
    }
    if (xmin <= 0 && xmax >= 0) {
        cy = (int)((q32_div(-xmin, span) * 24 + (1LL << 31)) >> 32);
        if (cy >= 0 && cy < 25)
            for (int r = 0; r < 9; r++) grid[r][cy] = '|';
    }
    if (rx >= 0 && rx < 9 && cy >= 0 && cy < 25) grid[rx][cy] = '+';

    for (int c = 0; c < 25; c++) {
        int idx = 2 * c;
        if (!gok[idx] || gy[idx] < pmin || gy[idx] > pmax) continue;
        int r = (int)((q32_mul(pmax - gy[idx], inv_yspan8) + (1LL << 31)) >> 32);
        if (r >= 0 && r < 9) grid[r][c] = '*';
    }

    for (int r = 0; r < 9; r++) {
        if ((r & 1) == 0 || r == rx) {
            q32_t ytick = (r == rx) ? 0 : (pmax - (yspan * r) / 8);
            fmt_q32(ytick, b1, sizeof b1, 2);
            int l = slen(b1);
            putc_rep(' ', l < 7 ? 7 - l : 0);
            puts_raw(b1);
        } else {
            putc_rep(' ', 7);
        }
        puts_raw(" |");
        int endc = 25;
        while (endc > 0 && grid[r][endc - 1] == ' ') endc--;
        for (int c = 0; c < endc; c++) rv_putchar(grid[r][c]);
        rv_putchar('\n');
    }
    // Bottom x-axis tick line
    fmt_q32(xmin, b1, sizeof b1, 2);
    fmt_q32((xmin >> 1) + (xmax >> 1), b2, sizeof b2, 2);
    fmt_q32(xmax, b3, sizeof b3, 2);
    puts_raw("        x=");
    puts_raw(b1);
    int used_w = 2 + slen(b1);
    int mid_target = 11 - (2 + slen(b2)) / 2;
    if (mid_target <= used_w) mid_target = used_w + 1;
    putc_rep(' ', mid_target - used_w);
    puts_raw("x=");
    puts_raw(b2);
    used_w = mid_target + 2 + slen(b2);
    int right_target = 25 - (2 + slen(b3));
    if (right_target <= used_w) right_target = used_w + 1;
    putc_rep(' ', right_target - used_w);
    puts_raw("x=");
    puts_raw(b3);
    rv_putchar('\n');
    return 0;
}

// ============================================================================
// Tool-Call Gathering, UTF-8 Normalization & Mode Dispatcher
// ============================================================================

static char raw_call[1024];
static char norm_call[1024];
static char fn_buf[32];
static char expr_buf[256];
static char xmin_buf[64];
static char xmax_buf[64];
static int has_xmin, has_xmax;

static void normalize_utf8(const char *src, char *dst, int cap) {
    int o = 0;
    for (int i = 0; src[i] && o < cap - 4; ) {
        unsigned char c0 = (unsigned char)src[i];
        unsigned char c1 = (unsigned char)src[i + 1];
        unsigned char c2 = c1 ? (unsigned char)src[i + 2] : 0;
        if (c0 == 0xC3 && c1 == 0x97) { dst[o++] = '*'; i += 2; }                      // ×
        else if (c0 == 0xC3 && c1 == 0xB7) { dst[o++] = '/'; i += 2; }                 // ÷
        else if (c0 == 0xC2 && c1 == 0xB2) { dst[o++] = '^'; dst[o++] = '2'; i += 2; } // ²
        else if (c0 == 0xC2 && c1 == 0xB3) { dst[o++] = '^'; dst[o++] = '3'; i += 2; } // ³
        else if (c0 == 0xCF && c1 == 0x80) { dst[o++] = 'p'; dst[o++] = 'i'; i += 2; } // π
        else if (c0 == 0xE2 && c1 == 0x88 && c2 == 0x92) { dst[o++] = '-'; i += 3; }   // −
        else if (c0 == 0xE2 && c1 == 0x80 && c2 == 0x93) { dst[o++] = '-'; i += 3; }   // –
        else if (c0 == '&' && starts_with_ci(src + i, "&lt;")) { dst[o++] = '<'; i += 4; }
        else if (c0 == '&' && starts_with_ci(src + i, "&gt;")) { dst[o++] = '>'; i += 4; }
        else { dst[o++] = src[i++]; }
    }
    dst[o] = 0;
}

static void parse_call_payload(void) {
    normalize_utf8(raw_call, norm_call, sizeof norm_call);
    fn_buf[0] = 0;
    expr_buf[0] = 0;
    xmin_buf[0] = 0;
    xmax_buf[0] = 0;
    has_xmin = 0;
    has_xmax = 0;

    int saw_param = 0;
    for (int i = 0; norm_call[i]; i++) {
        if (starts_with_ci(norm_call + i, "<function=")) {
            int s = i + 10, e = s;
            while (norm_call[e] && norm_call[e] != '>' && !is_space(norm_call[e])) e++;
            copy_trim(fn_buf, sizeof fn_buf, norm_call + s, e - s);
        } else if (starts_with_ci(norm_call + i, "<parameter=")) {
            saw_param = 1;
            int ks = i + 11, ke = ks;
            while (norm_call[ke] && norm_call[ke] != '>') ke++;
            if (!norm_call[ke]) break;
            char key[24];
            copy_trim(key, sizeof key, norm_call + ks, ke - ks);
            for (int k = 0; key[k]; k++) if (key[k] >= 'A' && key[k] <= 'Z') key[k] = (char)(key[k] + 32);
            int vs = ke + 1, ve = vs;
            while (norm_call[ve] && norm_call[ve] != '<') ve++;
            if (starts_with_ci(key, "expr") || starts_with_ci(key, "eq") ||
                starts_with_ci(key, "func") || starts_with_ci(key, "cmd") || streq(key, "f")) {
                copy_trim(expr_buf, sizeof expr_buf, norm_call + vs, ve - vs);
            } else if (streq(key, "xmin") || streq(key, "min") || streq(key, "a") || streq(key, "from")) {
                copy_trim(xmin_buf, sizeof xmin_buf, norm_call + vs, ve - vs);
                has_xmin = xmin_buf[0] != 0;
            } else if (streq(key, "xmax") || streq(key, "max") || streq(key, "b") || streq(key, "to")) {
                copy_trim(xmax_buf, sizeof xmax_buf, norm_call + vs, ve - vs);
                has_xmax = xmax_buf[0] != 0;
            } else if (streq(key, "mode")) {
                copy_trim(fn_buf, sizeof fn_buf, norm_call + vs, ve - vs);
            } else if (!expr_buf[0]) {
                copy_trim(expr_buf, sizeof expr_buf, norm_call + vs, ve - vs);
            }
            i = ve - 1;
        }
    }

    if (!saw_param) {
        // Fallback: JSON {"expression": "..."} or plain text expression
        int found_json = 0;
        for (int i = 0; norm_call[i]; i++) {
            if (starts_with_ci(norm_call + i, "\"expression\"") || starts_with_ci(norm_call + i, "\"equation\"")) {
                int j = i + 12;
                while (norm_call[j] && norm_call[j] != ':') j++;
                if (norm_call[j] == ':') j++;
                while (is_space(norm_call[j])) j++;
                if (norm_call[j] == '"') {
                    int vs = ++j;
                    while (norm_call[j] && norm_call[j] != '"') j++;
                    copy_trim(expr_buf, sizeof expr_buf, norm_call + vs, j - vs);
                    found_json = 1;
                    break;
                }
            }
        }
        if (!found_json) copy_trim(expr_buf, sizeof expr_buf, norm_call, slen(norm_call));
    }

    // Strip optional inline domain suffix: " on [a, b]" or " from x = a to b" or " from a to b"
    int n = slen(expr_buf);
    for (int i = 0; i < n; i++) {
        if (!has_xmin && (starts_with_ci(expr_buf + i, " on [") || starts_with_ci(expr_buf + i, " in ["))) {
            int lb = i + 5, cm = lb;
            while (expr_buf[cm] && expr_buf[cm] != ',' && expr_buf[cm] != ']') cm++;
            if (expr_buf[cm] == ',') {
                int rb = cm + 1;
                while (expr_buf[rb] && expr_buf[rb] != ']') rb++;
                if (expr_buf[rb] == ']') {
                    copy_trim(xmin_buf, sizeof xmin_buf, expr_buf + lb, cm - lb);
                    copy_trim(xmax_buf, sizeof xmax_buf, expr_buf + cm + 1, rb - cm - 1);
                    has_xmin = has_xmax = 1;
                    expr_buf[i] = 0;
                    copy_trim(expr_buf, sizeof expr_buf, expr_buf, i);
                    break;
                }
            }
        } else if (!has_xmin && starts_with_ci(expr_buf + i, " from ")) {
            int p = i + 6;
            while (is_space(expr_buf[p])) p++;
            if ((expr_buf[p] == 'x' || expr_buf[p] == 'X') && (expr_buf[p + 1] == '=' || expr_buf[p + 1] == ' ')) {
                p++;
                while (is_space(expr_buf[p]) || expr_buf[p] == '=') p++;
            }
            int to_pos = -1;
            for (int k = p; expr_buf[k]; k++) {
                if (starts_with_ci(expr_buf + k, " to ")) { to_pos = k; break; }
            }
            if (to_pos > p) {
                copy_trim(xmin_buf, sizeof xmin_buf, expr_buf + p, to_pos - p);
                copy_trim(xmax_buf, sizeof xmax_buf, expr_buf + to_pos + 4, slen(expr_buf) - (to_pos + 4));
                has_xmin = has_xmax = 1;
                expr_buf[i] = 0;
                copy_trim(expr_buf, sizeof expr_buf, expr_buf, i);
                break;
            }
        }
    }

    // Optional command prefix ("graph ...", "plot ...", "solve ...", "table ...")
    const char *prefixes[4] = {"graph ", "plot ", "solve ", "table "};
    for (int k = 0; k < 4; k++) {
        if (starts_with_ci(expr_buf, prefixes[k])) {
            int plen = slen(prefixes[k]);
            if (!fn_buf[0] || streq(fn_buf, "graphcalc"))
                copy_trim(fn_buf, sizeof fn_buf, prefixes[k], plen - 1);
            copy_trim(expr_buf, sizeof expr_buf, expr_buf + plen, slen(expr_buf) - plen);
            break;
        }
    }

    // Strip leading "y =" or "f(x) =" when there is no second '='
    int eq_cnt = 0, first_eq = -1;
    for (int i = 0; expr_buf[i]; i++) {
        if (expr_buf[i] == '=') {
            if (eq_cnt == 0) first_eq = i;
            eq_cnt++;
        }
    }
    if (eq_cnt == 1 && first_eq > 0) {
        char lhs_tmp[32];
        copy_trim(lhs_tmp, sizeof lhs_tmp, expr_buf, first_eq);
        if (streq(lhs_tmp, "y") || streq(lhs_tmp, "Y") ||
            starts_with_ci(lhs_tmp, "f(x)") || starts_with_ci(lhs_tmp, "g(x)")) {
            copy_trim(expr_buf, sizeof expr_buf, expr_buf + first_eq + 1, slen(expr_buf) - first_eq - 1);
        }
    }
}

static int has_alpha_outside_gcd(const char *s) {
    for (int i = 0; s[i]; ) {
        if (starts_with_ci(s + i, "gcd(")) { i += 4; continue; }
        if (is_alpha(s[i])) return 1;
        i++;
    }
    return 0;
}

static void dispatch_call(void) {
    parse_call_payload();
    if (!expr_buf[0]) {
        puts_raw("error: empty expression\n");
        return;
    }

    q32_t xmin = -(5LL << 32), xmax = 5LL << 32;
    int got_bounds = 0;
    if (has_xmin && parse_q32_num(xmin_buf, NULL, &xmin)) got_bounds = 1;
    if (has_xmax && parse_q32_num(xmax_buf, NULL, &xmax)) got_bounds = 1;

    // 1. Check for Equation (`LHS = RHS`)
    int eq_pos = -1;
    for (int i = 0; expr_buf[i]; i++) {
        if (expr_buf[i] == '=') { eq_pos = i; break; }
    }
    if (eq_pos > 0 && expr_buf[eq_pos + 1]) {
        char lhs_s[160], rhs_s[160];
        copy_trim(lhs_s, sizeof lhs_s, expr_buf, eq_pos);
        copy_trim(rhs_s, sizeof rhs_s, expr_buf + eq_pos + 1, slen(expr_buf) - eq_pos - 1);
        if (run_solver_mode(expr_buf, lhs_s, rhs_s, got_bounds, xmin, xmax) == 0) return;
        puts_raw("error: invalid equation syntax\n");
        return;
    }

    // 2. Pure Bignum / Exact Rational Arithmetic (no variable x, no trig/log/sqrt)
    int is_explicit_graph_fn = streq_ci(fn_buf, "graph") || streq_ci(fn_buf, "plot") || streq_ci(fn_buf, "table");
    if (!has_alpha_outside_gcd(expr_buf) && !got_bounds && !is_explicit_graph_fn) {
        if (run_bignum_engine(expr_buf) == 0) return;
    }

    // 3. Compile RPN for Scientific Numeric Eval or Function Graphing / Table
    RpnProg prog;
    if (!compile_rpn(expr_buf, &prog)) {
        puts_raw("error: syntax error in expression\n");
        return;
    }

    int force_graph = got_bounds || is_explicit_graph_fn;
    if (!prog.uses_x && !force_graph) {
        char step_log[512];
        int ok = 0;
        q32_t val = eval_rpn(&prog, 0, &ok, step_log, sizeof step_log);
        if (!ok) {
            puts_raw("error: math domain or overflow error\n");
            return;
        }
        puts_raw("expression: ");
        puts_raw(expr_buf);
        rv_putchar('\n');
        if (step_log[0]) puts_raw(step_log);
        char vbuf[48];
        fmt_q32_rich(val, vbuf, sizeof vbuf);
        puts_raw("value: ");
        puts_raw(vbuf);
        rv_putchar('\n');
        return;
    }

    int table_only = starts_with_ci(fn_buf, "table");
    run_graph_mode(expr_buf, &prog, xmin, xmax, table_only);
}

#ifndef QEMU
// P3 Bridge memory & CarbonText extended special tokens (256..264)
#define TEXT_BUF ((const volatile unsigned char *)0x8000)
#define CAP_ROWS 80

enum {
    TOK_TOOL_RESPONSE = 256,
    TOK_TOOL_RESPONSE_END = 257,
    TOK_THINK = 258,
    TOK_THINK_END = 259,
    TOK_IM_END = 260,
    TOK_IM_START = 261,
    TOK_USER = 262,
    TOK_ASSISTANT = 263,
    TOK_NL2 = 264
};

static void put_tok(int c) { rv_ecall2(1, c, 0, 0); }

static int gather_bridge_text(void) {
    int n = 0;
    for (int i = 1; i <= CAP_ROWS; i++) {
        const volatile unsigned char *r = TEXT_BUF + 16 * i;
        if (!r[0]) break;
        for (int k = 0; k < 16 && r[k]; k++) {
            if (n < (int)sizeof(raw_call) - 1) raw_call[n++] = (char)r[k];
        }
    }
    raw_call[n] = 0;
    return n;
}

int main(void) {
    for (;;) {
        // Park until Qwen emits </tool_call> (GO resumes right after this exit ecall)
        rv_ecall2(93, 0, 0, 0);
        int n = gather_bridge_text();

        // Native Qwen3.5 ChatML tool-response turn header:
        // <|im_end|>\n<|im_start|>user\n<tool_response>\n
        put_tok(TOK_IM_END);
        put_tok('\n');
        put_tok(TOK_IM_START);
        put_tok(TOK_USER);
        put_tok('\n');
        put_tok(TOK_TOOL_RESPONSE);
        put_tok('\n');

        if (n <= 0) {
            puts_raw("error: empty tool call\n");
        } else {
            dispatch_call();
        }

        // Native Qwen3.5 ChatML assistant continuation header:
        // </tool_response><|im_end|>\n<|im_start|>assistant\n<think>\n\n</think>\n\n##[HALT -> ' ']
        put_tok(TOK_TOOL_RESPONSE_END);
        put_tok(TOK_IM_END);
        put_tok('\n');
        put_tok(TOK_IM_START);
        put_tok(TOK_ASSISTANT);
        put_tok('\n');
        put_tok(TOK_THINK);
        put_tok(TOK_NL2);
        put_tok(TOK_THINK_END);
        put_tok(TOK_NL2);
        put_tok('#');
        put_tok('#');
    }
}
#else
int main(void) {
    int pos = 0;
    while (pos < (int)sizeof(raw_call) - 1) {
        long got = rv_ecall2(63, 0, (long)(raw_call + pos), (long)(sizeof(raw_call) - 1 - pos));
        if (got <= 0) break;
        pos += (int)got;
    }
    raw_call[pos] = 0;
    dispatch_call();
    return 0;
}
#endif