tawasuyu/hifas

1//! Aritmética **exacta** sobre 𝔽₂ (§7.1), bit-packed con `u64`.
2//!
3//! La versión "topológica pura": sin epsilon. Reinterpreta las entradas de una
4//! `Matrix` `f64` módulo 2 (las restricciones del MVP viven en `{-1,0,1}`, y
5//! `-1 ≡ 1 (mód 2)`), y hace eliminación gaussiana sobre 𝔽₂ con XOR de filas.
6
7use crate::linalg::{LinAlg, Matrix, Scalar};
8
9/// Matriz densa sobre 𝔽₂: cada fila son `words_per_row` palabras de 64 bits.
10#[derive(Debug, Clone)]
11struct BitMatrix {
12    rows: usize,
13    cols: usize,
14    words_per_row: usize,
15    /// `rows * words_per_row` palabras, fila mayor.
16    data: Vec<u64>,
17}
18
19impl BitMatrix {
20    fn new(rows: usize, cols: usize) -> Self {
21        let words_per_row = cols.div_ceil(64);
22        Self {
23            rows,
24            cols,
25            words_per_row,
26            data: vec![0; rows * words_per_row],
27        }
28    }
29
30    /// Construye la matriz binaria reduciendo cada entrada `f64` módulo 2.
31    fn from_matrix(m: &Matrix) -> Self {
32        let mut b = Self::new(m.rows, m.cols);
33        for r in 0..m.rows {
34            for c in 0..m.cols {
35                if (m.get(r, c).round() as i64).rem_euclid(2) == 1 {
36                    b.set(r, c);
37                }
38            }
39        }
40        b
41    }
42
43    #[inline]
44    fn get(&self, r: usize, c: usize) -> bool {
45        let idx = r * self.words_per_row + (c >> 6);
46        (self.data[idx] >> (c & 63)) & 1 == 1
47    }
48
49    #[inline]
50    fn set(&mut self, r: usize, c: usize) {
51        let idx = r * self.words_per_row + (c >> 6);
52        self.data[idx] |= 1 << (c & 63);
53    }
54
55    fn swap_rows(&mut self, a: usize, b: usize) {
56        if a == b {
57            return;
58        }
59        for w in 0..self.words_per_row {
60            self.data.swap(a * self.words_per_row + w, b * self.words_per_row + w);
61        }
62    }
63
64    /// `fila[dst] ^= fila[src]` (suma sobre 𝔽₂).
65    fn xor_row_into(&mut self, dst: usize, src: usize) {
66        for w in 0..self.words_per_row {
67            self.data[dst * self.words_per_row + w] ^= self.data[src * self.words_per_row + w];
68        }
69    }
70
71    /// Lleva la matriz a RREF sobre 𝔽₂ in situ; devuelve las columnas pivote.
72    fn rref(&mut self) -> Vec<usize> {
73        let mut pivots = Vec::new();
74        let mut row = 0;
75        for col in 0..self.cols {
76            if row >= self.rows {
77                break;
78            }
79            // Busca una fila (>= row) con bit en esta columna.
80            let Some(sel) = (row..self.rows).find(|&r| self.get(r, col)) else {
81                continue; // columna libre
82            };
83            self.swap_rows(row, sel);
84            // Elimina la columna del resto de filas (arriba y abajo → RREF).
85            for r in 0..self.rows {
86                if r != row && self.get(r, col) {
87                    self.xor_row_into(r, row);
88                }
89            }
90            pivots.push(col);
91            row += 1;
92        }
93        pivots
94    }
95
96    /// Base del espacio nulo `{x : self·x = 0}` sobre 𝔽₂, una por columna libre.
97    fn null_space(&self) -> Vec<Vec<Scalar>> {
98        let mut a = self.clone();
99        let pivots = a.rref();
100        let mut basis = Vec::new();
101        for free in 0..a.cols {
102            if pivots.contains(&free) {
103                continue;
104            }
105            let mut vec = vec![0.0; a.cols];
106            vec[free] = 1.0;
107            // Sobre 𝔽₂, el coeficiente de despeje es +rref[i][free] (−x ≡ x).
108            for (i, &pc) in pivots.iter().enumerate() {
109                if a.get(i, free) {
110                    vec[pc] = 1.0;
111                }
112            }
113            basis.push(vec);
114        }
115        basis
116    }
117}
118
119/// Backend de álgebra lineal **exacto** sobre 𝔽₂ (§7.1).
120#[derive(Debug, Default, Clone, Copy)]
121pub struct Gf2Backend;
122
123impl LinAlg for Gf2Backend {
124    fn rank(&self, m: &Matrix) -> usize {
125        BitMatrix::from_matrix(m).rref().len()
126    }
127
128    fn cokernel_basis(&self, m: &Matrix) -> Vec<Vec<Scalar>> {
129        // coker(m) ≅ ker(mᵀ) sobre 𝔽₂: vectores del espacio de aristas.
130        BitMatrix::from_matrix(&m.transpose()).null_space()
131    }
132}
133
134#[cfg(test)]
135mod tests {
136    use super::*;
137
138    fn mat(rows: &[&[Scalar]]) -> Matrix {
139        let mut m = Matrix::zeros(rows.len(), rows[0].len());
140        for (i, row) in rows.iter().enumerate() {
141            for (j, &v) in row.iter().enumerate() {
142                m.set(i, j, v);
143            }
144        }
145        m
146    }
147
148    #[test]
149    fn rank_exacto_sobre_gf2() {
150        // Sobre 𝔽₂ la tercera fila = fila0 XOR fila1, luego rango 2.
151        let m = mat(&[&[1.0, 1.0, 0.0], &[0.0, 1.0, 1.0], &[1.0, 0.0, 1.0]]);
152        assert_eq!(Gf2Backend.rank(&m), 2);
153    }
154
155    #[test]
156    fn menos_uno_es_uno_mod_dos() {
157        // Fila [1,-1] y [−1,1] son idénticas mód 2 → rango 1.
158        let m = mat(&[&[1.0, -1.0], &[-1.0, 1.0]]);
159        assert_eq!(Gf2Backend.rank(&m), 1);
160    }
161
162    #[test]
163    fn cokernel_del_triangulo() {
164        // Incidencia (no signada) del triángulo sobre 𝔽₂: cokernel dim 1 = [1,1,1].
165        let delta = mat(&[&[1.0, 1.0, 0.0], &[0.0, 1.0, 1.0], &[1.0, 0.0, 1.0]]);
166        let basis = Gf2Backend.cokernel_basis(&delta);
167        assert_eq!(basis.len(), 1);
168        assert!(basis[0].iter().all(|&x| x == 1.0), "el ciclo toca las 3 aristas");
169    }
170}