qualia_core_db/solvers/calculus/
exterior.rs

1//! Bounded exterior algebra for three-dimensional calculus and geometry.
2
3use super::analysis::{AnalysisError, LinearMap, Vector};
4
5pub const EXTERIOR3_COMPONENTS: usize = 8;
6
7#[repr(C)]
8#[derive(Debug, Clone, Copy, PartialEq)]
9pub struct Exterior3 {
10    /// Coefficients indexed by basis bit mask: 0, e1, e2, e12, e3, ...
11    pub coefficients: [f64; EXTERIOR3_COMPONENTS],
12}
13
14impl Exterior3 {
15    pub const fn zero() -> Self {
16        Self {
17            coefficients: [0.0; EXTERIOR3_COMPONENTS],
18        }
19    }
20
21    pub const fn scalar(value: f64) -> Self {
22        let mut coefficients = [0.0; EXTERIOR3_COMPONENTS];
23        coefficients[0] = value;
24        Self { coefficients }
25    }
26
27    pub const fn basis(mask: usize, value: f64) -> Self {
28        let mut coefficients = [0.0; EXTERIOR3_COMPONENTS];
29        if mask < EXTERIOR3_COMPONENTS {
30            coefficients[mask] = value;
31        }
32        Self { coefficients }
33    }
34
35    pub fn validate(&self) -> Result<(), AnalysisError> {
36        if self.coefficients.iter().all(|value| value.is_finite()) {
37            Ok(())
38        } else {
39            Err(AnalysisError::NonFinite)
40        }
41    }
42
43    pub fn grade(self, grade: u32) -> Self {
44        let mut result = Self::zero();
45        for mask in 0..EXTERIOR3_COMPONENTS {
46            if mask.count_ones() == grade {
47                result.coefficients[mask] = self.coefficients[mask];
48            }
49        }
50        result
51    }
52
53    pub fn add(self, other: Self) -> Self {
54        let mut result = Self::zero();
55        for mask in 0..EXTERIOR3_COMPONENTS {
56            result.coefficients[mask] = self.coefficients[mask] + other.coefficients[mask];
57        }
58        result
59    }
60
61    pub fn scale(self, scalar: f64) -> Result<Self, AnalysisError> {
62        if !scalar.is_finite() {
63            return Err(AnalysisError::NonFinite);
64        }
65        let mut result = self;
66        for value in &mut result.coefficients {
67            *value *= scalar;
68        }
69        Ok(result)
70    }
71
72    pub fn wedge(self, other: Self) -> Result<Self, AnalysisError> {
73        self.validate()?;
74        other.validate()?;
75        let mut result = Self::zero();
76        for left in 0..EXTERIOR3_COMPONENTS {
77            for right in 0..EXTERIOR3_COMPONENTS {
78                if left & right != 0 {
79                    continue;
80                }
81                let (mask, sign) = wedge_basis(left, right);
82                result.coefficients[mask] +=
83                    sign * self.coefficients[left] * other.coefficients[right];
84            }
85        }
86        Ok(result)
87    }
88
89    pub fn interior(self, vector: Vector<3>) -> Result<Self, AnalysisError> {
90        self.validate()?;
91        vector.validate()?;
92        let mut result = Self::zero();
93        for mask in 0..EXTERIOR3_COMPONENTS {
94            let coefficient = self.coefficients[mask];
95            for axis in 0..3 {
96                let bit = 1usize << axis;
97                if mask & bit == 0 {
98                    continue;
99                }
100                let lower = (mask & (bit - 1)).count_ones();
101                let sign = if lower & 1 == 0 { 1.0 } else { -1.0 };
102                result.coefficients[mask ^ bit] += sign * vector.data[axis] * coefficient;
103            }
104        }
105        Ok(result)
106    }
107
108    /// Euclidean, positively oriented Hodge star.
109    pub fn hodge_star(self) -> Result<Self, AnalysisError> {
110        self.validate()?;
111        let mut result = Self::zero();
112        for mask in 0..EXTERIOR3_COMPONENTS {
113            let complement = 0b111 ^ mask;
114            let (_, sign) = wedge_basis(mask, complement);
115            result.coefficients[complement] += sign * self.coefficients[mask];
116        }
117        Ok(result)
118    }
119}
120
121fn wedge_basis(left: usize, right: usize) -> (usize, f64) {
122    let mut inversions = 0_u32;
123    for left_axis in 0..3 {
124        if left & (1 << left_axis) == 0 {
125            continue;
126        }
127        inversions += (right & ((1 << left_axis) - 1)).count_ones();
128    }
129    (left | right, if inversions & 1 == 0 { 1.0 } else { -1.0 })
130}
131
132pub fn permutation_parity(permutation: &[usize]) -> Result<i8, AnalysisError> {
133    for (index, value) in permutation.iter().enumerate() {
134        if *value >= permutation.len()
135            || permutation[index + 1..].iter().any(|other| other == value)
136        {
137            return Err(AnalysisError::InvalidDomain);
138        }
139    }
140    let mut inversions = 0usize;
141    for left in 0..permutation.len() {
142        for right in left + 1..permutation.len() {
143            inversions += usize::from(permutation[left] > permutation[right]);
144        }
145    }
146    Ok(if inversions & 1 == 0 { 1 } else { -1 })
147}
148
149pub fn determinant_from_wedge(columns: [Vector<3>; 3]) -> Result<f64, AnalysisError> {
150    let mut forms = [Exterior3::zero(); 3];
151    for column in 0..3 {
152        for row in 0..3 {
153            forms[column].coefficients[1 << row] = columns[column].data[row];
154        }
155    }
156    Ok(forms[0].wedge(forms[1])?.wedge(forms[2])?.coefficients[0b111])
157}
158
159pub fn pullback_one_form(
160    derivative: &LinearMap<3, 3>,
161    one_form: Vector<3>,
162) -> Result<Vector<3>, AnalysisError> {
163    derivative.transpose().apply(one_form)
164}
165
166#[cfg(test)]
167mod tests {
168    use super::*;
169
170    #[test]
171    fn wedge_is_graded_commutative_and_associative() {
172        let e1 = Exterior3::basis(0b001, 1.0);
173        let e2 = Exterior3::basis(0b010, 1.0);
174        let e3 = Exterior3::basis(0b100, 1.0);
175        assert_eq!(e1.wedge(e1).unwrap(), Exterior3::zero());
176        assert_eq!(
177            e1.wedge(e2).unwrap(),
178            e2.wedge(e1).unwrap().scale(-1.0).unwrap()
179        );
180        assert_eq!(
181            e1.wedge(e2).unwrap().wedge(e3).unwrap(),
182            e1.wedge(e2.wedge(e3).unwrap()).unwrap()
183        );
184    }
185
186    #[test]
187    fn determinant_and_hodge_identities_hold() {
188        let columns = [
189            Vector::new([2.0, 0.0, 0.0]),
190            Vector::new([1.0, 3.0, 0.0]),
191            Vector::new([0.0, 2.0, 4.0]),
192        ];
193        assert_eq!(determinant_from_wedge(columns).unwrap(), 24.0);
194
195        for mask in 0..8 {
196            let blade = Exterior3::basis(mask, 1.0);
197            let twice = blade.hodge_star().unwrap().hodge_star().unwrap();
198            let grade = mask.count_ones();
199            let sign = if (grade * (3 - grade)) & 1 == 0 {
200                1.0
201            } else {
202                -1.0
203            };
204            assert_eq!(twice, blade.scale(sign).unwrap());
205        }
206    }
207
208    #[test]
209    fn permutation_and_pullback_are_correct() {
210        assert_eq!(permutation_parity(&[2, 0, 1]), Ok(1));
211        assert_eq!(permutation_parity(&[1, 0, 2]), Ok(-1));
212        assert_eq!(
213            permutation_parity(&[0, 0, 1]),
214            Err(AnalysisError::InvalidDomain)
215        );
216
217        let derivative = LinearMap::new([[2.0, 0.0, 0.0], [1.0, 3.0, 0.0], [0.0, 0.0, 4.0]]);
218        let pulled = pullback_one_form(&derivative, Vector::new([1.0, 2.0, 3.0])).unwrap();
219        assert_eq!(pulled, Vector::new([4.0, 6.0, 12.0]));
220    }
221}
qualia_core_db/solvers/calculus/exterior.rs

qualia_core_db/solvers/calculus/
exterior.rs
