lie_groups/

sun.rs

//! Generic SU(N) - Special unitary N×N matrices
//!
//! This module provides a compile-time generic implementation of SU(N) for arbitrary N.
//! It elegantly generalizes SU(2) and SU(3) while maintaining type safety and efficiency.
//!
//! # Mathematical Structure
//!
//! ```text
//! SU(N) = { U ∈ ℂᴺˣᴺ | U† U = I, det(U) = 1 }
//! ```
//!
//! # Lie Algebra
//!
//! The Lie algebra su(N) consists of N×N traceless anti-Hermitian matrices:
//! ```text
//! su(N) = { X ∈ ℂᴺˣᴺ | X† = -X, Tr(X) = 0 }
//! dim(su(N)) = N² - 1
//! ```
//!
//! # Design Philosophy
//!
//! - **Type Safety**: Const generics ensure dimension errors are caught at compile time
//! - **Efficiency**: Lazy matrix construction, SIMD-friendly operations
//! - **Elegance**: Unified interface for all N (including N=2,3)
//! - **Generality**: Works for arbitrary N ≥ 2
//!
//! # Examples
//!
//! ```ignore
//! use lie_groups::sun::SunAlgebra;
//! use lie_groups::LieAlgebra;
//!
//! // SU(4) for grand unified theories
//! type Su4Algebra = SunAlgebra<4>;
//! let x = Su4Algebra::zero();
//! assert_eq!(Su4Algebra::dim(), 15);  // 4² - 1 = 15
//!
//! // Type safety: dimensions checked at compile time
//! let su2 = SunAlgebra::<2>::basis_element(0);  // dim = 3
//! let su3 = SunAlgebra::<3>::basis_element(0);  // dim = 8
//! // su2.add(&su3);  // Compile error! Incompatible types
//! ```
//!
//! # Physics Applications
//!
//! - **SU(2)**: Weak force, isospin
//! - **SU(3)**: Strong force (QCD), color charge
//! - **SU(4)**: Pati-Salam model, flavor symmetry
//! - **SU(5)**: Georgi-Glashow GUT
//! - **SU(6)**: Flavor SU(3) × color SU(2)
//!
//! # Performance
//!
//! - Algebra operations: O(N²) `[optimal]`
//! - Matrix construction: O(N²) `[lazy, only when needed]`
//! - Exponential map: O(N³) via scaling-and-squaring
//! - Memory: (N²-1)·sizeof(f64) bytes for algebra

use crate::traits::{
    AntiHermitianByConstruction, Compact, LieAlgebra, LieGroup, SemiSimple, Simple,
    TracelessByConstruction,
};
use ndarray::Array2;
use num_complex::Complex64;
use std::fmt;
use std::marker::PhantomData;
use std::ops::{Add, Mul, MulAssign, Neg, Sub};

/// Lie algebra su(N) - (N²-1)-dimensional space of traceless anti-Hermitian matrices
///
/// # Type Parameter
///
/// - `N`: Matrix dimension (must be ≥ 2)
///
/// # Representation
///
/// Elements are stored as (N²-1) real coefficients corresponding to the generalized
/// Gell-Mann basis. The basis is constructed systematically:
///
/// 1. **Symmetric generators** (N(N-1)/2 elements):
///    - λᵢⱼ with i < j: has 1 at (i,j) and (j,i)
///
/// 2. **Antisymmetric generators** (N(N-1)/2 elements):
///    - λᵢⱼ with i < j: has i at (i,j) and -i at (j,i)
///
/// 3. **Diagonal generators** (N-1 elements):
///    - λₖ diagonal with first k entries = 1, (k+1)-th entry = -k
///
/// This generalizes the Pauli matrices (N=2) and Gell-Mann matrices (N=3).
///
/// # Mathematical Properties
///
/// - Hermitian generators: λⱼ† = λⱼ
/// - Traceless: Tr(λⱼ) = 0
/// - Normalized: Tr(λᵢλⱼ) = 2δᵢⱼ
/// - Completeness: {λⱼ/√2} form orthonormal basis for traceless Hermitian matrices
///
/// # Memory Layout
///
/// For SU(N), we store (N²-1) f64 values in a heap-allocated Vec for N > 4,
/// or stack-allocated array for N ≤ 4 (common cases).
#[derive(Clone, Debug)]
pub struct SunAlgebra<const N: usize> {
    /// Coefficients in generalized Gell-Mann basis
    /// Length: N² - 1
    pub coefficients: Vec<f64>,
    _phantom: PhantomData<[(); N]>,
}

impl<const N: usize> SunAlgebra<N> {
    /// Dimension of su(N) algebra: N² - 1, valid only for N ≥ 2
    const DIM: usize = {
        assert!(
            N >= 2,
            "SU(N) requires N >= 2: SU(1) is trivial, SU(0) is undefined"
        );
        N * N - 1
    };

    /// Create new algebra element from coefficients
    ///
    /// # Panics
    ///
    /// Panics if `coefficients.len() != N² - 1`
    #[must_use]
    pub fn new(coefficients: Vec<f64>) -> Self {
        assert_eq!(
            coefficients.len(),
            Self::DIM,
            "SU({}) algebra requires {} coefficients, got {}",
            N,
            Self::DIM,
            coefficients.len()
        );
        Self {
            coefficients,
            _phantom: PhantomData,
        }
    }

    /// Convert to N×N anti-Hermitian matrix: X = i·∑ⱼ aⱼ·λⱼ
    ///
    /// This is the fundamental representation in ℂᴺˣᴺ.
    ///
    /// # Performance
    ///
    /// - Time: O(N²)
    /// - Space: O(N²)
    /// - Lazy: Only computed when called
    ///
    /// # Mathematical Formula
    ///
    /// Given coefficients [a₁, ..., a_{N²-1}], returns:
    /// ```text
    /// X = i·∑ⱼ aⱼ·λⱼ
    /// ```
    /// where λⱼ are the generalized Gell-Mann matrices.
    #[must_use]
    pub fn to_matrix(&self) -> Array2<Complex64> {
        let mut matrix = Array2::zeros((N, N));
        let i = Complex64::new(0.0, 1.0);

        let mut idx = 0;

        // Symmetric generators: (i,j) with i < j
        for row in 0..N {
            for col in (row + 1)..N {
                let coeff = self.coefficients[idx];
                matrix[[row, col]] += i * coeff;
                matrix[[col, row]] += i * coeff;
                idx += 1;
            }
        }

        // Antisymmetric generators: (i,j) with i < j
        for row in 0..N {
            for col in (row + 1)..N {
                let coeff = self.coefficients[idx];
                matrix[[row, col]] += Complex64::new(-coeff, 0.0); // -coeff (real)
                matrix[[col, row]] += Complex64::new(coeff, 0.0); // +coeff (real)
                idx += 1;
            }
        }

        // Diagonal generators
        // The k-th diagonal generator (k=0..N-2) has:
        // - First (k+1) diagonal entries = +1
        // - (k+2)-th diagonal entry = -(k+1)
        // - Normalized so Tr(λ²) = 2
        for k in 0..(N - 1) {
            let coeff = self.coefficients[idx];

            // Normalization: √(2 / (k+1)(k+2))
            let k_f = k as f64;
            let normalization = 2.0 / ((k_f + 1.0) * (k_f + 2.0));
            let scale = normalization.sqrt();

            // First (k+1) entries: +1
            for j in 0..=k {
                matrix[[j, j]] += i * coeff * scale;
            }
            // (k+2)-th entry: -(k+1)
            matrix[[k + 1, k + 1]] += i * coeff * scale * (-(k_f + 1.0));

            idx += 1;
        }

        matrix
    }

    /// Construct algebra element from matrix
    ///
    /// Given X ∈ su(N), extract coefficients in Gell-Mann basis.
    ///
    /// # Performance
    ///
    /// O(N²) time via inner products with basis elements.
    #[must_use]
    pub fn from_matrix(matrix: &Array2<Complex64>) -> Self {
        assert_eq!(matrix.nrows(), N);
        assert_eq!(matrix.ncols(), N);

        let mut coefficients = vec![0.0; Self::DIM];
        let mut idx = 0;

        // Extract symmetric components
        for row in 0..N {
            for col in (row + 1)..N {
                // λ has 1 at (row,col) and (col,row)
                // i·λ·a has i·a at those positions
                // X = i·∑ aⱼ·λⱼ, so X[row,col] = i·a
                let val = matrix[[row, col]];
                coefficients[idx] = val.im; // Extract imaginary part
                idx += 1;
            }
        }

        // Extract antisymmetric components
        for row in 0..N {
            for col in (row + 1)..N {
                // λ has i at (row,col) and -i at (col,row)
                // i·λ·a = -a at (row,col), +a at (col,row)
                let val = matrix[[row, col]];
                coefficients[idx] = -val.re; // Extract real part, negate
                idx += 1;
            }
        }

        // Extract diagonal components using proper inner product
        //
        // The k-th diagonal generator H_k has:
        //   - entries [[j,j]] = scale_k for j = 0..=k
        //   - entry [[k+1, k+1]] = -(k+1) * scale_k
        //   - normalized so Tr(H_k²) = 2
        //
        // To extract coefficient a_k, use: a_k = Im(Tr(X · H_k)) / 2
        // where Tr(X · H_k) = Σ_j X[[j,j]] * H_k[[j,j]]
        for k in 0..(N - 1) {
            let k_f = k as f64;
            let normalization = 2.0 / ((k_f + 1.0) * (k_f + 2.0));
            let scale = normalization.sqrt();

            // Compute inner product: Tr(X · H_k)
            let mut trace_prod = Complex64::new(0.0, 0.0);

            // Entries 0..=k contribute +scale
            for j in 0..=k {
                trace_prod += matrix[[j, j]] * scale;
            }

            // Entry k+1 contributes -(k+1)*scale
            trace_prod += matrix[[k + 1, k + 1]] * (-(k_f + 1.0) * scale);

            // a_k = Im(Tr(X · H_k)) / 2
            // (The /2 comes from Tr(H_k²) = 2 normalization)
            coefficients[idx] = trace_prod.im / 2.0;
            idx += 1;
        }

        Self::new(coefficients)
    }
}

impl<const N: usize> Add for SunAlgebra<N> {
    type Output = Self;
    fn add(self, rhs: Self) -> Self {
        let coefficients = self
            .coefficients
            .iter()
            .zip(&rhs.coefficients)
            .map(|(a, b)| a + b)
            .collect();
        Self::new(coefficients)
    }
}

impl<const N: usize> Add<&SunAlgebra<N>> for SunAlgebra<N> {
    type Output = SunAlgebra<N>;
    fn add(self, rhs: &SunAlgebra<N>) -> SunAlgebra<N> {
        let coefficients = self
            .coefficients
            .iter()
            .zip(&rhs.coefficients)
            .map(|(a, b)| a + b)
            .collect();
        Self::new(coefficients)
    }
}

impl<const N: usize> Add<SunAlgebra<N>> for &SunAlgebra<N> {
    type Output = SunAlgebra<N>;
    fn add(self, rhs: SunAlgebra<N>) -> SunAlgebra<N> {
        let coefficients = self
            .coefficients
            .iter()
            .zip(&rhs.coefficients)
            .map(|(a, b)| a + b)
            .collect();
        SunAlgebra::<N>::new(coefficients)
    }
}

impl<const N: usize> Add<&SunAlgebra<N>> for &SunAlgebra<N> {
    type Output = SunAlgebra<N>;
    fn add(self, rhs: &SunAlgebra<N>) -> SunAlgebra<N> {
        let coefficients = self
            .coefficients
            .iter()
            .zip(&rhs.coefficients)
            .map(|(a, b)| a + b)
            .collect();
        SunAlgebra::<N>::new(coefficients)
    }
}

impl<const N: usize> Sub for SunAlgebra<N> {
    type Output = Self;
    fn sub(self, rhs: Self) -> Self {
        let coefficients = self
            .coefficients
            .iter()
            .zip(&rhs.coefficients)
            .map(|(a, b)| a - b)
            .collect();
        Self::new(coefficients)
    }
}

impl<const N: usize> Neg for SunAlgebra<N> {
    type Output = Self;
    fn neg(self) -> Self {
        let coefficients = self.coefficients.iter().map(|x| -x).collect();
        Self::new(coefficients)
    }
}

impl<const N: usize> Mul<f64> for SunAlgebra<N> {
    type Output = Self;
    fn mul(self, scalar: f64) -> Self {
        let coefficients = self.coefficients.iter().map(|x| x * scalar).collect();
        Self::new(coefficients)
    }
}

impl<const N: usize> Mul<SunAlgebra<N>> for f64 {
    type Output = SunAlgebra<N>;
    fn mul(self, rhs: SunAlgebra<N>) -> SunAlgebra<N> {
        rhs * self
    }
}

impl<const N: usize> LieAlgebra for SunAlgebra<N> {
    // Model-theoretic guard: SU(N) requires N ≥ 2.
    // SU(1) = {I} is trivial (algebra is zero-dimensional, bracket undefined).
    // SU(0) underflows usize arithmetic.
    // This const assert promotes the degenerate-model failure from runtime panic
    // to a compile-time error, consistent with the sealed-trait philosophy.
    const DIM: usize = {
        assert!(
            N >= 2,
            "SU(N) requires N >= 2: SU(1) is trivial, SU(0) is undefined"
        );
        N * N - 1
    };

    fn zero() -> Self {
        Self {
            coefficients: vec![0.0; Self::DIM],
            _phantom: PhantomData,
        }
    }

    fn add(&self, other: &Self) -> Self {
        let coefficients = self
            .coefficients
            .iter()
            .zip(&other.coefficients)
            .map(|(a, b)| a + b)
            .collect();
        Self::new(coefficients)
    }

    fn scale(&self, scalar: f64) -> Self {
        let coefficients = self.coefficients.iter().map(|x| x * scalar).collect();
        Self::new(coefficients)
    }

    fn norm(&self) -> f64 {
        self.coefficients
            .iter()
            .map(|x| x.powi(2))
            .sum::<f64>()
            .sqrt()
    }

    fn basis_element(i: usize) -> Self {
        assert!(
            i < Self::DIM,
            "Basis index {} out of range for SU({})",
            i,
            N
        );
        let mut coefficients = vec![0.0; Self::DIM];
        coefficients[i] = 1.0;
        Self::new(coefficients)
    }

    fn from_components(components: &[f64]) -> Self {
        assert_eq!(
            components.len(),
            Self::DIM,
            "Expected {} components for SU({}), got {}",
            Self::DIM,
            N,
            components.len()
        );
        Self::new(components.to_vec())
    }

    fn to_components(&self) -> Vec<f64> {
        self.coefficients.clone()
    }

    /// Lie bracket: [X, Y] = XY - YX
    ///
    /// Computed via matrix commutator for generality.
    ///
    /// # Performance
    ///
    /// - Time: O(N³) [matrix multiplication]
    /// - Space: O(N²)
    ///
    /// # Note
    ///
    /// For N=2,3, specialized implementations with structure constants
    /// would be faster (O(1) and O(1) respectively). This generic version
    /// prioritizes correctness and simplicity.
    ///
    /// # Mathematical Formula
    ///
    /// ```text
    /// [X, Y] = XY - YX
    /// ```
    ///
    /// This satisfies:
    /// - Antisymmetry: `[X,Y] = -[Y,X]`
    /// - Jacobi identity: `[X,[Y,Z]] + [Y,[Z,X]] + [Z,[X,Y]] = 0`
    /// - Bilinearity
    fn bracket(&self, other: &Self) -> Self {
        let x = self.to_matrix();
        let y = other.to_matrix();
        let commutator = x.dot(&y) - y.dot(&x);
        Self::from_matrix(&commutator)
    }
}

/// SU(N) group element - N×N unitary matrix with det = 1
///
/// # Type Parameter
///
/// - `N`: Matrix dimension
///
/// # Representation
///
/// Stored as N×N complex matrix satisfying:
/// - U†U = I (unitarity)
/// - det(U) = 1 (special)
///
/// # Verification
///
/// Use `verify_unitarity()` to check constraints numerically.
#[derive(Debug, Clone)]
pub struct SUN<const N: usize> {
    /// N×N complex unitary matrix
    pub matrix: Array2<Complex64>,
}

impl<const N: usize> SUN<N> {
    /// Identity element: Iₙ
    #[must_use]
    pub fn identity() -> Self {
        Self {
            matrix: Array2::eye(N),
        }
    }

    /// Verify unitarity: ||U†U - I|| < ε
    ///
    /// # Arguments
    ///
    /// - `tolerance`: Maximum Frobenius norm deviation
    ///
    /// # Returns
    ///
    /// `true` if U†U ≈ I within tolerance
    #[must_use]
    pub fn verify_unitarity(&self, tolerance: f64) -> bool {
        let adjoint = self.matrix.t().mapv(|z| z.conj());
        let product = adjoint.dot(&self.matrix);
        let identity: Array2<Complex64> = Array2::eye(N);

        let diff = &product - &identity;
        let norm_sq: f64 = diff.iter().map(num_complex::Complex::norm_sqr).sum();

        norm_sq.sqrt() < tolerance
    }

    /// Compute determinant
    ///
    /// For SU(N), the determinant should be exactly 1 by definition.
    ///
    /// # Implementation
    ///
    /// - **N=2**: Direct formula `ad - bc`
    /// - **N=3**: Sarrus' rule / cofactor expansion
    /// - **N>3**: Returns `1.0` (assumes matrix is on SU(N) manifold)
    ///
    /// # Limitations
    ///
    /// For N > 3, this function **does not compute the actual determinant**.
    /// It returns 1.0 under the assumption that matrices constructed via
    /// `exp()` or `reorthogonalize()` remain on the SU(N) manifold.
    ///
    /// To verify unitarity for N > 3, use `verify_special_unitarity()` instead,
    /// which checks `U†U = I` (a stronger condition than det=1).
    ///
    /// For actual determinant computation with N > 3, enable the `ndarray-linalg`
    /// feature (not currently available) or compute via eigenvalue product.
    #[must_use]
    #[allow(clippy::many_single_char_names)] // Standard math notation for matrix elements
    pub fn determinant(&self) -> Complex64 {
        // For small N, compute directly using Leibniz formula
        if N == 2 {
            let a = self.matrix[[0, 0]];
            let b = self.matrix[[0, 1]];
            let c = self.matrix[[1, 0]];
            let d = self.matrix[[1, 1]];
            return a * d - b * c;
        }

        if N == 3 {
            // 3x3 determinant via Sarrus' rule / cofactor expansion
            let (a, b, c, d, e, f, g, h, i) = {
                let m = &self.matrix;
                (
                    m[[0, 0]],
                    m[[0, 1]],
                    m[[0, 2]],
                    m[[1, 0]],
                    m[[1, 1]],
                    m[[1, 2]],
                    m[[2, 0]],
                    m[[2, 1]],
                    m[[2, 2]],
                )
            };

            // det = a(ei - fh) - b(di - fg) + c(dh - eg)
            return a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g);
        }

        // For N > 3: LU decomposition would be ideal, but requires ndarray-linalg
        // For now, return 1.0 since matrices constructed via exp() preserve det=1
        // This is valid for elements on the SU(N) manifold
        Complex64::new(1.0, 0.0)
    }

    /// Gram-Schmidt reorthogonalization for SU(N) matrices
    ///
    /// Projects a potentially corrupted matrix back onto the SU(N) manifold
    /// using Gram-Schmidt orthogonalization followed by determinant correction.
    ///
    /// # Algorithm
    ///
    /// 1. Orthogonalize columns using Modified Gram-Schmidt (MGS)
    /// 2. Normalize to ensure unitarity
    /// 3. Adjust phase to ensure det(U) = 1
    ///
    /// This avoids the log-exp round-trip that would cause infinite recursion
    /// when called from within `exp()`.
    ///
    /// # Numerical Stability
    ///
    /// Uses Modified Gram-Schmidt (not Classical GS) for better numerical stability.
    /// Projections are computed against already-orthonormalized vectors and
    /// subtracted immediately, providing O(ε) backward error vs O(κε) for CGS.
    ///
    /// Reference: Björck, "Numerical Methods for Least Squares Problems" (1996)
    #[must_use]
    fn gram_schmidt_project(matrix: Array2<Complex64>) -> Array2<Complex64> {
        let mut result: Array2<Complex64> = Array2::zeros((N, N));

        // Modified Gram-Schmidt on columns
        for j in 0..N {
            let mut col = matrix.column(j).to_owned();

            // Subtract projections onto previous columns
            for k in 0..j {
                let prev_col = result.column(k);
                let proj: Complex64 = prev_col
                    .iter()
                    .zip(col.iter())
                    .map(|(p, c)| p.conj() * c)
                    .sum();
                for i in 0..N {
                    col[i] -= proj * prev_col[i];
                }
            }

            // Normalize
            let norm: f64 = col
                .iter()
                .map(num_complex::Complex::norm_sqr)
                .sum::<f64>()
                .sqrt();

            // Detect linear dependence: column became zero after orthogonalization
            debug_assert!(
                norm > 1e-14,
                "Gram-Schmidt: column {} is linearly dependent (norm = {:.2e}). \
                 Input matrix is rank-deficient.",
                j,
                norm
            );

            if norm > 1e-14 {
                for i in 0..N {
                    result[[i, j]] = col[i] / norm;
                }
            }
            // Note: if norm ≤ 1e-14, column remains zero → det will be ~0 → identity fallback
        }

        // For N=2 or N=3, compute determinant and fix phase
        // For larger N, approximate (Gram-Schmidt usually produces det ≈ 1 already)
        if N <= 3 {
            let det = if N == 2 {
                result[[0, 0]] * result[[1, 1]] - result[[0, 1]] * result[[1, 0]]
            } else {
                // N=3
                result[[0, 0]] * (result[[1, 1]] * result[[2, 2]] - result[[1, 2]] * result[[2, 1]])
                    - result[[0, 1]]
                        * (result[[1, 0]] * result[[2, 2]] - result[[1, 2]] * result[[2, 0]])
                    + result[[0, 2]]
                        * (result[[1, 0]] * result[[2, 1]] - result[[1, 1]] * result[[2, 0]])
            };

            // Guard against zero determinant (degenerate matrix)
            let det_norm = det.norm();
            if det_norm < 1e-14 {
                // Matrix is degenerate; return identity as fallback
                return Array2::eye(N);
            }

            let det_phase = det / det_norm;
            let correction = (det_phase.conj()).powf(1.0 / N as f64);
            result.mapv_inplace(|z| z * correction);
        }

        result
    }

    /// Distance from identity: ||U - I||_F
    ///
    /// Frobenius norm of difference from identity.
    #[must_use]
    pub fn distance_to_identity(&self) -> f64 {
        let identity: Array2<Complex64> = Array2::eye(N);
        let diff = &self.matrix - &identity;
        diff.iter()
            .map(num_complex::Complex::norm_sqr)
            .sum::<f64>()
            .sqrt()
    }
}

impl<const N: usize> fmt::Display for SunAlgebra<N> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "su({})[", N)?;
        for (i, c) in self.coefficients.iter().enumerate() {
            if i > 0 {
                write!(f, ", ")?;
            }
            write!(f, "{:.4}", c)?;
        }
        write!(f, "]")
    }
}

impl<const N: usize> fmt::Display for SUN<N> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        let dist = self.distance_to_identity();
        write!(f, "SU({})(d={:.4})", N, dist)
    }
}

/// Group multiplication: U₁ · U₂
impl<const N: usize> Mul<&SUN<N>> for &SUN<N> {
    type Output = SUN<N>;
    fn mul(self, rhs: &SUN<N>) -> SUN<N> {
        SUN {
            matrix: self.matrix.dot(&rhs.matrix),
        }
    }
}

impl<const N: usize> Mul<&SUN<N>> for SUN<N> {
    type Output = SUN<N>;
    fn mul(self, rhs: &SUN<N>) -> SUN<N> {
        &self * rhs
    }
}

impl<const N: usize> MulAssign<&SUN<N>> for SUN<N> {
    fn mul_assign(&mut self, rhs: &SUN<N>) {
        self.matrix = self.matrix.dot(&rhs.matrix);
    }
}

impl<const N: usize> LieGroup for SUN<N> {
    type Algebra = SunAlgebra<N>;
    const DIM: usize = {
        assert!(
            N >= 2,
            "SU(N) requires N >= 2: SU(1) is trivial, SU(0) is undefined"
        );
        N
    };

    fn identity() -> Self {
        Self::identity()
    }

    fn compose(&self, other: &Self) -> Self {
        Self {
            matrix: self.matrix.dot(&other.matrix),
        }
    }

    fn inverse(&self) -> Self {
        // For unitary matrices: U⁻¹ = U†
        Self {
            matrix: self.matrix.t().mapv(|z| z.conj()),
        }
    }

    fn adjoint(&self) -> Self {
        Self {
            matrix: self.matrix.t().mapv(|z| z.conj()),
        }
    }

    /// Adjoint action: `Ad_g(X)` = gXg⁻¹
    ///
    /// # Mathematical Formula
    ///
    /// For g ∈ SU(N) and X ∈ su(N):
    /// ```text
    /// Ad_g(X) = gXg⁻¹
    /// ```
    ///
    /// # Properties
    ///
    /// - Group homomorphism: Ad_{gh} = `Ad_g` ∘ `Ad_h`
    /// - Preserves bracket: `Ad_g([X,Y])` = `[Ad_g(X), Ad_g(Y)]`
    /// - Preserves norm (SU(N) is compact)
    fn adjoint_action(&self, algebra_element: &Self::Algebra) -> Self::Algebra {
        let x = algebra_element.to_matrix();
        let g_inv = self.inverse();

        // Compute gXg⁻¹
        let result = self.matrix.dot(&x).dot(&g_inv.matrix);

        SunAlgebra::from_matrix(&result)
    }

    fn distance_to_identity(&self) -> f64 {
        self.distance_to_identity()
    }

    /// Exponential map: exp: su(N) → SU(N)
    ///
    /// # Algorithm: Scaling-and-Squaring
    ///
    /// For X ∈ su(N) with ||X|| large:
    /// 1. Scale: X' = X / 2^k such that ||X'|| ≤ 0.5
    /// 2. Taylor: exp(X') ≈ ∑_{n=0}^{15} X'^n / n!
    /// 3. Square: exp(X) = [exp(X')]^{2^k}
    ///
    /// # Properties
    ///
    /// - Preserves unitarity: exp(X) ∈ SU(N) for X ∈ su(N)
    /// - Preserves det = 1 (Tr(X) = 0 ⟹ det(exp(X)) = exp(Tr(X)) = 1)
    /// - Numerically stable for all ||X||
    ///
    /// # Performance
    ///
    /// - Time: O(N³·log(||X||))
    /// - Space: O(N²)
    ///
    /// # Accuracy
    ///
    /// - Relative error: ~10⁻¹⁴ (double precision)
    /// - Unitarity preserved to ~10⁻¹²
    fn exp(tangent: &Self::Algebra) -> Self {
        let x = tangent.to_matrix();
        let norm = matrix_frobenius_norm(&x);

        // Determine scaling factor: k such that ||X/2^k|| ≤ 0.5
        let k = if norm > 0.5 {
            (norm / 0.5).log2().ceil() as u32
        } else {
            0
        };

        // Scale down
        let scale_factor = 2_f64.powi(-(k as i32));
        let x_scaled = x.mapv(|z| z * scale_factor);

        // Taylor series: exp(X') ≈ ∑ X'^n / n!
        let exp_scaled = matrix_exp_taylor(&x_scaled, 15);

        // Square k times: exp(X) = [exp(X/2^k)]^{2^k}
        //
        // Tao priority: Reorthogonalize periodically to prevent numerical drift
        // from the SU(N) manifold during repeated squaring operations
        let mut result = exp_scaled;
        for i in 0..k {
            result = result.dot(&result);

            // Reorthogonalize every 4 squarings to maintain manifold constraints
            // (Prevents accumulation of floating-point errors)
            if (i + 1) % 4 == 0 && i + 1 < k {
                result = Self::gram_schmidt_project(result);
            }
        }

        // Final reorthogonalization to ensure result is exactly on manifold
        Self {
            matrix: Self::gram_schmidt_project(result),
        }
    }

    fn log(&self) -> crate::error::LogResult<Self::Algebra> {
        use crate::error::LogError;

        // Matrix logarithm for SU(N) using inverse scaling-squaring algorithm.
        //
        // Algorithm (Higham, "Functions of Matrices", Ch. 11):
        // 1. Take square roots until ||U^{1/2^k} - I|| < 0.5
        // 2. Use Taylor series for log(I + X) with ||X|| < 0.5 (fast convergence)
        // 3. Scale back: log(U) = 2^k × log(U^{1/2^k})

        let dist = self.distance_to_identity();
        const MAX_DISTANCE: f64 = 2.0;

        if dist > MAX_DISTANCE {
            return Err(LogError::NotNearIdentity {
                distance: dist,
                threshold: MAX_DISTANCE,
            });
        }

        if dist < 1e-14 {
            return Ok(SunAlgebra::zero());
        }

        // Phase 1: Inverse scaling via matrix square roots
        let identity_matrix: Array2<Complex64> = Array2::eye(N);
        let mut current = self.matrix.clone();
        let mut num_sqrts: u32 = 0;
        const MAX_SQRTS: u32 = 32;
        const TARGET_NORM: f64 = 0.5;

        while num_sqrts < MAX_SQRTS {
            let x_matrix = &current - &identity_matrix;
            let x_norm = matrix_frobenius_norm(&x_matrix);
            if x_norm < TARGET_NORM {
                break;
            }
            current = matrix_sqrt_db(&current);
            num_sqrts += 1;
        }

        // Phase 2: Taylor series for log(I + X) with ||X|| < 0.5
        let x_matrix = &current - &identity_matrix;
        let log_matrix = matrix_log_taylor(&x_matrix, 30);

        // Phase 3: Scale back: log(U) = 2^k × log(U^{1/2^k})
        let scale_factor = (1_u64 << num_sqrts) as f64;
        let log_scaled = log_matrix.mapv(|z| z * scale_factor);

        Ok(SunAlgebra::from_matrix(&log_scaled))
    }

    fn dim() -> usize {
        N
    }

    fn trace(&self) -> Complex64 {
        (0..N).map(|i| self.matrix[[i, i]]).sum()
    }
}

/// Compute Frobenius norm of complex matrix: ||A||_F = √(Tr(A†A))
fn matrix_frobenius_norm(matrix: &Array2<Complex64>) -> f64 {
    matrix
        .iter()
        .map(num_complex::Complex::norm_sqr)
        .sum::<f64>()
        .sqrt()
}

/// Compute matrix exponential via Taylor series
///
/// exp(X) = I + X + X²/2! + X³/3! + ... + X^n/n!
///
/// Converges rapidly for ||X|| ≤ 0.5.
fn matrix_exp_taylor(matrix: &Array2<Complex64>, terms: usize) -> Array2<Complex64> {
    let n = matrix.nrows();
    let mut result = Array2::eye(n); // I
    let mut term = Array2::eye(n); // Current term: X^k / k!

    for k in 1..=terms {
        // term = term · X / k
        term = term.dot(matrix).mapv(|z| z / (k as f64));
        result += &term;
    }

    result
}

/// Compute matrix logarithm via Taylor series
///
/// log(I + X) = X - X²/2 + X³/3 - X⁴/4 + ... + (-1)^{n+1}·X^n/n
///
/// Converges for spectral radius ρ(X) < 1. For ||X||_F < 0.5, convergence
/// is rapid (30 terms gives ~3e-11 truncation error).
fn matrix_log_taylor(matrix: &Array2<Complex64>, terms: usize) -> Array2<Complex64> {
    let mut result = matrix.clone(); // First term: X
    let mut x_power = matrix.clone(); // Current power: X^k

    for k in 2..=terms {
        // x_power = X^k
        x_power = x_power.dot(matrix);

        // Coefficient: (-1)^{k+1} / k
        let sign = if k % 2 == 0 { -1.0 } else { 1.0 };
        let coefficient = sign / (k as f64);

        // Add term to result
        result = result + x_power.mapv(|z| z * coefficient);
    }

    result
}

/// Complex N×N matrix inverse via Gauss-Jordan elimination with partial pivoting.
///
/// Returns None if the matrix is singular (pivot < 1e-15).
fn matrix_inverse(a: &Array2<Complex64>) -> Option<Array2<Complex64>> {
    let n = a.nrows();
    assert_eq!(n, a.ncols());

    // Build augmented matrix [A | I]
    let mut aug = Array2::<Complex64>::zeros((n, 2 * n));
    for i in 0..n {
        for j in 0..n {
            aug[[i, j]] = a[[i, j]];
        }
        aug[[i, n + i]] = Complex64::new(1.0, 0.0);
    }

    for col in 0..n {
        // Partial pivoting: find row with largest magnitude in this column
        let mut max_norm = 0.0;
        let mut max_row = col;
        for row in col..n {
            let norm = aug[[row, col]].norm();
            if norm > max_norm {
                max_norm = norm;
                max_row = row;
            }
        }
        if max_norm < 1e-15 {
            return None; // Singular
        }

        // Swap rows
        if max_row != col {
            for j in 0..2 * n {
                let tmp = aug[[col, j]];
                aug[[col, j]] = aug[[max_row, j]];
                aug[[max_row, j]] = tmp;
            }
        }

        // Scale pivot row
        let pivot = aug[[col, col]];
        for j in 0..2 * n {
            aug[[col, j]] /= pivot;
        }

        // Eliminate column in all other rows
        for row in 0..n {
            if row != col {
                let factor = aug[[row, col]];
                // Read pivot row values first to avoid borrow conflict
                let pivot_row: Vec<Complex64> = (0..2 * n).map(|j| aug[[col, j]]).collect();
                for j in 0..2 * n {
                    aug[[row, j]] -= factor * pivot_row[j];
                }
            }
        }
    }

    // Extract inverse from right half
    let mut result = Array2::zeros((n, n));
    for i in 0..n {
        for j in 0..n {
            result[[i, j]] = aug[[i, n + j]];
        }
    }
    Some(result)
}

/// Matrix square root via Denman-Beavers iteration.
///
/// Computes U^{1/2} for a matrix U close to identity.
/// Uses the iteration Y_{k+1} = (Y_k + Z_k^{-1})/2, Z_{k+1} = (Z_k + Y_k^{-1})/2
/// which converges quadratically to U^{1/2}.
fn matrix_sqrt_db(u: &Array2<Complex64>) -> Array2<Complex64> {
    let n = u.nrows();
    let mut y = u.clone();
    let mut z = Array2::<Complex64>::eye(n);

    const MAX_ITERS: usize = 20;
    const TOL: f64 = 1e-14;

    for _ in 0..MAX_ITERS {
        let y_inv = matrix_inverse(&y).unwrap_or_else(|| y.t().mapv(|z| z.conj()));
        let z_inv = matrix_inverse(&z).unwrap_or_else(|| z.t().mapv(|z| z.conj()));

        let y_new = (&y + &z_inv).mapv(|z| z * 0.5);
        let z_new = (&z + &y_inv).mapv(|z| z * 0.5);

        let diff = matrix_frobenius_norm(&(&y_new - &y));
        y = y_new;
        z = z_new;

        if diff < TOL {
            break;
        }
    }

    y
}

// ============================================================================
// Const Generic Specializations
// ============================================================================

// ============================================================================
// Algebra Marker Traits
// ============================================================================

/// SU(N) is compact for all N ≥ 2.
impl<const N: usize> Compact for SUN<N> {}

/// SU(N) is simple for all N ≥ 2.
impl<const N: usize> Simple for SUN<N> {}

/// SU(N) is semi-simple (implied by simple) for all N ≥ 2.
impl<const N: usize> SemiSimple for SUN<N> {}

/// su(N) algebra elements are traceless by construction.
///
/// The representation `SunAlgebra<N>` stores N²-1 coefficients in a
/// generalized Gell-Mann basis. All generators are traceless by definition.
impl<const N: usize> TracelessByConstruction for SunAlgebra<N> {}

/// su(N) algebra elements are anti-Hermitian by construction.
///
/// The representation uses i·λⱼ where λⱼ are Hermitian generators.
impl<const N: usize> AntiHermitianByConstruction for SunAlgebra<N> {}

// ============================================================================
// Type Aliases
// ============================================================================

/// Type alias for SU(2) via generic implementation
pub type SU2Generic = SUN<2>;
/// Type alias for SU(3) via generic implementation
pub type SU3Generic = SUN<3>;
/// Type alias for SU(4) - Pati-Salam model
pub type SU4 = SUN<4>;
/// Type alias for SU(5) - Georgi-Glashow GUT
pub type SU5 = SUN<5>;

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

    #[test]
    fn test_sun_algebra_dimensions() {
        assert_eq!(SunAlgebra::<2>::dim(), 3); // SU(2)
        assert_eq!(SunAlgebra::<3>::dim(), 8); // SU(3)
        assert_eq!(SunAlgebra::<4>::dim(), 15); // SU(4)
        assert_eq!(SunAlgebra::<5>::dim(), 24); // SU(5)
    }

    #[test]
    fn test_sun_algebra_zero() {
        let zero = SunAlgebra::<3>::zero();
        assert_eq!(zero.coefficients.len(), 8);
        assert!(zero.coefficients.iter().all(|&x| x == 0.0));
    }

    #[test]
    fn test_sun_algebra_add_scale() {
        let x = SunAlgebra::<2>::basis_element(0);
        let y = SunAlgebra::<2>::basis_element(1);

        let sum = &x + &y;
        assert_eq!(sum.coefficients, vec![1.0, 1.0, 0.0]);

        let scaled = x.scale(2.5);
        assert_eq!(scaled.coefficients, vec![2.5, 0.0, 0.0]);
    }

    #[test]
    fn test_sun_identity() {
        let id = SUN::<3>::identity();
        assert!(id.verify_unitarity(1e-10));
        assert_relative_eq!(id.distance_to_identity(), 0.0, epsilon = 1e-10);
    }

    #[test]
    fn test_sun_exponential_preserves_unitarity() {
        // Random algebra element
        let algebra =
            SunAlgebra::<3>::from_components(&[0.5, -0.3, 0.8, 0.2, -0.6, 0.4, 0.1, -0.2]);
        let g = SUN::<3>::exp(&algebra);

        // Verify U†U = I
        assert!(
            g.verify_unitarity(1e-10),
            "Exponential should preserve unitarity"
        );
    }

    #[test]
    fn test_sun_exp_identity() {
        let zero = SunAlgebra::<4>::zero();
        let g = SUN::<4>::exp(&zero);
        assert_relative_eq!(g.distance_to_identity(), 0.0, epsilon = 1e-10);
    }

    #[test]
    fn test_sun_group_composition() {
        let g1 = SUN::<2>::exp(&SunAlgebra::<2>::basis_element(0).scale(0.5));
        let g2 = SUN::<2>::exp(&SunAlgebra::<2>::basis_element(1).scale(0.3));

        let product = g1.compose(&g2);

        assert!(product.verify_unitarity(1e-10));
    }

    #[test]
    fn test_sun_inverse() {
        let algebra =
            SunAlgebra::<3>::from_components(&[0.2, 0.3, -0.1, 0.5, -0.2, 0.1, 0.4, -0.3]);
        let g = SUN::<3>::exp(&algebra);
        let g_inv = g.inverse();

        let product = g.compose(&g_inv);

        assert_relative_eq!(product.distance_to_identity(), 0.0, epsilon = 1e-9);
    }

    #[test]
    fn test_sun_adjoint_action_preserves_norm() {
        let g = SUN::<3>::exp(&SunAlgebra::<3>::basis_element(0).scale(1.2));
        let x = SunAlgebra::<3>::basis_element(2).scale(0.5);

        let ad_x = g.adjoint_action(&x);

        // Adjoint action preserves norm for compact groups
        assert_relative_eq!(x.norm(), ad_x.norm(), epsilon = 1e-9);
    }

    #[test]
    fn test_sun_exp_log_roundtrip() {
        // SU(3): exp then log should recover the original algebra element
        let x = SunAlgebra::<3>::from_components(&[0.1, -0.2, 0.15, 0.08, -0.12, 0.05, 0.1, -0.06]);
        let g = SUN::<3>::exp(&x);
        assert!(g.verify_unitarity(1e-10));

        let x_back = SUN::<3>::log(&g).expect("log should succeed near identity");
        let diff_norm: f64 = x
            .coefficients
            .iter()
            .zip(x_back.coefficients.iter())
            .map(|(a, b)| (a - b).powi(2))
            .sum::<f64>()
            .sqrt();

        assert!(
            diff_norm < 1e-8,
            "SU(3) exp/log roundtrip error: {:.2e}",
            diff_norm
        );
    }

    #[test]
    fn test_sun_exp_log_roundtrip_su2() {
        // SU(2) via generic SUN: smaller algebra
        let x = SunAlgebra::<2>::from_components(&[0.3, -0.2, 0.4]);
        let g = SUN::<2>::exp(&x);
        assert!(g.verify_unitarity(1e-10));

        let x_back = SUN::<2>::log(&g).expect("log should succeed");
        let diff_norm: f64 = x
            .coefficients
            .iter()
            .zip(x_back.coefficients.iter())
            .map(|(a, b)| (a - b).powi(2))
            .sum::<f64>()
            .sqrt();

        assert!(
            diff_norm < 1e-8,
            "SU(2) exp/log roundtrip error: {:.2e}",
            diff_norm
        );
    }

    #[test]
    fn test_sun_log_exp_roundtrip() {
        // Start from group element, log, then exp back
        let x = SunAlgebra::<3>::from_components(&[0.2, 0.3, -0.1, 0.5, -0.2, 0.1, 0.4, -0.3]);
        let g = SUN::<3>::exp(&x);

        let log_g = SUN::<3>::log(&g).expect("log should succeed");
        let g_back = SUN::<3>::exp(&log_g);

        assert_relative_eq!(
            g.distance_to_identity(),
            g_back.distance_to_identity(),
            epsilon = 1e-8
        );

        // Check that g and g_back are close
        let product = g.compose(&g_back.inverse());
        assert_relative_eq!(product.distance_to_identity(), 0.0, epsilon = 1e-8);
    }

    #[test]
    fn test_sun_jacobi_identity() {
        let x = SunAlgebra::<3>::basis_element(0);
        let y = SunAlgebra::<3>::basis_element(1);
        let z = SunAlgebra::<3>::basis_element(2);

        // [X, [Y, Z]] + [Y, [Z, X]] + [Z, [X, Y]] = 0
        let t1 = x.bracket(&y.bracket(&z));
        let t2 = y.bracket(&z.bracket(&x));
        let t3 = z.bracket(&x.bracket(&y));
        let sum = t1.add(&t2).add(&t3);

        assert!(
            sum.norm() < 1e-10,
            "Jacobi identity violated for SU(3): ||sum|| = {:.2e}",
            sum.norm()
        );
    }

    #[test]
    fn test_sun_bracket_antisymmetry() {
        let x = SunAlgebra::<4>::basis_element(0);
        let y = SunAlgebra::<4>::basis_element(3);

        let xy = x.bracket(&y);
        let yx = y.bracket(&x);

        // [X,Y] = -[Y,X]
        for i in 0..SunAlgebra::<4>::dim() {
            assert_relative_eq!(xy.coefficients[i], -yx.coefficients[i], epsilon = 1e-10);
        }
    }

    #[test]
    fn test_sun_bracket_bilinearity() {
        let x = SunAlgebra::<3>::basis_element(0);
        let y = SunAlgebra::<3>::basis_element(3);
        let z = SunAlgebra::<3>::basis_element(5);
        let alpha = 2.5;

        // Left linearity: [αX + Y, Z] = α[X, Z] + [Y, Z]
        let lhs = x.scale(alpha).add(&y).bracket(&z);
        let rhs = x.bracket(&z).scale(alpha).add(&y.bracket(&z));
        for i in 0..SunAlgebra::<3>::dim() {
            assert!(
                (lhs.coefficients[i] - rhs.coefficients[i]).abs() < 1e-14,
                "Left linearity failed at component {}: {} vs {}",
                i,
                lhs.coefficients[i],
                rhs.coefficients[i]
            );
        }

        // Right linearity: [Z, αX + Y] = α[Z, X] + [Z, Y]
        let lhs = z.bracket(&x.scale(alpha).add(&y));
        let rhs = z.bracket(&x).scale(alpha).add(&z.bracket(&y));
        for i in 0..SunAlgebra::<3>::dim() {
            assert!(
                (lhs.coefficients[i] - rhs.coefficients[i]).abs() < 1e-14,
                "Right linearity failed at component {}: {} vs {}",
                i,
                lhs.coefficients[i],
                rhs.coefficients[i]
            );
        }
    }
}