versorc_math/
multivector.rs

1// ============================================================================
2// VERSORC-MATH: GENERALIZED CLIFFORD MULTIVECTOR IMPLEMENTATION
3// License: Apache-2.0
4// ============================================================================
5
6use std::ops::{Add, AddAssign, Mul, Neg, Sub, SubAssign};
7use crate::blade::Blade;
8
9/// Fixed-size stack-allocated Multivector in Clifford Algebra Cl(p, q, r).
10/// - P, Q, R: Metric signature.
11/// - N: Algebra dimension (2^(P + Q + R)).
12#[derive(Debug, Clone, Copy, PartialEq)]
13pub struct Multivector<const P: usize, const Q: usize, const R: usize, const N: usize> {
14    pub coords: [f32; N],
15}
16
17impl<const P: usize, const Q: usize, const R: usize, const N: usize> Default
18    for Multivector<P, Q, R, N>
19{
20    #[inline(always)]
21    fn default() -> Self {
22        Self::zero()
23    }
24}
25
26impl<const P: usize, const Q: usize, const R: usize, const N: usize> Multivector<P, Q, R, N> {
27    /// Zero multivector.
28    #[inline(always)]
29    pub const fn zero() -> Self {
30        Self { coords: [0.0f32; N] }
31    }
32
33    /// Scalar multivector: alpha * 1.
34    #[inline]
35    pub fn from_scalar(alpha: f32) -> Self {
36        let mut mv = Self::zero();
37        mv.coords[0] = alpha;
38        mv
39    }
40
41    /// Creates a 1-vector from coordinate slice.
42    #[inline]
43    pub fn from_vector(vec: &[f32]) -> Self {
44        let dim = P + Q + R;
45        assert!(vec.len() <= dim, "Vector dimension exceeds algebra basis");
46        let mut mv = Self::zero();
47        for (i, &val) in vec.iter().enumerate() {
48            mv.coords[1 << i] = val;
49        }
50        mv
51    }
52
53    /// Basis element e_i (1-indexed: e_1, e_2, ...).
54    #[inline]
55    pub fn basis(i: usize) -> Self {
56        assert!(i >= 1 && i <= P + Q + R, "Basis index out of range");
57        let mut mv = Self::zero();
58        mv.coords[1 << (i - 1)] = 1.0;
59        mv
60    }
61
62    /// Scalar part <A>_0.
63    #[inline(always)]
64    pub fn scalar(&self) -> f32 {
65        self.coords[0]
66    }
67
68    /// Extracts the k-th grade projection <A>_k.
69    #[inline]
70    pub fn grade_projection(&self, grade: usize) -> Self {
71        let mut res = Self::zero();
72        for (mask, &coeff) in self.coords.iter().enumerate() {
73            if (mask as u32).count_ones() as usize == grade {
74                res.coords[mask] = coeff;
75            }
76        }
77        res
78    }
79
80    /// Clifford Reverse: A~ = sum (-1)^{k(k-1)/2} <A>_k.
81    #[inline]
82    pub fn reverse(&self) -> Self {
83        let mut res = Self::zero();
84        for (mask, &coeff) in self.coords.iter().enumerate() {
85            let sign = Blade::new(mask as u32).reverse_sign();
86            res.coords[mask] = coeff * (sign as f32);
87        }
88        res
89    }
90
91    /// Grade Involution: A^ = sum (-1)^k <A>_k.
92    #[inline]
93    pub fn grade_involution(&self) -> Self {
94        let mut res = Self::zero();
95        for (mask, &coeff) in self.coords.iter().enumerate() {
96            let sign = Blade::new(mask as u32).involution_sign();
97            res.coords[mask] = coeff * (sign as f32);
98        }
99        res
100    }
101
102    /// Clifford Conjugation: \bar{A} = sum (-1)^{k(k+1)/2} <A>_k.
103    #[inline]
104    pub fn conjugation(&self) -> Self {
105        let mut res = Self::zero();
106        for (mask, &coeff) in self.coords.iter().enumerate() {
107            let sign = Blade::new(mask as u32).conjugation_sign();
108            res.coords[mask] = coeff * (sign as f32);
109        }
110        res
111    }
112
113    /// Computes the Clifford Geometric Product: C = A * B.
114    #[inline]
115    pub fn geometric_product(&self, rhs: &Self) -> Self {
116        let mut res = Self::zero();
117        for (i, &a) in self.coords.iter().enumerate() {
118            if a == 0.0 {
119                continue;
120            }
121            let blade_a = Blade::new(i as u32);
122            for (j, &b) in rhs.coords.iter().enumerate() {
123                if b == 0.0 {
124                    continue;
125                }
126                let blade_b = Blade::new(j as u32);
127                if let Some((blade_c, sign)) = blade_a.geometric_product::<P, Q, R>(blade_b) {
128                    res.coords[blade_c.0 as usize] += a * b * (sign as f32);
129                }
130            }
131        }
132        res
133    }
134
135    /// Computes the Exterior Wedge Product: C = A ^ B.
136    #[inline]
137    pub fn wedge(&self, rhs: &Self) -> Self {
138        let mut res = Self::zero();
139        for (i, &a) in self.coords.iter().enumerate() {
140            if a == 0.0 {
141                continue;
142            }
143            let blade_a = Blade::new(i as u32);
144            for (j, &b) in rhs.coords.iter().enumerate() {
145                if b == 0.0 {
146                    continue;
147                }
148                let blade_b = Blade::new(j as u32);
149                if let Some((blade_c, sign)) = blade_a.wedge_product(blade_b) {
150                    res.coords[blade_c.0 as usize] += a * b * (sign as f32);
151                }
152            }
153        }
154        res
155    }
156
157    /// Computes the Left Contraction: C = A _| B.
158    #[inline]
159    pub fn left_contraction(&self, rhs: &Self) -> Self {
160        let mut res = Self::zero();
161        for (i, &a) in self.coords.iter().enumerate() {
162            if a == 0.0 {
163                continue;
164            }
165            let blade_a = Blade::new(i as u32);
166            for (j, &b) in rhs.coords.iter().enumerate() {
167                if b == 0.0 {
168                    continue;
169                }
170                let blade_b = Blade::new(j as u32);
171                if let Some((blade_c, sign)) = blade_a.left_contraction::<P, Q, R>(blade_b) {
172                    res.coords[blade_c.0 as usize] += a * b * (sign as f32);
173                }
174            }
175        }
176        res
177    }
178
179    /// Computes the Right Contraction: C = A |_ B.
180    #[inline]
181    pub fn right_contraction(&self, rhs: &Self) -> Self {
182        let mut res = Self::zero();
183        for (i, &a) in self.coords.iter().enumerate() {
184            if a == 0.0 {
185                continue;
186            }
187            let blade_a = Blade::new(i as u32);
188            for (j, &b) in rhs.coords.iter().enumerate() {
189                if b == 0.0 {
190                    continue;
191                }
192                let blade_b = Blade::new(j as u32);
193                if let Some((blade_c, sign)) = blade_a.right_contraction::<P, Q, R>(blade_b) {
194                    res.coords[blade_c.0 as usize] += a * b * (sign as f32);
195                }
196            }
197        }
198        res
199    }
200
201    /// Scalar product <A, B> = <A B~>_0.
202    #[inline]
203    pub fn scalar_product(&self, rhs: &Self) -> f32 {
204        self.geometric_product(&rhs.reverse()).scalar()
205    }
206
207    /// Squared magnitude: ||A||^2 = <A A~>_0.
208    #[inline]
209    pub fn norm_sq(&self) -> f32 {
210        self.scalar_product(self)
211    }
212
213    /// Magnitude / L2 norm.
214    #[inline]
215    pub fn norm(&self) -> f32 {
216        self.norm_sq().max(0.0).sqrt()
217    }
218
219    /// Normalized multivector: A / ||A||.
220    #[inline]
221    pub fn normalize(&self) -> Self {
222        let n = self.norm();
223        if n > 1e-12 {
224            *self * (1.0 / n)
225        } else {
226            *self
227        }
228    }
229
230    /// Versor sandwich product: X' = V X V~.
231    #[inline]
232    pub fn sandwich(&self, x: &Self) -> Self {
233        self.geometric_product(x).geometric_product(&self.reverse())
234    }
235
236    /// Returns the unit pseudoscalar I = e_1 ^ e_2 ^ ... ^ e_n.
237    #[inline]
238    pub fn pseudoscalar() -> Self {
239        let mut mv = Self::zero();
240        mv.coords[N - 1] = 1.0;
241        mv
242    }
243
244    /// Computes the pseudoscalar inverse I^{-1} for non-degenerate metrics (R == 0).
245    /// Under signature (P, Q, 0): I^2 = (-1)^{n(n-1)/2 + Q}, so I^{-1} = I * (-1)^{n(n-1)/2 + Q}.
246    #[inline]
247    pub fn pseudoscalar_inverse() -> Option<Self> {
248        if R > 0 {
249            return None; // Degenerate metric has I^2 == 0
250        }
251        let n = P + Q;
252        let sign_power = (n * (n - 1) / 2 + Q) % 2;
253        let sign = if sign_power == 0 { 1.0f32 } else { -1.0f32 };
254        let mut mv = Self::zero();
255        mv.coords[N - 1] = sign;
256        Some(mv)
257    }
258
259    /// Computes the Hodge Dual: A* = A I^{-1} (Macdonald Def 6.16, Hestenes Section 1.2).
260    #[inline]
261    pub fn dual(&self) -> Option<Self> {
262        let inv_i = Self::pseudoscalar_inverse()?;
263        Some(self.geometric_product(&inv_i))
264    }
265
266    /// Computes the Undual (reverse duality): A I.
267    #[inline]
268    pub fn undual(&self) -> Self {
269        self.geometric_product(&Self::pseudoscalar())
270    }
271
272    /// Computes the Lie Commutator Product: [A, B] = 1/2 (A B - B A) (Hestenes Eq 1.55).
273    /// Anticommutative: [A, B] = -[B, A]. Satisfies the Jacobi identity.
274    #[inline]
275    pub fn commutator(&self, rhs: &Self) -> Self {
276        (*self * *rhs - *rhs * *self) * 0.5
277    }
278
279    /// Computes the Jordan Anticommutator (symmetric product): A ⊙ B = 1/2 (A B + B A).
280    #[inline]
281    pub fn anticommutator(&self, rhs: &Self) -> Self {
282        (*self * *rhs + *rhs * *self) * 0.5
283    }
284
285    /// Computes the Regressive (Meet) Product: A ∨ B = (A* ∧ B*).undual()
286    /// In non-degenerate spaces, gives the geometric intersection of blade subspaces.
287    #[inline]
288    pub fn regressive_product(&self, rhs: &Self) -> Option<Self> {
289        let dual_a = self.dual()?;
290        let dual_b = rhs.dual()?;
291        let wedge_dual = dual_a.wedge(&dual_b);
292        Some(wedge_dual.undual())
293    }
294
295    /// Computes the Hestenes Inner Product (Fat Dot): A • B = sum_{r>=1, s>=1} < <A>_r <B>_s >_{|r - s|}.
296    /// As defined in Hestenes & Sobczyk (1984), Section 1.1, Eq. 1.8.
297    #[inline]
298    pub fn hestenes_inner_product(&self, rhs: &Self) -> Self {
299        let dim = P + Q + R;
300        let mut res = Self::zero();
301        for r in 1..=dim {
302            let a_r = self.grade_projection(r);
303            if a_r.norm_sq() == 0.0 {
304                continue;
305            }
306            for s in 1..=dim {
307                let b_s = rhs.grade_projection(s);
308                if b_s.norm_sq() == 0.0 {
309                    continue;
310                }
311                let target_grade = if r >= s { r - s } else { s - r };
312                let prod = a_r.geometric_product(&b_s);
313                res += prod.grade_projection(target_grade);
314            }
315        }
316        res
317    }
318
319    /// Computes the exact closed-form algebraic inverse A^{-1} for multivectors:
320    /// - For Cl(2, 0, 0): A^{-1} = \bar{A} / (A \bar{A})_0
321    /// - For Cl(3, 0, 0): A^{-1} = \bar{A} (A \bar{A})~ / det(A)
322    /// - General fallback for versors: A~ / ||A||^2
323    /// Returns None if multivector is singular (det(A) == 0).
324    #[inline]
325    pub fn inverse(&self) -> Option<Self> {
326        if N == 4 && P == 2 && Q == 0 && R == 0 {
327            // Cl(2) exact inverse
328            let conj = self.conjugation();
329            let denom = (*self * conj).scalar();
330            if denom.abs() > 1e-12 {
331                Some(conj * (1.0 / denom))
332            } else {
333                None
334            }
335        } else if N == 8 && P == 3 && Q == 0 && R == 0 {
336            // Cl(3) exact closed-form inverse
337            let conj = self.conjugation();
338            let c = *self * conj; // C in grade 0 and 3
339            let c_rev = c.reverse();
340            let det = (c * c_rev).scalar();
341            if det.abs() > 1e-12 {
342                let num = conj * c_rev;
343                Some(num * (1.0 / det))
344            } else {
345                None
346            }
347        } else {
348            // General versor fallback: if A A~ is a nonzero scalar
349            let rev = self.reverse();
350            let norm_sq = self.norm_sq();
351            if norm_sq.abs() > 1e-12 {
352                let candidate = rev * (1.0 / norm_sq);
353                let test = *self * candidate;
354                if (test.scalar() - 1.0).abs() < 1e-4 && (test - Self::from_scalar(1.0)).norm() < 1e-4 {
355                    return Some(candidate);
356                }
357            }
358            None
359        }
360    }
361}
362
363// Operator overloads
364impl<const P: usize, const Q: usize, const R: usize, const N: usize> Add for Multivector<P, Q, R, N> {
365    type Output = Self;
366    #[inline]
367    fn add(self, rhs: Self) -> Self {
368        let mut coords = [0.0f32; N];
369        for i in 0..N {
370            coords[i] = self.coords[i] + rhs.coords[i];
371        }
372        Self { coords }
373    }
374}
375
376impl<const P: usize, const Q: usize, const R: usize, const N: usize> AddAssign for Multivector<P, Q, R, N> {
377    #[inline]
378    fn add_assign(&mut self, rhs: Self) {
379        for i in 0..N {
380            self.coords[i] += rhs.coords[i];
381        }
382    }
383}
384
385impl<const P: usize, const Q: usize, const R: usize, const N: usize> Sub for Multivector<P, Q, R, N> {
386    type Output = Self;
387    #[inline]
388    fn sub(self, rhs: Self) -> Self {
389        let mut coords = [0.0f32; N];
390        for i in 0..N {
391            coords[i] = self.coords[i] - rhs.coords[i];
392        }
393        Self { coords }
394    }
395}
396
397impl<const P: usize, const Q: usize, const R: usize, const N: usize> SubAssign for Multivector<P, Q, R, N> {
398    #[inline]
399    fn sub_assign(&mut self, rhs: Self) {
400        for i in 0..N {
401            self.coords[i] -= rhs.coords[i];
402        }
403    }
404}
405
406impl<const P: usize, const Q: usize, const R: usize, const N: usize> Mul<f32> for Multivector<P, Q, R, N> {
407    type Output = Self;
408    #[inline]
409    fn mul(self, scalar: f32) -> Self {
410        let mut coords = [0.0f32; N];
411        for i in 0..N {
412            coords[i] = self.coords[i] * scalar;
413        }
414        Self { coords }
415    }
416}
417
418impl<const P: usize, const Q: usize, const R: usize, const N: usize> Mul for Multivector<P, Q, R, N> {
419    type Output = Self;
420    #[inline]
421    fn mul(self, rhs: Self) -> Self {
422        self.geometric_product(&rhs)
423    }
424}
425
426impl<const P: usize, const Q: usize, const R: usize, const N: usize> Neg for Multivector<P, Q, R, N> {
427    type Output = Self;
428    #[inline]
429    fn neg(self) -> Self {
430        let mut coords = [0.0f32; N];
431        for i in 0..N {
432            coords[i] = -self.coords[i];
433        }
434        Self { coords }
435    }
436}
437
438// Canonical Type Aliases
439pub type Cl2 = Multivector<2, 0, 0, 4>;
440pub type Cl3 = Multivector<3, 0, 0, 8>;
441pub type ClMinkowski = Multivector<1, 3, 0, 16>;
442pub type ClPGA3D = Multivector<3, 0, 1, 16>;
443pub type ClCGA3D = Multivector<4, 1, 0, 32>;
444pub type ClHamilton = Multivector<0, 2, 0, 4>;
versorc_math/multivector.rs

versorc_math/
multivector.rs
