compressor.c

/*
 * PMC v4: Bit-Level Context-Mixing Compressor
 *
 * PAQ-inspired architecture: multiple context models combined via adaptive
 * logistic mixing, APM correction, match model, BCJ filter, binary rANS.
 *
 * v4.8 changes vs v4.7:
 *   - Mixer learning rate fine-tuning: 24/24 → 28/28
 *   - Word-aware SSE (ISSE5) rate: 5 → 4 (faster adaptation for text)
 *     (~0.27% additional improvement across full Silesia corpus + specialized files)
 *
 * v4.7 changes vs v4.6:
 *   - Mixer learning rate tuning: divisors 10/12 → 24/24 for both mixers
 *     (~2.4% improvement across full Silesia corpus, zero regressions on files > 256 KB)
 *
 * v4.4 changes vs v4.3:
 *   - Word trigram context model (prev_prev_word + prev_word + position)
 *   - MTF word cache pre-filter for text files (replaces BCJ)
 *
 * v4.3 changes vs v4.2:
 *   - Block-level text detection with model gating (gates 36 noise models on text)
 *   - Previous-word context model for word-bigram prediction
 *
 * v4 changes vs v3:
 *   - 34 sparse context models using non-adjacent byte patterns (~90 KB gain)
 *   - Tuned mixer learning rates (10/12 vs 8/8) (~1.6 KB gain)
 *
 * Sparse models are the dominant source of compression gain. They capture
 * correlations between non-adjacent bytes (e.g. structured fields at regular
 * offsets in ELF/binary data). Each sparse model costs ~8 MB of memory
 * (4M tagged entries x 2 bytes) and yields diminishing returns:
 *   - First 10 models: ~3-7 KB each
 *   - Models 10-20:    ~1-3 KB each
 *   - Models 20-34:    ~250-600 bytes each
 *
 * To squeeze further when memory is not a concern, more sparse models can be
 * added by incrementing NUM_CM/NUM_MODELS, defining additional ctx[] entries
 * with new byte-skip patterns, and adding corresponding cm_init() calls.
 * Promising patterns to try: longer strides (pb[16]..pb[31] — requires
 * expanding the pb[] history buffer), higher-order sparse combinations
 * (4-5 non-adjacent bytes), or data-type-specific patterns. Each new model
 * adds ~8 MB memory for ~200-500 bytes of compression gain at this point.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <math.h>

/* ---- Configuration ---- */
#define BLOCK_SIZE      (1 << 23)   /* 8 MB */
#define PROB_BITS       12
#define PROB_SCALE      (1 << PROB_BITS)  /* 4096 */
#define PROB_HALF       (PROB_SCALE / 2)
#define RANS_L          (1u << 23)
#define MAX_BIT_EVENTS  (BLOCK_SIZE * 12)

/*
 * Context model table sizes.
 *
 * 9 order-N models (o0..o4, o6, o8, o12, o16), 1 word model, 1 prev_word
 * model, and NUM_SPARSE sparse models that use non-adjacent byte combinations.
 *
 * To add more sparse models for better compression (at ~8 MB each):
 *   1. Increase NUM_SPARSE
 *   2. Add the new ctx[] hash in process_block()
 *   3. Add the corresponding cm_init() call in do_compress/do_decompress
 * NUM_CM and NUM_MODELS are computed automatically.
 */
#define NUM_SPARSE      34
#define NUM_CM          (12 + NUM_SPARSE) /* 11 fixed (o0-o4, o6, o8, o12, word, prev_word, word_trigram) + sparseN + o16 */
#define NUM_MODELS      (NUM_CM + 6)      /* CM + match + ICM + delta + stride + wpred + lzp */
#define CM_O0_BITS      8     /* 256 bytes */
#define CM_O1_BITS      20
#define CM_O2_BITS      22
#define CM_O3_BITS      23
#define CM_O4_BITS      24
#define CM_O6_BITS      24
#define CM_O8_BITS      24
#define CM_O12_BITS     24    /* 16M tagged entries */
#define CM_WORD_BITS    22    /* 4M tagged entries */
#define CM_SPARSE_BITS  22    /* 4M tagged entries — used for all sparse models */
#define CM_O16_BITS     24    /* 16M tagged entries */

/* ICM (Indirect Context Model) */
#define ICM_BITS        24    /* 16M entries */
#define ICM_TAB_SIZE    (1u << ICM_BITS)
#define ICM_NCTR        256   /* shared adaptive counters */
#define ICM_RATE        5     /* counter += (target - counter) >> 5 */

/* Delta model */
#define DELTA_BITS      22    /* 4M entries */
#define DELTA_TAB_SIZE  (1u << DELTA_BITS)

/* Stride/Record model */
#define STRIDE_BITS     22    /* 4M entries */
#define STRIDE_TAB_SIZE (1u << STRIDE_BITS)
#define MAX_STRIDE      128

/* Match model */
#define MATCH_BITS      22   /* 4M entries */
#define MATCH_TAB_SIZE  (1u << MATCH_BITS)

/* LZP (Lempel-Ziv Prediction) model */
#define LZP_BITS        22    /* 4M entries */
#define LZP_SIZE        (1u << LZP_BITS)

/* Mixer 1: bp + prev_byte + prev2_byte + text + match => up to 4096 slots */
#define MIXER_SLOTS     4096
/* Mixer 2: bit_ctx */
#define MIXER2_SLOTS    256

/* ISSE chain: 4 refinement stages */
#define ISSE1_CTX_BITS  8    /* 256 ctx */
#define ISSE1_PROB_BITS 9    /* 512 bins */
#define ISSE2_CTX_BITS  8    /* 256 ctx */
#define ISSE2_PROB_BITS 9    /* 512 bins */
#define ISSE3_CTX_BITS  7    /* 128 ctx */
#define ISSE3_PROB_BITS 9    /* 512 bins */
#define ISSE4_CTX_BITS  10   /* 1024 ctx (with match) */
#define ISSE4_PROB_BITS 6    /* 64 bins */
#define ISSE5_CTX_BITS  18   /* 262144 ctx (word + bp) */
#define ISSE5_PROB_BITS 5    /* 32 bins */

/* Stretch/squash table range */
#define LOGIT_SCALE     256
#define LOGIT_MIN       (-2048)
#define LOGIT_MAX       2047
#define LOGIT_RANGE     (LOGIT_MAX - LOGIT_MIN + 1)

/* Interleaved rANS */
#define RANS_INTERLEAVE 4

static const uint8_t MAGIC[4] = {'P','M','C','4'};
#define FLAG_BCJ       1
#define FLAG_IMG_DELTA 2
#define FLAG_WAV_DELTA 4
#define FLAG_TEXT_MTF  8
#define FLAG_TEXT      16   /* file is text — enables text gating in blocks */
#define FLAG_IMG_RAW   32   /* raw image detected — 2D delta filter applied */

/* ---- Mixer learning rate divisors (tunable via -DMIX1_LR=N -DMIX2_LR=N) ---- */
#ifndef MIX1_LR
#define MIX1_LR 28
#endif
#ifndef MIX2_LR
#define MIX2_LR 28
#endif

/* ---- ISSE/APM learning rate shifts (tunable via -DISSE1_LR=N etc.) ---- */
#ifndef ISSE1_LR
#define ISSE1_LR 7
#endif
#ifndef ISSE2_LR
#define ISSE2_LR 5
#endif
#ifndef ISSE3_LR
#define ISSE3_LR 5
#endif
#ifndef ISSE4_LR
#define ISSE4_LR 4
#endif
#ifndef ISSE5_LR
#define ISSE5_LR 4
#endif
#ifndef WPRED_BASE_CONF
#define WPRED_BASE_CONF 8
#endif
#ifndef MATCH_CONF_DIV
#define MATCH_CONF_DIV 16
#endif

/* ---- Profiling (entropy map of predictions) ---- */
static int profile_mode = 0;

/* ---- Ablation (silence model families via bitmask) ---- */
/* bit 1=order_low, 2=order_high, 3=sparse, 4=linguistic,
   5=match, 6=icm, 7=delta, 8=stride, 9=wpred, 10=lzp */
static uint32_t ablate_mask = 0;

/* ---- Stretch / Squash Tables ---- */
static int16_t  stretch_tab[PROB_SCALE];
static uint16_t squash_tab[LOGIT_RANGE];

static void init_stretch_squash(void) {
    for (int i = 0; i < LOGIT_RANGE; i++) {
        double x = (double)(i + LOGIT_MIN) / LOGIT_SCALE;
        double p = 1.0 / (1.0 + exp(-x));
        int v = (int)(p * PROB_SCALE + 0.5);
        if (v < 1) v = 1;
        if (v > PROB_SCALE - 1) v = PROB_SCALE - 1;
        squash_tab[i] = (uint16_t)v;
    }
    for (int i = 0; i < PROB_SCALE; i++) {
        double p = ((double)i + 0.5) / PROB_SCALE;
        if (p < 0.001) p = 0.001;
        if (p > 0.999) p = 0.999;
        double logit = log(p / (1.0 - p)) * LOGIT_SCALE;
        int v = (int)(logit + (logit >= 0 ? 0.5 : -0.5));
        if (v < LOGIT_MIN) v = LOGIT_MIN;
        if (v > LOGIT_MAX) v = LOGIT_MAX;
        stretch_tab[i] = (int16_t)v;
    }
}

static inline int16_t stretch(uint16_t p) {
    return stretch_tab[p < PROB_SCALE ? p : PROB_SCALE - 1];
}

static inline uint16_t squash(int32_t logit) {
    if (logit < LOGIT_MIN) logit = LOGIT_MIN;
    if (logit > LOGIT_MAX) logit = LOGIT_MAX;
    return squash_tab[logit - LOGIT_MIN];
}

/* ---- PAQ-Style State Table ---- */
static uint8_t  next_state[256][2];
static uint16_t state_prob[256];
static int16_t  state_stretch[256];

static void init_state_table(void) {
    /* Improved state table: prioritize small counts for fast adaptation,
     * but allow higher counts for stable contexts. Uses (n0,n1) pairs
     * with a mapping that favors low counts and asymmetric distributions. */
    uint8_t sn0[256], sn1[256];
    int16_t state_map[42][42];
    memset(state_map, -1, sizeof(state_map));

    int idx = 0;
    /* First: allocate states by total count, favoring low counts */
    for (int sum = 0; sum <= 80 && idx < 256; sum++) {
        for (int n1 = 0; n1 <= sum && idx < 256; n1++) {
            int n0 = sum - n1;
            if (n0 > 41 || n1 > 41) continue;
            if (state_map[n0][n1] >= 0) continue;
            state_map[n0][n1] = idx;
            sn0[idx] = (uint8_t)n0;
            sn1[idx] = (uint8_t)n1;
            idx++;
        }
    }

    for (int s = 0; s < 256; s++) {
        for (int bit = 0; bit < 2; bit++) {
            int n0 = sn0[s], n1 = sn1[s];
            if (bit == 0) n0++; else n1++;
            while (n0 + n1 > 41) { n0 = (n0 + 1) >> 1; n1 = (n1 + 1) >> 1; }
            if (n0 > 41) n0 = 41;
            if (n1 > 41) n1 = 41;
            next_state[s][bit] = (state_map[n0][n1] >= 0) ?
                                 (uint8_t)state_map[n0][n1] : 0;
        }
    }

    for (int s = 0; s < 256; s++) {
        double p = (sn1[s] + 0.5) / (sn0[s] + sn1[s] + 1.0);
        int v = (int)(p * PROB_SCALE + 0.5);
        if (v < 1) v = 1;
        if (v > PROB_SCALE - 1) v = PROB_SCALE - 1;
        state_prob[s] = (uint16_t)v;
        state_stretch[s] = stretch_tab[v];
    }
}

/* ---- Context Model ---- */
typedef struct {
    union {
        uint8_t  *table8;
        uint16_t *table16;
    };
    uint32_t  mask;
    int       tagged;
} ContextModel;

static ContextModel models[NUM_CM];

static void cm_init(ContextModel *m, int bits, int tagged) {
    m->mask = (1u << bits) - 1;
    m->tagged = tagged;
    if (tagged) {
        m->table16 = (uint16_t *)calloc((size_t)(m->mask) + 1, sizeof(uint16_t));
        if (!m->table16) { fprintf(stderr, "OOM cm\n"); exit(1); }
    } else {
        m->table8 = (uint8_t *)calloc((size_t)(m->mask) + 1, 1);
        if (!m->table8) { fprintf(stderr, "OOM cm\n"); exit(1); }
    }
}

static void cm_free(ContextModel *m) {
    if (m->tagged)
        free(m->table16);
    else
        free(m->table8);
    m->table16 = NULL;
}

/* ---- Mixer ---- */
typedef struct {
    int32_t w[NUM_MODELS];
} MixerSlot;

#define MIX_W_SCALE  2048

static MixerSlot mixer[MIXER_SLOTS];
static MixerSlot mixer2[MIXER2_SLOTS];

static void mixer_init(void) {
    int init_w = MIX_W_SCALE / (NUM_MODELS - 4);
    for (int i = 0; i < MIXER_SLOTS; i++)
        for (int j = 0; j < NUM_MODELS; j++)
            mixer[i].w[j] = (j >= NUM_CM || j == 9 + NUM_SPARSE || j == 10 + NUM_SPARSE) ? 0 : init_w;
    for (int i = 0; i < MIXER2_SLOTS; i++)
        for (int j = 0; j < NUM_MODELS; j++)
            mixer2[i].w[j] = (j >= NUM_CM || j == 9 + NUM_SPARSE || j == 10 + NUM_SPARSE) ? 0 : init_w;
}

/* ---- APM (Adaptive Probability Map) ---- */
typedef struct {
    uint16_t *table;
    int ctx_count;
    int prob_bins;
} APM;

static APM isse1, isse2, isse3, isse4, isse5;

static void apm_init(APM *a, int ctx_bits, int prob_bits) {
    a->ctx_count = 1 << ctx_bits;
    a->prob_bins = 1 << prob_bits;
    a->table = (uint16_t *)malloc((size_t)a->ctx_count * a->prob_bins * sizeof(uint16_t));
    if (!a->table) { fprintf(stderr, "OOM apm\n"); exit(1); }
    for (int c = 0; c < a->ctx_count; c++)
        for (int p = 0; p < a->prob_bins; p++) {
            int v = (int)((p + 0.5) / a->prob_bins * PROB_SCALE);
            if (v < 1) v = 1;
            if (v > PROB_SCALE - 1) v = PROB_SCALE - 1;
            a->table[c * a->prob_bins + p] = (uint16_t)v;
        }
}

static void apm_free(APM *a) { free(a->table); a->table = NULL; }

static inline uint16_t apm_predict(APM *a, int ctx, uint16_t prob) {
    int scaled = (int)prob * (a->prob_bins - 1);
    int bin = scaled / PROB_SCALE;
    int frac = scaled % PROB_SCALE;
    if (bin >= a->prob_bins - 1) { bin = a->prob_bins - 2; frac = PROB_SCALE - 1; }
    int base = ctx * a->prob_bins;
    uint32_t p0 = a->table[base + bin];
    uint32_t p1 = a->table[base + bin + 1];
    uint32_t p = (p0 * (PROB_SCALE - frac) + p1 * frac) / PROB_SCALE;
    if (p < 1) p = 1;
    if (p > PROB_SCALE - 1) p = PROB_SCALE - 1;
    return (uint16_t)p;
}

static inline void apm_update_rate(APM *a, int ctx, uint16_t prob, int bit, int rate) {
    int scaled = (int)prob * (a->prob_bins - 1);
    int bin = scaled / PROB_SCALE;
    if (bin >= a->prob_bins - 1) bin = a->prob_bins - 2;
    int base = ctx * a->prob_bins;
    int target = bit ? PROB_SCALE : 0;
    for (int d = 0; d < 2; d++) {
        int idx = base + bin + d;
        int err = target - (int)a->table[idx];
        a->table[idx] = (uint16_t)((int)a->table[idx] + (err >> rate));
        if (a->table[idx] < 1) a->table[idx] = 1;
        if (a->table[idx] > PROB_SCALE - 1) a->table[idx] = PROB_SCALE - 1;
    }
}

/* ---- Match Model ---- */
#define MATCH_CHAIN_DEPTH 16

typedef struct {
    uint32_t *pos_table;   /* hash -> head of chain (position, 1-based) */
    uint32_t *prev;        /* prev[pos] = previous pos in chain with same hash */
    uint32_t match_pos;
    int match_len;         /* backward context match length (bytes), then +1 per matching bit */
    int match_active;
    uint32_t rep[3];
    int      rep_count;
} MatchModel;

static MatchModel mmatch;

static void match_init(void) {
    mmatch.pos_table = (uint32_t *)calloc(MATCH_TAB_SIZE, sizeof(uint32_t));
    mmatch.prev = (uint32_t *)calloc(BLOCK_SIZE + 1, sizeof(uint32_t));
    if (!mmatch.pos_table || !mmatch.prev) { fprintf(stderr, "OOM match\n"); exit(1); }
    mmatch.match_pos = 0;
    mmatch.match_len = 0;
    mmatch.match_active = 0;
    mmatch.rep[0] = mmatch.rep[1] = mmatch.rep[2] = 0;
    mmatch.rep_count = 0;
}

static void match_free(void) {
    free(mmatch.pos_table); mmatch.pos_table = NULL;
    free(mmatch.prev); mmatch.prev = NULL;
}

/* ---- Indirect Context Model (ICM) ---- */
static uint8_t  *icm_table;              /* context hash -> counter index */
static uint16_t  icm_ctr[ICM_NCTR];      /* adaptive probabilities */

static void icm_init(void) {
    icm_table = (uint8_t *)calloc(ICM_TAB_SIZE, 1);
    if (!icm_table) { fprintf(stderr, "OOM icm\n"); exit(1); }
    for (int i = 0; i < ICM_NCTR; i++) icm_ctr[i] = PROB_HALF;
}

static void icm_free(void) { free(icm_table); icm_table = NULL; }

/* ---- Delta Model ---- */
static uint16_t *delta_table;

static void delta_init(void) {
    delta_table = (uint16_t *)calloc(DELTA_TAB_SIZE, sizeof(uint16_t));
    if (!delta_table) { fprintf(stderr, "OOM delta\n"); exit(1); }
}
static void delta_free(void) { free(delta_table); delta_table = NULL; }

/* ---- Stride/Record Model ---- */
static uint16_t *stride_table;

static void stride_init(void) {
    stride_table = (uint16_t *)calloc(STRIDE_TAB_SIZE, sizeof(uint16_t));
    if (!stride_table) { fprintf(stderr, "OOM stride\n"); exit(1); }
}
static void stride_free(void) { free(stride_table); stride_table = NULL; }

/* ---- LZP (Lempel-Ziv Prediction) Model ---- */
static uint16_t *lzp_table;  /* tag(8) | predicted_byte(8) */

static void lzp_init(void) {
    lzp_table = (uint16_t *)calloc(LZP_SIZE, sizeof(uint16_t));
    if (!lzp_table) { fprintf(stderr, "OOM lzp\n"); exit(1); }
}
static void lzp_free(void) { free(lzp_table); lzp_table = NULL; }

static int detect_stride(const uint8_t *blk, uint32_t bl) {
    int best_p = 0, best_score = 0;
    uint32_t limit = bl < 65536 ? bl : 65536;
    for (int p = 2; p <= MAX_STRIDE && (uint32_t)p < bl; p++) {
        int score = 0;
        for (uint32_t i = p; i < limit; i += 7)
            if (blk[i] == blk[i - p]) score++;
        if (score > best_score) { best_score = score; best_p = p; }
    }
    int random_expect = (int)(limit / 7) / 256;
    return (best_score > random_expect * 4) ? best_p : 0;
}

/* ---- Word Predictor (Shadow Dictionary) ---- */
#define WPRED_BITS      18    /* 256K entries */
#define WPRED_SIZE      (1u << WPRED_BITS)
#define WPRED_MAX_WORD  24

typedef struct {
    uint16_t tag;
    uint8_t  len;
    uint8_t  word[WPRED_MAX_WORD]; /* lowercase predicted bytes */
} WPredEntry;

static WPredEntry *wpred_table;   /* unigram: prev_word → next_word */
static WPredEntry *wpred_table2;  /* bigram:  (prev_prev + prev) → next_word */

static void wpred_init(void) {
    wpred_table = (WPredEntry *)calloc(WPRED_SIZE, sizeof(WPredEntry));
    wpred_table2 = (WPredEntry *)calloc(WPRED_SIZE, sizeof(WPredEntry));
    if (!wpred_table || !wpred_table2) { fprintf(stderr, "OOM wpred\n"); exit(1); }
}

static void wpred_free(void) {
    free(wpred_table); wpred_table = NULL;
    free(wpred_table2); wpred_table2 = NULL;
}

/* ---- BCJ x86 Filter ---- */
static void bcj_encode(uint8_t *buf, uint32_t len) {
    if (len < 5) return;
    for (uint32_t i = 0; i + 4 < len; i++) {
        if (buf[i] == 0xE8) {
            uint32_t rel = (uint32_t)buf[i+1] | ((uint32_t)buf[i+2] << 8) |
                           ((uint32_t)buf[i+3] << 16) | ((uint32_t)buf[i+4] << 24);
            uint32_t abs_addr = rel + i + 5;
            buf[i+1] = (uint8_t)(abs_addr);
            buf[i+2] = (uint8_t)(abs_addr >> 8);
            buf[i+3] = (uint8_t)(abs_addr >> 16);
            buf[i+4] = (uint8_t)(abs_addr >> 24);
            i += 4;
        }
    }
}

static void bcj_decode(uint8_t *buf, uint32_t len) {
    if (len < 5) return;
    for (uint32_t i = 0; i + 4 < len; i++) {
        if (buf[i] == 0xE8) {
            uint32_t abs_addr = (uint32_t)buf[i+1] | ((uint32_t)buf[i+2] << 8) |
                                ((uint32_t)buf[i+3] << 16) | ((uint32_t)buf[i+4] << 24);
            uint32_t rel = abs_addr - i - 5;
            buf[i+1] = (uint8_t)(rel);
            buf[i+2] = (uint8_t)(rel >> 8);
            buf[i+3] = (uint8_t)(rel >> 16);
            buf[i+4] = (uint8_t)(rel >> 24);
            i += 4;
        }
    }
}

/* ---- BMP Vertical Delta Filter ---- */
static void bmp_delta_encode(uint8_t *buf, uint32_t len, int data_offset, int row_bytes) {
    if (row_bytes <= 0 || data_offset < 0) return;
    int num_rows = (len - data_offset) / row_bytes;
    if (num_rows < 2) return;
    /* Subtract row above, working backward to preserve source rows */
    for (int y = num_rows - 1; y >= 1; y--) {
        uint8_t *row = buf + data_offset + y * row_bytes;
        uint8_t *prev = row - row_bytes;
        for (int x = 0; x < row_bytes; x++)
            row[x] -= prev[x];
    }
}

static void bmp_delta_decode(uint8_t *buf, uint32_t len, int data_offset, int row_bytes) {
    if (row_bytes <= 0 || data_offset < 0) return;
    int num_rows = (len - data_offset) / row_bytes;
    if (num_rows < 2) return;
    /* Add row above, working forward */
    for (int y = 1; y < num_rows; y++) {
        uint8_t *row = buf + data_offset + y * row_bytes;
        uint8_t *prev = row - row_bytes;
        for (int x = 0; x < row_bytes; x++)
            row[x] += prev[x];
    }
}

/* Raw image delta filter — uses 16-bit arithmetic when stride is even,
 * byte-level otherwise. 16-bit subtraction handles carry across byte
 * boundaries correctly for 16-bit grayscale/depth images. */
static void raw_delta_encode(uint8_t *buf, uint32_t len, int row_bytes) {
    int num_rows = (int)(len / row_bytes);
    if (num_rows < 2) return;
    int use16 = (row_bytes % 2 == 0);
    for (int y = num_rows - 1; y >= 1; y--) {
        uint8_t *row = buf + y * row_bytes;
        uint8_t *prev = row - row_bytes;
        if (use16) {
            for (int x = 0; x < row_bytes; x += 2) {
                uint16_t cur = (uint16_t)(row[x] | (row[x + 1] << 8));
                uint16_t prv = (uint16_t)(prev[x] | (prev[x + 1] << 8));
                uint16_t d = cur - prv;
                row[x] = (uint8_t)d;
                row[x + 1] = (uint8_t)(d >> 8);
            }
        } else {
            for (int x = 0; x < row_bytes; x++)
                row[x] -= prev[x];
        }
    }
}

static void raw_delta_decode(uint8_t *buf, uint32_t len, int row_bytes) {
    int num_rows = (int)(len / row_bytes);
    if (num_rows < 2) return;
    int use16 = (row_bytes % 2 == 0);
    for (int y = 1; y < num_rows; y++) {
        uint8_t *row = buf + y * row_bytes;
        uint8_t *prev = row - row_bytes;
        if (use16) {
            for (int x = 0; x < row_bytes; x += 2) {
                uint16_t cur = (uint16_t)(row[x] | (row[x + 1] << 8));
                uint16_t prv = (uint16_t)(prev[x] | (prev[x + 1] << 8));
                uint16_t d = cur + prv;
                row[x] = (uint8_t)d;
                row[x + 1] = (uint8_t)(d >> 8);
            }
        } else {
            for (int x = 0; x < row_bytes; x++)
                row[x] += prev[x];
        }
    }
}

/* Detect raw image geometry by finding the row width that minimizes
 * vertical byte differences. Returns row_bytes (>= 128) or 0 if not an image. */
static int detect_image_geometry(const uint8_t *buf, uint32_t len) {
    if (len < 65536) return 0;
    uint32_t sample = len < 262144 ? len : 262144;

    /* Baseline: average |buf[i] - buf[i-1]| */
    int64_t base_diff = 0;
    int base_n = 0;
    for (uint32_t i = 1; i < sample; i += 7) {
        base_diff += abs((int)buf[i] - (int)buf[i - 1]);
        base_n++;
    }
    if (base_n == 0) return 0;
    int64_t base_avg = base_diff / base_n;
    if (base_avg < 8) return 0;  /* data already very smooth */

    int best_w = 0;
    int64_t best_avg = base_avg;

    /* Test candidate widths: multiples of 4 from 128 to 16384 */
    for (int w = 128; w <= 16384 && (uint32_t)w < sample / 4; w += 4) {
        int64_t diff = 0;
        int n = 0;
        for (uint32_t i = (uint32_t)w; i < sample; i += 31) {
            diff += abs((int)buf[i] - (int)buf[i - w]);
            n++;
        }
        if (n == 0) continue;
        int64_t avg = diff / n;
        if (avg < best_avg) {
            best_avg = avg;
            best_w = w;
        }
    }

    /* Accept only if vertical diff is dramatically lower (< 50% of horizontal) */
    if (best_w > 0 && best_avg * 2 < base_avg)
        return best_w;
    return 0;
}

/* Scan for "data" sub-chunk in RIFF/WAVE, return offset of audio data or -1 */
static int wav_find_data(const uint8_t *buf, uint32_t len) {
    if (len < 44) return -1;
    uint32_t off = 12; /* skip "RIFF" + size + "WAVE" */
    while (off + 8 <= len) {
        uint32_t chunk_size = buf[off+4] | (buf[off+5]<<8) | (buf[off+6]<<16) | (buf[off+7]<<24);
        if (buf[off]=='d' && buf[off+1]=='a' && buf[off+2]=='t' && buf[off+3]=='a')
            return (int)(off + 8);
        off += 8 + ((chunk_size + 1) & ~1); /* chunks are word-aligned */
    }
    return -1;
}

/* ---- WAV Per-Channel Delta Filter ---- */
static void wav_delta_encode(uint8_t *buf, uint32_t len, int data_offset, int num_channels) {
    if (num_channels < 1 || data_offset < 0 || (uint32_t)data_offset >= len) return;
    int total_samples = (len - data_offset) / 2;
    if (total_samples <= num_channels) return;
    int16_t *samples = (int16_t *)(buf + data_offset);
    for (int i = total_samples - 1; i >= num_channels; i--)
        samples[i] -= samples[i - num_channels];
}

static void wav_delta_decode(uint8_t *buf, uint32_t len, int data_offset, int num_channels) {
    if (num_channels < 1 || data_offset < 0 || (uint32_t)data_offset >= len) return;
    int total_samples = (len - data_offset) / 2;
    if (total_samples <= num_channels) return;
    int16_t *samples = (int16_t *)(buf + data_offset);
    for (int i = num_channels; i < total_samples; i++)
        samples[i] += samples[i - num_channels];
}

/* ---- MTF Word Cache Pre-filter ---- */
#define MTF_CACHE_SIZE   128
#define MTF_MIN_WORD     3
#define MTF_MAX_WORD     63

typedef struct {
    uint8_t words[MTF_CACHE_SIZE][MTF_MAX_WORD + 1];
    int     lens[MTF_CACHE_SIZE];
    int     count;
} MTFCache;

static void mtf_init(MTFCache *c) { c->count = 0; }

static int mtf_find(MTFCache *c, const uint8_t *word, int len) {
    for (int i = 0; i < c->count; i++)
        if (c->lens[i] == len && memcmp(c->words[i], word, len) == 0)
            return i;
    return -1;
}

static void mtf_promote(MTFCache *c, int idx) {
    if (idx == 0) return;
    uint8_t tmp[MTF_MAX_WORD + 1];
    int tmp_len = c->lens[idx];
    memcpy(tmp, c->words[idx], tmp_len);
    memmove(&c->words[1], &c->words[0], idx * sizeof(c->words[0]));
    memmove(&c->lens[1], &c->lens[0], idx * sizeof(c->lens[0]));
    memcpy(c->words[0], tmp, tmp_len);
    c->lens[0] = tmp_len;
}

static void mtf_insert(MTFCache *c, const uint8_t *word, int len) {
    if (c->count < MTF_CACHE_SIZE) c->count++;
    memmove(&c->words[1], &c->words[0], (c->count - 1) * sizeof(c->words[0]));
    memmove(&c->lens[1], &c->lens[0], (c->count - 1) * sizeof(c->lens[0]));
    memcpy(c->words[0], word, len);
    c->lens[0] = len;
}

static inline int is_word_char(uint8_t c) {
    return (c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z') ||
           (c >= '0' && c <= '9') || c == '_';
}

/* 0=mixed, 1=all lower, 2=title (first upper rest lower), 3=all upper */
static int detect_word_case(const uint8_t *w, int len) {
    int has_upper = 0, has_lower = 0, first_upper = 0;
    int upper_after_first = 0;
    for (int i = 0; i < len; i++) {
        if (w[i] >= 'A' && w[i] <= 'Z') {
            has_upper = 1;
            if (i == 0) first_upper = 1;
            else upper_after_first = 1;
        }
        if (w[i] >= 'a' && w[i] <= 'z') has_lower = 1;
    }
    if (!has_upper) return 1;                                   /* all lower */
    if (!has_lower) return 3;                                   /* all upper */
    if (first_upper && !upper_after_first) return 2;            /* title */
    return 0;                                                   /* mixed */
}

static void to_lower(const uint8_t *src, uint8_t *dst, int len) {
    for (int i = 0; i < len; i++)
        dst[i] = (src[i] >= 'A' && src[i] <= 'Z') ? (src[i] + 32) : src[i];
}

static uint32_t text_mtf_encode(const uint8_t *in, uint32_t in_len,
                                uint8_t *out, uint32_t out_cap) {
    MTFCache cache;
    mtf_init(&cache);
    uint32_t ip = 0, op = 0;
    uint8_t lower_buf[MTF_MAX_WORD + 1];

    while (ip < in_len) {
        /* Check for literal escape of 0x01-0x03 */
        if (in[ip] >= 0x01 && in[ip] <= 0x03) {
            if (op + 2 > out_cap) break;
            out[op++] = in[ip];
            out[op++] = 0x80;
            ip++;
            continue;
        }

        /* Not a word char — copy literally */
        if (!is_word_char(in[ip])) {
            if (op + 1 > out_cap) break;
            out[op++] = in[ip++];
            continue;
        }

        /* Accumulate word */
        uint32_t wstart = ip;
        while (ip < in_len && is_word_char(in[ip]) && (int)(ip - wstart) < MTF_MAX_WORD)
            ip++;
        int wlen = (int)(ip - wstart);

        if (wlen < MTF_MIN_WORD) {
            /* Too short — emit literally */
            if (op + (uint32_t)wlen > out_cap) break;
            memcpy(out + op, in + wstart, wlen);
            op += wlen;
            continue;
        }

        int wcase = detect_word_case(in + wstart, wlen);
        if (wcase == 0) {
            /* Mixed case — emit literally, do NOT update cache */
            if (op + (uint32_t)wlen > out_cap) break;
            memcpy(out + op, in + wstart, wlen);
            op += wlen;
            continue;
        }

        /* Lowercase for cache lookup */
        to_lower(in + wstart, lower_buf, wlen);
        int idx = mtf_find(&cache, lower_buf, wlen);

        if (idx >= 0) {
            /* Cache hit — emit reference */
            if (op + 2 > out_cap) break;
            out[op++] = (uint8_t)wcase;  /* 1=lower, 2=title, 3=upper */
            out[op++] = (uint8_t)idx;
            mtf_promote(&cache, idx);
        } else {
            /* Cache miss — emit literally, insert */
            if (op + (uint32_t)wlen > out_cap) break;
            memcpy(out + op, in + wstart, wlen);
            op += wlen;
            mtf_insert(&cache, lower_buf, wlen);
        }
    }
    return op;
}

static uint32_t text_mtf_decode(const uint8_t *in, uint32_t in_len,
                                uint8_t *out, uint32_t out_cap) {
    MTFCache cache;
    mtf_init(&cache);
    uint32_t ip = 0, op = 0;
    uint8_t lower_buf[MTF_MAX_WORD + 1];

    while (ip < in_len) {
        /* Check for escape bytes */
        if (in[ip] >= 0x01 && in[ip] <= 0x03) {
            if (ip + 1 >= in_len) break;
            uint8_t esc = in[ip];
            uint8_t next = in[ip + 1];
            if (next == 0x80) {
                /* Literal escape */
                if (op + 1 > out_cap) break;
                out[op++] = esc;
                ip += 2;
                continue;
            }
            /* Word reference */
            int idx = (int)next;
            if (idx >= cache.count) break;  /* invalid ref */
            int wlen = cache.lens[idx];
            if (op + (uint32_t)wlen > out_cap) break;

            /* Copy word with case applied */
            const uint8_t *w = cache.words[idx];
            if (esc == 1) {
                /* lowercase */
                memcpy(out + op, w, wlen);
            } else if (esc == 2) {
                /* title case */
                out[op] = (w[0] >= 'a' && w[0] <= 'z') ? (w[0] - 32) : w[0];
                if (wlen > 1) memcpy(out + op + 1, w + 1, wlen - 1);
            } else {
                /* all upper */
                for (int j = 0; j < wlen; j++)
                    out[op + j] = (w[j] >= 'a' && w[j] <= 'z') ? (w[j] - 32) : w[j];
            }
            op += wlen;
            mtf_promote(&cache, idx);
            ip += 2;
            continue;
        }

        /* Not a word char — copy literally */
        if (!is_word_char(in[ip])) {
            if (op + 1 > out_cap) break;
            out[op++] = in[ip++];
            continue;
        }

        /* Accumulate literal word */
        uint32_t wstart = ip;
        uint32_t ostart = op;
        while (ip < in_len && is_word_char(in[ip]) && (int)(ip - wstart) < MTF_MAX_WORD) {
            if (op + 1 > out_cap) break;
            out[op++] = in[ip++];
        }
        int wlen = (int)(ip - wstart);

        if (wlen >= MTF_MIN_WORD) {
            /* Insert into cache only if word has valid case (not mixed) —
             * must mirror the encoder, which skips mixed-case words */
            int wcase = detect_word_case(out + ostart, wlen);
            if (wcase != 0) {
                to_lower(out + ostart, lower_buf, wlen);
                if (mtf_find(&cache, lower_buf, wlen) < 0)
                    mtf_insert(&cache, lower_buf, wlen);
            }
        }
    }
    return op;
}

/* ---- Binary rANS ---- */
static inline uint32_t rans_enc_bit(uint32_t r, uint8_t **p, int bit, uint16_t prob) {
    uint16_t f = bit ? prob : (PROB_SCALE - prob);
    uint32_t xm = ((RANS_L >> PROB_BITS) << 8) * f;
    while (r >= xm) { *(--(*p)) = (uint8_t)(r & 0xFF); r >>= 8; }
    if (bit)
        return ((r / f) << PROB_BITS) + (r % f) + (PROB_SCALE - prob);
    else
        return ((r / f) << PROB_BITS) + (r % f);
}

static inline void rans_flush(uint32_t r, uint8_t **p) {
    *(--(*p)) = (uint8_t)(r >> 24);
    *(--(*p)) = (uint8_t)(r >> 16);
    *(--(*p)) = (uint8_t)(r >> 8);
    *(--(*p)) = (uint8_t)(r);
}

static inline uint32_t rans_dec_init(const uint8_t **p) {
    uint32_t r = (uint32_t)(*(*p)++);
    r |= (uint32_t)(*(*p)++) << 8;
    r |= (uint32_t)(*(*p)++) << 16;
    r |= (uint32_t)(*(*p)++) << 24;
    return r;
}

static inline int rans_dec_bit(uint32_t *r, const uint8_t **p, uint16_t prob) {
    uint32_t x = *r;
    uint32_t cumfreq = x & (PROB_SCALE - 1);
    int bit;
    uint16_t f;
    if (cumfreq >= (uint32_t)(PROB_SCALE - prob)) {
        bit = 1; f = prob;
        x = f * (x >> PROB_BITS) + cumfreq - (PROB_SCALE - prob);
    } else {
        bit = 0; f = PROB_SCALE - prob;
        x = f * (x >> PROB_BITS) + cumfreq;
    }
    while (x < RANS_L) x = (x << 8) | (uint32_t)(*(*p)++);
    *r = x;
    return bit;
}

/* ---- Event Buffer ---- */
static uint16_t *ev_buf;
static int32_t   ev_cnt;

static inline void ev_emit(int bit, uint16_t prob) {
    ev_buf[ev_cnt++] = (uint16_t)((bit << 15) | prob);
}

/* ---- Context Hashing ---- */
static inline uint32_t hmix(uint32_t h, uint32_t b) {
    return (h + b + 1) * 2654435761u;
}

/* Finalize hash to reduce clustering */
static inline uint32_t hfin(uint32_t h) {
    h ^= h >> 16;
    h *= 0x45d9f3bu;
    return h;
}

/* ---- File I/O ---- */
static void wr32(FILE *f, uint32_t v) {
    uint8_t b[4] = {v, v>>8, v>>16, v>>24};
    fwrite(b, 1, 4, f);
}

static uint32_t rd32(FILE *f) {
    uint8_t b[4]; fread(b, 1, 4, f);
    return (uint32_t)b[0]|((uint32_t)b[1]<<8)|((uint32_t)b[2]<<16)|((uint32_t)b[3]<<24);
}

/* ---- Shared Modeling Logic ---- */
static void process_block(uint8_t *blk, uint32_t bl, int mode,
                          uint32_t *rs, const uint8_t **dec_ptr, int32_t *ev_idx_p,
                          int stride, int force_text) {
    /* Don't reset context model tables — let them learn across blocks */

    match_init();

    /* Text mode is determined at file level and stored in the header (FLAG_TEXT),
     * ensuring encoder and decoder agree. Block-level detection was unreliable
     * because the decoder's buffer starts as zeros before decoding. */
    int is_text_block = force_text;

    /* Text blocks have no meaningful stride — the detector gives false positives
     * due to non-uniform byte distribution (spaces, common letters). */
    if (is_text_block) stride = 0;

    /* Deep match chain for text: literature has long-range repetitions */
    int chain_depth = is_text_block ? 512 : MATCH_CHAIN_DEPTH;

    /* Previous byte history for context building */
    uint8_t pb[16] = {0}; /* pb[0]=prev, pb[1]=prev2, ... pb[15]=prev16 */
    uint32_t word_hash = 0;      /* rolling hash of current word */
    uint32_t prev_word_hash = 0; /* hash of last completed word */
    uint32_t prev_prev_word_hash = 0; /* hash of word before last completed word */
    int word_pos = 0;            /* character index within current word */
    int32_t is_text = 0;     /* text/binary detector: positive = text, negative = binary */

    /* Word predictor state (shadow dictionary) */
    uint8_t wpred_word[WPRED_MAX_WORD]; /* predicted next word (lowercase) */
    int wpred_len = 0;    /* length of predicted word */
    int wpred_pos = 0;    /* current char index within prediction */
    int wpred_active = 0; /* whether prediction is active */
    int wpred_was_active = 0; /* was prediction active at start of this word? */
    int wpred_from_bigram = 0; /* was current prediction from bigram table? */
    uint8_t cur_word_lc[WPRED_MAX_WORD]; /* current word accumulator (lowercase) */
    int cur_word_len = 0;

    /* wpred diagnostic counters (only used with --profile) */
    uint32_t wp_lookups = 0, wp_bigram_hit = 0, wp_unigram_hit = 0, wp_miss = 0;
    uint32_t wp_bigram_collision = 0, wp_unigram_collision = 0;
    uint32_t wp_word_correct = 0, wp_word_partial = 0, wp_word_wrong = 0;
    uint32_t wp_bi_correct = 0, wp_bi_partial = 0, wp_bi_wrong = 0;
    uint32_t wp_uni_correct = 0, wp_uni_partial = 0, wp_uni_wrong = 0;

    if (mode == 0) ev_cnt = 0;

    for (uint32_t i = 0; i < bl; i++) {
        /* Build context hashes from previous bytes */
        uint32_t ctx[NUM_CM];
        ctx[0] = 0;  /* order-0: just bit_ctx, no byte context */
        ctx[1] = hfin(hmix(0, pb[0]));                                        /* order-1 */
        ctx[2] = hfin(hmix(hmix(0, pb[1]), pb[0]));                            /* order-2 */
        ctx[3] = hfin(hmix(hmix(hmix(0, pb[2]), pb[1]), pb[0]));               /* order-3 */
        ctx[4] = hfin(hmix(hmix(hmix(hmix(0, pb[3]), pb[2]), pb[1]), pb[0]));  /* order-4 */
        ctx[5] = i >= 5 ? hfin(hmix(hmix(hmix(hmix(hmix(hmix(0,
                    pb[5]), pb[4]), pb[3]), pb[2]), pb[1]), pb[0])) : ctx[4];   /* order-6 */
        ctx[6] = i >= 7 ? hfin(hmix(hmix(hmix(hmix(hmix(hmix(hmix(hmix(0,
                    pb[7]), pb[6]), pb[5]), pb[4]), pb[3]), pb[2]), pb[1]), pb[0])) : ctx[5]; /* order-8 */
        ctx[7] = i >= 11 ? hfin(hmix(hmix(hmix(hmix(ctx[6],
                    pb[8]), pb[9]), pb[10]), pb[11])) : ctx[6]; /* order-12 */
        ctx[8] = hfin(word_hash);  /* word context */

        /*
         * Sparse context models: non-adjacent byte combinations.
         * Each captures structural correlations that sequential order-N models
         * miss (e.g. fields at fixed offsets in ELF headers, instruction
         * patterns with register encodings at regular intervals).
         *
         * To add more sparse models when memory is not a concern:
         *   - Increment NUM_SPARSE in the defines
         *   - Add new ctx[9 + N] entries here with novel byte-skip patterns
         *   - Add matching cm_init() calls in do_compress() and do_decompress()
         *   - Expand pb[] if patterns need bytes beyond pb[15]
         * Each model costs ~8 MB and yields ~200-500 bytes at this point.
         */
        ctx[9]  = i >= 2 ? hfin(hmix(hmix(0, pb[2]), pb[0])) : 0;              /* sparse1:  skip-1 order-2 */
        ctx[10] = i >= 3 ? hfin(hmix(hmix(hmix(0, pb[3]), pb[1]), pb[0])) : 0;  /* sparse2:  (3,1,0) */
        ctx[11] = i >= 4 ? hfin(hmix(hmix(0, pb[4]), pb[0])) : 0;              /* sparse3:  skip-3 order-2 */
        ctx[12] = i >= 5 ? hfin(hmix(hmix(hmix(0, pb[5]), pb[2]), pb[0])) : 0;  /* sparse4:  stride-2+ order-3 */
        ctx[13] = i >= 7 ? hfin(hmix(hmix(hmix(0, pb[7]), pb[3]), pb[0])) : 0;  /* sparse5:  stride-3+ order-3 */
        ctx[14] = i >= 3 ? hfin(hmix(hmix(hmix(0, pb[3]), pb[2]), pb[0])) : 0;  /* sparse6:  (3,2,0) */
        ctx[15] = i >= 5 ? hfin(hmix(hmix(0, pb[5]), pb[1])) : 0;              /* sparse7:  (5,1) */
        ctx[16] = i >= 8 ? hfin(hmix(hmix(hmix(0, pb[8]), pb[4]), pb[0])) : 0;  /* sparse8:  stride-4 order-3 */
        ctx[17] = i >= 6 ? hfin(hmix(hmix(hmix(0, pb[6]), pb[3]), pb[0])) : 0;  /* sparse9:  stride-3 order-3 */
        ctx[18] = i >= 6 ? hfin(hmix(hmix(0, pb[6]), pb[2])) : 0;              /* sparse10: (6,2) */
        ctx[19] = i >= 4 ? hfin(hmix(hmix(hmix(hmix(0, pb[4]), pb[3]), pb[1]), pb[0])) : 0; /* sparse11: skip-1 order-4 */
        ctx[20] = i >= 10 ? hfin(hmix(hmix(hmix(0, pb[10]), pb[5]), pb[0])) : 0;/* sparse12: stride-5 order-3 */
        ctx[21] = i >= 3 ? hfin(hmix(hmix(0, pb[1]), pb[3])) : 0;              /* sparse13: (1,3) reverse */
        ctx[22] = i >= 8 ? hfin(hmix(hmix(0, pb[8]), pb[0])) : 0;              /* sparse14: skip-7 order-2 */
        ctx[23] = i >= 12 ? hfin(hmix(hmix(hmix(0, pb[12]), pb[6]), pb[0])) : 0;/* sparse15: stride-6 order-3 */
        ctx[24] = i >= 6 ? hfin(hmix(hmix(hmix(hmix(0, pb[6]), pb[4]), pb[2]), pb[0])) : 0; /* sparse16: stride-2 order-4 */
        ctx[25] = i >= 9 ? hfin(hmix(hmix(hmix(0, pb[9]), pb[3]), pb[0])) : 0;  /* sparse17: (9,3,0) */
        ctx[26] = i >= 5 ? hfin(hmix(hmix(hmix(0, pb[5]), pb[3]), pb[1])) : 0;  /* sparse18: odd positions */
        ctx[27] = i >= 7 ? hfin(hmix(hmix(0, pb[7]), pb[1])) : 0;              /* sparse19: (7,1) */
        ctx[28] = i >= 10 ? hfin(hmix(hmix(0, pb[10]), pb[2])) : 0;            /* sparse20: (10,2) */
        ctx[29] = i >= 4 ? hfin(hmix(hmix(hmix(0, pb[4]), pb[2]), pb[1])) : 0;  /* sparse21: (4,2,1) */
        ctx[30] = i >= 8 ? hfin(hmix(hmix(hmix(hmix(0, pb[8]), pb[6]), pb[4]), pb[2])) : 0; /* sparse22: even stride-2 */
        ctx[31] = i >= 12 ? hfin(hmix(hmix(hmix(hmix(0, pb[12]), pb[8]), pb[4]), pb[0])) : 0; /* sparse23: stride-4 order-4 */
        ctx[32] = i >= 9 ? hfin(hmix(hmix(0, pb[9]), pb[0])) : 0;              /* sparse24: (9,0) */
        ctx[33] = i >= 7 ? hfin(hmix(hmix(hmix(0, pb[7]), pb[5]), pb[3])) : 0;  /* sparse25: odd stride-2 */
        ctx[34] = i >= 6 ? hfin(hmix(hmix(hmix(hmix(0, pb[6]), pb[5]), pb[1]), pb[0])) : 0; /* sparse26: mixed */
        ctx[35] = i >= 11 ? hfin(hmix(hmix(0, pb[11]), pb[0])) : 0;            /* sparse27: (11,0) */
        ctx[36] = i >= 14 ? hfin(hmix(hmix(hmix(0, pb[14]), pb[7]), pb[0])) : 0;/* sparse28: stride-7 order-3 */
        ctx[37] = i >= 10 ? hfin(hmix(hmix(hmix(hmix(hmix(0, pb[10]), pb[8]), pb[6]), pb[4]), pb[2])) : 0; /* sparse29: even stride-2 deep */
        ctx[38] = i >= 9 ? hfin(hmix(hmix(hmix(hmix(0, pb[9]), pb[6]), pb[3]), pb[0])) : 0; /* sparse30: stride-3 order-4 */
        ctx[39] = i >= 13 ? hfin(hmix(hmix(0, pb[13]), pb[0])) : 0;            /* sparse31: (13,0) */
        ctx[40] = i >= 15 ? hfin(hmix(hmix(hmix(0, pb[15]), pb[10]), pb[5])) : 0;/* sparse32: stride-5 deep */
        ctx[41] = i >= 11 ? hfin(hmix(hmix(hmix(hmix(0, pb[11]), pb[8]), pb[5]), pb[2])) : 0; /* sparse33: stride-3 shifted */
        ctx[42] = i >= 4 ? hfin(hmix(hmix(hmix(0, pb[4]), pb[0]), pb[1])) : 0;  /* sparse34: nonlinear mix */

        /* Previous-word context: predicts current word from previous word + position */
        ctx[9 + NUM_SPARSE] = (prev_word_hash != 0)
            ? hfin(hmix(prev_word_hash, (uint32_t)(word_pos < 16 ? word_pos : 15)))
            : 0;

        /* Word trigram context: prev_prev_word + prev_word + position */
        ctx[10 + NUM_SPARSE] = (prev_prev_word_hash != 0 && prev_word_hash != 0)
            ? hfin(hmix(hmix(prev_prev_word_hash, prev_word_hash),
                        (uint32_t)(word_pos < 16 ? word_pos : 15)))
            : 0;

        /* order-16 (always last) */
        ctx[NUM_CM - 1] = i >= 15 ? hfin(hmix(hmix(hmix(hmix(ctx[7],
                    pb[12]), pb[13]), pb[14]), pb[15])) : ctx[7];

        /* Match model: deep hash chain search + rep-matches */
        mmatch.match_active = 0;
        if (i >= 4) {
            int best_ext = 0;
            uint32_t best_pos = 0;

            /* Try rep-matches with backward extension */
            for (int r = 0; r < mmatch.rep_count; r++) {
                uint32_t d = mmatch.rep[r];
                if (d == 0 || i < d + 4) continue;
                uint32_t cand = i - d;
                if (blk[cand-4]==pb[3] && blk[cand-3]==pb[2] &&
                    blk[cand-2]==pb[1] && blk[cand-1]==pb[0]) {
                    int ext = 4;
                    for (int k = 5; k <= 16 && (uint32_t)k <= cand && (uint32_t)k <= i; k++) {
                        if (blk[cand - k] != pb[k - 1]) break;
                        ext++;
                    }
                    if (ext > best_ext) { best_ext = ext; best_pos = cand; }
                }
            }

            /* Deep hash chain walk */
            uint32_t mctx = hfin(hmix(hmix(hmix(hmix(0,
                        pb[3]), pb[2]), pb[1]), pb[0]));
            uint32_t mhash = mctx & (MATCH_TAB_SIZE - 1);
            uint32_t cand = mmatch.pos_table[mhash];
            /* Skip self-referencing head (insert at i-1 sets pos_table[hash]=i) */
            if (cand == i) cand = mmatch.prev[cand];
            for (int depth = 0; depth < chain_depth && cand > 0 && cand < i && cand < bl; depth++) {
                if (cand >= 4) {
                    int ext = 0;
                    for (int k = 1; k <= 16 && (uint32_t)k <= cand && (uint32_t)k <= i; k++) {
                        if (blk[cand - k] != pb[k - 1]) break;
                        ext++;
                    }
                    if (ext > best_ext) { best_ext = ext; best_pos = cand; }
                    if (ext >= 15) break;
                }
                cand = mmatch.prev[cand];
            }
            if (best_pos > 0 && best_ext >= 4) {
                mmatch.match_pos = best_pos;
                mmatch.match_len = best_ext;
                mmatch.match_active = 1;
                /* Update rep queue */
                uint32_t dist = i - best_pos;
                int found = 0;
                for (int r = 0; r < mmatch.rep_count; r++) {
                    if (mmatch.rep[r] == dist) {
                        for (int j = r; j > 0; j--) mmatch.rep[j] = mmatch.rep[j-1];
                        mmatch.rep[0] = dist;
                        found = 1; break;
                    }
                }
                if (!found) {
                    for (int r = 2; r > 0; r--) mmatch.rep[r] = mmatch.rep[r-1];
                    mmatch.rep[0] = dist;
                    if (mmatch.rep_count < 3) mmatch.rep_count++;
                }
            }
        }

        /* LZP prediction: order-4 context → predicted byte */
        uint32_t lzp_hash = hfin(hmix(hmix(hmix(hmix(0x9e3779b9,
                    pb[3]), pb[2]), pb[1]), pb[0]));
        uint32_t lzp_idx = lzp_hash & (LZP_SIZE - 1);
        uint8_t lzp_tag = (uint8_t)(lzp_hash >> 16);
        uint16_t lzp_entry = lzp_table[lzp_idx];
        int lzp_valid = (i >= 4 && lzp_entry != 0 && (lzp_entry >> 8) == lzp_tag);
        uint8_t lzp_pred = lzp_valid ? (lzp_entry & 0xFF) : 0;

        int bit_ctx = 1;
        uint8_t byte_val = 0;
        /* Stride model vars (hoisted for update access) */
        uint32_t stride_idx = 0;
        uint8_t stride_tag = 0;
        uint8_t stride_st = 0;
        int stride_active = 0;

        for (int bp = 7; bp >= 0; bp--) {
            uint32_t bctx = (uint32_t)bit_ctx;

            /* Hash indices for each context model */
            uint32_t idx[NUM_CM];
            uint32_t full_hash[NUM_CM];
            full_hash[0] = bctx;
            idx[0] = bctx & models[0].mask;  /* order-0: just bit context */
            for (int m = 1; m < NUM_CM; m++) {
                full_hash[m] = ctx[m] ^ (bctx * (987654321u + m * 111111111u));
                idx[m] = full_hash[m] & models[m].mask;
            }

            /* State lookups and predictions */
            uint8_t st[NUM_CM];
            int16_t stv[NUM_MODELS];
            for (int m = 0; m < NUM_CM; m++) {
                if (models[m].tagged) {
                    uint16_t entry = models[m].table16[idx[m]];
                    uint8_t tag = (uint8_t)(full_hash[m] >> 16);
                    st[m] = (entry == 0 || (entry >> 8) == tag) ? (entry & 0xFF) : 0;
                } else {
                    st[m] = models[m].table8[idx[m]];
                }
                stv[m] = state_stretch[st[m]];
            }

            /* Match model prediction */
            stv[NUM_CM] = 0;
            if (mmatch.match_active && mmatch.match_pos < bl) {
                int matched_bit = (blk[mmatch.match_pos] >> bp) & 1;
                int conf = mmatch.match_len;
                if (conf > 15) conf = 15;
                int p_match;
                int mcdiv = is_text_block ? 20 : MATCH_CONF_DIV;
                if (matched_bit)
                    p_match = PROB_HALF + (PROB_HALF - 32) * (conf + 1) / mcdiv;
                else
                    p_match = PROB_HALF - (PROB_HALF - 32) * (conf + 1) / mcdiv;
                if (p_match < 1) p_match = 1;
                if (p_match > PROB_SCALE - 1) p_match = PROB_SCALE - 1;
                stv[NUM_CM] = stretch((uint16_t)p_match);
            }

            /* ICM prediction (using order-6 context + bit_ctx) */
            uint32_t icm_h = hfin(ctx[5] ^ (bctx * 0xCC9E2D51u)) & (ICM_TAB_SIZE - 1);
            uint8_t icm_si = icm_table[icm_h];
            stv[NUM_CM + 1] = stretch(icm_ctr[icm_si]);

            /* Delta model: linear prediction */
            uint8_t delta_pred = (uint8_t)((int)pb[0] * 2 - (int)pb[1]);
            uint32_t delta_full = hfin(hmix(hmix(0, delta_pred), bctx));
            uint32_t delta_idx = delta_full & (DELTA_TAB_SIZE - 1);
            uint8_t delta_tag = (uint8_t)(delta_full >> 16);
            uint16_t delta_entry = delta_table[delta_idx];
            uint8_t delta_st = (delta_entry == 0 || (delta_entry >> 8) == delta_tag)
                                ? (delta_entry & 0xFF) : 0;
            stv[NUM_CM + 2] = state_stretch[delta_st];

            /* Stride/Record model */
            stv[NUM_CM + 3] = 0;
            stride_active = 0;
            if (stride > 0 && i >= (uint32_t)stride) {
                uint8_t stride_byte = blk[i - stride];
                uint32_t stride_full = hfin(hmix(hmix(0, stride_byte), bctx));
                stride_idx = stride_full & (STRIDE_TAB_SIZE - 1);
                stride_tag = (uint8_t)(stride_full >> 16);
                uint16_t stride_entry = stride_table[stride_idx];
                stride_st = (stride_entry == 0 || (stride_entry >> 8) == stride_tag)
                             ? (stride_entry & 0xFF) : 0;
                stv[NUM_CM + 3] = state_stretch[stride_st];
                stride_active = 1;
            }

            /* Word predictor (shadow dictionary): predicts next word's bytes */
            stv[NUM_CM + 4] = 0;
            if (wpred_active && wpred_pos < wpred_len) {
                uint8_t pred_byte = wpred_word[wpred_pos];
                int pred_bit = (pred_byte >> bp) & 1;
                int conf = wpred_pos + WPRED_BASE_CONF;
                if (conf > 15) conf = 15;
                int p_wp;
                if (pred_bit)
                    p_wp = PROB_HALF + (PROB_HALF - 32) * (conf + 1) / 16;
                else
                    p_wp = PROB_HALF - (PROB_HALF - 32) * (conf + 1) / 16;
                if (p_wp < 1) p_wp = 1;
                if (p_wp > PROB_SCALE - 1) p_wp = PROB_SCALE - 1;
                stv[NUM_CM + 4] = stretch((uint16_t)p_wp);
            }

            /* LZP prediction */
            stv[NUM_CM + 5] = 0;
            if (lzp_valid) {
                int lzp_bit = (lzp_pred >> bp) & 1;
                int p_lzp = lzp_bit ? (PROB_SCALE * 3 / 4) : (PROB_SCALE / 4);
                stv[NUM_CM + 5] = stretch((uint16_t)p_lzp);
            }

            /* Gate noise models based on block type */
            if (is_text_block) {
                for (int m = 9; m < 9 + NUM_SPARSE; m++) stv[m] = 0;
                stv[NUM_CM + 2] = 0;  /* delta */
                stv[NUM_CM + 3] = 0;  /* stride */
            } else {
                stv[9 + NUM_SPARSE] = 0;   /* prev_word */
                stv[10 + NUM_SPARSE] = 0;  /* word_trigram */
                stv[NUM_CM + 4] = 0;       /* word predictor */
            }

            /* Ablation: zero out model families by bitmask */
            if (ablate_mask) {
                if (ablate_mask & (1u<<1))  for (int m = 0; m <= 4; m++) stv[m] = 0;
                if (ablate_mask & (1u<<2))  { for (int m = 5; m <= 7; m++) stv[m] = 0; stv[NUM_CM-1] = 0; }
                if (ablate_mask & (1u<<3))  for (int m = 9; m < 9+NUM_SPARSE; m++) stv[m] = 0;
                if (ablate_mask & (1u<<4))  { stv[8]=0; stv[9+NUM_SPARSE]=0; stv[10+NUM_SPARSE]=0; }
                if (ablate_mask & (1u<<5))  stv[NUM_CM] = 0;
                if (ablate_mask & (1u<<6))  stv[NUM_CM+1] = 0;
                if (ablate_mask & (1u<<7))  stv[NUM_CM+2] = 0;
                if (ablate_mask & (1u<<8))  stv[NUM_CM+3] = 0;
                if (ablate_mask & (1u<<9))  stv[NUM_CM+4] = 0;
                if (ablate_mask & (1u<<10)) stv[NUM_CM+5] = 0;
            }

            /* Mixer 1: combine all predictions */
            int mix_ctx = ((7 - bp) << 4) | (pb[0] >> 4);
            mix_ctx += (pb[1] >> 6) * 128;
            if (is_text > 0) mix_ctx += 512;
            if (mmatch.match_active) {
                int ml = mmatch.match_len;
                int mb = (ml < 8) ? 1 : (ml < 16) ? 2 : 3;
                mix_ctx += mb * 1024;
            }
            MixerSlot *ms = &mixer[mix_ctx];
            int32_t sum = 0;
            for (int m = 0; m < NUM_MODELS; m++)
                sum += (int32_t)ms->w[m] * stv[m];
            int32_t logit_mix1 = sum / MIX_W_SCALE;

            /* Mixer 2: different context (bit_ctx-based) */
            int mix2_ctx = (bit_ctx - 1) & (MIXER2_SLOTS - 1);
            MixerSlot *ms2 = &mixer2[mix2_ctx];
            int32_t sum2 = 0;
            for (int m = 0; m < NUM_MODELS; m++)
                sum2 += (int32_t)ms2->w[m] * stv[m];
            int32_t logit_mix2 = sum2 / MIX_W_SCALE;

            /* Average the two mixers in logit domain */
            int32_t logit_mix = (logit_mix1 + logit_mix2) / 2;
            uint16_t p_mix = squash(logit_mix);

            /* ISSE pipeline: serial APM chain with blend stabilization */
            /* Stage 1: prev_byte context (serial) */
            int isse1_ctx = pb[0];
            uint16_t p_i1 = apm_predict(&isse1, isse1_ctx, p_mix);

            /* Stage 2: bp + partial bits (serial from stage 1) */
            int isse2_ctx = ((7 - bp) | ((bit_ctx & 0x1F) << 3)) & 0xFF;
            uint16_t p_i2 = apm_predict(&isse2, isse2_ctx, p_i1);

            /* Stage 3: prev_byte_hi4 + bp (serial from stage 1) */
            int isse3_ctx = ((pb[0] >> 4) | ((7 - bp) << 4)) & 0x7F;
            uint16_t p_i3 = apm_predict(&isse3, isse3_ctx, p_i1);

            /* Blend stage 2 & 3 outputs with mixer as anchor */
            uint16_t p_blend = ((uint32_t)p_mix * 2 + (uint32_t)p_i2 + (uint32_t)p_i3) / 4;
            if (p_blend < 1) p_blend = 1;
            if (p_blend > PROB_SCALE - 1) p_blend = PROB_SCALE - 1;

            /* Stage 5: word-aware APM (text blocks only) */
            int isse5_ctx = 0;
            uint16_t p_blend_pre5 = p_blend;
            if (is_text_block) {
                isse5_ctx = (((hfin(hmix(word_hash, prev_word_hash)) >> 17)
                            & 0x7FFF) << 3) | (7 - bp); /* 15-bit word hash + 3-bit bp */
                uint16_t p_i5 = apm_predict(&isse5, isse5_ctx, p_blend);
                p_blend = ((uint32_t)p_blend * 3 + (uint32_t)p_i5) / 4;
                if (p_blend < 1) p_blend = 1;
                if (p_blend > PROB_SCALE - 1) p_blend = PROB_SCALE - 1;
            }

            /* Stage 4: SSE — quantized_p + last3 + match_bit */
            int p_q = ((int)p_blend * 31) / PROB_SCALE;
            int last3 = (bit_ctx & 7);
            int isse4_ctx = (p_q << 3) | last3;
            if (mmatch.match_active && mmatch.match_pos < bl) {
                int mb = (blk[mmatch.match_pos] >> bp) & 1;
                isse4_ctx += (mb + 1) * 256;
            }
            isse4_ctx &= ((1 << ISSE4_CTX_BITS) - 1);
            uint16_t p_sse = apm_predict(&isse4, isse4_ctx, p_blend);

            /* Final: blend SSE with pre-SSE */
            uint16_t p_final = ((uint32_t)p_blend * 3 + (uint32_t)p_sse) / 4;
            if (p_final < 1) p_final = 1;
            if (p_final > PROB_SCALE - 1) p_final = PROB_SCALE - 1;

            /* Get the actual bit */
            int bit;
            if (mode == 0) {
                bit = (blk[i] >> bp) & 1;
                ev_emit(bit, p_final);
            } else {
                int r_idx = (*ev_idx_p) % RANS_INTERLEAVE;
                bit = rans_dec_bit(&rs[r_idx], dec_ptr, p_final);
                (*ev_idx_p)++;
            }

            /* Update state tables */
            for (int m = 0; m < NUM_CM; m++) {
                uint8_t ns = next_state[st[m]][bit];
                if (models[m].tagged) {
                    uint8_t tag = (uint8_t)(full_hash[m] >> 16);
                    models[m].table16[idx[m]] = ((uint16_t)tag << 8) | ns;
                } else {
                    models[m].table8[idx[m]] = ns;
                }
            }

            /* Update ICM */
            {
                int tgt = bit ? PROB_SCALE : 0;
                icm_ctr[icm_si] += (uint16_t)((tgt - (int)icm_ctr[icm_si]) >> ICM_RATE);
                if (icm_ctr[icm_si] < 1) icm_ctr[icm_si] = 1;
                if (icm_ctr[icm_si] > PROB_SCALE - 1) icm_ctr[icm_si] = PROB_SCALE - 1;
                icm_table[icm_h] = next_state[icm_si][bit];
            }

            /* Update delta model */
            {
                uint8_t delta_ns = next_state[delta_st][bit];
                delta_table[delta_idx] = ((uint16_t)delta_tag << 8) | delta_ns;
            }

            /* Update stride model */
            if (stride_active) {
                uint8_t stride_ns = next_state[stride_st][bit];
                stride_table[stride_idx] = ((uint16_t)stride_tag << 8) | stride_ns;
            }

            /* Update mixer weights */
            int32_t err = (bit ? PROB_SCALE : 0) - (int32_t)p_mix;
            for (int m = 0; m < NUM_MODELS; m++) {
                ms->w[m]  += (err * stv[m]) / (MIX_W_SCALE * MIX1_LR);
                ms2->w[m] += (err * stv[m]) / (MIX_W_SCALE * MIX2_LR);
            }

            /* Update ISSE stages */
            apm_update_rate(&isse1, isse1_ctx, p_mix, bit, ISSE1_LR);
            apm_update_rate(&isse2, isse2_ctx, p_i1, bit, ISSE2_LR);
            apm_update_rate(&isse3, isse3_ctx, p_i1, bit, ISSE3_LR);
            apm_update_rate(&isse4, isse4_ctx, p_blend, bit, ISSE4_LR);
            if (is_text_block)
                apm_update_rate(&isse5, isse5_ctx, p_blend_pre5, bit, ISSE5_LR);

            /* Update match model tracking */
            if (mmatch.match_active && mmatch.match_pos < bl) {
                int matched_bit = (blk[mmatch.match_pos] >> bp) & 1;
                if (matched_bit == bit) mmatch.match_len++;
                else { mmatch.match_active = 0; mmatch.match_len = 0; }
            }

            byte_val = (byte_val << 1) | bit;
            bit_ctx = (bit_ctx << 1) | bit;
        }

        if (mode == 1) blk[i] = byte_val;

        if (mmatch.match_active) mmatch.match_pos++;

        /* Update LZP table */
        lzp_table[lzp_idx] = ((uint16_t)lzp_tag << 8) | byte_val;

        /* Update match hash chain */
        if (i >= 3) {
            uint32_t mhkey = hfin(hmix(hmix(hmix(hmix(0,
                        pb[2]), pb[1]), pb[0]), byte_val));
            uint32_t mhash = mhkey & (MATCH_TAB_SIZE - 1);
            uint32_t new_pos = i + 1;
            mmatch.prev[new_pos] = mmatch.pos_table[mhash];
            mmatch.pos_table[mhash] = new_pos;
        }

        /* Shift previous byte history */
        for (int j = 15; j > 0; j--) pb[j] = pb[j-1];
        pb[0] = (mode == 0) ? blk[i] : byte_val;

        /* Update text/binary detector */
        {
            uint8_t tc = pb[0];
            if ((tc >= 32 && tc <= 126) || tc == 10 || tc == 13 || tc == 9)
                is_text += (is_text < 64) ? 1 : 0;
            else
                is_text -= (is_text > -64) ? 1 : 0;
        }

        /* Update word hash, word predictor tracking */
        uint8_t c = pb[0];
        if ((c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z') ||
            (c >= '0' && c <= '9') || c == '_') {
            uint8_t lc = (c >= 'A' && c <= 'Z') ? (c + 32) : c;
            word_hash = hmix(word_hash, lc);
            word_pos++;
            /* Track current word for predictor */
            if (cur_word_len < WPRED_MAX_WORD)
                cur_word_lc[cur_word_len] = lc;
            cur_word_len++;
            /* Verify prediction match */
            if (wpred_active && wpred_pos < wpred_len) {
                if (lc == wpred_word[wpred_pos])
                    wpred_pos++;
                else
                    wpred_active = 0;
            }
        } else {
            if (word_hash != 0) {
                /* Store completed word in both tables (before hash update) */
                if (cur_word_len >= 3 && cur_word_len <= WPRED_MAX_WORD) {
                    /* Unigram: prev_word → this_word */
                    if (prev_word_hash != 0) {
                        uint32_t wkey = prev_word_hash;
                        uint32_t widx = wkey & (WPRED_SIZE - 1);
                        uint16_t wtag = (uint16_t)(wkey >> 16);
                        wpred_table[widx].tag = wtag;
                        wpred_table[widx].len = (uint8_t)cur_word_len;
                        memcpy(wpred_table[widx].word, cur_word_lc, cur_word_len);
                    }
                    /* Bigram: (prev_prev + prev) → this_word */
                    if (prev_word_hash != 0 && prev_prev_word_hash != 0) {
                        uint32_t wkey = hfin(hmix(prev_prev_word_hash, prev_word_hash));
                        uint32_t widx = wkey & (WPRED_SIZE - 1);
                        uint16_t wtag = (uint16_t)(wkey >> 16);
                        wpred_table2[widx].tag = wtag;
                        wpred_table2[widx].len = (uint8_t)cur_word_len;
                        memcpy(wpred_table2[widx].word, cur_word_lc, cur_word_len);
                    }
                }
                /* Track outcome of previous prediction */
                if (profile_mode && wpred_was_active) {
                    int outcome; /* 0=correct, 1=partial, 2=wrong */
                    if (wpred_pos >= wpred_len && cur_word_len == wpred_len)
                        { wp_word_correct++; outcome = 0; }
                    else if (wpred_pos > 0)
                        { wp_word_partial++; outcome = 1; }
                    else
                        { wp_word_wrong++; outcome = 2; }
                    if (wpred_from_bigram) {
                        if (outcome==0) wp_bi_correct++;
                        else if (outcome==1) wp_bi_partial++;
                        else wp_bi_wrong++;
                    } else {
                        if (outcome==0) wp_uni_correct++;
                        else if (outcome==1) wp_uni_partial++;
                        else wp_uni_wrong++;
                    }
                }

                prev_prev_word_hash = prev_word_hash;
                prev_word_hash = hfin(word_hash);
                /* Look up prediction: bigram first, fall back to unigram */
                wpred_active = 0;
                wpred_was_active = 0;
                /* Bigram lookup: (prev_prev + prev) → next_word */
                if (prev_prev_word_hash != 0) {
                    uint32_t wkey = hfin(hmix(prev_prev_word_hash, prev_word_hash));
                    uint32_t widx = wkey & (WPRED_SIZE - 1);
                    uint16_t wtag = (uint16_t)(wkey >> 16);
                    if (profile_mode) wp_lookups++;
                    if (wpred_table2[widx].len > 0 && wpred_table2[widx].tag == wtag) {
                        wpred_len = wpred_table2[widx].len;
                        memcpy(wpred_word, wpred_table2[widx].word, wpred_len);
                        wpred_pos = 0;
                        wpred_active = 1;
                        wpred_was_active = 1;
                        wpred_from_bigram = 1;
                        if (profile_mode) wp_bigram_hit++;
                    } else if (profile_mode && wpred_table2[widx].len > 0) {
                        wp_bigram_collision++;
                    }
                }
                /* Unigram fallback: prev → next_word */
                if (!wpred_active) {
                    uint32_t wkey = prev_word_hash;
                    uint32_t widx = wkey & (WPRED_SIZE - 1);
                    uint16_t wtag = (uint16_t)(wkey >> 16);
                    if (profile_mode && prev_prev_word_hash == 0) wp_lookups++;
                    if (wpred_table[widx].len > 0 && wpred_table[widx].tag == wtag) {
                        wpred_len = wpred_table[widx].len;
                        memcpy(wpred_word, wpred_table[widx].word, wpred_len);
                        wpred_pos = 0;
                        wpred_active = 1;
                        wpred_was_active = 1;
                        wpred_from_bigram = 0;
                        if (profile_mode) wp_unigram_hit++;
                    } else if (profile_mode) {
                        if (wpred_table[widx].len > 0) wp_unigram_collision++;
                        else wp_miss++;
                    }
                }
            }
            cur_word_len = 0;
            word_hash = 0;
            word_pos = 0;
        }
    }

    if (profile_mode && wp_lookups > 0) {
        uint32_t wp_hits = wp_bigram_hit + wp_unigram_hit;
        uint32_t wp_preds = wp_word_correct + wp_word_partial + wp_word_wrong;
        fprintf(stderr, "  WPRED STATS:\n");
        fprintf(stderr, "    lookups:     %u\n", wp_lookups);
        fprintf(stderr, "    bigram hit:  %u (%.1f%%)\n", wp_bigram_hit, 100.0*wp_bigram_hit/wp_lookups);
        fprintf(stderr, "    unigram hit: %u (%.1f%%)\n", wp_unigram_hit, 100.0*wp_unigram_hit/wp_lookups);
        fprintf(stderr, "    miss:        %u (%.1f%%)\n", wp_miss, 100.0*wp_miss/wp_lookups);
        fprintf(stderr, "    bigram coll: %u (%.1f%%)\n", wp_bigram_collision, 100.0*wp_bigram_collision/wp_lookups);
        fprintf(stderr, "    unigram col: %u (%.1f%%)\n", wp_unigram_collision, 100.0*wp_unigram_collision/wp_lookups);
        if (wp_preds > 0) {
            fprintf(stderr, "    predictions: %u (correct=%u %.1f%%, partial=%u %.1f%%, wrong=%u %.1f%%)\n",
                    wp_preds, wp_word_correct, 100.0*wp_word_correct/wp_preds,
                    wp_word_partial, 100.0*wp_word_partial/wp_preds,
                    wp_word_wrong, 100.0*wp_word_wrong/wp_preds);
            uint32_t bi_total = wp_bi_correct + wp_bi_partial + wp_bi_wrong;
            uint32_t uni_total = wp_uni_correct + wp_uni_partial + wp_uni_wrong;
            if (bi_total > 0)
                fprintf(stderr, "    bigram pred: %u (correct=%.1f%%, partial=%.1f%%, wrong=%.1f%%)\n",
                        bi_total, 100.0*wp_bi_correct/bi_total,
                        100.0*wp_bi_partial/bi_total, 100.0*wp_bi_wrong/bi_total);
            if (uni_total > 0)
                fprintf(stderr, "    unigram pred: %u (correct=%.1f%%, partial=%.1f%%, wrong=%.1f%%)\n",
                        uni_total, 100.0*wp_uni_correct/uni_total,
                        100.0*wp_uni_partial/uni_total, 100.0*wp_uni_wrong/uni_total);
        }
    }

    match_free();
}

/* ---- Compress ---- */
static int do_compress(const char *inp, const char *outp) {
    FILE *fi = fopen(inp, "rb");
    if (!fi) { perror("open"); return 1; }
    fseek(fi, 0, SEEK_END);
    uint32_t orig = (uint32_t)ftell(fi);
    fseek(fi, 0, SEEK_SET);
    uint8_t *buf = (uint8_t *)malloc(orig);
    if (!buf) { fprintf(stderr, "OOM\n"); return 1; }
    fread(buf, 1, orig, fi);
    fclose(fi);

    /* Auto-detect file type for pre-filtering */
    uint8_t flags = 0;
    uint32_t data_len = orig;
    uint16_t img_raw_row = 0;

    if (orig >= 54 && buf[0] == 'B' && buf[1] == 'M') {
        /* BMP detected — apply image delta filter, skip BCJ */
        int bpp = buf[28] | (buf[29] << 8);
        if (bpp == 24 || bpp == 32) {
            int width = (int32_t)(buf[18] | (buf[19]<<8) | (buf[20]<<16) | (buf[21]<<24));
            int data_offset = buf[10] | (buf[11]<<8) | (buf[12]<<16) | (buf[13]<<24);
            if (width > 0) {
                int bmp_row_bytes = ((width * (bpp / 8) + 3) & ~3);
                bmp_delta_encode(buf, orig, data_offset, bmp_row_bytes);
                flags = FLAG_IMG_DELTA;
            }
        }
    } else if (orig >= 44 && buf[0]=='R' && buf[1]=='I' && buf[2]=='F' && buf[3]=='F'
               && buf[8]=='W' && buf[9]=='A' && buf[10]=='V' && buf[11]=='E') {
        /* WAV detected — apply per-channel delta filter, skip BCJ */
        int audio_fmt = buf[20] | (buf[21] << 8);
        int num_channels = buf[22] | (buf[23] << 8);
        int bps = buf[34] | (buf[35] << 8);
        if (audio_fmt == 1 && bps == 16 && num_channels >= 1 && num_channels <= 8) {
            int data_offset = wav_find_data(buf, orig);
            if (data_offset > 0) {
                wav_delta_encode(buf, orig, data_offset, num_channels);
                flags = FLAG_WAV_DELTA;
            }
        }
    } else {
        /* Try raw image detection before falling back to text/binary */
        int img_row = detect_image_geometry(buf, orig);
        if (img_row > 0) {
            img_raw_row = (uint16_t)img_row;
            raw_delta_encode(buf, orig, img_row);
            flags = FLAG_IMG_RAW;
        } else {
            /* Detect text vs binary */
            int sample = orig < 4096 ? (int)orig : 4096;
            int non_text = 0;
            for (int i = 0; i < sample; i++) {
                uint8_t c = buf[i];
                if (c < 9 || (c > 13 && c < 32) || c == 127) non_text++;
            }
            bcj_encode(buf, orig);
            flags = FLAG_BCJ;
            if (non_text * 20 < sample) flags |= FLAG_TEXT;
        }
    }

    FILE *fo = fopen(outp, "wb");
    fwrite(MAGIC, 1, 4, fo);
    wr32(fo, orig);
    wr32(fo, BLOCK_SIZE);
    fwrite(&flags, 1, 1, fo);
    if (flags & FLAG_TEXT_MTF) wr32(fo, data_len);
    if (flags & FLAG_IMG_RAW) fwrite(&img_raw_row, 2, 1, fo);

    /* Order-N models (untagged for small tables, tagged for large) */
    cm_init(&models[0],  CM_O0_BITS,     0);
    cm_init(&models[1],  CM_O1_BITS,     0);
    cm_init(&models[2],  CM_O2_BITS,     0);
    cm_init(&models[3],  CM_O3_BITS,     0);
    cm_init(&models[4],  CM_O4_BITS,     1);
    cm_init(&models[5],  CM_O6_BITS,     1);
    cm_init(&models[6],  CM_O8_BITS,     1);
    cm_init(&models[7],  CM_O12_BITS,    1);
    cm_init(&models[8],  CM_WORD_BITS,   1);
    /* Sparse models (all same size — add more here to improve compression) */
    for (int s = 0; s < NUM_SPARSE; s++)
        cm_init(&models[9 + s], CM_SPARSE_BITS, 1);
    /* Previous-word context model */
    cm_init(&models[9 + NUM_SPARSE], CM_WORD_BITS, 1);
    /* Word trigram context model */
    cm_init(&models[10 + NUM_SPARSE], CM_WORD_BITS, 1);
    /* Order-16 (always last) */
    cm_init(&models[NUM_CM - 1], CM_O16_BITS, 1);

    ev_buf = (uint16_t *)malloc((size_t)MAX_BIT_EVENTS * sizeof(uint16_t));
    size_t rbsz = (size_t)BLOCK_SIZE * 16 + 65536;
    uint8_t *rb = (uint8_t *)malloc(rbsz);
    if (!ev_buf || !rb) { fprintf(stderr, "OOM\n"); return 1; }

    mixer_init();
    apm_init(&isse1, ISSE1_CTX_BITS, ISSE1_PROB_BITS);
    apm_init(&isse2, ISSE2_CTX_BITS, ISSE2_PROB_BITS);
    apm_init(&isse3, ISSE3_CTX_BITS, ISSE3_PROB_BITS);
    apm_init(&isse4, ISSE4_CTX_BITS, ISSE4_PROB_BITS);
    apm_init(&isse5, ISSE5_CTX_BITS, ISSE5_PROB_BITS);
    icm_init();
    delta_init();
    stride_init();
    wpred_init();
    lzp_init();

    int ft = (flags & (FLAG_TEXT_MTF | FLAG_TEXT)) ? 1 : 0;
    uint32_t pos = 0;
    int bn = 0;
    while (pos < data_len) {
        uint32_t bl = data_len - pos;
        if (bl > BLOCK_SIZE) bl = BLOCK_SIZE;
        uint8_t *blk = buf + pos;

        int stride = detect_stride(blk, bl);
        uint8_t stride_byte = (uint8_t)stride;

        process_block(blk, bl, 0, NULL, NULL, NULL, stride, ft);

        /* Entropy profiler: compute log-loss per 4KB region */
        if (profile_mode) {
            #define PROFILE_WINDOW 4096
            int32_t ev_i = 0;
            double block_total = 0.0;
            fprintf(stderr, "PROFILE block %d (offset %u, %u bytes):\n", bn+1, pos, bl);
            fprintf(stderr, "  %8s %10s %10s\n", "offset", "bits", "bpb");
            for (uint32_t wi = 0; wi < bl; wi += PROFILE_WINDOW) {
                uint32_t wend = wi + PROFILE_WINDOW;
                if (wend > bl) wend = bl;
                uint32_t wbytes = wend - wi;
                double window_bits = 0.0;
                int32_t window_events = (int32_t)wbytes * 8;
                for (int32_t e = 0; e < window_events && ev_i < ev_cnt; e++, ev_i++) {
                    int bit = (ev_buf[ev_i] >> 15) & 1;
                    uint16_t prob = ev_buf[ev_i] & 0x0FFF;
                    if (prob < 1) prob = 1;
                    if (prob > PROB_SCALE - 1) prob = PROB_SCALE - 1;
                    double p = bit ? (double)prob / PROB_SCALE
                                   : (double)(PROB_SCALE - prob) / PROB_SCALE;
                    if (p < 1e-9) p = 1e-9;
                    window_bits += -log2(p);
                }
                double bpb = window_bits / wbytes;
                block_total += window_bits;
                fprintf(stderr, "  %8u %10.1f %10.4f", pos + wi, window_bits, bpb);
                if (bpb > 7.0) fprintf(stderr, "  *** HIGH");
                else if (bpb < 1.0) fprintf(stderr, "  +++ LOW");
                fprintf(stderr, "\n");
            }
            fprintf(stderr, "  BLOCK TOTAL: %.1f bits (%.4f bpb)\n",
                    block_total, block_total / bl);

            /* Dump average mixer weights per model (text slots only: ctx >= 512) */
            const char *model_names[] = {
                "o0","o1","o2","o3","o4","o6","o8","o12","word",
                "sp0","sp1","sp2","sp3","sp4","sp5","sp6","sp7","sp8","sp9",
                "sp10","sp11","sp12","sp13","sp14","sp15","sp16","sp17","sp18","sp19",
                "sp20","sp21","sp22","sp23","sp24","sp25","sp26","sp27","sp28","sp29",
                "sp30","sp31","sp32","sp33",
                "prev_w","w_tri","o16",
                "match","icm","delta","stride","wpred","lzp"
            };
            fprintf(stderr, "  MIXER WEIGHTS (avg |w| across text slots 512-1023):\n");
            for (int m = 0; m < NUM_MODELS; m++) {
                int64_t sum_abs = 0;
                int cnt = 0;
                for (int s = 512; s < 1024 && s < MIXER_SLOTS; s++) {
                    int32_t w = mixer[s].w[m];
                    sum_abs += (w < 0 ? -w : w);
                    cnt++;
                }
                double avg = (double)sum_abs / cnt;
                if (avg > 1.0)
                    fprintf(stderr, "    %-8s %8.1f\n", model_names[m], avg);
            }
        }

        uint8_t *rp = rb + rbsz;
        uint32_t rs_enc[RANS_INTERLEAVE];
        for (int r = 0; r < RANS_INTERLEAVE; r++) rs_enc[r] = RANS_L;

        for (int32_t e = ev_cnt - 1; e >= 0; e--) {
            int r = e % RANS_INTERLEAVE;
            int bit = (ev_buf[e] >> 15) & 1;
            uint16_t prob = ev_buf[e] & 0x0FFF;
            rs_enc[r] = rans_enc_bit(rs_enc[r], &rp, bit, prob);
        }

        for (int r = RANS_INTERLEAVE - 1; r >= 0; r--)
            rans_flush(rs_enc[r], &rp);

        uint32_t csz = (uint32_t)(rb + rbsz - rp);
        wr32(fo, csz);
        wr32(fo, (uint32_t)ev_cnt);
        fwrite(&stride_byte, 1, 1, fo);
        fwrite(rp, 1, csz, fo);

        bn++;
        fprintf(stderr, "Block %d: %u -> %u (stride=%d)\n", bn, bl, csz, stride);
        pos += bl;
    }

    fclose(fo);
    free(buf);
    for (int m = 0; m < NUM_CM; m++) cm_free(&models[m]);

    apm_free(&isse1);
    apm_free(&isse2);
    apm_free(&isse3);
    apm_free(&isse4);
    apm_free(&isse5);
    icm_free();
    delta_free();
    stride_free();
    wpred_free();
    lzp_free();
    free(ev_buf);
    free(rb);

    fi = fopen(outp, "rb");
    fseek(fi, 0, SEEK_END);
    fprintf(stderr, "Total: %u -> %ld\n", orig, ftell(fi));
    fclose(fi);
    return 0;
}

/* ---- Decompress ---- */
static int do_decompress(const char *inp, const char *outp) {
    FILE *fi = fopen(inp, "rb");
    if (!fi) { perror("open"); return 1; }

    uint8_t mg[4]; fread(mg, 1, 4, fi);
    if (memcmp(mg, MAGIC, 4)) { fprintf(stderr, "Bad magic\n"); return 1; }

    uint32_t orig = rd32(fi);
    uint32_t block_size = rd32(fi);
    uint8_t flags;
    fread(&flags, 1, 1, fi);
    uint32_t data_len = orig;
    if (flags & FLAG_TEXT_MTF) data_len = rd32(fi);
    uint16_t img_raw_row = 0;
    if (flags & FLAG_IMG_RAW) fread(&img_raw_row, 2, 1, fi);

    uint8_t *obuf = (uint8_t *)calloc(data_len, 1);
    if (!obuf) { fprintf(stderr, "OOM\n"); return 1; }

    /* Order-N models (untagged for small tables, tagged for large) */
    cm_init(&models[0],  CM_O0_BITS,     0);
    cm_init(&models[1],  CM_O1_BITS,     0);
    cm_init(&models[2],  CM_O2_BITS,     0);
    cm_init(&models[3],  CM_O3_BITS,     0);
    cm_init(&models[4],  CM_O4_BITS,     1);
    cm_init(&models[5],  CM_O6_BITS,     1);
    cm_init(&models[6],  CM_O8_BITS,     1);
    cm_init(&models[7],  CM_O12_BITS,    1);
    cm_init(&models[8],  CM_WORD_BITS,   1);
    /* Sparse models (all same size — add more here to improve compression) */
    for (int s = 0; s < NUM_SPARSE; s++)
        cm_init(&models[9 + s], CM_SPARSE_BITS, 1);
    /* Previous-word context model */
    cm_init(&models[9 + NUM_SPARSE], CM_WORD_BITS, 1);
    /* Word trigram context model */
    cm_init(&models[10 + NUM_SPARSE], CM_WORD_BITS, 1);
    /* Order-16 (always last) */
    cm_init(&models[NUM_CM - 1], CM_O16_BITS, 1);

    size_t rbsz = (size_t)block_size * 16 + 65536;
    uint8_t *rb = (uint8_t *)malloc(rbsz);
    if (!rb) { fprintf(stderr, "OOM\n"); return 1; }

    mixer_init();
    apm_init(&isse1, ISSE1_CTX_BITS, ISSE1_PROB_BITS);
    apm_init(&isse2, ISSE2_CTX_BITS, ISSE2_PROB_BITS);
    apm_init(&isse3, ISSE3_CTX_BITS, ISSE3_PROB_BITS);
    apm_init(&isse4, ISSE4_CTX_BITS, ISSE4_PROB_BITS);
    apm_init(&isse5, ISSE5_CTX_BITS, ISSE5_PROB_BITS);
    icm_init();
    delta_init();
    stride_init();
    wpred_init();
    lzp_init();

    int ft = (flags & (FLAG_TEXT_MTF | FLAG_TEXT)) ? 1 : 0;
    uint32_t pos = 0;
    int bn = 0;
    while (pos < data_len) {
        uint32_t bl = data_len - pos;
        if (bl > block_size) bl = block_size;
        uint8_t *blk = obuf + pos;

        uint32_t csz = rd32(fi);
        uint32_t n_events = rd32(fi);
        (void)n_events;
        uint8_t stride_byte;
        fread(&stride_byte, 1, 1, fi);
        int stride = (int)stride_byte;
        fread(rb, 1, csz, fi);
        memset(rb + csz, 0, 64);

        const uint8_t *rp = rb;
        uint32_t rs_dec[RANS_INTERLEAVE];
        for (int r = 0; r < RANS_INTERLEAVE; r++)
            rs_dec[r] = rans_dec_init(&rp);

        int32_t ev_idx = 0;
        process_block(blk, bl, 1, rs_dec, &rp, &ev_idx, stride, ft);

        bn++;
        fprintf(stderr, "Block %d: %u bytes (stride=%d)\n", bn, bl, stride);
        pos += bl;
    }
    fclose(fi);

    if (flags & FLAG_TEXT_MTF) {
        uint8_t *decoded = (uint8_t *)malloc(orig);
        if (!decoded) { fprintf(stderr, "OOM\n"); return 1; }
        text_mtf_decode(obuf, data_len, decoded, orig);
        free(obuf);
        obuf = decoded;
    }
    if (flags & FLAG_BCJ) bcj_decode(obuf, orig);
    if (flags & FLAG_IMG_DELTA) {
        int bpp = obuf[28] | (obuf[29] << 8);
        int width = (int32_t)(obuf[18] | (obuf[19]<<8) | (obuf[20]<<16) | (obuf[21]<<24));
        int data_offset = obuf[10] | (obuf[11]<<8) | (obuf[12]<<16) | (obuf[13]<<24);
        int bmp_row_bytes = ((width * (bpp / 8) + 3) & ~3);
        bmp_delta_decode(obuf, orig, data_offset, bmp_row_bytes);
    }
    if (flags & FLAG_WAV_DELTA) {
        int num_channels = obuf[22] | (obuf[23] << 8);
        int data_offset = wav_find_data(obuf, orig);
        if (data_offset > 0)
            wav_delta_decode(obuf, orig, data_offset, num_channels);
    }
    if (flags & FLAG_IMG_RAW) {
        raw_delta_decode(obuf, orig, (int)img_raw_row);
    }

    FILE *fo = fopen(outp, "wb");
    fwrite(obuf, 1, orig, fo);
    fclose(fo);

    free(obuf);
    for (int m = 0; m < NUM_CM; m++) cm_free(&models[m]);

    apm_free(&isse1);
    apm_free(&isse2);
    apm_free(&isse3);
    apm_free(&isse4);
    apm_free(&isse5);
    icm_free();
    delta_free();
    stride_free();
    wpred_free();
    lzp_free();
    free(rb);
    fprintf(stderr, "Decompressed %u bytes\n", orig);
    return 0;
}

static int parse_ablate(const char *s) {
    const char *names[] = {
        "order_low","order_high","sparse","linguistic",
        "match","icm","delta","stride","wpred","lzp"
    };
    for (int i = 0; i < 10; i++)
        if (strcmp(s, names[i]) == 0) return i + 1;
    int v = atoi(s);
    return (v >= 1 && v <= 10) ? v : -1;
}

int main(int argc, char **argv) {
    if (argc < 4) {
        fprintf(stderr, "Usage: %s c|d in out [--ablate=FAMILY]... [--profile]\n", argv[0]);
        fprintf(stderr, "Families: order_low order_high sparse linguistic "
                        "match icm delta stride wpred lzp\n");
        return 1;
    }

    for (int a = 4; a < argc; a++) {
        if (strncmp(argv[a], "--ablate=", 9) == 0) {
            int fam = parse_ablate(argv[a] + 9);
            if (fam < 0) {
                fprintf(stderr, "Unknown ablation family: %s\n", argv[a] + 9);
                return 1;
            }
            ablate_mask |= (1u << fam);
            fprintf(stderr, "Ablation: silencing %s\n", argv[a] + 9);
        } else if (strcmp(argv[a], "--profile") == 0) {
            profile_mode = 1;
        }
    }

    init_stretch_squash();
    init_state_table();

    if (argv[1][0] == 'c') return do_compress(argv[2], argv[3]);
    if (argv[1][0] == 'd') return do_decompress(argv[2], argv[3]);
    fprintf(stderr, "Unknown mode\n");
    return 1;
}