17841719366

Committed 18 Sep 2025 09:21PM UTC coverage: 64.717% (-2.2%) from 66.912%

Build # 17841719366

Build Type

push

github

Committed by

vigna

Commit Message

Completed delegations for custom iterators; fixed infinite recursion when delegating traits to MemCase

Run Details

0 of 7 new or added lines in 3 files covered. (0.0%)

909 existing lines in 21 files now uncovered.

3674 of 5677 relevant lines covered (64.72%)

117564551.11 hits per line

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

75.86

/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
 */

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

use crate::{bit_vec, traits::Word};
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 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() {
                // SAFETY: to add the two equations, multiple references to the vector
                // of equations are needed, one of which is mutable
                let eq_j = unsafe { &*(&equations[j] as *const Modulo2Equation<W>) };
                let eq_i = &mut equations[i];

                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 {
                    // SAFETY: to add the equations, multiple references to the vector
                    // of equations are needed, one of which is mutable
                    let equation = unsafe { &*(&equations[first] as *const Modulo2Equation<W>) };

                    let pivot = equation
                        .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;
                    var_to_eqs[pivot as usize]
                        .as_ref()
                        .iter()
                        .filter(|&&eq_idx| eq_idx != first)
                        .for_each(|&eq| {
                            equations[eq].add(equation);

                            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() {
        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!(solution.is_ok());
        assert!(system.check(&solution.unwrap()));
    }

    #[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() {
        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!(solution.is_ok());
        assert!(system.check(&solution.unwrap()));
    }

    #[test]
    fn test_small_system() {
        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!(solution.is_ok());
        assert!(system.check(&solution.unwrap()));
    }

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

vigna / sux-rs / 17841719366

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