resonance_spectral_gap/
lib.rs

1#![doc = include_str!("../README.md")]
2#![forbid(unsafe_code)]
3#![allow(clippy::doc_markdown)]
4#![allow(clippy::large_stack_arrays)]
5#![allow(clippy::large_const_arrays)]
6#![allow(clippy::module_name_repetitions)]
7
8use num_rational::Ratio;
9use num_traits::{One, Zero};
10use std::collections::{BTreeMap, HashMap, HashSet, VecDeque};
11use std::fmt;
12
13/// Exact rational number type (no floating point).
14pub type Rational = Ratio<i64>;
15
16// ─────────────────────────────────────────────────────────────────────────────
17// Atlas Graph: 96-vertex graph of Resonance Classes
18// ─────────────────────────────────────────────────────────────────────────────
19
20/// Number of vertices in the Atlas (96 = 2⁵ × 3)
21const ATLAS_VERTEX_COUNT: usize = 96;
22
23/// Number of edges in the Atlas (256)
24const ATLAS_EDGE_COUNT: usize = 256;
25
26/// Canonical label for an Atlas vertex.
27///
28/// Each vertex is a 6-tuple `(e1, e2, e3, d45, e6, e7)` where:
29/// - `e1, e2, e3, e6, e7 ∈ {0, 1}` (binary coordinates)
30/// - `d45 ∈ {-1, 0, +1}` (canonicalized difference `e4 - e5`)
31///
32/// This gives `2⁵ × 3 = 96` vertices.
33#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
34pub struct Label {
35    /// e1 coordinate (0 or 1)
36    pub e1: u8,
37    /// e2 coordinate (0 or 1)
38    pub e2: u8,
39    /// e3 coordinate (0 or 1)
40    pub e3: u8,
41    /// d45 = e4 - e5, canonicalized to {-1, 0, +1}
42    pub d45: i8,
43    /// e6 coordinate (0 or 1)
44    pub e6: u8,
45    /// e7 coordinate (0 or 1) — flipped by mirror τ
46    pub e7: u8,
47}
48
49impl Label {
50    /// Create a new label with validated coordinates.
51    ///
52    /// # Panics
53    ///
54    /// Panics if binary coordinates are not in {0, 1} or d45 is not in {-1, 0, 1}.
55    #[must_use]
56    pub const fn new(e1: u8, e2: u8, e3: u8, d45: i8, e6: u8, e7: u8) -> Self {
57        assert!(
58            e1 <= 1 && e2 <= 1 && e3 <= 1 && e6 <= 1 && e7 <= 1,
59            "Binary coordinates must be 0 or 1"
60        );
61        assert!(d45 >= -1 && d45 <= 1, "d45 must be in {{-1, 0, 1}}");
62        Self { e1, e2, e3, d45, e6, e7 }
63    }
64
65    /// Apply mirror transformation τ (flip e7).
66    #[must_use]
67    pub const fn mirror(&self) -> Self {
68        Self {
69            e1: self.e1,
70            e2: self.e2,
71            e3: self.e3,
72            d45: self.d45,
73            e6: self.e6,
74            e7: 1 - self.e7,
75        }
76    }
77}
78
79/// The Atlas of Resonance Classes: a 96-vertex, 256-edge graph.
80///
81/// Constructed from first principles by enumerating all canonical labels
82/// and connecting vertices that differ by a single coordinate flip
83/// (Hamming distance 1, excluding the e7 mirror coordinate).
84///
85/// # Structure
86///
87/// - **96 vertices**: labeled by `(e1, e2, e3, d45, e6, e7)`
88/// - **256 edges**: Hamming-1 flips of `{e1, e2, e3, e4, e5, e6}` (not e7)
89/// - **Two hemispheres**: split by `e7 ∈ {0, 1}`, no cross-edges
90/// - **Mirror symmetry**: τ flips e7, giving H₀ ≅ H₁
91#[derive(Debug, Clone)]
92pub struct Atlas {
93    labels: Vec<Label>,
94    #[allow(dead_code)]
95    label_index: HashMap<Label, usize>,
96    adjacency: Vec<HashSet<usize>>,
97}
98
99impl Atlas {
100    /// Construct the Atlas graph from first principles.
101    ///
102    /// Generates all 96 canonical labels and builds the adjacency structure.
103    /// Verifies invariants (vertex count, edge count, regularity) on construction.
104    ///
105    /// # Panics
106    ///
107    /// Panics if any structural invariant fails.
108    #[must_use]
109    pub fn new() -> Self {
110        let labels = Self::generate_labels();
111        let label_index: HashMap<Label, usize> =
112            labels.iter().enumerate().map(|(i, &lab)| (lab, i)).collect();
113        let adjacency = Self::build_adjacency(&labels, &label_index);
114
115        let atlas = Self { labels, label_index, adjacency };
116        atlas.verify_invariants();
117        atlas
118    }
119
120    /// Number of vertices (always 96).
121    #[must_use]
122    pub fn num_vertices(&self) -> usize {
123        self.labels.len()
124    }
125
126    /// Number of edges (always 256).
127    #[must_use]
128    pub fn num_edges(&self) -> usize {
129        let sum: usize = self.adjacency.iter().map(HashSet::len).sum();
130        sum / 2
131    }
132
133    /// Degree of a vertex.
134    #[must_use]
135    pub fn degree(&self, vertex: usize) -> usize {
136        self.adjacency[vertex].len()
137    }
138
139    /// Neighbors of a vertex.
140    #[must_use]
141    pub fn neighbors(&self, vertex: usize) -> &HashSet<usize> {
142        &self.adjacency[vertex]
143    }
144
145    /// Label of a vertex.
146    #[must_use]
147    pub fn label(&self, vertex: usize) -> Label {
148        self.labels[vertex]
149    }
150
151    // ─── Private construction ───────────────────────────────────────
152
153    fn generate_labels() -> Vec<Label> {
154        let mut labels = Vec::with_capacity(ATLAS_VERTEX_COUNT);
155        for e1 in 0..=1u8 {
156            for e2 in 0..=1u8 {
157                for e3 in 0..=1u8 {
158                    for e6 in 0..=1u8 {
159                        for e7 in 0..=1u8 {
160                            for d45 in -1..=1i8 {
161                                labels.push(Label::new(e1, e2, e3, d45, e6, e7));
162                            }
163                        }
164                    }
165                }
166            }
167        }
168        assert_eq!(labels.len(), ATLAS_VERTEX_COUNT);
169        labels
170    }
171
172    fn build_adjacency(
173        labels: &[Label],
174        label_index: &HashMap<Label, usize>,
175    ) -> Vec<HashSet<usize>> {
176        let mut adjacency = vec![HashSet::new(); labels.len()];
177        for (i, label) in labels.iter().enumerate() {
178            for neighbor_label in Self::compute_neighbors(*label) {
179                if let Some(&j) = label_index.get(&neighbor_label) {
180                    if i != j {
181                        adjacency[i].insert(j);
182                    }
183                }
184            }
185        }
186        adjacency
187    }
188
189    fn compute_neighbors(label: Label) -> Vec<Label> {
190        let Label { e1, e2, e3, d45, e6, e7 } = label;
191        vec![
192            Label::new(1 - e1, e2, e3, d45, e6, e7),
193            Label::new(e1, 1 - e2, e3, d45, e6, e7),
194            Label::new(e1, e2, 1 - e3, d45, e6, e7),
195            Label::new(e1, e2, e3, d45, 1 - e6, e7),
196            Label::new(e1, e2, e3, Self::flip_d45_by_e4(d45), e6, e7),
197            Label::new(e1, e2, e3, Self::flip_d45_by_e5(d45), e6, e7),
198        ]
199    }
200
201    /// d45 after flipping e4: -1→0, 0→+1, +1→0
202    fn flip_d45_by_e4(d: i8) -> i8 {
203        match d {
204            -1 | 1 => 0,
205            0 => 1,
206            _ => panic!("d45 must be in {{-1, 0, 1}}"),
207        }
208    }
209
210    /// d45 after flipping e5: -1→0, 0→-1, +1→0
211    fn flip_d45_by_e5(d: i8) -> i8 {
212        match d {
213            -1 | 1 => 0,
214            0 => -1,
215            _ => panic!("d45 must be in {{-1, 0, 1}}"),
216        }
217    }
218
219    fn verify_invariants(&self) {
220        assert_eq!(self.labels.len(), ATLAS_VERTEX_COUNT, "Atlas must have 96 vertices");
221        assert_eq!(self.num_edges(), ATLAS_EDGE_COUNT, "Atlas must have 256 edges");
222    }
223}
224
225impl Default for Atlas {
226    fn default() -> Self {
227        Self::new()
228    }
229}
230
231// ─────────────────────────────────────────────────────────────────────────────
232// Constants
233// ─────────────────────────────────────────────────────────────────────────────
234
235/// Number of vertices per hemisphere (48 = 96/2)
236const HEMISPHERE_VERTICES: usize = 48;
237
238/// Number of vertices per Q₄ hypercube block (16 = 2⁴)
239const Q4_VERTICES: usize = 16;
240
241/// Degree of each vertex in a Q₄ hypercube
242const Q4_DEGREE: usize = 4;
243
244/// Number of edges in a Q₄ hypercube (32 = 16×4/2)
245const Q4_EDGES: usize = 32;
246
247/// Q₄ Laplacian eigenvalues with multiplicities: `(eigenvalue, C(4,k))`
248///
249/// For the 4-dimensional hypercube, `L = 4I - A` has eigenvalues
250/// `2k` for `k = 0,1,2,3,4` with multiplicities `C(4,k) = 1,4,6,4,1`.
251const Q4_SPECTRUM: [(i64, usize); 5] = [
252    (0, 1),  // C(4,0) = 1
253    (2, 4),  // C(4,1) = 4
254    (4, 6),  // C(4,2) = 6
255    (6, 4),  // C(4,3) = 4
256    (8, 1),  // C(4,4) = 1
257];
258
259// ─────────────────────────────────────────────────────────────────────────────
260// Dense Rational Matrix (internal utility)
261// ─────────────────────────────────────────────────────────────────────────────
262
263#[derive(Clone)]
264struct DenseRatMatrix {
265    n: usize,
266    data: Vec<Rational>,
267}
268
269impl DenseRatMatrix {
270    fn new(n: usize) -> Self {
271        Self { n, data: vec![Rational::zero(); n * n] }
272    }
273
274    fn get(&self, i: usize, j: usize) -> Rational {
275        self.data[i * self.n + j]
276    }
277
278    fn set(&mut self, i: usize, j: usize, val: Rational) {
279        self.data[i * self.n + j] = val;
280    }
281
282    fn trace(&self) -> Rational {
283        (0..self.n).map(|i| self.get(i, i)).fold(Rational::zero(), |a, b| a + b)
284    }
285
286    fn is_symmetric(&self) -> bool {
287        for i in 0..self.n {
288            for j in (i + 1)..self.n {
289                if self.get(i, j) != self.get(j, i) {
290                    return false;
291                }
292            }
293        }
294        true
295    }
296
297    fn row_sum(&self, i: usize) -> Rational {
298        (0..self.n).map(|j| self.get(i, j)).fold(Rational::zero(), |a, b| a + b)
299    }
300
301    fn frobenius_squared(&self) -> Rational {
302        self.data.iter().fold(Rational::zero(), |a, &v| a + v * v)
303    }
304
305    fn has_nonneg_diagonal(&self) -> bool {
306        (0..self.n).all(|i| self.get(i, i) >= Rational::zero())
307    }
308
309    fn has_nonpos_off_diagonal(&self) -> bool {
310        for i in 0..self.n {
311            for j in 0..self.n {
312                if i != j && self.get(i, j) > Rational::zero() {
313                    return false;
314                }
315            }
316        }
317        true
318    }
319
320    fn from_atlas_subgraph(atlas: &Atlas, vertices: &[usize]) -> Self {
321        let n = vertices.len();
322        let mut mat = Self::new(n);
323        let index_map: HashMap<usize, usize> = vertices
324            .iter()
325            .enumerate()
326            .map(|(local, &atlas_idx)| (atlas_idx, local))
327            .collect();
328
329        for (i, &v) in vertices.iter().enumerate() {
330            let mut degree = 0i64;
331            for &neighbor in atlas.neighbors(v) {
332                if let Some(&j) = index_map.get(&neighbor) {
333                    mat.set(i, j, -Rational::one());
334                    degree += 1;
335                }
336            }
337            mat.set(i, i, Rational::from_integer(degree));
338        }
339        mat
340    }
341}
342
343impl fmt::Debug for DenseRatMatrix {
344    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
345        write!(f, "DenseRatMatrix({}×{})", self.n, self.n)
346    }
347}
348
349// ─────────────────────────────────────────────────────────────────────────────
350// 3×3 Matrix Operations
351// ─────────────────────────────────────────────────────────────────────────────
352
353/// Compute determinant of a 3×3 rational matrix (exact).
354fn det_3x3(m: &[[Rational; 3]; 3]) -> Rational {
355    m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
356        - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0])
357        + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0])
358}
359
360/// Construct the 3×3 block tridiagonal matrix M_ν.
361///
362/// ```text
363/// M_ν = [ ν+1,  -1,   0  ]
364///       [ -1,   ν+2,  -1 ]
365///       [  0,   -1,  ν+1 ]
366/// ```
367fn block_tridiagonal_matrix(nu: i64) -> [[Rational; 3]; 3] {
368    let n = Rational::from_integer(nu);
369    let one = Rational::one();
370    let two = Rational::new(2, 1);
371    let neg = -one;
372    let zero = Rational::zero();
373    [
374        [n + one, neg, zero],
375        [neg, n + two, neg],
376        [zero, neg, n + one],
377    ]
378}
379
380/// Verify eigenvalues of M_ν by checking `det(M_ν - λI) = 0`.
381///
382/// Returns `{ν, ν+1, ν+3}` after verification.
383fn verify_block_eigenvalues(nu: i64) -> [Rational; 3] {
384    let m = block_tridiagonal_matrix(nu);
385    let nu_r = Rational::from_integer(nu);
386    let eigenvalues = [
387        nu_r,
388        nu_r + Rational::one(),
389        nu_r + Rational::new(3, 1),
390    ];
391    for &lambda in &eigenvalues {
392        let shifted = [
393            [m[0][0] - lambda, m[0][1], m[0][2]],
394            [m[1][0], m[1][1] - lambda, m[1][2]],
395            [m[2][0], m[2][1], m[2][2] - lambda],
396        ];
397        let det = det_3x3(&shifted);
398        assert!(det.is_zero(), "det(M_{nu} - {lambda}·I) must be 0, got {det}");
399    }
400    eigenvalues
401}
402
403// ─────────────────────────────────────────────────────────────────────────────
404// Hemisphere and Block Decomposition
405// ─────────────────────────────────────────────────────────────────────────────
406
407fn decompose_hemispheres(atlas: &Atlas) -> [Vec<usize>; 2] {
408    let mut h0 = Vec::with_capacity(HEMISPHERE_VERTICES);
409    let mut h1 = Vec::with_capacity(HEMISPHERE_VERTICES);
410    for v in 0..atlas.num_vertices() {
411        if atlas.label(v).e7 == 0 { h0.push(v); } else { h1.push(v); }
412    }
413    assert_eq!(h0.len(), HEMISPHERE_VERTICES);
414    assert_eq!(h1.len(), HEMISPHERE_VERTICES);
415    [h0, h1]
416}
417
418fn verify_no_cross_hemisphere_edges(atlas: &Atlas, hemispheres: &[Vec<usize>; 2]) {
419    let h0_set: HashSet<usize> = hemispheres[0].iter().copied().collect();
420    for &v in &hemispheres[0] {
421        for &n in atlas.neighbors(v) {
422            assert!(h0_set.contains(&n), "Cross-hemisphere edge: {v} ~ {n}");
423        }
424    }
425}
426
427fn decompose_blocks(atlas: &Atlas, hemisphere: &[usize]) -> [Vec<usize>; 3] {
428    let mut b_neg = Vec::with_capacity(Q4_VERTICES);
429    let mut b_zero = Vec::with_capacity(Q4_VERTICES);
430    let mut b_pos = Vec::with_capacity(Q4_VERTICES);
431    for &v in hemisphere {
432        match atlas.label(v).d45 {
433            -1 => b_neg.push(v),
434            0 => b_zero.push(v),
435            1 => b_pos.push(v),
436            d => panic!("Invalid d45 value: {d}"),
437        }
438    }
439    assert_eq!(b_neg.len(), Q4_VERTICES);
440    assert_eq!(b_zero.len(), Q4_VERTICES);
441    assert_eq!(b_pos.len(), Q4_VERTICES);
442    [b_neg, b_zero, b_pos]
443}
444
445fn verify_q4_block(atlas: &Atlas, block: &[usize]) {
446    assert_eq!(block.len(), Q4_VERTICES);
447    let block_set: HashSet<usize> = block.iter().copied().collect();
448
449    for &v in block {
450        let within_degree = atlas.neighbors(v).iter().filter(|n| block_set.contains(n)).count();
451        assert_eq!(within_degree, Q4_DEGREE, "Q4 vertex {v} degree {within_degree} != {Q4_DEGREE}");
452    }
453
454    // Connected via BFS
455    let mut visited = HashSet::new();
456    let mut queue = VecDeque::new();
457    visited.insert(block[0]);
458    queue.push_back(block[0]);
459    while let Some(v) = queue.pop_front() {
460        for &n in atlas.neighbors(v) {
461            if block_set.contains(&n) && visited.insert(n) {
462                queue.push_back(n);
463            }
464        }
465    }
466    assert_eq!(visited.len(), Q4_VERTICES, "Q4 block not connected");
467
468    // Bipartite via parity
469    for &v in block {
470        let vl = atlas.label(v);
471        let v_parity = (vl.e1 + vl.e2 + vl.e3 + vl.e6) % 2;
472        for &n in atlas.neighbors(v) {
473            if block_set.contains(&n) {
474                let nl = atlas.label(n);
475                let n_parity = (nl.e1 + nl.e2 + nl.e3 + nl.e6) % 2;
476                assert_ne!(v_parity, n_parity, "Q4 bipartite violation: {v} ~ {n}");
477            }
478        }
479    }
480}
481
482fn verify_inter_block_structure(atlas: &Atlas, blocks: &[Vec<usize>; 3]) {
483    let block_sets: [HashSet<usize>; 3] = [
484        blocks[0].iter().copied().collect(),
485        blocks[1].iter().copied().collect(),
486        blocks[2].iter().copied().collect(),
487    ];
488    verify_identity_matching(atlas, &blocks[0], &block_sets[1]);
489    verify_identity_matching(atlas, &blocks[2], &block_sets[1]);
490
491    // No skip edges (d45=-1 ↔ d45=+1)
492    for &v in &blocks[0] {
493        for &n in atlas.neighbors(v) {
494            assert!(!block_sets[2].contains(&n), "Skip edge: {v} ~ {n}");
495        }
496    }
497}
498
499fn verify_identity_matching(atlas: &Atlas, source: &[usize], target: &HashSet<usize>) {
500    for &v in source {
501        let inter: Vec<usize> = atlas.neighbors(v).iter().filter(|n| target.contains(n)).copied().collect();
502        assert_eq!(inter.len(), 1, "Vertex {v}: {} inter-block neighbors", inter.len());
503        let vl = atlas.label(v);
504        let nl = atlas.label(inter[0]);
505        assert_eq!(vl.e1, nl.e1);
506        assert_eq!(vl.e2, nl.e2);
507        assert_eq!(vl.e3, nl.e3);
508        assert_eq!(vl.e6, nl.e6);
509    }
510}
511
512#[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
513fn verify_q4_spectrum(atlas: &Atlas, block: &[usize]) {
514    let block_set: HashSet<usize> = block.iter().copied().collect();
515    let mut degree_sum: usize = 0;
516    let mut degree_sq_sum: usize = 0;
517    for &v in block {
518        let deg = atlas.neighbors(v).iter().filter(|n| block_set.contains(n)).count();
519        degree_sum += deg;
520        degree_sq_sum += deg * deg;
521    }
522    let edge_count = degree_sum / 2;
523    assert_eq!(degree_sum, Q4_VERTICES * Q4_DEGREE);
524    assert_eq!(edge_count, Q4_EDGES);
525    let trace = degree_sum;
526    let expected_trace: usize = Q4_SPECTRUM.iter().map(|&(ev, mult)| ev as usize * mult).sum();
527    assert_eq!(trace, expected_trace);
528    let trace_sq = degree_sq_sum + 2 * edge_count;
529    let expected_trace_sq: usize = Q4_SPECTRUM.iter().map(|&(ev, mult)| (ev as usize) * (ev as usize) * mult).sum();
530    assert_eq!(trace_sq, expected_trace_sq);
531}
532
533fn compute_hemisphere_spectrum() -> Vec<(Rational, usize)> {
534    let mut spectrum_map: BTreeMap<Rational, usize> = BTreeMap::new();
535    for &(nu, q4_mult) in &Q4_SPECTRUM {
536        let block_eigs = verify_block_eigenvalues(nu);
537        for &eig in &block_eigs {
538            *spectrum_map.entry(eig).or_insert(0) += q4_mult;
539        }
540    }
541    spectrum_map.into_iter().collect()
542}
543
544// ─────────────────────────────────────────────────────────────────────────────
545// SpectralAnalysis: Main Public Interface
546// ─────────────────────────────────────────────────────────────────────────────
547
548/// Spectral analysis of the Atlas graph Laplacian.
549///
550/// Computes the complete spectrum through a block tridiagonal decomposition,
551/// proving that the **spectral gap is exactly λ₁ = 1**.
552///
553/// # Construction
554///
555/// 1. Decompose the Atlas into two hemispheres (by e₇)
556/// 2. Decompose each hemisphere into three Q₄ hypercube blocks (by d₄₅)
557/// 3. Verify block tridiagonal structure and identity matchings
558/// 4. Compute eigenvalues of 3×3 block matrices M_ν
559/// 5. Assemble the complete spectrum
560/// 6. Cross-check against trace identities and explicit Laplacian
561///
562/// # Example
563///
564/// ```
565/// use resonance_spectral_gap::{Atlas, SpectralAnalysis, Rational};
566///
567/// let atlas = Atlas::new();
568/// let spectral = SpectralAnalysis::from_atlas(&atlas);
569///
570/// // Main result: spectral gap is exactly 1
571/// assert_eq!(spectral.spectral_gap(), Rational::from_integer(1));
572///
573/// // 11 distinct eigenvalues, 48 counting multiplicity
574/// let total: usize = spectral.hemisphere_spectrum().iter().map(|(_, m)| m).sum();
575/// assert_eq!(total, 48);
576/// ```
577#[derive(Debug, Clone)]
578pub struct SpectralAnalysis {
579    hemisphere_spectrum: Vec<(Rational, usize)>,
580    full_spectrum: Vec<(Rational, usize)>,
581    spectral_gap: Rational,
582    hemisphere_trace: Rational,
583    hemisphere_trace_squared: Rational,
584    max_eigenvalue: Rational,
585    #[allow(dead_code)]
586    hemisphere_laplacian: DenseRatMatrix,
587}
588
589impl SpectralAnalysis {
590    /// Construct spectral analysis from an Atlas graph.
591    ///
592    /// Performs the full block tridiagonal decomposition, verifying all
593    /// structural assumptions and computing the complete spectrum with
594    /// trace cross-checks.
595    ///
596    /// # Panics
597    ///
598    /// Panics if any structural verification fails.
599    #[must_use]
600    pub fn from_atlas(atlas: &Atlas) -> Self {
601        assert_eq!(atlas.num_vertices(), 96);
602        assert_eq!(atlas.num_edges(), 256);
603
604        let hemispheres = decompose_hemispheres(atlas);
605        verify_no_cross_hemisphere_edges(atlas, &hemispheres);
606
607        let blocks = decompose_blocks(atlas, &hemispheres[0]);
608        for block in &blocks {
609            verify_q4_block(atlas, block);
610            verify_q4_spectrum(atlas, block);
611        }
612        verify_inter_block_structure(atlas, &blocks);
613
614        let hemisphere_spectrum = compute_hemisphere_spectrum();
615        let total_mult: usize = hemisphere_spectrum.iter().map(|(_, m)| *m).sum();
616        assert_eq!(total_mult, HEMISPHERE_VERTICES);
617
618        let spectral_gap = Self::find_spectral_gap(&hemisphere_spectrum);
619        let hemisphere_trace = Self::compute_trace(&hemisphere_spectrum);
620        let hemisphere_trace_squared = Self::compute_trace_squared(&hemisphere_spectrum);
621        let max_eigenvalue = Self::find_max_eigenvalue(&hemisphere_spectrum);
622
623        let hemisphere_laplacian = DenseRatMatrix::from_atlas_subgraph(atlas, &hemispheres[0]);
624
625        // Cross-check trace against explicit Laplacian
626        assert_eq!(hemisphere_laplacian.trace(), hemisphere_trace);
627        assert_eq!(hemisphere_laplacian.frobenius_squared(), hemisphere_trace_squared);
628        assert!(hemisphere_laplacian.is_symmetric());
629        for i in 0..HEMISPHERE_VERTICES {
630            assert!(hemisphere_laplacian.row_sum(i).is_zero());
631        }
632        assert!(hemisphere_laplacian.has_nonneg_diagonal());
633        assert!(hemisphere_laplacian.has_nonpos_off_diagonal());
634
635        let full_spectrum: Vec<(Rational, usize)> =
636            hemisphere_spectrum.iter().map(|&(ev, mult)| (ev, mult * 2)).collect();
637
638        Self {
639            hemisphere_spectrum,
640            full_spectrum,
641            spectral_gap,
642            hemisphere_trace,
643            hemisphere_trace_squared,
644            max_eigenvalue,
645            hemisphere_laplacian,
646        }
647    }
648
649    /// The spectral gap: smallest nonzero eigenvalue.
650    ///
651    /// **Theorem**: λ₁ = 1, arising from the M₀ block (ν = 0, eigenvalues {0, 1, 3}).
652    #[must_use]
653    pub const fn spectral_gap(&self) -> Rational {
654        self.spectral_gap
655    }
656
657    /// Complete hemisphere spectrum as `(eigenvalue, multiplicity)` pairs.
658    ///
659    /// Sorted by eigenvalue. Total multiplicity = 48.
660    ///
661    /// | λ  | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 11 |
662    /// |----|---|---|---|---|---|---|---|---|---|---|----|
663    /// | mult | 1 | 1 | 4 | 5 | 6 | 10| 4 | 10| 1 | 5 | 1  |
664    #[must_use]
665    pub fn hemisphere_spectrum(&self) -> &[(Rational, usize)] {
666        &self.hemisphere_spectrum
667    }
668
669    /// Complete full Atlas spectrum (double the hemisphere multiplicities).
670    ///
671    /// Total multiplicity = 96 (two disconnected hemispheres).
672    #[must_use]
673    pub fn full_spectrum(&self) -> &[(Rational, usize)] {
674        &self.full_spectrum
675    }
676
677    /// Hemisphere Laplacian trace: `tr(L_H) = 2|E_H| = 256`.
678    #[must_use]
679    pub const fn hemisphere_trace(&self) -> Rational {
680        self.hemisphere_trace
681    }
682
683    /// Hemisphere `tr(L_H²) = 1632`.
684    #[must_use]
685    pub const fn hemisphere_trace_squared(&self) -> Rational {
686        self.hemisphere_trace_squared
687    }
688
689    /// Maximum eigenvalue (11, from M₈ block).
690    #[must_use]
691    pub const fn max_eigenvalue(&self) -> Rational {
692        self.max_eigenvalue
693    }
694
695    /// Number of distinct eigenvalues (11).
696    #[must_use]
697    pub fn num_distinct_eigenvalues(&self) -> usize {
698        self.hemisphere_spectrum.len()
699    }
700
701    /// Gershgorin upper bound: `2 × max_degree = 12`.
702    #[must_use]
703    pub fn gershgorin_upper_bound(&self) -> Rational {
704        Rational::from_integer(12)
705    }
706
707    /// Whether all eigenvalues are integers.
708    #[must_use]
709    pub fn all_eigenvalues_integer(&self) -> bool {
710        self.hemisphere_spectrum.iter().all(|(ev, _)| *ev.denom() == 1)
711    }
712
713    /// Get the 3×3 block tridiagonal matrix M_ν for a Q₄ eigenvalue.
714    ///
715    /// # Panics
716    ///
717    /// Panics if ν is not in {0, 2, 4, 6, 8}.
718    #[must_use]
719    pub fn block_matrix(&self, nu: i64) -> [[Rational; 3]; 3] {
720        assert!(Q4_SPECTRUM.iter().any(|&(ev, _)| ev == nu), "ν = {nu} is not a Q4 eigenvalue");
721        block_tridiagonal_matrix(nu)
722    }
723
724    fn find_spectral_gap(spectrum: &[(Rational, usize)]) -> Rational {
725        spectrum.iter().find(|(ev, _)| !ev.is_zero()).map(|(ev, _)| *ev).expect("no nonzero eigenvalue")
726    }
727
728    #[allow(clippy::cast_possible_wrap)]
729    fn compute_trace(spectrum: &[(Rational, usize)]) -> Rational {
730        spectrum.iter().fold(Rational::zero(), |acc, &(ev, mult)| acc + ev * Rational::from_integer(mult as i64))
731    }
732
733    #[allow(clippy::cast_possible_wrap)]
734    fn compute_trace_squared(spectrum: &[(Rational, usize)]) -> Rational {
735        spectrum.iter().fold(Rational::zero(), |acc, &(ev, mult)| acc + ev * ev * Rational::from_integer(mult as i64))
736    }
737
738    fn find_max_eigenvalue(spectrum: &[(Rational, usize)]) -> Rational {
739        spectrum.last().map_or_else(Rational::zero, |&(ev, _)| ev)
740    }
741}
742
743impl fmt::Display for SpectralAnalysis {
744    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
745        writeln!(f, "Atlas Spectral Analysis")?;
746        writeln!(f, "======================")?;
747        writeln!(f)?;
748        writeln!(f, "Spectral gap: λ₁ = {}", self.spectral_gap)?;
749        writeln!(f, "Max eigenvalue: λ_max = {}", self.max_eigenvalue)?;
750        writeln!(f, "Distinct eigenvalues: {}", self.num_distinct_eigenvalues())?;
751        writeln!(f)?;
752        writeln!(f, "Hemisphere spectrum (48 eigenvalues):")?;
753        for &(ev, mult) in &self.hemisphere_spectrum {
754            writeln!(f, "  λ = {ev:>3}  multiplicity = {mult}")?;
755        }
756        writeln!(f)?;
757        writeln!(f, "Trace:  {}", self.hemisphere_trace)?;
758        writeln!(f, "Trace²: {}", self.hemisphere_trace_squared)?;
759        writeln!(f, "Gershgorin bound: [0, {}]", self.gershgorin_upper_bound())?;
760        Ok(())
761    }
762}
763
764// ─────────────────────────────────────────────────────────────────────────────
765// Unit Tests
766// ─────────────────────────────────────────────────────────────────────────────
767
768#[cfg(test)]
769mod tests {
770    use super::*;
771
772    #[test]
773    fn test_atlas_construction() {
774        let atlas = Atlas::new();
775        assert_eq!(atlas.num_vertices(), 96);
776        assert_eq!(atlas.num_edges(), 256);
777    }
778
779    #[test]
780    fn test_hemisphere_decomposition() {
781        let atlas = Atlas::new();
782        let hemispheres = decompose_hemispheres(&atlas);
783        assert_eq!(hemispheres[0].len(), 48);
784        assert_eq!(hemispheres[1].len(), 48);
785        verify_no_cross_hemisphere_edges(&atlas, &hemispheres);
786    }
787
788    #[test]
789    fn test_q4_blocks() {
790        let atlas = Atlas::new();
791        let hemispheres = decompose_hemispheres(&atlas);
792        let blocks = decompose_blocks(&atlas, &hemispheres[0]);
793        for block in &blocks {
794            assert_eq!(block.len(), 16);
795            verify_q4_block(&atlas, block);
796            verify_q4_spectrum(&atlas, block);
797        }
798        verify_inter_block_structure(&atlas, &blocks);
799    }
800
801    #[test]
802    fn test_block_eigenvalues() {
803        for &(nu, _) in &Q4_SPECTRUM {
804            let eigs = verify_block_eigenvalues(nu);
805            let nu_r = Rational::from_integer(nu);
806            assert_eq!(eigs[0], nu_r);
807            assert_eq!(eigs[1], nu_r + Rational::one());
808            assert_eq!(eigs[2], nu_r + Rational::new(3, 1));
809        }
810    }
811
812    #[test]
813    fn test_hemisphere_spectrum() {
814        let spectrum = compute_hemisphere_spectrum();
815        let expected: Vec<(Rational, usize)> = vec![
816            (Rational::zero(), 1),
817            (Rational::from_integer(1), 1),
818            (Rational::from_integer(2), 4),
819            (Rational::from_integer(3), 5),
820            (Rational::from_integer(4), 6),
821            (Rational::from_integer(5), 10),
822            (Rational::from_integer(6), 4),
823            (Rational::from_integer(7), 10),
824            (Rational::from_integer(8), 1),
825            (Rational::from_integer(9), 5),
826            (Rational::from_integer(11), 1),
827        ];
828        assert_eq!(spectrum, expected);
829    }
830
831    #[test]
832    fn test_spectral_gap_is_one() {
833        let atlas = Atlas::new();
834        let spectral = SpectralAnalysis::from_atlas(&atlas);
835        assert_eq!(spectral.spectral_gap(), Rational::from_integer(1));
836        assert_eq!(*spectral.spectral_gap().numer(), 1);
837        assert_eq!(*spectral.spectral_gap().denom(), 1);
838    }
839
840    #[test]
841    fn test_traces() {
842        let atlas = Atlas::new();
843        let spectral = SpectralAnalysis::from_atlas(&atlas);
844        assert_eq!(spectral.hemisphere_trace(), Rational::from_integer(256));
845        assert_eq!(spectral.hemisphere_trace_squared(), Rational::from_integer(1632));
846    }
847
848    #[test]
849    fn test_max_eigenvalue() {
850        let atlas = Atlas::new();
851        let spectral = SpectralAnalysis::from_atlas(&atlas);
852        assert_eq!(spectral.max_eigenvalue(), Rational::from_integer(11));
853    }
854
855    #[test]
856    fn test_all_eigenvalues_integer() {
857        let atlas = Atlas::new();
858        let spectral = SpectralAnalysis::from_atlas(&atlas);
859        assert!(spectral.all_eigenvalues_integer());
860    }
861
862    #[test]
863    fn test_full_spectrum() {
864        let atlas = Atlas::new();
865        let spectral = SpectralAnalysis::from_atlas(&atlas);
866        let total: usize = spectral.full_spectrum().iter().map(|(_, m)| m).sum();
867        assert_eq!(total, 96);
868    }
869
870    #[test]
871    fn test_gershgorin() {
872        let atlas = Atlas::new();
873        let spectral = SpectralAnalysis::from_atlas(&atlas);
874        let bound = spectral.gershgorin_upper_bound();
875        for &(ev, _) in spectral.hemisphere_spectrum() {
876            assert!(ev >= Rational::zero());
877            assert!(ev <= bound);
878        }
879    }
880
881    #[test]
882    fn test_display() {
883        let atlas = Atlas::new();
884        let spectral = SpectralAnalysis::from_atlas(&atlas);
885        let output = format!("{spectral}");
886        assert!(output.contains("Spectral gap"));
887        assert!(output.contains("λ₁ = 1"));
888    }
889}
resonance_spectral_gap/lib.rs

resonance_spectral_gap/
lib.rs
