hekate_math/towers/

block64.rs

// SPDX-License-Identifier: Apache-2.0
// This file is part of the hekate-math project.
// Copyright (C) 2026 Andrei Kochergin <andrei@oumuamua.dev>
// Copyright (C) 2026 Oumuamua Labs <info@oumuamua.dev>. All rights reserved.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

//! GF(2^64) as `Block32[X] / (X^2 + X + τ)`, its flat basis
//! modulo x^64 + x^4 + x^3 + x + 1, and the aarch64 NEON kernels.

use crate::algebra::impl_binary_field_extras;
use crate::constants::FLAT_TO_TOWER_BIT_MASKS_64;
use crate::towers::bit::Bit;
use crate::towers::block8::Block8;
use crate::towers::block16::Block16;
use crate::towers::block32::Block32;
use crate::{
    BinaryFieldExtras, CanonicalDeserialize, CanonicalSerialize, Flat, FlatPromote, HardwareField,
    PackableField, PackedFlat, TowerField, constants,
};
use core::ops::{Add, AddAssign, BitXor, BitXorAssign, Mul, MulAssign, Sub, SubAssign};
use serde::{Deserialize, Serialize};
use zeroize::Zeroize;

#[cfg(not(feature = "table-math"))]
#[repr(align(64))]
struct CtConvertBasisU64<const N: usize>([u64; N]);

#[cfg(not(feature = "table-math"))]
static TOWER_TO_FLAT_BASIS_64: CtConvertBasisU64<64> =
    CtConvertBasisU64(constants::RAW_TOWER_TO_FLAT_64);

#[cfg(not(feature = "table-math"))]
static FLAT_TO_TOWER_BASIS_64: CtConvertBasisU64<64> =
    CtConvertBasisU64(constants::RAW_FLAT_TO_TOWER_64);

/// GF(2^64) in the tower basis: `lo + hi·X`
/// over [`Block32`], `lo` in the low 32 bits.
#[derive(Copy, Clone, Default, Debug, Eq, PartialEq, Serialize, Deserialize, Zeroize)]
#[repr(transparent)]
pub struct Block64(pub u64);

impl Block64 {
    /// τ of the next level, `Block64[X] / (X^2 + X + τ)`;
    /// equals [`TowerField::EXTENSION_TAU`].
    pub const TAU: Self = Block64(0x2000_0000_0000_0000);

    /// Builds `lo + hi·X` from its halves.
    pub fn new(lo: Block32, hi: Block32) -> Self {
        Self((hi.0 as u64) << 32 | (lo.0 as u64))
    }

    /// Returns the halves `(lo, hi)` of `lo + hi·X`.
    #[inline(always)]
    pub fn split(self) -> (Block32, Block32) {
        (Block32(self.0 as u32), Block32((self.0 >> 32) as u32))
    }

    #[inline(always)]
    pub(crate) fn mul_tau(self) -> Self {
        let (a0, a1) = self.split();
        let t = a1.mul_tau();

        Self::new(t.mul_tau(), (a0 + a1).mul_tau())
    }
}

impl TowerField for Block64 {
    const BITS: usize = 64;
    const ZERO: Self = Block64(0);
    const ONE: Self = Block64(1);

    const EXTENSION_TAU: Self = Self::TAU;

    fn invert(&self) -> Self {
        let (l, h) = self.split();

        let norm = h.square().mul_tau() + h * l + l.square();
        let norm_inv = norm.invert();

        let res_hi = h * norm_inv;
        let res_lo = (h + l) * norm_inv;

        Self::new(res_lo, res_hi)
    }

    fn from_uniform_bytes(bytes: &[u8; 32]) -> Self {
        let mut buf = [0u8; 8];
        buf.copy_from_slice(&bytes[0..8]);

        Self(u64::from_le_bytes(buf))
    }
}

impl Add for Block64 {
    type Output = Self;

    fn add(self, rhs: Self) -> Self {
        Self(self.0.bitxor(rhs.0))
    }
}

impl Sub for Block64 {
    type Output = Self;

    fn sub(self, rhs: Self) -> Self {
        self.add(rhs)
    }
}

impl Mul for Block64 {
    type Output = Self;

    #[inline]
    fn mul(self, rhs: Self) -> Self {
        #[cfg(pmull)]
        {
            mul_iso_64(self, rhs)
        }

        #[cfg(not(pmull))]
        {
            let (a0, a1) = self.split();
            let (b0, b1) = rhs.split();

            let v0 = a0 * b0;
            let v1 = a1 * b1;
            let vs = (a0 + a1) * (b0 + b1);

            Self::new(v0 + v1.mul_tau(), v0 + vs)
        }
    }
}

impl AddAssign for Block64 {
    fn add_assign(&mut self, rhs: Self) {
        self.0.bitxor_assign(rhs.0);
    }
}

impl SubAssign for Block64 {
    fn sub_assign(&mut self, rhs: Self) {
        self.0.bitxor_assign(rhs.0);
    }
}

impl MulAssign for Block64 {
    fn mul_assign(&mut self, rhs: Self) {
        *self = *self * rhs;
    }
}

impl CanonicalSerialize for Block64 {
    fn serialized_size(&self) -> usize {
        8
    }

    fn serialize(&self, writer: &mut [u8]) -> Result<(), ()> {
        if writer.len() < 8 {
            return Err(());
        }

        writer[..8].copy_from_slice(&self.0.to_le_bytes());

        Ok(())
    }
}

impl CanonicalDeserialize for Block64 {
    fn deserialize(bytes: &[u8]) -> Result<Self, ()> {
        if bytes.len() < 8 {
            return Err(());
        }

        let mut buf = [0u8; 8];
        buf.copy_from_slice(&bytes[0..8]);

        Ok(Self(u64::from_le_bytes(buf)))
    }
}

impl From<u8> for Block64 {
    #[inline(always)]
    fn from(val: u8) -> Self {
        Self(val as u64)
    }
}

impl From<u32> for Block64 {
    #[inline(always)]
    fn from(val: u32) -> Self {
        Self::from(val as u64)
    }
}

impl From<u64> for Block64 {
    #[inline(always)]
    fn from(val: u64) -> Self {
        Self(val)
    }
}

impl From<u128> for Block64 {
    #[inline(always)]
    fn from(val: u128) -> Self {
        Self(val as u64)
    }
}

// ========================================
// FIELD LIFTING
// ========================================

impl From<Bit> for Block64 {
    #[inline(always)]
    fn from(val: Bit) -> Self {
        Self(val.get() as u64)
    }
}

impl From<Block8> for Block64 {
    #[inline(always)]
    fn from(val: Block8) -> Self {
        Self(val.0 as u64)
    }
}

impl From<Block16> for Block64 {
    #[inline(always)]
    fn from(val: Block16) -> Self {
        Self(val.0 as u64)
    }
}

impl From<Block32> for Block64 {
    #[inline(always)]
    fn from(val: Block32) -> Self {
        Self(val.0 as u64)
    }
}

// ===================================
// PACKED BLOCK 64 (Width = 2)
// ===================================

/// Lanes in [`PackedBlock64`].
pub const PACKED_WIDTH_64: usize = 2;

/// Two [`Block64`] lanes; operators act lane by lane.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Default)]
#[repr(C, align(16))] // 128-bit alignment
pub struct PackedBlock64(pub [Block64; PACKED_WIDTH_64]);

impl PackedBlock64 {
    /// Returns the packed value with every lane zero.
    #[inline(always)]
    pub fn zero() -> Self {
        Self([Block64::ZERO; PACKED_WIDTH_64])
    }
}

impl PackableField for Block64 {
    type Packed = PackedBlock64;

    const WIDTH: usize = PACKED_WIDTH_64;

    #[inline(always)]
    fn pack(chunk: &[Self]) -> Self::Packed {
        assert!(
            chunk.len() >= PACKED_WIDTH_64,
            "PackableField::pack: input slice too short",
        );

        let mut arr = [Self::ZERO; PACKED_WIDTH_64];
        arr.copy_from_slice(&chunk[..PACKED_WIDTH_64]);

        PackedBlock64(arr)
    }

    #[inline(always)]
    fn unpack(packed: Self::Packed, output: &mut [Self]) {
        assert!(
            output.len() >= PACKED_WIDTH_64,
            "PackableField::unpack: output slice too short",
        );

        output[..PACKED_WIDTH_64].copy_from_slice(&packed.0);
    }
}

impl Add for PackedBlock64 {
    type Output = Self;

    #[inline(always)]
    fn add(self, rhs: Self) -> Self {
        let mut res = [Block64::ZERO; PACKED_WIDTH_64];
        for ((out, l), r) in res.iter_mut().zip(self.0.iter()).zip(rhs.0.iter()) {
            *out = *l + *r;
        }

        Self(res)
    }
}

impl AddAssign for PackedBlock64 {
    #[inline(always)]
    fn add_assign(&mut self, rhs: Self) {
        for (l, r) in self.0.iter_mut().zip(rhs.0.iter()) {
            *l += *r;
        }
    }
}

impl Sub for PackedBlock64 {
    type Output = Self;

    #[inline(always)]
    fn sub(self, rhs: Self) -> Self {
        self.add(rhs)
    }
}

impl SubAssign for PackedBlock64 {
    #[inline(always)]
    fn sub_assign(&mut self, rhs: Self) {
        self.add_assign(rhs);
    }
}

impl Mul for PackedBlock64 {
    type Output = Self;

    #[inline(always)]
    fn mul(self, rhs: Self) -> Self {
        let mut res = [Block64::ZERO; PACKED_WIDTH_64];
        for ((out, l), r) in res.iter_mut().zip(self.0.iter()).zip(rhs.0.iter()) {
            *out = *l * *r;
        }

        Self(res)
    }
}

impl MulAssign for PackedBlock64 {
    #[inline(always)]
    fn mul_assign(&mut self, rhs: Self) {
        for (l, r) in self.0.iter_mut().zip(rhs.0.iter()) {
            *l *= *r;
        }
    }
}

impl Mul<Block64> for PackedBlock64 {
    type Output = Self;

    #[inline(always)]
    fn mul(self, rhs: Block64) -> Self {
        let mut res = [Block64::ZERO; PACKED_WIDTH_64];
        for (out, v) in res.iter_mut().zip(self.0.iter()) {
            *out = *v * rhs;
        }

        Self(res)
    }
}

// ===================================
// Hardware Field
// ===================================

impl HardwareField for Block64 {
    #[inline(always)]
    fn to_hardware(self) -> Flat<Self> {
        #[cfg(feature = "table-math")]
        {
            Flat::from_raw(apply_matrix_64(self, &constants::TOWER_TO_FLAT_64))
        }

        #[cfg(not(feature = "table-math"))]
        {
            Flat::from_raw(Block64(map_ct_64(self.0, &TOWER_TO_FLAT_BASIS_64.0)))
        }
    }

    #[inline(always)]
    fn from_hardware(value: Flat<Self>) -> Self {
        let value = value.into_raw();

        #[cfg(feature = "table-math")]
        {
            apply_matrix_64(value, &constants::FLAT_TO_TOWER_64)
        }

        #[cfg(not(feature = "table-math"))]
        {
            Block64(map_ct_64(value.0, &FLAT_TO_TOWER_BASIS_64.0))
        }
    }

    #[inline(always)]
    fn add_hardware(lhs: Flat<Self>, rhs: Flat<Self>) -> Flat<Self> {
        Flat::from_raw(lhs.into_raw() + rhs.into_raw())
    }

    #[inline(always)]
    fn add_hardware_packed(lhs: PackedFlat<Self>, rhs: PackedFlat<Self>) -> PackedFlat<Self> {
        let lhs = lhs.into_raw();
        let rhs = rhs.into_raw();

        #[cfg(target_arch = "aarch64")]
        {
            PackedFlat::from_raw(neon::add_packed_64(lhs, rhs))
        }

        #[cfg(not(target_arch = "aarch64"))]
        {
            PackedFlat::from_raw(lhs + rhs)
        }
    }

    #[inline(always)]
    fn mul_hardware(lhs: Flat<Self>, rhs: Flat<Self>) -> Flat<Self> {
        let lhs = lhs.into_raw();
        let rhs = rhs.into_raw();

        #[cfg(pmull)]
        {
            Flat::from_raw(neon::mul_flat_64(lhs, rhs))
        }

        #[cfg(not(pmull))]
        {
            let a_tower = Self::from_hardware(Flat::from_raw(lhs));
            let b_tower = Self::from_hardware(Flat::from_raw(rhs));

            (a_tower * b_tower).to_hardware()
        }
    }

    #[inline(always)]
    fn mul_hardware_packed(lhs: PackedFlat<Self>, rhs: PackedFlat<Self>) -> PackedFlat<Self> {
        let lhs = lhs.into_raw();
        let rhs = rhs.into_raw();

        #[cfg(pmull)]
        {
            PackedFlat::from_raw(neon::mul_flat_packed_64(lhs, rhs))
        }

        #[cfg(not(pmull))]
        {
            let mut l = [Self::ZERO; <Self as PackableField>::WIDTH];
            let mut r = [Self::ZERO; <Self as PackableField>::WIDTH];
            let mut res = [Self::ZERO; <Self as PackableField>::WIDTH];

            Self::unpack(lhs, &mut l);
            Self::unpack(rhs, &mut r);

            for i in 0..<Self as PackableField>::WIDTH {
                res[i] = Self::mul_hardware(Flat::from_raw(l[i]), Flat::from_raw(r[i])).into_raw();
            }

            PackedFlat::from_raw(Self::pack(&res))
        }
    }

    #[inline(always)]
    fn mul_hardware_scalar_packed(lhs: PackedFlat<Self>, rhs: Flat<Self>) -> PackedFlat<Self> {
        let broadcasted = PackedBlock64([rhs.into_raw(); PACKED_WIDTH_64]);
        Self::mul_hardware_packed(lhs, PackedFlat::from_raw(broadcasted))
    }

    #[inline(always)]
    fn tower_bit_from_hardware(value: Flat<Self>, bit_idx: usize) -> u8 {
        let mask = FLAT_TO_TOWER_BIT_MASKS_64[bit_idx];

        // Parity of (x & mask) without popcount
        // Folds 64 bits down to 4,
        // then uses a lookup table.
        let mut v = value.into_raw().0 & mask;
        v ^= v >> 32;
        v ^= v >> 16;
        v ^= v >> 8;
        v ^= v >> 4;

        let idx = (v & 0xF) as u8;

        // Nibble parity lookup encoded
        // in a 16-bit constant (0x6996).
        ((0x6996u16 >> idx) & 1) as u8
    }
}

impl FlatPromote<Block8> for Block64 {
    #[inline(always)]
    fn promote_flat(val: Flat<Block8>) -> Flat<Self> {
        let val = val.into_raw();

        #[cfg(not(feature = "table-math"))]
        {
            Flat::from_raw(Block64(lift_ct_64::<8>(
                val.0 as u64,
                &constants::LIFT_BASIS_8_TO_64,
            )))
        }

        #[cfg(feature = "table-math")]
        {
            Flat::from_raw(Block64(constants::LIFT_TABLE_8_TO_64[val.0 as usize]))
        }
    }
}

impl FlatPromote<Block16> for Block64 {
    #[inline(always)]
    fn promote_flat(val: Flat<Block16>) -> Flat<Self> {
        let val = val.into_raw();

        #[cfg(not(feature = "table-math"))]
        {
            Flat::from_raw(Block64(lift_ct_64::<16>(
                val.0 as u64,
                &constants::LIFT_BASIS_16_TO_64,
            )))
        }

        #[cfg(feature = "table-math")]
        {
            let v = val.0;
            let res = constants::PROMOTE_16_BYTE_0_TO_64[(v & 0xFF) as usize]
                ^ constants::PROMOTE_16_BYTE_1_TO_64[(v >> 8) as usize];

            Flat::from_raw(Block64(res))
        }
    }
}

impl FlatPromote<Block32> for Block64 {
    #[inline(always)]
    fn promote_flat(val: Flat<Block32>) -> Flat<Self> {
        let val = val.into_raw();

        #[cfg(not(feature = "table-math"))]
        {
            Flat::from_raw(Block64(lift_ct_64::<32>(
                val.0 as u64,
                &constants::LIFT_BASIS_32_TO_64,
            )))
        }

        #[cfg(feature = "table-math")]
        {
            let v = val.0;
            let res = constants::PROMOTE_32_BYTE_0_TO_64[(v & 0xFF) as usize]
                ^ constants::PROMOTE_32_BYTE_1_TO_64[((v >> 8) & 0xFF) as usize]
                ^ constants::PROMOTE_32_BYTE_2_TO_64[((v >> 16) & 0xFF) as usize]
                ^ constants::PROMOTE_32_BYTE_3_TO_64[(v >> 24) as usize];

            Flat::from_raw(Block64(res))
        }
    }
}

// ===========================================
// Binary Field Extras
// ===========================================

impl_binary_field_extras!(Block64, map_ct_64, TRACE_MASK_64, SOLVE_QUADRATIC_BASIS_64);

// ===========================================
// UTILS
// ===========================================

/// Returns `a * b` through the flat basis, φ⁻¹(φ(a) · φ(b)), with
/// the PMULL kernel; aarch64 with the `aes` target feature only.
#[cfg(pmull)]
#[inline(always)]
pub fn mul_iso_64(a: Block64, b: Block64) -> Block64 {
    let a_flat = a.to_hardware();
    let b_flat = b.to_hardware();

    let c_flat = Flat::from_raw(neon::mul_flat_64(a_flat.into_raw(), b_flat.into_raw()));

    c_flat.to_tower()
}

#[cfg(feature = "table-math")]
#[inline(always)]
pub(crate) fn apply_matrix_64(val: Block64, table: &[u64; 2048]) -> Block64 {
    let mut res = 0u64;
    let v = val.0;

    // 8 lookups (8-bit window)
    for i in 0..8 {
        let byte = (v >> (i * 8)) & 0xFF;
        let idx = (i * 256) + (byte as usize);

        // SAFETY:
        // idx = i * 256 + byte < 8 * 256 = table.len().
        res ^= unsafe { *table.get_unchecked(idx) };
    }

    Block64(res)
}

#[cfg(not(feature = "table-math"))]
#[inline(always)]
fn lift_ct_64<const N: usize>(x: u64, basis: &[u64; N]) -> u64 {
    let mut acc = 0u64;
    let mut i = 0usize;

    while i < N {
        let bit = (x >> i) & 1;
        let mask = 0u64.wrapping_sub(bit);

        acc ^= basis[i] & mask;
        i += 1;
    }

    acc
}

#[inline(always)]
fn map_ct_64(x: u64, basis: &[u64; 64]) -> u64 {
    let mut acc = 0u64;
    let mut i = 0usize;

    while i < 64 {
        let bit = (x >> i) & 1;
        let mask = 0u64.wrapping_sub(bit);

        acc ^= basis[i] & mask;
        i += 1;
    }

    acc
}

// ===========================================
// SIMD INSTRUCTIONS
// ===========================================

#[cfg(target_arch = "aarch64")]
mod neon {
    use super::*;
    use core::arch::aarch64::*;
    use core::mem::transmute;

    const _: () = assert!(constants::POLY_64 == 0x1b, "verus twins hardcode R = 0x1b");

    #[inline(always)]
    pub fn add_packed_64(lhs: PackedBlock64, rhs: PackedBlock64) -> PackedBlock64 {
        unsafe {
            let l: uint8x16_t = transmute::<[Block64; PACKED_WIDTH_64], uint8x16_t>(lhs.0);
            let r: uint8x16_t = transmute::<[Block64; PACKED_WIDTH_64], uint8x16_t>(rhs.0);
            let res = veorq_u8(l, r);
            let out: [Block64; PACKED_WIDTH_64] =
                transmute::<uint8x16_t, [Block64; PACKED_WIDTH_64]>(res);

            PackedBlock64(out)
        }
    }

    #[cfg(pmull)]
    #[inline(always)]
    pub fn mul_flat_packed_64(lhs: PackedBlock64, rhs: PackedBlock64) -> PackedBlock64 {
        unsafe {
            let ap: poly64x2_t = transmute(lhs.0);
            let bp: poly64x2_t = transmute(rhs.0);

            // vdupq_n_p64, not [c; 2]: the repeat
            // form lowers to a memcpy libcall.
            let rv: poly64x2_t = vdupq_n_p64(constants::POLY_64);

            // Both lanes per stage: PMULL/PMULL2 pairs,
            // uzp1/uzp2 regroup [lo, hi] products by lane.
            let p0: uint64x2_t =
                transmute(vmull_p64(vgetq_lane_p64::<0>(ap), vgetq_lane_p64::<0>(bp)));
            let p1: uint64x2_t = transmute(vmull_high_p64(ap, bp));

            let los = vuzp1q_u64(p0, p1);
            let his = vuzp2q_u64(p0, p1);
            let hisp: poly64x2_t = transmute(his);

            let hr0: uint64x2_t = transmute(vmull_p64(
                vgetq_lane_p64::<0>(hisp),
                vgetq_lane_p64::<0>(rv),
            ));
            let hr1: uint64x2_t = transmute(vmull_high_p64(hisp, rv));

            let folded = vuzp1q_u64(hr0, hr1);
            let carries = vuzp2q_u64(hr0, hr1);
            let cp: poly64x2_t = transmute(carries);

            let cr0: uint64x2_t =
                transmute(vmull_p64(vgetq_lane_p64::<0>(cp), vgetq_lane_p64::<0>(rv)));
            let cr1: uint64x2_t = transmute(vmull_high_p64(cp, rv));
            let carry_red = vuzp1q_u64(cr0, cr1);

            let res = veorq_u64(veorq_u64(los, folded), carry_red);

            PackedBlock64(transmute::<uint64x2_t, [Block64; 2]>(res))
        }
    }

    #[cfg(pmull)]
    #[inline(always)]
    pub fn mul_flat_64(a: Block64, b: Block64) -> Block64 {
        unsafe {
            // vdupq_n_p64, not [c; 2]: the repeat
            // form lowers to a memcpy libcall.
            let rv: poly64x2_t = vdupq_n_p64(constants::POLY_64);

            // Reduction P(x) = x^64 + R(x), R(x) = 0x1b
            let prod: uint8x16_t = transmute(vmull_p64(a.0, b.0));
            let h_red: uint8x16_t = transmute(vmull_high_p64(
                transmute::<uint8x16_t, poly64x2_t>(prod),
                rv,
            ));

            // carry = hi lane of h_red < 2^3;
            // deg(carry·R) < 7: no third fold.
            let c_red: uint8x16_t = transmute(vmull_high_p64(
                transmute::<uint8x16_t, poly64x2_t>(h_red),
                rv,
            ));

            let res = veorq_u8(veorq_u8(prod, h_red), c_red);

            Block64(vgetq_lane_u64(vreinterpretq_u64_u8(res), 0))
        }
    }
}

// ==================================
// BLOCK 64 TESTS
// ==================================

#[cfg(test)]
mod tests {
    use super::*;
    use proptest::prelude::*;
    use rand::{RngExt, rng};

    // ==================================
    // BASIC
    // ==================================

    #[test]
    fn tower_constants() {
        // Check that tau is propagated correctly
        // For Block64, tau must be (0, 1) from Block32.
        let tau64 = Block64::EXTENSION_TAU;
        let (lo64, hi64) = tau64.split();
        assert_eq!(lo64, Block32::ZERO);
        assert_eq!(hi64, Block32::TAU);
    }

    #[test]
    fn add_truth() {
        let zero = Block64::ZERO;
        let one = Block64::ONE;

        assert_eq!(zero + zero, zero);
        assert_eq!(zero + one, one);
        assert_eq!(one + zero, one);
        assert_eq!(one + one, zero);
    }

    #[test]
    fn mul_truth() {
        let zero = Block64::ZERO;
        let one = Block64::ONE;

        assert_eq!(zero * zero, zero);
        assert_eq!(zero * one, zero);
        assert_eq!(one * one, one);
    }

    #[test]
    fn add() {
        // 5 ^ 3 = 6
        // 101 ^ 011 = 110
        assert_eq!(Block64(5) + Block64(3), Block64(6));
    }

    #[test]
    fn mul_simple() {
        // Check for prime numbers (without overflow)
        // x^1 * x^1 = x^2 (2 * 2 = 4)
        assert_eq!(Block64(2) * Block64(2), Block64(4));
    }

    #[test]
    fn mul_overflow() {
        // Reduction verification (AES test vectors)
        // Example from the AES specification:
        // 0x57 * 0x83 = 0xC1
        assert_eq!(Block64(0x57) * Block64(0x83), Block64(0xC1));
    }

    #[test]
    fn karatsuba_correctness() {
        // Let's check using Block64 as an example
        // Let A = X (hi=1, lo=0)
        // Let B = X (hi=1, lo=0)
        // A * B = X^2
        // According to the rule:
        // X^2 = X + tau
        // Where tau for Block32 = 0x2000_0000.
        // So the result should be:
        // hi=1 (X), lo=0x20 (tau)

        // Construct X manually
        let x = Block64::new(Block32::ZERO, Block32::ONE);
        let squared = x * x;

        // Verify result via splitting
        let (res_lo, res_hi) = squared.split();

        assert_eq!(res_hi, Block32::ONE, "X^2 should contain X component");
        assert_eq!(
            res_lo,
            Block32(0x2000_0000),
            "X^2 should contain tau component (0x2000_0000)"
        );
    }

    #[test]
    fn mul_tau_matches_mul() {
        let mut rng = rng();
        let edges = [Block64::ZERO, Block64::ONE, Block64::TAU, Block64(u64::MAX)];
        let random = (0..1000).map(|_| Block64(rng.random()));

        for x in edges.into_iter().chain(random) {
            assert_eq!(
                x.mul_tau(),
                x * Block64::EXTENSION_TAU,
                "Block64 mul_tau mismatch at {:#018x}",
                x.0
            );
        }
    }

    #[test]
    fn security_zeroize() {
        let mut secret_val = Block64::from(0xDEAD_BEEF_CAFE_BABE_u64);
        assert_ne!(secret_val, Block64::ZERO);

        secret_val.zeroize();

        assert_eq!(secret_val, Block64::ZERO);
        assert_eq!(secret_val.0, 0, "Block64 memory leak detected");
    }

    #[test]
    fn invert_zero() {
        // Zero check
        assert_eq!(
            Block64::ZERO.invert(),
            Block64::ZERO,
            "invert(0) must return 0"
        );
    }

    #[test]
    fn inversion_random() {
        let mut rng = rng();
        for _ in 0..1000 {
            let val = Block64(rng.random());
            if val != Block64::ZERO {
                let inv = val.invert();
                assert_eq!(
                    val * inv,
                    Block64::ONE,
                    "Inversion identity failed: a * a^-1 != 1"
                );
            }
        }
    }

    #[test]
    fn tower_embedding() {
        let mut rng = rng();
        for _ in 0..100 {
            let a = Block32(rng.random());
            let b = Block32(rng.random());

            // 1. Structure check
            let a_lifted: Block64 = a.into();
            let (lo, hi) = a_lifted.split();

            assert_eq!(lo, a, "Embedding structure failed: low part mismatch");
            assert_eq!(
                hi,
                Block32::ZERO,
                "Embedding structure failed: high part must be zero"
            );

            // 2. Addition Homomorphism
            let sum_sub = a + b;
            let sum_lifted: Block64 = sum_sub.into();
            let sum_in_super = Block64::from(a) + Block64::from(b);

            assert_eq!(sum_lifted, sum_in_super, "Homomorphism failed: add");

            // 3. Multiplication Homomorphism
            let prod_sub = a * b;
            let prod_lifted: Block64 = prod_sub.into();
            let prod_in_super = Block64::from(a) * Block64::from(b);

            assert_eq!(prod_lifted, prod_in_super, "Homomorphism failed: mul");
        }
    }

    // ==================================
    // HARDWARE
    // ==================================

    #[test]
    fn isomorphism_roundtrip() {
        let mut rng = rng();
        for _ in 0..1000 {
            let val = Block64(rng.random::<u64>());
            assert_eq!(val.to_hardware().to_tower(), val);
        }
    }

    #[test]
    fn flat_mul_homomorphism() {
        let mut rng = rng();
        for _ in 0..1000 {
            let a = Block64(rng.random());
            let b = Block64(rng.random());

            let expected_flat = (a * b).to_hardware();
            let actual_flat = a.to_hardware() * b.to_hardware();

            assert_eq!(
                actual_flat, expected_flat,
                "Block64 flat multiplication mismatch: (a*b)^H != a^H * b^H"
            );
        }
    }

    #[test]
    fn promote_block8_matches_embedding() {
        for v in 0u16..256 {
            let x = Block8(v as u8);
            assert_eq!(
                <Block64 as FlatPromote<Block8>>::promote_flat(x.to_hardware()),
                Block64::from(x).to_hardware(),
                "Block8 -> Block64 promote mismatch at {v:#04x}"
            );
        }
    }

    #[test]
    fn promote_block16_matches_embedding() {
        for v in 0u32..0x1_0000 {
            let x = Block16(v as u16);
            assert_eq!(
                <Block64 as FlatPromote<Block16>>::promote_flat(x.to_hardware()),
                Block64::from(x).to_hardware(),
                "Block16 -> Block64 promote mismatch at {v:#06x}"
            );
        }
    }

    #[test]
    fn promote_block32_matches_embedding() {
        let mut rng = rng();
        let edges = [Block32::ZERO, Block32::ONE, Block32::TAU, Block32(u32::MAX)];
        let random = (0..100_000).map(|_| Block32(rng.random()));

        for x in edges.into_iter().chain(random) {
            assert_eq!(
                <Block64 as FlatPromote<Block32>>::promote_flat(x.to_hardware()),
                Block64::from(x).to_hardware(),
                "Block32 -> Block64 promote mismatch at {:#010x}",
                x.0
            );
        }
    }

    #[test]
    fn packed_consistency() {
        let mut rng = rng();
        for _ in 0..100 {
            let a_vals = [Block64(rng.random()), Block64(rng.random())];
            let b_vals = [Block64(rng.random()), Block64(rng.random())];

            let a_flat_vals = a_vals.map(|x| x.to_hardware());
            let b_flat_vals = b_vals.map(|x| x.to_hardware());
            let a_packed = Flat::<Block64>::pack(&a_flat_vals);
            let b_packed = Flat::<Block64>::pack(&b_flat_vals);

            // 1. Test SIMD Add (XOR)
            let add_res = Block64::add_hardware_packed(a_packed, b_packed);

            let mut add_out = [Block64::ZERO.to_hardware(); 2];
            Flat::<Block64>::unpack(add_res, &mut add_out);

            assert_eq!(add_out[0], (a_vals[0] + b_vals[0]).to_hardware());
            assert_eq!(add_out[1], (a_vals[1] + b_vals[1]).to_hardware());

            // 2. Test SIMD Mul (Isomorphic/Flat basis)
            let mul_res = Block64::mul_hardware_packed(a_packed, b_packed);

            let mut mul_out = [Block64::ZERO.to_hardware(); 2];
            Flat::<Block64>::unpack(mul_res, &mut mul_out);

            assert_eq!(
                mul_out[0],
                (a_vals[0] * b_vals[0]).to_hardware(),
                "Block64 SIMD mul mismatch at index 0"
            );
            assert_eq!(
                mul_out[1],
                (a_vals[1] * b_vals[1]).to_hardware(),
                "Block64 SIMD mul mismatch at index 1"
            );
        }
    }

    // ==================================
    // PACKED
    // ==================================

    #[test]
    fn pack_unpack_roundtrip() {
        let mut rng = rng();
        let data = [Block64(rng.random()), Block64(rng.random())];

        let packed = Block64::pack(&data);
        let mut unpacked = [Block64::ZERO; 2];

        Block64::unpack(packed, &mut unpacked);
        assert_eq!(data, unpacked);
    }

    #[test]
    fn packed_add_consistency() {
        let mut rng = rng();
        let a_vals = [Block64(rng.random()), Block64(rng.random())];
        let b_vals = [Block64(rng.random()), Block64(rng.random())];

        let res_packed = Block64::pack(&a_vals) + Block64::pack(&b_vals);
        let mut res_unpacked = [Block64::ZERO; 2];
        Block64::unpack(res_packed, &mut res_unpacked);

        assert_eq!(res_unpacked[0], a_vals[0] + b_vals[0]);
        assert_eq!(res_unpacked[1], a_vals[1] + b_vals[1]);
    }

    #[test]
    fn packed_mul_consistency() {
        let mut rng = rng();

        for _ in 0..1000 {
            let mut a_arr = [Block64::ZERO; PACKED_WIDTH_64];
            let mut b_arr = [Block64::ZERO; PACKED_WIDTH_64];

            for i in 0..PACKED_WIDTH_64 {
                let val_a: u64 = rng.random();
                let val_b: u64 = rng.random();
                a_arr[i] = Block64(val_a);
                b_arr[i] = Block64(val_b);
            }

            let a_packed = PackedBlock64(a_arr);
            let b_packed = PackedBlock64(b_arr);

            // Perform SIMD multiplication
            let c_packed = a_packed * b_packed;

            // Verify against Scalar
            let mut c_expected = [Block64::ZERO; PACKED_WIDTH_64];
            for i in 0..PACKED_WIDTH_64 {
                c_expected[i] = a_arr[i] * b_arr[i];
            }

            assert_eq!(c_packed.0, c_expected, "SIMD Block64 mismatch!");
        }
    }

    proptest! {
        #[test]
        fn parity_masks_match_from_hardware(x_flat in any::<u64>()) {
            let tower = Block64::from_hardware(Flat::from_raw(Block64(x_flat))).0;

            for (k, &mask) in FLAT_TO_TOWER_BIT_MASKS_64.iter().enumerate() {
                // Ensure the static masks
                // themselves are correct.
                let parity = ((x_flat & mask).count_ones() & 1) as u8;
                let bit = ((tower >> k) & 1) as u8;
                prop_assert_eq!(parity, bit, "Block64 static mask mismatch at k={}", k);

                // Ensure XOR-tree implementation matches.
                let via_api = Flat::from_raw(Block64(x_flat)).tower_bit(k);
                prop_assert_eq!(
                    via_api, bit,
                    "Block64 tower_bit_from_hardware mismatch at x_flat={:#018x}, bit_idx={}",
                    x_flat, k
                );
            }
        }
    }
}