//! The narrow convolution kernels: 16-bit operands, 32-bit partial sums, 64-bit sums
//! (the parent module's documentation, *The narrow kernels*, states the rule).
//!
//! A convolution of the operator set forms each output's exact sum `Σ w·x`. Where every
//! weight has `|w| ≤ 2^15 − 1` and every input `|x| < 2^30`, this module forms the same
//! sum from 16-bit operands:
//!
//! - an input with `|x| ≤ 2^15 − 1` is one operand (`Limb::Whole`); otherwise each input
//!   is split into `hi = ⌊x / 2^15⌋` and `lo = x − hi·2^15 ∈ [0, 2^15)`, both 16-bit, and
//!   `Σ w·x = 2^15·Σ w·hi + Σ w·lo` (`Limb::Hi`, `Limb::Lo`);
//! - each output's products are summed in a 32-bit lane over a block of at most
//!   `kblk` taps, where `kblk` is chosen from the operands' largest magnitudes so that
//!   no partial sum of a block can leave `i32`, and each block's sum is then added,
//!   widened and shifted by its limb's `2^15` or `2^0`, into a 64-bit sum;
//! - the plan admits a convolution only where every such 64-bit sum, at every step,
//!   stays within `i64` too ([`plan`]).
//!
//! No value at any step leaves the type it is held in, so every step is exact and the
//! result is the exact integer `Σ w·x`: the same integer the interpreter forms in
//! `i128`, in whatever order the terms are added. On x86_64 two input channels share one
//! 32-bit word (`vpmaddwd` and `vpdpwssd` multiply 16-bit pairs and add both products
//! into one 32-bit lane); on aarch64 each channel has its own 16-bit lane (`smlal`). Rows
//! whose length is not a multiple of the vector width finish in scalar 64-bit arithmetic.

use super::{par_for_each, weight_at, ConvGeom, Isa, Layout, Tile};

/// Taps one pass over a row takes at most (fewer where the plan's block is shorter): the
/// kernels sweep the row's elements once per pass, so its weights and the input lines its
/// taps read stay in the first-level cache across the sweep.
pub(super) const TAPS_PER_PASS: usize = 128;

/// The largest operand magnitude a 16-bit lane takes here. With `|w| ≤ 2^15 − 1`, a sum
/// of two products of 16-bit operands stays below `2^31` (only `(−2^15)·(−2^15)` twice
/// reaches it).
const I16: u64 = i16::MAX as u64;

/// How an input value becomes a 16-bit operand.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
enum Limb {
    /// `x` itself, `|x| ≤ 2^15 − 1`.
    Whole,
    /// `x mod 2^15`, in `[0, 2^15)`.
    Lo,
    /// `⌊x / 2^15⌋`, in `[−2^15, 2^15)` for `|x| < 2^30`.
    Hi,
}

impl Limb {
    #[inline(always)]
    fn of(self, v: i64) -> i16 {
        match self {
            Limb::Whole => v as i16,
            Limb::Lo => (v & 0x7fff) as i16,
            Limb::Hi => (v >> 15) as i16,
        }
    }
    /// The power of two its sum is scaled by.
    fn shift(self) -> u32 {
        if self == Limb::Hi { 15 } else { 0 }
    }
}

/// The kernel family an [`Isa`] runs here, or none.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
enum Kernel {
    #[cfg(target_arch = "x86_64")]
    X86(super::x86::Kind),
    #[cfg(target_arch = "aarch64")]
    Neon,
}

impl Kernel {
    /// `isa` must already be [`Isa::effective`], so a set named here is one the processor runs.
    fn of(isa: Isa) -> Option<Kernel> {
        match isa {
            #[cfg(target_arch = "x86_64")]
            Isa::Avx2 => Some(Kernel::X86(super::x86::Kind::Avx2)),
            #[cfg(target_arch = "x86_64")]
            Isa::Avx2Vnni => Some(Kernel::X86(super::x86::Kind::Avx2Vnni)),
            #[cfg(target_arch = "x86_64")]
            Isa::Avx512 => Some(Kernel::X86(super::x86::Kind::Avx512)),
            #[cfg(target_arch = "x86_64")]
            Isa::Avx512Vnni => Some(Kernel::X86(super::x86::Kind::Avx512Vnni)),
            #[cfg(target_arch = "aarch64")]
            Isa::Neon => Some(Kernel::Neon),
            _ => None,
        }
    }
    /// Whether two channels share a 32-bit word.
    fn pairs(self) -> bool {
        match self {
            #[cfg(target_arch = "x86_64")]
            Kernel::X86(_) => true,
            #[cfg(target_arch = "aarch64")]
            Kernel::Neon => false,
        }
    }
    /// Output channels one tile holds in registers.
    fn block(self) -> usize {
        match self {
            #[cfg(target_arch = "x86_64")]
            Kernel::X86(k) => k.block(),
            #[cfg(target_arch = "aarch64")]
            Kernel::Neon => 8,
        }
    }
}

/// The limbs a convolution runs, each with the taps one 32-bit partial sum may take.
#[derive(Clone, Debug, PartialEq, Eq)]
struct Plan {
    limbs: Vec<(Limb, usize)>,
}

/// The plan for `g` with inputs `|x| ≤ xm` and weights `|w| ≤ wm`, or `None` where the
/// narrow kernels do not apply. `pairs`: a tap adds two products (two channels).
///
/// For each limb with operand bound `l`, one product is at most `p = wm·l` in magnitude,
/// so a block of `k` taps (`k` or `2k` products) has every partial sum within
/// `k·p` (`2k·p`), which `kblk` keeps within `i32::MAX`. Every 64-bit sum is a sum of
/// some of an output's at most `terms` products of each limb, scaled by its limb's power
/// of two, so it stays within `Σ_limbs terms·p·2^shift`, which must be within `i64::MAX`.
fn plan(g: &ConvGeom, xm: u64, wm: u64, pairs: bool) -> Option<Plan> {
    if wm > I16 || xm >= 1 << 30 {
        return None;
    }
    let terms = g.terms_per_output() as u128;
    let per_tap: u128 = if pairs { 2 } else { 1 };
    let limbs: &[(Limb, u64)] = if xm <= I16 { &[(Limb::Whole, xm)] } else { &[(Limb::Hi, (xm >> 15) + 1), (Limb::Lo, I16)] };
    let mut total: u128 = 0;
    let mut out = Vec::with_capacity(limbs.len());
    for &(l, lmax) in limbs {
        let p = wm as u128 * lmax as u128;
        let kblk = if p == 0 { usize::MAX } else { usize::try_from(i32::MAX as u128 / (per_tap * p)).unwrap_or(usize::MAX) };
        if kblk == 0 {
            return None;
        }
        total = total.checked_add(terms.checked_mul(p)? << l.shift())?;
        out.push((l, kblk));
    }
    (total <= i64::MAX as u128).then_some(Plan { limbs: out })
}

/// One limb's padded input in the kernel's form.
enum Planes {
    /// Two channels per 32-bit word: the first in the low half, both as 16-bit two's
    /// complement. Plane `grp·pg + q` holds channels `grp·cin_g + 2q` and `+ 1` (zero
    /// where the group has no such channel).
    #[cfg(target_arch = "x86_64")]
    Pairs(Vec<u32>),
    /// One channel per 16-bit lane, plane `ic` for channel `ic`.
    #[cfg(target_arch = "aarch64")]
    Single(Vec<i16>),
}

/// The pair planes of [`Planes::Pairs`] from the channel planes `xp` (`cin` planes of
/// `chan` values each).
#[cfg(target_arch = "x86_64")]
fn interleave(g: &ConvGeom, chan: usize, xp: &[i16], threads: usize) -> Vec<u32> {
    let cin_g = g.cin / g.groups;
    let pg = cin_g.div_ceil(2);
    let mut out = vec![0u32; g.groups * pg * chan];
    let threads = threads.min((out.len() / (1 << 16)).max(1));
    par_for_each(out.chunks_mut(chan).enumerate().collect(), threads, |(pl, dst): (usize, &mut [u32])| {
        let (grp, q) = (pl / pg, pl % pg);
        let a = grp * cin_g + 2 * q;
        let lo = &xp[a * chan..(a + 1) * chan];
        if 2 * q + 1 < cin_g {
            let hi = &xp[(a + 1) * chan..(a + 2) * chan];
            for ((d, &l), &h) in dst.iter_mut().zip(lo).zip(hi) {
                *d = (l as u16 as u32) | ((h as u16 as u32) << 16);
            }
        } else {
            for (d, &l) in dst.iter_mut().zip(lo) {
                *d = l as u16 as u32;
            }
        }
    });
    out
}

/// One tap table: each tap's weights for the tile's channels (`STRIDE` lanes per tap,
/// [`TapW`]), and its input offset.
struct Taps {
    w: Vec<TapW>,
    offs: Vec<usize>,
    /// The largest offset (0 for an empty table).
    max_off: usize,
}

impl Taps {
    fn new() -> Self {
        Taps { w: Vec::new(), offs: Vec::new(), max_off: 0 }
    }
    fn push(&mut self, w: &[TapW], off: usize) {
        self.w.extend_from_slice(w);
        self.offs.push(off);
        self.max_off = self.max_off.max(off);
    }
    fn is_empty(&self) -> bool {
        self.offs.is_empty()
    }
}

/// The narrow convolution of `g`, or `None` where it does not apply (an [`Isa`] with no
/// narrow kernels, or operands past [`plan`]'s bounds); the caller then runs the wide
/// kernels. `epi(oc, sum)` finishes each output from its exact sum.
#[allow(clippy::too_many_arguments)]
pub(super) fn conv<E>(g: &ConvGeom, x: &[i64], wt: &[i64], xm: u64, wm: u64, threads: usize, isa: Isa, epi: &E) -> Option<Vec<i64>>
where
    E: Fn(usize, i64) -> i64 + Sync,
{
    let k = Kernel::of(isa)?;
    let plan = plan(g, xm, wm, k.pairs())?;
    super::note_narrow();
    let lay = Layout::new(g).expect("ConvGeom::sized checked the layout");
    let chan = lay.hp * lay.rw;
    let planes: Vec<(Planes, usize, u32)> = plan
        .limbs
        .iter()
        .map(|&(l, kblk)| {
            let xp = lay.pad_with(g, x, threads, |v| l.of(v));
            #[cfg(target_arch = "x86_64")]
            let p = Planes::Pairs(interleave(g, chan, &xp, threads));
            #[cfg(target_arch = "aarch64")]
            let p = {
                let _ = chan;
                Planes::Single(xp)
            };
            (p, kblk, l.shift())
        })
        .collect();

    let plane = g.ho * g.wo;
    let mut out = vec![0i64; g.cout * plane];
    let block = k.block();
    let cout_g = g.cout / g.groups;
    let mut blocks = Vec::new();
    for grp in 0..g.groups {
        let mut oc = grp * cout_g;
        while oc < (grp + 1) * cout_g {
            let nb = block.min((grp + 1) * cout_g - oc);
            blocks.push((oc, nb));
            oc += nb;
        }
    }
    let nbands = if threads > 1 { threads.saturating_mul(4).div_ceil(blocks.len()).clamp(1, g.ho) } else { 1 };
    let rows = g.ho.div_ceil(nbands);
    let mut bands: Vec<std::slice::ChunksMut<'_, i64>> = out.chunks_mut(plane).map(|p| p.chunks_mut(rows * g.wo)).collect();
    let mut tiles = Vec::new();
    for &(oc0, nb) in &blocks {
        for b in 0..g.ho.div_ceil(rows) {
            let slices: Vec<&mut [i64]> = (oc0..oc0 + nb).map(|oc| bands[oc].next().unwrap()).collect();
            tiles.push(Tile { oc0, y0: b * rows, slices });
        }
    }
    let ctx = Ctx { k, g, lay: &lay, planes: &planes, wt };
    par_for_each(tiles, threads, |t| match t.slices.len() {
        1 => ctx.tile::<1, E>(epi, t),
        2 => ctx.tile::<2, E>(epi, t),
        3 => ctx.tile::<3, E>(epi, t),
        4 => ctx.tile::<4, E>(epi, t),
        5 => ctx.tile::<5, E>(epi, t),
        6 => ctx.tile::<6, E>(epi, t),
        7 => ctx.tile::<7, E>(epi, t),
        _ => ctx.tile::<8, E>(epi, t),
    });
    Some(out)
}

/// What every tile of one convolution reads.
struct Ctx<'a> {
    k: Kernel,
    g: &'a ConvGeom,
    lay: &'a Layout,
    planes: &'a [(Planes, usize, u32)],
    wt: &'a [i64],
}

/// One weight lane of a tap: a pair of 16-bit weights (x86_64), or one (aarch64).
#[cfg(target_arch = "x86_64")]
type TapW = i32;
#[cfg(target_arch = "aarch64")]
type TapW = i16;

impl Ctx<'_> {
    /// The weights of channels `oc0 .. oc0 + NB` for input channel `ic` (or, with pairs,
    /// channels `ic` and `ic + 1` of the group, `ic + 1` only where `second`) at kernel
    /// tap `(ky, kx)`; `None` where every one is zero, a tap the sums may omit.
    ///
    /// x86_64 lays a tap out as `NB` words, each a pair of 16-bit weights (channel `ic`
    /// in the low half); aarch64 as 8 16-bit weights, the lanes past `NB` zero, so that
    /// one vector load reads them.
    #[inline(always)]
    #[allow(clippy::too_many_arguments)]
    fn push_tap<const NB: usize>(&self, taps: &mut Taps, off: usize, oc0: usize, icl: usize, second: bool, ky: usize, kx: usize) {
        let (g, wt) = (self.g, self.wt);
        #[cfg(target_arch = "x86_64")]
        let w: [i32; NB] = std::array::from_fn(|j| {
            let a = weight_at(g, wt, oc0 + j, icl, ky, kx) as i16 as u16 as u32;
            let b = if second { weight_at(g, wt, oc0 + j, icl + 1, ky, kx) as i16 as u16 as u32 } else { 0 };
            (a | (b << 16)) as i32
        });
        #[cfg(target_arch = "aarch64")]
        let w: [i16; 8] = {
            let _ = second;
            std::array::from_fn(|j| if j < NB { weight_at(g, wt, oc0 + j, icl, ky, kx) as i16 } else { 0 })
        };
        if w.iter().any(|&v| v != 0) {
            taps.push(&w, off);
        }
    }

    /// Input channels per kernel lane step: two with pairs, else one.
    fn step(&self) -> usize {
        if self.k.pairs() { 2 } else { 1 }
    }

    fn tile<const NB: usize, E: Fn(usize, i64) -> i64 + Sync>(&self, epi: &E, mut t: Tile<'_>) {
        let (g, lay) = (self.g, self.lay);
        let rows = t.slices[0].len() / g.wo;
        let chan = lay.hp * lay.rw;
        let step = self.step();

        if !g.transposed {
            let cin_g = g.cin / g.groups;
            let grp = t.oc0 / (g.cout / g.groups);
            // the group's first plane: pair planes with pairs, channel planes without
            let p0 = if step == 2 { grp * cin_g.div_ceil(2) } else { grp * cin_g };
            let mut taps = Taps::new();
            for (q, icl) in (0..cin_g).step_by(step).enumerate() {
                for ky in 0..g.kh {
                    for kx in 0..g.kw {
                        let off = (p0 + q) * chan + ky * lay.rw + (kx % g.sw) * lay.pwl + kx / g.sw;
                        self.push_tap::<NB>(&mut taps, off, t.oc0, icl, icl + 1 < cin_g, ky, kx);
                    }
                }
            }
            let mut acc = vec![0i64; NB * g.wo];
            for r in 0..rows {
                let oy = t.y0 + r;
                self.sums::<NB>(&taps, oy * g.sh * lay.rw, g.wo, &mut acc);
                for (j, s) in t.slices.iter_mut().enumerate() {
                    let oc = t.oc0 + j;
                    for (d, &a) in s[r * g.wo..(r + 1) * g.wo].iter_mut().zip(&acc[j * g.wo..(j + 1) * g.wo]) {
                        *d = epi(oc, a);
                    }
                }
            }
            return;
        }

        // Transposed: the phase tables of the wide kernels' `tile`, over pairs of input
        // channels where the kernel takes pairs.
        let (sh, sw) = (g.sh, g.sw);
        let (nsig, nrho) = (rows.min(sh), g.wo.min(sw));
        let mut tables: Vec<(Taps, usize)> = Vec::with_capacity(nsig * nrho);
        for i in 0..nsig {
            let sigma = (t.y0 + i + g.ph) % sh;
            let nqy = g.kh.saturating_sub(sigma).div_ceil(sh);
            for b in 0..nrho {
                let rho = (b + g.pw) % sw;
                let nqx = g.kw.saturating_sub(rho).div_ceil(sw);
                let n = (g.wo - b).div_ceil(sw);
                let u0 = (b + g.pw) / sw;
                let mut taps = Taps::new();
                if nqy > 0 && nqx > 0 {
                    for (q, ic) in (0..g.cin).step_by(step).enumerate() {
                        for qy in (0..nqy).rev() {
                            let ky = sigma + qy * sh;
                            for qx in (0..nqx).rev() {
                                let kx = rho + qx * sw;
                                let off = q * chan + (lay.pt - qy) * lay.rw + (lay.pl - qx) + u0;
                                self.push_tap::<NB>(&mut taps, off, t.oc0, ic, ic + 1 < g.cin, ky, kx);
                            }
                        }
                    }
                }
                tables.push((taps, n));
            }
        }
        let mut acc = vec![0i64; NB * g.wo.div_ceil(sw)];
        for r in 0..rows {
            let v = (t.y0 + r + g.ph) / sh;
            let i = r % nsig;
            for (b, (taps, n)) in tables[i * nrho..(i + 1) * nrho].iter().enumerate() {
                let n = *n;
                self.sums::<NB>(taps, v * lay.rw, n, &mut acc[..NB * n]);
                for (j, s) in t.slices.iter_mut().enumerate() {
                    let oc = t.oc0 + j;
                    let row = &mut s[r * g.wo..(r + 1) * g.wo];
                    for (e, &a) in acc[j * n..(j + 1) * n].iter().enumerate() {
                        row[b + e * sw] = epi(oc, a);
                    }
                }
            }
        }
    }

    /// `acc[j·n + e] = Σ_t w_t[j] · x[offs_t + base + e]` over every limb, for the `n`
    /// elements of one row and the tile's `NB` channels.
    fn sums<const NB: usize>(&self, taps: &Taps, base: usize, n: usize, acc: &mut [i64]) {
        let acc = &mut acc[..NB * n];
        acc.fill(0);
        if n == 0 || taps.is_empty() {
            return;
        }
        for (p, kblk, sh) in self.planes {
            match (self.k, p) {
                #[cfg(target_arch = "x86_64")]
                (Kernel::X86(kind), Planes::Pairs(x)) => {
                    // Every element a kernel reads, `x[offs_t + base + e]` for `e < n`, lies
                    // within `x`: the raw loads below rely on this check.
                    assert!(taps.max_off + base + n <= x.len(), "a narrow row reads within its padded input");
                    // SAFETY: the processor runs `kind` (`Kernel::of` takes only an
                    // `Isa::effective` set), and the assertion bounds every read.
                    unsafe { super::x86::sums::<NB>(kind, &taps.w, &taps.offs, x, base, n, *kblk, *sh, acc) }
                }
                #[cfg(target_arch = "aarch64")]
                (Kernel::Neon, Planes::Single(x)) => {
                    assert!(taps.max_off + base + n <= x.len(), "a narrow row reads within its padded input");
                    // SAFETY: NEON is part of the aarch64 baseline, and the assertion bounds
                    // every read.
                    unsafe { neon::sums::<NB>(&taps.w, &taps.offs, x, base, n, *kblk, *sh, acc) }
                }
            }
        }
    }
}

/// Elements `e0 .. n` of [`Ctx::sums`] for one limb in scalar 64-bit arithmetic, with
/// pairs: `acc[j·n + e] += (Σ_t lo(w)·lo(x) + hi(w)·hi(x)) · 2^sh`.
#[cfg(target_arch = "x86_64")]
#[allow(clippy::too_many_arguments)]
pub(super) fn tail_pairs<const NB: usize>(ws: &[i32], offs: &[usize], x: &[u32], base: usize, e0: usize, n: usize, sh: u32, acc: &mut [i64]) {
    for e in e0..n {
        let mut s = [0i64; NB];
        for (w, &o) in ws.as_chunks::<NB>().0.iter().zip(offs) {
            let xv = x[o + base + e];
            let (xl, xh) = (xv as u16 as i16 as i64, (xv >> 16) as u16 as i16 as i64);
            for j in 0..NB {
                let wv = w[j] as u32;
                s[j] += (wv as u16 as i16 as i64) * xl + ((wv >> 16) as u16 as i16 as i64) * xh;
            }
        }
        for j in 0..NB {
            acc[j * n + e] += s[j] << sh;
        }
    }
}

/// [`tail_pairs`] with one channel per lane.
#[cfg(target_arch = "aarch64")]
#[allow(clippy::too_many_arguments)]
fn tail_single<const NB: usize>(ws: &[i16], offs: &[usize], x: &[i16], base: usize, e0: usize, n: usize, sh: u32, acc: &mut [i64]) {
    for e in e0..n {
        let mut s = [0i64; NB];
        for (w, &o) in ws.as_chunks::<8>().0.iter().zip(offs) {
            let xv = x[o + base + e] as i64;
            for j in 0..NB {
                s[j] += w[j] as i64 * xv;
            }
        }
        for j in 0..NB {
            acc[j * n + e] += s[j] << sh;
        }
    }
}

/// The NEON narrow kernel: `NB ≤ 8` channels × 8 adjacent elements in 32-bit lanes
/// (`smlal` / `smlal2` by weight lane), each block of at most `kblk` taps summed from
/// zero and then widened, shifted and added into the 64-bit sums.
#[cfg(target_arch = "aarch64")]
mod neon {
    use std::arch::aarch64::*;

    macro_rules! lanes {
        ($nb:expr, $a:ident, $xl:ident, $xv:ident, $wv:ident; $($j:literal)*) => {
            $(
                if $nb > $j {
                    $a[$j][0] = vmlal_laneq_s16::<$j>($a[$j][0], $xl, $wv);
                    $a[$j][1] = vmlal_high_laneq_s16::<$j>($a[$j][1], $xv, $wv);
                }
            )*
        };
    }

    /// See [`super::Ctx::sums`]. `kblk` taps' products fit a 32-bit lane ([`super::plan`]).
    ///
    /// # Safety
    /// `offs_t + base + n ≤ x.len()` for every tap `t`, and `acc.len() ≥ NB·n`.
    #[allow(clippy::too_many_arguments)]
    pub(super) unsafe fn sums<const NB: usize>(ws: &[i16], offs: &[usize], x: &[i16], base: usize, n: usize, kblk: usize, sh: u32, acc: &mut [i64]) {
        let taps = offs.len();
        let xp = x.as_ptr();
        let ap = acc.as_mut_ptr();
        // SAFETY: NEON is part of the aarch64 baseline.
        let shv = unsafe { vdupq_n_s64(sh as i64) };
        // Tap blocks outermost, so that one block's weights and the input lines its taps
        // read stay in the first-level cache while every 8-element step of the row
        // passes over them.
        let blk = kblk.min(super::TAPS_PER_PASS);
        let n8 = n - n % 8;
        let mut t0 = 0;
        while t0 < taps {
            let t1 = t0 + blk.min(taps - t0);
            let mut e = 0;
            while e < n8 {
                // SAFETY (all raw accesses below): each load reads `x[offs_t + base + e ..
                // + 8]` with `e + 8 ≤ n`, in bounds by the caller's contract; each weight
                // load reads `ws[8t .. 8t + 8]`; each store writes `acc[j·n + e .. + 8]`.
                unsafe {
                    let z = vdupq_n_s32(0);
                    let mut a = [[z; 2]; NB];
                    for t in t0..t1 {
                        let xv = vld1q_s16(xp.add(*offs.get_unchecked(t) + base + e));
                        let xl = vget_low_s16(xv);
                        let wv = vld1q_s16(ws.as_ptr().add(8 * t));
                        lanes!(NB, a, xl, xv, wv; 0 1 2 3 4 5 6 7);
                    }
                    for (j, aj) in a.iter().enumerate() {
                        let d = ap.add(j * n + e);
                        for (h, &v) in aj.iter().enumerate() {
                            let lo = vshlq_s64(vmovl_s32(vget_low_s32(v)), shv);
                            let hi = vshlq_s64(vmovl_high_s32(v), shv);
                            let p = d.add(4 * h);
                            vst1q_s64(p, vaddq_s64(vld1q_s64(p), lo));
                            vst1q_s64(p.add(2), vaddq_s64(vld1q_s64(p.add(2)), hi));
                        }
                    }
                }
                e += 8;
            }
            t0 = t1;
        }
        let e = n8;
        super::tail_single::<NB>(ws, offs, x, base, e, n, sh, acc);
    }
}

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

    #[test]
    fn a_limb_split_rebuilds_every_input() {
        for v in [-(1i64 << 30) + 1, -(1 << 15) - 1, -(1 << 15), -32767, -1, 0, 1, 32767, 32768, (1 << 30) - 1, 316_606, -316_606] {
            let (hi, lo) = (Limb::Hi.of(v) as i64, Limb::Lo.of(v) as i64);
            assert!((0..1 << 15).contains(&lo), "{v}");
            assert_eq!(hi * (1 << 15) + lo, v, "{v}");
        }
        for v in [-32767i64, -1, 0, 1, 32767] {
            assert_eq!(Limb::Whole.of(v) as i64, v);
        }
    }

    #[test]
    fn the_plan_keeps_every_partial_sum_in_its_type() {
        let g = ConvGeom::transposed(192, 128, 192, 192, 5, 5, 2, 2, 2, 2, 1, 1, 192 * 128 * 192, 192 * 192 * 25).unwrap();
        for pairs in [false, true] {
            let per = if pairs { 2u128 } else { 1 };
            for (xm, wm) in [(7622u64, 2047u64), (32767, 32767), (32768, 2047), (316_606, 2047), ((1 << 30) - 1, 32767), (0, 2047), (5, 0)] {
                let p = plan(&g, xm, wm, pairs).unwrap_or_else(|| panic!("{xm} {wm}"));
                let mut total = 0u128;
                for &(l, k) in &p.limbs {
                    let lmax = match l { Limb::Whole => xm, Limb::Hi => (xm >> 15) + 1, Limb::Lo => 32767 } as u128;
                    let prod = wm as u128 * lmax;
                    assert!(k >= 1 && (prod == 0 || per * k as u128 * prod <= i32::MAX as u128), "{xm} {wm} {l:?} {k}");
                    total += (g.terms_per_output() as u128 * prod) << l.shift();
                }
                assert!(total <= i64::MAX as u128);
                assert_eq!(p.limbs.len(), if xm <= 32767 { 1 } else { 2 });
            }
            // past either width, the wide kernels run
            assert_eq!(plan(&g, 1 << 30, 1, pairs), None);
            assert_eq!(plan(&g, 1, 32768, pairs), None);
        }
    }
}
