geoit 0.0.2

//! Multivector inversion.
//!
//! Three paths, tried in order:
//! 1. **Blade inverse**: if Mv is a single blade, `inv(B) = B / B²`.
//! 2. **Versor inverse**: if `M * rev(M)` is a nonzero scalar, `inv(M) = rev(M) / (M * rev(M))`.
//! 3. **Even/odd decomposition**: split `M = M⁺ + M⁻`, solve `M * X = 1` via
//!    conjugate structure. Works for all invertible Mvs in non-degenerate algebras.
//!
//! Returns `Err(InverseError)` for non-invertible Mvs.

use crate::algebra::mv::Mv;
use crate::algebra::ops;
use crate::algebra::signature::Signature;
use crate::scalar::Scalar;

/// Error from Mv inversion.
#[derive(Clone, Debug)]
pub enum InverseError {
    /// The Mv is zero.
    ZeroMv,
    /// The Mv has zero norm squared (null element — e.g., CGA null vector).
    NullElement,
    /// The Mv is not invertible in this algebra (degenerate case).
    NotInvertible,
    /// Iterative inversion did not converge.
    DidNotConverge { iterations: usize },
}

impl std::fmt::Display for InverseError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            InverseError::ZeroMv => write!(f, "cannot invert zero multivector"),
            InverseError::NullElement => write!(f, "null element (norm² = 0) is not invertible"),
            InverseError::NotInvertible => {
                write!(f, "multivector is not invertible in this algebra")
            }
            InverseError::DidNotConverge { iterations } => write!(
                f,
                "iterative inversion did not converge after {} iterations",
                iterations
            ),
        }
    }
}

/// Check if an Mv is a versor: a product of invertible vectors.
///
/// Test: `M * rev(M)` must be a nonzero scalar.
/// This is necessary but not sufficient in general, but covers
/// all versors (rotors, reflectors, etc.) which are the practical case.
pub fn is_versor(m: &Mv, sig: &Signature) -> bool {
    if m.is_zero() {
        return false;
    }
    let m_rev = ops::reverse(m);
    let product = ops::geometric(m, &m_rev, sig);
    // Must be a pure nonzero scalar
    if product.len() != 1 {
        return false;
    }
    let scalar_coeff = product.coefficient(0);
    !scalar_coeff.is_zero()
}

/// Invert a multivector. Tries blade → versor → general paths.
/// # Errors
/// Returns `InverseError` if the Mv is zero, null, or not invertible.
pub fn inverse(m: &Mv, sig: &Signature) -> Result<Mv, InverseError> {
    if m.is_zero() {
        return Err(InverseError::ZeroMv);
    }

    // Path 1: Single blade
    if m.len() == 1 {
        return blade_inverse(m, sig);
    }

    // Path 2: Versor (M * rev(M) = scalar)
    if let Some(inv) = try_versor_inverse(m, sig) {
        return Ok(inv);
    }

    // Path 3: Even/odd decomposition
    general_inverse(m, sig)
}

/// Blade inverse: B⁻¹ = B / B².
/// Only valid for single-blade Mvs with nonzero self-product.
fn blade_inverse(m: &Mv, sig: &Signature) -> Result<Mv, InverseError> {
    debug_assert_eq!(m.len(), 1);
    let (mask, coeff) = m.blades().next().unwrap();
    let (_, sign) = crate::algebra::product_new::blade_product(mask, mask, sig);
    if sign == 0 {
        return Err(InverseError::NullElement);
    }
    // B² = sign * coeff² (as scalar)
    let b_squared = coeff.clone() * coeff.clone() * Scalar::from(sign as i64);
    if b_squared.is_zero() {
        return Err(InverseError::NullElement);
    }
    // B⁻¹ = B / B² = (1/B²) * B
    let inv_b_sq = Scalar::from(1i64) / b_squared;
    Ok(m.scale(&inv_b_sq))
}

/// Try versor inverse: rev(M) / scalar_part(M * rev(M)).
/// Returns None if M * rev(M) is not a pure nonzero scalar.
fn try_versor_inverse(m: &Mv, sig: &Signature) -> Option<Mv> {
    let m_rev = ops::reverse(m);
    let product = ops::geometric(m, &m_rev, sig);

    // Check: must be a pure nonzero scalar
    if product.len() != 1 {
        return None;
    }
    let ns = product.coefficient(0);
    if ns.is_zero() {
        return None;
    }

    // inv(M) = rev(M) / norm²
    let inv_ns = Scalar::from(1i64) / ns;
    Some(m_rev.scale(&inv_ns))
}

/// General inverse via Clifford conjugation structure.
///
/// For an Mv `M` in a non-degenerate algebra, we use the identity:
///   M * conjugate(M) = M * grade_involution(reverse(M))
/// If this product is a scalar, then inv(M) = conj(M) / (M * conj(M)).
///
/// If that fails, try the Shirokov method:
///   For M in Cl(p,q), compute the sequence:
///     W₁ = M, N₁ = tr(W₁)
///     Wₖ = M * (Wₖ₋₁ - Nₖ₋₁), Nₖ = tr(Wₖ) / k
///   where tr(X) = scalar part of X.
///   After n steps (n = dimension of algebra), either Nₙ ≠ 0 and
///   inv(M) = (Wₙ₋₁ - Nₙ₋₁) / Nₙ, or M is not invertible.
fn general_inverse(m: &Mv, sig: &Signature) -> Result<Mv, InverseError> {
    // Try conjugate path first
    let conj = ops::conjugate(m);
    let m_conj = ops::geometric(m, &conj, sig);
    if m_conj.len() == 1 {
        let ns = m_conj.coefficient(0);
        if !ns.is_zero() {
            let inv_ns = Scalar::from(1i64) / ns;
            return Ok(conj.scale(&inv_ns));
        }
    }

    // Shirokov's method (Faddeev-LeVerrier for Clifford algebras)
    let n = sig.dimension() as usize;
    let max_iter = n.min(64); // cap iterations for large algebras
    let mut w = m.clone();
    let mut prev_adj = Mv::scalar(Scalar::from(1i64)); // identity

    for k in 1..=max_iter {
        let trace_k = w.coefficient(0); // scalar part
        let n_k = trace_k / Scalar::from(k as i64);

        // adj_k = W_k - N_k * I
        let adj_k = w.clone() - Mv::scalar(n_k.clone());

        if k == max_iter {
            // Final step: N_n should be the "determinant"
            // inv(M) = adj_{n-1} / N_n
            if n_k.is_zero() {
                return Err(InverseError::NotInvertible);
            }
            let inv_det = Scalar::from(1i64) / n_k;
            return Ok(prev_adj.scale(&inv_det));
        }

        prev_adj = adj_k.clone();
        w = ops::geometric(m, &adj_k, sig);
    }

    Err(InverseError::DidNotConverge {
        iterations: max_iter,
    })
}

/// Sandwich with proper normalization: `M * v * inv(M)`.
/// Unlike `sandwich(m, v) = m * v * rev(m)` which assumes unit norm,
/// this divides by `norm²(M)` to handle non-unit versors.
///
/// For unit versors (rotors), this equals `sandwich(m, v)`.
/// For non-unit versors, this is the correct transformation.
/// # Errors
/// Returns `InverseError` if `m` is not invertible.
pub fn sandwich_normalized(m: &Mv, v: &Mv, sig: &Signature) -> Result<Mv, InverseError> {
    let m_inv = inverse(m, sig)?;
    let mv = ops::geometric(m, v, sig);
    Ok(ops::geometric(&mv, &m_inv, sig))
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::algebra::blade_new::BladeMask;
    use crate::scalar::Rat;

    fn sig_vga3() -> Signature {
        Signature::new(0, 0, 3).unwrap()
    }
    fn sig_pga3() -> Signature {
        Signature::new(0, 1, 3).unwrap()
    }
    fn sig_cga2() -> Signature {
        Signature::new(1, 0, 3).unwrap()
    }
    fn e(k: u8) -> Mv {
        Mv::generator(k)
    }

    // ─── Blade inverse ───

    #[test]
    fn blade_inverse_vector() {
        let sig = sig_vga3();
        let v = e(0); // e₀, e₀² = 1
        let inv = inverse(&v, &sig).unwrap();
        let product = ops::geometric(&v, &inv, &sig);
        assert_eq!(product.coefficient(0), Scalar::from(1));
        assert_eq!(product.len(), 1);
    }

    #[test]
    fn blade_inverse_scaled_vector() {
        let sig = sig_vga3();
        let v = Mv::term(0b001, Scalar::from(3)); // 3e₀
        let inv = inverse(&v, &sig).unwrap();
        let product = ops::geometric(&v, &inv, &sig);
        assert_eq!(product.coefficient(0), Scalar::from(1));
    }

    #[test]
    fn blade_inverse_bivector() {
        let sig = sig_vga3();
        let bv = Mv::term(0b011, Scalar::from(1)); // e₀₁
        let inv = inverse(&bv, &sig).unwrap();
        let product = ops::geometric(&bv, &inv, &sig);
        assert_eq!(product.coefficient(0), Scalar::from(1));
        assert_eq!(product.len(), 1);
    }

    #[test]
    fn blade_inverse_negative_square() {
        // In Cl(1,0,0): e₀² = -1, so inv(e₀) = -e₀
        let sig = Signature::new(1, 0, 0).unwrap();
        let v = e(0);
        let inv = inverse(&v, &sig).unwrap();
        assert_eq!(inv.coefficient(0b1), Scalar::from(-1));
    }

    #[test]
    fn blade_inverse_degenerate_fails() {
        // In PGA: e₀² = 0
        let sig = sig_pga3();
        let v = e(0); // degenerate generator
        let result = inverse(&v, &sig);
        assert!(result.is_err());
    }

    // ─── Versor inverse ───

    #[test]
    fn versor_inverse_rotor_vga3() {
        // Rotor R = e₀e₁ (unit: R*rev(R) = e₀₁ * e₁₀ = e₀₁*(-e₀₁) = -(-1) = 1)
        let sig = sig_vga3();
        let r = ops::geometric(&e(0), &e(1), &sig); // e₀₁
        assert!(is_versor(&r, &sig));
        let inv = inverse(&r, &sig).unwrap();
        let product = ops::geometric(&r, &inv, &sig);
        assert_eq!(product.coefficient(0), Scalar::from(1));
        assert_eq!(product.len(), 1);
    }

    #[test]
    fn versor_check_vector() {
        let sig = sig_vga3();
        assert!(is_versor(&e(0), &sig));
    }

    #[test]
    fn versor_check_zero() {
        let sig = sig_vga3();
        assert!(!is_versor(&Mv::new(), &sig));
    }

    #[test]
    fn versor_check_null_vector() {
        // In CGA2: null vector e₀ + e₁ has norm² = 0
        let sig = sig_cga2();
        let null = Mv::from_rat_terms(&[(0b01, Rat::from(1)), (0b10, Rat::from(1))]);
        assert!(!is_versor(&null, &sig)); // null vectors are not versors
    }

    // ─── Inverse of multivectors ───

    #[test]
    fn inverse_scalar() {
        let sig = sig_vga3();
        let s = Mv::scalar(Scalar::from(5));
        let inv = inverse(&s, &sig).unwrap();
        let product = ops::geometric(&s, &inv, &sig);
        assert_eq!(product.coefficient(0), Scalar::from(1));
    }

    #[test]
    fn inverse_sum_of_vectors() {
        // v = e₀ + e₁ in VGA3: v² = 1+1 = 2, so inv(v) = v/2
        let sig = sig_vga3();
        let v = Mv::from_rat_terms(&[(0b001, Rat::from(1)), (0b010, Rat::from(1))]);
        let inv = inverse(&v, &sig).unwrap();
        let product = ops::geometric(&v, &inv, &sig);
        assert_eq!(product.coefficient(0), Scalar::from(1));
        assert_eq!(product.len(), 1);
    }

    #[test]
    fn inverse_zero_fails() {
        let sig = sig_vga3();
        assert!(matches!(
            inverse(&Mv::new(), &sig),
            Err(InverseError::ZeroMv)
        ));
    }

    #[test]
    fn inverse_null_vector_fails() {
        // CGA null vector: norm² = 0
        let sig = sig_cga2();
        let null = Mv::from_rat_terms(&[(0b01, Rat::from(1)), (0b10, Rat::from(1))]);
        assert!(inverse(&null, &sig).is_err());
    }

    // ─── Sandwich normalized ───

    #[test]
    fn sandwich_normalized_unit_versor() {
        // For a unit vector, sandwich_normalized should equal sandwich
        let sig = sig_vga3();
        let r = e(0);
        let v = e(1);
        let s1 = ops::sandwich(&r, &v, &sig);
        let s2 = sandwich_normalized(&r, &v, &sig).unwrap();
        assert_eq!(s1, s2);
    }

    #[test]
    fn sandwich_normalized_scaled_versor() {
        // 2e₀ is a versor with norm² = 4
        // sandwich(2e₀, e₁) = 2e₀ * e₁ * 2e₀ = 4 * (e₀*e₁*e₀) = -4e₁
        // sandwich_normalized(2e₀, e₁) = 2e₀ * e₁ * inv(2e₀) = -e₁
        let sig = sig_vga3();
        let r = Mv::term(0b001, Scalar::from(2)); // 2e₀
        let v = e(1);
        let result = sandwich_normalized(&r, &v, &sig).unwrap();
        assert_eq!(result.coefficient(0b010), Scalar::from(-1));
        assert_eq!(result.len(), 1);
    }

    // ─── Roundtrip stress ───

    #[test]
    fn inverse_roundtrip_all_single_blades_vga3() {
        let sig = sig_vga3();
        // Every single blade in VGA3 should be invertible
        for k in 0..sig.n() {
            let blade = e(k);
            let inv = inverse(&blade, &sig).unwrap();
            let product = ops::geometric(&blade, &inv, &sig);
            assert_eq!(
                product.coefficient(0),
                Scalar::from(1),
                "roundtrip failed for generator {}",
                k
            );
            assert_eq!(product.len(), 1);
        }
    }

    #[test]
    fn inverse_roundtrip_bivectors_vga3() {
        let sig = sig_vga3();
        let blades: Vec<BladeMask> = vec![0b011, 0b101, 0b110];
        for mask in blades {
            let bv = Mv::term(mask, Scalar::from(1));
            let inv = inverse(&bv, &sig).unwrap();
            let product = ops::geometric(&bv, &inv, &sig);
            assert_eq!(
                product.coefficient(0),
                Scalar::from(1),
                "roundtrip failed for blade {:#b}",
                mask
            );
        }
    }
}