cova-algebra 0.2.2

//! Implementation of Clifford algebras, also known as geometric algebras.
//!
//! A Clifford algebra is a unital associative algebra that generalizes the real numbers, complex
//! numbers, quaternions, and several other hypercomplex number systems. The key feature of a
//! Clifford algebra is that it is generated by a vector space equipped with a quadratic form.
//!
//! # Key Concepts
//!
//! - **Quadratic Form**: A function Q: V → F that satisfies Q(av) = a²Q(v) for all a ∈ F and v ∈ V.
//!   In this implementation, we assume the quadratic form is diagonal, represented by a vector of
//!   coefficients. The quadratic form is encoded at the type level using the
//!   [`QuadraticFormMarker`] trait, ensuring that only elements with the same quadratic form can be
//!   combined.
//!
//! - **Basis Blades**: The fundamental building blocks of a Clifford algebra. In an n-dimensional
//!   space, there are 2ⁿ basis blades, formed by taking the exterior product of basis vectors.
//!
//! - **Grade**: The number of basis vectors in a blade. A scalar has grade 0, a vector has grade 1,
//!   a bivector has grade 2, and so on.
//!
//! # Implementation Details
//!
//! The implementation uses a bit-based representation for basis blades:
//!
//! - Each basis blade is represented by a bit position in a 2ⁿ-dimensional vector
//! - The grade of a blade is determined by the number of set bits
//! - Multiplication is implemented using the geometric product, which combines the exterior product
//!   and the inner product defined by the quadratic form
//! - **Compile-time safety**: The quadratic form is encoded as a type parameter, so operations
//!   between elements with different quadratic forms result in compile-time errors
//!
//! # Examples
//!
//! ```
//! #![feature(generic_const_exprs)]
//! use cova_algebra::algebras::clifford::{CliffordAlgebra, Signature};
//!
//! // Create a 3D Clifford algebra with signature (+, +, -)
//! type MySignature = Signature<3, 1, 1, -1>;
//! let algebra = CliffordAlgebra::<f64, 3, MySignature>::new();
//!
//! // Create basis vectors
//! let e0 = algebra.blade([0]);
//! let e1 = algebra.blade([1]);
//! let e2 = algebra.blade([2]);
//!
//! // Example: e1 * e01 = -e0 (using the quadratic form)
//! let e01 = algebra.blade([0, 1]);
//! assert_eq!(format!("{}", e1 * e01), "-1e₀");
//!
//! // Example: e2 * e2 = -1 (using the quadratic form)
//! assert_eq!(format!("{}", e2.clone() * e2), "-1");
//! ```
//!
//! # Predefined Signatures
//!
//! Common quadratic forms are provided as type aliases:
//! - [`Euclidean2D`], [`Euclidean3D`], [`Euclidean4D`]: Standard Euclidean spaces
//! - [`Minkowski`]: Spacetime with signature (+, -, -, -)
//! - [`AntiMinkowski`]: Spacetime with signature (-, +, +, +)
//!
//! # Geometric Interpretation
//!
//! The geometric product of two vectors a and b can be decomposed as:
//!
//! a * b = a·b + a∧b
//!
//! where:
//! - a·b is the inner product (symmetric part)
//! - a∧b is the exterior product (antisymmetric part)
//!
//! This decomposition is fundamental to geometric algebra and allows for a unified treatment of
//! various geometric operations.
//!
//! # Common Use Cases
//!
//! - **3D Euclidean Space**: Used for 3D rotations and transformations
//! - **Spacetime Physics**: Used in special relativity with signature (+, -, -, -)
//! - **Computer Graphics**: Used for efficient representation of transformations
//! - **Robotics**: Used for representing and manipulating rigid body motions
//!
//! # References
//!
//! - [Geometric Algebra for Computer Science](https://www.geometricalgebra.net/)
//! - [Clifford Algebra](https://en.wikipedia.org/wiki/Clifford_algebra)
//! - [Geometric Algebra](https://en.wikipedia.org/wiki/Geometric_algebra)

use std::{
  fmt::{Debug, Display, Formatter},
  marker::PhantomData,
};

use num_traits::One;

use super::*;
use crate::{
  algebras::Algebra,
  groups::{AbelianGroup, Group},
  modules::{LeftModule, RightModule, TwoSidedModule},
  rings::Field,
  tensors::SVector,
};

/// A marker trait for quadratic forms on a vector space.
///
/// This trait enables compile-time verification that Clifford algebra elements
/// have compatible quadratic forms. By encoding the quadratic form as a type parameter,
/// operations between elements with different forms result in compile-time errors.
///
/// A quadratic form Q is a function that satisfies Q(av) = a²Q(v) for all scalars a and vectors v.
/// In this implementation, we assume the quadratic form is diagonal, meaning it can be represented
/// by a vector of coefficients where Q(v) = Σᵢ cᵢvᵢ².
///
/// # Implementing Custom Quadratic Forms
///
/// ```
/// use cova_algebra::{algebras::clifford::QuadraticFormMarker, rings::Field, tensors::SVector};
///
/// // Define a custom quadratic form for 2D split-complex numbers
/// #[derive(Copy, Clone, Debug)]
/// pub struct SplitComplex;
///
/// impl<F: Field + From<f64>> QuadraticFormMarker<F, 2> for SplitComplex {
///   fn coefficients() -> SVector<F, 2> { SVector::from_row_slice(&[F::from(1.0), F::from(-1.0)]) }
/// }
/// ```
pub trait QuadraticFormMarker<F: Field, const N: usize>: 'static + Copy + Debug {
  /// Returns the diagonal coefficients of the quadratic form.
  ///
  /// For a quadratic form Q(v) = Σᵢ cᵢvᵢ², this returns the vector [c₀, c₁, ..., cₙ₋₁].
  fn coefficients() -> SVector<F, N>;

  /// Evaluates the quadratic form at a given vector.
  ///
  /// # Arguments
  ///
  /// * `v` - The vector at which to evaluate the quadratic form
  ///
  /// # Returns
  ///
  /// The value of Q(v) = Σᵢ cᵢvᵢ²
  fn evaluate(v: &SVector<F, N>) -> F {
    let coeffs = Self::coefficients();
    let mut result = F::zero();
    for i in 0..N {
      result += coeffs[i] * v[i] * v[i];
    }
    result
  }
}

/// A generic signature type for defining quadratic forms via const generics.
///
/// This type allows defining quadratic forms using const parameters. For simple cases
/// with up to 4 dimensions using ±1 coefficients, use the predefined type aliases.
///
/// # Type Parameters
///
/// * `N` - The dimension of the vector space
/// * `C0`, `C1`, `C2`, `C3` - Coefficients as i64 values (typically 1, -1, or 0)
///
/// # Examples
///
/// ```
/// use cova_algebra::algebras::clifford::Signature;
///
/// // Define a signature for Cl(2,1) - 2D space with one negative direction
/// type MySig = Signature<3, 1, 1, -1>;
/// ```
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct Signature<
  const N: usize,
  const C0: i64 = 1,
  const C1: i64 = 1,
  const C2: i64 = 1,
  const C3: i64 = 1,
>;

impl<const N: usize, const C0: i64, const C1: i64, const C2: i64, const C3: i64>
  QuadraticFormMarker<f64, N> for Signature<N, C0, C1, C2, C3>
{
  fn coefficients() -> SVector<f64, N> {
    let mut coeffs = SVector::<f64, N>::zeros();
    if N > 0 {
      coeffs[0] = C0 as f64;
    }
    if N > 1 {
      coeffs[1] = C1 as f64;
    }
    if N > 2 {
      coeffs[2] = C2 as f64;
    }
    if N > 3 {
      coeffs[3] = C3 as f64;
    }
    // For N > 4, remaining coefficients default to 0 (could extend pattern)
    coeffs
  }
}

impl<const N: usize, const C0: i64, const C1: i64, const C2: i64, const C3: i64>
  QuadraticFormMarker<f32, N> for Signature<N, C0, C1, C2, C3>
{
  fn coefficients() -> SVector<f32, N> {
    let mut coeffs = SVector::<f32, N>::zeros();
    if N > 0 {
      coeffs[0] = C0 as f32;
    }
    if N > 1 {
      coeffs[1] = C1 as f32;
    }
    if N > 2 {
      coeffs[2] = C2 as f32;
    }
    if N > 3 {
      coeffs[3] = C3 as f32;
    }
    coeffs
  }
}

// ============================================================================
// Predefined Signatures for Common Algebras
// ============================================================================

/// Euclidean 2D signature: (+, +)
pub type Euclidean2D = Signature<2, 1, 1>;

/// Euclidean 3D signature: (+, +, +)
pub type Euclidean3D = Signature<3, 1, 1, 1>;

/// Euclidean 4D signature: (+, +, +, +)
pub type Euclidean4D = Signature<4, 1, 1, 1, 1>;

/// Minkowski spacetime signature: (+, -, -, -)
pub type Minkowski = Signature<4, 1, -1, -1, -1>;

/// Anti-Minkowski spacetime signature: (-, +, +, +)
pub type AntiMinkowski = Signature<4, -1, 1, 1, 1>;

/// Split-complex numbers signature: (+, -)
pub type SplitComplex = Signature<2, 1, -1>;

/// Complex numbers as Cl(0,1): (-)
pub type ComplexSig = Signature<1, -1>;

/// Quaternions as Cl(0,2): (-, -)
pub type QuaternionSig = Signature<2, -1, -1>;

/// A Clifford algebra over a vector space with a given quadratic form.
///
/// A Clifford algebra is a unital associative algebra that generalizes the real numbers,
/// complex numbers, and quaternions. It is generated by a vector space equipped with a
/// quadratic form.
///
/// # Type Parameters
///
/// * `F` - The field over which the algebra is defined
/// * `N` - The dimension of the underlying vector space
/// * `Q` - The quadratic form marker type (ensures compile-time compatibility)
///
/// # Examples
///
/// ```
/// #![feature(generic_const_exprs)]
/// use cova_algebra::algebras::clifford::{CliffordAlgebra, Euclidean3D};
///
/// // Create a 3D Euclidean Clifford algebra
/// let algebra = CliffordAlgebra::<f64, 3, Euclidean3D>::new();
/// let e0 = algebra.blade([0]);
/// let e1 = algebra.blade([1]);
/// ```
pub struct CliffordAlgebra<F: Field, const N: usize, Q: QuadraticFormMarker<F, N>> {
  _marker: PhantomData<(F, Q)>,
}

impl<F: Field + Copy, const N: usize, Q: QuadraticFormMarker<F, N>> CliffordAlgebra<F, N, Q>
where [(); 1 << N]:
{
  /// Creates a new Clifford algebra with the quadratic form specified by the type parameter.
  ///
  /// # Examples
  ///
  /// ```
  /// #![feature(generic_const_exprs)]
  /// use cova_algebra::algebras::clifford::{CliffordAlgebra, Euclidean3D, Minkowski};
  ///
  /// // Create a 3D Euclidean Clifford algebra
  /// let euclidean = CliffordAlgebra::<f64, 3, Euclidean3D>::new();
  ///
  /// // Create a Minkowski spacetime algebra
  /// let spacetime = CliffordAlgebra::<f64, 4, Minkowski>::new();
  /// ```
  pub const fn new() -> Self { Self { _marker: PhantomData } }

  /// Creates a new element in the algebra from a vector of coefficients.
  ///
  /// # Arguments
  ///
  /// * `value` - A vector of coefficients for each basis blade
  ///
  /// # Examples
  ///
  /// ```
  /// #![feature(generic_const_exprs)]
  /// use cova_algebra::{
  ///   algebras::clifford::{CliffordAlgebra, Euclidean3D},
  ///   tensors::SVector,
  /// };
  ///
  /// let algebra = CliffordAlgebra::<f64, 3, Euclidean3D>::new();
  /// let element =
  ///   algebra.element(SVector::<f64, 8>::from_row_slice(&[1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]));
  /// ```
  pub fn element(&self, value: SVector<F, { 1 << N }>) -> CliffordAlgebraElement<F, N, Q> {
    CliffordAlgebraElement { value, _marker: PhantomData }
  }

  /// Creates a basis blade with the given indices.
  ///
  /// A basis blade is a product of basis vectors, represented by their indices.
  /// The indices must be sorted and unique.
  ///
  /// # Arguments
  ///
  /// * `indices` - The indices of the basis vectors that form the blade
  ///
  /// # Panics
  ///
  /// Panics if the indices are not sorted or if any index is out of range.
  ///
  /// # Examples
  ///
  /// ```
  /// #![feature(generic_const_exprs)]
  /// use cova_algebra::algebras::clifford::{CliffordAlgebra, Euclidean3D};
  ///
  /// let algebra = CliffordAlgebra::<f64, 3, Euclidean3D>::new();
  /// let e1 = algebra.blade([1]);
  /// let e01 = algebra.blade([0, 1]);
  /// ```
  pub fn blade<const I: usize>(&self, indices: [usize; I]) -> CliffordAlgebraElement<F, N, Q> {
    // Validate indices are in range and sorted
    for i in 1..I {
      assert!(
        indices[i - 1] < indices[i] && indices[i] < N,
        "Indices must be sorted and in
    range"
      );
    }

    // Convert indices to bit position using our helper function
    let bit_position = Self::blade_indices_to_bit(&indices);

    let mut value = SVector::<F, { 1 << N }>::zeros();
    value[bit_position] = <F as One>::one();

    CliffordAlgebraElement { value, _marker: PhantomData }
  }

  /// Maps a set of basis blade indices back to a bit position.
  /// This is the inverse of [`bit_to_blade_indices`].
  fn blade_indices_to_bit(indices: &[usize]) -> usize {
    let grade = indices.len();
    let mut bit_position = 0;

    // Add up all the positions from lower grades
    for g in 0..grade {
      bit_position += binomial(N, g);
    }

    // Now add the position within this grade
    let mut remaining_bits = 0;
    for (i, &idx) in indices.iter().enumerate() {
      // For each index, add the number of combinations that come before it
      for j in if i == 0 { 0 } else { indices[i - 1] + 1 }..idx {
        remaining_bits += binomial(N - j - 1, grade - i - 1);
      }
    }

    bit_position + remaining_bits
  }
}

impl<F: Field + Copy, const N: usize, Q: QuadraticFormMarker<F, N>> Default
  for CliffordAlgebra<F, N, Q>
where [(); 1 << N]:
{
  fn default() -> Self { Self::new() }
}

/// Computes the binomial coefficient C(n, k)
fn binomial(n: usize, k: usize) -> usize {
  if k > n {
    return 0;
  }
  if k == 0 || k == n {
    return 1;
  }
  let k = std::cmp::min(k, n - k);
  let mut result = 1;
  for i in 1..=k {
    result = result * (n - k + i) / i;
  }
  result
}

/// An element of a Clifford algebra.
///
/// An element is represented as a linear combination of basis blades, where each
/// basis blade is a product of basis vectors. The quadratic form is encoded in the
/// type parameter `Q`, ensuring compile-time compatibility between elements.
///
/// # Type Parameters
///
/// * `F` - The field over which the algebra is defined
/// * `N` - The dimension of the underlying vector space
/// * `Q` - The quadratic form marker type (ensures compile-time compatibility)
///
/// # Examples
///
/// ```
/// #![feature(generic_const_exprs)]
/// use cova_algebra::algebras::clifford::{CliffordAlgebra, Euclidean3D};
///
/// let algebra = CliffordAlgebra::<f64, 3, Euclidean3D>::new();
/// let e1 = algebra.blade([1]);
/// let e2 = algebra.blade([2]);
/// let sum = e1 + e2; // Works: same quadratic form
/// ```
///
/// ```compile_fail
/// #![feature(generic_const_exprs)]
/// use cova_algebra::algebras::clifford::{CliffordAlgebra, Euclidean3D, Minkowski};
///
/// let euclidean = CliffordAlgebra::<f64, 3, Euclidean3D>::new();
/// let minkowski = CliffordAlgebra::<f64, 4, Minkowski>::new();
/// let e1 = euclidean.blade([1]);
/// let f1 = minkowski.blade([1]);
/// let sum = e1 + f1; // Compile error: different quadratic forms!
/// ```
#[derive(Clone, Copy, Debug)]
pub struct CliffordAlgebraElement<F: Field, const N: usize, Q: QuadraticFormMarker<F, N>>
where [(); 1 << N]: {
  value:   SVector<F, { 1 << N }>,
  _marker: PhantomData<Q>,
}

impl<F: Field, const N: usize, Q: QuadraticFormMarker<F, N>> PartialEq
  for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
  fn eq(&self, other: &Self) -> bool { self.value == other.value }
}

// Compile-time safety: operations only work between elements with the same Q type parameter

impl<F: Field, const N: usize, Q: QuadraticFormMarker<F, N>> Add for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
  type Output = Self;

  fn add(self, other: Self) -> Self::Output {
    Self { value: self.value + other.value, _marker: PhantomData }
  }
}

impl<F: Field, const N: usize, Q: QuadraticFormMarker<F, N>> AddAssign
  for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
  fn add_assign(&mut self, rhs: Self) { self.value += rhs.value; }
}

impl<F: Field, const N: usize, Q: QuadraticFormMarker<F, N>> Neg for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
  type Output = Self;

  fn neg(self) -> Self::Output { Self { value: -self.value, _marker: PhantomData } }
}

impl<F: Field, const N: usize, Q: QuadraticFormMarker<F, N>> Sub for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
  type Output = Self;

  fn sub(self, other: Self) -> Self::Output { self + -other }
}

impl<F: Field, const N: usize, Q: QuadraticFormMarker<F, N>> SubAssign
  for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
  fn sub_assign(&mut self, rhs: Self) { self.value -= rhs.value; }
}

impl<F: Field, const N: usize, Q: QuadraticFormMarker<F, N>> Mul for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
  type Output = Self;

  fn mul(self, other: Self) -> Self::Output {
    // Get the quadratic form coefficients from the type parameter
    let quadratic_coeffs = Q::coefficients();

    let mut result = SVector::<F, { 1 << N }>::zeros();

    // For each non-zero component in the first element
    for i in 0..(1 << N) {
      if self.value[i].is_zero() {
        continue;
      }

      // For each non-zero component in the second element
      for j in 0..(1 << N) {
        if other.value[j].is_zero() {
          continue;
        }

        // Get the indices for both basis blades
        let left_indices = bit_to_blade_indices::<N>(i);
        let right_indices = bit_to_blade_indices::<N>(j);

        // Calculate the sign and product indices
        let (sign, product_indices) = multiply_blades::<N>(&left_indices, &right_indices);

        // Calculate the coefficient
        let mut coefficient = self.value[i] * other.value[j];

        // Apply sign
        coefficient = match sign {
          Sign::Positive => coefficient,
          Sign::Negative => -coefficient,
        };

        // Apply quadratic form for any repeated indices
        let mut repeated_indices = Vec::new();
        let mut i = 0;
        let mut j = 0;
        while i < left_indices.len() && j < right_indices.len() {
          match left_indices[i].cmp(&right_indices[j]) {
            std::cmp::Ordering::Equal => {
              repeated_indices.push(left_indices[i]);
              i += 1;
              j += 1;
            },
            std::cmp::Ordering::Less => {
              i += 1;
            },
            std::cmp::Ordering::Greater => {
              j += 1;
            },
          }
        }

        // Multiply by quadratic form coefficients for each repeated index
        for &idx in &repeated_indices {
          coefficient *= quadratic_coeffs[idx];
        }

        // Add to the result
        let product_bit = blade_indices_to_bit::<N>(&product_indices);
        result[product_bit] += coefficient;
      }
    }

    Self { value: result, _marker: PhantomData }
  }
}

impl<F: Field, const N: usize, Q: QuadraticFormMarker<F, N>> MulAssign
  for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
  fn mul_assign(&mut self, rhs: Self) {
    let result = *self * rhs;
    *self = result;
  }
}

impl<F: Field, const N: usize, Q: QuadraticFormMarker<F, N>> Mul<F>
  for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
  type Output = Self;

  fn mul(self, rhs: F) -> Self::Output { Self { value: self.value * rhs, _marker: PhantomData } }
}

impl<F: Field, const N: usize, Q: QuadraticFormMarker<F, N>> Multiplicative
  for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
}

impl<F: Field, const N: usize, Q: QuadraticFormMarker<F, N>> Zero
  for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
  fn zero() -> Self { Self { value: SVector::<F, { 1 << N }>::zeros(), _marker: PhantomData } }

  fn is_zero(&self) -> bool { self.value.iter().all(num_traits::Zero::is_zero) }
}

impl<F: Field, const N: usize, Q: QuadraticFormMarker<F, N>> Additive
  for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
}

impl<F: Field, const N: usize, Q: QuadraticFormMarker<F, N>> Group
  for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
  fn identity() -> Self { Self { value: SVector::<F, { 1 << N }>::zeros(), _marker: PhantomData } }

  fn inverse(&self) -> Self { Self { value: -self.value, _marker: PhantomData } }
}

impl<F: Field, const N: usize, Q: QuadraticFormMarker<F, N>> AbelianGroup
  for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
}

impl<F: Field + Mul<Self>, const N: usize, Q: QuadraticFormMarker<F, N>> LeftModule
  for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
  type Ring = F;
}

impl<F: Field + Mul<Self>, const N: usize, Q: QuadraticFormMarker<F, N>> RightModule
  for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
  type Ring = F;
}

impl<F: Field + Mul<Self>, const N: usize, Q: QuadraticFormMarker<F, N>> TwoSidedModule
  for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
  type Ring = F;
}

impl<F: Field + Mul<Self>, const N: usize, Q: QuadraticFormMarker<F, N>> VectorSpace
  for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
}

impl<F: Field + Mul<Self>, const N: usize, Q: QuadraticFormMarker<F, N>> Algebra
  for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
}

impl<F: Field + Display, const N: usize, Q: QuadraticFormMarker<F, N>> Display
  for CliffordAlgebraElement<F, N, Q>
where [(); 1 << N]:
{
  fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
    let mut first = true;

    // Helper function to write basis element
    let write_basis = |f: &mut Formatter<'_>, indices: &[usize]| -> std::fmt::Result {
      if indices.is_empty() {
        return Ok(());
      }

      write!(f, "e")?;
      for (i, &index) in indices.iter().enumerate() {
        // Convert number to Unicode subscript by converting each digit
        let num = index;
        for digit in num.to_string().chars() {
          let subscript = match digit {
            '0' => "₀",
            '1' => "₁",
            '2' => "₂",
            '3' => "₃",
            '4' => "₄",
            '5' => "₅",
            '6' => "₆",
            '7' => "₇",
            '8' => "₈",
            '9' => "₉",
            _ => panic!("Invalid digit"),
          };
          write!(f, "{subscript}")?;
        }

        if i < indices.len() - 1 {
          write!(f, "‚")?;
        }
      }
      Ok(())
    };

    // Print each grade in order
    for grade in 0..=N {
      // Generate all possible combinations of indices for this grade
      let mut indices = Vec::with_capacity(grade);
      let mut combinations = Vec::new();
      generate_combinations(0, N, grade, &mut indices, &mut combinations);

      // Sort combinations by their bit position
      combinations.sort_by_key(|indices| CliffordAlgebra::<F, N, Q>::blade_indices_to_bit(indices));

      // Print each combination in order
      for indices in combinations {
        let bit_position = CliffordAlgebra::<F, N, Q>::blade_indices_to_bit(&indices);
        if !self.value[bit_position].is_zero() {
          if !first {
            write!(f, " + ")?;
          }
          write!(f, "{}", self.value[bit_position])?;
          if grade > 0 {
            write_basis(f, &indices)?;
          }
          first = false;
        }
      }
    }

    if first {
      write!(f, "0")?;
    }
    Ok(())
  }
}

/// Helper function to generate all combinations of indices for a given grade
fn generate_combinations(
  start: usize,
  n: usize,
  k: usize,
  current: &mut Vec<usize>,
  result: &mut Vec<Vec<usize>>,
) {
  if k == 0 {
    result.push(current.clone());
    return;
  }

  for i in start..n {
    current.push(i);
    generate_combinations(i + 1, n, k - 1, current, result);
    current.pop();
  }
}

/// A macro that implements scalar multiplication for Clifford algebra elements.
///
/// This macro generates implementations of the `Mul` trait for scalar types (like `f32` and `f64`)
/// to allow multiplication with `CliffordAlgebraElement`s. It enables syntax like `2.0 * element`
/// in addition to `element * 2.0`.
///
/// # Examples
///
/// ```
/// #![feature(generic_const_exprs)]
/// use cova_algebra::algebras::clifford::{CliffordAlgebra, Euclidean3D};
///
/// let algebra = CliffordAlgebra::<f64, 3, Euclidean3D>::new();
/// let e1 = algebra.blade([1]);
///
/// // Both of these are valid thanks to this macro
/// let doubled1 = 2.0 * e1.clone();
/// let doubled2 = e1 * 2.0;
/// ```
/// This allows scalar multiplication to be commutative, delegating to the existing
/// right-multiplication implementation.
#[macro_export]
macro_rules! impl_mul_scalar_clifford {
  ($($t:ty)*) => ($(
    impl<const N: usize, Q: QuadraticFormMarker<$t, N>> Mul<CliffordAlgebraElement<$t, N, Q>> for $t
    where [(); 1 << N]:
    {
      type Output = CliffordAlgebraElement<$t, N, Q>;

      fn mul(self, rhs: CliffordAlgebraElement<$t, N, Q>) -> Self::Output { rhs * self }
    }
  )*)
}

impl_mul_scalar_clifford!(f32);
impl_mul_scalar_clifford!(f64);

/// Represents the sign of a term in the algebra.
///
/// Used internally to track sign changes during multiplication.
#[derive(Debug)]
pub enum Sign {
  /// Positive sign
  Positive,
  /// Negative sign
  Negative,
}

/// Helper function to multiply two basis blades
fn multiply_blades<const N: usize>(left: &[usize], right: &[usize]) -> (Sign, Vec<usize>) {
  let mut result_indices = Vec::new();
  let mut sign = Sign::Positive;

  // Handle the case where either blade is empty (scalar)
  if left.is_empty() {
    return (Sign::Positive, right.to_vec());
  }
  if right.is_empty() {
    return (Sign::Positive, left.to_vec());
  }

  // Merge the indices while keeping track of sign changes
  let mut i = 0;
  let mut j = 0;

  while i < left.len() && j < right.len() {
    match left[i].cmp(&right[j]) {
      std::cmp::Ordering::Equal => {
        // Same index: apply quadratic form
        i += 1;
        j += 1;
      },
      std::cmp::Ordering::Less => {
        // Left index comes first: no sign change
        result_indices.push(left[i]);
        i += 1;
      },
      std::cmp::Ordering::Greater => {
        // Right index comes first: count swaps for sign
        result_indices.push(right[j]);
        if (left.len() - i) % 2 == 1 {
          sign = match sign {
            Sign::Positive => Sign::Negative,
            Sign::Negative => Sign::Positive,
          };
        }
        j += 1;
      },
    }
  }

  // Add remaining indices
  result_indices.extend_from_slice(&left[i..]);
  result_indices.extend_from_slice(&right[j..]);

  (sign, result_indices)
}

/// Maps a bit position to the corresponding basis blade indices.
fn bit_to_blade_indices<const N: usize>(bits: usize) -> Vec<usize> {
  let mut remaining_bits = bits;
  let mut grade = 0;

  // Find which grade this element belongs to
  while grade <= N {
    let grade_size = binomial(N, grade);
    if remaining_bits < grade_size {
      break;
    }
    remaining_bits -= grade_size;
    grade += 1;
  }

  // Now we know the grade, we need to find which combination of indices
  // corresponds to the remaining_bits position within that grade
  let mut indices = Vec::with_capacity(grade);
  let mut current = 0;

  for _ in 0..grade {
    // Find the next index to include
    while current < N {
      let remaining_combinations = binomial(N - current - 1, grade - indices.len() - 1);
      if remaining_bits < remaining_combinations {
        indices.push(current);
        current += 1;
        break;
      }
      remaining_bits -= remaining_combinations;
      current += 1;
    }
  }

  indices
}

/// Maps a set of basis blade indices back to a bit position.
fn blade_indices_to_bit<const N: usize>(indices: &[usize]) -> usize {
  let grade = indices.len();
  let mut bit_position = 0;

  // Add up all the positions from lower grades
  for g in 0..grade {
    bit_position += binomial(N, g);
  }

  // Now add the position within this grade
  let mut remaining_bits = 0;
  for (i, &idx) in indices.iter().enumerate() {
    for j in if i == 0 { 0 } else { indices[i - 1] + 1 }..idx {
      remaining_bits += binomial(N - j - 1, grade - i - 1);
    }
  }

  bit_position + remaining_bits
}

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

  /// Non-Euclidean signature: (+, +, -)
  type NonEuclidean3D = Signature<3, 1, 1, -1>;

  fn clifford_algebra_non_euclidean() -> CliffordAlgebra<f64, 3, NonEuclidean3D> {
    CliffordAlgebra::new()
  }

  fn clifford_algebra_euclidean() -> CliffordAlgebra<f64, 3, Euclidean3D> { CliffordAlgebra::new() }

  #[test]
  fn test_display_order() {
    let algebra = clifford_algebra_non_euclidean();
    let one =
      algebra.element(SVector::<f64, 8>::from_row_slice(&[1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]));
    let e0 =
      algebra.element(SVector::<f64, 8>::from_row_slice(&[0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]));
    let e1 =
      algebra.element(SVector::<f64, 8>::from_row_slice(&[0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0]));
    let e2 =
      algebra.element(SVector::<f64, 8>::from_row_slice(&[0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0]));
    let e01 =
      algebra.element(SVector::<f64, 8>::from_row_slice(&[0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0]));
    let e02 =
      algebra.element(SVector::<f64, 8>::from_row_slice(&[0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0]));
    let e12 =
      algebra.element(SVector::<f64, 8>::from_row_slice(&[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0]));
    let e012 =
      algebra.element(SVector::<f64, 8>::from_row_slice(&[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]));

    let sum = one + 2.0 * e0 + 3.0 * e1 + 4.0 * e2 + 5.0 * e01 + 6.0 * e02 + 7.0 * e12 + 8.0 * e012;
    assert_eq!(format!("{sum}"), "1 + 2e₀ + 3e₁ + 4e₂ + 5e₀‚₁ + 6e₀‚₂ + 7e₁‚₂ + 8e₀‚₁‚₂");
  }

  #[test]
  fn test_blade_indices_to_bit() {
    // Test scalar (empty indices)
    assert_eq!(CliffordAlgebra::<f64, 3, Euclidean3D>::blade_indices_to_bit(&[]), 0);

    // Test vectors
    assert_eq!(CliffordAlgebra::<f64, 3, Euclidean3D>::blade_indices_to_bit(&[0]), 1);
    assert_eq!(CliffordAlgebra::<f64, 3, Euclidean3D>::blade_indices_to_bit(&[1]), 2);
    assert_eq!(CliffordAlgebra::<f64, 3, Euclidean3D>::blade_indices_to_bit(&[2]), 3);

    // Test bivectors
    assert_eq!(CliffordAlgebra::<f64, 3, Euclidean3D>::blade_indices_to_bit(&[0, 1]), 4);
    assert_eq!(CliffordAlgebra::<f64, 3, Euclidean3D>::blade_indices_to_bit(&[0, 2]), 5);
    assert_eq!(CliffordAlgebra::<f64, 3, Euclidean3D>::blade_indices_to_bit(&[1, 2]), 6);

    // Test trivector
    assert_eq!(CliffordAlgebra::<f64, 3, Euclidean3D>::blade_indices_to_bit(&[0, 1, 2]), 7);
  }

  #[test]
  fn test_blade() {
    let algebra = clifford_algebra_non_euclidean();
    let e1 = algebra.blade([1]);
    assert_eq!(format!("{e1}"), "1e₁");
  }

  #[test]
  fn test_add() {
    let algebra = clifford_algebra_non_euclidean();
    let e1 = algebra.blade([1]);
    let e2 = algebra.blade([2]);
    let sum = e1 + e2;
    assert_eq!(format!("{sum}"), "1e₁ + 1e₂");
  }

  #[test]
  fn test_mul_basic() {
    let algebra = clifford_algebra_non_euclidean();
    let e1 = algebra.blade([1]);
    let e2 = algebra.blade([2]);
    let product = e1 * e2;
    assert_eq!(format!("{product}"), "1e₁‚₂");
  }

  #[test]
  fn test_mul_with_quadratic_form() {
    let algebra = clifford_algebra_non_euclidean();
    let e1 = algebra.blade([1]);
    let e01 = algebra.blade([0, 1]);
    let product = e1 * e01;
    assert_eq!(format!("{product}"), "-1e₀");
  }

  #[test]
  fn test_mul_euclidean() {
    let algebra = clifford_algebra_euclidean();
    let e1 = algebra.blade([1]);
    let e01 = algebra.blade([0, 1]);
    let product = e1 * e01;
    assert_eq!(format!("{product}"), "-1e₀");
  }

  #[test]
  fn test_mul_anti_commutativity() {
    let algebra = clifford_algebra_non_euclidean();
    let e0 = algebra.blade([0]);
    let e1 = algebra.blade([1]);
    let e2 = algebra.blade([2]);

    // Test anti-commutativity of basis vectors (clone to avoid move)
    assert_eq!(format!("{}", e0 * e1), "1e₀‚₁");
    assert_eq!(format!("{}", e1 * e0), "-1e₀‚₁");
    assert_eq!(format!("{}", e1 * e2), "1e₁‚₂");
    assert_eq!(format!("{}", e2 * e1), "-1e₁‚₂");
  }

  #[test]
  fn test_mul_scalar() {
    let algebra = clifford_algebra_non_euclidean();
    let one =
      algebra.element(SVector::<f64, 8>::from_row_slice(&[1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]));
    let e1 = algebra.blade([1]);
    let e2 = algebra.blade([2]);

    // Test scalar multiplication (clone to avoid move)
    assert_eq!(format!("{}", one * e1), "1e₁");
    assert_eq!(format!("{}", e1 * one), "1e₁");
    assert_eq!(format!("{}", one * e2), "1e₂");
    assert_eq!(format!("{}", e2 * one), "1e₂");
  }

  #[test]
  fn test_mul_quadratic_form_application() {
    let algebra = clifford_algebra_non_euclidean();
    let e0 = algebra.blade([0]);
    let e1 = algebra.blade([1]);
    let e2 = algebra.blade([2]);

    // Test quadratic form application
    assert_eq!(format!("{}", e0 * e0), "1"); // Q(e0) = 1
    assert_eq!(format!("{}", e1 * e1), "1"); // Q(e1) = 1
    assert_eq!(format!("{}", e2 * e2), "-1"); // Q(e2) = -1
  }

  #[test]
  fn test_mul_higher_grade() {
    let algebra = clifford_algebra_non_euclidean();
    let e01 = algebra.blade([0, 1]);
    let e12 = algebra.blade([1, 2]);
    let e02 = algebra.blade([0, 2]);

    // Test multiplication of bivectors
    assert_eq!(format!("{}", e01 * e12), "1e₀‚₂"); // e01 * e12 = e0 * e1 * e1 * e2 = e0 * Q(e1) * e2 = e0 * 1 * e2 = e0 * e2 = e02

    assert_eq!(format!("{}", e12 * e01), "1e₀‚₂");
    assert_eq!(format!("{}", e02 * e12), "1e₀‚₁"); // e02 * e12 = e0 * e2 * e1 * e2 = -e0 *
    // e1 * e2
    // * e2 = -e0 * e1 * Q(e2) = -e0 * e1 * -1 = e0
    // * e1 = e01
  }

  #[test]
  fn test_mul_trivector() {
    let algebra = clifford_algebra_non_euclidean();
    let e01 = algebra.blade([0, 1]);
    let e2 = algebra.blade([2]);
    let e012 = algebra.blade([0, 1, 2]);

    // Test multiplication with trivector
    assert_eq!(format!("{}", e01 * e2), "1e₀‚₁‚₂");
    assert_eq!(format!("{}", e2 * e01), "1e₀‚₁‚₂");
    assert_eq!(format!("{}", e012 * e2), "-1e₀‚₁"); // e012 * e2 = e0 * e1 * e2 * e2 = e0 *
    // e1 *
    // Q(e2) = e0 * e1 * -1 = -e0 * e1 =
    // -e01
  }

  /// Test that different quadratic forms cannot be mixed (compile-time safety).
  /// This test exists to document the expected behavior - mixing types should not compile.
  #[test]
  fn test_type_safety() {
    // These two algebras have different quadratic forms encoded in their types
    let euclidean = CliffordAlgebra::<f64, 3, Euclidean3D>::new();
    let non_euclidean = clifford_algebra_non_euclidean();

    let e1_euclidean = euclidean.blade([1]);
    let _e1_non_euclidean = non_euclidean.blade([1]);

    // The following would NOT compile (uncomment to verify):
    // let sum = e1_euclidean + e1_non_euclidean; // Error: mismatched types

    // But same-type operations work fine
    let e2_euclidean = euclidean.blade([2]);
    let _sum = e1_euclidean + e2_euclidean; // OK: same quadratic form
  }
}