breen deligne homology

main.rs

// [package]
// name = "breen-deligne-homology"
// version = "0.1.0"
// authors = ["Mario Carneiro <[email protected]>"]
// edition = "2018"

// [dependencies]
// modinverse = "0.1"
// rand = "0.8"

use std::{rc::Rc, hash::Hash};
use std::collections::{HashMap, hash_map::Entry};

#[derive(Copy, Clone, Hash, PartialEq, Eq)]
struct Fp(bool);
const B: u8 = 2; // This version has been specialized to B = 2. See earlier revisions
const I: u8 = 4;

const RANDOM: bool = true;

impl std::ops::Add<Fp> for Fp {
  type Output = Fp;
  fn add(self, rhs: Fp) -> Fp {
    Fp(self.0 ^ rhs.0)
  }
}
impl std::ops::AddAssign<Fp> for Fp {
  fn add_assign(&mut self, rhs: Fp) { *self = *self + rhs }
}
impl std::ops::Sub<Fp> for Fp {
  type Output = Fp;
  fn sub(self, rhs: Fp) -> Fp {
    Fp(self.0 ^ rhs.0)
  }
}
impl std::ops::Neg for Fp {
  type Output = Fp;
  fn neg(self) -> Fp { self }
}
impl std::ops::SubAssign<Fp> for Fp {
  fn sub_assign(&mut self, rhs: Fp) { *self = *self - rhs }
}
impl std::ops::Mul<Fp> for Fp {
  type Output = Fp;
  fn mul(self, rhs: Fp) -> Fp {
    Fp(self.0 & rhs.0)
  }
}
impl std::ops::MulAssign<Fp> for Fp {
  fn mul_assign(&mut self, rhs: Fp) { *self = *self * rhs }
}

impl Fp {
  fn inv(self) -> Option<Fp> {
    if self.0 { Some(self) } else { None }
    // modinverse::modinverse(self.0 as i16, B as i16).map(|i| Fp(i as _))
  }
}

fn sel(m: &[bool], t: &[Fp]) -> Box<[Fp]> {
  m.iter().zip(t).filter(|p| *p.0).map(|p| *p.1).collect()
}

fn iadd(mut a: Box<[Fp]>, b: &[Fp]) -> Box<[Fp]> {
  a.iter_mut().zip(b).for_each(|(a, b)| *a += *b);
  a
}

#[derive(Clone, Hash, PartialEq, Eq)]
enum NaiveMat {
  One(Box<[bool]>),
  Two(Box<[bool]>, Box<[bool]>),
}

impl NaiveMat {
  fn double(&self) -> Self {
    let f = |m: &Box<[bool]>| m.iter().chain(&**m).copied().collect();
    match self {
      NaiveMat::One(m) => NaiveMat::One(f(m)),
      NaiveMat::Two(m1, m2) => NaiveMat::Two(f(m1), f(m2)),
    }
  }
}

#[derive(Clone)]
struct FreeAb<K=Rc<[Fp]>>(HashMap<K, i64>);

impl<K> Default for FreeAb<K> {
  fn default() -> Self { Self(HashMap::new()) }
}
impl<K: Hash + Eq> FreeAb<K> {
  fn single(t: impl Into<K>) -> Self {
    let mut m = HashMap::new();
    m.insert(t.into(), 1);
    Self(m)
  }

  fn add_kv(&mut self, k: K, v: i64) {
    match self.0.entry(k) {
      Entry::Occupied(mut e) => *e.get_mut() += v,
      Entry::Vacant(e) => {e.insert(v);}
    }
  }

  fn dmap(&self, f: impl Fn(&K) -> K) -> Self {
    let mut m = Self(HashMap::new());
    for (k, &v) in &self.0 { m.add_kv(f(k), v) }
    m
  }

  fn add_smul(&mut self, n: i64, rhs: &FreeAb<K>) where K: Clone {
    for (k, &v) in &rhs.0 { self.add_kv(k.clone(), n * v) }
  }
}

impl FreeAb {
  fn ev_nv(&self, nm: &NaiveMat) -> Self {
    match nm {
      NaiveMat::One(m) => self.dmap(|k| sel(m, k).into()),
      NaiveMat::Two(m1, m2) => self.dmap(|k| iadd(sel(m1, k), &sel(m2, k)).into()),
    }
  }

  fn ev(&self, nms: &FreeAb<NaiveMat>) -> Self {
    let mut r = Self::default();
    for (nm, &n) in &nms.0 {
      r.add_smul(n, &self.ev_nv(nm));
    }
    r
  }
}

impl FreeAb<NaiveMat> {
  fn double(&self) -> Self { self.dmap(NaiveMat::double) }
}

impl<'a, K: Hash + Eq + Clone> std::ops::AddAssign<&'a FreeAb<K>> for FreeAb<K> {
  fn add_assign(&mut self, rhs: &'a FreeAb<K>) {
    for (k, &v) in &rhs.0 { self.add_kv(k.clone(), v) }
  }
}

impl<'a, K: Hash + Eq + Clone> std::ops::Add<&'a FreeAb<K>> for FreeAb<K> {
  type Output = Self;
  fn add(mut self, rhs: &'a FreeAb<K>) -> Self { self += rhs; self }
}

impl<'a, K: Hash + Eq + Clone> std::ops::SubAssign<&'a FreeAb<K>> for FreeAb<K> {
  fn sub_assign(&mut self, rhs: &'a FreeAb<K>) {
    for (k, &v) in &rhs.0 { self.add_kv(k.clone(), -v) }
  }
}
impl<'a, K: Hash + Eq + Clone> std::ops::Sub<&'a FreeAb<K>> for FreeAb<K> {
  type Output = Self;
  fn sub(mut self, rhs: &'a FreeAb<K>) -> Self { self -= rhs; self }
}

impl<'a, K> std::ops::MulAssign<i64> for FreeAb<K> {
  fn mul_assign(&mut self, rhs: i64) {
    for v in self.0.values_mut() { *v *= rhs }
  }
}

const D: usize = 1 << I;
const W: usize = (B as usize).pow(D as u32);
const S: usize = 1 << (I + 1);

fn encode(v: &[Fp]) -> usize {
  v.iter().fold(0, |r, &x| r * B as usize + x.0 as usize)
}

fn enc_ev(d_n: &FreeAb<NaiveMat>, t: &[Fp; S], r: &mut [Fp; W]) {
  let dt = FreeAb::<Rc<[Fp]>>::single(t.to_vec()).ev(&d_n);
  for (k, v) in dt.0 {
    let i = encode(&k);
    r[i] += Fp(v & 1 != 0);
  }
}

fn d(n: u8) -> FreeAb<NaiveMat> {
  let m = 1 << n;
  let mut pi1_n = vec![false; 2 * m].into_boxed_slice();
  let mut pi2_n = vec![false; 2 * m].into_boxed_slice();
  for i in 0..m { pi1_n[i] = true; pi2_n[m + i] = true; }
  let pi_n =
    FreeAb::<NaiveMat>::single(NaiveMat::One(pi1_n.clone())) +
    &FreeAb::single(NaiveMat::One(pi2_n.clone()));
  let sigma_n = FreeAb::single(NaiveMat::Two(pi1_n, pi2_n));
  let r = pi_n - &sigma_n;
  if n == 0 { r } else { r - &d(n-1).double() }
}

struct BasisElement {
  pivot: usize,
  values: Box<[Fp; W]>,
}

fn gaussian_eliminate_one(basis: &mut Vec<BasisElement>, r: &mut [Fp; W]) -> bool {
  let mut next = 0;
  while let Some(i) = r[next..].iter().position(|&x| x != Fp(false)).map(|a| a+next) {
    match basis.binary_search_by_key(&i, |b| b.pivot) {
      Ok(a) => {
        let b = &basis[a].values;
        for j in i..W { r[j] -= b[j] }
        next = i;
      }
      Err(a) => {
        let z = r[i].inv().unwrap();
        r[i] = Fp(true);
        for j in i+1..W { r[j] *= z }
        // for BasisElement {values, ..} in &mut basis[..a] {
        //   let x = -values[i];
        //   if x != Fp(0) {
        //     for j in i..W { values[j] += x * r[j] }
        //   }
        // }
        basis.insert(a, BasisElement {pivot: i, values: Box::new(r.clone())});
        return true
      }
    }
  }
  false
}

fn run_one(d_n: &FreeAb<NaiveMat>, t: &[Fp; S], basis: &mut Vec<BasisElement>, r: &mut [Fp; W]) {
  enc_ev(d_n, t, r);
  if gaussian_eliminate_one(basis, r) {
    *r = [Fp(false); W];
    println!("{:?}", basis.len());
  }
}

const MAX: usize = (B as usize).pow(S as u32);

fn exhaustive(d_n: &FreeAb<NaiveMat>, basis: &mut Vec<BasisElement>) {
  let mut iteration = 0;
  let mut per_10000 = 0;
  let mut t = [Fp(false); S];
  let mut r = [Fp(false); W];
  'a: loop {
    run_one(d_n, &t, basis, &mut r);
    iteration += 1;
    if iteration == MAX / 10000 {
      iteration = 0;
      per_10000 += 1;
      println!("==>      {}/10000", per_10000);
    }
    for i in t.iter_mut() {
      i.0 = !i.0;
      if !i.0 { continue 'a }
    }
    break
  }
}

fn random(d_n: &FreeAb<NaiveMat>, basis: &mut Vec<BasisElement>) {
  let between = rand::distributions::Uniform::from(0..B);
  let mut rng = rand::thread_rng();
  let mut iteration = 0;
  let mut per_10000 = 0;
  let mut t = [Fp(false); S];
  let mut r = [Fp(false); W];
  loop {
    for i in &mut t { i.0 = rand::distributions::Distribution::sample(&between, &mut rng) != 0 }
    run_one(d_n, &t, basis, &mut r);
    iteration += 1;
    if iteration == MAX / 10000 {
      iteration = 0;
      per_10000 += 1;
      println!("==>      {}/10000", per_10000);
    }
  }
}

fn main() {
  let d_n = d(I);
  let mut basis: Vec<BasisElement> = vec![];
  if RANDOM {
    random(&d_n, &mut basis);
  } else {
    exhaustive(&d_n, &mut basis);
  }
  println!("{:?}", basis.len());
}