zxc_huffman.c

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-3-Clause
/*
 * ZXC - High-performance lossless compression
 *
 * Copyright (c) Bertrand Lebonnois and contributors.
 *
 */
 
// Function Multi-Versioning Support
// With ZXC_FUNCTION_SUFFIX defined (e.g. _avx2), each variant TU gets its own
// copy under a unique symbol, which the runtime dispatcher routes to. The
// defines precede zxc_internal.h so the header's prototypes take the same
// suffix as the definitions below.
#ifdef ZXC_FUNCTION_SUFFIX
#define ZXC_CAT_IMPL(x, y) x##y
#define ZXC_CAT(x, y) ZXC_CAT_IMPL(x, y)
#define zxc_huf_build_code_lengths ZXC_CAT(zxc_huf_build_code_lengths, ZXC_FUNCTION_SUFFIX)
#define zxc_huf_pack_lengths ZXC_CAT(zxc_huf_pack_lengths, ZXC_FUNCTION_SUFFIX)
#define zxc_huf_unpack_lengths ZXC_CAT(zxc_huf_unpack_lengths, ZXC_FUNCTION_SUFFIX)
#define zxc_huf_calc_size ZXC_CAT(zxc_huf_calc_size, ZXC_FUNCTION_SUFFIX)
#define zxc_huf_calc_size_dict ZXC_CAT(zxc_huf_calc_size_dict, ZXC_FUNCTION_SUFFIX)
#define zxc_huf_encode_section ZXC_CAT(zxc_huf_encode_section, ZXC_FUNCTION_SUFFIX)
#define zxc_huf_encode_section_dict ZXC_CAT(zxc_huf_encode_section_dict, ZXC_FUNCTION_SUFFIX)
#define zxc_huf_decode_section ZXC_CAT(zxc_huf_decode_section, ZXC_FUNCTION_SUFFIX)
#define zxc_huf_decode_section_dict ZXC_CAT(zxc_huf_decode_section_dict, ZXC_FUNCTION_SUFFIX)
#endif
 
// Mark the primary variant (only _default, or a no-suffix build) so ISA-
// independent cold code compiles once, not in every per-ISA copy. Keyed off the
// suffix value, so every build gets it with no extra flag.
#ifdef ZXC_FUNCTION_SUFFIX
#define ZXC_PRIMARY__default 1
#define ZXC_PRIMARY__avx2 0
#define ZXC_PRIMARY__avx512 0
#define ZXC_PRIMARY__neon32 0
#define ZXC_PRIMARY_CAT_(a, b) a##b
#define ZXC_PRIMARY_CAT(a, b) ZXC_PRIMARY_CAT_(a, b)
#if ZXC_PRIMARY_CAT(ZXC_PRIMARY_, ZXC_FUNCTION_SUFFIX)
#define ZXC_VARIANT_PRIMARY 1
#endif
#else
#define ZXC_VARIANT_PRIMARY 1
#endif
 
#include "../../include/zxc_error.h"
#include "zxc_internal.h"
#include "zxc_pivco_tables.h"
 
// ===========================================================================
// Length-limited Huffman: boundary package-merge
// ===========================================================================
//
// Builds optimal length-limited Huffman code lengths (max length
// ZXC_HUF_MAX_CODE_LEN_ULTRA) on 256-symbol alphabets. Package-merge is run for
// ZXC_HUF_MAX_CODE_LEN_ULTRA levels; each level holds up to 2N items (leaves +
// paired packages). Selection of the cheapest 2N - 2 items at level
// ZXC_HUF_MAX_CODE_LEN_ULTRA gives the appearance count of each leaf, which is
// its code length.
 
typedef zxc_huf_pm_item_t pm_item_t;
 
typedef zxc_huf_pm_leaf_t pm_leaf_t;
 
typedef zxc_huf_pm_frame_t frame_t;
 
static void pm_leaves_sort(pm_leaf_t* RESTRICT leaves, const int n) {
    // One bucket per possible value of floor(log2(weight)) for a 32-bit
    // weight, i.e. 32 buckets.
    enum { NUM_BUCKETS = 32 };
    int count[NUM_BUCKETS];
    int offset[NUM_BUCKETS + 1]; /* +1 sentinel = n, avoids end-of-bucket branch. */
    uint8_t bucket_of[ZXC_HUF_NUM_SYMBOLS];
    pm_leaf_t tmp[ZXC_HUF_NUM_SYMBOLS];
 
    ZXC_MEMSET(count, 0, sizeof(count));
    for (int i = 0; i < n; i++) {
        const unsigned b = zxc_log2_u32(leaves[i].w);
        bucket_of[i] = (uint8_t)b;
        count[b]++;
    }
 
    int acc = 0;
    for (int b = 0; b < NUM_BUCKETS; b++) {
        offset[b] = acc;
        acc += count[b];
    }
    offset[NUM_BUCKETS] = n;
 
    int pos[NUM_BUCKETS];
    ZXC_MEMCPY(pos, offset, sizeof(pos));
    for (int i = 0; i < n; i++) {
        tmp[pos[bucket_of[i]]++] = leaves[i];
    }
 
    for (int b = 0; b < NUM_BUCKETS; b++) {
        if (count[b] < 2) continue;
        const int s = offset[b];
        const int e = offset[b + 1];
        for (int i = s + 1; i < e; i++) {
            const pm_leaf_t key = tmp[i];
            int j = i - 1;
            while (j >= s && (tmp[j].w > key.w || (tmp[j].w == key.w && tmp[j].sym > key.sym))) {
                tmp[j + 1] = tmp[j];
                j--;
            }
            tmp[j + 1] = key;
        }
    }
 
    ZXC_MEMCPY(leaves, tmp, (size_t)n * sizeof(pm_leaf_t));
}
 
int zxc_huf_build_code_lengths(const uint32_t* RESTRICT freq, uint8_t* RESTRICT code_len,
                               void* RESTRICT scratch, const int max_code_len) {
    ZXC_MEMSET(code_len, 0, ZXC_HUF_NUM_SYMBOLS);
 
    pm_leaf_t leaves[ZXC_HUF_NUM_SYMBOLS];
    int n = 0;
    for (int i = 0; i < ZXC_HUF_NUM_SYMBOLS; i++) {
        if (freq[i] > 0) {
            leaves[n].w = freq[i];
            leaves[n].sym = (int16_t)i;
            n++;
        }
    }
    if (UNLIKELY(n == 0)) return ZXC_ERROR_CORRUPT_DATA;
    if (n == 1) {
        code_len[leaves[0].sym] = 1;
        return ZXC_OK;
    }
 
    pm_leaves_sort(leaves, n);
 
    // Callers pass max_code_len >= 8, and n <= 256 <= 2^max_code_len, so the
    // length limit is always feasible (a full 8-bit code covers all 256 symbols).
    const int max_per_level = 2 * n;
 
    // Working buffers: either carve from caller-provided scratch (sized for
    // the worst-case alphabet) or fall back to per-call malloc/free.
    pm_item_t* items;
    int* counts;
    frame_t* stack;
    pm_item_t* owned_items = NULL;
    int* owned_counts = NULL;
    frame_t* owned_stack = NULL;
    if (scratch) {
        uint8_t* p = (uint8_t*)scratch;
        items = (pm_item_t*)p;
        p +=
            (size_t)ZXC_HUF_MAX_CODE_LEN_ULTRA * (size_t)ZXC_HUF_PM_LEVEL_BOUND * sizeof(pm_item_t);
        p = (uint8_t*)(((uintptr_t)p + 7U) & ~(uintptr_t)7U);
        counts = (int*)p;
        ZXC_MEMSET(counts, 0, (size_t)ZXC_HUF_MAX_CODE_LEN_ULTRA * sizeof(int));
        p += (size_t)ZXC_HUF_MAX_CODE_LEN_ULTRA * sizeof(int);
        p = (uint8_t*)(((uintptr_t)p + 7U) & ~(uintptr_t)7U);
        stack = (frame_t*)p;
    } else {
        owned_items = (pm_item_t*)ZXC_MALLOC((size_t)ZXC_HUF_MAX_CODE_LEN_ULTRA *
                                             (size_t)max_per_level * sizeof(pm_item_t));
        owned_counts = (int*)ZXC_CALLOC((size_t)ZXC_HUF_MAX_CODE_LEN_ULTRA, sizeof(int));
        owned_stack = (frame_t*)ZXC_MALLOC((size_t)ZXC_HUF_MAX_CODE_LEN_ULTRA *
                                           (size_t)max_per_level * sizeof(frame_t));
        if (UNLIKELY(!owned_items || !owned_counts || !owned_stack)) {
            ZXC_FREE(owned_items);
            ZXC_FREE(owned_counts);
            ZXC_FREE(owned_stack);
            return ZXC_ERROR_MEMORY;
        }
        items = owned_items;
        counts = owned_counts;
        stack = owned_stack;
    }
#define ITEM(k, i) items[(size_t)(k) * (size_t)max_per_level + (size_t)(i)]
 
    // Level 0 (logical level 1): the leaves themselves, already sorted.
    for (int i = 0; i < n; i++) {
        ITEM(0, i).weight = leaves[i].w;
        ITEM(0, i).left = -1;
        ITEM(0, i).right = -1;
        ITEM(0, i).sym = leaves[i].sym;
    }
    counts[0] = n;
 
    // Levels 1..max_code_len-1: merge sorted leaves with sorted packages from the previous
    // level.
    for (int k = 1; k < max_code_len; k++) {
        const int prev = counts[k - 1];
        const int packs = prev / 2;
        int li = 0;
        int pi = 0;
        int n_lvl = 0;
        while (li < n || pi < packs) {
            const uint32_t wl = (li < n) ? leaves[li].w : UINT32_MAX;
            const uint32_t wp =
                (pi < packs)
                    ? (uint32_t)(ITEM(k - 1, 2 * pi).weight + ITEM(k - 1, 2 * pi + 1).weight)
                    : UINT32_MAX;
            if (wl <= wp && li < n) {
                ITEM(k, n_lvl).weight = wl;
                ITEM(k, n_lvl).left = -1;
                ITEM(k, n_lvl).right = -1;
                ITEM(k, n_lvl).sym = leaves[li].sym;
                li++;
            } else {
                ITEM(k, n_lvl).weight = wp;
                ITEM(k, n_lvl).left = (int16_t)(2 * pi);
                ITEM(k, n_lvl).right = (int16_t)(2 * pi + 1);
                ITEM(k, n_lvl).sym = -1;
                pi++;
            }
            n_lvl++;
        }
        counts[k] = n_lvl;
    }
 
    // Step 3: take first 2n-2 items at level max_code_len-1; trace back, counting leaf
    // appearances.
    int n_take = 2 * n - 2;
    if (n_take > counts[max_code_len - 1]) n_take = counts[max_code_len - 1];
 
    // Worst case stack depth: (max_code_len * n_take) frames; bounded by
    // max_code_len * 2n <= ZXC_HUF_MAX_CODE_LEN_ULTRA * 2n. `stack` was set up earlier from
    // scratch (or the local malloc fallback), sized for the ceiling.
    int sp = 0;
    for (int i = 0; i < n_take; i++) {
        stack[sp].lvl = (int8_t)(max_code_len - 1);
        stack[sp].idx = (int16_t)i;
        sp++;
    }
    while (sp > 0) {
        frame_t f = stack[--sp];
        const pm_item_t* it = &ITEM(f.lvl, f.idx);
        if (it->sym >= 0) {
            code_len[it->sym]++;
        } else {
            stack[sp].lvl = (int8_t)(f.lvl - 1);
            stack[sp].idx = it->left;
            sp++;
            stack[sp].lvl = (int8_t)(f.lvl - 1);
            stack[sp].idx = it->right;
            sp++;
        }
    }
 
    if (owned_items) {
        ZXC_FREE(owned_items);
        ZXC_FREE(owned_counts);
        ZXC_FREE(owned_stack);
    }
#undef ITEM
    return ZXC_OK;
}
 
// ===========================================================================
// Encoder-side joint flat/length nudge (idea: pivco-huffman issue #20)
// ===========================================================================
//
// Package-merge minimizes bits alone; the PivCo decoder's cost also depends on
// the SHAPE of the length histogram (pass count = max_depth + 1, and maximal
// complete subtrees collapse into single unpacks -- see the decode section).
// The functions below reshape the class counts toward power-of-two boundaries
// and shallower caps under an adoption guard, trading a bounded ratio loss for
// fewer modeled decode level-touches.
//
// Everything is exact integer arithmetic on the canonical code space. With
// leaves assigned in (length, symbol) ascending order -- the canonical rule of
// zxc_pivco_tree_build -- the depth-l leaves occupy one contiguous ALIGNED run
// [S_{l-1}, S_l) of the code space at 2^-11 slot granularity, where
// S_l = sum_{j<=l} bl_count[j] << (11-j) (each prior term is divisible by
// 2^(12-l), so runs start on even slot boundaries). The maximal dyadic (buddy)
// blocks of each run are EXACTLY the decoder's maximal flat subtrees, so any
// candidate histogram is priced without building a tree.
//
// Compiled once in the primary variant: this is ISA-independent cold encoder
// policy, and a single copy guarantees cross-ISA identical archives. For an
// A/B build without the nudge, set ZXC_HUF_NUDGE_MERGE_Q8 to 0 in
// zxc_internal.h: the guard then rejects every candidate.
#if defined(ZXC_VARIANT_PRIMARY)
 
typedef struct {
    uint64_t bits;    
    uint64_t touches; 
} zxc_huf_nudge_cost_t;
 
static int zxc_huf_nudge_classes(const uint8_t* RESTRICT code_len, uint32_t* RESTRICT blc) {
    ZXC_MEMSET(blc, 0, (ZXC_HUF_MAX_CODE_LEN_ULTRA + 1) * sizeof(uint32_t));
    int n = 0;
    for (int s = 0; s < ZXC_HUF_NUM_SYMBOLS; s++) {
        const int l = code_len[s];
        if (!l) continue;
        blc[l]++;
        n++;
    }
    return n;
}
 
static void zxc_huf_nudge_pf_canonical(const uint8_t* RESTRICT code_len,
                                       const uint32_t* RESTRICT freq, const uint32_t* RESTRICT blc,
                                       uint64_t* RESTRICT pf) {
    uint32_t pos[ZXC_HUF_MAX_CODE_LEN_ULTRA + 1];
    uint32_t acc = 0;
    for (int l = 1; l <= ZXC_HUF_MAX_CODE_LEN_ULTRA; l++) {
        pos[l] = acc;
        acc += blc[l];
    }
    uint32_t val[ZXC_HUF_NUM_SYMBOLS];
    for (int s = 0; s < ZXC_HUF_NUM_SYMBOLS; s++) {
        const int l = code_len[s];
        if (l) val[pos[l]++] = freq[s];
    }
    pf[0] = 0;
    for (uint32_t i = 0; i < acc; i++) pf[i + 1] = pf[i] + val[i];
}
 
static void zxc_huf_nudge_eval(const uint32_t* RESTRICT blc, const uint64_t* RESTRICT pf,
                               zxc_huf_nudge_cost_t* RESTRICT out) {
    enum { LU = ZXC_HUF_MAX_CODE_LEN_ULTRA };
    uint64_t bits = 0;
    uint64_t touches = 0;
    int max_len = 0;
    uint32_t S = 0;    /* code-space cursor, in 2^-LU slots */
    uint32_t base = 0; /* leaf index where the current run starts */
    for (int l = 1; l <= LU; l++) {
        if (!blc[l]) continue;
        max_len = l;
        bits += (uint64_t)l * (pf[base + blc[l]] - pf[base]);
        const uint32_t w = 1U << (LU - l);
        const uint32_t end = S + blc[l] * w;
        uint32_t x = S;
        while (x < end) {
            const uint32_t wx = x ? (x & (0U - x)) : (1U << LU);
            const uint32_t wr = 1U << zxc_log2_u32(end - x);
            const uint32_t W = wx < wr ? wx : wr;
            const int D = (int)zxc_log2_u32(W) - (LU - l);
            const uint32_t i0 = base + ((x - S) >> (LU - l));
            const uint64_t mass = pf[i0 + (W >> (LU - l))] - pf[i0];
            int t;
            if (D == 0) {
                t = l + 1; /* lone leaf: one memset + merges above */
            } else if (D == 1) {
                t = l; /* leaf pair: parent XOR-blend at depth l-1 */
            } else {
                t = (l - D) + 1; /* flat root: one unpack + merges above */
                if (D > ZXC_PIVCO_UNPACK_FLAT_SIMD_MAX) t += ZXC_HUF_NUDGE_DEEP_FLAT_PENALTY;
            }
            touches += mass * (uint64_t)t;
            x += W;
        }
        S = end;
        base += blc[l];
    }
    touches += (uint64_t)ZXC_HUF_NUDGE_LEVEL_COST * (uint64_t)(max_len + 1);
    out->bits = bits;
    out->touches = touches;
}
 
static ZXC_ALWAYS_INLINE uint64_t zxc_huf_nudge_j(const zxc_huf_nudge_cost_t* c) {
    return 256U * c->bits + (uint64_t)ZXC_HUF_NUDGE_LAMBDA_Q8 * c->touches;
}
 
static uint32_t zxc_huf_nudge_clamp(const uint32_t want, const uint32_t s, const uint32_t n_rem,
                                    const int l, const int cap) {
    if (n_rem <= s || l >= cap) return n_rem; /* forced finish */
    const uint32_t m = 1U << (cap - l);
    const uint32_t lo = (2 * s > n_rem) ? 2 * s - n_rem : 0;
    uint32_t hi = s - 1;
    const uint32_t cap_hi = (s * m - n_rem) / (m - 1);
    if (cap_hi < hi) hi = cap_hi;
    const uint32_t c = want < lo ? lo : want;
    return c > hi ? hi : c;
}
 
static void zxc_huf_nudge_complete(uint32_t* RESTRICT blc, const int l, uint32_t s, uint32_t n_rem,
                                   const uint32_t* RESTRICT blc0, const int cap) {
    for (int j = l + 1; j <= cap && n_rem; j++) {
        const uint32_t c = zxc_huf_nudge_clamp(blc0[j], s, n_rem, j, cap);
        blc[j] = c;
        n_rem -= c;
        s = 2 * (s - c);
    }
}
 
static void zxc_huf_nudge_walk(const uint32_t* RESTRICT blc0, const uint64_t* RESTRICT pf_rank,
                               const int n, const int cap, uint32_t* RESTRICT out_blc) {
    enum { LU = ZXC_HUF_MAX_CODE_LEN_ULTRA };
    ZXC_MEMSET(out_blc, 0, (LU + 1) * sizeof(uint32_t));
    uint32_t s = 2; /* the root's two slots at level 1 (n >= 4 => root is internal) */
    uint32_t n_rem = (uint32_t)n;
    for (int l = 1; l <= cap && n_rem; l++) {
        struct {
            uint32_t c;
            int flat_d; /* > 0: complete as one uniform run at depth l + flat_d */
        } cand[6];
        int n_cand = 0;
        const uint32_t c_base = zxc_huf_nudge_clamp(blc0[l], s, n_rem, l, cap);
        cand[n_cand].c = c_base;
        cand[n_cand++].flat_d = 0;
        if (c_base < n_rem) {
            // Shape the internal-node count i = s - c toward few set bits.
            const uint32_t i_base = s - c_base;
            const uint32_t i_dn = 1U << zxc_log2_u32(i_base);
            const uint32_t rest = i_base - i_dn;
            const uint32_t shapes[3] = {i_dn, i_dn << 1,
                                        i_dn | (rest ? 1U << zxc_log2_u32(rest) : 0U)};
            for (int k = 0; k < 3; k++) {
                const uint32_t want = s > shapes[k] ? s - shapes[k] : 0;
                const uint32_t c = zxc_huf_nudge_clamp(want, s, n_rem, l, cap);
                int dup = 0;
                for (int p = 0; p < n_cand; p++) dup |= (cand[p].c == c);
                if (!dup) {
                    cand[n_cand].c = c;
                    cand[n_cand++].flat_d = 0;
                }
            }
            // Uniform finish: c leaves here, the rest as i complete depth-D
            // subtrees (c*(2^D - 1) == s*2^D - n_rem must divide exactly).
            for (int d = 1; d <= cap - l && n_cand < 6; d++) {
                const uint32_t den = (1U << d) - 1;
                const int32_t num = (int32_t)(s << d) - (int32_t)n_rem;
                if (num < 0 || (uint32_t)num % den) continue;
                const uint32_t c = (uint32_t)num / den;
                if (c >= s || c >= n_rem) continue; /* need i >= 1 and a remainder */
                cand[n_cand].c = c;
                cand[n_cand++].flat_d = d;
                break; /* the shallowest exact finish is the aggressive one */
            }
        }
        uint64_t best_j = UINT64_MAX;
        uint32_t best_c = c_base;
        for (int k = 0; k < n_cand; k++) {
            uint32_t tmp[LU + 1];
            ZXC_MEMCPY(tmp, out_blc, sizeof(tmp));
            tmp[l] = cand[k].c;
            const uint32_t rem = n_rem - cand[k].c;
            if (rem) {
                if (cand[k].flat_d)
                    tmp[l + cand[k].flat_d] = rem;
                else
                    zxc_huf_nudge_complete(tmp, l, 2 * (s - cand[k].c), rem, blc0, cap);
            }
            zxc_huf_nudge_cost_t cc;
            zxc_huf_nudge_eval(tmp, pf_rank, &cc);
            const uint64_t j = zxc_huf_nudge_j(&cc);
            if (j < best_j) {
                best_j = j;
                best_c = cand[k].c;
            }
        }
        out_blc[l] = best_c;
        n_rem -= best_c;
        s = 2 * (s - best_c);
    }
}
 
static uint64_t zxc_huf_nudge_dp_run_j(const int lu, const int lc, const int g_log2,
                                       const uint32_t s, const uint32_t c,
                                       const uint64_t* RESTRICT pfg, const uint32_t k) {
    if (!c) return 0;
    const int lr = lc + g_log2; /* real code length of these leaves */
    const uint32_t w = 1U << (lu - lc);
    const uint32_t S = (1U << lu) - s * w;
    const uint32_t end = S + c * w;
    const uint64_t bits = (uint64_t)lr * (pfg[k + c] - pfg[k]);
    uint64_t touches = 0;
    uint32_t x = S;
    while (x < end) {
        const uint32_t wx = x ? (x & (0U - x)) : (1U << lu);
        const uint32_t wr = 1U << zxc_log2_u32(end - x);
        const uint32_t W = wx < wr ? wx : wr;
        const int d = (int)zxc_log2_u32(W >> (lu - lc)) + g_log2;
        const uint32_t i0 = k + ((x - S) >> (lu - lc));
        const uint64_t mass = pfg[i0 + (W >> (lu - lc))] - pfg[i0];
        int t;
        if (d == 0) {
            t = lr + 1;
        } else if (d == 1) {
            t = lr;
        } else {
            t = (lr - d) + 1;
            if (d > ZXC_PIVCO_UNPACK_FLAT_SIMD_MAX) t += ZXC_HUF_NUDGE_DEEP_FLAT_PENALTY;
        }
        touches += mass * (uint64_t)t;
        x += W;
    }
    return 256U * bits + (uint64_t)ZXC_HUF_NUDGE_LAMBDA_Q8 * touches;
}
 
static int zxc_huf_nudge_dp_solve(zxc_huf_nudge_ws_t* ws, const uint64_t* RESTRICT pfg, const int m,
                                  const int cap_c, const int lu, const int g_log2,
                                  uint32_t* RESTRICT out_cblc) {
    if (UNLIKELY(m < 2 || m > ZXC_HUF_NUDGE_DP_M || cap_c < 1 ||
                 cap_c > ZXC_HUF_MAX_CODE_LEN_ULTRA))
        return 0;
 
    const size_t plane = (size_t)(m + 1) * (size_t)(m + 1);
    const size_t arrive_cnt = (size_t)(cap_c + 1) * plane;
    uint64_t* jcur = ws->jcur;
    uint64_t* jnxt = ws->jnxt;
    uint8_t* const arrive = ws->arrive;
    int ok = 0;
 
#define ZXC_NUDGE_DP_IDX(k, s) ((size_t)(k) * (size_t)(m + 1) + (size_t)(s))
    for (size_t i = 0; i < plane; i++) jcur[i] = UINT64_MAX;
    jcur[ZXC_NUDGE_DP_IDX(0, 2)] = 0;
    uint64_t best_j = UINT64_MAX;
    int best_l = 0;
    int best_k = 0;
    int best_s = 0;
    for (int lc = 1; lc <= cap_c; lc++) {
        for (size_t i = 0; i < plane; i++) jnxt[i] = UINT64_MAX;
        for (int k = 0; k < m; k++) {
            const uint32_t n_rem = (uint32_t)(m - k);
            for (uint32_t s = 1; s <= n_rem; s++) {
                const size_t from = ZXC_NUDGE_DP_IDX(k, s);
                if (UNLIKELY(from >= plane)) continue;
                const uint64_t j0 = jcur[from];
                if (j0 == UINT64_MAX) continue;
                if (s == n_rem) {
                    // Forced finish: fill every slot, close the tree here.
                    const uint64_t j =
                        j0 + zxc_huf_nudge_dp_run_j(lu, lc, g_log2, s, s, pfg, (uint32_t)k) +
                        (uint64_t)ZXC_HUF_NUDGE_LAMBDA_Q8 * (uint64_t)ZXC_HUF_NUDGE_LEVEL_COST *
                            (uint64_t)(lc + g_log2 + 1);
                    if (j < best_j) {
                        best_j = j;
                        best_l = lc;
                        best_k = k;
                        best_s = (int)s;
                    }
                    continue;
                }
                if (lc == cap_c) continue; /* must finish at the cap */
                const uint32_t mm = 1U << (cap_c - lc);
                const uint32_t lo = (2 * s > n_rem) ? 2 * s - n_rem : 0;
                uint32_t hi = s - 1;
                const uint32_t cap_hi = (s * mm - n_rem) / (mm - 1);
                if (cap_hi < hi) hi = cap_hi;
                for (uint32_t c = lo; c <= hi; c++) {
                    const uint64_t j =
                        j0 + zxc_huf_nudge_dp_run_j(lu, lc, g_log2, s, c, pfg, (uint32_t)k);
                    const size_t to = ZXC_NUDGE_DP_IDX(k + (int)c, 2 * (s - c));
                    const size_t arr = (size_t)(lc + 1) * plane + to;
 
                    if (UNLIKELY(to >= plane || arr >= arrive_cnt)) continue;
                    if (j < jnxt[to]) {
                        jnxt[to] = j;
                        arrive[arr] = (uint8_t)c;
                    }
                }
            }
        }
        uint64_t* tmp = jcur;
        jcur = jnxt;
        jnxt = tmp;
    }
    if (best_j == UINT64_MAX) goto done;
    ZXC_MEMSET(out_cblc, 0, (ZXC_HUF_MAX_CODE_LEN_ULTRA + 1) * sizeof(uint32_t));
    out_cblc[best_l] = (uint32_t)best_s;
    {
        int k = best_k;
        int s = best_s;
        for (int lc = best_l; lc > 1; lc--) {
            const uint32_t c = arrive[(size_t)lc * plane + ZXC_NUDGE_DP_IDX(k, s)];
            out_cblc[lc - 1] = c;
            s = s / 2 + (int)c;
            k -= (int)c;
        }
        if (UNLIKELY(k != 0 || s != 2)) goto done; /* backtrack corruption: drop candidate */
    }
    ok = 1;
done:
#undef ZXC_NUDGE_DP_IDX
    return ok;
}
 
static void zxc_huf_nudge_assign(const uint32_t* RESTRICT blc, const int16_t* RESTRICT sym_order,
                                 uint8_t* RESTRICT cand_len) {
    ZXC_MEMSET(cand_len, 0, ZXC_HUF_NUM_SYMBOLS);
    int r = 0;
    for (int l = 1; l <= ZXC_HUF_MAX_CODE_LEN_ULTRA; l++)
        for (uint32_t k = 0; k < blc[l]; k++) cand_len[sym_order[r++]] = (uint8_t)l;
}
 
void zxc_huf_nudge_cost(const uint8_t* RESTRICT code_len, const uint32_t* RESTRICT freq,
                        uint64_t* RESTRICT bits, uint64_t* RESTRICT touches) {
    uint32_t blc[ZXC_HUF_MAX_CODE_LEN_ULTRA + 1];
    (void)zxc_huf_nudge_classes(code_len, blc);
    uint64_t pf[ZXC_HUF_NUM_SYMBOLS + 1];
    zxc_huf_nudge_pf_canonical(code_len, freq, blc, pf);
    zxc_huf_nudge_cost_t c;
    zxc_huf_nudge_eval(blc, pf, &c);
    *bits = c.bits;
    *touches = c.touches;
}
 
int zxc_huf_nudge_code_lengths(const uint32_t* RESTRICT freq, uint8_t* RESTRICT code_len,
                               void* RESTRICT scratch, const size_t scratch_cap,
                               const int max_code_len) {
    enum { LU = ZXC_HUF_MAX_CODE_LEN_ULTRA };
    uint32_t blc0[LU + 1];
    const int n = zxc_huf_nudge_classes(code_len, blc0);
    // Degenerate alphabets have no shape freedom (and n == 1 is format-pinned
    // to code length exactly 1): leave them untouched.
    if (n < 4) return 0;
 
    // Tables and DP planes live past the rebuilds' region of the scratch when it
    // has the room (opt_scratch does), else in one allocation: no 10 KiB frame.
    const int owned = !scratch || scratch_cap < ZXC_HUF_NUDGE_SCRATCH_SIZE;
    zxc_huf_nudge_ws_t* const ws =
        owned ? (zxc_huf_nudge_ws_t*)ZXC_MALLOC(sizeof(*ws))
              : (zxc_huf_nudge_ws_t*)((uint8_t*)scratch + ZXC_HUF_NUDGE_SCRATCH_OFF);
    if (UNLIKELY(!ws)) return 0;
 
    // Baseline cost, exact (canonical-order frequency weighting).
    uint64_t* const pf = ws->pf;
    zxc_huf_nudge_pf_canonical(code_len, freq, blc0, pf);
    zxc_huf_nudge_cost_t c0;
    zxc_huf_nudge_eval(blc0, pf, &c0);
 
    // Frequency-rank order shared by every candidate (rank 0 = most frequent;
    // pm_leaves_sort is ascending, ties on symbol, so read it backwards).
    pm_leaf_t* const leaves = ws->leaves;
    int k = 0;
    for (int s = 0; s < ZXC_HUF_NUM_SYMBOLS; s++) {
        if (!freq[s]) continue;
        leaves[k].w = freq[s];
        leaves[k].sym = (int16_t)s;
        k++;
    }
    pm_leaves_sort(leaves, n);
    int16_t* const sym_order = ws->sym_order;
    uint64_t* const pf_rank = ws->pf_rank;
    pf_rank[0] = 0;
    for (int r = 0; r < n; r++) {
        sym_order[r] = leaves[n - 1 - r].sym;
        pf_rank[r + 1] = pf_rank[r] + leaves[n - 1 - r].w;
    }
 
    // Candidates: the greedy ledger walk, package-merge rebuilt at reduced caps
    // (the pass loop runs max_depth + 1 times, so depth cuts attack the decoder's
    // biggest fixed cost), and the slot-ledger DP below. All Kraft-exact.
    uint8_t (*const cand)[ZXC_HUF_NUM_SYMBOLS] = ws->cand;
    int n_cand = 0;
    {
        uint32_t blc_w[LU + 1];
        zxc_huf_nudge_walk(blc0, pf_rank, n, max_code_len, blc_w);
        zxc_huf_nudge_assign(blc_w, sym_order, cand[n_cand]);
        n_cand++;
    }
    int max_len0 = 0;
    for (int l = LU; l >= 1; l--) {
        if (blc0[l]) {
            max_len0 = l;
            break;
        }
    }
    if (max_len0 >= 2) {
        for (int cut = 1; cut <= 2; cut++) {
            const int cap2 = max_len0 - cut;
            if (cap2 < 2 || ((uint32_t)1 << cap2) < (uint32_t)n) break;
            if (zxc_huf_build_code_lengths(freq, cand[n_cand], scratch, cap2) != ZXC_OK) break;
            n_cand++;
        }
    }
 
    // Slot-ledger DP: optimal class counts under the rank-weighted model, over
    // groups of 2^G frequency-adjacent symbols (G = 0 is exact; coarser tiers
    // keep the plane cache-resident, bounding the call to a few hundred us).
    // When G does not divide n the lightest group is padded with zero-freq
    // ghost leaves on unused byte values - wire-legal, empty runs, and the
    // evaluator prices the burnt code space like anything else.
    {
        int g_log2 = 0; /* coarsen until the groups fit the DP plane */
        while (((n + (1 << g_log2) - 1) >> g_log2) > ZXC_HUF_NUDGE_DP_M) g_log2++;
        const int g = 1 << g_log2;
        const int m = (n + g - 1) / g;
        const int cap_c = max_code_len - g_log2;
        if (m >= 2 && cap_c >= 1 && m <= (1 << cap_c)) {
            uint64_t* const pfg = ws->pfg;
            for (int j2 = 0; j2 <= m; j2++) {
                int r = j2 * g;
                if (r > n) r = n;
                pfg[j2] = pf_rank[r];
            }
            uint32_t cblc[LU + 1];
            if (zxc_huf_nudge_dp_solve(ws, pfg, m, cap_c, LU - g_log2, g_log2, cblc)) {
                uint8_t* const cl = cand[n_cand];
                ZXC_MEMSET(cl, 0, ZXC_HUF_NUM_SYMBOLS);
                int r = 0;
                int ghosts = 0;
                uint8_t ghost_len = 0;
                for (int lc = 1; lc <= cap_c; lc++) {
                    for (uint32_t q = 0; q < cblc[lc]; q++) {
                        for (int e = 0; e < g; e++, r++) {
                            if (r < n) {
                                cl[sym_order[r]] = (uint8_t)(lc + g_log2);
                            } else {
                                ghost_len = (uint8_t)(lc + g_log2);
                                ghosts++;
                            }
                        }
                    }
                }
                for (int s = 0; s < ZXC_HUF_NUM_SYMBOLS && ghosts; s++) {
                    if (freq[s] == 0 && cl[s] == 0) {
                        cl[s] = ghost_len;
                        ghosts--;
                    }
                }
                n_cand++;
            }
        }
    }
 
    // Adopt the cheapest candidate clearing both guard rails and strictly
    // beating the baseline J; otherwise leave code_len byte-for-byte alone, so a
    // rejection never depends on package-merge tie assignment.
    const uint64_t j0 = zxc_huf_nudge_j(&c0);
    uint64_t best_j = j0;
    int best = -1;
    for (int ci = 0; ci < n_cand; ci++) {
        // Every live symbol must keep a code; ghost-padded candidates carry
        // MORE nonzero lengths than n, which is fine.
        int valid = 1;
        for (int s = 0; s < ZXC_HUF_NUM_SYMBOLS; s++) {
            if (freq[s] != 0 && cand[ci][s] == 0) {
                valid = 0;
                break;
            }
        }
        if (!valid) continue;
        uint32_t blc[LU + 1];
        (void)zxc_huf_nudge_classes(cand[ci], blc);
        zxc_huf_nudge_pf_canonical(cand[ci], freq, blc, pf);
        zxc_huf_nudge_cost_t c1;
        zxc_huf_nudge_eval(blc, pf, &c1);
        if (c1.bits * 1000 > c0.bits * ZXC_HUF_NUDGE_BITS_PERMIL) continue;
        if (c1.touches * 256 > c0.touches * ZXC_HUF_NUDGE_MERGE_Q8) continue;
        const uint64_t j = zxc_huf_nudge_j(&c1);
        if (j < best_j) {
            best_j = j;
            best = ci;
        }
    }
    if (best >= 0) ZXC_MEMCPY(code_len, cand[best], ZXC_HUF_NUM_SYMBOLS);
    if (owned) ZXC_FREE(ws);
    return best >= 0;
}
#endif /* ZXC_VARIANT_PRIMARY */
 
// ===========================================================================
// 128-byte length header: 256 x 4-bit lengths, low nibble first.
// =========================================================================
 
void zxc_huf_pack_lengths(const uint8_t* RESTRICT code_len, uint8_t* RESTRICT out) {
    for (int i = 0; i < ZXC_HUF_NUM_SYMBOLS; i += 2) {
        const uint8_t lo = code_len[i] & 0x0F;
        const uint8_t hi = code_len[i + 1] & 0x0F;
        out[i >> 1] = (uint8_t)(lo | (hi << 4));
    }
}
 
int zxc_huf_unpack_lengths(const uint8_t* RESTRICT in, uint8_t* RESTRICT code_len) {
    int max_len = 0;
    int n_present = 0;
    for (int i = 0; i < ZXC_HUF_NUM_SYMBOLS; i += 2) {
        const uint8_t b = in[i >> 1];
        const uint8_t lo = b & 0x0F;
        const uint8_t hi = (uint8_t)(b >> 4);
        code_len[i] = lo;
        code_len[i + 1] = hi;
        if (lo > max_len) max_len = lo;
        if (hi > max_len) max_len = hi;
        if (lo) n_present++;
        if (hi) n_present++;
    }
    if (UNLIKELY(max_len > ZXC_HUF_MAX_CODE_LEN_ULTRA || n_present == 0))
        return ZXC_ERROR_CORRUPT_DATA;
    return ZXC_OK;
}
 
// ===========================================================================
// PivCo-Huffman section codec (enc 2/3)
// ===========================================================================
//
// Level-ordered layout and the merge-based decode are from PivCo-Huffman
// by Marcin Zukowski (https://github.com/MarcinZukowski/pivco-huffman).
// Implemented independently here.
//
// Layout: [128-byte packed code lengths (literal sections only)] then, for
// every EMITTING node of the canonical code tree in BFS order, that node's
// run: one branch bit per symbol routed through the node (0 = left child,
// 1 = right), in sequence order -- or, for flat subtree roots, D packed bits
// per symbol (bit j = branch at relative depth j). Runs are LSB-first within
// bytes and byte-padded; descendants of a flat root emit nothing.
//
// The decoder recovers each node's symbol count for free: the root handles
// n symbols, and popcounting a node's bits yields its right child's count.
// Reconstruction runs bottom-up, one level at a time: leaves are runs of a
// single symbol; an internal node MERGES its two children's sequences under
// the control of its bits. Merges are branch-free shuffles (16 outputs per
// step on AArch64 via a two-register TBL, 8 on ARMv7 via VTBL4, 64 on
// AVX-512-VBMI2 via vpexpandb), which is what makes this layout
// decode faster than the serial bit-chain of the classic 4-stream layout on
// any target with a 16-byte shuffle.
//
// Level buffers ping-pong between `scratch` (odd depths) and `dst` (even
// depths, final output at depth 0), so one n-sized scratch suffices. Both
// buffers need read slack past `n` (speculative 16-byte kernel loads):
// ZXC_PAD_SIZE on dst (already true for lit/token buffers), and
// ZXC_PIVCO_SCRATCH_PAD on scratch.
 
static ZXC_ALWAYS_INLINE int zxc_pivco_popcnt32(const uint32_t v) { return (int)ZXC_POPCOUNT32(v); }
 
static ZXC_ALWAYS_INLINE int zxc_pivco_popcnt64(const uint64_t v) { return (int)ZXC_POPCOUNT64(v); }
 
static ZXC_ALWAYS_INLINE int zxc_pivco_leaves(const zxc_pivco_tree_t* RESTRICT t, const int d) {
    return t->leaf_base[d + 1] - t->leaf_base[d];
}
 
static ZXC_ALWAYS_INLINE int zxc_pivco_flat_leaves(const zxc_pivco_tree_t* RESTRICT t, const int g,
                                                   const int d, const int D) {
    const int e = d + D;
    return t->leaf_base[e] + t->base[e] - t->node_base[e] + ((g - t->base[d]) << D);
}
 
static void zxc_pivco_bitrev_permute(uint8_t* RESTRICT s, const int D) {
    const unsigned m = 1U << D;
    unsigned c = 0; /* r bit-reversed */
    for (unsigned r = 0; r < m; r++) {
        if (r < c) {
            const uint8_t x = s[r];
            s[r] = s[c];
            s[c] = x;
        }
        unsigned bit = m >> 1; /* reversed increment */
        for (; c & bit; bit >>= 1) c ^= bit;
        c |= bit;
    }
}
 
static int zxc_pivco_tree_build(const uint8_t* RESTRICT code_len, zxc_pivco_tree_t* RESTRICT t) {
    uint16_t bl_count[ZXC_HUF_MAX_CODE_LEN_ULTRA + 1] = {0};
    int n_present = 0;
    int max_depth = 0;
    for (int s = 0; s < ZXC_HUF_NUM_SYMBOLS; s++) {
        const int l = code_len[s];
        if (l == 0) continue;
        if (UNLIKELY(l > ZXC_HUF_MAX_CODE_LEN_ULTRA)) return -1;
        bl_count[l]++;
        n_present++;
        if (l > max_depth) max_depth = l;
    }
    if (UNLIKELY(n_present == 0)) return -1;
    if (n_present >= 2) {
        uint32_t kraft = 0;
        for (int l = 1; l <= ZXC_HUF_MAX_CODE_LEN_ULTRA; l++)
            kraft += (uint32_t)bl_count[l] << (ZXC_HUF_MAX_CODE_LEN_ULTRA - l);
        if (UNLIKELY(kraft != (1U << ZXC_HUF_MAX_CODE_LEN_ULTRA))) return -1;
    } else if (UNLIKELY(bl_count[1] != 1)) {
        return -1; /* the lone symbol is pinned at length 1 */
    }
    t->max_depth = max_depth;
 
    // Two children per internal node above, one under the lone symbol's root.
    // The d-bit prefixes start at `first`.
    unsigned width = 1;
    unsigned first = 0;
    t->node_base[0] = 0;
    t->leaf_base[0] = 0;
    for (int d = 0; d <= max_depth; d++) {
        t->base[d] = (int16_t)(t->node_base[d] - (int)first);
        t->node_base[d + 1] = (uint16_t)(t->node_base[d] + width);
        t->leaf_base[d + 1] = (uint16_t)(t->leaf_base[d] + bl_count[d]);
        first = 2 * (first + bl_count[d]);
        width = (d == 0 && n_present == 1) ? 1U : 2 * (width - bl_count[d]);
    }
 
    uint16_t pos[ZXC_HUF_MAX_CODE_LEN_ULTRA + 1];
    for (int d = 0; d <= max_depth; d++) pos[d] = t->leaf_base[d];
    for (int s = 0; s < ZXC_HUF_NUM_SYMBOLS; s++)
        if (code_len[s]) t->syms[pos[code_len[s]]++] = (uint8_t)s;
 
    // FORMAT RULE: a flat root is a MAXIMAL complete subtree of depth D >= 2
    // (D = 1 is a leaf pair). Bottom-up, flat_d holds complete heights: 0 at
    // leaves and the slot, INCOMPLETE across depths.
    enum { INCOMPLETE = 0xFE };
    t->flat_d[t->node_base[max_depth + 1]] = 0;
    for (int d = max_depth; d >= 0; d--) {
        uint8_t* const h = t->flat_d + t->node_base[d];
        const uint8_t* const hc = t->flat_d + t->node_base[d + 1];
        const int L = zxc_pivco_leaves(t, d);
        const int N = t->node_base[d + 1] - t->node_base[d];
        for (int i = 0; i < L; i++) h[i] = 0;
        for (int i = L; i < N; i++) {
            const uint8_t a = hc[2 * (i - L)];
            const uint8_t b = hc[2 * (i - L) + 1];
            h[i] = (a == b && a != INCOMPLETE) ? (uint8_t)(a + 1) : (uint8_t)INCOMPLETE;
        }
    }
    // Top-down: the first height >= 2 on a path is a flat root, the internal
    // nodes below are covered. The wire packs a path LSB first, so leaf r has
    // packed code r bit-reversed: store the leaves in that order.
    for (int d = 0; d < max_depth; d++) {
        uint8_t* const fd = t->flat_d + t->node_base[d];
        const uint8_t* const up =
            d ? t->flat_d + t->node_base[d - 1] + zxc_pivco_leaves(t, d - 1) : NULL;
        const int N = t->node_base[d + 1] - t->node_base[d];
        for (int i = zxc_pivco_leaves(t, d); i < N; i++) {
            if (up && up[i >> 1]) {
                fd[i] = ZXC_PIVCO_COVERED;
            } else if (fd[i] < 2 || fd[i] == INCOMPLETE) {
                fd[i] = 0;
            } else {
                const int g = t->node_base[d] + i;
                zxc_pivco_bitrev_permute(t->syms + zxc_pivco_flat_leaves(t, g, d, fd[i]), fd[i]);
            }
        }
    }
    return 0;
}
 
static void zxc_pivco_counts(const zxc_pivco_tree_t* RESTRICT t, const uint32_t* RESTRICT freq,
                             uint32_t* RESTRICT count) {
    count[t->node_base[t->max_depth + 1]] = 0;
    for (int d = t->max_depth; d >= 0; d--) {
        const int L = zxc_pivco_leaves(t, d);
        const int N = t->node_base[d + 1] - t->node_base[d];
        uint32_t* const cnt = count + t->node_base[d];
        const uint32_t* const ccnt = count + t->node_base[d + 1];
        const uint8_t* const sym = t->syms + t->leaf_base[d];
        for (int i = 0; i < L; i++) cnt[i] = freq[sym[i]];
        for (int i = L; i < N; i++) cnt[i] = ccnt[2 * (i - L)] + ccnt[2 * (i - L) + 1];
    }
}
 
static void zxc_pivco_codes(const zxc_pivco_tree_t* RESTRICT t, uint32_t* RESTRICT codes) {
    for (int d = 1; d <= t->max_depth; d++)
        for (int i = 0; i < zxc_pivco_leaves(t, d); i++)
            codes[t->syms[t->leaf_base[d] + i]] = (uint32_t)(t->node_base[d] + i - t->base[d]);
}
 
static ZXC_ALWAYS_INLINE size_t zxc_pivco_run_bytes(const uint32_t count, const uint8_t flat_d) {
    return flat_d ? ((size_t)count * flat_d + 7) / 8 : ((size_t)count + 7) / 8;
}
 
static size_t zxc_pivco_layout(const zxc_pivco_tree_t* RESTRICT t, const uint32_t* RESTRICT count,
                               uint32_t* RESTRICT cur) {
    size_t total = 0;
    for (int d = 0; d < t->max_depth; d++) {
        for (int g = t->node_base[d] + zxc_pivco_leaves(t, d); g < t->node_base[d + 1]; g++) {
            if (t->flat_d[g] == ZXC_PIVCO_COVERED) continue;
            if (cur) cur[g] = (uint32_t)(total << 3);
            total += zxc_pivco_run_bytes(count[g], t->flat_d[g]);
        }
    }
    return total;
}
 
size_t zxc_huf_calc_size(const uint32_t* RESTRICT freq, const uint8_t* RESTRICT code_len,
                         const int with_header) {
    zxc_pivco_tree_t t;
    if (UNLIKELY(zxc_pivco_tree_build(code_len, &t) != 0)) return SIZE_MAX;
    const size_t body = zxc_huf_calc_size_dict(freq, code_len, &t);
    if (UNLIKELY(body == SIZE_MAX)) return SIZE_MAX;
    return body + (with_header ? (size_t)ZXC_HUF_TABLE_SIZE : 0);
}
 
size_t zxc_huf_calc_size_dict(const uint32_t* RESTRICT freq, const uint8_t* RESTRICT code_len,
                              const zxc_pivco_tree_t* RESTRICT tree) {
    // Encodability is part of the estimate: a histogram symbol the table has no
    // code for would be silently ignored below and undercount. Report the
    // candidate unencodable instead.
    for (int k = 0; k < ZXC_HUF_NUM_SYMBOLS; k++)
        if (UNLIKELY(freq[k] != 0 && code_len[k] == 0)) return SIZE_MAX;
    uint32_t count[ZXC_PIVCO_NODE_SLOTS];
    zxc_pivco_counts(tree, freq, count);
    return zxc_pivco_layout(tree, count, NULL);
}
 
static int zxc_pivco_encode_core(const uint8_t* RESTRICT literals, const size_t n_literals,
                                 const uint32_t* RESTRICT freq, const uint8_t* RESTRICT code_len,
                                 const zxc_pivco_tree_t* RESTRICT t, const uint32_t* RESTRICT codes,
                                 uint8_t* RESTRICT dst, const size_t dst_cap,
                                 const int with_header) {
    if (UNLIKELY(n_literals == 0)) return ZXC_ERROR_CORRUPT_DATA;
    // Every literal present in the histogram must have a code (a dict table
    // may lack codes for symbols unseen in training: the caller then falls
    // back to a per-block table, but a direct call must fail loudly).
    for (int k = 0; k < ZXC_HUF_NUM_SYMBOLS; k++)
        if (UNLIKELY(freq[k] != 0 && code_len[k] == 0)) return ZXC_ERROR_CORRUPT_DATA;
 
    uint32_t count[ZXC_PIVCO_NODE_SLOTS];
    zxc_pivco_counts(t, freq, count);
    uint32_t cur[ZXC_PIVCO_MAX_NODES]; /* per-node write cursor, in payload bits */
    const size_t payload = zxc_pivco_layout(t, count, cur);
    const size_t hdr = with_header ? (size_t)ZXC_HUF_TABLE_SIZE : 0;
    // +2: the packed-code emitter uses a 3-byte read-modify-write that may
    // touch up to 2 bytes past the payload end.
    if (UNLIKELY(hdr + payload + 2 > dst_cap)) return ZXC_ERROR_DST_TOO_SMALL;
 
    if (with_header) zxc_huf_pack_lengths(code_len, dst);
    uint8_t* const out = dst + hdr;
    ZXC_MEMSET(out, 0, payload + 2);
 
    // Emit MSB-first while descending. At a flat root, emit the code's last D
    // bits, its packed code, in one shot and stop. Batching the per-bit
    // read-modify-write buys nothing: accumulators would still touch memory
    // once per bit, and flat roots absorb the dense levels.
    for (size_t i = 0; i < n_literals; i++) {
        const uint8_t s = literals[i];
        const uint32_t c = codes[s];
        const int l = code_len[s];
        for (int d = 0; d < l; d++) {
            const int g = t->base[d] + (int)(c >> (l - d));
            const int fd = t->flat_d[g];
            if (fd) {
                const uint32_t p = cur[g];
                cur[g] += (uint32_t)fd;
                uint8_t* q = out + (p >> 3);
                const uint32_t sh = p & 7U;
                uint32_t w = (uint32_t)q[0] | ((uint32_t)q[1] << 8) | ((uint32_t)q[2] << 16);
                w |= (c & ((1U << fd) - 1U)) << sh; /* fd <= 8, sh <= 7: fits 24 bits */
                q[0] = (uint8_t)w;
                q[1] = (uint8_t)(w >> 8);
                q[2] = (uint8_t)(w >> 16);
                break;
            }
            const uint32_t bit = (c >> (l - 1 - d)) & 1U;
            const uint32_t p = cur[g]++;
            out[p >> 3] |= (uint8_t)(bit << (p & 7));
        }
    }
    return (int)(hdr + payload);
}
 
int zxc_huf_encode_section(const uint8_t* RESTRICT literals, const size_t n_literals,
                           const uint32_t* RESTRICT freq, const uint8_t* RESTRICT code_len,
                           uint8_t* RESTRICT dst, const size_t dst_cap) {
    zxc_pivco_tree_t t;
    uint32_t codes[ZXC_HUF_NUM_SYMBOLS];
    if (UNLIKELY(zxc_pivco_tree_build(code_len, &t) != 0)) return ZXC_ERROR_CORRUPT_DATA;
    zxc_pivco_codes(&t, codes);
    return zxc_pivco_encode_core(literals, n_literals, freq, code_len, &t, codes, dst, dst_cap, 1);
}
 
int zxc_huf_encode_section_dict(const uint8_t* RESTRICT literals, const size_t n_literals,
                                const uint32_t* RESTRICT freq, const uint8_t* RESTRICT code_len,
                                const zxc_pivco_tree_t* RESTRICT tree,
                                const uint32_t* RESTRICT codes, uint8_t* RESTRICT dst,
                                const size_t dst_cap) {
    return zxc_pivco_encode_core(literals, n_literals, freq, code_len, tree, codes, dst, dst_cap,
                                 0);
}
 
// ISA-independent cold dict setup
#if defined(ZXC_VARIANT_PRIMARY)
 
int zxc_huf_dict_tree_build(const uint8_t* RESTRICT packed_lengths, zxc_pivco_tree_t* RESTRICT tree,
                            uint32_t* RESTRICT codes, uint8_t* RESTRICT code_len) {
    const int rc = zxc_huf_unpack_lengths(packed_lengths, code_len);
    if (UNLIKELY(rc != ZXC_OK)) return rc;
    if (UNLIKELY(zxc_pivco_tree_build(code_len, tree) != 0)) return ZXC_ERROR_CORRUPT_DATA;
    zxc_pivco_codes(tree, codes);
    return ZXC_OK;
}
#endif /* ZXC_VARIANT_PRIMARY */
 
#define ZXC_PIVCO_MERGE_STEP(w, pc) \
    do {                            \
        rp += (size_t)(pc);         \
        lp += (size_t)((w) - (pc)); \
        i += (w);                   \
    } while (0)
 
static ZXC_ALWAYS_INLINE void zxc_pivco_merge(uint8_t* RESTRICT out, const uint8_t* RESTRICT src,
                                              const size_t nl, const size_t nr,
                                              const uint8_t* RESTRICT bits) {
    const uint8_t* L = src;
    const uint8_t* R = src + nl;
    const size_t n = nl + nr;
    size_t i = 0;
    size_t lp = 0;
    size_t rp = 0;
#if defined(ZXC_USE_NEON64)
    // 16 outputs per step: two-register TBL over {L[0..15], R[0..15]} with the
    // signed-index tables (see zxc_pivco_tables.h).
    while (i + 16 <= n) {
        const uint8_t b0 = bits[i >> 3];
        const uint8_t b1 = bits[(i >> 3) + 1];
        uint8x16x2_t tb;
        tb.val[0] = vld1q_u8(L + lp);
        tb.val[1] = vld1q_u8(R + rp);
        const int pc0 = zxc_pivco_popcnt32(b0);
        const uint8x8_t ia = vld1_u8(zxc_pivco_idxa_u8[b0]);
        const uint8x8_t ib = vld1_u8(zxc_pivco_idxb_pre[pc0][b1]);
        vst1q_u8(out + i, vqtbl2q_u8(tb, vcombine_u8(ia, ib)));
        const int pc = pc0 + zxc_pivco_popcnt32(b1);
        ZXC_PIVCO_MERGE_STEP(16, pc);
    }
#else /* x86 tiers */
#if defined(ZXC_USE_AVX512) && defined(__AVX512VBMI2__)
    // 64 outputs per step: the merge IS a pair of byte expands - L's bytes fill
    // the control word's 0-bit lanes, R's the 1-bit lanes. The masked expand
    // loads only touch the bytes they consume, so the speculative reads stay
    // fault-safe without extra buffer slack.
    const __m512i zero = _mm512_setzero_si512();
    while (i + 64 <= n) {
        uint64_t ctrl;
        ZXC_MEMCPY(&ctrl, bits + (i >> 3), 8);
        const __m512i vl =
            _mm512_mask_expandloadu_epi8(zero, (__mmask64)~ctrl, (const void*)(L + lp));
        const __m512i outv =
            _mm512_mask_expandloadu_epi8(vl, (__mmask64)ctrl, (const void*)(R + rp));
        _mm512_storeu_si512((void*)(out + i), outv);
        const int pc = zxc_pivco_popcnt64(ctrl);
        ZXC_PIVCO_MERGE_STEP(64, pc);
    }
#endif
#if defined(ZXC_USE_AVX512) || defined(ZXC_USE_AVX2) || \
    (defined(ZXC_USE_SSE2) && defined(__SSSE3__))
    // 16 outputs/step: pshufb(L, ix+0x70) | pshufb(R, ix-16) selects L[ix] or
    // R[ix-16] (out-of-range lanes zero out); 32/step unrolls it twice.
#if defined(ZXC_USE_AVX2)
    // Same thing at 256 bits, 32 outputs/step. Both selects take the index
    // unbiased: table lanes all sit in [0, 31], so pshufb masks to 4 bits and
    // bit 4 is the L/R selector - shifted into the sign it drives the blend,
    // and both bias vectors go away. L and R stay two 128-bit loads: the
    // second cursor is lp + (16 - pcA), never lp + 16.
    while (i + 32 <= n) {
        const uint8_t* cb = bits + (i >> 3);
        const int pcA0 = zxc_pivco_popcnt32(cb[0]);
        const int pcA = pcA0 + zxc_pivco_popcnt32(cb[1]);
        const int pcB0 = zxc_pivco_popcnt32(cb[2]);
        const int pcB = pcB0 + zxc_pivco_popcnt32(cb[3]);
        const size_t lpB = lp + (size_t)(16 - pcA);
        const size_t rpB = rp + (size_t)pcA;
        const __m256i vl = _mm256_inserti128_si256(
            _mm256_castsi128_si256(_mm_loadu_si128((const __m128i*)(const void*)(L + lp))),
            _mm_loadu_si128((const __m128i*)(const void*)(L + lpB)), 1);
        const __m256i vr = _mm256_inserti128_si256(
            _mm256_castsi128_si256(_mm_loadu_si128((const __m128i*)(const void*)(R + rp))),
            _mm_loadu_si128((const __m128i*)(const void*)(R + rpB)), 1);
        const __m128i ixA = _mm_unpacklo_epi64(
            _mm_loadl_epi64((const __m128i*)(const void*)zxc_pivco_idxa_u8[cb[0]]),
            _mm_loadl_epi64((const __m128i*)(const void*)zxc_pivco_idxb_pre[pcA0][cb[1]]));
        const __m128i ixB = _mm_unpacklo_epi64(
            _mm_loadl_epi64((const __m128i*)(const void*)zxc_pivco_idxa_u8[cb[2]]),
            _mm_loadl_epi64((const __m128i*)(const void*)zxc_pivco_idxb_pre[pcB0][cb[3]]));
        const __m256i ix = _mm256_inserti128_si256(_mm256_castsi128_si256(ixA), ixB, 1);
        _mm256_storeu_si256(
            (__m256i*)(void*)(out + i),
            _mm256_blendv_epi8(_mm256_shuffle_epi8(vl, ix), _mm256_shuffle_epi8(vr, ix),
                               _mm256_slli_epi16(ix, 3)));
        lp = lpB + (size_t)(16 - pcB);
        rp = rpB + (size_t)pcB;
        i += 32;
    }
#else
    while (i + 32 <= n) {
        const uint8_t* cb = bits + (i >> 3);
        const int pcA0 = zxc_pivco_popcnt32(cb[0]);
        const int pcA = pcA0 + zxc_pivco_popcnt32(cb[1]);
        const int pcB0 = zxc_pivco_popcnt32(cb[2]);
        const int pcB = pcB0 + zxc_pivco_popcnt32(cb[3]);
        const size_t lpB = lp + (size_t)(16 - pcA);
        const size_t rpB = rp + (size_t)pcA;
        const __m128i vlA = _mm_loadu_si128((const __m128i*)(const void*)(L + lp));
        const __m128i vrA = _mm_loadu_si128((const __m128i*)(const void*)(R + rp));
        const __m128i vlB = _mm_loadu_si128((const __m128i*)(const void*)(L + lpB));
        const __m128i vrB = _mm_loadu_si128((const __m128i*)(const void*)(R + rpB));
        const __m128i iaA = _mm_loadl_epi64((const __m128i*)(const void*)zxc_pivco_idxa_u8[cb[0]]);
        const __m128i ibA =
            _mm_loadl_epi64((const __m128i*)(const void*)zxc_pivco_idxb_pre[pcA0][cb[1]]);
        const __m128i iaB = _mm_loadl_epi64((const __m128i*)(const void*)zxc_pivco_idxa_u8[cb[2]]);
        const __m128i ibB =
            _mm_loadl_epi64((const __m128i*)(const void*)zxc_pivco_idxb_pre[pcB0][cb[3]]);
        const __m128i ixA = _mm_unpacklo_epi64(iaA, ibA);
        const __m128i ixB = _mm_unpacklo_epi64(iaB, ibB);
        const __m128i selA =
            _mm_or_si128(_mm_shuffle_epi8(vlA, _mm_add_epi8(ixA, _mm_set1_epi8(0x70))),
                         _mm_shuffle_epi8(vrA, _mm_sub_epi8(ixA, _mm_set1_epi8(16))));
        const __m128i selB =
            _mm_or_si128(_mm_shuffle_epi8(vlB, _mm_add_epi8(ixB, _mm_set1_epi8(0x70))),
                         _mm_shuffle_epi8(vrB, _mm_sub_epi8(ixB, _mm_set1_epi8(16))));
        _mm_storeu_si128((__m128i*)(void*)(out + i), selA);
        _mm_storeu_si128((__m128i*)(void*)(out + i + 16), selB);
        lp = lpB + (size_t)(16 - pcB);
        rp = rpB + (size_t)pcB;
        i += 32;
    }
#endif
    while (i + 16 <= n) {
        const uint8_t b0 = bits[i >> 3];
        const uint8_t b1 = bits[(i >> 3) + 1];
        const __m128i vl = _mm_loadu_si128((const __m128i*)(const void*)(L + lp));
        const __m128i vr = _mm_loadu_si128((const __m128i*)(const void*)(R + rp));
        const int pc0 = zxc_pivco_popcnt32(b0);
        const __m128i ia = _mm_loadl_epi64((const __m128i*)(const void*)zxc_pivco_idxa_u8[b0]);
        const __m128i ib =
            _mm_loadl_epi64((const __m128i*)(const void*)zxc_pivco_idxb_pre[pc0][b1]);
        const __m128i ix = _mm_unpacklo_epi64(ia, ib);
        const __m128i sell = _mm_shuffle_epi8(vl, _mm_add_epi8(ix, _mm_set1_epi8(0x70)));
        const __m128i selr = _mm_shuffle_epi8(vr, _mm_sub_epi8(ix, _mm_set1_epi8(16)));
        _mm_storeu_si128((__m128i*)(void*)(out + i), _mm_or_si128(sell, selr));
        const int pc = pc0 + zxc_pivco_popcnt32(b1);
        ZXC_PIVCO_MERGE_STEP(16, pc);
    }
#elif defined(ZXC_USE_NEON32)
    // 8 outputs per step: four-register VTBL over {L[0..7], -, R[0..7], -}. The
    // shared index tables address a 32-lane view, and an 8-output step only
    // touches lanes 0..7 and 16..23, so the two spare d-registers stay zero.
    while (i + 8 <= n) {
        const uint8_t b = bits[i >> 3];
        uint8x8x4_t tb;
        tb.val[0] = vld1_u8(L + lp);
        tb.val[1] = vdup_n_u8(0);
        tb.val[2] = vld1_u8(R + rp);
        tb.val[3] = vdup_n_u8(0);
        vst1_u8(out + i, vtbl4_u8(tb, vld1_u8(zxc_pivco_idxa_u8[b])));
        const int pc = zxc_pivco_popcnt32(b);
        ZXC_PIVCO_MERGE_STEP(8, pc);
    }
#endif
#endif /* x86 tiers */
    // Portable 8-output path: byte selects from a 24-byte scratch laid out
    // as the low half of the 32-lane view the shared index tables assume
    // (L at offset 0, R at offset 16).
    while (i + 8 <= n) {
        const uint8_t b = bits[i >> 3];
        uint8_t comb[24];
        ZXC_MEMCPY(comb, L + lp, 8);
        ZXC_MEMCPY(comb + 16, R + rp, 8);
        const uint8_t* ix = zxc_pivco_idxa_u8[b];
        for (int j = 0; j < 8; j++) out[i + (size_t)j] = comb[ix[j]];
        const int pc = zxc_pivco_popcnt32(b);
        ZXC_PIVCO_MERGE_STEP(8, pc);
    }
    while (i < n) {
        const int bit = (bits[i >> 3] >> (i & 7)) & 1;
        out[i++] = bit ? R[rp++] : L[lp++];
    }
}
 
#undef ZXC_PIVCO_MERGE_STEP
 
static void zxc_pivco_unpack_flat(uint8_t* RESTRICT out, const size_t n, const int D,
                                  const uint8_t* RESTRICT bits, const uint8_t* RESTRICT c2s) {
    size_t i = 0;
#if defined(ZXC_USE_NEON64)
    if (D == 2) {
        uint8_t c2s16[16];
        for (int k = 0; k < 16; k++) c2s16[k] = c2s[k & 3];
        const uint8x16_t vc2s = vld1q_u8(c2s16);
        static const uint8_t rep4[16] = {0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3};
        static const int8_t sh2[16] = {0, -2, -4, -6, 0, -2, -4, -6, 0, -2, -4, -6, 0, -2, -4, -6};
        const uint8x16_t vrep = vld1q_u8(rep4);
        const int8x16_t vsh = vld1q_s8(sh2);
        const uint8x16_t vmask = vdupq_n_u8(0x03);
        while (i + 16 <= n) {
            uint8x16_t raw = vdupq_n_u8(0);
            raw = vreinterpretq_u8_u32(
                vsetq_lane_u32(zxc_le32(bits + ((i * 2) >> 3)), vreinterpretq_u32_u8(raw), 0));
            const uint8x16_t rep = vqtbl1q_u8(raw, vrep);
            const uint8x16_t codes = vandq_u8(vshlq_u8(rep, vsh), vmask);
            vst1q_u8(out + i, vqtbl1q_u8(vc2s, codes));
            i += 16;
        }
    } else if (D == 3) {
        // Odd D straddle bytes: build a 16-bit window {byte_j,byte_j+1} per lane
        // (byte_j = Di>>3), right-shift by Di&7, mask D bits. 16/step, two u16x8
        // halves sharing the shift vector. D=3: 6 source bytes.
        uint8_t c2s16[16];
        for (int k = 0; k < 16; k++) c2s16[k] = c2s[k & 7];
        const uint8x16_t vc2s = vld1q_u8(c2s16); /* 8 entries, padded to 16 */
        static const uint8_t idxA[16] = {0, 1, 0, 1, 0, 1, 1, 2, 1, 2, 1, 2, 2, 3, 2, 3};
        static const uint8_t idxB[16] = {3, 4, 3, 4, 3, 4, 4, 5, 4, 5, 4, 5, 5, 6, 5, 6};
        static const int16_t sh8[8] = {0, -3, -6, -1, -4, -7, -2, -5};
        const uint8x16_t vidxA = vld1q_u8(idxA);
        const uint8x16_t vidxB = vld1q_u8(idxB);
        const int16x8_t vsh = vld1q_s16(sh8);
        const uint16x8_t vmask = vdupq_n_u16(0x07);
        while (i + 16 <= n) {
            const size_t off = (i * 3) >> 3;
            uint32_t w0;
            uint16_t w1;
            ZXC_MEMCPY(&w0, bits + off, 4);
            ZXC_MEMCPY(&w1, bits + off + 4, 2); /* exactly 6 bytes read */
            uint8x16_t raw = vdupq_n_u8(0);
            raw = vreinterpretq_u8_u32(vsetq_lane_u32(w0, vreinterpretq_u32_u8(raw), 0));
            raw = vreinterpretq_u8_u16(vsetq_lane_u16(w1, vreinterpretq_u16_u8(raw), 2));
            const uint16x8_t winA = vreinterpretq_u16_u8(vqtbl1q_u8(raw, vidxA));
            const uint16x8_t winB = vreinterpretq_u16_u8(vqtbl1q_u8(raw, vidxB));
            const uint8x8_t cA = vmovn_u16(vandq_u16(vshlq_u16(winA, vsh), vmask));
            const uint8x8_t cB = vmovn_u16(vandq_u16(vshlq_u16(winB, vsh), vmask));
            vst1q_u8(out + i, vqtbl1q_u8(vc2s, vcombine_u8(cA, cB)));
            i += 16;
        }
    } else if (D == 4) {
        const uint8x16_t vc2s = vld1q_u8(c2s); /* 16 entries */
        static const uint8_t rep2[16] = {0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7};
        static const int8_t sh4[16] = {0, -4, 0, -4, 0, -4, 0, -4, 0, -4, 0, -4, 0, -4, 0, -4};
        const uint8x16_t vrep = vld1q_u8(rep2);
        const int8x16_t vsh = vld1q_s8(sh4);
        const uint8x16_t vmask = vdupq_n_u8(0x0F);
        while (i + 16 <= n) {
            uint8x16_t raw = vdupq_n_u8(0);
            raw = vreinterpretq_u8_u64(
                vsetq_lane_u64(zxc_le64(bits + ((i * 4) >> 3)), vreinterpretq_u64_u8(raw), 0));
            const uint8x16_t rep = vqtbl1q_u8(raw, vrep);
            const uint8x16_t codes = vandq_u8(vshlq_u8(rep, vsh), vmask);
            vst1q_u8(out + i, vqtbl1q_u8(vc2s, codes));
            i += 16;
        }
    } else if (D == 5) {
        // Same window scheme as D=3; c2s has 32 entries (vqtbl2q). 10 bytes.
        const uint8x16x2_t vc2s = vld1q_u8_x2(c2s);
        static const uint8_t idxA[16] = {0, 1, 0, 1, 1, 2, 1, 2, 2, 3, 3, 4, 3, 4, 4, 5};
        static const uint8_t idxB[16] = {5, 6, 5, 6, 6, 7, 6, 7, 7, 8, 8, 9, 8, 9, 9, 10};
        static const int16_t sh8[8] = {0, -5, -2, -7, -4, -1, -6, -3};
        const uint8x16_t vidxA = vld1q_u8(idxA);
        const uint8x16_t vidxB = vld1q_u8(idxB);
        const int16x8_t vsh = vld1q_s16(sh8);
        const uint16x8_t vmask = vdupq_n_u16(0x1F);
        while (i + 16 <= n) {
            const size_t off = (i * 5) >> 3;
            uint64_t w0;
            uint16_t w1;
            ZXC_MEMCPY(&w0, bits + off, 8);
            ZXC_MEMCPY(&w1, bits + off + 8, 2);
            uint8x16_t raw = vdupq_n_u8(0);
            raw = vreinterpretq_u8_u64(vsetq_lane_u64(w0, vreinterpretq_u64_u8(raw), 0));
            raw = vreinterpretq_u8_u16(vsetq_lane_u16(w1, vreinterpretq_u16_u8(raw), 4));
            const uint16x8_t winA = vreinterpretq_u16_u8(vqtbl1q_u8(raw, vidxA));
            const uint16x8_t winB = vreinterpretq_u16_u8(vqtbl1q_u8(raw, vidxB));
            const uint8x8_t cA = vmovn_u16(vandq_u16(vshlq_u16(winA, vsh), vmask));
            const uint8x8_t cB = vmovn_u16(vandq_u16(vshlq_u16(winB, vsh), vmask));
            vst1q_u8(out + i, vqtbl2q_u8(vc2s, vcombine_u8(cA, cB)));
            i += 16;
        }
    } else if (D == 6) {
        // Same window scheme as D=3; c2s has 64 entries (vqtbl4q). 12 bytes.
        const uint8x16x4_t vc2s = vld1q_u8_x4(c2s);
        static const uint8_t idxA[16] = {0, 1, 0, 1, 1, 2, 2, 3, 3, 4, 3, 4, 4, 5, 5, 6};
        static const uint8_t idxB[16] = {6, 7, 6, 7, 7, 8, 8, 9, 9, 10, 9, 10, 10, 11, 11, 12};
        static const int16_t sh8[8] = {0, -6, -4, -2, 0, -6, -4, -2};
        const uint8x16_t vidxA = vld1q_u8(idxA);
        const uint8x16_t vidxB = vld1q_u8(idxB);
        const int16x8_t vsh = vld1q_s16(sh8);
        const uint16x8_t vmask = vdupq_n_u16(0x3F);
        while (i + 16 <= n) {
            const size_t off = (i * 6) >> 3;
            uint64_t w0;
            uint32_t w1;
            ZXC_MEMCPY(&w0, bits + off, 8);
            ZXC_MEMCPY(&w1, bits + off + 8, 4);
            uint8x16_t raw = vdupq_n_u8(0);
            raw = vreinterpretq_u8_u64(vsetq_lane_u64(w0, vreinterpretq_u64_u8(raw), 0));
            raw = vreinterpretq_u8_u32(vsetq_lane_u32(w1, vreinterpretq_u32_u8(raw), 2));
            const uint16x8_t winA = vreinterpretq_u16_u8(vqtbl1q_u8(raw, vidxA));
            const uint16x8_t winB = vreinterpretq_u16_u8(vqtbl1q_u8(raw, vidxB));
            const uint8x8_t cA = vmovn_u16(vandq_u16(vshlq_u16(winA, vsh), vmask));
            const uint8x8_t cB = vmovn_u16(vandq_u16(vshlq_u16(winB, vsh), vmask));
            vst1q_u8(out + i, vqtbl4q_u8(vc2s, vcombine_u8(cA, cB)));
            i += 16;
        }
    }
#elif defined(ZXC_USE_NEON32)
    // d-register mirrors of the AArch64 kernels: 8 codes per step, c2s lookup
    // via VTBL (16-entry table = two d-registers).
    if (D == 2) {
        uint8_t c2s8[8];
        for (int k = 0; k < 8; k++) c2s8[k] = c2s[k & 3];
        const uint8x8_t vc2s = vld1_u8(c2s8);
        static const uint8_t rep4[8] = {0, 0, 0, 0, 1, 1, 1, 1};
        static const int8_t sh2[8] = {0, -2, -4, -6, 0, -2, -4, -6};
        const uint8x8_t vrep = vld1_u8(rep4);
        const int8x8_t vsh = vld1_s8(sh2);
        const uint8x8_t vmask = vdup_n_u8(0x03);
        while (i + 8 <= n) {
            const uint16_t w = (uint16_t)((uint16_t)bits[(i * 2) >> 3] |
                                          ((uint16_t)bits[((i * 2) >> 3) + 1] << 8));
            const uint8x8_t raw = vreinterpret_u8_u16(vdup_n_u16(w));
            const uint8x8_t rep = vtbl1_u8(raw, vrep);
            const uint8x8_t codes = vand_u8(vshl_u8(rep, vsh), vmask);
            vst1_u8(out + i, vtbl1_u8(vc2s, codes));
            i += 8;
        }
    } else if (D == 3) {
        // 8 codes/step. 3-bit codes straddle bytes, so build a 16-bit window
        // {byte_j, byte_j+1} per lane via two VTBLs, then extract bits [s,s+3)
        // by multiply-by-2^(13-s) then >>13 - ARMv7 has no per-lane u16 variable
        // shift. 8 codes fit in 3 source bytes.
        const uint8x8_t vc2s = vld1_u8(c2s); /* 8 entries (2^3) */
        static const uint8_t idxlo[8] = {0, 1, 0, 1, 0, 1, 1, 2};
        static const uint8_t idxhi[8] = {1, 2, 1, 2, 2, 3, 2, 3};
        static const uint16_t mul8[8] = {8192, 1024, 128, 4096, 512, 64, 2048, 256};
        const uint8x8_t vidxlo = vld1_u8(idxlo);
        const uint8x8_t vidxhi = vld1_u8(idxhi);
        const uint16x8_t vmul = vld1q_u16(mul8);
        const uint16x8_t vmask = vdupq_n_u16(0x07);
        while (i + 8 <= n) {
            const size_t off = (i * 3) >> 3;
            const uint32_t w = (uint32_t)bits[off] | ((uint32_t)bits[off + 1] << 8) |
                               ((uint32_t)bits[off + 2] << 16); /* exactly 3 bytes */
            const uint8x8_t raw = vreinterpret_u8_u32(vdup_n_u32(w));
            const uint8x8_t wlo = vtbl1_u8(raw, vidxlo);
            const uint8x8_t whi = vtbl1_u8(raw, vidxhi);
            const uint16x8_t win = vreinterpretq_u16_u8(vcombine_u8(wlo, whi));
            const uint16x8_t codes16 = vandq_u16(vshrq_n_u16(vmulq_u16(win, vmul), 13), vmask);
            vst1_u8(out + i, vtbl1_u8(vc2s, vmovn_u16(codes16)));
            i += 8;
        }
    } else if (D == 4) {
        /* Register by register: GCC has no AArch32 vld1_u8_x2/_x4 before 14.
         * Loaded once per call, so the tuple form bought nothing anyway. */
        uint8x8x2_t vc2s;
        vc2s.val[0] = vld1_u8(c2s);
        vc2s.val[1] = vld1_u8(c2s + 8);
        static const uint8_t rep2[8] = {0, 0, 1, 1, 2, 2, 3, 3};
        static const int8_t sh4[8] = {0, -4, 0, -4, 0, -4, 0, -4};
        const uint8x8_t vrep = vld1_u8(rep2);
        const int8x8_t vsh = vld1_s8(sh4);
        const uint8x8_t vmask = vdup_n_u8(0x0F);
        while (i + 8 <= n) {
            uint32_t w;
            ZXC_MEMCPY(&w, bits + ((i * 4) >> 3), 4);
            const uint8x8_t raw = vreinterpret_u8_u32(vdup_n_u32(w));
            const uint8x8_t rep = vtbl1_u8(raw, vrep);
            const uint8x8_t codes = vand_u8(vshl_u8(rep, vsh), vmask);
            vst1_u8(out + i, vtbl2_u8(vc2s, codes));
            i += 8;
        }
    } else if (D == 5) {
        // Same window scheme as D=3; c2s has 32 entries (VTBL4). 5 source bytes.
        // Register by register: see D == 4.
        uint8x8x4_t vc2s;
        vc2s.val[0] = vld1_u8(c2s);
        vc2s.val[1] = vld1_u8(c2s + 8);
        vc2s.val[2] = vld1_u8(c2s + 16);
        vc2s.val[3] = vld1_u8(c2s + 24);
        static const uint8_t idxlo[8] = {0, 1, 0, 1, 1, 2, 1, 2};
        static const uint8_t idxhi[8] = {2, 3, 3, 4, 3, 4, 4, 5};
        static const uint16_t mul8[8] = {2048, 64, 512, 16, 128, 1024, 32, 256};
        const uint8x8_t vidxlo = vld1_u8(idxlo);
        const uint8x8_t vidxhi = vld1_u8(idxhi);
        const uint16x8_t vmul = vld1q_u16(mul8);
        const uint16x8_t vmask = vdupq_n_u16(0x1F);
        while (i + 8 <= n) {
            uint64_t w = 0;
            ZXC_MEMCPY(&w, bits + ((i * 5) >> 3), 5);
            const uint8x8_t raw = vcreate_u8(w);
            const uint16x8_t win =
                vreinterpretq_u16_u8(vcombine_u8(vtbl1_u8(raw, vidxlo), vtbl1_u8(raw, vidxhi)));
            const uint8x8_t codes =
                vmovn_u16(vandq_u16(vshrq_n_u16(vmulq_u16(win, vmul), 11), vmask));
            vst1_u8(out + i, vtbl4_u8(vc2s, codes));
            i += 8;
        }
    } else if (D == 6) {
        // Same window scheme; c2s has 64 entries: two VTBL4 halves, select by
        // code bit 5. 6 source bytes.
        uint8x8x4_t lo, hi;  // register by register: see D == 4
        lo.val[0] = vld1_u8(c2s);
        lo.val[1] = vld1_u8(c2s + 8);
        lo.val[2] = vld1_u8(c2s + 16);
        lo.val[3] = vld1_u8(c2s + 24);
        hi.val[0] = vld1_u8(c2s + 32);
        hi.val[1] = vld1_u8(c2s + 40);
        hi.val[2] = vld1_u8(c2s + 48);
        hi.val[3] = vld1_u8(c2s + 56);
        static const uint8_t idxlo[8] = {0, 1, 0, 1, 1, 2, 2, 3};
        static const uint8_t idxhi[8] = {3, 4, 3, 4, 4, 5, 5, 6};
        static const uint16_t mul8[8] = {1024, 16, 64, 256, 1024, 16, 64, 256};
        const uint8x8_t vidxlo = vld1_u8(idxlo);
        const uint8x8_t vidxhi = vld1_u8(idxhi);
        const uint16x8_t vmul = vld1q_u16(mul8);
        const uint16x8_t vmask = vdupq_n_u16(0x3F);
        while (i + 8 <= n) {
            uint64_t w = 0;
            ZXC_MEMCPY(&w, bits + ((i * 6) >> 3), 6);
            const uint8x8_t raw = vcreate_u8(w);
            const uint16x8_t win =
                vreinterpretq_u16_u8(vcombine_u8(vtbl1_u8(raw, vidxlo), vtbl1_u8(raw, vidxhi)));
            const uint8x8_t codes =
                vmovn_u16(vandq_u16(vshrq_n_u16(vmulq_u16(win, vmul), 10), vmask));
            const uint8x8_t sel =
                vbsl_u8(vcge_u8(codes, vdup_n_u8(32)), vtbl4_u8(hi, vsub_u8(codes, vdup_n_u8(32))),
                        vtbl4_u8(lo, codes));
            vst1_u8(out + i, sel);
            i += 8;
        }
    }
#elif defined(ZXC_USE_AVX512) || defined(ZXC_USE_AVX2) || \
    (defined(ZXC_USE_SSE2) && defined(__SSE4_1__))
#if defined(ZXC_USE_AVX512) && defined(__AVX512VBMI__)
    /* vpmultishiftqb reads an 8-bit field at an arbitrary bit offset: the
     * unaligned D-bit read, in one instruction. Two facts carry the rest -
     * 64 codes are exactly 8*D bytes and byte-aligned, so the tables are the
     * same every step; and only the low D bits of a field survive, so bytes
     * past 8*D are irrelevant and both accesses can be clamped to the run.
     * Clamping the store also ends the run, hence nothing narrower below. */
    if (D >= 2 && D <= 6) {
        uint8_t arr[64];
        uint8_t shf[64];
        for (int k = 0; k < 8; k++) {
            for (int j = 0; j < 8; j++) {
                arr[k * 8 + j] = (uint8_t)(k * D + j);
                shf[k * 8 + j] = (uint8_t)(j * D);
            }
        }
        const __m512i varr = _mm512_loadu_si512((const void*)arr);
        const __m512i vshf = _mm512_loadu_si512((const void*)shf);
        const __m512i vand = _mm512_set1_epi8((char)((1 << D) - 1));
        /* c2s holds 2^D entries and nothing is guaranteed past them. */
        const __mmask64 tm = (D == 6) ? ~(__mmask64)0 : (((__mmask64)1 << (1 << D)) - 1);
        const __m512i vc2s = _mm512_maskz_loadu_epi8(tm, (const void*)c2s);
        const size_t nbytes = zxc_pivco_run_bytes((uint32_t)n, (uint8_t)D);
        while (i < n) {
            const size_t off = ((size_t)i * (size_t)D) >> 3;
            const size_t avail = nbytes - off;
            const size_t left = n - i;
            const __mmask64 lm = (avail >= 64) ? ~(__mmask64)0 : (((__mmask64)1 << avail) - 1);
            const __mmask64 sm = (left >= 64) ? ~(__mmask64)0 : (((__mmask64)1 << left) - 1);
            const __m512i raw = _mm512_maskz_loadu_epi8(lm, (const void*)(bits + off));
            const __m512i q = _mm512_permutexvar_epi8(varr, raw);
            const __m512i codes = _mm512_and_si512(_mm512_multishift_epi64_epi8(vshf, q), vand);
            _mm512_mask_storeu_epi8((void*)(out + i), sm, _mm512_permutexvar_epi8(codes, vc2s));
            i += 64;
        }
        if (i > n) i = n;
    }
#else
#if defined(ZXC_USE_AVX2)
    /* Codes 16..31 start exactly 2*D bytes in, so the upper half runs the same
     * pattern over a shifted source: one shuffle vector serves both, and packus
     * leaves each half's codes in order in its own lane - no cross-lane fixup.
     * The bound keeps both 16-byte reads inside the run; the exact-length
     * 16-wide loops below take the tail. */
    if (D >= 2 && D <= 6) {
        uint8_t sa[16];
        uint8_t sb[16];
        uint16_t mu[8];
        for (int q = 0; q < 8; q++) {
            const int ba = (q * D) >> 3;
            sa[2 * q] = (uint8_t)ba;
            sa[2 * q + 1] = (uint8_t)(ba + 1);
            sb[2 * q] = (uint8_t)(ba + D); /* code q+8 sits D bytes further */
            sb[2 * q + 1] = (uint8_t)(ba + D + 1);
            mu[q] = (uint16_t)(1U << (16 - D - ((q * D) & 7)));
        }
        const __m256i vsa = _mm256_broadcastsi128_si256(_mm_loadu_si128((const __m128i*)(void*)sa));
        const __m256i vsb = _mm256_broadcastsi128_si256(_mm_loadu_si128((const __m128i*)(void*)sb));
        const __m256i vmu = _mm256_broadcastsi128_si256(_mm_loadu_si128((const __m128i*)(void*)mu));
        const __m256i vand = _mm256_set1_epi8((char)((1 << D) - 1));
        uint8_t tab[64];
        for (int e = 0; e < 64; e++) tab[e] = c2s[e & ((1 << D) - 1)];
        const __m256i t0 = _mm256_broadcastsi128_si256(_mm_loadu_si128((const __m128i*)(void*)tab));
        const __m256i t1 =
            _mm256_broadcastsi128_si256(_mm_loadu_si128((const __m128i*)(void*)(tab + 16)));
        const __m256i t2 =
            _mm256_broadcastsi128_si256(_mm_loadu_si128((const __m128i*)(void*)(tab + 32)));
        const __m256i t3 =
            _mm256_broadcastsi128_si256(_mm_loadu_si128((const __m128i*)(void*)(tab + 48)));
        const size_t nbytes = zxc_pivco_run_bytes((uint32_t)n, (uint8_t)D);
        while (i + 32 <= n && ((size_t)i * (size_t)D >> 3) + (size_t)(2 * D) + 16 <= nbytes) {
            const uint8_t* src = bits + ((size_t)i * (size_t)D >> 3);
            const __m256i raw = _mm256_inserti128_si256(
                _mm256_castsi128_si256(_mm_loadu_si128((const __m128i*)(const void*)src)),
                _mm_loadu_si128((const __m128i*)(const void*)(src + 2 * D)), 1);
            const __m256i ca =
                _mm256_srli_epi16(_mm256_mullo_epi16(_mm256_shuffle_epi8(raw, vsa), vmu), 16 - D);
            const __m256i cb =
                _mm256_srli_epi16(_mm256_mullo_epi16(_mm256_shuffle_epi8(raw, vsb), vmu), 16 - D);
            const __m256i codes = _mm256_and_si256(_mm256_packus_epi16(ca, cb), vand);
            __m256i sym;
            if (D <= 4) {
                sym = _mm256_shuffle_epi8(t0, codes);
            } else {
                const __m256i b4 = _mm256_slli_epi16(codes, 3); /* code bit 4 -> sign */
                const __m256i lo = _mm256_blendv_epi8(_mm256_shuffle_epi8(t0, codes),
                                                      _mm256_shuffle_epi8(t1, codes), b4);
                if (D == 5) {
                    sym = lo;
                } else {
                    const __m256i hi = _mm256_blendv_epi8(_mm256_shuffle_epi8(t2, codes),
                                                          _mm256_shuffle_epi8(t3, codes), b4);
                    sym = _mm256_blendv_epi8(lo, hi, _mm256_slli_epi16(codes, 2));
                }
            }
            _mm256_storeu_si256((__m256i*)(void*)(out + i), sym);
            i += 32;
        }
    }
#endif
    if (D == 4) {
        const __m128i vc2s = _mm_loadu_si128((const __m128i*)(const void*)c2s);
        const __m128i vrep = _mm_setr_epi8(0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7);
        const __m128i vmask = _mm_set1_epi8(0x0F);
        // even lanes keep low nibble, odd lanes take the high nibble: shift the
        // replicated bytes right by 4 via a 16-bit shift + odd-lane select.
        const __m128i vodd = _mm_set1_epi16(0x0F00 - 0x0F00 + (short)0xFF00);
        while (i + 16 <= n) {
            const __m128i raw =
                _mm_loadl_epi64((const __m128i*)(const void*)(bits + ((i * 4) >> 3)));
            const __m128i rep = _mm_shuffle_epi8(raw, vrep);
            const __m128i hi = _mm_and_si128(_mm_srli_epi16(rep, 4), _mm_set1_epi8(0x0F));
            const __m128i lo = _mm_and_si128(rep, vmask);
            const __m128i codes = _mm_blendv_epi8(lo, hi, vodd);
            _mm_storeu_si128((__m128i*)(void*)(out + i), _mm_shuffle_epi8(vc2s, codes));
            i += 16;
        }
    } else if (D == 2) {
        // 16 codes from 4 bytes. SSE has no per-lane byte shift, so shift by
        // multiply: (b << (6 - 2j)) >> 6 keeps bits [2j+1:2j]. Bytes spread to
        // u16 lanes, times {64,16,4,1}, shifted, then packed back.
        uint8_t c2s16[16];
        for (int k = 0; k < 16; k++) c2s16[k] = c2s[k & 3];
        const __m128i vc2s = _mm_loadu_si128((const __m128i*)(const void*)c2s16);
        const __m128i vspreadA =
            _mm_setr_epi8(0, (char)0x80, 0, (char)0x80, 0, (char)0x80, 0, (char)0x80, 1, (char)0x80,
                          1, (char)0x80, 1, (char)0x80, 1, (char)0x80);
        const __m128i vspreadB =
            _mm_setr_epi8(2, (char)0x80, 2, (char)0x80, 2, (char)0x80, 2, (char)0x80, 3, (char)0x80,
                          3, (char)0x80, 3, (char)0x80, 3, (char)0x80);
        const __m128i vmul = _mm_setr_epi16(64, 16, 4, 1, 64, 16, 4, 1);
        while (i + 16 <= n) {
            const __m128i raw = _mm_cvtsi32_si128(
                (int)((uint32_t)bits[(i * 2) >> 3] | ((uint32_t)bits[((i * 2) >> 3) + 1] << 8) |
                      ((uint32_t)bits[((i * 2) >> 3) + 2] << 16) |
                      ((uint32_t)bits[((i * 2) >> 3) + 3] << 24)));
            const __m128i a16 =
                _mm_srli_epi16(_mm_mullo_epi16(_mm_shuffle_epi8(raw, vspreadA), vmul), 6);
            const __m128i b16 =
                _mm_srli_epi16(_mm_mullo_epi16(_mm_shuffle_epi8(raw, vspreadB), vmul), 6);
            const __m128i codes = _mm_and_si128(_mm_packus_epi16(a16, b16), _mm_set1_epi8(3));
            _mm_storeu_si128((__m128i*)(void*)(out + i), _mm_shuffle_epi8(vc2s, codes));
            i += 16;
        }
    } else if (D == 3) {
        // 16-bit window {byte_j,byte_j+1} per lane (byte_j = Di>>3); SSE has no
        // per-lane u16 shift, so extract bits [s,s+D) of window w by w*2^(16-D-s)
        // (field lands in the top D bits) then >>(16-D); s_j = Dj&7. D=3: 6 bytes.
        uint8_t c2s16[16];
        for (int k = 0; k < 16; k++) c2s16[k] = c2s[k & 7];
        const __m128i vc2s = _mm_loadu_si128((const __m128i*)(const void*)c2s16);
        const __m128i shufA = _mm_setr_epi8(0, 1, 0, 1, 0, 1, 1, 2, 1, 2, 1, 2, 2, 3, 2, 3);
        const __m128i shufB = _mm_setr_epi8(3, 4, 3, 4, 3, 4, 4, 5, 4, 5, 4, 5, 5, 6, 5, 6);
        const __m128i vmul = _mm_setr_epi16(8192, 1024, 128, 4096, 512, 64, 2048, 256);
        while (i + 16 <= n) {
            const size_t off = (i * 3) >> 3;
            uint32_t w0;
            uint16_t w1;
            ZXC_MEMCPY(&w0, bits + off, 4);
            ZXC_MEMCPY(&w1, bits + off + 4, 2); /* exactly 6 bytes read */
            const __m128i raw = _mm_insert_epi16(_mm_cvtsi32_si128((int)w0), w1, 2);
            const __m128i cA =
                _mm_srli_epi16(_mm_mullo_epi16(_mm_shuffle_epi8(raw, shufA), vmul), 13);
            const __m128i cB =
                _mm_srli_epi16(_mm_mullo_epi16(_mm_shuffle_epi8(raw, shufB), vmul), 13);
            const __m128i codes = _mm_and_si128(_mm_packus_epi16(cA, cB), _mm_set1_epi8(7));
            _mm_storeu_si128((__m128i*)(void*)(out + i), _mm_shuffle_epi8(vc2s, codes));
            i += 16;
        }
    } else if (D == 5) {
        // c2s has 32 entries: pshufb does 16, so select lo/hi by code bit 4.
        const __m128i vlo = _mm_loadu_si128((const __m128i*)(const void*)c2s);
        const __m128i vhi = _mm_loadu_si128((const __m128i*)(const void*)(c2s + 16));
        const __m128i shufA = _mm_setr_epi8(0, 1, 0, 1, 1, 2, 1, 2, 2, 3, 3, 4, 3, 4, 4, 5);
        const __m128i shufB = _mm_setr_epi8(5, 6, 5, 6, 6, 7, 6, 7, 7, 8, 8, 9, 8, 9, 9, 10);
        const __m128i vmul = _mm_setr_epi16(2048, 64, 512, 16, 128, 1024, 32, 256);
        while (i + 16 <= n) {
            const size_t off = (i * 5) >> 3;
            uint16_t w1;
            ZXC_MEMCPY(&w1, bits + off + 8, 2); /* 10 source bytes */
            const __m128i raw =
                _mm_insert_epi16(_mm_loadl_epi64((const __m128i*)(const void*)(bits + off)), w1, 4);
            const __m128i cA =
                _mm_srli_epi16(_mm_mullo_epi16(_mm_shuffle_epi8(raw, shufA), vmul), 11);
            const __m128i cB =
                _mm_srli_epi16(_mm_mullo_epi16(_mm_shuffle_epi8(raw, shufB), vmul), 11);
            const __m128i codes = _mm_and_si128(_mm_packus_epi16(cA, cB), _mm_set1_epi8(0x1F));
            const __m128i sel =
                _mm_blendv_epi8(_mm_shuffle_epi8(vlo, codes), _mm_shuffle_epi8(vhi, codes),
                                _mm_slli_epi16(codes, 3));
            _mm_storeu_si128((__m128i*)(void*)(out + i), sel);
            i += 16;
        }
    } else if (D == 6) {
        // c2s has 64 entries: 4 pshufb sub-tables, select by code bits 4-5.
        const __m128i t0 = _mm_loadu_si128((const __m128i*)(const void*)c2s);
        const __m128i t1 = _mm_loadu_si128((const __m128i*)(const void*)(c2s + 16));
        const __m128i t2 = _mm_loadu_si128((const __m128i*)(const void*)(c2s + 32));
        const __m128i t3 = _mm_loadu_si128((const __m128i*)(const void*)(c2s + 48));
        const __m128i shufA = _mm_setr_epi8(0, 1, 0, 1, 1, 2, 2, 3, 3, 4, 3, 4, 4, 5, 5, 6);
        const __m128i shufB = _mm_setr_epi8(6, 7, 6, 7, 7, 8, 8, 9, 9, 10, 9, 10, 10, 11, 11, 12);
        const __m128i vmul = _mm_setr_epi16(1024, 16, 64, 256, 1024, 16, 64, 256);
        while (i + 16 <= n) {
            const size_t off = (i * 6) >> 3;
            uint32_t w1;
            ZXC_MEMCPY(&w1, bits + off + 8, 4); /* 12 source bytes */
            const __m128i raw = _mm_insert_epi32(
                _mm_loadl_epi64((const __m128i*)(const void*)(bits + off)), (int)w1, 2);
            const __m128i cA =
                _mm_srli_epi16(_mm_mullo_epi16(_mm_shuffle_epi8(raw, shufA), vmul), 10);
            const __m128i cB =
                _mm_srli_epi16(_mm_mullo_epi16(_mm_shuffle_epi8(raw, shufB), vmul), 10);
            const __m128i codes = _mm_and_si128(_mm_packus_epi16(cA, cB), _mm_set1_epi8(0x3F));
            const __m128i m4 = _mm_slli_epi16(codes, 3);
            const __m128i lo =
                _mm_blendv_epi8(_mm_shuffle_epi8(t0, codes), _mm_shuffle_epi8(t1, codes), m4);
            const __m128i hi =
                _mm_blendv_epi8(_mm_shuffle_epi8(t2, codes), _mm_shuffle_epi8(t3, codes), m4);
            _mm_storeu_si128((__m128i*)(void*)(out + i),
                             _mm_blendv_epi8(lo, hi, _mm_slli_epi16(codes, 2)));
            i += 16;
        }
    }
#endif /* AVX-512 VBMI, else AVX2 + SSE */
#endif
    // Generic scalar bit-reader (any D).
    {
        uint64_t bitpos = i * (uint64_t)D;
        const uint32_t m = (1U << D) - 1U;
        const size_t fbytes = (size_t)(((uint64_t)n * (uint64_t)D + 7U) >> 3);
        while (i < n) {
            const size_t byte = (size_t)(bitpos >> 3);
            uint32_t w = bits[byte];
            if (byte + 1 < fbytes) w |= (uint32_t)bits[byte + 1] << 8;
            if (D > 8 && byte + 2 < fbytes) w |= (uint32_t)bits[byte + 2] << 16;
            out[i++] = c2s[(w >> (bitpos & 7U)) & m];
            bitpos += (uint64_t)D;
        }
    }
}
 
static ZXC_ALWAYS_INLINE void zxc_pivco_emit_leaf_pair(uint8_t* RESTRICT out, const size_t n,
                                                       const uint8_t sym0, const uint8_t sym1,
                                                       const uint8_t* RESTRICT bits) {
    size_t i = 0;
#if defined(ZXC_USE_NEON64)
    const uint8x16_t vsym0 = vdupq_n_u8(sym0);
    const uint8x16_t vdelta = vdupq_n_u8((uint8_t)(sym0 ^ sym1));
    static const uint8_t rep_idx[16] = {0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1};
    static const uint8_t bit_sel[16] = {1, 2, 4, 8, 16, 32, 64, 128, 1, 2, 4, 8, 16, 32, 64, 128};
    const uint8x16_t vrep = vld1q_u8(rep_idx);
    const uint8x16_t vsel = vld1q_u8(bit_sel);
    while (i + 16 <= n) {
        uint8x16_t ctrl = vdupq_n_u8(0);
        ctrl = vsetq_lane_u8(bits[i >> 3], ctrl, 0);
        ctrl = vsetq_lane_u8(bits[(i >> 3) + 1], ctrl, 1);
        const uint8x16_t rep = vqtbl1q_u8(ctrl, vrep);
        const uint8x16_t mask = vtstq_u8(rep, vsel);
        vst1q_u8(out + i, veorq_u8(vsym0, vandq_u8(vdelta, mask)));
        i += 16;
    }
#elif defined(ZXC_USE_NEON32)
    const uint8x8_t vsym0 = vdup_n_u8(sym0);
    const uint8x8_t vdelta = vdup_n_u8((uint8_t)(sym0 ^ sym1));
    static const uint8_t bit_sel8[8] = {1, 2, 4, 8, 16, 32, 64, 128};
    const uint8x8_t vsel = vld1_u8(bit_sel8);
    while (i + 8 <= n) {
        const uint8x8_t mask = vtst_u8(vdup_n_u8(bits[i >> 3]), vsel);
        vst1_u8(out + i, veor_u8(vsym0, vand_u8(vdelta, mask)));
        i += 8;
    }
#elif defined(ZXC_USE_AVX512) || defined(ZXC_USE_AVX2) || \
    (defined(ZXC_USE_SSE2) && defined(__SSSE3__))
#if defined(ZXC_USE_AVX512)
    /* The control word is already the blend mask; the masked store takes the
     * tail, hence nothing narrower below. */
    {
        const __m512i z0 = _mm512_set1_epi8((char)sym0);
        const __m512i z1 = _mm512_set1_epi8((char)sym1);
        while (i < n) {
            const size_t left = n - i;
            const size_t cb = (left + 7) >> 3;
            uint64_t ctrl = 0;
            ZXC_MEMCPY(&ctrl, bits + (i >> 3), cb < 8 ? cb : 8); /* not past the run */
            const __mmask64 sm = (left >= 64) ? ~(__mmask64)0 : (((__mmask64)1 << left) - 1);
            _mm512_mask_storeu_epi8((void*)(out + i), sm,
                                    _mm512_mask_blend_epi8((__mmask64)ctrl, z0, z1));
            i += 64;
        }
        if (i > n) i = n;
    }
#else
#if defined(ZXC_USE_AVX2)
    /* Broadcasting the control dword lets each 128-bit half pick its own two
     * bytes out of it. */
    {
        const __m256i y0 = _mm256_set1_epi8((char)sym0);
        const __m256i ydelta = _mm256_set1_epi8((char)(uint8_t)(sym0 ^ sym1));
        const __m256i yrep = _mm256_setr_epi8(0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2,
                                              2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3);
        const __m256i ysel =
            _mm256_setr_epi8(1, 2, 4, 8, 16, 32, 64, (char)128, 1, 2, 4, 8, 16, 32, 64, (char)128,
                             1, 2, 4, 8, 16, 32, 64, (char)128, 1, 2, 4, 8, 16, 32, 64, (char)128);
        while (i + 32 <= n) {
            uint32_t c4;
            ZXC_MEMCPY(&c4, bits + (i >> 3), 4);
            const __m256i rep = _mm256_shuffle_epi8(_mm256_set1_epi32((int)c4), yrep);
            const __m256i mask = _mm256_cmpeq_epi8(_mm256_and_si256(rep, ysel), ysel);
            _mm256_storeu_si256((__m256i*)(void*)(out + i),
                                _mm256_xor_si256(y0, _mm256_and_si256(ydelta, mask)));
            i += 32;
        }
    }
#endif
    const __m128i vsym0 = _mm_set1_epi8((char)sym0);
    const __m128i vdelta = _mm_set1_epi8((char)(uint8_t)(sym0 ^ sym1));
    const __m128i vrep = _mm_setr_epi8(0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1);
    const __m128i vsel =
        _mm_setr_epi8(1, 2, 4, 8, 16, 32, 64, (char)128, 1, 2, 4, 8, 16, 32, 64, (char)128);
    while (i + 16 <= n) {
        const __m128i ctrl = _mm_cvtsi32_si128(bits[i >> 3] | ((int)bits[(i >> 3) + 1] << 8));
        const __m128i rep = _mm_shuffle_epi8(ctrl, vrep);
        // lanes where the selected bit is set -> 0xFF
        const __m128i mask = _mm_cmpeq_epi8(_mm_and_si128(rep, vsel), vsel);
        _mm_storeu_si128((__m128i*)(void*)(out + i),
                         _mm_xor_si128(vsym0, _mm_and_si128(vdelta, mask)));
        i += 16;
    }
#endif /* AVX-512, else AVX2 + SSE */
#endif
    while (i < n) {
        const int bit = (bits[i >> 3] >> (i & 7)) & 1;
        out[i++] = bit ? sym1 : sym0;
    }
}
 
static int zxc_pivco_decode_core(const uint8_t* RESTRICT payload, const size_t payload_size,
                                 uint8_t* RESTRICT dst, const size_t n,
                                 const zxc_pivco_tree_t* RESTRICT t, uint8_t* RESTRICT scratch) {
    if (UNLIKELY(n == 0)) return ZXC_ERROR_CORRUPT_DATA;
 
    // Pass 1, in wire order: node counts from the partition bits, and where
    // each depth's runs start. Leaf-pair children (emitted by their parent) and
    // a flat root's inside get count 0: empty runs, nothing to flag.
    uint32_t count[ZXC_PIVCO_NODE_SLOTS];
    const uint8_t* lvl_run[ZXC_HUF_MAX_CODE_LEN_ULTRA];
    const int slot = t->node_base[t->max_depth + 1];
    count[0] = (uint32_t)n;
    count[slot] = 0;
    const uint8_t* p = payload;
    const uint8_t* const pend = payload + payload_size;
    for (int d = 0; d < t->max_depth; d++) {
        lvl_run[d] = p;
        const int L = zxc_pivco_leaves(t, d);
        const int N = t->node_base[d + 1] - t->node_base[d];
        const int Lc = zxc_pivco_leaves(t, d + 1);
        const uint32_t* const cnt = count + t->node_base[d];
        uint32_t* const ccnt = count + t->node_base[d + 1];
        const uint8_t* const fd = t->flat_d + t->node_base[d];
        for (int i = L; i < N; i++) {
            const uint32_t c = cnt[i];
            const int k2 = 2 * (i - L); /* first child */
            if (fd[i]) {
                // Packed-code run: no partition below, nothing to popcount.
                const size_t fbytes = zxc_pivco_run_bytes(c, fd[i]);
                if (UNLIKELY((size_t)(pend - p) < fbytes)) return ZXC_ERROR_CORRUPT_DATA;
                p += fbytes;
                ccnt[k2] = 0;
                ccnt[k2 + 1] = 0;
                continue;
            }
            const size_t nbytes = zxc_pivco_run_bytes(c, 0); /* merge node: c partition bits */
            if (UNLIKELY((size_t)(pend - p) < nbytes)) return ZXC_ERROR_CORRUPT_DATA;
            // popcount of the c valid bits = the right child's count. All but the
            // last byte counted fast, the last one with its padding masked - if the
            // 8-byte loop swallowed it unmasked, padding bits would inflate `ones`.
            uint32_t ones = 0;
            if (nbytes) {
                const size_t full = nbytes - 1;
                size_t k = 0;
                for (; k + 8 <= full; k += 8) {
                    uint64_t w;
                    ZXC_MEMCPY(&w, p + k, 8);
                    ones += (uint32_t)zxc_pivco_popcnt64(w);
                }
                for (; k < full; k++) ones += (uint32_t)zxc_pivco_popcnt32(p[k]);
                uint8_t last = p[full];
                if (c & 7U) last &= (uint8_t)((1U << (c & 7U)) - 1U);
                ones += (uint32_t)zxc_pivco_popcnt32(last);
            }
            p += nbytes;
            if (UNLIKELY(ones > c)) return ZXC_ERROR_CORRUPT_DATA;
            const int leaf_pair = k2 + 1 < Lc;
            ccnt[k2] = leaf_pair ? 0 : c - ones;
            ccnt[k2 + 1] = leaf_pair ? 0 : ones;
        }
    }
    // The slot stands for the lone-symbol root's missing child: no 1 bits.
    if (UNLIKELY(count[slot] != 0)) return ZXC_ERROR_CORRUPT_DATA;
 
    // Pass 2, bottom-up. Runs and children follow position order, so two
    // running offsets replace per-node tables.
    for (int d = t->max_depth; d >= 0; d--) {
        uint8_t* const buf_d = (d & 1) ? scratch : dst;
        const uint8_t* const buf_c = (d & 1) ? dst : scratch; /* children at d+1 (read-only) */
        const int L = zxc_pivco_leaves(t, d);
        const int N = t->node_base[d + 1] - t->node_base[d];
        const uint32_t* const cnt = count + t->node_base[d];
        const uint8_t* const sym = t->syms + t->leaf_base[d];
        size_t off = 0;
        for (int i = 0; i < L; i++) {
            if (cnt[i]) ZXC_MEMSET(buf_d + off, sym[i], cnt[i]);
            off += cnt[i];
        }
        if (d == t->max_depth) continue; /* leaves only, as in pass 1 */
        const uint8_t* run = lvl_run[d];
        const uint32_t* const ccnt = count + t->node_base[d + 1];
        const uint8_t* const csym = t->syms + t->leaf_base[d + 1];
        const int Lc = zxc_pivco_leaves(t, d + 1);
        const uint8_t* const fd = t->flat_d + t->node_base[d];
        size_t coff = 0;
        for (int i = L; i < N; i++) {
            const uint32_t c = cnt[i];
            const int k2 = 2 * (i - L);
            if (c) {
                if (fd[i]) {
                    const int D = fd[i];
                    const int leaves = zxc_pivco_flat_leaves(t, t->node_base[d] + i, d, D);
                    zxc_pivco_unpack_flat(buf_d + off, c, D, run, t->syms + leaves);
                } else if (k2 + 1 < Lc) {
                    zxc_pivco_emit_leaf_pair(buf_d + off, c, csym[k2], csym[k2 + 1], run);
                } else {
                    const uint32_t nl = ccnt[k2];
                    zxc_pivco_merge(buf_d + off, buf_c + coff, nl, c - nl, run);
                }
            }
            run += zxc_pivco_run_bytes(c, fd[i]);
            off += c;
            coff += ccnt[k2] + ccnt[k2 + 1];
        }
    }
    return ZXC_OK;
}
 
int zxc_huf_decode_section(const uint8_t* RESTRICT payload, const size_t payload_size,
                           uint8_t* RESTRICT dst, const size_t n, uint8_t* RESTRICT scratch) {
    if (UNLIKELY(payload_size < (size_t)ZXC_HUF_TABLE_SIZE)) return ZXC_ERROR_CORRUPT_DATA;
    uint8_t code_len[ZXC_HUF_NUM_SYMBOLS];
    const int rc = zxc_huf_unpack_lengths(payload, code_len);
    if (UNLIKELY(rc != ZXC_OK)) return rc;
    zxc_pivco_tree_t t;
    if (UNLIKELY(zxc_pivco_tree_build(code_len, &t) != 0)) return ZXC_ERROR_CORRUPT_DATA;
    return zxc_pivco_decode_core(payload + ZXC_HUF_TABLE_SIZE, payload_size - ZXC_HUF_TABLE_SIZE,
                                 dst, n, &t, scratch);
}
 
int zxc_huf_decode_section_dict(const uint8_t* RESTRICT payload, const size_t payload_size,
                                uint8_t* RESTRICT dst, const size_t n,
                                const zxc_pivco_tree_t* RESTRICT tree, uint8_t* RESTRICT scratch) {
    return zxc_pivco_decode_core(payload, payload_size, dst, n, tree, scratch);
}