Skip to main content

versorc_math/
blade.rs

1// ============================================================================
2// VERSORC-MATH: BASIS BLADES & CAYLEY MULTIPLICATION SIGN CALCULATOR
3// License: Apache-2.0
4// ============================================================================
5
6use crate::signature::MetricSignature;
7
8/// Basis blade represented as a bitmask of generator vectors.
9/// E.g., scalar = 0b0000, e_1 = 0b0001, e_2 = 0b0010, e_12 = 0b0011, e_123 = 0b0111.
10#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Default)]
11pub struct Blade(pub u32);
12
13impl Blade {
14    pub const SCALAR: Self = Self(0);
15
16    #[inline(always)]
17    pub const fn new(mask: u32) -> Self {
18        Self(mask)
19    }
20
21    #[inline(always)]
22    pub const fn generator(i: usize) -> Self {
23        Self(1 << i)
24    }
25
26    /// Grade of the blade (number of basis vectors spanning it).
27    #[inline(always)]
28    pub const fn grade(self) -> usize {
29        self.0.count_ones() as usize
30    }
31
32    /// Reverse sign of the blade: (-1)^{k(k-1)/2}.
33    #[inline(always)]
34    pub const fn reverse_sign(self) -> i8 {
35        let k = self.grade();
36        if k < 2 {
37            1
38        } else if (k * (k - 1) / 2) % 2 == 0 {
39            1
40        } else {
41            -1
42        }
43    }
44
45    /// Grade involution sign of the blade: (-1)^k.
46    #[inline(always)]
47    pub const fn involution_sign(self) -> i8 {
48        if self.grade() % 2 == 0 {
49            1
50        } else {
51            -1
52        }
53    }
54
55    /// Clifford conjugation sign of the blade: (-1)^{k(k+1)/2}.
56    #[inline(always)]
57    pub const fn conjugation_sign(self) -> i8 {
58        let k = self.grade();
59        if (k * (k + 1) / 2) % 2 == 0 {
60            1
61        } else {
62            -1
63        }
64    }
65
66    /// Computes the geometric product of two basis blades e_I * e_J under metric Cl(p, q, r).
67    /// Returns Some((resulting_blade, sign)) or None if the result is zero (due to degenerate r-generators).
68    #[inline]
69    pub fn geometric_product<const P: usize, const Q: usize, const R: usize>(
70        self,
71        rhs: Self,
72    ) -> Option<(Self, i8)> {
73        let mask_a = self.0;
74        let mask_b = rhs.0;
75
76        // 1. Calculate permutation transpositions (swaps) to order generators
77        let mut swaps = 0u32;
78        let mut temp_b = mask_b;
79        while temp_b != 0 {
80            let j = temp_b.trailing_zeros();
81            // Count generators in A that are strictly greater than generator j
82            swaps += (mask_a >> (j + 1)).count_ones();
83            temp_b &= temp_b - 1; // Clear lowest set bit
84        }
85
86        let mut sign: i8 = if swaps % 2 == 0 { 1 } else { -1 };
87
88        // 2. Evaluate metric contraction for common generators: e_k * e_k = eta_kk
89        let common = mask_a & mask_b;
90        if common != 0 {
91            let mut temp_common = common;
92            while temp_common != 0 {
93                let k = temp_common.trailing_zeros() as usize;
94                let metric_sq = MetricSignature::<P, Q, R>::basis_square(k);
95                if metric_sq == 0 {
96                    return None; // Degenerate generator squared to zero
97                }
98                sign *= metric_sq;
99                temp_common &= temp_common - 1;
100            }
101        }
102
103        let resulting_mask = mask_a ^ mask_b;
104        Some((Self(resulting_mask), sign))
105    }
106
107    /// Computes the exterior wedge product e_I ^ e_J.
108    /// Nonzero if and only if blades share no common generators (I & J == 0).
109    #[inline]
110    pub fn wedge_product(self, rhs: Self) -> Option<(Self, i8)> {
111        if (self.0 & rhs.0) != 0 {
112            return None; // Linearly dependent generators wedge to zero
113        }
114
115        // Swaps identical to geometric product when common == 0
116        let mut swaps = 0u32;
117        let mut temp_b = rhs.0;
118        while temp_b != 0 {
119            let j = temp_b.trailing_zeros();
120            swaps += (self.0 >> (j + 1)).count_ones();
121            temp_b &= temp_b - 1;
122        }
123
124        let sign: i8 = if swaps % 2 == 0 { 1 } else { -1 };
125        Some((Self(self.0 | rhs.0), sign))
126    }
127
128    /// Computes the left contraction e_I _| e_J under metric Cl(p, q, r).
129    /// Nonzero only if grade(result) == grade(J) - grade(I).
130    #[inline]
131    pub fn left_contraction<const P: usize, const Q: usize, const R: usize>(
132        self,
133        rhs: Self,
134    ) -> Option<(Self, i8)> {
135        let r = self.grade();
136        let s = rhs.grade();
137        if r > s {
138            return None;
139        }
140        let (prod_blade, sign) = self.geometric_product::<P, Q, R>(rhs)?;
141        if prod_blade.grade() == s - r {
142            Some((prod_blade, sign))
143        } else {
144            None
145        }
146    }
147
148    /// Computes the right contraction e_I |_ e_J under metric Cl(p, q, r).
149    /// Nonzero only if grade(result) == grade(I) - grade(J).
150    #[inline]
151    pub fn right_contraction<const P: usize, const Q: usize, const R: usize>(
152        self,
153        rhs: Self,
154    ) -> Option<(Self, i8)> {
155        let r = self.grade();
156        let s = rhs.grade();
157        if r < s {
158            return None;
159        }
160        let (prod_blade, sign) = self.geometric_product::<P, Q, R>(rhs)?;
161        if prod_blade.grade() == r - s {
162            Some((prod_blade, sign))
163        } else {
164            None
165        }
166    }
167}