21985226363

Committed 13 Feb 2026 11:29AM UTC coverage: 78.4% (-0.1%) from 78.499%

Build # 21985226363

Build Type

push

github

Committed by

vigna

Commit Message

fmt

Coverage Stats

3 of 3 new or added lines in 1 file covered. (100.0%)

94 existing lines in 5 files now uncovered.

4773 of 6088 relevant lines covered (78.4%)

120432792.49 hits per line

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

95.0

/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](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
    ///
    /// The caller must ensure that the variables are sorted.
    pub unsafe fn from_parts(vars: Vec<u32>, c: W) -> Self {
        debug_assert!(vars.iter().is_sorted());
        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
        }
    }

    pub fn add(&mut self, other: &Modulo2Equation<W>) {
        let left_range = self.vars.as_ptr_range();
        let left = left_range.start;
        let left_end = left_range.end;
        let right_range = other.vars.as_ptr_range();
        let right = right_range.start;
        let right_end = right_range.end;
        let mut vars = Vec::with_capacity(self.vars.len() + other.vars.len());
        let dst = vars.as_mut_ptr();

        unsafe {
            let copied =
                Self::add_ptr(left, left_end, right, right_end, dst).offset_from(dst) as usize;
            vars.set_len(copied);
        }

        self.vars = vars;
        self.c ^= other.c;
    }

    /// 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 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
    ///
    /// The caller must ensure that variables in each equation match the number
    /// of variables in the system.
    pub unsafe fn from_parts(num_vars: usize, equations: Vec<Modulo2Equation<W>>) -> Self {
        Modulo2System {
            num_vars,
            equations,
        }
    }

    /// Adds an equation to the system.
    pub fn push(&mut self, equation: Modulo2Equation<W>) {
        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(());
        }
        'main: for i in 0..equations.len() - 1 {
            ensure!(!equations[i].vars.is_empty());
            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];

                let first_var_j = eq_j.vars[0];

                if eq_i.vars[0] == first_var_j {
                    eq_i.add(eq_j);
                    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()?;
        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](https://doi.org/10.1016/j.ic.2020.104517).
    ///
    /// 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.
    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](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 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(&right[0]);
                        } else {
                            right[0].add(&left[lo]);
                        }

                        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]
    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_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](struct.Modulo2Equation.html).
34	#[derive(Clone, Debug)]
35	pub struct Modulo2System<W: Word = usize> {
36	/// The number of variables.
37	num_vars: usize,
38	/// The equations in the system.
39	equations: Vec<Modulo2Equation<W>>,
40	}
41
42	impl<W: Word> Modulo2Equation<W> {
43	/// Creates a new equation with given variables and constant term.
44	///
45	/// # Safety
46	///
47	/// The caller must ensure that the variables are sorted.
48	pub unsafe fn from_parts(vars: Vec<u32>, c: W) -> Self {	646,958✔
49	debug_assert!(vars.iter().is_sorted());	1,940,874✔
50	Modulo2Equation { vars, c }
51	}
52
53	#[inline(always)]
54	unsafe fn add_ptr(	1,727,775✔
55	mut left: *const u32,
56	left_end: *const u32,
57	mut right: *const u32,
58	right_end: *const u32,
59	mut dst: *mut u32,
60	) -> *mut u32 {
61	unsafe {
62	while left != left_end && right != right_end {	2,109,473,469✔
63	let less = left <= right;	2,107,357,158✔
64	let more = left >= right;	2,107,357,158✔
65
66	let src = if less { left } else { right };	2,147,483,647✔
67	dst = src;	1,053,678,579✔
68
69	left = left.add(less as usize);	2,107,357,158✔
70	right = right.add(more as usize);	2,107,357,158✔
71	dst = dst.add((less ^ more) as usize);	2,147,483,647✔
72	}
73
74	let rem_left = left_end.offset_from(left) as usize;	5,183,325✔
75	ptr::copy_nonoverlapping(left, dst, rem_left);	6,911,100✔
76	dst = dst.add(rem_left);	3,455,550✔
77	let rem_right = right_end.offset_from(right) as usize;	5,183,325✔
78	ptr::copy_nonoverlapping(right, dst, rem_right);	6,911,100✔
79	dst = dst.add(rem_right);	3,455,550✔
80	dst	1,727,775✔
81	}
82	}
83
84	pub fn add(&mut self, other: &Modulo2Equation<W>) {	1,727,774✔
85	let left_range = self.vars.as_ptr_range();	3,455,548✔
86	let left = left_range.start;	3,455,548✔
87	let left_end = left_range.end;	3,455,548✔
88	let right_range = other.vars.as_ptr_range();	3,455,548✔
89	let right = right_range.start;	3,455,548✔
90	let right_end = right_range.end;	3,455,548✔
91	let mut vars = Vec::with_capacity(self.vars.len() + other.vars.len());	8,638,870✔
92	let dst = vars.as_mut_ptr();	5,183,322✔
93
94	unsafe {
95	let copied =	3,455,548✔
96	Self::add_ptr(left, left_end, right, right_end, dst).offset_from(dst) as usize;	13,822,192✔
97	vars.set_len(copied);	3,455,548✔
98	}
99
100	self.vars = vars;	3,455,548✔
101	self.c ^= other.c;	1,727,774✔
102	}
103
104	/// Checks whether the equation is unsolvable.
105	fn is_unsolvable(&self) -> bool {	675,502✔
106	self.vars.is_empty() && self.c != W::ZERO	1,351,148✔
107	}
108
109	/// Checks whether the equation is an identity.
110	fn is_identity(&self) -> bool {	693,796✔
111	self.vars.is_empty() && self.c == W::ZERO	1,387,593✔
112	}
113
114	/// Evaluates the XOR of the values associated to the given variables.
115	fn eval_vars(vars: impl AsRef<[u32]>, values: impl AsRef<[W]>) -> W {	599,477✔
116	let mut sum = W::ZERO;	1,198,954✔
117	for &var in vars.as_ref() {	609,298,840✔
118	sum ^= values.as_ref()[var as usize];	912,149,829✔
119	}
120	sum	599,477✔
121	}
122	}
123
124	impl<W: Word> Modulo2System<W> {
125	pub fn new(num_vars: usize) -> Self {	4✔
126	Modulo2System {
127	num_vars,
128	equations: vec![],	4✔
129	}
130	}
131
132	pub fn num_vars(&self) -> usize {	×
133	self.num_vars	×
134	}
135
136	pub fn num_equations(&self) -> usize {	27✔
137	self.equations.len()	54✔
138	}
139
140	/// Creates a new `Modulo2System` from existing equations.
141	///
142	/// # Safety
143	///
144	/// The caller must ensure that variables in each equation match the number
145	/// of variables in the system.
146	pub unsafe fn from_parts(num_vars: usize, equations: Vec<Modulo2Equation<W>>) -> Self {	237✔
147	Modulo2System {
148	num_vars,
149	equations,
150	}
151	}
152
153	/// Adds an equation to the system.
154	pub fn push(&mut self, equation: Modulo2Equation<W>) {	11✔
155	self.equations.push(equation);	33✔
156	}
157
158	/// Checks if a given solution satisfies the system of equations.
159	pub fn check(&self, solution: &[W]) -> bool {	3✔
160	assert_eq!(	3✔
161	solution.len(),	6✔
162	self.num_vars,	×
163	"The number of variables in the solution ({}) does not match the number of variables in the system ({})",	×
164	solution.len(),	×
165	self.num_vars	×
166	);
167	self.equations	3✔
168	.iter()
169	.all(\|eq\| eq.c == Modulo2Equation::<W>::eval_vars(&eq.vars, solution))	30✔
170	}
171
172	/// Puts the system into echelon form.
173	fn echelon_form(&mut self) -> Result<()> {	67✔
174	let equations = &mut self.equations;	134✔
175	if equations.is_empty() {	134✔
176	return Ok(());	2✔
177	}
178	'main: for i in 0..equations.len() - 1 {	18,925✔
179	ensure!(!equations[i].vars.is_empty());	37,720✔
180	for j in i + 1..equations.len() {	14,553,983✔
181	let (left, right) = equations.split_at_mut(j);	58,065,052✔
182	let eq_i = &mut left[i];	29,032,526✔
183	let eq_j = &right[0];	29,032,526✔
184
185	let first_var_j = eq_j.vars[0];	29,032,526✔
186
187	if eq_i.vars[0] == first_var_j {	14,516,263✔
188	eq_i.add(eq_j);	1,963,554✔
189	if eq_i.is_unsolvable() {	1,309,036✔
190	bail!("System is unsolvable");	37✔
191	}
192	if eq_i.is_identity() {	1,308,962✔
UNCOV 193	continue 'main;	×
194	}
195	}
196
197	if eq_i.vars[0] > first_var_j {	14,516,226✔
198	equations.swap(i, j)	29,196,972✔
199	}
200	}
201	}
202	Ok(())	28✔
203	}
204
205	/// Solves the system using Gaussian elimination.
206	pub fn gaussian_elimination(&mut self) -> Result<Vec<W>> {	67✔
207	self.echelon_form()?;	171✔
208	let mut solution = vec![W::ZERO; self.num_vars];	90✔
209	self.equations	30✔
210	.iter()
211	.rev()
212	.filter(\|eq\| !eq.is_identity())	36,904✔
213	.for_each(\|eq\| {	18,467✔
214	solution[eq.vars[0] as usize] =	55,311✔
215	eq.c ^ Modulo2Equation::<W>::eval_vars(&eq.vars, &solution);	55,311✔
216	});
217	Ok(solution)	30✔
218	}
219
220	/// Builds the data structures needed for [lazy Gaussian
221	/// elimination](https://doi.org/10.1016/j.ic.2020.104517).
222	///
223	/// This method returns the variable-to-equation mapping, the weight of each
224	/// variable (the number of equations in which it appears), and the priority
225	/// of each equation (the number of variables in it). The
226	/// variable-to-equation mapping is materialized as a vector of mutable
227	/// slices, each of which points inside a provided backing vector. This
228	/// approach greatly reduces the number of allocations.
229	fn setup<'a>(	174✔
230	&self,
231	backing: &'a mut Vec<usize>,
232	) -> (Vec<&'a mut [usize]>, Vec<usize>, Vec<usize>) {
233	let mut weight = vec![0; self.num_vars];	522✔
234	let mut priority = vec![0; self.equations.len()];	696✔
235
236	let mut total_vars = 0;	348✔
237	for (i, equation) in self.equations.iter().enumerate() {	1,294,258✔
238	priority[i] = equation.vars.len();	1,293,910✔
239	total_vars += equation.vars.len();	646,955✔
240	for &var in &equation.vars {	4,528,665✔
241	weight[var as usize] += 1;	1,940,855✔
242	}
243	}
244
245	backing.resize(total_vars, 0);	522✔
246	let mut var_to_eq: Vec<&mut [usize]> = Vec::with_capacity(self.num_vars);	696✔
247
248	backing.arbitrary_chunks_mut(&weight).for_each(\|chunk\| {	1,238,636✔
249	var_to_eq.push(chunk);	3,714,342✔
250	});
251
252	let mut pos = vec![0usize; self.num_vars];	522✔
253	for (i, equation) in self.equations.iter().enumerate() {	1,294,258✔
254	for &var in &equation.vars {	4,528,665✔
255	let var = var as usize;	5,822,565✔
256	var_to_eq[var][pos[var]] = i;	5,822,565✔
257	pos[var] += 1;	1,940,855✔
258	}
259	}
260
261	(var_to_eq, weight, priority)	348✔
262	}
263
264	/// Solves the system using [lazy Gaussian
265	/// elimination](https://doi.org/10.1016/j.ic.2020.104517).
266	pub fn lazy_gaussian_elimination(&mut self) -> Result<Vec<W>> {	173✔
267	let num_vars = self.num_vars;	346✔
268	let num_equations = self.equations.len();	519✔
269
270	if num_equations == 0 {	173✔
UNCOV 271	return Ok(vec![W::ZERO; num_vars]);	×
272	}
273
274	let mut backing = vec![];	346✔
UNCOV 275	let (var_to_eqs, mut weight, mut priority);	×
276	(var_to_eqs, weight, priority) = self.setup(&mut backing);	865✔
277
278	let mut variables = vec![0; num_vars];	519✔
279	{
280	let mut count = vec![0; num_equations + 1];	519✔
281
282	for x in 0..num_vars {	1,238,276✔
283	count[weight[x]] += 1	2,476,206✔
284	}
285	for i in 1..count.len() {	647,122✔
286	count[i] += count[i - 1]	1,293,898✔
287	}
288	for i in (0..num_vars).rev() {	2,476,552✔
289	count[weight[i]] -= 1;	3,714,309✔
290	variables[count[weight[i]]] = i;	3,714,309✔
291	}
292	}
293
294	let mut equation_list: Vec<usize> = (0..priority.len())	519✔
295	.rev()
296	.filter(\|&x\| priority[x] <= 1)	647,122✔
297	.collect();
298
299	let mut dense = vec![];	346✔
300	let mut solved = vec![];	346✔
301	let mut pivots = vec![];	346✔
302
303	let equations = &mut self.equations;	346✔
304	let mut idle = bit_vec![true; num_vars];	519✔
305
306	let mut remaining = equations.len();	519✔
307
308	while remaining != 0 {	677,553✔
309	if equation_list.is_empty() {	1,354,972✔
310	let mut var = variables.pop().unwrap();	162,148✔
311	while weight[var] == 0 {	612,392✔
312	var = variables.pop().unwrap()	1,143,710✔
313	}
314	idle.set(var, false);	121,611✔
315	var_to_eqs[var].as_ref().iter().for_each(\|&eq\| {	347,751✔
316	priority[eq] -= 1;	226,140✔
317	if priority[eq] == 1 {	226,140✔
318	equation_list.push(eq)	126,240✔
319	}
320	});
321	} else {
322	remaining -= 1;	636,949✔
323	let first = equation_list.pop().unwrap();	2,547,796✔
324
325	if priority[first] == 0 {	636,949✔
326	let equation = &mut equations[first];	41,968✔
327	if equation.is_unsolvable() {	41,968✔
328	bail!("System is unsolvable")	106✔
329	}
330	if equation.is_identity() {	41,756✔
331	continue;	1✔
332	}
333	dense.push(equation.to_owned());	83,508✔
334	} else if priority[first] == 1 {	615,965✔
335	let pivot = equations[first]	1,231,930✔
336	.vars	615,965✔
337	.iter()
338	.copied()
339	.find(\|x\| idle.get(*x as usize))	440,004,560✔
340	.expect("Missing expected idle variable in equation");
341	pivots.push(pivot as usize);	1,847,895✔
342	solved.push(first);	1,847,895✔
343	weight[pivot as usize] = 0;	615,965✔
344	for &eq in var_to_eqs[pivot as usize]	1,689,220✔
345	.as_ref()	615,965✔
346	.iter()	615,965✔
347	.filter(\|&&eq_idx\| eq_idx != first)	3,994,405✔
348	{
349	let lo = eq.min(first);	4,293,020✔
350	let hi = eq.max(first);	4,293,020✔
351	let (left, right) = equations.split_at_mut(hi);	4,293,020✔
352	if eq < first {	1,702,202✔
353	left[lo].add(&right[0]);	1,257,894✔
354	} else {
355	right[0].add(&left[lo]);	888,616✔
356	}
357
358	priority[eq] -= 1;	1,073,255✔
359	if priority[eq] == 1 {	1,073,255✔
360	equation_list.push(eq)	1,785,948✔
361	}
362	}
363	}
364	}
365	}
366
367	// SAFETY: the usage of the method is safe, as the equations have the
368	// right number of variables
369	let mut dense_system = unsafe { Modulo2System::from_parts(num_vars, dense) };	268✔
370	let mut solution = dense_system.gaussian_elimination()?;	201✔
371
372	for i in 0..solved.len() {	581,061✔
373	let eq = &equations[solved[i]];	1,743,093✔
374	let pivot = pivots[i];	1,162,062✔
375	assert!(solution[pivot] == W::ZERO);	1,162,062✔
376	solution[pivot] = eq.c ^ Modulo2Equation::<W>::eval_vars(&eq.vars, &solution);	2,324,124✔
377	}
378
379	Ok(solution)	30✔
380	}
381	}
382
383	#[cfg(test)]
384	mod tests {
385	use super::*;
386
387	#[test]
388	fn test_add_ptr() {
389	let a = [0, 1, 3, 4];
390	let b = [0, 2, 4, 5];
391	let mut c = Vec::with_capacity(a.len() + b.len());
392
393	let ra = a.as_ptr_range();
394	let rb = b.as_ptr_range();
395	let mut dst = c.as_mut_ptr();
396	unsafe {
397	dst = Modulo2Equation::<u32>::add_ptr(ra.start, ra.end, rb.start, rb.end, dst);
398	assert_eq!(dst.offset_from(c.as_ptr()), 4);
399	c.set_len(4);
400	assert_eq!(c, vec![1, 2, 3, 5]);
401	}
402	}
403
404	#[test]
405	fn test_equation_builder() {
406	let eq = unsafe { Modulo2Equation::from_parts(vec![0, 1, 2], 3_usize) };
407	assert_eq!(eq.vars.len(), 3);
408	assert_eq!(eq.vars.to_vec(), vec![0, 1, 2]);
409	}
410
411	#[test]
412	fn test_equation_addition() {
413	let mut eq0 = unsafe { Modulo2Equation::from_parts(vec![1, 4, 9], 3_usize) };
414	let eq1 = unsafe { Modulo2Equation::from_parts(vec![1, 4, 10], 3_usize) };
415	eq0.add(&eq1);
416	assert_eq!(eq0.vars, vec![9, 10]);
417	assert_eq!(eq0.c, 0);
418	}
419
420	#[test]
421	fn test_system_one_equation() -> anyhow::Result<()> {
422	let mut system = Modulo2System::<usize>::new(2);
423	let eq = unsafe { Modulo2Equation::from_parts(vec![0], 3_usize) };
424	system.push(eq);
425	let solution = system.lazy_gaussian_elimination()?;
426	assert!(system.check(&solution));
427	Ok(())
428	}
429
430	#[test]
431	fn test_impossible_system() {
432	let mut system = Modulo2System::<usize>::new(1);
433	let eq0 = unsafe { Modulo2Equation::from_parts(vec![0], 0_usize) };
434	system.push(eq0);
435	let eq1 = unsafe { Modulo2Equation::from_parts(vec![0], 1_usize) };
436	system.push(eq1);
437	let solution = system.lazy_gaussian_elimination();
438	assert!(solution.is_err());
439	}
440
441	#[test]
442	fn test_redundant_system() -> anyhow::Result<()> {
443	let mut system = Modulo2System::<usize>::new(1);
444	let eq0 = unsafe { Modulo2Equation::from_parts(vec![0], 0_usize) };
445	system.push(eq0);
446	let eq1 = unsafe { Modulo2Equation::from_parts(vec![0], 0_usize) };
447	system.push(eq1);
448	let solution = system.lazy_gaussian_elimination()?;
449	assert!(system.check(&solution));
450	Ok(())
451	}
452
453	#[test]
454	fn test_small_system() -> anyhow::Result<()> {
455	let mut system = Modulo2System::<usize>::new(11);
456	let mut eq = unsafe { Modulo2Equation::from_parts(vec![1, 4, 10], 0) };
457	system.push(eq);
458	eq = unsafe { Modulo2Equation::from_parts(vec![1, 4, 9], 2) };
459	system.push(eq);
460	eq = unsafe { Modulo2Equation::from_parts(vec![0, 6, 8], 0) };
461	system.push(eq);
462	eq = unsafe { Modulo2Equation::from_parts(vec![0, 6, 9], 1) };
463	system.push(eq);
464	eq = unsafe { Modulo2Equation::from_parts(vec![2, 4, 8], 2) };
465	system.push(eq);
466	eq = unsafe { Modulo2Equation::from_parts(vec![2, 6, 10], 0) };
467	system.push(eq);
468
469	let solution = system.lazy_gaussian_elimination()?;
470	assert!(system.check(&solution));
471	Ok(())
472	}
473
474	#[test]
475	fn test_var_to_vec_builder() {
476	let system = unsafe {
477	Modulo2System::from_parts(
478	11,
479	vec![
480	Modulo2Equation::from_parts(vec![1, 4, 10], 0_usize),
481	Modulo2Equation::from_parts(vec![1, 4, 9], 1),
482	Modulo2Equation::from_parts(vec![0, 6, 8], 2),
483	Modulo2Equation::from_parts(vec![0, 6, 9], 3),
484	Modulo2Equation::from_parts(vec![2, 4, 8], 4),
485	Modulo2Equation::from_parts(vec![2, 6, 10], 5),
486	],
487	)
488	};
489	let mut backing: Vec<usize> = vec![];
490	let (var_to_eqs, _weight, _priority) = system.setup(&mut backing);
491	let expected_res = vec![
492	vec![2, 3],
493	vec![0, 1],
494	vec![4, 5],
495	vec![],
496	vec![0, 1, 4],
497	vec![],
498	vec![2, 3, 5],
499	vec![],
500	vec![2, 4],
501	vec![1, 3],
502	vec![0, 5],
503	];
504	var_to_eqs
505	.iter()
506	.zip(expected_res.iter())
507	.for_each(\|(v, e)\| v.iter().zip(e.iter()).for_each(\|(x, y)\| assert_eq!(x, y)));
508	}
509	}

vigna / sux-rs / 21985226363

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