alice/groups/
so.rs

1use alloc::vec::Vec;
2use core::marker::PhantomData;
3use crate::core::ops::{Group, Magma, Monoid, Semigroup, TopologicalGroup};
4use crate::groups::matrix_group::{GroupError, MatrixGroup};
5use crate::topology::manifold::{Atlas, Dim, Manifold};
6use crate::core::scalar::FiniteF64;
7use crate::maps::exp_log::HasExpMap;
8
9const DET_TOLERANCE: f64 = 1e-10;
10const ORTHO_TOLERANCE: f64 = 1e-8;
11const SKEW_TOLERANCE: f64 = 1e-10;
12
13#[derive(Debug, Clone, PartialEq, Eq)]
14pub struct So<D: Dim> {
15    inner: MatrixGroup,
16    _dim: PhantomData<D>,
17}
18
19#[derive(Debug, Clone, PartialEq, Eq)]
20pub struct SoAlgebra<D: Dim> {
21    inner: MatrixGroup,
22    _dim: PhantomData<D>,
23}
24
25impl<D: Dim> So<D> {
26    pub fn new(data: Vec<f64>, dim: D) -> Result<Self, GroupError> {
27        let n = dim.value();
28        let inner = MatrixGroup::new(data, n)?;
29        if libm::fabs(inner.det() - 1.0) > DET_TOLERANCE {
30            return Err(GroupError::ConstraintViolated("det must equal 1"));
31        }
32        let t = inner.transpose();
33        let rrt = inner.mul(&t).expect("multiplication must succeed");
34        for i in 0..n {
35            for j in 0..n {
36                let expected = if i == j { 1.0 } else { 0.0 };
37                if libm::fabs(rrt.get(i, j) - expected) > ORTHO_TOLERANCE {
38                    return Err(GroupError::ConstraintViolated("R*R^T must equal identity"));
39                }
40            }
41        }
42        Ok(Self { inner, _dim: PhantomData })
43    }
44
45    pub fn identity(dim: D) -> Self {
46        Self { inner: MatrixGroup::identity(dim.value()), _dim: PhantomData }
47    }
48
49    pub fn n(&self) -> usize { self.inner.n() }
50    pub fn data(&self) -> Vec<f64> { self.inner.data() }
51
52    pub fn transpose(&self) -> Self {
53        Self { inner: self.inner.transpose(), _dim: PhantomData }
54    }
55}
56
57impl<D: Dim> SoAlgebra<D> {
58    pub fn new(data: Vec<f64>, dim: D) -> Result<Self, GroupError> {
59        let n = dim.value();
60        let inner = MatrixGroup::new(data, n)?;
61        for i in 0..n {
62            for j in 0..n {
63                if libm::fabs(inner.get(i, j) + inner.get(j, i)) > SKEW_TOLERANCE {
64                    return Err(GroupError::ConstraintViolated("matrix must be skew-symmetric"));
65                }
66            }
67        }
68        Ok(Self { inner, _dim: PhantomData })
69    }
70
71    pub fn zero(dim: D) -> Self {
72        Self {
73            inner: MatrixGroup::new(alloc::vec![0.0; dim.value() * dim.value()], dim.value()).unwrap(),
74            _dim: PhantomData,
75        }
76    }
77
78    pub fn data(&self) -> Vec<f64> { self.inner.data() }
79    pub fn n(&self) -> usize { self.inner.n() }
80
81    pub fn bracket(&self, other: &Self) -> Self {
82        let n = self.inner.n();
83        let mut result = alloc::vec![0.0f64; n * n];
84        for i in 0..n {
85            for j in 0..n {
86                let mut val = 0.0f64;
87                for k in 0..n {
88                    val += self.inner.get(i, k) * other.inner.get(k, j)
89                         - other.inner.get(i, k) * self.inner.get(k, j);
90                }
91                result[i * n + j] = val;
92            }
93        }
94        Self {
95            inner: MatrixGroup::new(result, n).unwrap(),
96            _dim: PhantomData,
97        }
98    }
99}
100
101impl<D: Dim> Magma for So<D> {
102    fn op(&self, other: &Self) -> Self {
103        Self {
104            inner: self.inner.mul(&other.inner).expect("SO multiplication must succeed"),
105            _dim: PhantomData,
106        }
107    }
108}
109
110impl<D: Dim> Semigroup for So<D> {}
111
112impl<D: Dim> Monoid for So<D> {
113    fn identity() -> Self {
114        panic!("SO identity requires dimension; use So::identity(dim)")
115    }
116}
117
118impl<D: Dim> Group for So<D> {
119    fn inverse(&self) -> Self {
120        Self { inner: self.inner.transpose(), _dim: PhantomData }
121    }
122}
123
124impl<D: Dim> TopologicalGroup for So<D> {}
125
126impl<D: Dim> Manifold for So<D> {
127    type Scalar = FiniteF64;
128    fn dim(&self) -> usize { let n = self.inner.n(); n * (n - 1) / 2 }
129    fn atlas(&self) -> &Atlas<FiniteF64> { unimplemented!("SO atlas not yet constructed") }
130}
131
132impl<D: Dim> HasExpMap for So<D> {
133    type Algebra = SoAlgebra<D>;
134    fn exp(x: &SoAlgebra<D>) -> Self {
135        Self { inner: x.inner.exp(), _dim: PhantomData }
136    }
137    fn log(&self) -> Option<SoAlgebra<D>> {
138        Some(SoAlgebra { inner: self.inner.log()?, _dim: PhantomData })
139    }
140}
alice/groups/so.rs

alice/groups/
so.rs
