fgf 2.0.0

//! Precomputed per-coefficient multiplication tables.
//!
//! Two families, matching the two ways SIMD hardware can multiply in a binary
//! field:
//!
//! - [`ScaleTable`] — split-nibble shuffle tables for one GF(2^8) coefficient.
//!   `c * x` is decomposed as `c * (x & 0xf) ^ c * (x & 0xf0)`, each half a
//!   16-entry lookup that `PSHUFB`/`TBL` performs on every lane at once.
//! - [`TowerCoeff`] — the two broadcast words that turn one GF(2^16) tower
//!   multiply into two byte-wide multiplies of the source and its
//!   adjacent-swapped self. No table at all, which is what makes the GFNI and
//!   `PMULL` paths cheap.
//!
//! The GF(2^8) table bank is shared and built once; all 256 coefficients cost
//! 8 KiB and stay resident in L1/L2. GF(2^16) has 65536 coefficients, so a
//! full bank would be ~9 MiB and thrash cache. A GF(2^16) coefficient is
//! instead resolved per call into its four base-field factors (two base
//! multiplies, [`TowerCoeff::new`]) and, on shuffle backends, four table
//! copies out of the shared bank ([`TowerTables::new`]) — amortized over the
//! whole buffer, or hoisted out entirely with `Coeff`/`CoeffVec`.

use crate::field::fan_paar::{fp8, fp16};
use crate::field::{gf8b, gf8d, gf16};

/// Split-nibble multiplication tables for one GF(2^8) coefficient.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct ScaleTable {
    /// The coefficient these tables multiply by.
    pub coeff: gf8b::Elem,
    /// `lo[i] = coeff * i` for the low nibble.
    pub lo: [u8; 16],
    /// `hi[i] = coeff * (i << 4)` for the high nibble.
    pub hi: [u8; 16],
}

impl ScaleTable {
    /// Build the nibble tables for `coeff`.
    #[must_use]
    // Loop counters are bounded by the array sizes (16, 256), so every cast
    // below is exact; `const fn` rules out `try_into`.
    #[allow(clippy::cast_possible_truncation)]
    pub const fn new(coeff: gf8b::Elem) -> Self {
        let mut lo = [0u8; 16];
        let mut hi = [0u8; 16];
        let mut i = 0;
        while i < 16 {
            lo[i] = gf8b::Elem(i as u8).mul(coeff).0;
            hi[i] = gf8b::Elem((i as u8) << 4).mul(coeff).0;
            i += 1;
        }
        Self { coeff, lo, hi }
    }
}

/// Shared table bank, one entry per GF(2^8) coefficient.
///
/// A `static` rather than a `const`: a `const` array would be materialized
/// onto the stack at every use site instead of being borrowed from rodata.
static SCALE_TABLE_BANK: [ScaleTable; 256] = build_bank();

#[allow(clippy::cast_possible_truncation)]
const fn build_bank() -> [ScaleTable; 256] {
    let mut bank = [ScaleTable::new(gf8b::Elem(0)); 256];
    let mut i = 0;
    while i < 256 {
        bank[i] = ScaleTable::new(gf8b::Elem(i as u8));
        i += 1;
    }
    bank
}

/// Return the shared nibble tables for a GF(2^8) coefficient.
#[inline]
#[must_use]
pub fn scale_table(coeff: gf8b::Elem) -> &'static ScaleTable {
    &SCALE_TABLE_BANK[coeff.0 as usize]
}

/// GF(2^8)/`0x11D` nibble-table bank, one entry per coefficient.
///
/// Same layout and 8 KiB budget as the AES [`SCALE_TABLE_BANK`]; only the
/// products differ. Because the shuffle kernels read nothing but `lo`/`hi`,
/// one implementation serves both byte fields — this bank just carries the
/// `0x11D` products.
static SCALE_TABLE_BANK_8D: [ScaleTable; 256] = build_bank_8d();

#[allow(clippy::cast_possible_truncation)]
const fn build_bank_8d() -> [ScaleTable; 256] {
    let mut bank = [ScaleTable {
        coeff: gf8b::Elem(0),
        lo: [0; 16],
        hi: [0; 16],
    }; 256];
    let mut i = 0;
    while i < 256 {
        let c = gf8d::Elem(i as u8);
        let mut lo = [0u8; 16];
        let mut hi = [0u8; 16];
        let mut j = 0;
        while j < 16 {
            lo[j] = gf8d::Elem(j as u8).mul(c).0;
            hi[j] = gf8d::Elem((j as u8) << 4).mul(c).0;
            j += 1;
        }
        // `coeff` stores the coefficient byte; its AES-typed wrapper is inert
        // storage the field-agnostic shuffle kernels never read.
        bank[i] = ScaleTable {
            coeff: gf8b::Elem(i as u8),
            lo,
            hi,
        };
        i += 1;
    }
    bank
}

/// Return the shared `0x11D` nibble tables for a coefficient.
#[inline]
#[must_use]
pub fn scale_table_8d(coeff: gf8d::Elem) -> &'static ScaleTable {
    &SCALE_TABLE_BANK_8D[coeff.0 as usize]
}

/// Experimental GF(2^8)/`0x11B` affine-map bank, one
/// `VGF2P8AFFINEQB` matrix qword per coefficient.
///
/// This uses the same instruction convention as [`affine_8d`], but derives
/// each column from [`gf8b::Elem::mul`]. The bank supports experimental
/// affine kernels; dispatched `Gf8B` kernels continue to use native
/// `GF2P8MULB`.
#[allow(dead_code)]
static AFFINE_BANK_8B: [u64; 256] = build_affine_bank_8b();

#[allow(clippy::cast_possible_truncation)]
#[allow(dead_code)]
const fn build_affine_bank_8b() -> [u64; 256] {
    let mut bank = [0u64; 256];
    let mut i = 0;
    while i < 256 {
        let c = gf8b::Elem(i as u8);
        let mut map = 0u64;
        let mut r = 0u8;
        while r < 8 {
            let mut row = 0u8;
            let mut k = 0u8;
            while k < 8 {
                let column = c.mul(gf8b::Elem(1 << k)).0;
                row |= ((column >> (7 - r)) & 1) << k;
                k += 1;
            }
            map |= (row as u64) << (8 * r);
            r += 1;
        }
        bank[i] = map;
        i += 1;
    }
    bank
}

/// Return the experimental `VGF2P8AFFINEQB` matrix qword that multiplies by
/// `coeff` under `0x11B`.
#[inline]
#[must_use]
#[allow(dead_code)]
pub fn affine_8b(coeff: gf8b::Elem) -> u64 {
    AFFINE_BANK_8B[coeff.0 as usize]
}

/// GF(2^8)/`0x11D` affine-map bank, one `VGF2P8AFFINEQB` matrix qword per
/// coefficient.
///
/// Multiplication by a fixed coefficient in any GF(2^8) is an 8×8 GF(2)
/// linear map, and `VGF2P8AFFINEQB` applies an arbitrary such map per byte
/// lane — unlike `GF2P8MULB`, which bakes in the AES polynomial. The
/// instruction computes `dst.bit[b] = parity(map.byte[7 - b] & x)` per lane,
/// so the qword for coefficient `c` stores the product map transposed and
/// row-reversed: byte `r`, bit `k` is bit `7 - r` of the column `c * x^k`.
/// Every entry derives from the crate's own [`gf8d::Elem::mul`] oracle —
/// never from another library's tables — and the derivation is validated
/// exhaustively: the bank against a scalar model of the instruction in this
/// module's tests, the hardware kernels against the scalar oracle over all
/// 65 536 products in `kernel::tests`.
///
/// 2 KiB of rodata. Only the `Gf8D` GFNI kernels read the maps, but
/// `Gf8D::prepare` fills them on every target.
static AFFINE_BANK_8D: [u64; 256] = build_affine_bank_8d();

#[allow(clippy::cast_possible_truncation)]
const fn build_affine_bank_8d() -> [u64; 256] {
    let mut bank = [0u64; 256];
    let mut i = 0;
    while i < 256 {
        let c = gf8d::Elem(i as u8);
        let mut map = 0u64;
        // Row `r` of the instruction's matrix carries, across bits `k`, bit
        // `7 - r` of every column `c * x^k` of the multiplication map.
        let mut r = 0u8;
        while r < 8 {
            let mut row = 0u8;
            let mut k = 0u8;
            while k < 8 {
                let column = c.mul(gf8d::Elem(1 << k)).0;
                row |= ((column >> (7 - r)) & 1) << k;
                k += 1;
            }
            map |= (row as u64) << (8 * r);
            r += 1;
        }
        bank[i] = map;
        i += 1;
    }
    bank
}

/// Return the `VGF2P8AFFINEQB` matrix qword that multiplies by `coeff` under
/// `0x11D`.
#[inline]
#[must_use]
pub fn affine_8d(coeff: gf8d::Elem) -> u64 {
    AFFINE_BANK_8D[coeff.0 as usize]
}

/// Broadcast factors that express one GF(2^16) tower multiply as two
/// byte-wide GF(2^8) multiplies.
///
/// Interleaved source bytes are `[a, b]` meaning `a + b*u`. Multiplying by
/// `c0 + c1*u` gives
///
/// ```text
/// (c0*a + DELTA*c1*b) + (c1*a + (c0 + c1)*b) * u.
/// ```
///
/// The first term of each component multiplies the source in place; the
/// second multiplies the source with adjacent bytes swapped. So a single
/// alternating-coefficient byte multiply of `src` by [`TowerCoeff::same`],
/// `XORed` with the same of `swap16(src)` by [`TowerCoeff::cross`], produces
/// both components with no planar de-interleave.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct TowerCoeff {
    /// The GF(2^16) coefficient.
    pub coeff: gf16::Elem,
    /// Little-endian `[c0, c0 + c1]`, applied to the source directly.
    pub same: u16,
    /// Little-endian `[DELTA*c1, c1]`, applied to the byte-swapped source.
    pub cross: u16,
}

impl TowerCoeff {
    /// Derive the broadcast factors for `coeff`. Two base multiplies.
    #[inline]
    #[must_use]
    pub const fn new(coeff: gf16::Elem) -> Self {
        let (c0, c1) = coeff.to_components();
        let same = u16::from_le_bytes([c0.0, c0.add(c1).0]);
        let cross = u16::from_le_bytes([gf16::DELTA.mul(c1).0, c1.0]);
        Self { coeff, same, cross }
    }

    /// The four base-field coefficients in `[same.0, same.1, cross.0,
    /// cross.1]` order, i.e. `[c0, c0+c1, DELTA*c1, c1]`.
    #[inline]
    #[must_use]
    pub const fn factors(self) -> [gf8b::Elem; 4] {
        let [s0, s1] = self.same.to_le_bytes();
        let [x0, x1] = self.cross.to_le_bytes();
        [
            gf8b::Elem(s0),
            gf8b::Elem(s1),
            gf8b::Elem(x0),
            gf8b::Elem(x1),
        ]
    }
}

/// Nibble tables for the four base-field factors of a GF(2^16) coefficient.
///
/// Used by the shuffle backends (AVX2, SSSE3, NEON), which have no byte-wide
/// field multiply instruction and must emulate one per factor.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct TowerTables {
    /// The GF(2^16) coefficient.
    pub coeff: gf16::Elem,
    /// Tables for `[c0, c0+c1, DELTA*c1, c1]`. Entries 0 and 1 apply to the
    /// source's even and odd byte lanes; entries 2 and 3 apply to the
    /// adjacent-swapped source's even and odd lanes.
    pub factors: [ScaleTable; 4],
}

impl TowerTables {
    /// Build the four nibble tables for `coeff`.
    #[inline]
    #[must_use]
    pub fn new(coeff: gf16::Elem) -> Self {
        let f = TowerCoeff::new(coeff).factors();
        Self {
            coeff,
            factors: [
                *scale_table(f[0]),
                *scale_table(f[1]),
                *scale_table(f[2]),
                *scale_table(f[3]),
            ],
        }
    }
}

/// Split-nibble multiplication tables for one Fan–Paar `fp8` coefficient.
///
/// The canonical Fan–Paar byte field is *not* the AES field, so `GF2P8MULB`
/// cannot multiply its bytes; every backend emulates the byte product with a
/// 16-entry nibble lookup, exactly as [`ScaleTable`] does for GF(2^8). The
/// two share a shape and a 256-entry static bank, but the multiplication
/// used to *fill* them differs (fp8 tower vs. AES reduction).
#[allow(dead_code)]
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct FpScaleTable {
    /// `lo[i] = coeff * i` for the low nibble.
    pub lo: [u8; 16],
    /// `hi[i] = coeff * (i << 4)` for the high nibble.
    pub hi: [u8; 16],
}

impl FpScaleTable {
    /// Build the nibble tables for `coeff` in the Fan–Paar byte field.
    #[must_use]
    #[allow(clippy::cast_possible_truncation)]
    #[allow(dead_code)]
    pub const fn new(coeff: fp8::Elem) -> Self {
        let mut lo = [0u8; 16];
        let mut hi = [0u8; 16];
        let mut i = 0;
        while i < 16 {
            lo[i] = fp8::Elem(i as u8).mul(coeff).0;
            hi[i] = fp8::Elem((i as u8) << 4).mul(coeff).0;
            i += 1;
        }
        Self { lo, hi }
    }
}

/// Shared nibble-table bank, one entry per Fan–Paar `fp8` coefficient.
///
/// A separate bank from [`SCALE_TABLE_BANK`]: fp8 is a different field, so its
/// products are different bytes. 8 KiB, resident in L1/L2, touched only by
/// Fan–Paar kernels.
#[allow(dead_code)]
static FP_SCALE_TABLE_BANK: [FpScaleTable; 256] = build_fp_bank();

#[allow(clippy::cast_possible_truncation)]
#[allow(dead_code)]
const fn build_fp_bank() -> [FpScaleTable; 256] {
    let mut bank = [FpScaleTable::new(fp8::Elem(0)); 256];
    let mut i = 0;
    while i < 256 {
        bank[i] = FpScaleTable::new(fp8::Elem(i as u8));
        i += 1;
    }
    bank
}

/// Return the shared fp8 nibble tables for a Fan–Paar byte coefficient.
#[inline]
#[must_use]
#[allow(dead_code)]
pub fn fp_scale_table(coeff: fp8::Elem) -> &'static FpScaleTable {
    &FP_SCALE_TABLE_BANK[coeff.0 as usize]
}

/// Four nibble-table factors a shuffle kernel consumes lane-by-lane.
///
/// Both GF(2^16) over GF(2^8) and Fan–Paar GF(2^16) over `fp8` reduce a
/// 16-bit fixed-coefficient scale to four byte-wide multiplies — factors 0/1
/// on the source's even/odd bytes, factors 2/3 on the adjacent-swapped
/// source's. The multiply core in `kernel::x86::gf16` reads only the nibble
/// tables, not the field they encode, so this trait lets the two fields share
/// that core. Entries 0 and 1 apply to the source's even and odd byte lanes;
/// 2 and 3 to the swapped source's.
#[cfg(all(feature = "simd", any(target_arch = "x86", target_arch = "x86_64")))]
pub(crate) trait NibbleFactors {
    /// The low-nibble table of factor `i` (`0 <= i < 4`).
    fn lo(&self, i: usize) -> &[u8; 16];
    /// The high-nibble table of factor `i` (`0 <= i < 4`).
    fn hi(&self, i: usize) -> &[u8; 16];
}

#[cfg(all(feature = "simd", any(target_arch = "x86", target_arch = "x86_64")))]
impl NibbleFactors for TowerTables {
    #[inline]
    fn lo(&self, i: usize) -> &[u8; 16] {
        &self.factors[i].lo
    }
    #[inline]
    fn hi(&self, i: usize) -> &[u8; 16] {
        &self.factors[i].hi
    }
}

/// Nibble tables for the four `fp8` factors of a Fan–Paar GF(2^16)
/// coefficient.
///
/// The Fan–Paar tower relation is `X² + alpha·X + 1 = 0` with `alpha` the
/// `fp8` tower generator, so for `c = c0 + c1·X` and `x = x0 + x1·X`:
///
/// ```text
/// r0 = c0·x0 ^ c1·x1
/// r1 = (c0 ^ alpha·c1)·x1 ^ c1·x0
/// ```
///
/// `alpha·(c1·x1) = (alpha·c1)·x1` because `alpha`, `c1`, `x1` all lie in the
/// `fp8` subfield and the subfield commutes — so the `mul_alpha` fold lands in
/// coefficient preparation, not in the kernel. The four factors are therefore
/// `[c0, c0 ^ alpha·c1, c1, c1]`: 0/1 scale the source, 2/3 (both `c1`) the
/// swapped source. This is the same four-table shape as [`TowerTables`], which
/// is why the shuffle multiply core is shared.
#[allow(dead_code)]
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct FpTowerTables {
    /// The Fan–Paar GF(2^16) coefficient, recovered for the portable scalar
    /// tail of a vector kernel.
    pub coeff: fp16::Elem,
    /// Tables for `[c0, c0 ^ alpha·c1, c1, c1]`.
    pub factors: [FpScaleTable; 4],
}

impl FpTowerTables {
    /// Build the four `fp8` nibble tables for `coeff`.
    #[inline]
    #[must_use]
    #[allow(dead_code)]
    pub fn new(coeff: fp16::Elem) -> Self {
        let (c0, c1) = coeff.to_components();
        // `c1.mul_alpha()` is `alpha·c1` in the fp8 subfield.
        let b = c0.add(c1.mul_alpha());
        let [f0, f1, f2, f3] = [c0, b, c1, c1];
        Self {
            coeff,
            factors: [
                *fp_scale_table(f0),
                *fp_scale_table(f1),
                *fp_scale_table(f2),
                *fp_scale_table(f3),
            ],
        }
    }
}
#[cfg(all(feature = "simd", any(target_arch = "x86", target_arch = "x86_64")))]
impl NibbleFactors for FpTowerTables {
    #[inline]
    fn lo(&self, i: usize) -> &[u8; 16] {
        &self.factors[i].lo
    }
    #[inline]
    fn hi(&self, i: usize) -> &[u8; 16] {
        &self.factors[i].hi
    }
}

/// Period-2 alternating subfield coefficients for a two-level tower multiply.
///
/// Generalizes [`TowerCoeff`] one field level up. For a level-2 element
/// `c = c0 + c1·u` over `u² + u + DELTA`, the level-2 multiply is
///
/// ```text
/// c·x = lane_mul(x,      same)  ^  lane_mul(swap(x), cross)
/// ```
///
/// where `lane_mul` multiplies subfield elements under alternating `same`/
/// `cross` coefficients and `swap` exchanges the two subfield halves of each
/// element. `same = [c0, c0 + c1]` applies to the source's even and odd
/// subfield lanes; `cross = [DELTA·c1, c1]` to the half-swapped source. This
/// is exactly the identity the GF(2^16) kernel builds at byte granularity
/// via [`TowerCoeff`], one level up.
///
/// The form carries subfield elements, not byte broadcasts, so it is
/// representation-agnostic: GF(2^32) builds it over `gf16::Elem` and GF(2^64)
/// over `gf32::Elem`, each with its own `DELTA`.
/// This is a representation-agnostic seam: every level-2 SIMD tower — the
/// GFNI x86 path today, the aarch64/wasm paths to come — derives it from the
/// subfield element, so it is dead weight only on targets with no such
/// kernel.
#[allow(dead_code)]
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct Tower2Coeff<E> {
    /// `[c0, c0 + c1]`, applied to the even and odd subfield lanes of the
    /// source.
    pub same: [E; 2],
    /// `[DELTA·c1, c1]`, applied to the even and odd lanes of the half-swapped
    /// source.
    pub cross: [E; 2],
}

impl<E: crate::field::Elem> Tower2Coeff<E> {
    /// Derive the alternating subfield coefficient pair for `c0 + c1·u` over
    /// `u² + u + delta`.
    #[inline]
    #[must_use]
    #[allow(dead_code)]
    pub fn derive(c0: E, c1: E, delta: E) -> Self {
        Self {
            same: [c0, c0.add(c1)],
            cross: [delta.mul(c1), c1],
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::field::gf8b::Elem as E8;

    #[test]
    fn nibble_tables_reconstruct_the_product() {
        for c in 0..=u8::MAX {
            let table = scale_table(E8(c));
            for x in 0..=u8::MAX {
                let split = table.lo[(x & 0x0f) as usize] ^ table.hi[(x >> 4) as usize];
                assert_eq!(split, E8(x).mul(E8(c)).0, "coeff {c:#04x} value {x:#04x}");
            }
        }
    }

    #[test]
    #[allow(clippy::cast_possible_truncation)]
    fn affine_bank_matches_the_instruction_model() {
        // The semantics `VGF2P8AFFINEQB` documents, modeled in scalar code:
        // `dst.bit[b] = parity(map.byte[7 - b] & x)`. This proves both const
        // banks use the hardware bit-order convention.
        fn apply(map: u64, x: u8) -> u8 {
            let mut out = 0u8;
            for b in 0..8u8 {
                let row = (map >> (8 * (7 - b))) as u8;
                let parity = (row & x).count_ones() & 1;
                out |= u8::try_from(parity).unwrap() << b;
            }
            out
        }

        for c in 0..=u8::MAX {
            let aes_map = affine_8b(gf8b::Elem(c));
            let rs_map = affine_8d(gf8d::Elem(c));
            for x in 0..=u8::MAX {
                assert_eq!(
                    apply(aes_map, x),
                    gf8b::Elem(c).mul(gf8b::Elem(x)).0,
                    "0x11B coeff {c:#04x} value {x:#04x}"
                );
                assert_eq!(
                    apply(rs_map, x),
                    gf8d::Elem(c).mul(gf8d::Elem(x)).0,
                    "0x11D coeff {c:#04x} value {x:#04x}"
                );
            }
        }
    }

    #[test]
    fn tower_factors_reconstruct_the_product() {
        // Exercise the identity the SIMD kernels rely on, scalar-side.
        for coeff in [0u16, 1, 0x0108, 0x2000, 0xffff, 0x1234] {
            let coeff = gf16::Elem(coeff);
            let tc = TowerCoeff::new(coeff);
            let [f0, f1, f2, f3] = tc.factors();
            for value in [0u16, 1, 0x00ff, 0xff00, 0xbeef] {
                let value = gf16::Elem(value);
                let [a, b] = value.to_bytes();
                let (a, b) = (E8(a), E8(b));
                // even lane: same.0 * a  ^  cross.0 * b   (b is a's swap partner)
                // odd  lane: same.1 * b  ^  cross.1 * a
                let even = f0.mul(a).add(f2.mul(b));
                let odd = f1.mul(b).add(f3.mul(a));
                let got = gf16::Elem::from_bytes([even.0, odd.0]);
                assert_eq!(got, value.mul(coeff), "{value:?} * {coeff:?}");
            }
        }
    }

    #[test]
    fn const_built_banks_match_their_runtime_builders() {
        // The banks are `const`-evaluated into rodata; running the same
        // builders at runtime must reproduce them bit for bit. This is the
        // only place the builder bodies execute, and it doubles as a
        // differential between compile-time and runtime evaluation.
        assert_eq!(build_bank(), SCALE_TABLE_BANK);
        assert_eq!(build_bank_8d(), SCALE_TABLE_BANK_8D);
        assert_eq!(build_affine_bank_8d(), AFFINE_BANK_8D);
        assert_eq!(build_fp_bank(), FP_SCALE_TABLE_BANK);
        for c in 0..=u8::MAX {
            let runtime = ScaleTable::new(E8(c));
            assert_eq!(&runtime, scale_table(E8(c)), "0x11B coeff {c:#04x}");
            let runtime = FpScaleTable::new(fp8::Elem(c));
            assert_eq!(&runtime, fp_scale_table(fp8::Elem(c)), "fp8 coeff {c:#04x}");
        }
    }

    #[test]
    fn experimental_affine_bank_matches_its_runtime_builder() {
        assert_eq!(build_affine_bank_8b(), AFFINE_BANK_8B);
    }
}