amari_core/

verified.rs

//! Formally verified geometric algebra with phantom types and Creusot annotations
//!
//! This module provides type-safe, mathematically verified implementations of
//! geometric algebra operations using phantom types for compile-time invariants
//! and Creusot annotations for formal verification.

use num_traits::{Float, One, Zero};
use std::marker::PhantomData;

#[cfg(feature = "formal-verification")]
use creusot_contracts::macros::{ensures, requires};

/// Phantom type encoding the metric signature of a Clifford algebra Cl(p,q,r)
/// - P: number of positive basis vectors (e_i² = +1)
/// - Q: number of negative basis vectors (e_i² = -1)
/// - R: number of null basis vectors (e_i² = 0)
pub struct Signature<const P: usize, const Q: usize, const R: usize>;

/// Phantom type for compile-time grade tracking
pub struct Grade<const G: usize>;

/// Type-level dimension marker for vectors
pub struct Dim<const D: usize>;

/// A verified multivector with compile-time signature guarantees
///
/// This structure ensures at the type level that:
/// 1. The signature is fixed and consistent
/// 2. Operations preserve algebraic laws
/// 3. Dimension bounds are respected
#[derive(Debug, Clone)]
pub struct VerifiedMultivector<T, const P: usize, const Q: usize, const R: usize>
where
    T: Float + Zero + One,
{
    /// Coefficients in lexicographic basis blade order
    pub(crate) coefficients: Vec<T>,
    /// Phantom marker for signature
    _signature: PhantomData<Signature<P, Q, R>>,
}

impl<T, const P: usize, const Q: usize, const R: usize> VerifiedMultivector<T, P, Q, R>
where
    T: Float + Zero + One,
{
    /// The total dimension of the algebra
    pub const DIM: usize = P + Q + R;

    /// The number of basis blades (2^n)
    pub const BASIS_SIZE: usize = 1 << (P + Q + R);

    /// Create a new verified multivector from coefficients
    ///
    /// # Type Invariants
    /// - Coefficients array must have exactly 2^(P+Q+R) elements
    /// - Signature is encoded at type level and cannot be violated
    #[cfg_attr(feature = "formal-verification",
        requires(coefficients.len() == Self::BASIS_SIZE))]
    pub fn new(coefficients: Vec<T>) -> Result<Self, &'static str> {
        if coefficients.len() != Self::BASIS_SIZE {
            return Err("Coefficient array size must equal 2^(P+Q+R)");
        }

        Ok(Self {
            coefficients,
            _signature: PhantomData,
        })
    }

    /// Create a scalar multivector
    #[cfg_attr(feature = "formal-verification",
        ensures(result.is_scalar()),
        ensures(result.coefficients[0] == value),
        ensures(result.coefficients.len() == Self::BASIS_SIZE))]
    pub fn scalar(value: T) -> Self {
        let mut coefficients = vec![T::zero(); Self::BASIS_SIZE];
        coefficients[0] = value;

        Self {
            coefficients,
            _signature: PhantomData,
        }
    }

    /// Create a basis vector e_i
    ///
    /// # Type Safety
    /// The index is bounds-checked against the signature dimensions
    #[cfg_attr(feature = "formal-verification",
        requires(index < Self::DIM),
        ensures(result.grade() == 1))]
    pub fn basis_vector(index: usize) -> Result<Self, &'static str> {
        if index >= Self::DIM {
            return Err("Basis vector index exceeds dimension");
        }

        let mut coefficients = vec![T::zero(); Self::BASIS_SIZE];
        coefficients[1 << index] = T::one();

        Ok(Self {
            coefficients,
            _signature: PhantomData,
        })
    }

    /// Check if this is a pure scalar
    #[cfg_attr(feature = "formal-verification",
        ensures(result == (self.coefficients[1..].iter().all(|c| c.is_zero()))))]
    pub fn is_scalar(&self) -> bool {
        self.coefficients[1..].iter().all(|c| c.is_zero())
    }

    /// Compute the grade (highest non-zero grade component)
    #[cfg_attr(feature = "formal-verification",
        ensures(result <= Self::DIM))]
    pub fn grade(&self) -> usize {
        // Count bits to determine grade of each basis blade
        for (i, coeff) in self.coefficients.iter().enumerate().rev() {
            if !coeff.is_zero() {
                return i.count_ones() as usize;
            }
        }
        0
    }

    /// Addition of multivectors (same signature enforced by types)
    #[cfg_attr(feature = "formal-verification",
        ensures(result.coefficients.len() == self.coefficients.len()))]
    pub fn add(&self, other: &Self) -> Self {
        let coefficients: Vec<T> = self
            .coefficients
            .iter()
            .zip(&other.coefficients)
            .map(|(a, b)| *a + *b)
            .collect();

        Self {
            coefficients,
            _signature: PhantomData,
        }
    }

    /// Geometric product with compile-time signature matching
    ///
    /// # Mathematical Properties (verified by Creusot when enabled)
    /// - Associativity: (AB)C = A(BC)
    /// - Distributivity: A(B+C) = AB + AC
    /// - Identity: 1*A = A*1 = A
    #[cfg_attr(feature = "formal-verification",
        requires(self.coefficients.len() == Self::BASIS_SIZE),
        requires(other.coefficients.len() == Self::BASIS_SIZE),
        ensures(result.coefficients.len() == Self::BASIS_SIZE))]
    pub fn geometric_product(&self, other: &Self) -> Self {
        // Simplified geometric product implementation
        // In practice, this would use the Cayley table for the signature
        let mut result = vec![T::zero(); Self::BASIS_SIZE];

        for i in 0..Self::BASIS_SIZE {
            for j in 0..Self::BASIS_SIZE {
                let sign = self.compute_product_sign(i, j);
                let target_index = i ^ j; // XOR gives the resulting basis blade
                result[target_index] = result[target_index]
                    + self.coefficients[i] * other.coefficients[j] * T::from(sign).unwrap();
            }
        }

        Self {
            coefficients: result,
            _signature: PhantomData,
        }
    }

    /// Compute the sign (and zero factor) of the product of two basis blades
    ///
    /// Accounts for the metric signature (P,Q,R):
    /// - Basis vectors 0..P square to +1
    /// - Basis vectors P..P+Q square to -1
    /// - Basis vectors P+Q..P+Q+R square to 0
    fn compute_product_sign(&self, blade1: usize, blade2: usize) -> i32 {
        // Step 1: Count transposition swaps for reorder sign
        let swaps = self.count_swaps(blade1, blade2);
        let mut sign: i32 = if swaps.is_multiple_of(2) { 1 } else { -1 };

        // Step 2: Apply metric signature for shared basis vectors
        // When e_i appears in both blades, e_i * e_i = metric(i)
        let shared = blade1 & blade2;
        for i in 0..Self::DIM {
            if shared & (1 << i) != 0 {
                if i >= P + Q {
                    // Null basis vector (R range): e_i² = 0
                    return 0;
                } else if i >= P {
                    // Negative basis vector (Q range): e_i² = -1
                    sign = -sign;
                }
                // Positive basis vector (P range): e_i² = +1, no change
            }
        }

        sign
    }

    /// Count the number of swaps needed to reorder basis vectors
    fn count_swaps(&self, blade1: usize, blade2: usize) -> usize {
        // Count inversions when combining blade indices
        let mut count = 0;
        for i in 0..Self::DIM {
            if blade1 & (1 << i) != 0 {
                for j in 0..i {
                    if blade2 & (1 << j) != 0 {
                        count += 1;
                    }
                }
            }
        }
        count
    }
}

/// Type-safe k-vector (homogeneous grade element)
///
/// This type guarantees at compile-time that the multivector
/// contains only elements of grade K
pub struct KVector<T, const K: usize, const P: usize, const Q: usize, const R: usize>
where
    T: Float + Zero + One,
{
    multivector: VerifiedMultivector<T, P, Q, R>,
    _grade: PhantomData<Grade<K>>,
}

impl<T, const K: usize, const P: usize, const Q: usize, const R: usize> KVector<T, K, P, Q, R>
where
    T: Float + Zero + One,
{
    /// Create a k-vector from a general multivector by grade projection
    #[cfg_attr(feature = "formal-verification",
        requires(K <= P + Q + R),
        ensures(result.grade() == K))]
    pub fn from_multivector(mv: VerifiedMultivector<T, P, Q, R>) -> Self {
        let mut coefficients = vec![T::zero(); VerifiedMultivector::<T, P, Q, R>::BASIS_SIZE];

        // Extract only grade-K components
        for (i, coeff) in mv.coefficients.iter().enumerate() {
            if i.count_ones() as usize == K {
                coefficients[i] = *coeff;
            }
        }

        Self {
            multivector: VerifiedMultivector {
                coefficients,
                _signature: PhantomData,
            },
            _grade: PhantomData,
        }
    }

    /// Get the grade (always returns K due to type constraint)
    pub const fn grade(&self) -> usize {
        K
    }

    /// Inner product with another k-vector (produces a scalar)
    pub fn inner_product(&self, other: &Self) -> T {
        // Inner product of k-vectors of same grade
        self.multivector
            .coefficients
            .iter()
            .zip(&other.multivector.coefficients)
            .map(|(a, b)| *a * *b)
            .fold(T::zero(), |acc, x| acc + x)
    }
}

// Trait for outer product with specific grade combinations
// This works around const generic arithmetic limitations
pub trait OuterProduct<T, const J: usize, const P: usize, const Q: usize, const R: usize>
where
    T: Float + Zero + One,
{
    type Output;
    fn outer_product(&self, other: &KVector<T, J, P, Q, R>) -> Self::Output;
}

// Implement outer product for specific grade combinations
// Grade 1 ∧ Grade 1 = Grade 2 (vector ∧ vector = bivector)
impl<T, const P: usize, const Q: usize, const R: usize> OuterProduct<T, 1, P, Q, R>
    for KVector<T, 1, P, Q, R>
where
    T: Float + Zero + One,
{
    type Output = KVector<T, 2, P, Q, R>;

    fn outer_product(&self, other: &KVector<T, 1, P, Q, R>) -> Self::Output {
        let mut coefficients = vec![T::zero(); VerifiedMultivector::<T, P, Q, R>::BASIS_SIZE];

        // Compute outer product for grade-1 vectors
        for i in 0..VerifiedMultivector::<T, P, Q, R>::DIM {
            for j in i + 1..VerifiedMultivector::<T, P, Q, R>::DIM {
                let blade_i = 1 << i;
                let blade_j = 1 << j;
                let blade_ij = blade_i | blade_j;

                coefficients[blade_ij] = self.multivector.coefficients[blade_i]
                    * other.multivector.coefficients[blade_j]
                    - self.multivector.coefficients[blade_j]
                        * other.multivector.coefficients[blade_i];
            }
        }

        KVector {
            multivector: VerifiedMultivector {
                coefficients,
                _signature: PhantomData,
            },
            _grade: PhantomData,
        }
    }
}

// Grade 1 ∧ Grade 2 = Grade 3 (vector ∧ bivector = trivector)
impl<T, const P: usize, const Q: usize, const R: usize> OuterProduct<T, 2, P, Q, R>
    for KVector<T, 1, P, Q, R>
where
    T: Float + Zero + One,
{
    type Output = KVector<T, 3, P, Q, R>;

    fn outer_product(&self, other: &KVector<T, 2, P, Q, R>) -> Self::Output {
        let dim = VerifiedMultivector::<T, P, Q, R>::DIM;
        let basis_size = VerifiedMultivector::<T, P, Q, R>::BASIS_SIZE;
        let mut coefficients = vec![T::zero(); basis_size];

        // For each grade-1 blade (vectors) and grade-2 blade (bivectors),
        // compute their wedge product contributing to grade-3 (trivectors)
        for i in 0..dim {
            let blade_i = 1usize << i;
            for j in 0..basis_size {
                if j.count_ones() != 2 {
                    continue;
                }
                // Only wedge if basis vectors don't overlap
                if blade_i & j == 0 {
                    let target = blade_i | j;
                    // Compute sign from reordering e_i into canonical position within target
                    let mut swaps = 0;
                    for k in 0..i {
                        if j & (1 << k) != 0 {
                            swaps += 1;
                        }
                    }
                    let sign = if swaps % 2 == 0 { T::one() } else { -T::one() };
                    coefficients[target] = coefficients[target]
                        + sign
                            * self.multivector.coefficients[blade_i]
                            * other.multivector.coefficients[j];
                }
            }
        }

        KVector {
            multivector: VerifiedMultivector {
                coefficients,
                _signature: PhantomData,
            },
            _grade: PhantomData,
        }
    }
}

// Grade 2 ∧ Grade 1 = Grade 3 (bivector ∧ vector = trivector)
impl<T, const P: usize, const Q: usize, const R: usize> OuterProduct<T, 1, P, Q, R>
    for KVector<T, 2, P, Q, R>
where
    T: Float + Zero + One,
{
    type Output = KVector<T, 3, P, Q, R>;

    fn outer_product(&self, other: &KVector<T, 1, P, Q, R>) -> Self::Output {
        let dim = VerifiedMultivector::<T, P, Q, R>::DIM;
        let basis_size = VerifiedMultivector::<T, P, Q, R>::BASIS_SIZE;
        let mut coefficients = vec![T::zero(); basis_size];

        // For each grade-2 blade (bivectors) and grade-1 blade (vectors),
        // compute their wedge product contributing to grade-3 (trivectors)
        for i in 0..basis_size {
            if i.count_ones() != 2 {
                continue;
            }
            for j in 0..dim {
                let blade_j = 1usize << j;
                // Only wedge if basis vectors don't overlap
                if i & blade_j == 0 {
                    let target = i | blade_j;
                    // Compute sign: count how many bits in i are above j
                    let mut swaps = 0;
                    for k in (j + 1)..dim {
                        if i & (1 << k) != 0 {
                            swaps += 1;
                        }
                    }
                    let sign = if swaps % 2 == 0 { T::one() } else { -T::one() };
                    coefficients[target] = coefficients[target]
                        + sign
                            * self.multivector.coefficients[i]
                            * other.multivector.coefficients[blade_j];
                }
            }
        }

        KVector {
            multivector: VerifiedMultivector {
                coefficients,
                _signature: PhantomData,
            },
            _grade: PhantomData,
        }
    }
}

/// A verified rotor with compile-time guarantees
///
/// Rotors are even-grade multivectors with unit norm that
/// represent rotations in the geometric algebra
pub struct VerifiedRotor<T, const P: usize, const Q: usize, const R: usize>
where
    T: Float + Zero + One,
{
    multivector: VerifiedMultivector<T, P, Q, R>,
    /// Phantom data ensuring this is a valid rotor
    _rotor_invariant: PhantomData<()>,
}

impl<T, const P: usize, const Q: usize, const R: usize> VerifiedRotor<T, P, Q, R>
where
    T: Float + Zero + One,
{
    /// Create a rotor from a multivector
    ///
    /// # Invariants (verified at runtime, would be proven with Creusot)
    /// - Must be even grade
    /// - Must have unit norm
    #[cfg_attr(feature = "formal-verification",
        requires(mv.is_even_grade()),
        requires((mv.norm() - T::one()).abs() < T::from(0.0001).unwrap()),
        ensures(result.is_ok()))]
    pub fn new(mv: VerifiedMultivector<T, P, Q, R>) -> Result<Self, &'static str> {
        if !Self::is_even_grade(&mv) {
            return Err("Rotor must have even grade");
        }

        let norm = Self::compute_norm(&mv);
        if (norm - T::one()).abs() > T::from(0.0001).unwrap() {
            return Err("Rotor must have unit norm");
        }

        Ok(Self {
            multivector: mv,
            _rotor_invariant: PhantomData,
        })
    }

    /// Check if multivector has only even grade components
    fn is_even_grade(mv: &VerifiedMultivector<T, P, Q, R>) -> bool {
        for (i, coeff) in mv.coefficients.iter().enumerate() {
            let grade = i.count_ones() as usize;
            if !grade.is_multiple_of(2) && !coeff.is_zero() {
                return false;
            }
        }
        true
    }

    /// Compute the norm of a multivector
    fn compute_norm(mv: &VerifiedMultivector<T, P, Q, R>) -> T {
        mv.coefficients
            .iter()
            .map(|c| *c * *c)
            .fold(T::zero(), |acc, x| acc + x)
            .sqrt()
    }

    /// Compose two rotors (rotation composition)
    ///
    /// # Type Safety
    /// - Result is guaranteed to be a valid rotor
    /// - Unit norm is preserved
    #[cfg_attr(feature = "formal-verification",
        ensures(Self::compute_norm(&result.multivector) == T::one()))]
    pub fn compose(&self, other: &Self) -> Self {
        let composed = self.multivector.geometric_product(&other.multivector);

        // Normalize to ensure unit norm (should already be close due to properties)
        let norm = Self::compute_norm(&composed);
        let normalized = VerifiedMultivector {
            coefficients: composed.coefficients.iter().map(|c| *c / norm).collect(),
            _signature: PhantomData,
        };

        Self {
            multivector: normalized,
            _rotor_invariant: PhantomData,
        }
    }

    /// Apply rotor to a multivector (rotation/reflection)
    ///
    /// Computes R v R† where R† is the reverse of the rotor.
    pub fn apply_to_multivector(
        &self,
        v: &VerifiedMultivector<T, P, Q, R>,
    ) -> VerifiedMultivector<T, P, Q, R> {
        // Compute R† (reverse): negate grades where k*(k-1)/2 is odd
        let mut rev_coeffs = self.multivector.coefficients.clone();
        for (i, coeff) in rev_coeffs.iter_mut().enumerate() {
            let grade = i.count_ones() as usize;
            if grade >= 2 && (grade * (grade - 1) / 2) % 2 == 1 {
                *coeff = -*coeff;
            }
        }
        let r_rev = VerifiedMultivector::<T, P, Q, R> {
            coefficients: rev_coeffs,
            _signature: PhantomData,
        };

        // Compute R * v * R†
        let rv = self.multivector.geometric_product(v);
        rv.geometric_product(&r_rev)
    }
}

/// Type-safe vector with compile-time dimension
pub struct Vector<T, const D: usize>
where
    T: Float + Zero + One,
{
    pub(crate) data: Vec<T>,
    _dim: PhantomData<Dim<D>>,
}

impl<T, const D: usize> Vector<T, D>
where
    T: Float + Zero + One,
{
    /// Create a new vector (dimension checked at compile time)
    #[cfg_attr(feature = "formal-verification",
        requires(data.len() == D),
        ensures(result.dimension() == D))]
    pub fn new(data: Vec<T>) -> Result<Self, &'static str> {
        if data.len() != D {
            return Err("Vector data length must match dimension");
        }

        Ok(Self {
            data,
            _dim: PhantomData,
        })
    }

    /// Get the dimension (always D due to type parameter)
    pub const fn dimension(&self) -> usize {
        D
    }

    /// Dot product (only defined for same dimension - enforced by types)
    #[cfg_attr(feature = "formal-verification",
        ensures(result >= T::zero()))] // For positive-definite metrics
    pub fn dot(&self, other: &Self) -> T {
        self.data
            .iter()
            .zip(&other.data)
            .map(|(a, b)| *a * *b)
            .fold(T::zero(), |acc, x| acc + x)
    }

    /// Add two vectors (same dimension enforced at compile time)
    pub fn add(&self, other: &Self) -> Self {
        Self {
            data: self
                .data
                .iter()
                .zip(&other.data)
                .map(|(a, b)| *a + *b)
                .collect(),
            _dim: PhantomData,
        }
    }
}

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

    #[test]
    fn test_phantom_type_dimension_safety() {
        // This would fail to compile if dimensions don't match
        let v1 = Vector::<f64, 3>::new(vec![1.0, 2.0, 3.0]).unwrap();
        let v2 = Vector::<f64, 3>::new(vec![4.0, 5.0, 6.0]).unwrap();
        let _v3 = v1.add(&v2); // OK - same dimensions

        // The following would not compile:
        // let v4 = Vector::<f64, 2>::new(vec![1.0, 2.0]).unwrap();
        // let v5 = v1.add(&v4); // ERROR: type mismatch
    }

    #[test]
    fn test_signature_type_safety() {
        // Clifford algebra Cl(3,0,0) - 3D Euclidean space
        let mv1 = VerifiedMultivector::<f64, 3, 0, 0>::scalar(2.0);
        let mv2 = VerifiedMultivector::<f64, 3, 0, 0>::scalar(3.0);
        let _mv3 = mv1.add(&mv2); // OK - same signature

        // The following would not compile:
        // let mv4 = VerifiedMultivector::<f64, 2, 1, 0>::scalar(1.0); // Different signature
        // let mv5 = mv1.add(&mv4); // ERROR: type mismatch
    }

    #[test]
    fn test_grade_preservation() {
        // Create a bivector (grade 2)
        let bivector =
            KVector::<f64, 2, 3, 0, 0>::from_multivector(VerifiedMultivector::scalar(1.0));

        assert_eq!(bivector.grade(), 2);
        // Grade is guaranteed at compile time
    }

    #[test]
    fn test_outer_product_type_safety() {
        use super::OuterProduct;

        // Create two vectors (grade 1) in 3D Euclidean space
        let v1 = KVector::<f64, 1, 3, 0, 0>::from_multivector(
            VerifiedMultivector::basis_vector(0).unwrap(),
        );
        let v2 = KVector::<f64, 1, 3, 0, 0>::from_multivector(
            VerifiedMultivector::basis_vector(1).unwrap(),
        );

        // Outer product of two vectors gives a bivector (grade 2)
        let bivector: KVector<f64, 2, 3, 0, 0> = v1.outer_product(&v2);

        // The type system guarantees this is grade 2
        assert_eq!(bivector.grade(), 2);

        // The following would not compile due to type mismatch:
        // let wrong: KVector<f64, 3, 3, 0, 0> = v1.outer_product(&v2);
        // Error: expected KVector<_, 3, _, _, _>, found KVector<_, 2, _, _, _>
    }
}