diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f208129..402a7a6 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -71,6 +71,7 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@11bd71901bbe5b1630ceea73d27597364c9af683 # v4 + - run: cargo generate-lockfile - uses: rustsec/audit-check@69366f33c96575abad1ee0dba8212993eecbe998 # v2.0.0 with: token: ${{ secrets.GITHUB_TOKEN }} @@ -90,4 +91,4 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@11bd71901bbe5b1630ceea73d27597364c9af683 # v4 - - uses: EmbarkStudios/cargo-deny-action@34899fc7ba81ca6268d5947a7a16b4649013fea1 # v2 + - uses: EmbarkStudios/cargo-deny-action@44db170f6a7d12a6e90340e9e0fca1f650d34b14 # v2.0.15 diff --git a/Cargo.toml b/Cargo.toml index 5c9e07c..2b0a90d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "lie-groups" -version = "0.2.0" +version = "0.3.0" edition = "2021" rust-version = "1.75" authors = ["Gestell Labs LLC"] @@ -28,3 +28,14 @@ proptest = { version = "1", optional = true } [dev-dependencies] rand = "0.8" rand_distr = "0.4" + +[[example]] +name = "basics" +required-features = ["rand"] + +[[example]] +name = "random_sampling" +required-features = ["rand"] + +[[example]] +name = "bch_vs_exp" diff --git a/deny.toml b/deny.toml index 4c86b80..34485b3 100644 --- a/deny.toml +++ b/deny.toml @@ -1,11 +1,11 @@ [advisories] -vulnerability = "deny" -unmaintained = "warn" -yanked = "warn" -notice = "warn" +version = 2 +ignore = [ + "RUSTSEC-2024-0436", # paste: unmaintained but stable, transitive via nalgebra +] [licenses] -unlicensed = "deny" +version = 2 allow = [ "MIT", "Apache-2.0", diff --git a/examples/basics.rs b/examples/basics.rs new file mode 100644 index 0000000..6739615 --- /dev/null +++ b/examples/basics.rs @@ -0,0 +1,56 @@ +//! Basic Lie group operations in 30 lines. +//! +//! Run: `cargo run --example basics` + +use lie_groups::prelude::*; + +fn main() { + // Group elements + let g = SU2::rotation_x(0.5); + let h = SU2::rotation_y(0.3); + println!("g = {}", g); + println!("h = {}", h); + + // Composition via * operator + let product = &g * &h; + println!("g·h = {}", product); + + // Inverse + let g_inv = g.inverse(); + let should_be_identity = g_inv.compose(&g); + println!( + "g⁻¹·g = {} (distance from I: {:.2e})", + should_be_identity, + should_be_identity.distance_to_identity() + ); + + // Exponential map: algebra → group + let x = Su2Algebra::new([0.1, 0.2, 0.3]); + let exp_x = SU2::exp(&x); + println!("exp({:?}) = {}", x, exp_x); + + // Logarithm: group → algebra + let log_g = g.log().unwrap(); + println!("log(g) = {:?}", log_g); + + // Compose a path of rotations + let path: Vec = (0..10).map(|i| SU2::rotation_z(0.1 * i as f64)).collect(); + let holonomy: SU2 = path.iter().product(); + println!("holonomy = {}", holonomy); + + // BCH formula: log(exp(X)·exp(Y)) ≈ X + Y + ½[X,Y] + ... + let y = Su2Algebra::new([0.05, -0.1, 0.0]); + let bch = bch_checked::(&x, &y, 5).unwrap(); + println!("BCH(X,Y) = {:?}", bch); + + // Works generically over any Lie group + fn parallel_transport std::iter::Product<&'a G>>(path: &[G]) -> G { + path.iter().product() + } + let u1_path = vec![ + U1::from_angle(0.1), + U1::from_angle(0.2), + U1::from_angle(0.3), + ]; + println!("U(1) holonomy = {}", parallel_transport(&u1_path)); +} diff --git a/examples/bch_vs_exp.rs b/examples/bch_vs_exp.rs new file mode 100644 index 0000000..6da1f16 --- /dev/null +++ b/examples/bch_vs_exp.rs @@ -0,0 +1,50 @@ +//! Verify the Baker-Campbell-Hausdorff formula against direct exp·exp. +//! +//! The BCH formula gives: log(exp(X)·exp(Y)) = X + Y + ½[X,Y] + ... +//! This example compares the truncated series against the exact answer. +//! +//! Run: `cargo run --example bch_vs_exp` + +use lie_groups::prelude::*; + +fn main() { + println!("=== BCH formula vs direct computation ===\n"); + + // Small algebra elements (BCH converges well) + let x = Su2Algebra::new([0.1, 0.0, 0.0]); + let y = Su2Algebra::new([0.0, 0.1, 0.0]); + + // Direct: log(exp(X) · exp(Y)) + let direct = SU2::exp(&x) * SU2::exp(&y); + let z_exact = direct.log().unwrap(); + + // BCH approximation at various orders + for order in [2, 3, 4, 5] { + let z_bch = lie_groups::bch_safe::(&x, &y, order).unwrap(); + let error = z_exact.add(&z_bch.scale(-1.0)).norm(); + println!("Order {}: error = {:.2e}", order, error); + } + + // bch_checked with highest order + let z_best = lie_groups::bch_safe::(&x, &y, 5).unwrap(); + let error_best = z_exact.add(&z_best.scale(-1.0)).norm(); + println!("Best(5): error = {:.2e}", error_best); + + // Show the bracket structure + println!("\n=== Lie bracket structure constants ===\n"); + for i in 0..3 { + for j in (i + 1)..3 { + let ei = Su2Algebra::basis_element(i); + let ej = Su2Algebra::basis_element(j); + let bracket = ei.bracket(&ej); + println!("[e{}, e{}] = {:?}", i, j, bracket.to_components()); + } + } + + // SU(3) bracket + println!("\n=== SU(3) brackets ===\n"); + let t1 = Su3Algebra::basis_element(0); + let t2 = Su3Algebra::basis_element(1); + let bracket_12 = t1.bracket(&t2); + println!("[T1, T2] = {:?}", bracket_12.to_components()); +} diff --git a/examples/random_sampling.rs b/examples/random_sampling.rs new file mode 100644 index 0000000..a0bc563 --- /dev/null +++ b/examples/random_sampling.rs @@ -0,0 +1,55 @@ +//! Sample random group elements and verify statistical properties. +//! +//! Run: `cargo run --example random_sampling` + +use lie_groups::prelude::*; +use rand::SeedableRng; + +fn main() { + let mut rng = rand::rngs::StdRng::seed_from_u64(42); + + // Sample from each compact group + println!("=== Random Haar sampling ===\n"); + + // SU(2) + let g2 = SU2::random_haar(&mut rng); + println!("SU(2): {}", g2); + println!(" unitary: {}", g2.verify_unitarity(1e-10)); + + // SO(3) + let r = SO3::random_haar(&mut rng); + println!("SO(3): {}", r); + println!(" orthogonal: {}", r.verify_orthogonality(1e-10)); + + // SU(3) — the QCD gauge group + let g3 = SU3::random_haar(&mut rng); + println!("SU(3): {}", g3); + println!(" unitary: {}", g3.verify_unitarity(1e-10)); + + // SU(4) — Pati-Salam model + let g4 = SU4::random_haar(&mut rng); + println!("SU(4): {}", g4); + println!(" unitary: {}", g4.verify_unitarity(1e-10)); + + // Verify Haar measure: average trace should be 0 for SU(N), N≥2 + println!("\n=== Haar measure check: → 0 ===\n"); + let n_samples = 10_000; + + let avg_trace_su2: f64 = (0..n_samples) + .map(|_| SU2::random_haar(&mut rng).trace().re) + .sum::() + / n_samples as f64; + println!("SU(2): = {:.6} (expect ~0)", avg_trace_su2); + + let avg_trace_su3: f64 = (0..n_samples) + .map(|_| SU3::random_haar(&mut rng).trace().re) + .sum::() + / n_samples as f64; + println!("SU(3): = {:.6} (expect ~0)", avg_trace_su3); + + let avg_trace_sun4: f64 = (0..n_samples) + .map(|_| SUN::<4>::random_haar(&mut rng).trace().re) + .sum::() + / n_samples as f64; + println!("SU(4): = {:.6} (expect ~0)", avg_trace_sun4); +} diff --git a/src/bch.rs b/src/bch.rs index fd9c825..eb78d74 100644 --- a/src/bch.rs +++ b/src/bch.rs @@ -1237,8 +1237,10 @@ mod tests { /// BCH vs exp·log on SUN<3> with non-axis-aligned inputs. #[test] fn test_bch_vs_exp_log_sun3() { - let x = SunAlgebra::<3>::from_components(&[0.05, -0.03, 0.02, 0.04, -0.01, 0.03, -0.02, 0.01]); - let y = SunAlgebra::<3>::from_components(&[-0.02, 0.04, -0.03, 0.01, 0.05, -0.04, 0.02, -0.03]); + let x = + SunAlgebra::<3>::from_components(&[0.05, -0.03, 0.02, 0.04, -0.01, 0.03, -0.02, 0.01]); + let y = + SunAlgebra::<3>::from_components(&[-0.02, 0.04, -0.03, 0.01, 0.05, -0.04, 0.02, -0.03]); let direct = SUN::<3>::exp(&x).compose(&SUN::<3>::exp(&y)); let bch_z = bch_fifth_order(&x, &y); @@ -1281,8 +1283,10 @@ mod tests { #[test] fn test_bch_su3_vs_sun3_consistency() { // Build algebra elements via SU3, convert to SUN<3> via matrix roundtrip - let x_su3 = Su3Algebra::from_components(&[0.05, -0.03, 0.02, 0.04, -0.01, 0.03, -0.02, 0.01]); - let y_su3 = Su3Algebra::from_components(&[-0.02, 0.04, -0.03, 0.01, 0.05, -0.04, 0.02, -0.03]); + let x_su3 = + Su3Algebra::from_components(&[0.05, -0.03, 0.02, 0.04, -0.01, 0.03, -0.02, 0.01]); + let y_su3 = + Su3Algebra::from_components(&[-0.02, 0.04, -0.03, 0.01, 0.05, -0.04, 0.02, -0.03]); let x_sun = SunAlgebra::<3>::from_matrix(&x_su3.to_matrix()); let y_sun = SunAlgebra::<3>::from_matrix(&y_su3.to_matrix()); @@ -1301,7 +1305,10 @@ mod tests { assert!( (m_su3[(r, c)] - m_sun[(r, c)]).norm() < 1e-12, "BCH matrix disagrees at ({},{}): SU3={}, SUN<3>={}", - r, c, m_su3[(r, c)], m_sun[(r, c)] + r, + c, + m_su3[(r, c)], + m_sun[(r, c)] ); } } diff --git a/src/conversions.rs b/src/conversions.rs new file mode 100644 index 0000000..2808a22 --- /dev/null +++ b/src/conversions.rs @@ -0,0 +1,460 @@ +//! Type-safe conversions between equivalent Lie group representations. +//! +//! This module provides [`From`] implementations connecting specialized +//! implementations to their generic counterparts: +//! +//! - [`SU2`] ↔ [`SUN<2>`]: Same 2×2 unitary matrices, different types +//! - [`Su2Algebra`] ↔ [`SunAlgebra<2>`]: Same Pauli basis, same ordering +//! - [`SU3`] ↔ [`SUN<3>`]: Same 3×3 unitary matrices, different types +//! - [`Su3Algebra`] ↔ [`SunAlgebra<3>`]: Gell-Mann basis with reordering +//! - [`UnitQuaternion`] ↔ [`SU2`]: Quaternion ↔ matrix isomorphism +//! +//! # Mathematical Guarantee +//! +//! These conversions are **structure-preserving** (Lie group homomorphisms): +//! +//! ```text +//! from(g · h) = from(g) · from(h) // preserves composition +//! from(exp(X)) = exp(from(X)) // commutes with exp +//! from(g⁻¹) = from(g)⁻¹ // preserves inverse +//! ``` +//! +//! These properties are verified by roundtrip tests in this module. +//! +//! # Basis Ordering: SU(3) ↔ SU(N=3) +//! +//! SU(3) uses the standard Gell-Mann ordering (λ₁..λ₈), while SU(N=3) groups +//! generators by type (symmetric, antisymmetric, diagonal). The conversion +//! applies the permutation: +//! +//! ```text +//! SU3 index: [0, 1, 2, 3, 4, 5, 6, 7] (Gell-Mann: λ₁..λ₈) +//! SUN<3> index: [0, 3, 6, 1, 4, 2, 5, 7] (grouped by type) +//! ``` + +use crate::quaternion::UnitQuaternion; +use crate::su2::{Su2Algebra, SU2}; +use crate::su3::{Su3Algebra, SU3}; +use crate::sun::{SunAlgebra, SUN}; + +// ============================================================================ +// SU2 ↔ SUN<2> (matrix representations share identical basis ordering) +// ============================================================================ + +impl From for SUN<2> { + fn from(g: SU2) -> Self { + Self { + matrix: g.matrix.clone(), + } + } +} + +impl From> for SU2 { + fn from(g: SUN<2>) -> Self { + Self { + matrix: g.matrix.clone(), + } + } +} + +impl From for SunAlgebra<2> { + fn from(x: Su2Algebra) -> Self { + Self::new(x.components().to_vec()) + } +} + +impl From> for Su2Algebra { + fn from(x: SunAlgebra<2>) -> Self { + let c = x.coefficients(); + Self::new([c[0], c[1], c[2]]) + } +} + +// ============================================================================ +// SU3 ↔ SUN<3> (matrix identical, algebra basis needs permutation) +// ============================================================================ + +/// SU3 Gell-Mann index → SUN<3> generalized Gell-Mann index. +/// +/// SU3 uses interleaved Gell-Mann ordering (λ₁..λ₈). +/// SUN<3> groups by type: symmetric, antisymmetric, diagonal. +const SU3_TO_SUN3: [usize; 8] = [0, 3, 6, 1, 4, 2, 5, 7]; + +/// SUN<3> generalized Gell-Mann index → SU3 Gell-Mann index (inverse permutation). +const SUN3_TO_SU3: [usize; 8] = [0, 3, 5, 1, 4, 6, 2, 7]; + +impl From for SUN<3> { + fn from(g: SU3) -> Self { + Self { + matrix: g.matrix().clone(), + } + } +} + +impl From> for SU3 { + fn from(g: SUN<3>) -> Self { + Self { + matrix: g.matrix().clone(), + } + } +} + +impl From for SunAlgebra<3> { + fn from(x: Su3Algebra) -> Self { + let c = x.components(); + let mut sun_coeffs = vec![0.0; 8]; + for i in 0..8 { + sun_coeffs[SU3_TO_SUN3[i]] = c[i]; + } + Self::new(sun_coeffs) + } +} + +impl From> for Su3Algebra { + fn from(x: SunAlgebra<3>) -> Self { + let c = x.coefficients(); + let mut su3_coeffs = [0.0; 8]; + for i in 0..8 { + su3_coeffs[SUN3_TO_SU3[i]] = c[i]; + } + Self::new(su3_coeffs) + } +} + +// ============================================================================ +// UnitQuaternion ↔ SU2 (the SU(2) ≅ S³ isomorphism) +// ============================================================================ + +impl From for SU2 { + /// Convert quaternion q = w + xi + yj + zk to the SU(2) matrix: + /// + /// ```text + /// U = [[ w + ix, -y + iz ], + /// [ y + iz, w - ix ]] + /// ``` + /// + /// **Convention note:** This is an anti-homomorphism with respect to + /// multiplication: `U(q₁·q₂) = U(q₂)·U(q₁)`. This is standard in + /// physics — quaternion left-action on vectors corresponds to matrix + /// right-multiplication. + fn from(q: UnitQuaternion) -> Self { + let m = q.to_matrix(); + Self::from_matrix(m) + } +} + +impl From for UnitQuaternion { + /// Extract quaternion from SU(2) matrix: + /// + /// ```text + /// w = Re(U₀₀), x = Im(U₀₀), y = Re(U₁₀) with sign flip, z = Im(U₁₀) + /// ``` + fn from(g: SU2) -> Self { + let m = g.to_matrix(); + UnitQuaternion::from_matrix(m) + } +} + +// ============================================================================ +// Tests: Structure-preserving roundtrip verification +// ============================================================================ + +#[cfg(test)] +mod tests { + use super::*; + use crate::traits::{LieAlgebra, LieGroup}; + + // ==================================================================== + // SU2 ↔ SUN<2> roundtrips + // ==================================================================== + + #[test] + fn test_su2_sun2_group_roundtrip() { + let g = SU2::rotation_x(0.7); + let sun: SUN<2> = g.clone().into(); + let back: SU2 = sun.into(); + assert!( + g.distance(&back) < 1e-14, + "SU2 → SUN<2> → SU2 roundtrip failed" + ); + } + + #[test] + fn test_su2_sun2_algebra_roundtrip() { + let x = Su2Algebra::new([0.3, -0.7, 1.2]); + let sun: SunAlgebra<2> = x.into(); + let back: Su2Algebra = sun.into(); + assert_eq!( + x, back, + "Su2Algebra → SunAlgebra<2> → Su2Algebra roundtrip failed" + ); + } + + #[test] + fn test_su2_sun2_compose_homomorphism() { + let g = SU2::rotation_x(0.3); + let h = SU2::rotation_y(0.7); + let product = g.compose(&h); + + let g_sun: SUN<2> = g.into(); + let h_sun: SUN<2> = h.into(); + let product_sun: SUN<2> = product.into(); + + let composed_sun = g_sun.compose(&h_sun); + assert!( + product_sun.distance(&composed_sun) < 1e-13, + "from(g·h) ≠ from(g)·from(h): {}", + product_sun.distance(&composed_sun) + ); + } + + #[test] + fn test_su2_sun2_exp_commutes() { + let x = Su2Algebra::new([0.3, -0.2, 0.4]); + let g = SU2::exp(&x); + + let x_sun: SunAlgebra<2> = x.into(); + let g_sun: SUN<2> = g.into(); + let exp_sun = SUN::<2>::exp(&x_sun); + + assert!( + g_sun.distance(&exp_sun) < 1e-13, + "from(exp(X)) ≠ exp(from(X)): {}", + g_sun.distance(&exp_sun) + ); + } + + #[test] + fn test_su2_sun2_inverse_commutes() { + let g = SU2::rotation_z(1.5); + let g_inv = g.inverse(); + + let g_sun: SUN<2> = g.into(); + let g_inv_sun: SUN<2> = g_inv.into(); + + assert!( + g_sun.inverse().distance(&g_inv_sun) < 1e-14, + "from(g⁻¹) ≠ from(g)⁻¹" + ); + } + + #[test] + fn test_su2_sun2_bracket_commutes() { + let x = Su2Algebra::new([1.0, 0.0, 0.0]); + let y = Su2Algebra::new([0.0, 1.0, 0.0]); + let bracket = x.bracket(&y); + + let x_sun: SunAlgebra<2> = x.into(); + let y_sun: SunAlgebra<2> = y.into(); + let bracket_sun = x_sun.bracket(&y_sun); + + let bracket_converted: SunAlgebra<2> = bracket.into(); + let diff: f64 = bracket_converted + .coefficients() + .iter() + .zip(bracket_sun.coefficients().iter()) + .map(|(a, b)| (a - b).abs()) + .sum(); + assert!( + diff < 1e-13, + "from([X,Y]) ≠ [from(X),from(Y)]: diff={}", + diff + ); + } + + // ==================================================================== + // SU3 ↔ SUN<3> roundtrips + // ==================================================================== + + #[test] + fn test_su3_sun3_group_roundtrip() { + let x = Su3Algebra::new([0.1, 0.2, 0.3, 0.0, 0.0, 0.0, 0.0, 0.0]); + let g = SU3::exp(&x); + let sun: SUN<3> = g.clone().into(); + let back: SU3 = sun.into(); + assert!( + g.distance(&back) < 1e-13, + "SU3 → SUN<3> → SU3 roundtrip failed: {}", + g.distance(&back) + ); + } + + #[test] + fn test_su3_sun3_algebra_roundtrip() { + let x = Su3Algebra::new([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]); + let sun: SunAlgebra<3> = x.into(); + let back: Su3Algebra = sun.into(); + assert_eq!( + x, back, + "Su3Algebra → SunAlgebra<3> → Su3Algebra roundtrip failed" + ); + } + + #[test] + fn test_su3_sun3_compose_homomorphism() { + let x = Su3Algebra::new([0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]); + let y = Su3Algebra::new([0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]); + let g = SU3::exp(&x); + let h = SU3::exp(&y); + let product = g.compose(&h); + + let g_sun: SUN<3> = g.into(); + let h_sun: SUN<3> = h.into(); + let product_sun: SUN<3> = product.into(); + + let composed_sun = g_sun.compose(&h_sun); + assert!( + product_sun.distance(&composed_sun) < 1e-12, + "from(g·h) ≠ from(g)·from(h): {}", + product_sun.distance(&composed_sun) + ); + } + + #[test] + fn test_su3_sun3_exp_commutes() { + let x = Su3Algebra::new([0.1, 0.2, -0.1, 0.05, 0.0, 0.0, 0.0, 0.0]); + let g = SU3::exp(&x); + + let x_sun: SunAlgebra<3> = x.into(); + let g_sun: SUN<3> = g.into(); + let exp_sun = SUN::<3>::exp(&x_sun); + + assert!( + g_sun.distance(&exp_sun) < 1e-12, + "from(exp(X)) ≠ exp(from(X)): {}", + g_sun.distance(&exp_sun) + ); + } + + #[test] + fn test_su3_sun3_bracket_commutes() { + // Bracket in SU3 coordinates, convert, compare to bracket in SUN<3> coordinates + let x = Su3Algebra::new([1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]); + let y = Su3Algebra::new([0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]); + + // Compute bracket in SU3, then convert to SUN<3> + let bracket_su3 = x.bracket(&y); + let bracket_as_sun: SunAlgebra<3> = bracket_su3.into(); + + // Convert to SUN<3>, then compute bracket + let x_sun: SunAlgebra<3> = x.into(); + let y_sun: SunAlgebra<3> = y.into(); + let bracket_sun = x_sun.bracket(&y_sun); + + // Compare via matrix representation (avoids coefficient ordering issues) + let m1 = bracket_as_sun.to_matrix(); + let m2 = bracket_sun.to_matrix(); + for r in 0..3 { + for c in 0..3 { + assert!( + (m1[(r, c)] - m2[(r, c)]).norm() < 1e-12, + "from([X,Y]) ≠ [from(X),from(Y)] at ({},{}): {} vs {}", + r, + c, + m1[(r, c)], + m2[(r, c)] + ); + } + } + } + + // ==================================================================== + // Permutation self-consistency + // ==================================================================== + + #[test] + fn test_su3_sun3_permutation_is_inverse() { + for i in 0..8 { + assert_eq!( + super::SUN3_TO_SU3[super::SU3_TO_SUN3[i]], + i, + "SUN3_TO_SU3[SU3_TO_SUN3[{}]] ≠ {}", + i, + i + ); + assert_eq!( + super::SU3_TO_SUN3[super::SUN3_TO_SU3[i]], + i, + "SU3_TO_SUN3[SUN3_TO_SU3[{}]] ≠ {}", + i, + i + ); + } + } + + // ==================================================================== + // UnitQuaternion ↔ SU2 roundtrips + // ==================================================================== + + #[test] + fn test_quaternion_su2_group_roundtrip() { + let q = UnitQuaternion::exp([0.3, -0.5, 0.2]); + let g: SU2 = q.into(); + let back: UnitQuaternion = g.into(); + + // Compare quaternion components (up to global sign — q and -q are the same SU2 element) + let same = (q.w() - back.w()).abs() < 1e-14 + && (q.x() - back.x()).abs() < 1e-14 + && (q.y() - back.y()).abs() < 1e-14 + && (q.z() - back.z()).abs() < 1e-14; + let negated = (q.w() + back.w()).abs() < 1e-14 + && (q.x() + back.x()).abs() < 1e-14 + && (q.y() + back.y()).abs() < 1e-14 + && (q.z() + back.z()).abs() < 1e-14; + assert!( + same || negated, + "Quaternion roundtrip failed: q={:?}, back={:?}", + q, + back + ); + } + + #[test] + fn test_quaternion_su2_compose_anti_homomorphism() { + // The standard SU(2) parameterization U = [[α, -β*], [β, α*]] + // gives an anti-homomorphism: U(q₁·q₂) = U(q₂)·U(q₁). + // + // This is the standard physics convention — quaternion left-multiplication + // corresponds to matrix right-multiplication. Both represent the same + // rotation; the ordering reflects the column-vector vs row-vector convention. + let q1 = UnitQuaternion::exp([0.1, 0.2, 0.3]); + let q2 = UnitQuaternion::exp([0.4, -0.1, 0.2]); + let q_product = q1 * q2; + + let g1: SU2 = q1.into(); + let g2: SU2 = q2.into(); + let g_product: SU2 = q_product.into(); + + // Anti-homomorphism: from(q1*q2) = from(q2).compose(&from(q1)) + let g_anti = g2.compose(&g1); + assert!( + g_product.distance(&g_anti) < 1e-7, + "from(q1·q2) ≠ from(q2)·from(q1): {}", + g_product.distance(&g_anti) + ); + } + + #[test] + fn test_quaternion_su2_inverse_commutes() { + let q = UnitQuaternion::exp([0.5, -0.3, 0.8]); + let q_inv = q.inverse(); + + let g: SU2 = q.into(); + let g_inv: SU2 = q_inv.into(); + + assert!( + g.inverse().distance(&g_inv) < 1e-14, + "from(q⁻¹) ≠ from(q)⁻¹" + ); + } + + #[test] + fn test_quaternion_su2_identity() { + let q = UnitQuaternion::identity(); + let g: SU2 = q.into(); + assert!( + g.is_near_identity(1e-14), + "Quaternion identity should map to SU2 identity" + ); + } +} diff --git a/src/lib.rs b/src/lib.rs index a2b94a3..087e6bc 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -39,7 +39,9 @@ //! - **Numerical stability**: Conditioned logarithms, scaling-and-squaring exp pub mod bch; +pub mod conversions; pub mod error; +pub mod prelude; pub mod quaternion; pub mod representation; pub mod root_systems; diff --git a/src/prelude.rs b/src/prelude.rs new file mode 100644 index 0000000..aed5b5c --- /dev/null +++ b/src/prelude.rs @@ -0,0 +1,20 @@ +//! Convenience re-exports for common usage. +//! +//! ``` +//! use lie_groups::prelude::*; +//! +//! let g = SU2::rotation_x(0.5); +//! let h = SU2::rotation_y(0.3); +//! let product = &g * &h; +//! ``` + +pub use crate::bch::{bch_checked, BchError, BchMethod}; +pub use crate::error::{LogError, LogResult}; +pub use crate::quaternion::UnitQuaternion; +pub use crate::rplus::{RPlus, RPlusAlgebra}; +pub use crate::so3::{So3Algebra, SO3}; +pub use crate::su2::{Su2Algebra, SU2}; +pub use crate::su3::{Su3Algebra, SU3}; +pub use crate::sun::{SunAlgebra, SU4, SU5, SUN}; +pub use crate::traits::{LieAlgebra, LieGroup}; +pub use crate::u1::{U1Algebra, U1}; diff --git a/src/rplus.rs b/src/rplus.rs index 6cc530a..1f27131 100644 --- a/src/rplus.rs +++ b/src/rplus.rs @@ -377,12 +377,31 @@ impl Mul<&RPlus> for RPlus { } } +impl Mul for RPlus { + type Output = RPlus; + fn mul(self, rhs: RPlus) -> RPlus { + &self * &rhs + } +} + impl MulAssign<&RPlus> for RPlus { fn mul_assign(&mut self, rhs: &RPlus) { *self = self.compose(rhs); } } +impl std::iter::Product for RPlus { + fn product>(iter: I) -> Self { + iter.fold(Self::from_value(1.0), |acc, g| acc * g) + } +} + +impl<'a> std::iter::Product<&'a RPlus> for RPlus { + fn product>(iter: I) -> Self { + iter.fold(Self::from_value(1.0), |acc, g| acc * g) + } +} + impl LieGroup for RPlus { const MATRIX_DIM: usize = 1; diff --git a/src/so3.rs b/src/so3.rs index 97c6e44..da8eb6f 100644 --- a/src/so3.rs +++ b/src/so3.rs @@ -402,6 +402,50 @@ impl SO3 { (product - identity).norm() < tolerance } + /// Random SO(3) element uniformly distributed according to Haar measure. + /// + /// Uses the SU(2) → SO(3) double cover: sample a random quaternion on S³, + /// convert to rotation matrix via the adjoint representation. + /// + /// # Examples + /// + /// ```rust + /// use lie_groups::SO3; + /// use rand::SeedableRng; + /// + /// let mut rng = rand::rngs::StdRng::seed_from_u64(42); + /// let r = SO3::random_haar(&mut rng); + /// assert!(r.verify_orthogonality(1e-10)); + /// ``` + #[cfg(feature = "rand")] + #[must_use] + pub fn random_haar(rng: &mut R) -> Self { + use crate::UnitQuaternion; + use rand_distr::{Distribution, StandardNormal}; + + const MIN_NORM: f64 = 1e-10; + loop { + let a: f64 = StandardNormal.sample(rng); + let b: f64 = StandardNormal.sample(rng); + let c: f64 = StandardNormal.sample(rng); + let d: f64 = StandardNormal.sample(rng); + + let norm = (a * a + b * b + c * c + d * d).sqrt(); + if norm < MIN_NORM { + continue; + } + + let q = UnitQuaternion::new(a / norm, b / norm, c / norm, d / norm); + let rot = q.to_rotation_matrix(); + return Self { + matrix: Matrix3::new( + rot[0][0], rot[0][1], rot[0][2], rot[1][0], rot[1][1], rot[1][2], rot[2][0], + rot[2][1], rot[2][2], + ), + }; + } + } + /// Matrix inverse (equals transpose for orthogonal matrices) #[must_use] pub fn inverse(&self) -> Self { @@ -559,8 +603,27 @@ impl fmt::Display for So3Algebra { impl fmt::Display for SO3 { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - let dist = self.distance_to_identity(); - write!(f, "SO(3)(θ={:.4})", dist) + let angle = self.distance_to_identity(); + if angle.abs() < 1e-12 { + write!(f, "SO(3)(I)") + } else if let Ok(log) = self.log() { + let c = log.components(); + let n = (c[0] * c[0] + c[1] * c[1] + c[2] * c[2]).sqrt(); + if n > 1e-12 { + write!( + f, + "SO(3)(θ={:.4}, n̂=[{:.3}, {:.3}, {:.3}])", + angle, + c[0] / n, + c[1] / n, + c[2] / n + ) + } else { + write!(f, "SO(3)(θ={:.4})", angle) + } + } else { + write!(f, "SO(3)(θ={:.4})", angle) + } } } @@ -581,12 +644,31 @@ impl Mul<&SO3> for SO3 { } } +impl Mul for SO3 { + type Output = SO3; + fn mul(self, rhs: SO3) -> SO3 { + &self * &rhs + } +} + impl MulAssign<&SO3> for SO3 { fn mul_assign(&mut self, rhs: &SO3) { self.matrix *= rhs.matrix; } } +impl std::iter::Product for SO3 { + fn product>(iter: I) -> Self { + iter.fold(Self::identity(), |acc, g| acc * g) + } +} + +impl<'a> std::iter::Product<&'a SO3> for SO3 { + fn product>(iter: I) -> Self { + iter.fold(Self::identity(), |acc, g| &acc * g) + } +} + impl LieGroup for SO3 { const MATRIX_DIM: usize = 3; diff --git a/src/su2.rs b/src/su2.rs index 962c7b2..b18ebd4 100644 --- a/src/su2.rs +++ b/src/su2.rs @@ -887,12 +887,32 @@ impl Mul<&SU2> for SU2 { } } +impl Mul for SU2 { + type Output = SU2; + + fn mul(self, rhs: SU2) -> SU2 { + &self * &rhs + } +} + impl MulAssign<&SU2> for SU2 { fn mul_assign(&mut self, rhs: &SU2) { self.matrix = self.matrix.dot(&rhs.matrix); } } +impl std::iter::Product for SU2 { + fn product>(iter: I) -> Self { + iter.fold(Self::identity(), |acc, g| acc * g) + } +} + +impl<'a> std::iter::Product<&'a SU2> for SU2 { + fn product>(iter: I) -> Self { + iter.fold(Self::identity(), |acc, g| &acc * g) + } +} + impl approx::AbsDiffEq for Su2Algebra { type Epsilon = f64; @@ -938,9 +958,16 @@ impl fmt::Display for Su2Algebra { impl fmt::Display for SU2 { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - // Show rotation angle and distance to identity - let dist = self.distance_to_identity(); - write!(f, "SU(2)(θ={:.4})", dist) + let (angle, axis) = self.angle_and_axis(); + if angle.abs() < 1e-12 { + write!(f, "SU(2)(I)") + } else { + write!( + f, + "SU(2)(θ={:.4}, n̂=[{:.3}, {:.3}, {:.3}])", + angle, axis[0], axis[1], axis[2] + ) + } } } diff --git a/src/su3.rs b/src/su3.rs index 97101c4..d734b03 100644 --- a/src/su3.rs +++ b/src/su3.rs @@ -668,6 +668,27 @@ impl SU3 { norm < tolerance } + /// Random SU(3) element uniformly distributed according to Haar measure. + /// + /// Delegates to the generic `SUN<3>::random_haar` implementation using + /// the Mezzadri (2007) algorithm, then converts via `From>`. + /// + /// # Examples + /// + /// ```rust + /// use lie_groups::SU3; + /// use rand::SeedableRng; + /// + /// let mut rng = rand::rngs::StdRng::seed_from_u64(42); + /// let g = SU3::random_haar(&mut rng); + /// assert!(g.verify_unitarity(1e-10)); + /// ``` + #[cfg(feature = "rand")] + #[must_use] + pub fn random_haar(rng: &mut R) -> Self { + crate::sun::SUN::<3>::random_haar(rng).into() + } + /// Matrix inverse (equals conjugate transpose for unitary matrices) #[must_use] pub fn inverse(&self) -> Self { @@ -979,7 +1000,12 @@ impl fmt::Display for Su3Algebra { impl fmt::Display for SU3 { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { let dist = self.distance_to_identity(); - write!(f, "SU(3)(d={:.4})", dist) + let tr = self.trace(); + if dist < 1e-12 { + write!(f, "SU(3)(I)") + } else { + write!(f, "SU(3)(d={:.4}, tr={:.3}{:+.3}i)", dist, tr.re, tr.im) + } } } @@ -1000,12 +1026,31 @@ impl Mul<&SU3> for SU3 { } } +impl Mul for SU3 { + type Output = SU3; + fn mul(self, rhs: SU3) -> SU3 { + &self * &rhs + } +} + impl MulAssign<&SU3> for SU3 { fn mul_assign(&mut self, rhs: &SU3) { self.matrix = self.matrix.dot(&rhs.matrix); } } +impl std::iter::Product for SU3 { + fn product>(iter: I) -> Self { + iter.fold(Self::identity(), |acc, g| acc * g) + } +} + +impl<'a> std::iter::Product<&'a SU3> for SU3 { + fn product>(iter: I) -> Self { + iter.fold(Self::identity(), |acc, g| &acc * g) + } +} + impl LieGroup for SU3 { const MATRIX_DIM: usize = 3; diff --git a/src/sun.rs b/src/sun.rs index 1cf1de3..3f2966e 100644 --- a/src/sun.rs +++ b/src/sun.rs @@ -556,6 +556,146 @@ impl SUN { norm_sq.sqrt() < tolerance } + /// Random SU(N) element uniformly distributed according to Haar measure. + /// + /// Uses the standard algorithm: generate an N×N matrix of i.i.d. complex + /// Gaussians, compute the QR decomposition, and adjust phases so the + /// result has determinant 1. + /// + /// # Algorithm (Mezzadri 2007) + /// + /// 1. Sample Z: N×N matrix with Z_{ij} ~ N(0,1) + i·N(0,1) + /// 2. QR decompose Z = Q·R + /// 3. Adjust Q: multiply columns by sign(R_{ii}/|R_{ii}|) to make + /// the decomposition unique (positive diagonal R) + /// 4. Fix determinant: multiply first column by det(Q)* to get det=1 + /// + /// # References + /// + /// - Mezzadri, F.: "How to generate random matrices from the classical + /// compact groups" (Notices AMS, 2007) + /// + /// # Examples + /// + /// ```rust + /// use lie_groups::SUN; + /// use rand::SeedableRng; + /// + /// let mut rng = rand::rngs::StdRng::seed_from_u64(42); + /// let g = SUN::<3>::random_haar(&mut rng); + /// assert!(g.verify_unitarity(1e-10)); + /// ``` + #[cfg(feature = "rand")] + #[must_use] + pub fn random_haar(rng: &mut R) -> Self { + use ndarray::Array2; + use num_complex::Complex64; + use rand_distr::{Distribution, StandardNormal}; + + // Step 1: Generate random complex Gaussian matrix + let mut z = Array2::::zeros((N, N)); + for i in 0..N { + for j in 0..N { + let re: f64 = StandardNormal.sample(rng); + let im: f64 = StandardNormal.sample(rng); + z[[i, j]] = Complex64::new(re, im); + } + } + + // Step 2: Modified Gram-Schmidt QR decomposition + let mut q = z; + let mut diag_r = vec![Complex64::new(0.0, 0.0); N]; + for j in 0..N { + // Orthogonalize column j against previous columns + for k in 0..j { + let mut dot = Complex64::new(0.0, 0.0); + for i in 0..N { + dot += q[[i, k]].conj() * q[[i, j]]; + } + // Copy column k to avoid borrow conflict + let col_k: Vec = (0..N).map(|i| q[[i, k]]).collect(); + for i in 0..N { + q[[i, j]] -= dot * col_k[i]; + } + } + // Normalize + let mut norm_sq = 0.0; + for i in 0..N { + norm_sq += q[[i, j]].norm_sqr(); + } + let norm = norm_sq.sqrt(); + diag_r[j] = Complex64::new(norm, 0.0); + if norm > 1e-15 { + for i in 0..N { + q[[i, j]] /= norm; + } + } + } + + // Step 3: Adjust phases (Mezzadri correction) + // Multiply each column by the phase of the diagonal R element + // This makes the QR decomposition unique + for j in 0..N { + let r_jj = diag_r[j]; + if r_jj.norm() > 1e-15 { + let phase = r_jj / r_jj.norm(); + for i in 0..N { + q[[i, j]] *= phase.conj(); + } + } + } + + // Step 4: Fix determinant to 1 + // Compute det(Q) via product of diagonal after triangularization + // For small N, we use the direct formula or cofactor expansion + let det = Self::compute_det(&q); + if det.norm() > 1e-15 { + let phase = det / det.norm(); + // Multiply first column by phase* to set det = 1 + for i in 0..N { + q[[i, 0]] *= phase.conj(); + } + } + + Self { matrix: q } + } + + /// Compute determinant of an N×N complex matrix. + /// + /// Uses LU-style elimination for general N. + fn compute_det(m: &Array2) -> Complex64 { + let n = m.nrows(); + if n == 1 { + return m[[0, 0]]; + } + if n == 2 { + return m[[0, 0]] * m[[1, 1]] - m[[0, 1]] * m[[1, 0]]; + } + if n == 3 { + return m[[0, 0]] * (m[[1, 1]] * m[[2, 2]] - m[[1, 2]] * m[[2, 1]]) + - m[[0, 1]] * (m[[1, 0]] * m[[2, 2]] - m[[1, 2]] * m[[2, 0]]) + + m[[0, 2]] * (m[[1, 0]] * m[[2, 1]] - m[[1, 1]] * m[[2, 0]]); + } + // General case: cofactor expansion along first row + let mut det = Complex64::new(0.0, 0.0); + for j in 0..n { + let mut minor = Array2::::zeros((n - 1, n - 1)); + for ii in 1..n { + let mut col = 0; + for jj in 0..n { + if jj == j { + continue; + } + minor[[ii - 1, col]] = m[[ii, jj]]; + col += 1; + } + } + let sign = if j % 2 == 0 { 1.0 } else { -1.0 }; + det += Complex64::new(sign, 0.0) * m[[0, j]] * Self::compute_det(&minor); + } + det + } + /// Compute determinant /// /// For SU(N), the determinant should be exactly 1 by definition. @@ -774,7 +914,12 @@ impl fmt::Display for SunAlgebra { impl fmt::Display for SUN { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { let dist = self.distance_to_identity(); - write!(f, "SU({})(d={:.4})", N, dist) + let tr = self.trace(); + if dist < 1e-12 { + write!(f, "SU({})(I)", N) + } else { + write!(f, "SU({})(d={:.4}, tr={:.3}{:+.3}i)", N, dist, tr.re, tr.im) + } } } @@ -795,12 +940,31 @@ impl Mul<&SUN> for SUN { } } +impl Mul> for SUN { + type Output = SUN; + fn mul(self, rhs: SUN) -> SUN { + &self * &rhs + } +} + impl MulAssign<&SUN> for SUN { fn mul_assign(&mut self, rhs: &SUN) { self.matrix = self.matrix.dot(&rhs.matrix); } } +impl std::iter::Product for SUN { + fn product>(iter: I) -> Self { + iter.fold(Self::identity(), |acc, g| acc * g) + } +} + +impl<'a, const N: usize> std::iter::Product<&'a SUN> for SUN { + fn product>(iter: I) -> Self { + iter.fold(Self::identity(), |acc, g| &acc * g) + } +} + impl LieGroup for SUN { type Algebra = SunAlgebra; const MATRIX_DIM: usize = { @@ -1465,7 +1629,9 @@ mod tests { let x = Su3Algebra::from_components(&[0.3, -0.2, 0.1, 0.15, -0.1, 0.25, -0.15, 0.05]); let y = Su3Algebra::from_components(&[-0.1, 0.25, -0.15, 0.2, 0.05, -0.1, 0.3, -0.2]); - let g = SU3::exp(&Su3Algebra::from_components(&[0.4, -0.3, 0.2, 0.1, -0.2, 0.15, -0.1, 0.3])); + let g = SU3::exp(&Su3Algebra::from_components(&[ + 0.4, -0.3, 0.2, 0.1, -0.2, 0.15, -0.1, 0.3, + ])); let lhs = g.adjoint_action(&x.bracket(&y)); let rhs = g.adjoint_action(&x).bracket(&g.adjoint_action(&y)); @@ -1508,7 +1674,10 @@ mod tests { assert!( (m_su3[(r, c)] - m_sun[(r, c)]).norm() < 1e-12, "Bracket matrix disagrees at ({},{}): SU3={}, SUN<3>={}", - r, c, m_su3[(r, c)], m_sun[(r, c)] + r, + c, + m_su3[(r, c)], + m_sun[(r, c)] ); } } @@ -1520,7 +1689,8 @@ mod tests { assert!( (m_su3[(r, c)] - roundtrip[(r, c)]).norm() < 1e-12, "from_matrix roundtrip failed at ({},{})", - r, c + r, + c ); } } diff --git a/src/u1.rs b/src/u1.rs index c02a453..6174543 100644 --- a/src/u1.rs +++ b/src/u1.rs @@ -519,12 +519,31 @@ impl Mul<&U1> for U1 { } } +impl Mul for U1 { + type Output = U1; + fn mul(self, rhs: U1) -> U1 { + &self * &rhs + } +} + impl MulAssign<&U1> for U1 { fn mul_assign(&mut self, rhs: &U1) { *self = self.compose(rhs); } } +impl std::iter::Product for U1 { + fn product>(iter: I) -> Self { + iter.fold(Self::from_angle(0.0), |acc, g| acc * g) + } +} + +impl<'a> std::iter::Product<&'a U1> for U1 { + fn product>(iter: I) -> Self { + iter.fold(Self::from_angle(0.0), |acc, g| acc * g) + } +} + impl LieGroup for U1 { const MATRIX_DIM: usize = 1;