29426297655

Committed 15 Jul 2026 03:02PM UTC coverage: 74.502%. First build

Build # 29426297655

Build Type

Pull #111

github

Committed by

web-flow

Commit Message

Merge 94de354d0 into 4d383f6b9

Pull Request Pull Request #111: Harden succinct structures and construction contracts

Coverage Stats

850 of 1056 new or added lines in 40 files covered. (80.49%)

9046 of 12142 relevant lines covered (74.5%)

18306116.8 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

93.94

/src/utils/mod2_sys.rs

/*
 *
 * SPDX-FileCopyrightText: 2025 Dario Moschetti
 * SPDX-FileCopyrightText: 2025 Sebastiano Vigna
 *
 * SPDX-License-Identifier: Apache-2.0 OR LGPL-2.1-or-later
 */

//! Fast solution of random systems of linear equations with coefficients in
//! **F**₂.

#![allow(unexpected_cfgs)]
#![allow(clippy::comparison_chain)]
use std::ptr;

use crate::{
    bit_vec,
    traits::Word,
    traits::{BitVecOps, BitVecOpsMut},
};
use anyhow::{Result, bail, ensure};
use arbitrary_chunks::ArbitraryChunks;

/// An equation on `W::BITS`-dimensional vectors with coefficients in **F**₂.
#[derive(Clone, Debug)]
pub struct Modulo2Equation<W: Word = usize> {
    /// The variables in increasing order.
    vars: Vec<u32>,
    /// The constant term
    c: W,
}

/// A system of [equations].
///
/// [equations]: struct.Modulo2Equation.html
#[derive(Clone, Debug)]
pub struct Modulo2System<W: Word = usize> {
    /// The number of variables.
    num_vars: usize,
    /// The equations in the system.
    equations: Vec<Modulo2Equation<W>>,
}

impl<W: Word> Modulo2Equation<W> {
    /// Creates a new equation with given variables and constant term.
    ///
    /// # Safety
    ///
    /// `vars` must be strictly increasing. In particular, duplicate variables
    /// are invalid: repeated coefficients cancel in F2 and must be removed by
    /// the caller before construction.
    pub unsafe fn from_parts(vars: Vec<u32>, c: W) -> Self {
        debug_assert!(vars.windows(2).all(|pair| pair[0] < pair[1]));
        Modulo2Equation { vars, c }
    }

    #[inline(always)]
    unsafe fn add_ptr(
        mut left: *const u32,
        left_end: *const u32,
        mut right: *const u32,
        right_end: *const u32,
        mut dst: *mut u32,
    ) -> *mut u32 {
        unsafe {
            while left != left_end && right != right_end {
                let less = *left <= *right;
                let more = *left >= *right;

                let src = if less { left } else { right };
                *dst = *src;

                left = left.add(less as usize);
                right = right.add(more as usize);
                dst = dst.add((less ^ more) as usize);
            }

            let rem_left = left_end.offset_from(left) as usize;
            ptr::copy_nonoverlapping(left, dst, rem_left);
            dst = dst.add(rem_left);
            let rem_right = right_end.offset_from(right) as usize;
            ptr::copy_nonoverlapping(right, dst, rem_right);
            dst = dst.add(rem_right);
            dst
        }
    }

    /// Adds `other` into `self`, using `scratch` as a working buffer to avoid
    /// per-call allocation.
    ///
    /// After the call, `scratch` holds what used to be `self.vars` (with its
    /// original capacity intact), ready to be reused on the next call. Over
    /// many calls, `scratch`'s capacity grows to the peak merge size and
    /// further calls allocate nothing.
    pub(crate) fn add_into(&mut self, other: &Modulo2Equation<W>, scratch: &mut Vec<u32>) {
        scratch.clear();
        scratch.reserve(self.vars.len() + other.vars.len());
        let dst = scratch.as_mut_ptr();

        let left_range = self.vars.as_ptr_range();
        let right_range = other.vars.as_ptr_range();

        unsafe {
            let end = Self::add_ptr(
                left_range.start,
                left_range.end,
                right_range.start,
                right_range.end,
                dst,
            );
            scratch.set_len(end.offset_from(dst) as usize);
        }

        self.c ^= other.c;
        std::mem::swap(&mut self.vars, scratch);
    }

    pub fn add(&mut self, other: &Modulo2Equation<W>) {
        let mut scratch = Vec::new();
        self.add_into(other, &mut scratch);
    }

    /// Checks whether the equation is unsolvable.
    fn is_unsolvable(&self) -> bool {
        self.vars.is_empty() && self.c != W::ZERO
    }

    /// Checks whether the equation is an identity.
    fn is_identity(&self) -> bool {
        self.vars.is_empty() && self.c == W::ZERO
    }

    /// Evaluates the XOR of the values associated to the given variables.
    fn eval_vars(vars: impl AsRef<[u32]>, values: impl AsRef<[W]>) -> W {
        let mut sum = W::ZERO;
        for &var in vars.as_ref() {
            sum ^= values.as_ref()[var as usize];
        }
        sum
    }
}

impl<W: Word> Modulo2System<W> {
    pub fn new(num_vars: usize) -> Self {
        Modulo2System {
            num_vars,
            equations: vec![],
        }
    }

    pub const fn num_vars(&self) -> usize {
        self.num_vars
    }

    pub fn num_equations(&self) -> usize {
        self.equations.len()
    }

    /// Creates a new `Modulo2System` from existing equations.
    ///
    /// # Safety
    ///
    /// Every variable index in every equation must be strictly less than
    /// `num_vars`. Equation variable lists must also satisfy
    /// [`Modulo2Equation::from_parts`]'s strictly-increasing invariant.
    pub unsafe fn from_parts(num_vars: usize, equations: Vec<Modulo2Equation<W>>) -> Self {
        Modulo2System {
            num_vars,
            equations,
        }
    }

    /// Adds an equation to the system.
    ///
    /// # Panics
    ///
    /// Panics if the equation contains a variable outside this system.
    pub fn push(&mut self, equation: Modulo2Equation<W>) {
        if let Some(&last_var) = equation.vars.last() {
            let last_var = usize::try_from(last_var).expect("variable index does not fit in usize");
            assert!(
                last_var < self.num_vars,
                "variable {last_var} is out of bounds for a system with {} variables",
                self.num_vars
            );
        }
        self.equations.push(equation);
    }

    /// Checks if a given solution satisfies the system of equations.
    pub fn check(&self, solution: &[W]) -> bool {
        assert_eq!(
            solution.len(),
            self.num_vars,
            "The number of variables in the solution ({}) does not match the number of variables in the system ({})",
            solution.len(),
            self.num_vars
        );
        self.equations
            .iter()
            .all(|eq| eq.c == Modulo2Equation::<W>::eval_vars(&eq.vars, solution))
    }

    /// Puts the system into echelon form.
    fn echelon_form(&mut self) -> Result<()> {
        let equations = &mut self.equations;
        if equations.is_empty() {
            return Ok(());
        }
        let mut scratch: Vec<u32> = Vec::new();
        'main: for i in 0..equations.len() - 1 {
            // A pivot row reduced to [] = c: unsolvable if c != 0, otherwise
            // a redundant identity that takes no part in elimination.
            if equations[i].vars.is_empty() {
                ensure!(!equations[i].is_unsolvable(), "System is unsolvable");
                continue;
            }
            for j in i + 1..equations.len() {
                let (left, right) = equations.split_at_mut(j);
                let eq_i = &mut left[i];
                let eq_j = &right[0];

                // A later row reduced to [] = c: unsolvable if c != 0, else a
                // redundant identity to skip. Without this, eq_j.vars[0] below
                // would panic on the empty variable list.
                if eq_j.vars.is_empty() {
                    ensure!(!eq_j.is_unsolvable(), "System is unsolvable");
                    continue;
                }

                let first_var_j = eq_j.vars[0];

                if eq_i.vars[0] == first_var_j {
                    eq_i.add_into(eq_j, &mut scratch);
                    if eq_i.is_unsolvable() {
                        bail!("System is unsolvable");
                    }
                    if eq_i.is_identity() {
                        continue 'main;
                    }
                }

                if eq_i.vars[0] > first_var_j {
                    equations.swap(i, j)
                }
            }
        }
        Ok(())
    }

    /// Solves the system using Gaussian elimination.
    pub fn gaussian_elimination(&mut self) -> Result<Vec<W>> {
        self.echelon_form()?;
        // An equation reduced to [] = c with c != 0 has no solution; the
        // back-substitution below would index eq.vars[0] on an empty variable
        // list. echelon_form only checks emptiness for rows before the last.
        for eq in &self.equations {
            ensure!(!eq.is_unsolvable(), "System is unsolvable");
        }
        let mut solution = vec![W::ZERO; self.num_vars];
        self.equations
            .iter()
            .rev()
            .filter(|eq| !eq.is_identity())
            .for_each(|eq| {
                solution[eq.vars[0] as usize] =
                    eq.c ^ Modulo2Equation::<W>::eval_vars(&eq.vars, &solution);
            });
        Ok(solution)
    }

    /// Builds the data structures needed for [lazy Gaussian elimination].
    ///
    /// This method returns the variable-to-equation mapping, the weight of each
    /// variable (the number of equations in which it appears), and the priority
    /// of each equation (the number of variables in it). The
    /// variable-to-equation mapping is materialized as a vector of mutable
    /// slices, each of which points inside a provided backing vector. This
    /// approach greatly reduces the number of allocations.
    ///
    /// [lazy Gaussian elimination]: https://doi.org/10.1016/j.ic.2020.104517
    fn setup<'a>(
        &self,
        backing: &'a mut Vec<usize>,
    ) -> (Vec<&'a mut [usize]>, Vec<usize>, Vec<usize>) {
        let mut weight = vec![0; self.num_vars];
        let mut priority = vec![0; self.equations.len()];

        let mut total_vars = 0;
        for (i, equation) in self.equations.iter().enumerate() {
            priority[i] = equation.vars.len();
            total_vars += equation.vars.len();
            for &var in &equation.vars {
                weight[var as usize] += 1;
            }
        }

        backing.resize(total_vars, 0);
        let mut var_to_eq: Vec<&mut [usize]> = Vec::with_capacity(self.num_vars);

        backing.arbitrary_chunks_mut(&weight).for_each(|chunk| {
            var_to_eq.push(chunk);
        });

        let mut pos = vec![0usize; self.num_vars];
        for (i, equation) in self.equations.iter().enumerate() {
            for &var in &equation.vars {
                let var = var as usize;
                var_to_eq[var][pos[var]] = i;
                pos[var] += 1;
            }
        }

        (var_to_eq, weight, priority)
    }

    /// Solves the system using [lazy Gaussian elimination].
    ///
    /// [lazy Gaussian elimination]: https://doi.org/10.1016/j.ic.2020.104517
    pub fn lazy_gaussian_elimination(&mut self) -> Result<Vec<W>> {
        let num_vars = self.num_vars;
        let num_equations = self.equations.len();

        if num_equations == 0 {
            return Ok(vec![W::ZERO; num_vars]);
        }

        let mut backing = vec![];
        let (var_to_eqs, mut weight, mut priority);
        (var_to_eqs, weight, priority) = self.setup(&mut backing);

        let mut variables = vec![0; num_vars];
        {
            let mut count = vec![0; num_equations + 1];

            for x in 0..num_vars {
                count[weight[x]] += 1
            }
            for i in 1..count.len() {
                count[i] += count[i - 1]
            }
            for i in (0..num_vars).rev() {
                count[weight[i]] -= 1;
                variables[count[weight[i]]] = i;
            }
        }

        let mut equation_list: Vec<usize> = (0..priority.len())
            .rev()
            .filter(|&x| priority[x] <= 1)
            .collect();

        let mut dense = vec![];
        let mut solved = vec![];
        let mut pivots = vec![];

        let equations = &mut self.equations;
        let mut idle = bit_vec![true; num_vars];
        let mut scratch: Vec<u32> = Vec::new();

        let mut remaining = equations.len();

        while remaining != 0 {
            if equation_list.is_empty() {
                let mut var = variables.pop().unwrap();
                while weight[var] == 0 {
                    var = variables.pop().unwrap()
                }
                idle.set(var, false);
                var_to_eqs[var].as_ref().iter().for_each(|&eq| {
                    priority[eq] -= 1;
                    if priority[eq] == 1 {
                        equation_list.push(eq)
                    }
                });
            } else {
                remaining -= 1;
                let first = equation_list.pop().unwrap();

                if priority[first] == 0 {
                    let equation = &mut equations[first];
                    if equation.is_unsolvable() {
                        bail!("System is unsolvable")
                    }
                    if equation.is_identity() {
                        continue;
                    }
                    dense.push(equation.to_owned());
                } else if priority[first] == 1 {
                    let pivot = equations[first]
                        .vars
                        .iter()
                        .copied()
                        .find(|x| idle.get(*x as usize))
                        .expect("Missing expected idle variable in equation");
                    pivots.push(pivot as usize);
                    solved.push(first);
                    weight[pivot as usize] = 0;
                    for &eq in var_to_eqs[pivot as usize]
                        .as_ref()
                        .iter()
                        .filter(|&&eq_idx| eq_idx != first)
                    {
                        let lo = eq.min(first);
                        let hi = eq.max(first);
                        let (left, right) = equations.split_at_mut(hi);
                        if eq < first {
                            left[lo].add_into(&right[0], &mut scratch);
                        } else {
                            right[0].add_into(&left[lo], &mut scratch);
                        }

                        priority[eq] -= 1;
                        if priority[eq] == 1 {
                            equation_list.push(eq)
                        }
                    }
                }
            }
        }

        // SAFETY: the usage of the method is safe, as the equations have the
        // right number of variables
        let mut dense_system = unsafe { Modulo2System::from_parts(num_vars, dense) };
        let mut solution = dense_system.gaussian_elimination()?;

        for i in 0..solved.len() {
            let eq = &equations[solved[i]];
            let pivot = pivots[i];
            assert!(solution[pivot] == W::ZERO);
            solution[pivot] = eq.c ^ Modulo2Equation::<W>::eval_vars(&eq.vars, &solution);
        }

        Ok(solution)
    }
}

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

    #[test]
    fn test_add_ptr() {
        let a = [0, 1, 3, 4];
        let b = [0, 2, 4, 5];
        let mut c = Vec::with_capacity(a.len() + b.len());

        let ra = a.as_ptr_range();
        let rb = b.as_ptr_range();
        let mut dst = c.as_mut_ptr();
        unsafe {
            dst = Modulo2Equation::<u32>::add_ptr(ra.start, ra.end, rb.start, rb.end, dst);
            assert_eq!(dst.offset_from(c.as_ptr()), 4);
            c.set_len(4);
            assert_eq!(c, vec![1, 2, 3, 5]);
        }
    }

    #[test]
    fn test_equation_builder() {
        let eq = unsafe { Modulo2Equation::from_parts(vec![0, 1, 2], 3_usize) };
        assert_eq!(eq.vars.len(), 3);
        assert_eq!(eq.vars.to_vec(), vec![0, 1, 2]);
    }

    #[test]
    fn test_equation_addition() {
        let mut eq0 = unsafe { Modulo2Equation::from_parts(vec![1, 4, 9], 3_usize) };
        let eq1 = unsafe { Modulo2Equation::from_parts(vec![1, 4, 10], 3_usize) };
        eq0.add(&eq1);
        assert_eq!(eq0.vars, vec![9, 10]);
        assert_eq!(eq0.c, 0);
    }

    #[test]
    fn test_system_one_equation() -> anyhow::Result<()> {
        let mut system = Modulo2System::<usize>::new(2);
        let eq = unsafe { Modulo2Equation::from_parts(vec![0], 3_usize) };
        system.push(eq);
        let solution = system.lazy_gaussian_elimination()?;
        assert!(system.check(&solution));
        Ok(())
    }

    #[test]
    #[should_panic(expected = "out of bounds")]
    fn test_push_rejects_out_of_range_variable() {
        let mut system = Modulo2System::<usize>::new(1);
        // SAFETY: a one-element variable list is strictly increasing.
        let equation = unsafe { Modulo2Equation::from_parts(vec![1], 0_usize) };
        system.push(equation);
    }

    #[test]
    fn test_impossible_system() {
        let mut system = Modulo2System::<usize>::new(1);
        let eq0 = unsafe { Modulo2Equation::from_parts(vec![0], 0_usize) };
        system.push(eq0);
        let eq1 = unsafe { Modulo2Equation::from_parts(vec![0], 1_usize) };
        system.push(eq1);
        let solution = system.lazy_gaussian_elimination();
        assert!(solution.is_err());
    }

    #[test]
    fn test_gaussian_elimination_last_unsolvable_equation() {
        let mut system = Modulo2System::<usize>::new(1);
        // SAFETY: variable 0 is in bounds for this one-variable system, and
        // the variable list is sorted.
        let eq0 = unsafe { Modulo2Equation::from_parts(vec![0], 0_usize) };
        system.push(eq0);
        // SAFETY: the empty variable list is sorted and contains no
        // out-of-bounds variable indices.
        let eq1 = unsafe { Modulo2Equation::from_parts(vec![], 1_usize) };
        system.push(eq1);

        let solution = system.gaussian_elimination();
        assert!(solution.is_err());
    }

    #[test]
    fn test_gaussian_elimination_leading_identity_row() {
        let mut system = Modulo2System::<usize>::new(1);
        // SAFETY: the empty variable list is sorted and contains no
        // out-of-bounds variable indices.
        let eq0 = unsafe { Modulo2Equation::from_parts(vec![], 0_usize) };
        system.push(eq0);
        // SAFETY: variable 0 is in bounds for this one-variable system, and
        // the variable list is sorted.
        let eq1 = unsafe { Modulo2Equation::from_parts(vec![0], 1_usize) };
        system.push(eq1);

        let solution = system.gaussian_elimination().unwrap();
        assert!(system.check(&solution));
    }

    #[test]
    fn test_gaussian_elimination_leading_unsolvable_row() {
        let mut system = Modulo2System::<usize>::new(1);
        // SAFETY: the empty variable list is sorted and contains no
        // out-of-bounds variable indices.
        let eq0 = unsafe { Modulo2Equation::from_parts(vec![], 1_usize) };
        system.push(eq0);
        // SAFETY: variable 0 is in bounds for this one-variable system, and
        // the variable list is sorted.
        let eq1 = unsafe { Modulo2Equation::from_parts(vec![0], 1_usize) };
        system.push(eq1);

        assert!(system.gaussian_elimination().is_err());
    }

    #[test]
    fn test_redundant_system() -> anyhow::Result<()> {
        let mut system = Modulo2System::<usize>::new(1);
        let eq0 = unsafe { Modulo2Equation::from_parts(vec![0], 0_usize) };
        system.push(eq0);
        let eq1 = unsafe { Modulo2Equation::from_parts(vec![0], 0_usize) };
        system.push(eq1);
        let solution = system.lazy_gaussian_elimination()?;
        assert!(system.check(&solution));
        Ok(())
    }

    #[test]
    fn test_small_system() -> anyhow::Result<()> {
        let mut system = Modulo2System::<usize>::new(11);
        let mut eq = unsafe { Modulo2Equation::from_parts(vec![1, 4, 10], 0) };
        system.push(eq);
        eq = unsafe { Modulo2Equation::from_parts(vec![1, 4, 9], 2) };
        system.push(eq);
        eq = unsafe { Modulo2Equation::from_parts(vec![0, 6, 8], 0) };
        system.push(eq);
        eq = unsafe { Modulo2Equation::from_parts(vec![0, 6, 9], 1) };
        system.push(eq);
        eq = unsafe { Modulo2Equation::from_parts(vec![2, 4, 8], 2) };
        system.push(eq);
        eq = unsafe { Modulo2Equation::from_parts(vec![2, 6, 10], 0) };
        system.push(eq);

        let solution = system.lazy_gaussian_elimination()?;
        assert!(system.check(&solution));
        Ok(())
    }

    #[test]
    fn test_var_to_vec_builder() {
        let system = unsafe {
            Modulo2System::from_parts(
                11,
                vec![
                    Modulo2Equation::from_parts(vec![1, 4, 10], 0_usize),
                    Modulo2Equation::from_parts(vec![1, 4, 9], 1),
                    Modulo2Equation::from_parts(vec![0, 6, 8], 2),
                    Modulo2Equation::from_parts(vec![0, 6, 9], 3),
                    Modulo2Equation::from_parts(vec![2, 4, 8], 4),
                    Modulo2Equation::from_parts(vec![2, 6, 10], 5),
                ],
            )
        };
        let mut backing: Vec<usize> = vec![];
        let (var_to_eqs, _weight, _priority) = system.setup(&mut backing);
        let expected_res = vec![
            vec![2, 3],
            vec![0, 1],
            vec![4, 5],
            vec![],
            vec![0, 1, 4],
            vec![],
            vec![2, 3, 5],
            vec![],
            vec![2, 4],
            vec![1, 3],
            vec![0, 5],
        ];
        var_to_eqs
            .iter()
            .zip(expected_res.iter())
            .for_each(|(v, e)| v.iter().zip(e.iter()).for_each(|(x, y)| assert_eq!(x, y)));
    }
}

1	/*
2	*
3	* SPDX-FileCopyrightText: 2025 Dario Moschetti
4	* SPDX-FileCopyrightText: 2025 Sebastiano Vigna
5	*
6	* SPDX-License-Identifier: Apache-2.0 OR LGPL-2.1-or-later
7	*/
8
9	//! Fast solution of random systems of linear equations with coefficients in
10	//! F₂.
11
12	#![allow(unexpected_cfgs)]
13	#![allow(clippy::comparison_chain)]
14	use std::ptr;
15
16	use crate::{
17	bit_vec,
18	traits::Word,
19	traits::{BitVecOps, BitVecOpsMut},
20	};
21	use anyhow::{Result, bail, ensure};
22	use arbitrary_chunks::ArbitraryChunks;
23
24	/// An equation on `W::BITS`-dimensional vectors with coefficients in F₂.
25	#[derive(Clone, Debug)]
26	pub struct Modulo2Equation<W: Word = usize> {
27	/// The variables in increasing order.
28	vars: Vec<u32>,
29	/// The constant term
30	c: W,
31	}
32
33	/// A system of [equations].
34	///
35	/// [equations]: struct.Modulo2Equation.html
36	#[derive(Clone, Debug)]
37	pub struct Modulo2System<W: Word = usize> {
38	/// The number of variables.
39	num_vars: usize,
40	/// The equations in the system.
41	equations: Vec<Modulo2Equation<W>>,
42	}
43
44	impl<W: Word> Modulo2Equation<W> {
45	/// Creates a new equation with given variables and constant term.
46	///
47	/// # Safety
48	///
49	/// `vars` must be strictly increasing. In particular, duplicate variables
50	/// are invalid: repeated coefficients cancel in F2 and must be removed by
51	/// the caller before construction.
52	pub unsafe fn from_parts(vars: Vec<u32>, c: W) -> Self {	352,391✔
53	debug_assert!(vars.windows(2).all(\|pair\| pair[0] < pair[1]));	2,466,689✔
54	Modulo2Equation { vars, c }
55	}
56
57	#[inline(always)]
58	unsafe fn add_ptr(	718,803✔
59	mut left: *const u32,
60	left_end: *const u32,
61	mut right: *const u32,
62	right_end: *const u32,
63	mut dst: *mut u32,
64	) -> *mut u32 {
65	unsafe {
66	while left != left_end && right != right_end {	437,822,446✔
67	let less = left <= right;	436,970,970✔
68	let more = left >= right;	436,970,970✔
69
70	let src = if less { left } else { right };	655,456,455✔
71	dst = src;	218,485,485✔
72
73	left = left.add(less as usize);	436,970,970✔
74	right = right.add(more as usize);	436,970,970✔
75	dst = dst.add((less ^ more) as usize);	655,456,455✔
76	}
77
78	let rem_left = left_end.offset_from(left) as usize;	2,156,409✔
79	ptr::copy_nonoverlapping(left, dst, rem_left);	2,875,212✔
80	dst = dst.add(rem_left);	1,437,606✔
81	let rem_right = right_end.offset_from(right) as usize;	2,156,409✔
82	ptr::copy_nonoverlapping(right, dst, rem_right);	2,875,212✔
83	dst = dst.add(rem_right);	1,437,606✔
84	dst	718,803✔
85	}
86	}
87
88	/// Adds `other` into `self`, using `scratch` as a working buffer to avoid
89	/// per-call allocation.
90	///
91	/// After the call, `scratch` holds what used to be `self.vars` (with its
92	/// original capacity intact), ready to be reused on the next call. Over
93	/// many calls, `scratch`'s capacity grows to the peak merge size and
94	/// further calls allocate nothing.
95	pub(crate) fn add_into(&mut self, other: &Modulo2Equation<W>, scratch: &mut Vec<u32>) {	718,802✔
96	scratch.clear();	1,437,604✔
97	scratch.reserve(self.vars.len() + other.vars.len());	3,594,010✔
98	let dst = scratch.as_mut_ptr();	2,156,406✔
99
100	let left_range = self.vars.as_ptr_range();	1,437,604✔
101	let right_range = other.vars.as_ptr_range();	1,437,604✔
102
103	unsafe {
104	let end = Self::add_ptr(
105	left_range.start,	1,437,604✔
106	left_range.end,	1,437,604✔
107	right_range.start,	1,437,604✔
108	right_range.end,	1,437,604✔
109	dst,	1,437,604✔
110	);
111	scratch.set_len(end.offset_from(dst) as usize);	2,875,208✔
112	}
113
114	self.c ^= other.c;	718,802✔
115	std::mem::swap(&mut self.vars, scratch);	2,156,406✔
116	}
117
118	pub fn add(&mut self, other: &Modulo2Equation<W>) {	1✔
119	let mut scratch = Vec::new();	2✔
120	self.add_into(other, &mut scratch);	4✔
121	}
122
123	/// Checks whether the equation is unsolvable.
124	fn is_unsolvable(&self) -> bool {	123,738✔
125	self.vars.is_empty() && self.c != W::ZERO	247,631✔
126	}
127
128	/// Checks whether the equation is an identity.
129	fn is_identity(&self) -> bool {	123,587✔
130	self.vars.is_empty() && self.c == W::ZERO	247,178✔
131	}
132
133	/// Evaluates the XOR of the values associated to the given variables.
134	fn eval_vars(vars: impl AsRef<[u32]>, values: impl AsRef<[W]>) -> W {	348,900✔
135	let mut sum = W::ZERO;	697,800✔
136	for &var in vars.as_ref() {	263,854,600✔
137	sum ^= values.as_ref()[var as usize];	394,735,200✔
138	}
139	sum	348,900✔
140	}
141	}
142
143	impl<W: Word> Modulo2System<W> {
144	pub fn new(num_vars: usize) -> Self {	8✔
145	Modulo2System {
146	num_vars,
147	equations: vec![],	8✔
148	}
149	}
150
151	pub const fn num_vars(&self) -> usize {	×
152	self.num_vars	×
153	}
154
155	pub fn num_equations(&self) -> usize {	42✔
156	self.equations.len()	84✔
157	}
158
159	/// Creates a new `Modulo2System` from existing equations.
160	///
161	/// # Safety
162	///
163	/// Every variable index in every equation must be strictly less than
164	/// `num_vars`. Equation variable lists must also satisfy
165	/// [`Modulo2Equation::from_parts`]'s strictly-increasing invariant.
166	pub unsafe fn from_parts(num_vars: usize, equations: Vec<Modulo2Equation<W>>) -> Self {	236✔
167	Modulo2System {
168	num_vars,
169	equations,
170	}
171	}
172
173	/// Adds an equation to the system.
174	///
175	/// # Panics
176	///
177	/// Panics if the equation contains a variable outside this system.
178	pub fn push(&mut self, equation: Modulo2Equation<W>) {	18✔
179	if let Some(&last_var) = equation.vars.last() {	33✔
180	let last_var = usize::try_from(last_var).expect("variable index does not fit in usize");	75✔
181	assert!(	15✔
182	last_var < self.num_vars,	15✔
NEW 183	"variable {last_var} is out of bounds for a system with {} variables",	×
NEW 184	self.num_vars	×
185	);
186	}
187	self.equations.push(equation);	51✔
188	}
189
190	/// Checks if a given solution satisfies the system of equations.
191	pub fn check(&self, solution: &[W]) -> bool {	4✔
192	assert_eq!(	4✔
193	solution.len(),	8✔
194	self.num_vars,	×
195	"The number of variables in the solution ({}) does not match the number of variables in the system ({})",	×
196	solution.len(),	×
197	self.num_vars	×
198	);
199	self.equations	4✔
200	.iter()
201	.all(\|eq\| eq.c == Modulo2Equation::<W>::eval_vars(&eq.vars, solution))	37✔
202	}
203
204	/// Puts the system into echelon form.
205	fn echelon_form(&mut self) -> Result<()> {	49✔
206	let equations = &mut self.equations;	98✔
207	if equations.is_empty() {	98✔
208	return Ok(());	2✔
209	}
210	let mut scratch: Vec<u32> = Vec::new();	141✔
211	'main: for i in 0..equations.len() - 1 {	8,092✔
212	// A pivot row reduced to [] = c: unsolvable if c != 0, otherwise
213	// a redundant identity that takes no part in elimination.
214	if equations[i].vars.is_empty() {	16,090✔
215	ensure!(!equations[i].is_unsolvable(), "System is unsolvable");	4✔
216	continue;	1✔
217	}
218	for j in i + 1..equations.len() {	2,793,421✔
219	let (left, right) = equations.split_at_mut(j);	11,109,340✔
220	let eq_i = &mut left[i];	5,554,670✔
221	let eq_j = &right[0];	5,554,670✔
222
223	// A later row reduced to [] = c: unsolvable if c != 0, else a
224	// redundant identity to skip. Without this, eq_j.vars[0] below
225	// would panic on the empty variable list.
226	if eq_j.vars.is_empty() {	5,554,670✔
227	ensure!(!eq_j.is_unsolvable(), "System is unsolvable");	2✔
228	continue;	×
229	}
230
231	let first_var_j = eq_j.vars[0];	5,554,668✔
232
233	if eq_i.vars[0] == first_var_j {	2,777,334✔
234	eq_i.add_into(eq_j, &mut scratch);	429,196✔
235	if eq_i.is_unsolvable() {	214,598✔
236	bail!("System is unsolvable");	1✔
237	}
238	if eq_i.is_identity() {	214,596✔
239	continue 'main;	×
240	}
241	}
242
243	if eq_i.vars[0] > first_var_j {	2,777,333✔
244	equations.swap(i, j)	7,221,213✔
245	}
246	}
247	}
248	Ok(())	44✔
249	}
250
251	/// Solves the system using Gaussian elimination.
252	pub fn gaussian_elimination(&mut self) -> Result<Vec<W>> {	49✔
253	self.echelon_form()?;	101✔
254	// An equation reduced to [] = c with c != 0 has no solution; the
255	// back-substitution below would index eq.vars[0] on an empty variable
256	// list. echelon_form only checks emptiness for rows before the last.
257	for eq in &self.equations {	8,130✔
258	ensure!(!eq.is_unsolvable(), "System is unsolvable");	16,168✔
259	}
260	let mut solution = vec![W::ZERO; self.num_vars];	138✔
261	self.equations	46✔
262	.iter()
263	.rev()
264	.filter(\|eq\| !eq.is_identity())	16,214✔
265	.for_each(\|eq\| {	8,129✔
266	solution[eq.vars[0] as usize] =	24,249✔
267	eq.c ^ Modulo2Equation::<W>::eval_vars(&eq.vars, &solution);	24,249✔
268	});
269	Ok(solution)	46✔
270	}
271
272	/// Builds the data structures needed for [lazy Gaussian elimination].
273	///
274	/// This method returns the variable-to-equation mapping, the weight of each
275	/// variable (the number of equations in which it appears), and the priority
276	/// of each equation (the number of variables in it). The
277	/// variable-to-equation mapping is materialized as a vector of mutable
278	/// slices, each of which points inside a provided backing vector. This
279	/// approach greatly reduces the number of allocations.
280	///
281	/// [lazy Gaussian elimination]: https://doi.org/10.1016/j.ic.2020.104517
282	fn setup<'a>(	194✔
283	&self,
284	backing: &'a mut Vec<usize>,
285	) -> (Vec<&'a mut [usize]>, Vec<usize>, Vec<usize>) {
286	let mut weight = vec![0; self.num_vars];	582✔
287	let mut priority = vec![0; self.equations.len()];	776✔
288
289	let mut total_vars = 0;	388✔
290	for (i, equation) in self.equations.iter().enumerate() {	705,150✔
291	priority[i] = equation.vars.len();	704,762✔
292	total_vars += equation.vars.len();	352,381✔
293	for &var in &equation.vars {	2,466,647✔
294	weight[var as usize] += 1;	1,057,133✔
295	}
296	}
297
298	backing.resize(total_vars, 0);	582✔
299	let mut var_to_eq: Vec<&mut [usize]> = Vec::with_capacity(self.num_vars);	776✔
300
301	backing.arbitrary_chunks_mut(&weight).for_each(\|chunk\| {	695,824✔
302	var_to_eq.push(chunk);	2,085,726✔
303	});
304
305	let mut pos = vec![0usize; self.num_vars];	582✔
306	for (i, equation) in self.equations.iter().enumerate() {	705,150✔
307	for &var in &equation.vars {	2,466,647✔
308	let var = var as usize;	3,171,399✔
309	var_to_eq[var][pos[var]] = i;	3,171,399✔
310	pos[var] += 1;	1,057,133✔
311	}
312	}
313
314	(var_to_eq, weight, priority)	388✔
315	}
316
317	/// Solves the system using [lazy Gaussian elimination].
318	///
319	/// [lazy Gaussian elimination]: https://doi.org/10.1016/j.ic.2020.104517
320	pub fn lazy_gaussian_elimination(&mut self) -> Result<Vec<W>> {	193✔
321	let num_vars = self.num_vars;	386✔
322	let num_equations = self.equations.len();	579✔
323
324	if num_equations == 0 {	193✔
325	return Ok(vec![W::ZERO; num_vars]);	×
326	}
327
328	let mut backing = vec![];	386✔
329	let (var_to_eqs, mut weight, mut priority);	×
330	(var_to_eqs, weight, priority) = self.setup(&mut backing);	965✔
331
332	let mut variables = vec![0; num_vars];	579✔
333	{
334	let mut count = vec![0; num_equations + 1];	579✔
335
336	for x in 0..num_vars {	695,424✔
337	count[weight[x]] += 1	1,390,462✔
338	}
339	for i in 1..count.len() {	352,568✔
340	count[i] += count[i - 1]	704,750✔
341	}
342	for i in (0..num_vars).rev() {	1,390,848✔
343	count[weight[i]] -= 1;	2,085,693✔
344	variables[count[weight[i]]] = i;	2,085,693✔
345	}
346	}
347
348	let mut equation_list: Vec<usize> = (0..priority.len())	579✔
349	.rev()
350	.filter(\|&x\| priority[x] <= 1)	352,568✔
351	.collect();
352
353	let mut dense = vec![];	386✔
354	let mut solved = vec![];	386✔
355	let mut pivots = vec![];	386✔
356
357	let equations = &mut self.equations;	386✔
358	let mut idle = bit_vec![true; num_vars];	579✔
359	let mut scratch: Vec<u32> = Vec::new();	579✔
360
361	let mut remaining = equations.len();	579✔
362
363	while remaining != 0 {	369,611✔
364	if equation_list.is_empty() {	739,130✔
365	let mut var = variables.pop().unwrap();	71,056✔
366	while weight[var] == 0 {	336,038✔
367	var = variables.pop().unwrap()	636,548✔
368	}
369	idle.set(var, false);	53,292✔
370	var_to_eqs[var].as_ref().iter().for_each(\|&eq\| {	154,225✔
371	priority[eq] -= 1;	100,933✔
372	if priority[eq] == 1 {	100,933✔
373	equation_list.push(eq)	53,265✔
374	}
375	});
376	} else {
377	remaining -= 1;	351,801✔
378	let first = equation_list.pop().unwrap();	1,407,204✔
379
380	if priority[first] == 0 {	351,801✔
381	let equation = &mut equations[first];	16,704✔
382	if equation.is_unsolvable() {	16,704✔
383	bail!("System is unsolvable")	147✔
384	}
385	if equation.is_identity() {	16,410✔
386	continue;	3✔
387	}
388	dense.push(equation.to_owned());	32,808✔
389	} else if priority[first] == 1 {	343,449✔
390	let pivot = equations[first]	686,898✔
391	.vars	343,449✔
392	.iter()
393	.copied()
394	.find(\|x\| idle.get(*x as usize))	23,532,690✔
395	.expect("Missing expected idle variable in equation");
396	pivots.push(pivot as usize);	1,030,347✔
397	solved.push(first);	1,030,347✔
398	weight[pivot as usize] = 0;	343,449✔
399	for &eq in var_to_eqs[pivot as usize]	954,951✔
400	.as_ref()	343,449✔
401	.iter()	343,449✔
402	.filter(\|&&eq_idx\| eq_idx != first)	2,253,351✔
403	{
404	let lo = eq.min(first);	2,446,008✔
405	let hi = eq.max(first);	2,446,008✔
406	let (left, right) = equations.split_at_mut(hi);	2,446,008✔
407	if eq < first {	1,010,251✔
408	left[lo].add_into(&right[0], &mut scratch);	1,196,247✔
409	} else {
410	right[0].add_into(&left[lo], &mut scratch);	638,259✔
411	}
412
413	priority[eq] -= 1;	611,502✔
414	if priority[eq] == 1 {	611,502✔
415	equation_list.push(eq)	1,002,312✔
416	}
417	}
418	}
419	}
420	}
421
422	// SAFETY: the usage of the method is safe, as the equations have the
423	// right number of variables
424	let mut dense_system = unsafe { Modulo2System::from_parts(num_vars, dense) };	184✔
425	let mut solution = dense_system.gaussian_elimination()?;	138✔
426
427	for i in 0..solved.len() {	340,851✔
428	let eq = &equations[solved[i]];	1,022,418✔
429	let pivot = pivots[i];	681,612✔
430	assert!(solution[pivot] == W::ZERO);	681,612✔
431	solution[pivot] = eq.c ^ Modulo2Equation::<W>::eval_vars(&eq.vars, &solution);	1,363,224✔
432	}
433
434	Ok(solution)	45✔
435	}
436	}
437
438	#[cfg(test)]
439	mod tests {
440	use super::*;
441
442	#[test]
443	fn test_add_ptr() {
444	let a = [0, 1, 3, 4];
445	let b = [0, 2, 4, 5];
446	let mut c = Vec::with_capacity(a.len() + b.len());
447
448	let ra = a.as_ptr_range();
449	let rb = b.as_ptr_range();
450	let mut dst = c.as_mut_ptr();
451	unsafe {
452	dst = Modulo2Equation::<u32>::add_ptr(ra.start, ra.end, rb.start, rb.end, dst);
453	assert_eq!(dst.offset_from(c.as_ptr()), 4);
454	c.set_len(4);
455	assert_eq!(c, vec![1, 2, 3, 5]);
456	}
457	}
458
459	#[test]
460	fn test_equation_builder() {
461	let eq = unsafe { Modulo2Equation::from_parts(vec![0, 1, 2], 3_usize) };
462	assert_eq!(eq.vars.len(), 3);
463	assert_eq!(eq.vars.to_vec(), vec![0, 1, 2]);
464	}
465
466	#[test]
467	fn test_equation_addition() {
468	let mut eq0 = unsafe { Modulo2Equation::from_parts(vec![1, 4, 9], 3_usize) };
469	let eq1 = unsafe { Modulo2Equation::from_parts(vec![1, 4, 10], 3_usize) };
470	eq0.add(&eq1);
471	assert_eq!(eq0.vars, vec![9, 10]);
472	assert_eq!(eq0.c, 0);
473	}
474
475	#[test]
476	fn test_system_one_equation() -> anyhow::Result<()> {
477	let mut system = Modulo2System::<usize>::new(2);
478	let eq = unsafe { Modulo2Equation::from_parts(vec![0], 3_usize) };
479	system.push(eq);
480	let solution = system.lazy_gaussian_elimination()?;
481	assert!(system.check(&solution));
482	Ok(())
483	}
484
485	#[test]
486	#[should_panic(expected = "out of bounds")]
487	fn test_push_rejects_out_of_range_variable() {
488	let mut system = Modulo2System::<usize>::new(1);
489	// SAFETY: a one-element variable list is strictly increasing.
490	let equation = unsafe { Modulo2Equation::from_parts(vec![1], 0_usize) };
491	system.push(equation);
492	}
493
494	#[test]
495	fn test_impossible_system() {
496	let mut system = Modulo2System::<usize>::new(1);
497	let eq0 = unsafe { Modulo2Equation::from_parts(vec![0], 0_usize) };
498	system.push(eq0);
499	let eq1 = unsafe { Modulo2Equation::from_parts(vec![0], 1_usize) };
500	system.push(eq1);
501	let solution = system.lazy_gaussian_elimination();
502	assert!(solution.is_err());
503	}
504
505	#[test]
506	fn test_gaussian_elimination_last_unsolvable_equation() {
507	let mut system = Modulo2System::<usize>::new(1);
508	// SAFETY: variable 0 is in bounds for this one-variable system, and
509	// the variable list is sorted.
510	let eq0 = unsafe { Modulo2Equation::from_parts(vec![0], 0_usize) };
511	system.push(eq0);
512	// SAFETY: the empty variable list is sorted and contains no
513	// out-of-bounds variable indices.
514	let eq1 = unsafe { Modulo2Equation::from_parts(vec![], 1_usize) };
515	system.push(eq1);
516
517	let solution = system.gaussian_elimination();
518	assert!(solution.is_err());
519	}
520
521	#[test]
522	fn test_gaussian_elimination_leading_identity_row() {
523	let mut system = Modulo2System::<usize>::new(1);
524	// SAFETY: the empty variable list is sorted and contains no
525	// out-of-bounds variable indices.
526	let eq0 = unsafe { Modulo2Equation::from_parts(vec![], 0_usize) };
527	system.push(eq0);
528	// SAFETY: variable 0 is in bounds for this one-variable system, and
529	// the variable list is sorted.
530	let eq1 = unsafe { Modulo2Equation::from_parts(vec![0], 1_usize) };
531	system.push(eq1);
532
533	let solution = system.gaussian_elimination().unwrap();
534	assert!(system.check(&solution));
535	}
536
537	#[test]
538	fn test_gaussian_elimination_leading_unsolvable_row() {
539	let mut system = Modulo2System::<usize>::new(1);
540	// SAFETY: the empty variable list is sorted and contains no
541	// out-of-bounds variable indices.
542	let eq0 = unsafe { Modulo2Equation::from_parts(vec![], 1_usize) };
543	system.push(eq0);
544	// SAFETY: variable 0 is in bounds for this one-variable system, and
545	// the variable list is sorted.
546	let eq1 = unsafe { Modulo2Equation::from_parts(vec![0], 1_usize) };
547	system.push(eq1);
548
549	assert!(system.gaussian_elimination().is_err());
550	}
551
552	#[test]
553	fn test_redundant_system() -> anyhow::Result<()> {
554	let mut system = Modulo2System::<usize>::new(1);
555	let eq0 = unsafe { Modulo2Equation::from_parts(vec![0], 0_usize) };
556	system.push(eq0);
557	let eq1 = unsafe { Modulo2Equation::from_parts(vec![0], 0_usize) };
558	system.push(eq1);
559	let solution = system.lazy_gaussian_elimination()?;
560	assert!(system.check(&solution));
561	Ok(())
562	}
563
564	#[test]
565	fn test_small_system() -> anyhow::Result<()> {
566	let mut system = Modulo2System::<usize>::new(11);
567	let mut eq = unsafe { Modulo2Equation::from_parts(vec![1, 4, 10], 0) };
568	system.push(eq);
569	eq = unsafe { Modulo2Equation::from_parts(vec![1, 4, 9], 2) };
570	system.push(eq);
571	eq = unsafe { Modulo2Equation::from_parts(vec![0, 6, 8], 0) };
572	system.push(eq);
573	eq = unsafe { Modulo2Equation::from_parts(vec![0, 6, 9], 1) };
574	system.push(eq);
575	eq = unsafe { Modulo2Equation::from_parts(vec![2, 4, 8], 2) };
576	system.push(eq);
577	eq = unsafe { Modulo2Equation::from_parts(vec![2, 6, 10], 0) };
578	system.push(eq);
579
580	let solution = system.lazy_gaussian_elimination()?;
581	assert!(system.check(&solution));
582	Ok(())
583	}
584
585	#[test]
586	fn test_var_to_vec_builder() {
587	let system = unsafe {
588	Modulo2System::from_parts(
589	11,
590	vec![
591	Modulo2Equation::from_parts(vec![1, 4, 10], 0_usize),
592	Modulo2Equation::from_parts(vec![1, 4, 9], 1),
593	Modulo2Equation::from_parts(vec![0, 6, 8], 2),
594	Modulo2Equation::from_parts(vec![0, 6, 9], 3),
595	Modulo2Equation::from_parts(vec![2, 4, 8], 4),
596	Modulo2Equation::from_parts(vec![2, 6, 10], 5),
597	],
598	)
599	};
600	let mut backing: Vec<usize> = vec![];
601	let (var_to_eqs, _weight, _priority) = system.setup(&mut backing);
602	let expected_res = vec![
603	vec![2, 3],
604	vec![0, 1],
605	vec![4, 5],
606	vec![],
607	vec![0, 1, 4],
608	vec![],
609	vec![2, 3, 5],
610	vec![],
611	vec![2, 4],
612	vec![1, 3],
613	vec![0, 5],
614	];
615	var_to_eqs
616	.iter()
617	.zip(expected_res.iter())
618	.for_each(\|(v, e)\| v.iter().zip(e.iter()).for_each(\|(x, y)\| assert_eq!(x, y)));
619	}
620	}

vigna / sux-rs / 29426297655

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous