clifford 0.3.0

//! Basis blades: the fundamental building blocks of geometric algebra.
//!
//! # What is a Blade?
//!
//! A **blade** is a fundamental geometric object that represents an oriented
//! subspace of a given dimension. Think of blades as the "atoms" from which
//! all geometric algebra elements are built.
//!
//! ## Grade 0: Scalars
//!
//! The simplest blade is the **scalar** (grade 0). Scalars are just numbers
//! with no geometric direction—they represent pure magnitude.
//!
//! ## Grade 1: Vectors
//!
//! **Vectors** (grade 1) are blades that represent directed line segments.
//! In 3D, we have three basis vectors: `e₁`, `e₂`, `e₃`, pointing along
//! the x, y, and z axes respectively. Any vector can be written as a
//! linear combination: `v = a·e₁ + b·e₂ + c·e₃`.
//!
//! ## Grade 2: Bivectors
//!
//! **Bivectors** (grade 2) represent oriented plane segments. They are formed
//! by the **wedge product** (∧) of two vectors:
//!
//! ```text
//! e₁ ∧ e₂ = e₁₂  (the xy-plane)
//! e₂ ∧ e₃ = e₂₃  (the yz-plane)
//! e₁ ∧ e₃ = e₁₃  (the xz-plane)
//! ```
//!
//! The bivector `e₁₂` represents a unit area in the xy-plane with a specific
//! orientation (counterclockwise when viewed from +z). Bivectors encode
//! rotations: a rotation in the xy-plane is generated by the bivector `e₁₂`.
//!
//! ## Grade 3: Trivectors (and higher)
//!
//! **Trivectors** represent oriented volumes. In 3D, there's only one basis
//! trivector: `e₁₂₃ = e₁ ∧ e₂ ∧ e₃`, called the **pseudoscalar**. It
//! represents a unit cube with a specific handedness.
//!
//! In higher dimensions, we get quadvectors (grade 4) and so on.
//!
//! # The Geometric Product
//!
//! Blades multiply using the **geometric product**, which combines:
//!
//! 1. **Wedge product** (∧): Creates higher-grade blades when vectors are
//!    independent. `e₁ ∧ e₂ = e₁₂`
//!
//! 2. **Contraction**: When a vector appears twice, it "contracts" using the
//!    metric. In Euclidean space: `e₁ · e₁ = e₁² = 1`
//!
//! The geometric product of two vectors `a` and `b` is:
//! ```text
//! ab = a·b + a∧b  (dot product + wedge product)
//! ```
//!
//! ## Why Vectors Anticommute
//!
//! For distinct basis vectors, swapping order gives a minus sign:
//! ```text
//! e₁e₂ = e₁₂
//! e₂e₁ = -e₁₂  (one swap → one minus sign)
//! ```
//!
//! This is because the geometric product is associative but not commutative.
//! The rule is: each swap of adjacent basis vectors introduces a factor of -1.
//!
//! # Representation in This Library
//!
//! We represent a basis blade by its **index**, a bitmask where bit `i` is set
//! if basis vector `eᵢ` is present:
//!
//! | Index | Binary | Blade | Grade | Geometric Meaning |
//! |-------|--------|-------|-------|-------------------|
//! | 0 | `000` | `1` | 0 | Scalar (no direction) |
//! | 1 | `001` | `e₁` | 1 | Vector along x |
//! | 2 | `010` | `e₂` | 1 | Vector along y |
//! | 3 | `011` | `e₁₂` | 2 | xy-plane bivector |
//! | 4 | `100` | `e₃` | 1 | Vector along z |
//! | 5 | `101` | `e₁₃` | 2 | xz-plane bivector |
//! | 6 | `110` | `e₂₃` | 2 | yz-plane bivector |
//! | 7 | `111` | `e₁₂₃` | 3 | Pseudoscalar (volume) |
//!
//! The grade equals the number of set bits (population count).

use core::fmt;

/// A basis blade in a Clifford algebra.
///
/// A blade represents an oriented subspace: scalars (grade 0), vectors (grade 1),
/// bivectors (grade 2), trivectors (grade 3), and so on. This type represents
/// a single basis blade, identified by which basis vectors it contains.
///
/// # Representation
///
/// Internally, a blade is stored as a bitmask index. Bit `i` is set if and only
/// if basis vector `eᵢ` is part of this blade. For example:
///
/// - `Blade::scalar()` has index `0b000` = 0
/// - `Blade::basis_vector(0)` (e₁) has index `0b001` = 1
/// - The bivector e₁₂ has index `0b011` = 3
///
/// # Example
///
/// ```
/// use clifford::basis::Blade;
///
/// // Create some basis blades
/// let scalar = Blade::scalar();
/// let e1 = Blade::basis_vector(0);
/// let e2 = Blade::basis_vector(1);
/// let e12 = Blade::from_index(0b011);
///
/// // Check their grades
/// assert_eq!(scalar.grade(), 0);
/// assert_eq!(e1.grade(), 1);
/// assert_eq!(e12.grade(), 2);
///
/// // Multiply blades (using Euclidean metric where all eᵢ² = +1)
/// let euclidean = |_: usize| 1i8;
///
/// // e₁ * e₂ = e₁₂
/// let (sign, result) = e1.product(&e2, euclidean);
/// assert_eq!(sign, 1);
/// assert_eq!(result, e12);
///
/// // e₂ * e₁ = -e₁₂ (anticommutes!)
/// let (sign, result) = e2.product(&e1, euclidean);
/// assert_eq!(sign, -1);
/// assert_eq!(result, e12);
/// ```
#[derive(Clone, Copy, PartialEq, Eq, Hash)]
pub struct Blade {
    /// Bitmask index: bit i is set iff basis vector eᵢ is present
    index: usize,
}

impl Blade {
    /// Creates a blade from its bitmask index.
    ///
    /// The index encodes which basis vectors are present: bit `i` set means
    /// basis vector `eᵢ` is part of this blade.
    ///
    /// # Example
    ///
    /// ```
    /// use clifford::basis::Blade;
    ///
    /// let e12 = Blade::from_index(0b011); // e₁ and e₂ → e₁₂
    /// assert_eq!(e12.grade(), 2);
    /// ```
    #[inline]
    pub const fn from_index(index: usize) -> Self {
        Self { index }
    }

    /// Returns the scalar blade (grade 0, the multiplicative identity).
    ///
    /// The scalar blade represents pure magnitude with no direction.
    /// It is the identity element for the geometric product: `1 * A = A * 1 = A`.
    ///
    /// # Example
    ///
    /// ```
    /// use clifford::basis::Blade;
    ///
    /// let one = Blade::scalar();
    /// assert_eq!(one.grade(), 0);
    /// assert_eq!(one.index(), 0);
    /// ```
    #[inline]
    pub const fn scalar() -> Self {
        Self { index: 0 }
    }

    /// Creates a basis vector blade (grade 1).
    ///
    /// Basis vector `eᵢ` is created by `basis_vector(i)`. These are the
    /// fundamental directional elements from which all other blades are built.
    ///
    /// # Arguments
    ///
    /// * `i` - The index of the basis vector (0-indexed)
    ///
    /// # Example
    ///
    /// ```
    /// use clifford::basis::Blade;
    ///
    /// let e1 = Blade::basis_vector(0); // e₁
    /// let e2 = Blade::basis_vector(1); // e₂
    /// let e3 = Blade::basis_vector(2); // e₃
    ///
    /// assert_eq!(e1.grade(), 1);
    /// assert_eq!(e1.index(), 0b001);
    /// assert_eq!(e2.index(), 0b010);
    /// assert_eq!(e3.index(), 0b100);
    /// ```
    #[inline]
    pub const fn basis_vector(i: usize) -> Self {
        Self { index: 1 << i }
    }

    /// Returns the bitmask index of this blade.
    ///
    /// Bit `i` is set if basis vector `eᵢ` is part of this blade.
    #[inline]
    pub const fn index(&self) -> usize {
        self.index
    }

    /// Returns the grade (dimension of the subspace) of this blade.
    ///
    /// - Grade 0: Scalar
    /// - Grade 1: Vector (line)
    /// - Grade 2: Bivector (plane)
    /// - Grade 3: Trivector (volume)
    /// - Grade k: k-vector (k-dimensional subspace)
    ///
    /// The grade equals the number of basis vectors in the blade, which
    /// is the population count of the index.
    ///
    /// # Example
    ///
    /// ```
    /// use clifford::basis::Blade;
    ///
    /// assert_eq!(Blade::scalar().grade(), 0);
    /// assert_eq!(Blade::basis_vector(0).grade(), 1);
    /// assert_eq!(Blade::from_index(0b111).grade(), 3); // e₁₂₃
    /// ```
    #[inline]
    pub const fn grade(&self) -> usize {
        self.index.count_ones() as usize
    }

    /// Checks if this blade contains basis vector `eᵢ`.
    ///
    /// # Example
    ///
    /// ```
    /// use clifford::basis::Blade;
    ///
    /// let e12 = Blade::from_index(0b011); // e₁₂
    /// assert!(e12.contains_vector(0));    // contains e₁
    /// assert!(e12.contains_vector(1));    // contains e₂
    /// assert!(!e12.contains_vector(2));   // does not contain e₃
    /// ```
    #[inline]
    pub const fn contains_vector(&self, i: usize) -> bool {
        (self.index & (1 << i)) != 0
    }

    /// Computes the geometric product of two basis blades.
    ///
    /// Returns `(sign, result)` where:
    /// - `sign` is -1, 0, or +1
    /// - `result` is the resulting blade
    ///
    /// A sign of 0 indicates the product is zero (from a null basis vector).
    ///
    /// # The Geometric Product
    ///
    /// When multiplying `e_A * e_B`:
    /// 1. Vectors in both A and B "contract" using the metric
    /// 2. Remaining vectors combine via the wedge product
    /// 3. Sign changes occur from reordering vectors to canonical form
    ///
    /// # Arguments
    ///
    /// * `other` - The blade to multiply with
    /// * `metric` - Function returning eᵢ² for each basis vector i
    ///
    /// # Example
    ///
    /// ```
    /// use clifford::basis::Blade;
    ///
    /// let euclidean = |_: usize| 1i8; // All eᵢ² = +1
    ///
    /// let e1 = Blade::basis_vector(0);
    /// let e2 = Blade::basis_vector(1);
    ///
    /// // e₁ * e₁ = 1 (vector squared = scalar)
    /// let (sign, result) = e1.product(&e1, euclidean);
    /// assert_eq!(sign, 1);
    /// assert_eq!(result, Blade::scalar());
    ///
    /// // e₁ * e₂ = e₁₂ (independent vectors → bivector)
    /// let (sign, result) = e1.product(&e2, euclidean);
    /// assert_eq!(sign, 1);
    /// assert_eq!(result.grade(), 2);
    /// ```
    #[inline]
    pub fn product<F>(&self, other: &Self, metric: F) -> (i8, Self)
    where
        F: Fn(usize) -> i8,
    {
        let result_index = self.index ^ other.index;
        let sign = compute_sign(self.index, other.index, metric);
        (
            sign,
            Self {
                index: result_index,
            },
        )
    }

    /// Checks if this blade anticommutes with another.
    ///
    /// Two blades anticommute if `A * B = -B * A`. This happens when
    /// the product of their grades is odd.
    ///
    /// # Example
    ///
    /// ```
    /// use clifford::basis::Blade;
    ///
    /// let e1 = Blade::basis_vector(0);
    /// let e2 = Blade::basis_vector(1);
    ///
    /// // Two vectors anticommute: e₁e₂ = -e₂e₁
    /// assert!(e1.anticommutes_with(&e2));
    ///
    /// // Scalar commutes with everything
    /// assert!(!Blade::scalar().anticommutes_with(&e1));
    /// ```
    #[inline]
    pub const fn anticommutes_with(&self, other: &Self) -> bool {
        let grade_product = self.grade() * other.grade();
        grade_product % 2 == 1
    }
}

impl fmt::Debug for Blade {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        if self.index == 0 {
            write!(f, "Blade(1)")
        } else {
            write!(f, "Blade(e")?;
            for i in 0..usize::BITS as usize {
                if self.contains_vector(i) {
                    write!(f, "{}", i + 1)?; // Use 1-indexed for display
                }
            }
            write!(f, ")")
        }
    }
}

impl fmt::Display for Blade {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        if self.index == 0 {
            write!(f, "1")
        } else {
            write!(f, "e")?;
            for i in 0..usize::BITS as usize {
                if self.contains_vector(i) {
                    // Use subscript digits for nice display
                    let subscripts = ['₁', '₂', '₃', '₄', '₅', '₆', '₇', '₈', '₉'];
                    if i < subscripts.len() {
                        write!(f, "{}", subscripts[i])?;
                    } else {
                        write!(f, "_{}", i + 1)?;
                    }
                }
            }
            Ok(())
        }
    }
}

// ============================================================================
// Internal sign computation
// ============================================================================

/// Computes the sign when multiplying basis blades a and b.
#[inline]
fn compute_sign<F>(a: usize, b: usize, metric: F) -> i8
where
    F: Fn(usize) -> i8,
{
    let common = a & b;

    // Check for null vectors (zero metric → zero product)
    let mut temp_common = common;
    while temp_common != 0 {
        let i = temp_common.trailing_zeros() as usize;
        if metric(i) == 0 {
            return 0;
        }
        temp_common &= temp_common - 1;
    }

    let swaps = count_swaps(a, b);
    let metric_sign = compute_metric_sign(common, &metric);

    if swaps.is_multiple_of(2) {
        metric_sign
    } else {
        -metric_sign
    }
}

/// Counts swaps needed to reorder the concatenated basis vectors.
#[inline]
fn count_swaps(a: usize, b: usize) -> usize {
    let mut swaps = 0;
    let mut b_remaining = b;

    while b_remaining != 0 {
        let j = b_remaining.trailing_zeros() as usize;
        let mask_above_j = !((1usize << (j + 1)) - 1);
        swaps += (a & mask_above_j).count_ones() as usize;
        b_remaining &= b_remaining - 1;
    }

    swaps
}

/// Computes the sign from metric contractions.
#[inline]
fn compute_metric_sign<F>(common: usize, metric: F) -> i8
where
    F: Fn(usize) -> i8,
{
    let mut sign: i8 = 1;
    let mut remaining = common;

    while remaining != 0 {
        let i = remaining.trailing_zeros() as usize;
        sign *= metric(i);
        remaining &= remaining - 1;
    }

    sign
}

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

    fn euclidean_metric(_i: usize) -> i8 {
        1
    }

    // ========================================================================
    // Blade struct tests
    // ========================================================================

    #[test]
    fn blade_constructors() {
        assert_eq!(Blade::scalar().index(), 0);
        assert_eq!(Blade::basis_vector(0).index(), 1);
        assert_eq!(Blade::basis_vector(1).index(), 2);
        assert_eq!(Blade::basis_vector(2).index(), 4);
    }

    #[test]
    fn blade_grade() {
        assert_eq!(Blade::scalar().grade(), 0);
        assert_eq!(Blade::basis_vector(0).grade(), 1);
        assert_eq!(Blade::from_index(0b011).grade(), 2);
        assert_eq!(Blade::from_index(0b111).grade(), 3);
    }

    #[test]
    fn blade_contains_vector() {
        let e12 = Blade::from_index(0b011);
        assert!(e12.contains_vector(0));
        assert!(e12.contains_vector(1));
        assert!(!e12.contains_vector(2));
    }

    #[test]
    fn blade_display() {
        assert_eq!(format!("{}", Blade::scalar()), "1");
        assert_eq!(format!("{}", Blade::basis_vector(0)), "e₁");
        assert_eq!(format!("{}", Blade::from_index(0b011)), "e₁₂");
        assert_eq!(format!("{}", Blade::from_index(0b111)), "e₁₂₃");
    }

    #[test]
    fn blade_debug() {
        assert_eq!(format!("{:?}", Blade::scalar()), "Blade(1)");
        assert_eq!(format!("{:?}", Blade::basis_vector(0)), "Blade(e1)");
    }

    // ========================================================================
    // Product tests using Blade struct
    // ========================================================================

    #[test]
    fn vector_products_blade() {
        let e1 = Blade::basis_vector(0);
        let e2 = Blade::basis_vector(1);

        // e₁ * e₁ = 1
        let (sign, result) = e1.product(&e1, euclidean_metric);
        assert_eq!(sign, 1);
        assert_eq!(result, Blade::scalar());

        // e₁ * e₂ = e₁₂
        let (sign, result) = e1.product(&e2, euclidean_metric);
        assert_eq!(sign, 1);
        assert_eq!(result, Blade::from_index(0b011));

        // e₂ * e₁ = -e₁₂
        let (sign, result) = e2.product(&e1, euclidean_metric);
        assert_eq!(sign, -1);
        assert_eq!(result, Blade::from_index(0b011));
    }

    #[test]
    fn bivector_products_blade() {
        let e2 = Blade::basis_vector(1);
        let e12 = Blade::from_index(0b011);

        // e₁₂ * e₁₂ = -1
        let (sign, result) = e12.product(&e12, euclidean_metric);
        assert_eq!(sign, -1);
        assert_eq!(result, Blade::scalar());

        // e₁₂ * e₂ = e₁
        let (sign, result) = e12.product(&e2, euclidean_metric);
        assert_eq!(sign, 1);
        assert_eq!(result, Blade::basis_vector(0));
    }

    // ========================================================================
    // Property-based tests
    // ========================================================================

    proptest! {
        #[test]
        fn result_is_xor(a in 0usize..64, b in 0usize..64) {
            let blade_a = Blade::from_index(a);
            let blade_b = Blade::from_index(b);
            let (_, result) = blade_a.product(&blade_b, euclidean_metric);
            prop_assert_eq!(result.index(), a ^ b);
        }

        #[test]
        fn scalar_identity(a in 0usize..64) {
            let one = Blade::scalar();
            let blade = Blade::from_index(a);

            let (sign, result) = one.product(&blade, euclidean_metric);
            prop_assert_eq!(sign, 1);
            prop_assert_eq!(result, blade);

            let (sign, result) = blade.product(&one, euclidean_metric);
            prop_assert_eq!(sign, 1);
            prop_assert_eq!(result, blade);
        }

        #[test]
        fn sign_associativity(a in 0usize..16, b in 0usize..16, c in 0usize..16) {
            let blade_a = Blade::from_index(a);
            let blade_b = Blade::from_index(b);
            let blade_c = Blade::from_index(c);

            let (s_bc, bc) = blade_b.product(&blade_c, euclidean_metric);
            let (s_a_bc, _) = blade_a.product(&bc, euclidean_metric);

            let (s_ab, ab) = blade_a.product(&blade_b, euclidean_metric);
            let (s_ab_c, _) = ab.product(&blade_c, euclidean_metric);

            prop_assert_eq!(
                (s_a_bc as i32) * (s_bc as i32),
                (s_ab_c as i32) * (s_ab as i32)
            );
        }

        #[test]
        fn grade_is_popcount(index in 0usize..256) {
            let blade = Blade::from_index(index);
            prop_assert_eq!(blade.grade(), index.count_ones() as usize);
        }
    }

    // ========================================================================
    // Metric variation tests
    // ========================================================================

    #[test]
    fn minkowski_metric() {
        let minkowski = |i: usize| if i == 0 { 1 } else { -1 };
        let e0 = Blade::basis_vector(0);
        let e1 = Blade::basis_vector(1);

        // e₀² = +1 (time-like)
        let (sign, result) = e0.product(&e0, minkowski);
        assert_eq!(sign, 1);
        assert_eq!(result, Blade::scalar());

        // e₁² = -1 (space-like)
        let (sign, result) = e1.product(&e1, minkowski);
        assert_eq!(sign, -1);
        assert_eq!(result, Blade::scalar());
    }

    #[test]
    fn null_metric() {
        let pga_metric = |i: usize| if i == 0 { 0 } else { 1 };
        let e0 = Blade::basis_vector(0);

        // e₀² = 0 (null/degenerate)
        let (sign, _) = e0.product(&e0, pga_metric);
        assert_eq!(sign, 0);
    }
}