gwseq-io 0.2.1

//! rANS N×16 — CRAM 3.1's general-purpose entropy codec.
//!
//! Reference: `docs/cram_codecs_v3.1.md` §3.
//!
//! Four things are bundled into one stream here, and the leading flag byte says
//! which are present: an order-0 or order-1 rANS coder, run-length encoding,
//! bit packing, and a four-way byte stripe. They compose in a fixed order —
//! meta-data is read pack-then-RLE, and the transforms are undone
//! RLE-then-unpack — which is the part that is easy to get backwards and
//! impossible to notice, because a stream decoded in the wrong order is still
//! bytes.
//!
//! # Two corrections to the pseudocode
//!
//! §3.7's `RansDecodeNx16` reads `if flags AND NoSize != 0 then len <-
//! ReadUint7()`, which is inverted: the field table one page earlier says the
//! length is present *unless* `NoSize` is set, and `NoSize` exists precisely so
//! that a stripe's sub-streams can leave it out. Read the other way, every
//! ordinary block would take its length from whatever byte followed the flags.
//!
//! §1.4's five-byte `ReadITF8` ends `v + (ReadUint8() >> 4)`, taking the high
//! nibble of the last byte. It is the low nibble — the specification's own EOF
//! container encodes `-1` as `ff ff ff ff 0f`, which is `-1` only if the last
//! four bits are the ones kept. See [`super::super::container::read_itf8`].

use crate::error::{Error, Result};

use super::{ByteReader, PackMeta, MAX_CODEC_LEN};

/// The flag byte, §3.7.
mod flag {
    pub const ORDER: u8 = 1;
    pub const N32: u8 = 4;
    pub const STRIPE: u8 = 8;
    pub const NO_SIZE: u8 = 16;
    pub const CAT: u8 = 32;
    pub const RLE: u8 = 64;
    pub const PACK: u8 = 128;
}

/// The rANS lower bound, and the point renormalisation refills at.
const RANS_L: u32 = 1 << 15;

/// How deep a stripe may nest.
///
/// A stripe's sub-streams are themselves complete rANS N×16 streams, so a
/// crafted block can ask for a stripe of stripes of stripes. Real data never
/// nests: the encoder stripes once, over four byte positions of a 32-bit word.
const MAX_STRIPE_DEPTH: u32 = 4;

/// Decode a complete rANS N×16 block.
pub fn decode(data: &[u8], path: &str, offset: u64) -> Result<Vec<u8>> {
    let mut reader = ByteReader::new(data, path, offset);
    decode_stream(&mut reader, None, 0)
}

/// One stream, which may be a stripe's sub-stream.
///
/// `known_len` is the length a stripe computed for this sub-stream; it is used
/// only when the `NoSize` flag says the stream does not carry its own.
fn decode_stream(
    reader: &mut ByteReader<'_>,
    known_len: Option<usize>,
    depth: u32,
) -> Result<Vec<u8>> {
    let flags = reader.u8()?;

    let mut len = if flags & flag::NO_SIZE == 0 {
        reader.length()?
    } else {
        known_len.ok_or_else(|| {
            Error::corrupt(
                reader.path(),
                reader.offset(),
                "a rans4x16 stream sets NoSize but nothing outside it knows the size",
            )
        })?
    };

    if flags & flag::STRIPE != 0 {
        return decode_stripe(reader, len, depth);
    }

    let n_states = if flags & flag::N32 != 0 { 32 } else { 4 };

    // Meta-data is read pack-then-RLE, and each one rewrites `len` to the
    // length of what the *next* stage down produces. The two saved lengths are
    // what the transforms expand back out to on the way up.
    let mut pack = None;
    if flags & flag::PACK != 0 {
        let unpacked_len = len;
        let meta = PackMeta::read(reader)?;
        len = meta.packed_len;
        pack = Some((meta, unpacked_len));
    }
    let mut rle = None;
    if flags & flag::RLE != 0 {
        let expanded_len = len;
        let meta = RleMeta::read(reader, n_states)?;
        len = meta.pre_expansion_len;
        rle = Some((meta, expanded_len));
    }

    if len > MAX_CODEC_LEN {
        return Err(Error::corrupt(
            reader.path(),
            reader.offset(),
            format!("a rans4x16 stream declares {len} bytes, past this reader's ceiling"),
        ));
    }

    let mut data = if flags & flag::CAT != 0 {
        reader.take(len)?.to_vec()
    } else if flags & flag::ORDER != 0 {
        decode_order_1(reader, len, n_states)?
    } else {
        decode_order_0(reader, len, n_states)?
    };

    if let Some((meta, expanded_len)) = rle {
        data = meta.expand(&data, expanded_len, reader.path(), reader.offset())?;
    }
    if let Some((meta, unpacked_len)) = pack {
        data = meta.unpack(&data, unpacked_len, reader.path(), reader.offset())?;
    }
    Ok(data)
}

/// §3.6: `N` sub-streams holding every `N`th byte, transposed back on the way
/// out.
fn decode_stripe(reader: &mut ByteReader<'_>, len: usize, depth: u32) -> Result<Vec<u8>> {
    if depth >= MAX_STRIPE_DEPTH {
        return Err(Error::corrupt(
            reader.path(),
            reader.offset(),
            "rans4x16 stripes nested past this reader's limit",
        ));
    }
    let n = reader.u8()? as usize;
    if n == 0 {
        return Err(Error::corrupt(
            reader.path(),
            reader.offset(),
            "a rans4x16 stripe of zero sub-streams",
        ));
    }
    let mut lengths = Vec::with_capacity(n);
    for _ in 0..n {
        lengths.push(reader.length()?);
    }

    let mut out = vec![0u8; len];
    for (j, compressed_len) in lengths.into_iter().enumerate() {
        // The last sub-streams are one byte shorter when the total is not a
        // multiple of N.
        let sub_len = len / n + usize::from(len % n > j);
        let bytes = reader.take(compressed_len)?;
        let mut sub = ByteReader::new(bytes, reader.path(), reader.offset());
        let decoded = decode_stream(&mut sub, Some(sub_len), depth + 1)?;
        if decoded.len() != sub_len {
            return Err(Error::corrupt(
                reader.path(),
                reader.offset(),
                format!(
                    "a rans4x16 stripe sub-stream gave {} bytes where {sub_len} were due",
                    decoded.len()
                ),
            ));
        }
        for (i, byte) in decoded.into_iter().enumerate() {
            out[i * n + j] = byte;
        }
    }
    Ok(out)
}

/// A symbol table: frequencies, cumulative frequencies, and the reverse lookup
/// from a cumulative frequency back to its symbol.
struct SymbolTable {
    freq: [u32; 256],
    cumulative: [u32; 256],
    /// `1 << bits` entries, so finding a symbol is an index rather than a
    /// search. Built once per context that is actually used — for base or
    /// quality data that is a handful of the 256, not all of them.
    lookup: Vec<u8>,
}

impl SymbolTable {
    /// Normalise to `bits`, accumulate, and build the reverse lookup.
    fn build(mut freq: [u32; 256], bits: u32, path: &str, offset: u64) -> Result<Self> {
        // Summed wide: each frequency is a `uint7` a file chose, so 256 of them
        // overflow a `u32` long before they reach anything meaningful.
        let total: u64 = freq.iter().map(|f| u64::from(*f)).sum();
        let target = 1u64 << bits;
        if total == 0 {
            return Err(Error::corrupt(
                path,
                offset,
                "a rans4x16 frequency table whose frequencies are all zero",
            ));
        }
        if total > target {
            return Err(Error::corrupt(
                path,
                offset,
                format!("rans4x16 frequencies summing to {total}, past the {target} they must fit"),
            ));
        }
        // §3.1: frequencies must sum to a power of two, but not necessarily to
        // this one, so they are scaled up until they do.
        let mut scaled = total;
        let mut shift = 0;
        while scaled < target {
            scaled *= 2;
            shift += 1;
        }
        if scaled != target {
            return Err(Error::corrupt(
                path,
                offset,
                format!("rans4x16 frequencies summing to {total}, which is not a power of two"),
            ));
        }
        // Every frequency is at most `total`, and `total << shift` is `target`,
        // so this cannot overflow — but only because of the two checks above.
        let target = target as u32;
        for f in freq.iter_mut() {
            *f <<= shift;
        }

        let mut cumulative = [0u32; 256];
        let mut lookup = vec![0u8; target as usize];
        let mut running = 0u32;
        for symbol in 0..256 {
            cumulative[symbol] = running;
            let f = freq[symbol] as usize;
            lookup[running as usize..running as usize + f].fill(symbol as u8);
            running += freq[symbol];
        }
        Ok(Self {
            freq,
            cumulative,
            lookup,
        })
    }
}

/// §3.1's `ReadAlphabet`: an ascending symbol list, run-length encoded where
/// the symbols are consecutive.
///
/// The symbol is an `i32` and not a `u8`, and that is the whole subtlety of
/// this function. The list is terminated by reading a symbol of zero, and a
/// consecutive symbol is spotted by `s == last + 1` — so when the alphabet runs
/// up to 255, the terminating zero must compare against **256** and not wrap
/// round to match it. Read as bytes, an alphabet ending at 255 swallows the
/// byte after its terminator as a run length, and every frequency after that is
/// off by one position. The result still decodes; it just decodes to the wrong
/// bytes. Found on a real 3.1 file whose `YS:c` tag alphabet ends at 255.
fn read_alphabet(reader: &mut ByteReader<'_>) -> Result<Vec<u8>> {
    let mut alphabet = Vec::new();
    let mut symbol = i32::from(reader.u8()?);
    let mut last = symbol;
    let mut run = 0u32;
    loop {
        if symbol > 255 {
            // A run that walked off the end of the alphabet. Nothing more can
            // be read, and the symbols gathered so far are the alphabet.
            break;
        }
        alphabet.push(symbol as u8);
        if alphabet.len() > 256 {
            return Err(Error::corrupt(
                reader.path(),
                reader.offset(),
                "a rans4x16 alphabet of more than 256 symbols",
            ));
        }
        if run > 0 {
            run -= 1;
            symbol += 1;
        } else {
            symbol = i32::from(reader.u8()?);
            if symbol == last + 1 {
                run = u32::from(reader.u8()?);
            }
        }
        last = symbol;
        if symbol == 0 {
            break;
        }
    }
    Ok(alphabet)
}

/// §3.1's `ReadFrequenciesNx16_0`, at a fixed twelve bits.
fn read_frequencies_0(reader: &mut ByteReader<'_>) -> Result<SymbolTable> {
    let alphabet = read_alphabet(reader)?;
    let mut freq = [0u32; 256];
    for symbol in alphabet {
        freq[symbol as usize] = reader.uint7()?;
    }
    SymbolTable::build(freq, 12, reader.path(), reader.offset())
}

/// §3.1's `ReadFrequenciesNx16_1`: one alphabet for every context, a
/// zero-run-length-encoded table, and the whole thing optionally packed into an
/// order-0 stream of its own.
fn read_frequencies_1(reader: &mut ByteReader<'_>) -> Result<(Vec<Option<SymbolTable>>, u32)> {
    let comp = reader.u8()?;
    let bits = u32::from(comp >> 4);
    if bits != 10 && bits != 12 {
        return Err(Error::corrupt(
            reader.path(),
            reader.offset(),
            format!("a rans4x16 order-1 frequency table scaled to {bits} bits, not 10 or 12"),
        ));
    }

    // The table may itself be an order-0 rANS stream — a raw one, with no flag
    // byte of its own — in which case everything below reads from that.
    let tables = if comp & 1 != 0 {
        let raw_len = reader.length()?;
        let compressed_len = reader.length()?;
        let bytes = reader.take(compressed_len)?;
        let mut inner = ByteReader::new(bytes, reader.path(), reader.offset());
        let table_bytes = decode_order_0(&mut inner, raw_len, 4)?;
        let mut table_reader = ByteReader::new(&table_bytes, reader.path(), reader.offset());
        read_frequencies_1_tables(&mut table_reader, bits)?
    } else {
        read_frequencies_1_tables(reader, bits)?
    };
    Ok((tables, bits))
}

fn read_frequencies_1_tables(
    reader: &mut ByteReader<'_>,
    bits: u32,
) -> Result<Vec<Option<SymbolTable>>> {
    let alphabet = read_alphabet(reader)?;
    let mut tables: Vec<Option<SymbolTable>> = (0..256).map(|_| None).collect();
    for &context in &alphabet {
        let mut freq = [0u32; 256];
        let mut run = 0u32;
        for &symbol in &alphabet {
            if run > 0 {
                run -= 1;
                continue;
            }
            let f = reader.uint7()?;
            freq[symbol as usize] = f;
            if f == 0 {
                run = u32::from(reader.u8()?);
            }
        }
        // A context that never occurs has an all-zero row, and the alphabet is
        // one list for every context, so this is ordinary rather than corrupt:
        // a symbol that only ever appears last is in the alphabet and follows
        // nothing. It stays `None`, and reaching it is what is an error.
        if freq.iter().any(|f| *f > 0) {
            tables[context as usize] = Some(SymbolTable::build(
                freq,
                bits,
                reader.path(),
                reader.offset(),
            )?);
        }
    }
    Ok(tables)
}

/// Renormalise one state: below the lower bound it takes on sixteen more bits.
#[inline]
fn renorm(state: u32, reader: &mut ByteReader<'_>) -> Result<u32> {
    if state < RANS_L {
        return Ok((state << 16) + u32::from(reader.u16()?));
    }
    Ok(state)
}

/// §3.2: `N` states, each taking every `N`th byte of the output.
fn decode_order_0(reader: &mut ByteReader<'_>, len: usize, n: usize) -> Result<Vec<u8>> {
    let table = read_frequencies_0(reader)?;
    let mut states = Vec::with_capacity(n);
    for _ in 0..n {
        states.push(reader.u32()?);
    }
    let mask = (1u32 << 12) - 1;
    let mut out = Vec::with_capacity(len.min(1 << 20));
    for i in 0..len {
        let j = i % n;
        let state = states[j];
        let symbol = table.lookup[(state & mask) as usize];
        let f = table.freq[symbol as usize];
        let c = table.cumulative[symbol as usize];
        out.push(symbol);
        states[j] = renorm(f * (state >> 12) + (state & mask) - c, reader)?;
    }
    Ok(out)
}

/// §3.3: `N` states over *contiguous* regions rather than interleaved bytes, so
/// each keeps its own previous symbol as its context.
fn decode_order_1(reader: &mut ByteReader<'_>, len: usize, n: usize) -> Result<Vec<u8>> {
    let (tables, bits) = read_frequencies_1(reader)?;
    let mut states = Vec::with_capacity(n);
    for _ in 0..n {
        states.push(reader.u32()?);
    }
    let mut contexts = vec![0u8; n];
    let mask = (1u32 << bits) - 1;

    let missing = |context: u8, reader: &ByteReader<'_>| {
        Error::corrupt(
            reader.path(),
            reader.offset(),
            format!("a rans4x16 order-1 stream reaching context {context}, which its table omits"),
        )
    };

    let mut out = vec![0u8; len];
    let stride = len / n;
    for i in 0..stride {
        for j in 0..n {
            let table = tables[contexts[j] as usize]
                .as_ref()
                .ok_or_else(|| missing(contexts[j], reader))?;
            let state = states[j];
            let symbol = table.lookup[(state & mask) as usize];
            let f = table.freq[symbol as usize];
            let c = table.cumulative[symbol as usize];
            out[i + j * stride] = symbol;
            states[j] = renorm(f * (state >> bits) + (state & mask) - c, reader)?;
            contexts[j] = symbol;
        }
    }
    // The tail the stride does not cover belongs to the last state, which
    // carries on from the symbol it left off at.
    let last = n - 1;
    for slot in out.iter_mut().take(len).skip(stride * n) {
        let table = tables[contexts[last] as usize]
            .as_ref()
            .ok_or_else(|| missing(contexts[last], reader))?;
        let state = states[last];
        let symbol = table.lookup[(state & mask) as usize];
        let f = table.freq[symbol as usize];
        let c = table.cumulative[symbol as usize];
        *slot = symbol;
        states[last] = renorm(f * (state >> bits) + (state & mask) - c, reader)?;
        contexts[last] = symbol;
    }
    Ok(out)
}

/// §3.4's run-length transform: which symbols carry a run, and the runs.
struct RleMeta {
    /// Whether a symbol is followed by a run length.
    has_run: [bool; 256],
    /// The run lengths, `uint7` encoded, read as the expansion walks.
    runs: Vec<u8>,
    /// How many bytes the entropy coder below produces.
    pre_expansion_len: usize,
}

impl RleMeta {
    fn read(reader: &mut ByteReader<'_>, n: usize) -> Result<Self> {
        let meta_len = reader.length()?;
        let pre_expansion_len = reader.length()?;
        // The bottom bit says whether the meta-data is stored plainly; the rest
        // is its length doubled.
        let raw_meta_len = meta_len / 2;
        let meta = if meta_len & 1 != 0 {
            reader.take(raw_meta_len)?.to_vec()
        } else {
            let compressed_len = reader.length()?;
            let bytes = reader.take(compressed_len)?;
            let mut inner = ByteReader::new(bytes, reader.path(), reader.offset());
            decode_order_0(&mut inner, raw_meta_len, n)?
        };

        let mut meta = ByteReader::new(&meta, reader.path(), reader.offset());
        let count = meta.u8()?;
        // Zero means every symbol carries a run.
        let count = if count == 0 { 256 } else { count as usize };
        let mut has_run = [false; 256];
        for _ in 0..count {
            has_run[meta.u8()? as usize] = true;
        }
        Ok(Self {
            has_run,
            runs: meta.take(meta.remaining())?.to_vec(),
            pre_expansion_len,
        })
    }

    fn expand(&self, data: &[u8], len: usize, path: &str, offset: u64) -> Result<Vec<u8>> {
        let mut runs = ByteReader::new(&self.runs, path, offset);
        let mut out = Vec::with_capacity(len.min(1 << 20));
        for &symbol in data {
            if out.len() >= len {
                break;
            }
            if self.has_run[symbol as usize] {
                // A run of `run + 1`, the symbol itself included.
                let run = runs.uint7()? as usize;
                let take = (run + 1).min(len - out.len());
                out.resize(out.len() + take, symbol);
            } else {
                out.push(symbol);
            }
        }
        if out.len() != len {
            return Err(Error::corrupt(
                path,
                offset,
                format!(
                    "a rans4x16 run-length stream expanded to {} bytes where {len} were due",
                    out.len()
                ),
            ));
        }
        Ok(out)
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    /// An order-0 encoder, so the decoder round-trips against something rather
    /// than against a second reading of the same prose.
    ///
    /// rANS encodes backwards — the decoder's first symbol is the encoder's
    /// last — so this fills a buffer in reverse and reverses it at the end.
    fn encode_order_0(data: &[u8], n: usize) -> Vec<u8> {
        let mut freq = [0u32; 256];
        for &b in data {
            freq[b as usize] += 1;
        }
        normalise(&mut freq, 12);
        let mut cumulative = [0u32; 256];
        let mut running = 0;
        for s in 0..256 {
            cumulative[s] = running;
            running += freq[s];
        }

        let mut words: Vec<u16> = Vec::new();
        let mut states = vec![RANS_L; n];
        for i in (0..data.len()).rev() {
            let j = i % n;
            let s = data[i] as usize;
            let (f, c) = (freq[s], cumulative[s]);
            let max = ((RANS_L >> 12) << 16) * f;
            if states[j] >= max {
                words.push((states[j] & 0xffff) as u16);
                states[j] >>= 16;
            }
            states[j] = ((states[j] / f) << 12) + (states[j] % f) + c;
        }

        let mut out = Vec::new();
        // The flag byte has to say which interleaving was used: encoding with
        // 32 states and declaring 4 decodes to plausible garbage.
        out.push(if n == 32 { flag::N32 } else { 0 });
        write_uint7(&mut out, data.len() as u32);
        write_frequencies(&mut out, &freq);
        for state in &states {
            out.extend_from_slice(&state.to_le_bytes());
        }
        for word in words.iter().rev() {
            out.extend_from_slice(&word.to_le_bytes());
        }
        out
    }

    /// The same for order-1, whose states run over contiguous regions and carry
    /// the previous symbol as context.
    fn encode_order_1(data: &[u8], n: usize) -> Vec<u8> {
        let len = data.len();
        let stride = len / n;
        assert!(stride > 0, "order-1 wants at least one byte per state");

        // Which symbol precedes each position, per the decoder's own layout.
        let context_of = |i: usize| -> u8 {
            if i >= stride * n {
                // The tail continues the last state.
                if i == stride * n {
                    if stride == 0 {
                        0
                    } else {
                        data[(n - 1) * stride + stride - 1]
                    }
                } else {
                    data[i - 1]
                }
            } else if i % stride == 0 {
                0
            } else {
                data[i - 1]
            }
        };

        let mut freq = [[0u32; 256]; 256];
        for i in 0..len {
            freq[context_of(i) as usize][data[i] as usize] += 1;
        }
        // The alphabet is one list for every context, so it holds every symbol
        // that occurs as a value as well as every one that occurs as a context.
        let used: Vec<u8> = freq
            .iter()
            .enumerate()
            .filter(|(symbol, row)| row.iter().any(|f| *f > 0) || data.contains(&(*symbol as u8)))
            .map(|(symbol, _)| symbol as u8)
            .collect();
        for row in freq.iter_mut() {
            if row.iter().any(|f| *f > 0) {
                normalise(row, 12);
            }
        }
        let mut cumulative = [[0u32; 256]; 256];
        for s in 0..256 {
            let mut running = 0;
            for t in 0..256 {
                cumulative[s][t] = running;
                running += freq[s][t];
            }
        }

        // Encode each state's region backwards, states in reverse order, which
        // mirrors the decoder walking them forwards.
        let mut words: Vec<u16> = Vec::new();
        let mut states = vec![RANS_L; n];
        let push = |i: usize, states: &mut Vec<u32>, words: &mut Vec<u16>, j: usize| {
            let ctx = context_of(i) as usize;
            let s = data[i] as usize;
            let (f, c) = (freq[ctx][s], cumulative[ctx][s]);
            let max = ((RANS_L >> 12) << 16) * f;
            if states[j] >= max {
                words.push((states[j] & 0xffff) as u16);
                states[j] >>= 16;
            }
            states[j] = ((states[j] / f) << 12) + (states[j] % f) + c;
        };
        for i in (stride * n..len).rev() {
            push(i, &mut states, &mut words, n - 1);
        }
        for i in (0..stride).rev() {
            for j in (0..n).rev() {
                push(i + j * stride, &mut states, &mut words, j);
            }
        }

        let mut out = Vec::new();
        out.push(flag::ORDER | if n == 32 { flag::N32 } else { 0 });
        write_uint7(&mut out, len as u32);
        out.push(12 << 4); // twelve-bit frequencies, table not itself compressed
        write_alphabet(&mut out, &used);
        for &ctx in &used {
            let mut run: usize = 0;
            let mut pending: Vec<u8> = Vec::new();
            for (k, &sym) in used.iter().enumerate() {
                if run > 0 {
                    run -= 1;
                    continue;
                }
                let f = freq[ctx as usize][sym as usize];
                write_uint7(&mut pending, f);
                if f == 0 {
                    // How many of the following symbols are also zero.
                    let mut zeros = 0u8;
                    for &next in &used[k + 1..] {
                        if freq[ctx as usize][next as usize] == 0 && zeros < 255 {
                            zeros += 1;
                        } else {
                            break;
                        }
                    }
                    pending.push(zeros);
                    run = zeros as usize;
                }
            }
            out.extend_from_slice(&pending);
        }
        for state in &states {
            out.extend_from_slice(&state.to_le_bytes());
        }
        for word in words.iter().rev() {
            out.extend_from_slice(&word.to_le_bytes());
        }
        out
    }

    fn normalise(freq: &mut [u32; 256], bits: u32) {
        let target = 1u32 << bits;
        let total: u32 = freq.iter().sum();
        if total == 0 {
            return;
        }
        // Scale to the target, giving every present symbol at least one count.
        let mut running = 0u32;
        let mut last = 0usize;
        for (symbol, count) in freq.iter_mut().enumerate() {
            if *count == 0 {
                continue;
            }
            let scaled = ((*count as u64 * target as u64) / total as u64).max(1) as u32;
            *count = scaled;
            running += scaled;
            last = symbol;
        }
        // Push the rounding error onto the last symbol present.
        while running > target {
            let take = (running - target).min(freq[last] - 1);
            freq[last] -= take;
            running -= take;
            if freq[last] == 1 {
                last = (0..256).rev().find(|s| freq[*s] > 1).unwrap_or(last);
            }
        }
        if running < target {
            freq[last] += target - running;
        }
    }

    fn write_uint7(out: &mut Vec<u8>, value: u32) {
        let mut groups = Vec::new();
        let mut value = value;
        loop {
            groups.push((value & 0x7f) as u8);
            value >>= 7;
            if value == 0 {
                break;
            }
        }
        for (i, group) in groups.iter().enumerate().rev() {
            out.push(if i == 0 { *group } else { group | 0x80 });
        }
    }

    fn write_alphabet(out: &mut Vec<u8>, symbols: &[u8]) {
        // The plain form: every symbol written, no consecutive run ever
        // declared, terminated by a zero. Legal, and what the decoder's
        // `rle == 0` path walks.
        let mut i = 0;
        while i < symbols.len() {
            out.push(symbols[i]);
            // A consecutive pair would be read as the start of a run, so the
            // run length has to be written when one occurs.
            if i + 1 < symbols.len() && symbols[i + 1] == symbols[i].wrapping_add(1) {
                let mut run = 0u8;
                while run < 255
                    && i + 1 + (run as usize) < symbols.len()
                    && symbols[i + 1 + run as usize] == symbols[i].wrapping_add(1).wrapping_add(run)
                {
                    run += 1;
                }
                out.push(symbols[i + 1]);
                out.push(run - 1);
                i += 1 + run as usize;
            } else {
                i += 1;
            }
        }
        out.push(0);
    }

    fn write_frequencies(out: &mut Vec<u8>, freq: &[u32; 256]) {
        let symbols: Vec<u8> = (0..256).filter(|s| freq[*s] > 0).map(|s| s as u8).collect();
        write_alphabet(out, &symbols);
        for &s in &symbols {
            write_uint7(out, freq[s as usize]);
        }
    }

    fn roundtrip_0(data: &[u8], n: usize) {
        let encoded = encode_order_0(data, n);
        let decoded = decode(&encoded, "test", 0).expect("decodes");
        assert_eq!(decoded, data, "order-0, N={n}");
    }

    fn roundtrip_1(data: &[u8], n: usize) {
        let encoded = encode_order_1(data, n);
        let decoded = decode(&encoded, "test", 0).expect("decodes");
        assert_eq!(decoded, data, "order-1, N={n}");
    }

    #[test]
    fn order_0_round_trips_at_both_interleavings() {
        let data: Vec<u8> = (0..5000u32).map(|i| (i * 7 % 41) as u8).collect();
        roundtrip_0(&data, 4);
        roundtrip_0(&data, 32);
    }

    #[test]
    fn order_0_round_trips_text() {
        let data =
            b"the quick brown fox jumps over the lazy dog, repeatedly and at length. ".repeat(40);
        roundtrip_0(&data, 4);
        roundtrip_0(&data, 32);
    }

    #[test]
    fn order_1_round_trips_at_both_interleavings() {
        let data = b"ACGTACGTTTTTACGNNNNNNACGTACGTACGGGGTTTACGT".repeat(120);
        roundtrip_1(&data, 4);
        roundtrip_1(&data, 32);
    }

    /// The length not dividing by the state count is the case the tail loop
    /// exists for, and the one an off-by-one hides in.
    #[test]
    fn order_1_round_trips_a_length_that_is_not_a_multiple_of_n() {
        let base = b"ACGTNACGTNACGTTTTACG".repeat(30);
        for extra in 0..8 {
            let data = &base[..base.len() - extra];
            roundtrip_1(data, 4);
        }
    }

    /// The alphabet's terminator is a symbol of zero, and consecutive symbols
    /// are spotted with `s == last + 1`. An alphabet reaching 255 puts those two
    /// rules in each other's way: read as bytes, `last + 1` wraps to zero and
    /// the terminator looks like the start of a run, so the next byte is eaten
    /// as a run length and every frequency after it shifts by one.
    ///
    /// The stream still decodes — to the wrong bytes. This is the shape of bug
    /// the whole module is careful about, and it was found on a real file's
    /// `YS:c` tag rather than by reading the specification again.
    #[test]
    fn an_alphabet_that_reaches_255_does_not_eat_the_byte_after_its_terminator() {
        // Symbol 0 and the run 226..=255, which is what the file had.
        let mut symbols = vec![0u8];
        symbols.extend(226..=255);
        let mut out = Vec::new();
        write_alphabet(&mut out, &symbols);
        let mut reader = ByteReader::new(&out, "test", 0);
        assert_eq!(read_alphabet(&mut reader).expect("an alphabet"), symbols);
        // And the terminator was the last thing consumed: nothing beyond it.
        assert_eq!(reader.remaining(), 0);

        // The same, round-tripped through a real stream.
        let data: Vec<u8> = (0..4000).map(|i| symbols[i % symbols.len()]).collect();
        roundtrip_0(&data, 4);
        roundtrip_1(&data, 4);
    }

    /// The other end of the same rule: an alphabet of exactly `[0]` terminates
    /// at once, and one of `[0, 1]` does declare a run.
    #[test]
    fn an_alphabet_starting_at_zero_round_trips() {
        for symbols in [vec![0u8], vec![0, 1], vec![0, 1, 2], vec![0, 5], vec![255]] {
            let mut out = Vec::new();
            write_alphabet(&mut out, &symbols);
            let mut reader = ByteReader::new(&out, "test", 0);
            assert_eq!(
                read_alphabet(&mut reader).expect("an alphabet"),
                symbols,
                "{symbols:?}"
            );
        }
    }

    #[test]
    fn a_single_symbol_stream_round_trips() {
        let data = vec![b'Q'; 1000];
        roundtrip_0(&data, 4);
        roundtrip_1(&data, 4);
    }

    /// `Cat` is the escape hatch for data too short to be worth coding.
    #[test]
    fn the_cat_flag_gives_the_bytes_back_verbatim() {
        let mut stream = vec![flag::CAT];
        write_uint7(&mut stream, 5);
        stream.extend_from_slice(b"hello");
        assert_eq!(decode(&stream, "test", 0).expect("decodes"), b"hello");
    }

    #[test]
    fn a_pack_of_one_symbol_expands_without_reading_any_data() {
        // nsym=1 means every value is the same and nothing was packed at all.
        let mut stream = vec![flag::PACK | flag::CAT];
        write_uint7(&mut stream, 6); // unpacked length
        stream.push(1); // one symbol
        stream.push(b'N');
        write_uint7(&mut stream, 0); // no packed bytes
        assert_eq!(decode(&stream, "test", 0).expect("decodes"), b"NNNNNN");
    }

    #[test]
    fn a_two_symbol_pack_unpacks_eight_values_to_the_byte() {
        let mut stream = vec![flag::PACK | flag::CAT];
        write_uint7(&mut stream, 8);
        stream.push(2);
        stream.extend_from_slice(b"AB");
        write_uint7(&mut stream, 1);
        // §3.5: "bits are unpacked from low to high", so 0b0000_1101 reads as
        // 1,0,1,1,0,0,0,0 and not as the byte written left to right.
        stream.push(0b0000_1101);
        assert_eq!(decode(&stream, "test", 0).expect("decodes"), b"BABBAAAA");
    }

    #[test]
    fn a_run_length_stream_expands_its_runs() {
        // Meta: one symbol carries runs ('N'), then the run lengths.
        let mut meta = vec![1u8, b'N'];
        write_uint7(&mut meta, 3); // a run of 3+1 Ns
        let mut stream = vec![flag::RLE | flag::CAT];
        write_uint7(&mut stream, 6); // expanded length: "ANNNNA"
        write_uint7(&mut stream, (meta.len() * 2 + 1) as u32); // stored plainly
        write_uint7(&mut stream, 3); // pre-expansion length: "ANA"
        stream.extend_from_slice(&meta);
        stream.extend_from_slice(b"ANA");
        assert_eq!(decode(&stream, "test", 0).expect("decodes"), b"ANNNNA");
    }

    #[test]
    fn a_stripe_transposes_its_sub_streams_back() {
        // Two sub-streams over five bytes: "acbdc" and "ABC" interleave to
        // "aAbBcC" — except the total is 5, so the last of stream 0 has no
        // partner.
        let sub0 = {
            let mut s = vec![flag::CAT | flag::NO_SIZE];
            s.extend_from_slice(b"abc");
            s
        };
        let sub1 = {
            let mut s = vec![flag::CAT | flag::NO_SIZE];
            s.extend_from_slice(b"AB");
            s
        };
        let mut stream = vec![flag::STRIPE];
        write_uint7(&mut stream, 5);
        stream.push(2);
        write_uint7(&mut stream, sub0.len() as u32);
        write_uint7(&mut stream, sub1.len() as u32);
        stream.extend_from_slice(&sub0);
        stream.extend_from_slice(&sub1);
        assert_eq!(decode(&stream, "test", 0).expect("decodes"), b"aAbBc");
    }

    #[test]
    fn a_stripe_nested_past_the_limit_is_refused_rather_than_recursed() {
        // A stripe whose sub-stream is a stripe, all the way down.
        let mut stream = Vec::new();
        for _ in 0..MAX_STRIPE_DEPTH + 2 {
            let mut outer = vec![flag::STRIPE];
            write_uint7(&mut outer, 4);
            outer.push(1);
            write_uint7(&mut outer, stream.len() as u32);
            outer.extend_from_slice(&stream);
            stream = outer;
        }
        assert!(decode(&stream, "test", 0).is_err());
    }

    #[test]
    fn truncated_streams_are_errors_rather_than_short_output() {
        let data: Vec<u8> = (0..2000u32).map(|i| (i % 37) as u8).collect();
        let encoded = encode_order_0(&data, 4);
        for cut in [1, 5, 20, encoded.len() / 2, encoded.len() - 1] {
            assert!(
                decode(&encoded[..cut], "test", 0).is_err(),
                "a stream cut to {cut} bytes decoded"
            );
        }
    }
}