//! Aritmética **exacta** sobre 𝔽₂ (§7.1), bit-packed con `u64`.
//!
//! La versión "topológica pura": sin epsilon. Reinterpreta las entradas de una
//! `Matrix` `f64` módulo 2 (las restricciones del MVP viven en `{-1,0,1}`, y
//! `-1 ≡ 1 (mód 2)`), y hace eliminación gaussiana sobre 𝔽₂ con XOR de filas.

use crate::linalg::{LinAlg, Matrix, Scalar};

/// Matriz densa sobre 𝔽₂: cada fila son `words_per_row` palabras de 64 bits.
#[derive(Debug, Clone)]
struct BitMatrix {
    rows: usize,
    cols: usize,
    words_per_row: usize,
    /// `rows * words_per_row` palabras, fila mayor.
    data: Vec<u64>,
}

impl BitMatrix {
    fn new(rows: usize, cols: usize) -> Self {
        let words_per_row = cols.div_ceil(64);
        Self {
            rows,
            cols,
            words_per_row,
            data: vec![0; rows * words_per_row],
        }
    }

    /// Construye la matriz binaria reduciendo cada entrada `f64` módulo 2.
    fn from_matrix(m: &Matrix) -> Self {
        let mut b = Self::new(m.rows, m.cols);
        for r in 0..m.rows {
            for c in 0..m.cols {
                if (m.get(r, c).round() as i64).rem_euclid(2) == 1 {
                    b.set(r, c);
                }
            }
        }
        b
    }

    #[inline]
    fn get(&self, r: usize, c: usize) -> bool {
        let idx = r * self.words_per_row + (c >> 6);
        (self.data[idx] >> (c & 63)) & 1 == 1
    }

    #[inline]
    fn set(&mut self, r: usize, c: usize) {
        let idx = r * self.words_per_row + (c >> 6);
        self.data[idx] |= 1 << (c & 63);
    }

    fn swap_rows(&mut self, a: usize, b: usize) {
        if a == b {
            return;
        }
        for w in 0..self.words_per_row {
            self.data.swap(a * self.words_per_row + w, b * self.words_per_row + w);
        }
    }

    /// `fila[dst] ^= fila[src]` (suma sobre 𝔽₂).
    fn xor_row_into(&mut self, dst: usize, src: usize) {
        for w in 0..self.words_per_row {
            self.data[dst * self.words_per_row + w] ^= self.data[src * self.words_per_row + w];
        }
    }

    /// Lleva la matriz a RREF sobre 𝔽₂ in situ; devuelve las columnas pivote.
    fn rref(&mut self) -> Vec<usize> {
        let mut pivots = Vec::new();
        let mut row = 0;
        for col in 0..self.cols {
            if row >= self.rows {
                break;
            }
            // Busca una fila (>= row) con bit en esta columna.
            let Some(sel) = (row..self.rows).find(|&r| self.get(r, col)) else {
                continue; // columna libre
            };
            self.swap_rows(row, sel);
            // Elimina la columna del resto de filas (arriba y abajo → RREF).
            for r in 0..self.rows {
                if r != row && self.get(r, col) {
                    self.xor_row_into(r, row);
                }
            }
            pivots.push(col);
            row += 1;
        }
        pivots
    }

    /// Base del espacio nulo `{x : self·x = 0}` sobre 𝔽₂, una por columna libre.
    fn null_space(&self) -> Vec<Vec<Scalar>> {
        let mut a = self.clone();
        let pivots = a.rref();
        let mut basis = Vec::new();
        for free in 0..a.cols {
            if pivots.contains(&free) {
                continue;
            }
            let mut vec = vec![0.0; a.cols];
            vec[free] = 1.0;
            // Sobre 𝔽₂, el coeficiente de despeje es +rref[i][free] (−x ≡ x).
            for (i, &pc) in pivots.iter().enumerate() {
                if a.get(i, free) {
                    vec[pc] = 1.0;
                }
            }
            basis.push(vec);
        }
        basis
    }
}

/// Backend de álgebra lineal **exacto** sobre 𝔽₂ (§7.1).
#[derive(Debug, Default, Clone, Copy)]
pub struct Gf2Backend;

impl LinAlg for Gf2Backend {
    fn rank(&self, m: &Matrix) -> usize {
        BitMatrix::from_matrix(m).rref().len()
    }

    fn cokernel_basis(&self, m: &Matrix) -> Vec<Vec<Scalar>> {
        // coker(m) ≅ ker(mᵀ) sobre 𝔽₂: vectores del espacio de aristas.
        BitMatrix::from_matrix(&m.transpose()).null_space()
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    fn mat(rows: &[&[Scalar]]) -> Matrix {
        let mut m = Matrix::zeros(rows.len(), rows[0].len());
        for (i, row) in rows.iter().enumerate() {
            for (j, &v) in row.iter().enumerate() {
                m.set(i, j, v);
            }
        }
        m
    }

    #[test]
    fn rank_exacto_sobre_gf2() {
        // Sobre 𝔽₂ la tercera fila = fila0 XOR fila1, luego rango 2.
        let m = mat(&[&[1.0, 1.0, 0.0], &[0.0, 1.0, 1.0], &[1.0, 0.0, 1.0]]);
        assert_eq!(Gf2Backend.rank(&m), 2);
    }

    #[test]
    fn menos_uno_es_uno_mod_dos() {
        // Fila [1,-1] y [−1,1] son idénticas mód 2 → rango 1.
        let m = mat(&[&[1.0, -1.0], &[-1.0, 1.0]]);
        assert_eq!(Gf2Backend.rank(&m), 1);
    }

    #[test]
    fn cokernel_del_triangulo() {
        // Incidencia (no signada) del triángulo sobre 𝔽₂: cokernel dim 1 = [1,1,1].
        let delta = mat(&[&[1.0, 1.0, 0.0], &[0.0, 1.0, 1.0], &[1.0, 0.0, 1.0]]);
        let basis = Gf2Backend.cokernel_basis(&delta);
        assert_eq!(basis.len(), 1);
        assert!(basis[0].iter().all(|&x| x == 1.0), "el ciclo toca las 3 aristas");
    }
}