embeddenator_vsa/

ternary_vec.rs

//! Packed ternary vector representation.
//!
//! This module implements the bitsliced balanced ternary representation that serves as
//! the internal substrate for fast dot/bind/bundle operations.
//!
//! # Bitsliced Design
//!
//! See `docs/BITSLICED_TERNARY_DESIGN.md` for comprehensive documentation.
//!
//! ## Representation: 2 bits per dimension
//! - 00 = Z (0)
//! - 01 = P (+1)
//! - 10 = N (-1)
//! - 11 = unused (treated as Z)
//!
//! ## Key Features
//!
//! - **Word-level parallelism**: 32 trits per u64 word
//! - **Bitplane separation**: Even bits (P plane), odd bits (N plane)
//! - **SIMD-friendly**: Branchless operations, vectorization-ready
//! - **GPU-ready**: Coalesced memory access, no divergence
//!
//! ## Performance
//!
//! - Dot product: 32 trits per word operation
//! - Bind: Pure bitwise operations (AND, OR, shift)
//! - Bundle: Saturating addition via bitwise logic
//! - Memory: 2 bits per trit (optimal for ternary encoding)
//!
//! ## When to Use
//!
//! Use `PackedTritVec` when:
//! - Vector density ≥ 25%
//! - Performing bulk operations (dot, bind, bundle)
//! - SIMD/GPU acceleration is needed
//! - Throughput over latency
//!
//! Use `SparseVec` when:
//! - Vector density < 25%
//! - Random access patterns
//! - Incremental construction
//! - Memory is constrained
//!
//! ## Future Enhancements
//!
//! - Explicit SIMD (AVX2, AVX-512, NEON) - Phase 2
//! - GPU acceleration (CUDA, OpenCL) - Phase 5
//! - See `IMPLEMENTATION_PLAN.md` for roadmap

use crate::ternary::Trit;
use crate::vsa::SparseVec;

// SIMD feature detection for runtime dispatch
#[cfg(target_arch = "x86_64")]
fn has_avx512_vpopcntdq() -> bool {
    // Check for AVX-512F (foundation) + AVX-512VPOPCNTDQ (native popcount)
    std::arch::is_x86_feature_detected!("avx512f")
        && std::arch::is_x86_feature_detected!("avx512vpopcntdq")
}

#[cfg(target_arch = "x86_64")]
fn has_avx2() -> bool {
    std::arch::is_x86_feature_detected!("avx2")
}

#[derive(Clone, Debug, PartialEq, Eq)]
pub struct PackedTritVec {
    len: usize,
    data: Vec<u64>,
}

impl PackedTritVec {
    const MASK_EVEN_BITS: u64 = 0x5555_5555_5555_5555;

    #[inline]
    fn ensure_len_and_clear(&mut self, len: usize) {
        self.len = len;
        let words = Self::word_count_for_len(len);
        if self.data.len() != words {
            self.data.resize(words, 0u64);
        }
        self.data.fill(0u64);
    }

    pub fn new_zero(len: usize) -> Self {
        let bits = len.saturating_mul(2);
        let words = bits.div_ceil(64);
        Self {
            len,
            data: vec![0u64; words],
        }
    }

    #[inline]
    fn word_count_for_len(len: usize) -> usize {
        let bits = len.saturating_mul(2);
        bits.div_ceil(64)
    }

    #[inline]
    fn last_word_mask(len: usize) -> u64 {
        let lanes_in_last = len % 32;
        if lanes_in_last == 0 {
            !0u64
        } else {
            let used_bits = lanes_in_last * 2;
            if used_bits >= 64 {
                !0u64
            } else {
                (1u64 << used_bits) - 1
            }
        }
    }

    pub fn len(&self) -> usize {
        self.len
    }

    pub fn is_empty(&self) -> bool {
        self.len == 0
    }

    #[inline]
    fn word_bit_index(i: usize) -> (usize, usize) {
        let bit = i * 2;
        (bit / 64, bit % 64)
    }

    pub fn get(&self, i: usize) -> Trit {
        if i >= self.len {
            return Trit::Z;
        }
        let (word, bit) = Self::word_bit_index(i);
        let w = self.data.get(word).copied().unwrap_or(0);
        let v = (w >> bit) & 0b11;
        match v {
            0b01 => Trit::P,
            0b10 => Trit::N,
            _ => Trit::Z,
        }
    }

    pub fn set(&mut self, i: usize, t: Trit) {
        if i >= self.len {
            return;
        }
        let (word, bit) = Self::word_bit_index(i);
        if let Some(w) = self.data.get_mut(word) {
            *w &= !(0b11u64 << bit);
            let enc = match t {
                Trit::Z => 0b00u64,
                Trit::P => 0b01u64,
                Trit::N => 0b10u64,
            };
            *w |= enc << bit;
        }
    }

    pub fn from_sparsevec(vec: &SparseVec, len: usize) -> Self {
        let mut out = Self::new_zero(len);
        out.fill_from_sparsevec(vec, len);
        out
    }

    /// Fill this packed vector from a SparseVec, reusing existing allocation.
    ///
    /// This is a hot-path helper for PackedTritVec operations to avoid repeated allocations.
    pub fn fill_from_sparsevec(&mut self, vec: &SparseVec, len: usize) {
        self.ensure_len_and_clear(len);

        // Fast set: output is already zeroed, so we can OR lane encodings.
        // 01 => P (+1) sets even bit; 10 => N (-1) sets odd bit.
        for &idx in &vec.pos {
            if idx < len {
                let bit = idx * 2;
                let word = bit / 64;
                let shift = bit % 64;
                self.data[word] |= 1u64 << shift;
            }
        }

        for &idx in &vec.neg {
            if idx < len {
                let bit = idx * 2;
                let word = bit / 64;
                let shift = bit % 64;
                self.data[word] |= 1u64 << (shift + 1);
            }
        }

        if !self.data.is_empty() {
            let last = self.data.len() - 1;
            self.data[last] &= Self::last_word_mask(self.len);
        }
    }

    pub fn to_sparsevec(&self) -> SparseVec {
        let mut pos: Vec<usize> = Vec::new();
        let mut neg: Vec<usize> = Vec::new();

        // Word-wise extraction: each u64 holds 32 trits (2 bits each).
        for (word_idx, &word_raw) in self.data.iter().enumerate() {
            let mut word = word_raw;
            if word_idx + 1 == self.data.len() {
                word &= Self::last_word_mask(self.len);
            }

            // P lanes have the even bit set; N lanes have the odd bit set.
            // Shift odd bits down to even positions to get per-lane masks.
            let pos_bits = word & Self::MASK_EVEN_BITS;
            let neg_bits = (word >> 1) & Self::MASK_EVEN_BITS;

            // CRITICAL FIX: Detect conflicting trits (0b11 = both P and N set)
            // This can occur if a SparseVec with overlapping pos/neg is converted to PackedTritVec
            // When both bits are set, treat as 0 (cancel out) to maintain invariant
            let conflict_bits = pos_bits & neg_bits;
            let clean_pos = pos_bits & !conflict_bits;
            let clean_neg = neg_bits & !conflict_bits;

            // Extract indices for P lanes.
            let mut m = clean_pos;
            while m != 0 {
                let tz = m.trailing_zeros() as usize;
                let lane = tz / 2;
                let idx = word_idx * 32 + lane;
                if idx < self.len {
                    pos.push(idx);
                }
                m &= m - 1;
            }

            // Extract indices for N lanes.
            let mut n = clean_neg;
            while n != 0 {
                let tz = n.trailing_zeros() as usize;
                let lane = tz / 2;
                let idx = word_idx * 32 + lane;
                if idx < self.len {
                    neg.push(idx);
                }
                n &= n - 1;
            }
        }

        SparseVec { pos, neg }
    }

    /// Sparse ternary dot product: sum over i of a_i * b_i.
    ///
    /// Automatically dispatches to the best available SIMD implementation:
    /// - AVX-512 with VPOPCNTDQ (servers, workstations)
    /// - AVX2 (most modern x86_64, including Intel 14th gen)
    /// - Scalar fallback (auto-vectorizes well with LLVM)
    pub fn dot(&self, other: &Self) -> i32 {
        let n = self.len.min(other.len);
        if n == 0 {
            return 0;
        }

        #[cfg(target_arch = "x86_64")]
        {
            let words = Self::word_count_for_len(n)
                .min(self.data.len())
                .min(other.data.len());

            // AVX-512 with native popcount (best, but rare on consumer CPUs)
            if words >= 16 && has_avx512_vpopcntdq() {
                return unsafe { self.dot_avx512(other, n) };
            }

            // AVX2 (available on 14700K and most modern CPUs)
            if words >= 8 && has_avx2() {
                return unsafe { self.dot_avx2(other, n) };
            }
        }

        // Scalar fallback (auto-vectorizes well)
        self.dot_scalar(other, n)
    }

    /// Scalar implementation of dot product (auto-vectorizes well with LLVM).
    #[inline]
    fn dot_scalar(&self, other: &Self, n: usize) -> i32 {
        let words = Self::word_count_for_len(n)
            .min(self.data.len())
            .min(other.data.len());

        let mut acc: i32 = 0;
        for w in 0..words {
            let mut a = self.data[w];
            let mut b = other.data[w];
            if w + 1 == words {
                let mask = Self::last_word_mask(n);
                a &= mask;
                b &= mask;
            }

            let a_pos = a & Self::MASK_EVEN_BITS;
            let a_neg = (a >> 1) & Self::MASK_EVEN_BITS;
            let b_pos = b & Self::MASK_EVEN_BITS;
            let b_neg = (b >> 1) & Self::MASK_EVEN_BITS;

            let pp = (a_pos & b_pos).count_ones() as i32;
            let nn = (a_neg & b_neg).count_ones() as i32;
            let pn = (a_pos & b_neg).count_ones() as i32;
            let np = (a_neg & b_pos).count_ones() as i32;

            acc += (pp + nn) - (pn + np);
        }

        acc
    }

    /// AVX-512 accelerated dot product using native VPOPCNTDQ.
    ///
    /// Processes 8 × u64 = 256 trits per iteration.
    ///
    /// # Safety
    /// Caller must ensure AVX-512F and AVX-512VPOPCNTDQ are available.
    #[cfg(target_arch = "x86_64")]
    #[target_feature(enable = "avx512f", enable = "avx512vpopcntdq")]
    unsafe fn dot_avx512(&self, other: &Self, n: usize) -> i32 {
        use std::arch::x86_64::*;

        let words = Self::word_count_for_len(n)
            .min(self.data.len())
            .min(other.data.len());

        // Broadcast the even-bits mask to all 8 lanes
        let mask_even = _mm512_set1_epi64(Self::MASK_EVEN_BITS as i64);

        // Accumulators for positive and negative contributions
        let mut acc_pos = _mm512_setzero_si512(); // pp + nn
        let mut acc_neg = _mm512_setzero_si512(); // pn + np

        // Process 8 words (256 trits) at a time
        let chunks = words / 8;
        for chunk in 0..chunks {
            let base = chunk * 8;

            // Load 8 × u64 from each vector
            let va = _mm512_loadu_si512(self.data[base..].as_ptr() as *const __m512i);
            let vb = _mm512_loadu_si512(other.data[base..].as_ptr() as *const __m512i);

            // Separate bitplanes: even bits = P plane, odd bits = N plane
            let a_pos = _mm512_and_si512(va, mask_even);
            let a_neg = _mm512_and_si512(_mm512_srli_epi64(va, 1), mask_even);
            let b_pos = _mm512_and_si512(vb, mask_even);
            let b_neg = _mm512_and_si512(_mm512_srli_epi64(vb, 1), mask_even);

            // Compute intersections
            let pp = _mm512_and_si512(a_pos, b_pos); // Both positive
            let nn = _mm512_and_si512(a_neg, b_neg); // Both negative
            let pn = _mm512_and_si512(a_pos, b_neg); // Opposite signs
            let np = _mm512_and_si512(a_neg, b_pos); // Opposite signs

            // Native popcount on 64-bit lanes (AVX-512 VPOPCNTDQ)
            let pp_cnt = _mm512_popcnt_epi64(pp);
            let nn_cnt = _mm512_popcnt_epi64(nn);
            let pn_cnt = _mm512_popcnt_epi64(pn);
            let np_cnt = _mm512_popcnt_epi64(np);

            // Accumulate: positive contributions (matching signs)
            acc_pos = _mm512_add_epi64(acc_pos, pp_cnt);
            acc_pos = _mm512_add_epi64(acc_pos, nn_cnt);

            // Accumulate: negative contributions (opposing signs)
            acc_neg = _mm512_add_epi64(acc_neg, pn_cnt);
            acc_neg = _mm512_add_epi64(acc_neg, np_cnt);
        }

        // Horizontal reduction: sum all 8 lanes
        let pos_sum = _mm512_reduce_add_epi64(acc_pos);
        let neg_sum = _mm512_reduce_add_epi64(acc_neg);
        let mut acc = (pos_sum - neg_sum) as i32;

        // Handle remaining words with scalar
        let remainder_start = chunks * 8;
        for w in remainder_start..words {
            let mut a = self.data[w];
            let mut b = other.data[w];
            if w + 1 == words {
                let mask = Self::last_word_mask(n);
                a &= mask;
                b &= mask;
            }

            let a_pos = a & Self::MASK_EVEN_BITS;
            let a_neg = (a >> 1) & Self::MASK_EVEN_BITS;
            let b_pos = b & Self::MASK_EVEN_BITS;
            let b_neg = (b >> 1) & Self::MASK_EVEN_BITS;

            let pp = (a_pos & b_pos).count_ones() as i32;
            let nn = (a_neg & b_neg).count_ones() as i32;
            let pn = (a_pos & b_neg).count_ones() as i32;
            let np = (a_neg & b_pos).count_ones() as i32;

            acc += (pp + nn) - (pn + np);
        }

        acc
    }

    /// AVX2 accelerated dot product with popcount emulation.
    ///
    /// Processes 4 × u64 = 128 trits per iteration.
    /// Uses the Wilkes-Wheeler-Gill algorithm for popcount via PSHUFB lookup.
    ///
    /// # Safety
    /// Caller must ensure AVX2 is available.
    #[cfg(target_arch = "x86_64")]
    #[target_feature(enable = "avx2")]
    unsafe fn dot_avx2(&self, other: &Self, n: usize) -> i32 {
        use std::arch::x86_64::*;

        let words = Self::word_count_for_len(n)
            .min(self.data.len())
            .min(other.data.len());

        // Broadcast masks for bitplane separation
        let mask_even = _mm256_set1_epi64x(Self::MASK_EVEN_BITS as i64);

        // Nibble lookup table for popcount (0-15 -> 0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4)
        let popcount_lut = _mm256_setr_epi8(
            0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4, 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2,
            3, 3, 4,
        );
        let low_nibble_mask = _mm256_set1_epi8(0x0F);

        // Accumulators (use 64-bit to avoid overflow)
        let mut acc_pos = _mm256_setzero_si256();
        let mut acc_neg = _mm256_setzero_si256();

        // Process 4 words (128 trits) at a time
        let chunks = words / 4;
        for chunk in 0..chunks {
            let base = chunk * 4;

            // Load 4 × u64 from each vector
            let va = _mm256_loadu_si256(self.data[base..].as_ptr() as *const __m256i);
            let vb = _mm256_loadu_si256(other.data[base..].as_ptr() as *const __m256i);

            // Separate bitplanes: even bits = P plane, odd bits = N plane
            let a_pos = _mm256_and_si256(va, mask_even);
            let a_neg = _mm256_and_si256(_mm256_srli_epi64(va, 1), mask_even);
            let b_pos = _mm256_and_si256(vb, mask_even);
            let b_neg = _mm256_and_si256(_mm256_srli_epi64(vb, 1), mask_even);

            // Compute intersections
            let pp = _mm256_and_si256(a_pos, b_pos);
            let nn = _mm256_and_si256(a_neg, b_neg);
            let pn = _mm256_and_si256(a_pos, b_neg);
            let np = _mm256_and_si256(a_neg, b_pos);

            // Popcount using PSHUFB lookup (Wilkes-Wheeler-Gill algorithm)
            let pp_lo = _mm256_shuffle_epi8(popcount_lut, _mm256_and_si256(pp, low_nibble_mask));
            let pp_hi = _mm256_shuffle_epi8(
                popcount_lut,
                _mm256_and_si256(_mm256_srli_epi16(pp, 4), low_nibble_mask),
            );
            let nn_lo = _mm256_shuffle_epi8(popcount_lut, _mm256_and_si256(nn, low_nibble_mask));
            let nn_hi = _mm256_shuffle_epi8(
                popcount_lut,
                _mm256_and_si256(_mm256_srli_epi16(nn, 4), low_nibble_mask),
            );

            let pn_lo = _mm256_shuffle_epi8(popcount_lut, _mm256_and_si256(pn, low_nibble_mask));
            let pn_hi = _mm256_shuffle_epi8(
                popcount_lut,
                _mm256_and_si256(_mm256_srli_epi16(pn, 4), low_nibble_mask),
            );
            let np_lo = _mm256_shuffle_epi8(popcount_lut, _mm256_and_si256(np, low_nibble_mask));
            let np_hi = _mm256_shuffle_epi8(
                popcount_lut,
                _mm256_and_si256(_mm256_srli_epi16(np, 4), low_nibble_mask),
            );

            // Sum byte counts
            let pos_bytes =
                _mm256_add_epi8(_mm256_add_epi8(pp_lo, pp_hi), _mm256_add_epi8(nn_lo, nn_hi));
            let neg_bytes =
                _mm256_add_epi8(_mm256_add_epi8(pn_lo, pn_hi), _mm256_add_epi8(np_lo, np_hi));

            // Horizontal sum of bytes to 64-bit (SAD against zero)
            let pos_sad = _mm256_sad_epu8(pos_bytes, _mm256_setzero_si256());
            let neg_sad = _mm256_sad_epu8(neg_bytes, _mm256_setzero_si256());

            acc_pos = _mm256_add_epi64(acc_pos, pos_sad);
            acc_neg = _mm256_add_epi64(acc_neg, neg_sad);
        }

        // Horizontal reduction: sum all 4 × 64-bit lanes
        let pos_lo = _mm256_castsi256_si128(acc_pos);
        let pos_hi = _mm256_extracti128_si256(acc_pos, 1);
        let pos_sum128 = _mm_add_epi64(pos_lo, pos_hi);

        let neg_lo = _mm256_castsi256_si128(acc_neg);
        let neg_hi = _mm256_extracti128_si256(acc_neg, 1);
        let neg_sum128 = _mm_add_epi64(neg_lo, neg_hi);

        let pos_final = _mm_extract_epi64(pos_sum128, 0) + _mm_extract_epi64(pos_sum128, 1);
        let neg_final = _mm_extract_epi64(neg_sum128, 0) + _mm_extract_epi64(neg_sum128, 1);

        let mut acc = (pos_final - neg_final) as i32;

        // Handle remaining words with scalar
        let remainder_start = chunks * 4;
        for w in remainder_start..words {
            let mut a = self.data[w];
            let mut b = other.data[w];
            if w + 1 == words {
                let mask = Self::last_word_mask(n);
                a &= mask;
                b &= mask;
            }

            let a_pos = a & Self::MASK_EVEN_BITS;
            let a_neg = (a >> 1) & Self::MASK_EVEN_BITS;
            let b_pos = b & Self::MASK_EVEN_BITS;
            let b_neg = (b >> 1) & Self::MASK_EVEN_BITS;

            let pp = (a_pos & b_pos).count_ones() as i32;
            let nn = (a_neg & b_neg).count_ones() as i32;
            let pn = (a_pos & b_neg).count_ones() as i32;
            let np = (a_neg & b_pos).count_ones() as i32;

            acc += (pp + nn) - (pn + np);
        }

        acc
    }

    /// Element-wise ternary multiplication (bind primitive).
    ///
    /// Automatically dispatches to the best available SIMD implementation:
    /// - AVX-512F (servers, workstations)
    /// - AVX2 (most modern x86_64)
    /// - Scalar fallback
    pub fn bind(&self, other: &Self) -> Self {
        let n = self.len.min(other.len);
        if n == 0 {
            return Self::new_zero(0);
        }

        let words = Self::word_count_for_len(n)
            .min(self.data.len())
            .min(other.data.len());
        let mut out = Self::new_zero(n);

        #[cfg(target_arch = "x86_64")]
        {
            // AVX-512F (processes 8 words = 256 trits per iteration)
            if words >= 8 && std::arch::is_x86_feature_detected!("avx512f") {
                unsafe { self.bind_avx512(other, n, &mut out) };
                return out;
            }

            // AVX2 (processes 4 words = 128 trits per iteration)
            if words >= 4 && has_avx2() {
                unsafe { self.bind_avx2(other, n, &mut out) };
                return out;
            }
        }

        // Scalar fallback
        self.bind_scalar(other, n, &mut out);
        out
    }

    /// Scalar bind implementation
    #[inline]
    fn bind_scalar(&self, other: &Self, n: usize, out: &mut Self) {
        let words = Self::word_count_for_len(n)
            .min(self.data.len())
            .min(other.data.len())
            .min(out.data.len());

        for w in 0..words {
            let mut a = self.data[w];
            let mut b = other.data[w];
            if w + 1 == words {
                let mask = Self::last_word_mask(n);
                a &= mask;
                b &= mask;
            }

            let a_pos = a & Self::MASK_EVEN_BITS;
            let a_neg = (a >> 1) & Self::MASK_EVEN_BITS;
            let b_pos = b & Self::MASK_EVEN_BITS;
            let b_neg = (b >> 1) & Self::MASK_EVEN_BITS;

            let same = (a_pos & b_pos) | (a_neg & b_neg);
            let opp = (a_pos & b_neg) | (a_neg & b_pos);

            out.data[w] = same | (opp << 1);
        }

        // Ensure any unused tail bits stay zero.
        if !out.data.is_empty() {
            let last = out.data.len() - 1;
            out.data[last] &= Self::last_word_mask(out.len);
        }
    }

    /// AVX-512 accelerated bind operation.
    ///
    /// Processes 8 × u64 = 256 trits per iteration using 512-bit registers.
    ///
    /// # Safety
    /// Caller must ensure AVX-512F is available.
    #[cfg(target_arch = "x86_64")]
    #[target_feature(enable = "avx512f")]
    unsafe fn bind_avx512(&self, other: &Self, n: usize, out: &mut Self) {
        use std::arch::x86_64::*;

        let words = Self::word_count_for_len(n)
            .min(self.data.len())
            .min(other.data.len())
            .min(out.data.len());

        let mask_even = _mm512_set1_epi64(Self::MASK_EVEN_BITS as i64);

        // Process 8 words at a time
        let chunks = words / 8;
        for chunk in 0..chunks {
            let base = chunk * 8;

            let va = _mm512_loadu_si512(self.data[base..].as_ptr() as *const __m512i);
            let vb = _mm512_loadu_si512(other.data[base..].as_ptr() as *const __m512i);

            // Separate bitplanes
            let a_pos = _mm512_and_si512(va, mask_even);
            let a_neg = _mm512_and_si512(_mm512_srli_epi64(va, 1), mask_even);
            let b_pos = _mm512_and_si512(vb, mask_even);
            let b_neg = _mm512_and_si512(_mm512_srli_epi64(vb, 1), mask_even);

            // Bind: same signs -> P, opposite signs -> N
            let same = _mm512_or_si512(
                _mm512_and_si512(a_pos, b_pos),
                _mm512_and_si512(a_neg, b_neg),
            );
            let opp = _mm512_or_si512(
                _mm512_and_si512(a_pos, b_neg),
                _mm512_and_si512(a_neg, b_pos),
            );

            let result = _mm512_or_si512(same, _mm512_slli_epi64(opp, 1));
            _mm512_storeu_si512(out.data[base..].as_mut_ptr() as *mut __m512i, result);
        }

        // Handle remaining words with scalar
        let remainder_start = chunks * 8;
        for w in remainder_start..words {
            let mut a = self.data[w];
            let mut b = other.data[w];
            if w + 1 == words {
                let mask = Self::last_word_mask(n);
                a &= mask;
                b &= mask;
            }

            let a_pos = a & Self::MASK_EVEN_BITS;
            let a_neg = (a >> 1) & Self::MASK_EVEN_BITS;
            let b_pos = b & Self::MASK_EVEN_BITS;
            let b_neg = (b >> 1) & Self::MASK_EVEN_BITS;

            let same = (a_pos & b_pos) | (a_neg & b_neg);
            let opp = (a_pos & b_neg) | (a_neg & b_pos);

            out.data[w] = same | (opp << 1);
        }

        if !out.data.is_empty() {
            let last = out.data.len() - 1;
            out.data[last] &= Self::last_word_mask(out.len);
        }
    }

    /// AVX2 accelerated bind operation.
    ///
    /// Processes 4 × u64 = 128 trits per iteration using 256-bit registers.
    ///
    /// # Safety
    /// Caller must ensure AVX2 is available.
    #[cfg(target_arch = "x86_64")]
    #[target_feature(enable = "avx2")]
    unsafe fn bind_avx2(&self, other: &Self, n: usize, out: &mut Self) {
        use std::arch::x86_64::*;

        let words = Self::word_count_for_len(n)
            .min(self.data.len())
            .min(other.data.len())
            .min(out.data.len());

        let mask_even = _mm256_set1_epi64x(Self::MASK_EVEN_BITS as i64);

        // Process 4 words at a time
        let chunks = words / 4;
        for chunk in 0..chunks {
            let base = chunk * 4;

            let va = _mm256_loadu_si256(self.data[base..].as_ptr() as *const __m256i);
            let vb = _mm256_loadu_si256(other.data[base..].as_ptr() as *const __m256i);

            // Separate bitplanes
            let a_pos = _mm256_and_si256(va, mask_even);
            let a_neg = _mm256_and_si256(_mm256_srli_epi64(va, 1), mask_even);
            let b_pos = _mm256_and_si256(vb, mask_even);
            let b_neg = _mm256_and_si256(_mm256_srli_epi64(vb, 1), mask_even);

            // Bind: same signs -> P, opposite signs -> N
            let same = _mm256_or_si256(
                _mm256_and_si256(a_pos, b_pos),
                _mm256_and_si256(a_neg, b_neg),
            );
            let opp = _mm256_or_si256(
                _mm256_and_si256(a_pos, b_neg),
                _mm256_and_si256(a_neg, b_pos),
            );

            let result = _mm256_or_si256(same, _mm256_slli_epi64(opp, 1));
            _mm256_storeu_si256(out.data[base..].as_mut_ptr() as *mut __m256i, result);
        }

        // Handle remaining words with scalar
        let remainder_start = chunks * 4;
        for w in remainder_start..words {
            let mut a = self.data[w];
            let mut b = other.data[w];
            if w + 1 == words {
                let mask = Self::last_word_mask(n);
                a &= mask;
                b &= mask;
            }

            let a_pos = a & Self::MASK_EVEN_BITS;
            let a_neg = (a >> 1) & Self::MASK_EVEN_BITS;
            let b_pos = b & Self::MASK_EVEN_BITS;
            let b_neg = (b >> 1) & Self::MASK_EVEN_BITS;

            let same = (a_pos & b_pos) | (a_neg & b_neg);
            let opp = (a_pos & b_neg) | (a_neg & b_pos);

            out.data[w] = same | (opp << 1);
        }

        if !out.data.is_empty() {
            let last = out.data.len() - 1;
            out.data[last] &= Self::last_word_mask(out.len);
        }
    }

    /// Element-wise ternary multiplication into an existing output buffer.
    pub fn bind_into(&self, other: &Self, out: &mut Self) {
        let n = self.len.min(other.len);
        out.ensure_len_and_clear(n);
        if n == 0 {
            return;
        }

        let words = Self::word_count_for_len(n)
            .min(self.data.len())
            .min(other.data.len())
            .min(out.data.len());

        for w in 0..words {
            let mut a = self.data[w];
            let mut b = other.data[w];
            if w + 1 == words {
                let mask = Self::last_word_mask(n);
                a &= mask;
                b &= mask;
            }

            let a_pos = a & Self::MASK_EVEN_BITS;
            let a_neg = (a >> 1) & Self::MASK_EVEN_BITS;
            let b_pos = b & Self::MASK_EVEN_BITS;
            let b_neg = (b >> 1) & Self::MASK_EVEN_BITS;

            let same = (a_pos & b_pos) | (a_neg & b_neg);
            let opp = (a_pos & b_neg) | (a_neg & b_pos);

            out.data[w] = same | (opp << 1);
        }

        if !out.data.is_empty() {
            let last = out.data.len() - 1;
            out.data[last] &= Self::last_word_mask(out.len);
        }
    }

    /// Element-wise saturating ternary addition (bundle primitive for two vectors).
    ///
    /// Automatically dispatches to the best available SIMD implementation:
    /// - AVX-512F (servers, workstations)
    /// - AVX2 (most modern x86_64)
    /// - Scalar fallback
    pub fn bundle(&self, other: &Self) -> Self {
        let n = self.len.min(other.len);
        if n == 0 {
            return Self::new_zero(0);
        }

        let words = Self::word_count_for_len(n)
            .min(self.data.len())
            .min(other.data.len());
        let mut out = Self::new_zero(n);

        #[cfg(target_arch = "x86_64")]
        {
            // AVX-512F (processes 8 words = 256 trits per iteration)
            if words >= 8 && std::arch::is_x86_feature_detected!("avx512f") {
                unsafe { self.bundle_avx512(other, n, &mut out) };
                return out;
            }

            // AVX2 (processes 4 words = 128 trits per iteration)
            if words >= 4 && has_avx2() {
                unsafe { self.bundle_avx2(other, n, &mut out) };
                return out;
            }
        }

        // Scalar fallback
        self.bundle_scalar(other, n, &mut out);
        out
    }

    /// Scalar bundle implementation
    #[inline]
    fn bundle_scalar(&self, other: &Self, n: usize, out: &mut Self) {
        let words = Self::word_count_for_len(n)
            .min(self.data.len())
            .min(other.data.len())
            .min(out.data.len());

        for w in 0..words {
            let mut a = self.data[w];
            let mut b = other.data[w];
            if w + 1 == words {
                let mask = Self::last_word_mask(n);
                a &= mask;
                b &= mask;
            }

            let a_pos = a & Self::MASK_EVEN_BITS;
            let a_neg = (a >> 1) & Self::MASK_EVEN_BITS;
            let b_pos = b & Self::MASK_EVEN_BITS;
            let b_neg = (b >> 1) & Self::MASK_EVEN_BITS;

            let mask = Self::MASK_EVEN_BITS;
            let not_b_neg = (!b_neg) & mask;
            let not_a_neg = (!a_neg) & mask;
            let not_b_pos = (!b_pos) & mask;
            let not_a_pos = (!a_pos) & mask;

            let pos = (a_pos & not_b_neg) | (b_pos & not_a_neg);
            let neg = (a_neg & not_b_pos) | (b_neg & not_a_pos);

            out.data[w] = pos | (neg << 1);
        }

        if !out.data.is_empty() {
            let last = out.data.len() - 1;
            out.data[last] &= Self::last_word_mask(out.len);
        }
    }

    /// AVX-512 accelerated bundle operation.
    ///
    /// Processes 8 × u64 = 256 trits per iteration using 512-bit registers.
    ///
    /// # Safety
    /// Caller must ensure AVX-512F is available.
    #[cfg(target_arch = "x86_64")]
    #[target_feature(enable = "avx512f")]
    unsafe fn bundle_avx512(&self, other: &Self, n: usize, out: &mut Self) {
        use std::arch::x86_64::*;

        let words = Self::word_count_for_len(n)
            .min(self.data.len())
            .min(other.data.len())
            .min(out.data.len());

        let mask_even = _mm512_set1_epi64(Self::MASK_EVEN_BITS as i64);

        // Process 8 words at a time
        let chunks = words / 8;
        for chunk in 0..chunks {
            let base = chunk * 8;

            let va = _mm512_loadu_si512(self.data[base..].as_ptr() as *const __m512i);
            let vb = _mm512_loadu_si512(other.data[base..].as_ptr() as *const __m512i);

            // Separate bitplanes
            let a_pos = _mm512_and_si512(va, mask_even);
            let a_neg = _mm512_and_si512(_mm512_srli_epi64(va, 1), mask_even);
            let b_pos = _mm512_and_si512(vb, mask_even);
            let b_neg = _mm512_and_si512(_mm512_srli_epi64(vb, 1), mask_even);

            // Bundle: saturating add
            // pos = (a_pos & !b_neg) | (b_pos & !a_neg)
            // neg = (a_neg & !b_pos) | (b_neg & !a_pos)
            let not_b_neg = _mm512_andnot_si512(b_neg, mask_even);
            let not_a_neg = _mm512_andnot_si512(a_neg, mask_even);
            let not_b_pos = _mm512_andnot_si512(b_pos, mask_even);
            let not_a_pos = _mm512_andnot_si512(a_pos, mask_even);

            let pos = _mm512_or_si512(
                _mm512_and_si512(a_pos, not_b_neg),
                _mm512_and_si512(b_pos, not_a_neg),
            );
            let neg = _mm512_or_si512(
                _mm512_and_si512(a_neg, not_b_pos),
                _mm512_and_si512(b_neg, not_a_pos),
            );

            let result = _mm512_or_si512(pos, _mm512_slli_epi64(neg, 1));
            _mm512_storeu_si512(out.data[base..].as_mut_ptr() as *mut __m512i, result);
        }

        // Handle remaining words with scalar
        let remainder_start = chunks * 8;
        for w in remainder_start..words {
            let mut a = self.data[w];
            let mut b = other.data[w];
            if w + 1 == words {
                let mask = Self::last_word_mask(n);
                a &= mask;
                b &= mask;
            }

            let a_pos = a & Self::MASK_EVEN_BITS;
            let a_neg = (a >> 1) & Self::MASK_EVEN_BITS;
            let b_pos = b & Self::MASK_EVEN_BITS;
            let b_neg = (b >> 1) & Self::MASK_EVEN_BITS;

            let mask = Self::MASK_EVEN_BITS;
            let not_b_neg = (!b_neg) & mask;
            let not_a_neg = (!a_neg) & mask;
            let not_b_pos = (!b_pos) & mask;
            let not_a_pos = (!a_pos) & mask;

            let pos = (a_pos & not_b_neg) | (b_pos & not_a_neg);
            let neg = (a_neg & not_b_pos) | (b_neg & not_a_pos);

            out.data[w] = pos | (neg << 1);
        }

        if !out.data.is_empty() {
            let last = out.data.len() - 1;
            out.data[last] &= Self::last_word_mask(out.len);
        }
    }

    /// AVX2 accelerated bundle operation.
    ///
    /// Processes 4 × u64 = 128 trits per iteration using 256-bit registers.
    ///
    /// # Safety
    /// Caller must ensure AVX2 is available.
    #[cfg(target_arch = "x86_64")]
    #[target_feature(enable = "avx2")]
    unsafe fn bundle_avx2(&self, other: &Self, n: usize, out: &mut Self) {
        use std::arch::x86_64::*;

        let words = Self::word_count_for_len(n)
            .min(self.data.len())
            .min(other.data.len())
            .min(out.data.len());

        let mask_even = _mm256_set1_epi64x(Self::MASK_EVEN_BITS as i64);

        // Process 4 words at a time
        let chunks = words / 4;
        for chunk in 0..chunks {
            let base = chunk * 4;

            let va = _mm256_loadu_si256(self.data[base..].as_ptr() as *const __m256i);
            let vb = _mm256_loadu_si256(other.data[base..].as_ptr() as *const __m256i);

            // Separate bitplanes
            let a_pos = _mm256_and_si256(va, mask_even);
            let a_neg = _mm256_and_si256(_mm256_srli_epi64(va, 1), mask_even);
            let b_pos = _mm256_and_si256(vb, mask_even);
            let b_neg = _mm256_and_si256(_mm256_srli_epi64(vb, 1), mask_even);

            // Bundle: saturating add
            let not_b_neg = _mm256_andnot_si256(b_neg, mask_even);
            let not_a_neg = _mm256_andnot_si256(a_neg, mask_even);
            let not_b_pos = _mm256_andnot_si256(b_pos, mask_even);
            let not_a_pos = _mm256_andnot_si256(a_pos, mask_even);

            let pos = _mm256_or_si256(
                _mm256_and_si256(a_pos, not_b_neg),
                _mm256_and_si256(b_pos, not_a_neg),
            );
            let neg = _mm256_or_si256(
                _mm256_and_si256(a_neg, not_b_pos),
                _mm256_and_si256(b_neg, not_a_pos),
            );

            let result = _mm256_or_si256(pos, _mm256_slli_epi64(neg, 1));
            _mm256_storeu_si256(out.data[base..].as_mut_ptr() as *mut __m256i, result);
        }

        // Handle remaining words with scalar
        let remainder_start = chunks * 4;
        for w in remainder_start..words {
            let mut a = self.data[w];
            let mut b = other.data[w];
            if w + 1 == words {
                let mask = Self::last_word_mask(n);
                a &= mask;
                b &= mask;
            }

            let a_pos = a & Self::MASK_EVEN_BITS;
            let a_neg = (a >> 1) & Self::MASK_EVEN_BITS;
            let b_pos = b & Self::MASK_EVEN_BITS;
            let b_neg = (b >> 1) & Self::MASK_EVEN_BITS;

            let mask = Self::MASK_EVEN_BITS;
            let not_b_neg = (!b_neg) & mask;
            let not_a_neg = (!a_neg) & mask;
            let not_b_pos = (!b_pos) & mask;
            let not_a_pos = (!a_pos) & mask;

            let pos = (a_pos & not_b_neg) | (b_pos & not_a_neg);
            let neg = (a_neg & not_b_pos) | (b_neg & not_a_pos);

            out.data[w] = pos | (neg << 1);
        }

        if !out.data.is_empty() {
            let last = out.data.len() - 1;
            out.data[last] &= Self::last_word_mask(out.len);
        }
    }

    /// Element-wise saturating ternary addition into an existing output buffer.
    pub fn bundle_into(&self, other: &Self, out: &mut Self) {
        let n = self.len.min(other.len);
        out.ensure_len_and_clear(n);
        if n == 0 {
            return;
        }

        let words = Self::word_count_for_len(n)
            .min(self.data.len())
            .min(other.data.len())
            .min(out.data.len());

        for w in 0..words {
            let mut a = self.data[w];
            let mut b = other.data[w];
            if w + 1 == words {
                let mask = Self::last_word_mask(n);
                a &= mask;
                b &= mask;
            }

            let a_pos = a & Self::MASK_EVEN_BITS;
            let a_neg = (a >> 1) & Self::MASK_EVEN_BITS;
            let b_pos = b & Self::MASK_EVEN_BITS;
            let b_neg = (b >> 1) & Self::MASK_EVEN_BITS;

            let mask = Self::MASK_EVEN_BITS;
            let not_b_neg = (!b_neg) & mask;
            let not_a_neg = (!a_neg) & mask;
            let not_b_pos = (!b_pos) & mask;
            let not_a_pos = (!a_pos) & mask;

            let pos = (a_pos & not_b_neg) | (b_pos & not_a_neg);
            let neg = (a_neg & not_b_pos) | (b_neg & not_a_pos);

            out.data[w] = pos | (neg << 1);
        }

        if !out.data.is_empty() {
            let last = out.data.len() - 1;
            out.data[last] &= Self::last_word_mask(out.len);
        }
    }
}