Skip to main content

versorc_math/lattice/
e8.rs

1// ============================================================================
2// VERSORC-MATH: GOSSET E8 ROOT LATTICE & CONWAY-SLOANE FAST DECODER
3// License: Apache-2.0
4//
5// E8 is the unique even unimodular lattice in dimension 8.
6// Composed of the union of two D8 cosets: E8 = D8 U (D8 + 1/2 * 1).
7// Possesses exactly 240 minimal root versors with norm squared ||r||^2 = 2.0.
8// ============================================================================
9
10pub const E8_DIMENSION: usize = 8;
11pub const E8_NUM_ROOTS: usize = 240;
12pub const E8_CODEBOOK_SIZE: usize = 1024; // Exactly 1 KiB
13pub const E8_NUM_CENTROIDS: usize = E8_CODEBOOK_SIZE / E8_DIMENSION; // 128 centroids in 8D
14
15/// E8 Lattice Codebook of 1 KiB.
16/// Explicitly aligned to 64 bytes to guarantee residence in exactly 16 contiguous L1d cache lines.
17#[repr(C, align(64))]
18#[derive(Debug, Clone, PartialEq)]
19pub struct E8LatticeCodebook {
20    /// 128 centroids x 8 integer coordinates (scaled int8) = 1024 bytes
21    pub table: [i8; E8_CODEBOOK_SIZE],
22}
23
24const _: () = assert!(
25    std::mem::size_of::<E8LatticeCodebook>() == E8_CODEBOOK_SIZE,
26    "E8LatticeCodebook must measure exactly 1024 bytes (1 KiB)"
27);
28
29impl Default for E8LatticeCodebook {
30    fn default() -> Self {
31        Self::new()
32    }
33}
34
35impl E8LatticeCodebook {
36    /// Initializes E8P codebook with canonical E8 half-integer roots scaled by 2 (coordinates in {+1, -1}).
37    /// Every 8D centroid has an even number of minus signs, satisfying E8 membership.
38    pub fn new() -> Self {
39        let mut table = [0i8; E8_CODEBOOK_SIZE];
40        let mut idx = 0;
41        for c in 0..E8_NUM_CENTROIDS {
42            let mut parity = 0;
43            for d in 0..7 {
44                let sign = if ((c >> d) & 1) == 0 { 1i8 } else { -1i8 };
45                if sign == -1 {
46                    parity ^= 1;
47                }
48                table[idx] = sign;
49                idx += 1;
50            }
51            // 8th coordinate enforces even number of negative signs
52            table[idx] = if parity == 0 { 1i8 } else { -1i8 };
53            idx += 1;
54        }
55        Self { table }
56    }
57
58    /// Decodes a 4-bit nibble [0..15] back into its 2D constituent weights (w0, w1).
59    /// Invertible bijection with project_and_encode_to_nibbles. Unbiased symmetric levels in [-1.5, +1.5].
60    #[inline(always)]
61    pub fn decode_nibble(nibble: u8) -> (f32, f32) {
62        let q0 = (nibble & 0x03) as f32;
63        let q1 = ((nibble >> 2) & 0x03) as f32;
64        let w0 = q0 - 1.5;
65        let w1 = q1 - 1.5;
66        (w0, w1)
67    }
68
69    /// Decodes 4 4-bit nibbles (2 bytes) back to an 8-dimensional E8 vector.
70    #[inline(always)]
71    pub fn decode_4_nibbles(nibbles: &[u8; 4], scale: f32) -> [f32; E8_DIMENSION] {
72        let mut out = [0.0f32; E8_DIMENSION];
73        for k in 0..4 {
74            let (w0, w1) = Self::decode_nibble(nibbles[k]);
75            out[2 * k] = w0 * scale;
76            out[2 * k + 1] = w1 * scale;
77        }
78        out
79    }
80
81    /// Projects an 8D continuous vector onto the closest vector in the Conway-Sloane half-integer E8 coset (D8 + 1/2*1),
82    /// then discretizes each 2D sub-vector to 1 of 16 symmetric zero-mean states (4 bits: 2 bits per coordinate).
83    /// Levels per coordinate: {-1.5, -0.5, +0.5, +1.5} * scale. Strictly unbiased (mean = 0.0).
84    pub fn project_and_encode_to_nibbles(x: &[f32; E8_DIMENSION], scale: f32) -> [u8; 4] {
85        let inv_scale = if scale > 1e-7 { 1.0 / scale } else { 1.0 };
86        let mut shifted_x = [0.0f32; E8_DIMENSION];
87        for d in 0..E8_DIMENSION {
88            shifted_x[d] = x[d] * inv_scale - 0.5;
89        }
90
91        let p_int = E8Lattice::project_to_d8(&shifted_x);
92        let mut nibbles = [0u8; 4];
93        for k in 0..4 {
94            let w0 = p_int[2 * k] as f32 + 0.5;
95            let w1 = p_int[2 * k + 1] as f32 + 0.5;
96            let q0 = ((w0 + 1.5).round() as i32).clamp(0, 3) as u8;
97            let q1 = ((w1 + 1.5).round() as i32).clamp(0, 3) as u8;
98            nibbles[k] = (q0 & 0x03) | ((q1 & 0x03) << 2);
99        }
100        nibbles
101    }
102}
103
104/// The 8-dimensional Gosset E8 lattice structure and Conway-Sloane closest-point decoder.
105pub struct E8Lattice;
106
107impl E8Lattice {
108    /// Generates all 240 minimal root vectors of the E8 lattice in R^8.
109    /// - 112 roots in D8: permutations of (+-1, +-1, 0, 0, 0, 0, 0, 0).
110    /// - 128 roots in D8 + 1/2*1: (+-1/2, ..., +-1/2) with an even number of minus signs.
111    /// Every root satisfies ||r||^2 = 2.0.
112    pub fn generate_240_roots() -> [[f32; E8_DIMENSION]; E8_NUM_ROOTS] {
113        let mut roots = [[0.0f32; E8_DIMENSION]; E8_NUM_ROOTS];
114        let mut idx = 0;
115
116        // 1. Generate 112 integer roots in D8: (+-1, +-1, 0^6)
117        for i in 0..E8_DIMENSION {
118            for j in (i + 1)..E8_DIMENSION {
119                for &s_i in &[1.0f32, -1.0f32] {
120                    for &s_j in &[1.0f32, -1.0f32] {
121                        let mut root = [0.0f32; E8_DIMENSION];
122                        root[i] = s_i;
123                        root[j] = s_j;
124                        roots[idx] = root;
125                        idx += 1;
126                    }
127                }
128            }
129        }
130        assert_eq!(idx, 112, "Group 1 must contain exactly 112 roots");
131
132        // 2. Generate 128 half-integer roots in D8 + 1/2*1: (+-1/2)^8 with even negative count
133        for mask in 0..256u32 {
134            // Parity of negative signs: count number of set bits in mask
135            if mask.count_ones() % 2 == 0 {
136                let mut root = [0.0f32; E8_DIMENSION];
137                for d in 0..E8_DIMENSION {
138                    let sign = if ((mask >> d) & 1) == 0 { 0.5f32 } else { -0.5f32 };
139                    root[d] = sign;
140                }
141                roots[idx] = root;
142                idx += 1;
143            }
144        }
145        assert_eq!(idx, 240, "Total E8 minimal roots must equal exactly 240");
146
147        roots
148    }
149
150    /// Fast Conway & Sloane closest vector algorithm for the D8 checkerboard lattice.
151    /// In D8, coordinates are integers and their sum is even.
152    pub fn project_to_d8(x: &[f32; E8_DIMENSION]) -> [i32; E8_DIMENSION] {
153        let mut f = [0i32; E8_DIMENSION];
154        let mut errors = [0.0f32; E8_DIMENSION];
155        let mut sum_f = 0i32;
156
157        for i in 0..E8_DIMENSION {
158            f[i] = x[i].round() as i32;
159            errors[i] = (x[i] - f[i] as f32).abs();
160            sum_f += f[i];
161        }
162
163        // If sum of rounded integers is odd, flip the coordinate with the largest error
164        if sum_f % 2 != 0 {
165            let mut max_err_idx = 0;
166            let mut max_err = errors[0];
167            for i in 1..E8_DIMENSION {
168                if errors[i] > max_err {
169                    max_err = errors[i];
170                    max_err_idx = i;
171                }
172            }
173
174            if x[max_err_idx] >= f[max_err_idx] as f32 {
175                f[max_err_idx] += 1;
176            } else {
177                f[max_err_idx] -= 1;
178            }
179        }
180
181        f
182    }
183
184    /// Fast Conway & Sloane Closest Vector Algorithm for E8:
185    /// E8 = D8 U (D8 + 1/2 * 1)
186    /// Finds the globally closest lattice point in strictly linear O(8) operations.
187    pub fn project_to_e8(x: &[f32; E8_DIMENSION]) -> [f32; E8_DIMENSION] {
188        // Candidate 1: Nearest point in D8
189        let p1_int = Self::project_to_d8(x);
190        let mut p1 = [0.0f32; E8_DIMENSION];
191        let mut dist_sq1 = 0.0f32;
192        for i in 0..E8_DIMENSION {
193            p1[i] = p1_int[i] as f32;
194            let diff = x[i] - p1[i];
195            dist_sq1 += diff * diff;
196        }
197
198        // Candidate 2: Nearest point in shifted coset D8 + 1/2 * 1
199        let mut shifted_x = [0.0f32; E8_DIMENSION];
200        for i in 0..E8_DIMENSION {
201            shifted_x[i] = x[i] - 0.5;
202        }
203        let p2_int = Self::project_to_d8(&shifted_x);
204        let mut p2 = [0.0f32; E8_DIMENSION];
205        let mut dist_sq2 = 0.0f32;
206        for i in 0..E8_DIMENSION {
207            p2[i] = p2_int[i] as f32 + 0.5;
208            let diff = x[i] - p2[i];
209            dist_sq2 += diff * diff;
210        }
211
212        // Return candidate with minimal Euclidean distance
213        if dist_sq1 <= dist_sq2 {
214            p1
215        } else {
216            p2
217        }
218    }
219}