fgf 2.0.0

//! Throughput harness for the vector kernels.
//!
//! No criterion: these kernels are memory-bound straight-line loops with no
//! branchy tail, so a warm loop and a wall-clock median is enough to see the
//! effects that matter — backend, buffer size relative to cache, and whether
//! the blocked multi-row shapes actually beat repeated single-row AXPY.
//!
//! ```sh
//! cargo bench --features internals --bench kernels
//! SIMD_BACKEND=v3     cargo bench --features internals --bench kernels # compare backends
//! SIMD_BACKEND=scalar cargo bench --features internals --bench kernels
//! ```

// Toolchain-drift lint (not in the MSRV); see `src/lib.rs`.
#![allow(unknown_lints, clippy::chunks_exact_to_as_chunks)]
use std::hint::black_box;
use std::time::{Duration, Instant};

use fgf::{
    FanPaar16, FanPaar32, FanPaar64, Gf8B, Gf16, Gf32, Gf64, Goldilocks, Mersenne31,
    QuadMersenne31, backend, fan_paar, gf8b, gf16, gf32, gf64, goldilocks, mersenne31, ops,
    quad_mersenne31,
};

fn noise(len: usize, seed: u64) -> Vec<u8> {
    let mut state = seed | 1;
    (0..len)
        .map(|_| {
            state = state
                .wrapping_mul(6_364_136_223_846_793_005)
                .wrapping_add(1);
            (state >> 33) as u8
        })
        .collect()
}

/// Format a nanosecond count the way `Duration`'s debug output scales its
/// unit, but from an `f64`: `Duration` cannot hold sub-nanosecond values, so
/// dividing one by a repetition count quantizes every estimate to whole
/// nanoseconds.
fn fmt_ns(ns: f64) -> String {
    if ns < 1_000.0 {
        format!("{ns:.2}ns")
    } else if ns < 1_000_000.0 {
        format!("{:.2}µs", ns / 1_000.0)
    } else {
        format!("{:.2}ms", ns / 1_000_000.0)
    }
}

/// Run `body` until it has been timed enough times to trust the median, and
/// report bytes per second over `bytes` of logical traffic per iteration.
///
/// Returns the median per-iteration time in nanoseconds so callers can print
/// a ratio between two shapes measured back to back in the same process.
fn bench(label: &str, bytes: usize, mut body: impl FnMut()) -> f64 {
    // Warm caches, branch predictors, and the lazy backend detection.
    for _ in 0..16 {
        body();
    }

    let mut samples = Vec::with_capacity(64);
    let deadline = Instant::now() + Duration::from_millis(250);
    while Instant::now() < deadline && samples.len() < 64 {
        let start = Instant::now();
        let reps = 32;
        for _ in 0..reps {
            body();
        }
        samples.push(start.elapsed().as_secs_f64() * 1e9 / reps as f64);
    }
    samples.sort_unstable_by(f64::total_cmp);
    let median = samples[samples.len() / 2];

    let gib_per_sec = bytes as f64 / (median / 1e9) / (1024.0 * 1024.0 * 1024.0);
    println!(
        "  {label:<44} {:>9}  {gib_per_sec:>7.2} GiB/s",
        fmt_ns(median)
    );
    median
}

/// Row lengths for the preparation crossover sweep: below one lane through
/// well past the point where the byte loop dominates table derivation.
const CROSSOVER_LENGTHS: &[usize] = &[
    16, 32, 64, 128, 256, 512, 1_024, 2_048, 4_096, 16_384, 65_536,
];

/// Where does preparing a coefficient stop paying for itself?
///
/// One-shot `mul_add` derives the backend's coefficient form on every call —
/// two broadcast words on GFNI, four nibble tables on a shuffle backend. The
/// `_with` form derives it once. The ratio printed here is the whole reason
/// `Coeff`/`CoeffVec` exist, and it is the number `BENCHMARKS.md` refers to.
fn bench_preparation_crossover() {
    println!("preparation crossover — one-shot vs prepared, by row length:");
    for &len in CROSSOVER_LENGTHS {
        let src = noise(len, 0xa00 + len as u64);
        let mut dst = noise(len, 0xb00 + len as u64);
        let coeff8 = gf8b::Elem::from_raw(0x53);
        let coeff16 = gf16::Elem::from_raw(0x53a7);
        let prepared8 = ops::Coeff::<Gf8B>::new(coeff8);
        let prepared16 = ops::Coeff::<Gf16>::new(coeff16);

        println!("  row {len} B:");
        let one8 = bench("  mul_add one-shot           gf8", len, || {
            ops::mul_add::<Gf8B>(black_box(&mut dst), coeff8, black_box(&src));
        });
        let with8 = bench("  mul_add prepared           gf8", len, || {
            ops::mul_add_with::<Gf8B>(black_box(&mut dst), &prepared8, black_box(&src));
        });
        let one16 = bench("  mul_add one-shot          gf16", len, || {
            ops::mul_add::<Gf16>(black_box(&mut dst), coeff16, black_box(&src));
        });
        let with16 = bench("  mul_add prepared          gf16", len, || {
            ops::mul_add_with::<Gf16>(black_box(&mut dst), &prepared16, black_box(&src));
        });
        println!(
            "    one-shot/prepared: gf8 {:.2}x  gf16 {:.2}x",
            one8 / with8,
            one16 / with16
        );
    }
    println!();
}

/// Payload lengths used by network-facing consumers.
///
/// The run around 1,200 bytes alternates exact 32-byte lanes with 16-byte
/// remainders so SIMD tail cliffs stay visible.
const NETWORK_LENGTHS: &[usize] = &[
    64, 256, 512, 1_152, 1_168, 1_184, 1_200, 1_216, 1_232, 1_248, 1_400,
];

fn bench_network_payloads() {
    println!("network-size GF(2^8) payloads:");
    for &len in NETWORK_LENGTHS {
        println!("  payload {len} B:");
        let src = noise(len, 0x700 + len as u64);
        let mut dst = noise(len, 0x800 + len as u64);

        bench("xor", len, || {
            ops::add_assign::<Gf8B>(black_box(&mut dst), black_box(&src));
        });
        bench("mul_add", len, || {
            ops::mul_add::<Gf8B>(
                black_box(&mut dst),
                gf8b::Elem::from_raw(0x53),
                black_box(&src),
            );
        });
        bench("mul_assign", len, || {
            ops::mul_assign::<Gf8B>(black_box(&mut dst), gf8b::Elem::from_raw(0x53));
        });

        for nrows in [4usize, 16] {
            let coeffs: Vec<_> = (0..nrows)
                .map(|row| gf8b::Elem::from_raw((row as u8).wrapping_mul(37).wrapping_add(2)))
                .collect();
            let mut rows = noise(len * nrows, 0x900 + nrows as u64);
            let label = format!("scatter ({nrows} rows)");
            bench(&label, len * nrows, || {
                ops::mul_add_scatter::<Gf8B>(
                    black_box(&mut rows),
                    len,
                    black_box(&coeffs),
                    black_box(&src),
                );
            });
        }
    }
    println!();
}

/// Long batches isolate the payload loop from the clock and scheduler.
fn bench_network(label: &str, bytes: usize, mut body: impl FnMut()) -> f64 {
    for _ in 0..64 {
        body();
    }
    let mut samples = [0.0; 64];
    for sample in &mut samples {
        let start = Instant::now();
        for _ in 0..512 {
            body();
        }
        *sample = start.elapsed().as_secs_f64() * 1e9 / 512.0;
    }
    samples.sort_unstable_by(f64::total_cmp);
    let median = samples[samples.len() / 2];
    println!(
        "  {label:<44} {:>9}  {:>7.2} GiB/s",
        fmt_ns(median),
        bytes as f64 / (median / 1e9) / (1024.0 * 1024.0 * 1024.0)
    );
    median
}

/// Pinned diagnostic: fixed allocation alignment separates row pitch from
/// the allocator's base address.
fn bench_network_alignment() {
    for &len in &[1152usize, 1168, 1184, 1200] {
        let src_storage = noise(len + 128, 0x700 + len as u64);
        let mut dst_storage = noise(len + 128, 0x800 + len as u64);
        for &off in &[0usize, 16, 32, 48] {
            let s = src_storage.as_ptr().align_offset(64) + off;
            let d = dst_storage.as_ptr().align_offset(64) + off;
            let src = &src_storage[s..s + len];
            let dst = &mut dst_storage[d..d + len];
            println!(
                "network diagnostic {len} B off={off} dst={}",
                (dst.as_ptr() as usize) & 63
            );
            bench_network("xor", len, || {
                ops::add_assign::<Gf8B>(black_box(&mut *dst), black_box(src))
            });
            bench_network("mul_add", len, || {
                ops::mul_add::<Gf8B>(
                    black_box(&mut *dst),
                    gf8b::Elem::from_raw(0x53),
                    black_box(src),
                )
            });
            bench_network("mul_assign", len, || {
                ops::mul_assign::<Gf8B>(black_box(&mut *dst), gf8b::Elem::from_raw(0x53))
            });
            for nrows in [4usize, 16] {
                let coeffs: Vec<_> = (0..nrows)
                    .map(|row| gf8b::Elem::from_raw((row as u8).wrapping_mul(37).wrapping_add(2)))
                    .collect();
                let mut storage = noise(len * nrows + 128, 0x900 + nrows as u64);
                let r = storage.as_ptr().align_offset(64) + off;
                println!(
                    "rows={nrows} base={}",
                    (storage[r..].as_ptr() as usize) & 63
                );
                let rows = &mut storage[r..r + len * nrows];
                bench_network(&format!("scatter {nrows}"), len * nrows, || {
                    ops::mul_add_scatter::<Gf8B>(
                        black_box(&mut *rows),
                        len,
                        black_box(&coeffs),
                        black_box(src),
                    )
                });
            }
        }
    }
}

/// Blocked multi-row GF(2^16) kernels measured against repeated single-row
/// AXPY, bypassing dispatch.
///
/// The SSSE3 and GFNI tiers dispatch the blocked gather. This harness times
/// each against repeated AXPY in one process through `internals`
/// token-proven entries.
#[cfg(all(feature = "simd", any(target_arch = "x86", target_arch = "x86_64")))]
fn bench_blocked_vs_axpy() {
    use fgf::internals::kernel::tables::{TowerCoeff, TowerTables};
    use fgf::internals::kernel::x86;
    use fgf::internals::kernel::{SimdToken, X64V2Token, X64V3GfniCryptoToken};

    // Genuine capability tokens, summoned once: each tier that the host
    // cannot prove is skipped rather than SIGILLed.
    let Some(ssse3) = X64V2Token::summon() else {
        return;
    };
    let gfni = X64V3GfniCryptoToken::summon();

    println!("blocked vs AXPY — direct GF(2^16) kernel calls (dispatch bypassed):");
    for &row_len in &[4 * 1024usize, 16 * 1024, 64 * 1024] {
        for &nsrc in &[2usize, 4, 8, 16] {
            let sources: Vec<Vec<u8>> = (0..nsrc)
                .map(|t| noise(row_len, 0xc00 + t as u64))
                .collect();
            let srcs: Vec<&[u8]> = sources.iter().map(Vec::as_slice).collect();
            let coeffs: Vec<gf16::Elem> = (0..nsrc)
                .map(|t| gf16::Elem::from_raw(((t * 7919) as u16).wrapping_add(3)))
                .collect();
            let mut dst = noise(row_len, 0xd00);
            let traffic = row_len * nsrc;

            println!("  gather {nsrc} sources x {} KiB:", row_len / 1024);
            let blocked = bench("  gather blocked           ssse3", traffic, || {
                x86::gf16::mul_add_gather_ssse3(
                    ssse3,
                    black_box(&mut dst),
                    &coeffs,
                    black_box(&srcs),
                );
            });
            let axpy = bench("  gather AXPY              ssse3", traffic, || {
                for (&coeff, &src) in coeffs.iter().zip(&srcs) {
                    x86::gf16::mul_add_ssse3(
                        ssse3,
                        black_box(&mut dst),
                        &TowerTables::new(coeff),
                        black_box(src),
                    );
                }
            });
            println!("    ssse3 blocked/AXPY: {:.2}x", axpy / blocked);

            if let Some(gfni) = gfni {
                let blocked = bench("  gather blocked            gfni", traffic, || {
                    x86::gf16::mul_add_gather_gfni(
                        gfni,
                        black_box(&mut dst),
                        &coeffs,
                        black_box(&srcs),
                    );
                });
                let axpy = bench("  gather AXPY               gfni", traffic, || {
                    for (&coeff, &src) in coeffs.iter().zip(&srcs) {
                        x86::gf16::mul_add_gfni(
                            gfni,
                            black_box(&mut dst),
                            TowerCoeff::new(coeff),
                            black_box(src),
                        );
                    }
                });
                println!("    gfni  blocked/AXPY: {:.2}x", axpy / blocked);
            }
        }
    }
    println!();
}

/// Destination sizes straddling the non-temporal store threshold
/// (`x86::NT_STORE_MIN`, 2 MiB): one below it as a control, two above.
const LARGE_DST_LENGTHS: &[usize] = &[1 << 20, 8 << 20, 32 << 20];

/// `mul_into` over destinations large enough for the store policy to matter.
///
/// `mul_into` never reads its destination, so an ordinary store still fetches
/// every line it overwrites. Past 2 MiB the x86 kernels switch to `vmovntdq`,
/// which skips that fetch and does not keep the destination cached. The
/// read-back row is the workload that trade pessimizes: encode, then consume
/// the result immediately. `mul_add` is the control — it reads its
/// destination, so it is unaffected either way.
fn bench_large_destination() {
    println!("large destinations — mul_into store policy:");
    for &len in LARGE_DST_LENGTHS {
        let src = noise(len, 0xd00);
        let mut dst = noise(len, 0xd01);
        let mib = len / (1024 * 1024);

        bench(&format!("{mib:3} MiB mul_into          gf8"), len, || {
            ops::mul_into::<Gf8B>(
                black_box(&mut dst),
                gf8b::Elem::from_raw(0x53),
                black_box(&src),
            );
        });
        bench(&format!("{mib:3} MiB mul_into         gf16"), len, || {
            ops::mul_into::<Gf16>(
                black_box(&mut dst),
                gf16::Elem::from_raw(0x53a7),
                black_box(&src),
            );
        });
        bench(&format!("{mib:3} MiB mul_into+read     gf8"), len, || {
            ops::mul_into::<Gf8B>(
                black_box(&mut dst),
                gf8b::Elem::from_raw(0x53),
                black_box(&src),
            );
            let mut acc = 0u64;
            for word in dst.chunks_exact(8) {
                acc ^= u64::from_le_bytes(word.try_into().unwrap());
            }
            black_box(acc);
        });
        bench(&format!("{mib:3} MiB mul_add           gf8"), len, || {
            ops::mul_add::<Gf8B>(
                black_box(&mut dst),
                gf8b::Elem::from_raw(0x53),
                black_box(&src),
            );
        });
    }
    println!();
}

/// Multi-row scatter with the destination at each 32-byte residue.
///
/// A 32-byte access at an odd multiple of 32 touches two cache lines, and a
/// scatter loads and stores every row of a group once per source window, so
/// the destination side is where that doubles. The GFNI and AVX2 kernels peel
/// up to 31 bytes per row group to align the rest of the pass; the skew rows
/// below are what set that policy, and an allocator-fresh buffer alone will
/// not show it because its residue is fixed.
fn bench_destination_alignment() {
    println!("destination alignment — multi-row scatter:");
    let nrows = 8;
    for &row_len in &[64 * 1024usize, 256 * 1024] {
        let src = noise(row_len, 0xe00);
        let coeffs8: Vec<_> = (0..nrows)
            .map(|j| gf8b::Elem::from_raw((j as u8).wrapping_mul(37).wrapping_add(2)))
            .collect();
        let coeffs16: Vec<_> = (0..nrows)
            .map(|j| gf16::Elem::from_raw((j as u16).wrapping_mul(9871).wrapping_add(2)))
            .collect();
        let mut backing = noise(row_len * nrows + 32, 0xe01);
        let kib = row_len / 1024;
        for skew in [0usize, 16] {
            let rows = &mut backing[skew..skew + row_len * nrows];
            let traffic = row_len * nrows;
            bench(
                &format!("{kib:4} KiB rows, skew {skew:2}   gf8"),
                traffic,
                || {
                    ops::mul_add_scatter::<Gf8B>(
                        black_box(rows),
                        row_len,
                        &coeffs8,
                        black_box(&src),
                    );
                },
            );
            bench(
                &format!("{kib:4} KiB rows, skew {skew:2}  gf16"),
                traffic,
                || {
                    ops::mul_add_scatter::<Gf16>(
                        black_box(rows),
                        row_len,
                        &coeffs16,
                        black_box(&src),
                    );
                },
            );
        }
    }
    println!();
}

/// Row lengths where a GF(2^16) coefficient's four nibble tables are a
/// visible share of the work, plus one length where they are not.
const SMALL_ROW_LENGTHS: &[usize] = &[64, 256, 1_024, 16_384];

/// Multi-row GF(2^16) shapes at row lengths where preparation dominates.
///
/// A GF(2^16) coefficient costs four derived nibble tables. Multi-row kernels
/// that rebuild them per `(term, row)` — as the wasm kernels used to — pay
/// that per row instead of once, which only shows up when the row is short
/// enough that the byte loop cannot hide it. Coefficient sets deliberately
/// include zeros and ones, the case the zero/one specializations exist for.
fn bench_small_row_shapes() {
    println!("small-row GF(2^16) multi-row shapes — preparation-dominated:");
    let nrows = 16;
    let nsrc = 8;
    for &row_len in SMALL_ROW_LENGTHS {
        let sources: Vec<Vec<u8>> = (0..nsrc)
            .map(|t| noise(row_len, 0xf00 + t as u64))
            .collect();
        let srcs: Vec<&[u8]> = sources.iter().map(Vec::as_slice).collect();
        // Every fourth coefficient is zero and every fifth is one.
        let coeff_at = |i: usize| match i % 5 {
            0 => gf16::Elem::ZERO,
            1 => gf16::Elem::ONE,
            _ => gf16::Elem::from_raw(((i * 7919) as u16) | 0x0100),
        };
        let row_coeffs: Vec<gf16::Elem> = (0..nrows).map(coeff_at).collect();
        let coeff_sets: Vec<Vec<gf16::Elem>> = (0..nsrc)
            .map(|t| (0..nrows).map(|j| coeff_at(t * 3 + j)).collect())
            .collect();
        let terms: Vec<(&[gf16::Elem], &[u8])> = coeff_sets
            .iter()
            .zip(&sources)
            .map(|(c, s)| (c.as_slice(), s.as_slice()))
            .collect();
        let mut rows = noise(row_len * nrows, 0x1000);
        let mut dst = noise(row_len, 0x1100);

        println!("  row {row_len} B, {nrows} rows, {nsrc} sources:");
        bench("  scatter                   gf16", row_len * nrows, || {
            ops::mul_add_scatter::<Gf16>(
                black_box(&mut rows),
                row_len,
                &row_coeffs,
                black_box(srcs[0]),
            );
        });
        bench("  gather                    gf16", row_len * nsrc, || {
            ops::mul_add_gather::<Gf16>(black_box(&mut dst), &row_coeffs[..nsrc], black_box(&srcs));
        });
        bench(
            "  matrix                    gf16",
            row_len * nrows * nsrc,
            || {
                ops::mul_add_matrix::<Gf16>(black_box(&mut rows), row_len, nrows, &terms);
            },
        );

        // Fused `mul_into` against the copy-then-scale it replaces: the
        // fused kernel touches the destination once.
        let coeff = gf16::Elem::from_raw(0x53a7);
        let fused = bench("  mul_into fused            gf16", row_len, || {
            ops::mul_into::<Gf16>(black_box(&mut dst), coeff, black_box(srcs[0]));
        });
        let copied = bench("  mul_into copy+scale       gf16", row_len, || {
            let dst = black_box(&mut dst);
            dst.copy_from_slice(black_box(srcs[0]));
            ops::mul_assign::<Gf16>(dst, coeff);
        });
        println!("    copy+scale/fused: {:.2}x", copied / fused);
    }
    println!();
}

/// Bit-packed GF(2) surface: xor (dispatched byte XOR), the portable AND
/// word loop, the popcount folds, and a masked range — against the
/// byte-per-element GF(2^8) add as the density control. The packing itself
/// is the expected win (8x the elements per byte); the interesting number
/// is whether the portable loops sit at the same memory ceiling as the
/// dispatched kernel.
fn bench_gf2_bits() {
    println!("bit-packed GF(2) — packed vs byte-per-element control:");
    // L1-resident, L2-resident, and DRAM-resident buffer sizes.
    for &bytes in &[4 * 1024usize, 256 * 1024, 8 * 1024 * 1024, 64 * 1024 * 1024] {
        let bits = 8 * bytes;
        let human = if bytes >= 1024 * 1024 {
            format!("{} MiB", bytes / (1024 * 1024))
        } else {
            format!("{} KiB", bytes / 1024)
        };
        let a = noise(bytes, 0x2100 + bytes as u64);
        let b = noise(bytes, 0x2200 + bytes as u64);
        let mut dst = noise(bytes, 0x2300 + bytes as u64);
        println!("  {human} buffers ({bits} elements):");
        bench("  bits::xor_assign              packed", bytes, || {
            fgf::bits::xor_assign(black_box(&mut dst), black_box(&a));
        });
        bench("  bits::and_into         packed", bytes, || {
            fgf::bits::and_into(black_box(&mut dst), black_box(&a), black_box(&b));
        });
        bench("  bits::weight           packed", bytes, || {
            black_box(fgf::bits::weight(black_box(&a), bits));
        });
        bench("  bits::dot_product       packed", bytes, || {
            black_box(fgf::bits::dot_product(black_box(&a), black_box(&b), bits));
        });
        let from = bits / 4;
        let to = bits - bits / 8;
        bench("  bits::xor_range  5/8    packed", bytes, || {
            fgf::bits::xor_range(black_box(&mut dst), black_box(&a), bits, from, to);
        });
        let plan = fgf::bits::XorRange::new(bits, from, to);
        bench("  bits::XorRange   5/8    packed", bytes, || {
            fgf::bits::xor_range_with(black_box(&mut dst), black_box(&a), black_box(&plan));
        });
        bench("  ops::add_assign  gf8 control", bytes, || {
            ops::add_assign::<Gf8B>(black_box(&mut dst), black_box(&a));
        });
    }
    println!();
}

/// The GF(2) elimination's per-call shape: short rows, sub-word ranges.
///
/// An elimination XORs an eight-bit range of one row into another, once
/// per row per pivot. Bandwidth is irrelevant at this size — the cost is
/// the per-call surface — so this section compares the checked one-shot
/// `xor_range` against the prepared `XorRange` + `xor_range_with` split
/// over the same buffer pairs. The range `[61, 69)` straddles a byte
/// boundary, the worst two-byte window.
fn bench_gf2_short_rows() {
    const ROWS: usize = 512;
    const ROW_BYTES: usize = 16; // 128 elements: the 128-column row shape
    let mut dst = noise(ROWS * ROW_BYTES, 0x5100);
    let src = noise(ROWS * ROW_BYTES, 0x5200);
    let bits = ROW_BYTES * 8;
    let pairs = ROWS / 2;
    println!("  {pairs} row pairs, 16-byte rows, bits 61..69:");
    let one_shot = bench("  bits::xor_range  8-bit row", ROW_BYTES * pairs, || {
        for i in 0..pairs {
            let off = i * 2 * ROW_BYTES;
            fgf::bits::xor_range(
                black_box(&mut dst[off..off + ROW_BYTES]),
                black_box(&src[off..off + ROW_BYTES]),
                bits,
                61,
                69,
            );
        }
    });
    let plan = fgf::bits::XorRange::new(bits, 61, 69);
    let prepared = bench("  bits::XorRange   8-bit row", ROW_BYTES * pairs, || {
        for i in 0..pairs {
            let off = i * 2 * ROW_BYTES;
            fgf::bits::xor_range_with(
                black_box(&mut dst[off..off + ROW_BYTES]),
                black_box(&src[off..off + ROW_BYTES]),
                black_box(&plan),
            );
        }
    });
    let per_call = |d: f64| d / pairs as f64;
    println!(
        "  per call: one-shot {:.2} ns, prepared {:.2} ns ({:.2}x)\n",
        per_call(one_shot),
        per_call(prepared),
        one_shot / prepared
    );
}

/// The row-interleaved addition matrix: flat `add_assign` against
/// `add_assign_rows` and a literal row-by-row loop, over the row counts and
/// row lengths where four-stream interleaving should pay (and its edges,
/// where it must not regress). Rerun under each `SIMD_BACKEND` tier for the
/// per-backend numbers; the crossover thresholds in `kernel::xor_rows` were
/// set from this matrix on the reference host.
fn bench_add_assign_rows<F: fgf::FieldKernels>(name: &str) {
    println!("add_assign_rows — {name}:");
    for &row_len in &[64usize, 1024, 64 * 1024] {
        for &rows in &[1usize, 2, 4, 8, 16, 32] {
            let len = row_len * rows;
            let src = noise(len, 0x2a00 + row_len as u64);
            let mut dst = noise(len, 0x2b00 + rows as u64);

            let flat = bench("  flat add_assign", len, || {
                ops::add_assign::<F>(black_box(&mut dst), black_box(&src));
            });
            let interleaved = bench("  add_assign_rows", len, || {
                ops::add_assign_rows::<F>(black_box(&mut dst), row_len, black_box(&src));
            });
            let looped = bench("  add_assign per row", len, || {
                let dst = black_box(&mut dst);
                let src = black_box(&src);
                for (d, s) in dst.chunks_exact_mut(row_len).zip(src.chunks_exact(row_len)) {
                    ops::add_assign::<F>(d, s);
                }
            });
            println!(
                "    rows/flat: {:.2}x, rows/loop: {:.2}x  ({rows} rows x {row_len} B)",
                flat / interleaved,
                looped / interleaved,
            );
        }
    }
    println!();
}

fn main() {
    // Family flags `--gf`, `--gdl`, `--m31`, and `--gf2` run only that
    // family's panels. Unknown flags are ignored so wrapper arguments pass
    // through unchanged, and no family flag runs every panel.
    let args: Vec<String> = std::env::args().collect();
    let has = |flag: &str| args.iter().any(|arg| arg.as_str() == flag);
    let named = has("--gf") || has("--gdl") || has("--m31") || has("--gf2");
    let gf = !named || has("--gf");
    let gdl = !named || has("--gdl");
    let m31 = !named || has("--m31");
    let gf2 = !named || has("--gf2");

    println!("fgf kernel benchmark — backend: {}", backend().name());
    if has("--network-diagnostic") {
        bench_network_alignment();
        return;
    }
    println!("  (override with SIMD_BACKEND=v3_gfni_crypto|v3|v2|neon|scalar)\n");

    if gf {
        bench_preparation_crossover();
        bench_small_row_shapes();
        bench_add_assign_rows::<Gf8B>("gf8");
        bench_add_assign_rows::<Gf16>("gf16");
        bench_large_destination();
        bench_destination_alignment();
        #[cfg(all(feature = "simd", any(target_arch = "x86", target_arch = "x86_64")))]
        bench_blocked_vs_axpy();
        bench_network_payloads();
    }
    if m31 {
        bench_add_assign_rows::<Mersenne31>("m31 (prime control, flat default)");
    }
    if gf2 {
        bench_gf2_bits();
        bench_gf2_short_rows();
    }

    // L1-resident, L2-resident, and DRAM-resident. The panel's rows all
    // belong to the gf and prime families, so a gf2-only run skips it.
    for &len in &[4 * 1024usize, 256 * 1024, 8 * 1024 * 1024] {
        if !(gf || m31 || gdl) {
            continue;
        }
        let human = if len >= 1024 * 1024 {
            format!("{} MiB", len / (1024 * 1024))
        } else {
            format!("{} KiB", len / 1024)
        };
        println!("buffer {human}:");

        let src = noise(len, 1);
        let mut dst = noise(len, 2);
        let prepared16 = ops::Coeff::<Gf16>::new(gf16::Elem::from_raw(0x53a7));
        let rhs = noise(len, 0x602);
        let mut product = vec![0; len];

        if gf {
            bench("xor                       gf8", len, || {
                ops::add_assign::<Gf8B>(black_box(&mut dst), black_box(&src));
            });
            bench("mul_add                   gf8", len, || {
                ops::mul_add::<Gf8B>(
                    black_box(&mut dst),
                    gf8b::Elem::from_raw(0x53),
                    black_box(&src),
                );
            });
            bench("mul_add                  gf16", len, || {
                ops::mul_add::<Gf16>(
                    black_box(&mut dst),
                    gf16::Elem::from_raw(0x53a7),
                    black_box(&src),
                );
            });
            bench("mul_add prepared         gf16", len, || {
                ops::mul_add_with::<Gf16>(black_box(&mut dst), &prepared16, black_box(&src));
            });
            bench("mul_into                  gf8", len, || {
                ops::mul_into::<Gf8B>(
                    black_box(&mut product),
                    gf8b::Elem::from_raw(0x53),
                    black_box(&src),
                );
            });
            bench("mul_into                 gf16", len, || {
                ops::mul_into::<Gf16>(
                    black_box(&mut product),
                    gf16::Elem::from_raw(0x53a7),
                    black_box(&src),
                );
            });
            bench("mul_assign                gf8", len, || {
                ops::mul_assign::<Gf8B>(black_box(&mut dst), gf8b::Elem::from_raw(0x53));
            });
            bench("elementwise                gf8", len, || {
                ops::mul_elementwise::<Gf8B>(
                    black_box(&mut product),
                    black_box(&src),
                    black_box(&rhs),
                );
            });
            bench("elementwise               gf16", len, || {
                ops::mul_elementwise::<Gf16>(
                    black_box(&mut product),
                    black_box(&src),
                    black_box(&rhs),
                );
            });
            bench("elementwise_assign        gf8", len, || {
                ops::mul_elementwise_assign::<Gf8B>(black_box(&mut product), black_box(&src));
            });
            bench("elementwise_assign       gf16", len, || {
                ops::mul_elementwise_assign::<Gf16>(black_box(&mut product), black_box(&src));
            });
            bench("mul_assign               gf16", len, || {
                ops::mul_assign::<Gf16>(black_box(&mut dst), gf16::Elem::from_raw(0x53a7));
            });
            bench("add_assign                gf16", len, || {
                ops::add_assign::<Gf16>(black_box(&mut dst), black_box(&src));
            });
            bench("add_assign_scalar          gf8", len, || {
                ops::add_assign_scalar::<Gf8B>(black_box(&mut dst), gf8b::Elem::from_raw(0x53));
            });
            bench("sub_assign_scalar          gf8", len, || {
                ops::sub_assign_scalar::<Gf8B>(black_box(&mut dst), gf8b::Elem::from_raw(0x53));
            });
            bench("add_assign_scalar         gf16", len, || {
                ops::add_assign_scalar::<Gf16>(black_box(&mut dst), gf16::Elem::from_raw(0x53a7));
            });
            bench("sub_assign_scalar         gf16", len, || {
                ops::sub_assign_scalar::<Gf16>(black_box(&mut dst), gf16::Elem::from_raw(0x53a7));
            });
        }

        // Prime-field elementwise and broadcast scalar ops. `add_assign` and
        // `elementwise` are the unchanged controls for the new assign forms.
        if m31 {
            bench("add_assign                   m31", len, || {
                ops::add_assign::<Mersenne31>(black_box(&mut dst), black_box(&src));
            });
            bench("add_assign_scalar            m31", len, || {
                ops::add_assign_scalar::<Mersenne31>(
                    black_box(&mut dst),
                    mersenne31::Elem::from_raw(0x1234_5678),
                );
            });
            bench("sub_assign_scalar            m31", len, || {
                ops::sub_assign_scalar::<Mersenne31>(
                    black_box(&mut dst),
                    mersenne31::Elem::from_raw(0x1234_5678),
                );
            });
            bench("elementwise                  m31", len, || {
                ops::mul_elementwise::<Mersenne31>(
                    black_box(&mut product),
                    black_box(&src),
                    black_box(&rhs),
                );
            });
            bench("elementwise_assign           m31", len, || {
                ops::mul_elementwise_assign::<Mersenne31>(black_box(&mut product), black_box(&src));
            });
        }
        if gdl {
            bench("add_assign                   gld", len, || {
                ops::add_assign::<Goldilocks>(black_box(&mut dst), black_box(&src));
            });
            bench("add_assign_scalar            gld", len, || {
                ops::add_assign_scalar::<Goldilocks>(
                    black_box(&mut dst),
                    goldilocks::Elem::from_raw(0x1234_5678_9abc_def0),
                );
            });
            bench("sub_assign_scalar            gld", len, || {
                ops::sub_assign_scalar::<Goldilocks>(
                    black_box(&mut dst),
                    goldilocks::Elem::from_raw(0x1234_5678_9abc_def0),
                );
            });
            bench("elementwise                  gld", len, || {
                ops::mul_elementwise::<Goldilocks>(
                    black_box(&mut product),
                    black_box(&src),
                    black_box(&rhs),
                );
            });
            bench("elementwise_assign           gld", len, || {
                ops::mul_elementwise_assign::<Goldilocks>(black_box(&mut product), black_box(&src));
            });
        }
        if m31 {
            bench("add_assign                  qm31", len, || {
                ops::add_assign::<QuadMersenne31>(black_box(&mut dst), black_box(&src));
            });
            bench("add_assign_scalar           qm31", len, || {
                ops::add_assign_scalar::<QuadMersenne31>(
                    black_box(&mut dst),
                    quad_mersenne31::Elem::from_raw(0x1234_5678, 0x0987_6543),
                );
            });
            bench("sub_assign_scalar           qm31", len, || {
                ops::sub_assign_scalar::<QuadMersenne31>(
                    black_box(&mut dst),
                    quad_mersenne31::Elem::from_raw(0x1234_5678, 0x0987_6543),
                );
            });
            bench("elementwise                 qm31", len, || {
                ops::mul_elementwise::<QuadMersenne31>(
                    black_box(&mut product),
                    black_box(&src),
                    black_box(&rhs),
                );
            });
            bench("elementwise_assign          qm31", len, || {
                ops::mul_elementwise_assign::<QuadMersenne31>(
                    black_box(&mut product),
                    black_box(&src),
                );
            });
        }
        println!();
    }

    // The remaining panels cover the binary tower fields only.
    if !gf {
        return;
    }

    // Multi-row shapes. Geometry is a realistic erasure code: k data rows
    // folded into m parity rows, each row a 64 KiB symbol.
    let row_len = 64 * 1024;
    for &nrows in &[2usize, 4, 8, 16] {
        println!("scatter/matrix — {nrows} rows x {} KiB:", row_len / 1024);

        let src = noise(row_len, 3);
        let mut rows = noise(row_len * nrows, 4);
        let coeffs8: Vec<_> = (0..nrows)
            .map(|j| gf8b::Elem::from_raw((j as u8).wrapping_mul(37).wrapping_add(2)))
            .collect();
        let coeffs16: Vec<_> = (0..nrows)
            .map(|j| gf16::Elem::from_raw((j as u16).wrapping_mul(9871).wrapping_add(2)))
            .collect();
        let scatter_coeffs8 = ops::CoeffVec::<Gf8B>::new(&coeffs8);
        let scatter_coeffs16 = ops::CoeffVec::<Gf16>::new(&coeffs16);

        let traffic = row_len * nrows;
        bench("scatter                   gf8", traffic, || {
            ops::mul_add_scatter::<Gf8B>(black_box(&mut rows), row_len, &coeffs8, black_box(&src));
        });
        bench("scatter                  gf16", traffic, || {
            ops::mul_add_scatter::<Gf16>(black_box(&mut rows), row_len, &coeffs16, black_box(&src));
        });
        bench("scatter prepared          gf8", traffic, || {
            ops::mul_add_scatter_with::<Gf8B>(
                black_box(&mut rows),
                row_len,
                scatter_coeffs8.as_ref(),
                black_box(&src),
            );
        });
        bench("scatter prepared         gf16", traffic, || {
            ops::mul_add_scatter_with::<Gf16>(
                black_box(&mut rows),
                row_len,
                scatter_coeffs16.as_ref(),
                black_box(&src),
            );
        });
        bench("scatter (unblocked)       gf8", traffic, || {
            for (row, &coeff) in rows.chunks_exact_mut(row_len).zip(&coeffs8) {
                ops::mul_add::<Gf8B>(black_box(row), coeff, black_box(&src));
            }
        });
        bench("scatter (unblocked)      gf16", traffic, || {
            for (row, &coeff) in rows.chunks_exact_mut(row_len).zip(&coeffs16) {
                ops::mul_add::<Gf16>(black_box(row), coeff, black_box(&src));
            }
        });

        // The blocked-vs-unblocked comparison: 8 sources into `nrows` rows,
        // as one matrix call versus eight scatter calls. Same arithmetic,
        // different destination memory traffic.
        let sources: Vec<Vec<u8>> = (0..16).map(|t| noise(row_len, 100 + t as u64)).collect();
        let coeff_sets: Vec<Vec<gf8b::Elem>> = (0..8)
            .map(|t| {
                (0..nrows)
                    .map(|j| gf8b::Elem::from_raw(((t * 31 + j * 17) as u8).wrapping_add(1)))
                    .collect()
            })
            .collect();
        let coeff_sets16: Vec<Vec<gf16::Elem>> = (0..8)
            .map(|t| {
                (0..nrows)
                    .map(|j| gf16::Elem::from_raw(((t * 7919 + j * 613) as u16).wrapping_add(1)))
                    .collect()
            })
            .collect();
        let terms: Vec<(&[gf8b::Elem], &[u8])> = coeff_sets
            .iter()
            .zip(&sources)
            .map(|(c, s)| (c.as_slice(), s.as_slice()))
            .collect();
        let terms16: Vec<(&[gf16::Elem], &[u8])> = coeff_sets16
            .iter()
            .zip(&sources)
            .map(|(c, s)| (c.as_slice(), s.as_slice()))
            .collect();
        let matrix_coeffs8: Vec<_> = coeff_sets.iter().flatten().copied().collect();
        let matrix_coeffs16: Vec<_> = coeff_sets16.iter().flatten().copied().collect();
        let matrix8 = ops::CoeffMatrix::<Gf8B>::from_source_major(8, nrows, &matrix_coeffs8);
        let matrix16 = ops::CoeffMatrix::<Gf16>::from_source_major(8, nrows, &matrix_coeffs16);
        let matrix_srcs: Vec<&[u8]> = sources.iter().take(8).map(Vec::as_slice).collect();

        let traffic = row_len * nrows * 8;
        bench("matrix (selected)          gf8", traffic, || {
            ops::mul_add_matrix::<Gf8B>(black_box(&mut rows), row_len, nrows, &terms);
        });
        bench("matrix (unblocked AXPY)   gf8", traffic, || {
            for &(coeffs, src) in &terms {
                for (row, &coeff) in rows.chunks_exact_mut(row_len).zip(coeffs) {
                    ops::mul_add::<Gf8B>(black_box(row), coeff, src);
                }
            }
        });
        bench("matrix (selected)         gf16", traffic, || {
            ops::mul_add_matrix::<Gf16>(black_box(&mut rows), row_len, nrows, &terms16);
        });
        bench("matrix prepared           gf8", traffic, || {
            ops::mul_add_matrix_with::<Gf8B>(
                black_box(&mut rows),
                row_len,
                &matrix8,
                black_box(&matrix_srcs),
            );
        });
        bench("matrix prepared          gf16", traffic, || {
            ops::mul_add_matrix_with::<Gf16>(
                black_box(&mut rows),
                row_len,
                &matrix16,
                black_box(&matrix_srcs),
            );
        });
        bench("matrix (unblocked AXPY)  gf16", traffic, || {
            for &(coeffs, src) in &terms16 {
                for (row, &coeff) in rows.chunks_exact_mut(row_len).zip(coeffs) {
                    ops::mul_add::<Gf16>(black_box(row), coeff, src);
                }
            }
        });

        let gather_srcs: Vec<&[u8]> = sources.iter().take(nrows).map(Vec::as_slice).collect();
        let mut gathered = noise(row_len, 5);
        let gather_coeffs8 = ops::CoeffVec::<Gf8B>::new(&coeffs8);
        let gather_coeffs16 = ops::CoeffVec::<Gf16>::new(&coeffs16);
        let gather_traffic = row_len * nrows;
        bench("gather (selected)          gf8", gather_traffic, || {
            ops::mul_add_gather::<Gf8B>(
                black_box(&mut gathered),
                &coeffs8,
                black_box(&gather_srcs),
            );
        });
        bench("gather (unblocked)        gf8", gather_traffic, || {
            for (&coeff, &source) in coeffs8.iter().zip(&gather_srcs) {
                ops::mul_add::<Gf8B>(black_box(&mut gathered), coeff, black_box(source));
            }
        });
        bench("gather (selected)         gf16", gather_traffic, || {
            ops::mul_add_gather::<Gf16>(
                black_box(&mut gathered),
                &coeffs16,
                black_box(&gather_srcs),
            );
        });
        bench("gather prepared           gf8", gather_traffic, || {
            ops::mul_add_gather_with::<Gf8B>(
                black_box(&mut gathered),
                gather_coeffs8.as_ref(),
                black_box(&gather_srcs),
            );
        });
        bench("gather prepared          gf16", gather_traffic, || {
            ops::mul_add_gather_with::<Gf16>(
                black_box(&mut gathered),
                gather_coeffs16.as_ref(),
                black_box(&gather_srcs),
            );
        });
        bench("gather (unblocked)       gf16", gather_traffic, || {
            for (&coeff, &source) in coeffs16.iter().zip(&gather_srcs) {
                ops::mul_add::<Gf16>(black_box(&mut gathered), coeff, black_box(source));
            }
        });
        println!();
    }

    // Same byte volume means the same information volume: one GF(2^32)
    // symbol occupies exactly two GF(2^16) symbols. This makes the cost of
    // choosing a wider field directly visible rather than hiding it behind a
    // symbol-count comparison.
    let tier3_len = 256 * 1024;
    let tier3_src = noise(tier3_len, 0x900);
    let mut tier3_dst = noise(tier3_len, 0x901);
    println!("Tier 3 field cost — {tier3_len} bytes:");
    bench("mul_add polynomial tower     gf16", tier3_len, || {
        ops::mul_add::<Gf16>(
            black_box(&mut tier3_dst),
            gf16::Elem::from_raw(0x53a7),
            black_box(&tier3_src),
        );
    });
    bench("mul_add polynomial tower     gf32", tier3_len, || {
        ops::mul_add::<Gf32>(
            black_box(&mut tier3_dst),
            gf32::Elem::from_raw(0xdead_beef),
            black_box(&tier3_src),
        );
    });
    bench("mul_add polynomial tower     gf64", tier3_len, || {
        ops::mul_add::<Gf64>(
            black_box(&mut tier3_dst),
            gf64::Elem::from_raw(0x0123_4567_89ab_cdef),
            black_box(&tier3_src),
        );
    });
    bench("mul_add canonical Fan-Paar  fp16", tier3_len, || {
        ops::mul_add::<FanPaar16>(
            black_box(&mut tier3_dst),
            fan_paar::fp16::Elem::from_raw(0xe2de),
            black_box(&tier3_src),
        );
    });
    bench("mul_add canonical Fan-Paar  fp32", tier3_len, || {
        ops::mul_add::<FanPaar32>(
            black_box(&mut tier3_dst),
            fan_paar::fp32::Elem::from_raw(0x03e2_1cea),
            black_box(&tier3_src),
        );
    });
    bench("mul_add canonical Fan-Paar  fp64", tier3_len, || {
        ops::mul_add::<FanPaar64>(
            black_box(&mut tier3_dst),
            fan_paar::fp64::Elem::from_raw(0x070f_870d_cd9c_1d88),
            black_box(&tier3_src),
        );
    });
}