rustiq_core/routines/
f2_linalg.rs

1use crate::structures::{CliffordCircuit, CliffordGate};
2pub type Matrix = Vec<Vec<bool>>;
3
4pub fn xor_vec(a: &mut [bool], b: &[bool]) {
5    for i in 0..a.len() {
6        a[i] ^= b[i];
7    }
8}
9
10pub fn rowop(table: &mut Matrix, i: usize, j: usize) {
11    if !table.is_empty() {
12        for k in 0..table.first().unwrap().len() {
13            table[j][k] ^= table[i][k];
14        }
15    }
16}
17
18pub fn colop(table: &mut Matrix, i: usize, j: usize) {
19    for row in table.iter_mut() {
20        row[j] ^= row[i];
21    }
22}
23
24pub fn row_echelon(table: &mut Matrix, k: usize) {
25    let mut rank = 0;
26    for i in 0..table.first().unwrap().len() {
27        let mut pivot = None;
28        for (j, row) in table.iter().enumerate().take(k).skip(rank) {
29            if row[i] {
30                pivot = Some(j);
31                break;
32            }
33        }
34        if let Some(pivot) = pivot {
35            if pivot != rank {
36                rowop(table, pivot, rank);
37            }
38            for j in 0..table.len() {
39                if table[j][i] && j != rank {
40                    rowop(table, rank, j);
41                }
42            }
43            rank += 1;
44        }
45    }
46}
47
48pub fn diagonalize(table: &mut Matrix, friend: &mut Matrix, rank: usize) {
49    let n = table.first().unwrap().len();
50    for i in 0..rank {
51        let mut pivot = None;
52        for j in i..n {
53            if table[i][j] {
54                pivot = Some(j);
55                break;
56            }
57        }
58        if let Some(pivot) = pivot {
59            if pivot != i {
60                colop(table, pivot, i);
61                colop(friend, pivot, i);
62            }
63            for j in 0..n {
64                if table[i][j] && j != i {
65                    colop(table, i, j);
66                    colop(friend, i, j);
67                }
68            }
69        } else {
70            panic!("This is not gooood!");
71        }
72    }
73}
74
75pub fn f2_rank(table: &Matrix) -> usize {
76    let mut rank = 0;
77    let mut table = table.clone();
78    let nrows = table.len();
79    let ncols = table.first().unwrap().len();
80    for i in 0..ncols {
81        let mut pivot = None;
82        for (j, row) in table.iter().enumerate().take(nrows).skip(rank) {
83            if row[i] {
84                pivot = Some(j);
85                break;
86            }
87        }
88        if let Some(pivot) = pivot {
89            if pivot != rank {
90                rowop(&mut table, pivot, rank);
91            }
92            for j in rank + 1..nrows {
93                if table[j][i] {
94                    rowop(&mut table, rank, j);
95                }
96            }
97            rank += 1;
98        }
99    }
100
101    rank
102}
103
104pub fn inverse_f2(table: &Matrix) -> Matrix {
105    let n = table.len();
106    let mut table = table.clone();
107    let mut friend = vec![vec![false; n]; n];
108    for (i, row) in friend.iter_mut().enumerate().take(n) {
109        row[i] = true;
110    }
111    for i in 0..n {
112        let mut pivot = None;
113        for (j, row) in table.iter().enumerate().take(n).skip(i) {
114            if row[i] {
115                pivot = Some(j);
116                break;
117            }
118        }
119        if let Some(pivot) = pivot {
120            table.swap(i, pivot);
121            friend.swap(i, pivot);
122            for j in 0..n {
123                if j != i && table[j][i] {
124                    rowop(&mut table, i, j);
125                    rowop(&mut friend, i, j)
126                }
127            }
128        }
129    }
130    friend
131}
132
133pub fn mult_f2(a: &Matrix, b: &Matrix) -> Matrix {
134    let (k, l, m) = (a.len(), a.first().unwrap().len(), b.first().unwrap().len());
135    assert_eq!(l, b.len());
136    let mut result = vec![vec![false; m]; k];
137    for i in 0..k {
138        for j in 0..m {
139            for (y, b_row) in b.iter().enumerate().take(l) {
140                result[i][j] ^= a[i][y] & b_row[j];
141            }
142        }
143    }
144    result
145}
146
147pub fn transpose(table: &Matrix) -> Matrix {
148    let n = table.len();
149    let k = table.first().unwrap().len();
150    let mut result = vec![vec![false; n]; k];
151    for (i, row) in result.iter_mut().enumerate().take(k) {
152        for (j, table_row) in table.iter().enumerate().take(n) {
153            row[j] = table_row[i];
154        }
155    }
156    result
157}
158
159pub fn plu_facto(table: &Matrix) -> (Matrix, Matrix, Matrix) {
160    let n = table.len();
161    let k = table.first().unwrap().len();
162    let mut u = table.clone();
163    let mut l = vec![vec![false; k]; n];
164    let mut p = vec![vec![false; n]; n];
165    for (i, row) in p.iter_mut().enumerate().take(n) {
166        row[i] = true;
167    }
168    for i in 0..k {
169        let mut pivot = None;
170        for (j, row) in u.iter().enumerate().take(n).skip(i) {
171            if row[i] {
172                pivot = Some(j);
173                break;
174            }
175        }
176
177        if let Some(pivot) = pivot {
178            if i != pivot {
179                u.swap(i, pivot);
180                p.swap(i, pivot);
181                l.swap(i, pivot);
182            }
183            l[i][i] = true;
184            for j in i + 1..n {
185                if u[j][i] {
186                    l[j][i] = true;
187                    rowop(&mut u, i, j);
188                }
189            }
190        }
191    }
192
193    (p, l, u.drain(..k).collect())
194}
195
196pub fn lu_facto(table: &Matrix) -> (Matrix, Matrix, Matrix, CliffordCircuit) {
197    let mut table = table.clone();
198    let (c, ops) = non_zero_leading_principal_minors(&table);
199    let mut output = CliffordCircuit::new(table.len());
200    for (a, b) in ops.iter() {
201        rowop(&mut table, *a, *b);
202        output.gates.push(CliffordGate::CNOT(*b, *a));
203    }
204    let (p, l, u) = plu_facto(&table);
205    for (i, row) in p.iter().enumerate().take(table.len()) {
206        assert!(row[i]);
207    }
208    (l, u, c, output)
209}
210
211fn f2_rank_square(matrix: &Matrix) -> usize {
212    let matrix: Matrix = matrix
213        .iter()
214        .map(|v| v.clone().into_iter().take(matrix.len()).collect())
215        .collect();
216    f2_rank(&matrix)
217}
218pub fn print_matrix(matrix: &Matrix) {
219    for row in matrix.iter() {
220        for elem in row.iter() {
221            if *elem {
222                print!("1");
223            } else {
224                print!("0");
225            }
226        }
227        println!();
228    }
229}
230/// Finds a sequence of row operations that makes sure that all the leading principal
231/// minors are 1
232/// In other words, makes sure that matrix[:i, :i] is invertible
233pub fn non_zero_leading_principal_minors(matrix: &Matrix) -> (Matrix, Vec<(usize, usize)>) {
234    let mut piece: Matrix = Vec::new();
235    let mut moves = Vec::new();
236    for i in 0..f2_rank(matrix) {
237        piece.push(matrix[i].clone());
238        let mut current_rk = f2_rank_square(&piece);
239        while current_rk == i {
240            for (k, row) in matrix.iter().enumerate().skip(i + 1) {
241                xor_vec(&mut piece[i], row);
242                let new_rank = f2_rank_square(&piece);
243                if new_rank == i + 1 {
244                    moves.push((k, i));
245                    current_rk = new_rank;
246                    break;
247                } else {
248                    xor_vec(&mut piece[i], row);
249                }
250            }
251            if current_rk == i + 1 {
252                break;
253            }
254        }
255    }
256    let mut c = vec![vec![false; matrix.len()]; matrix.len()];
257    for (i, row) in c.iter_mut().enumerate() {
258        row[i] = true;
259    }
260    for (a, b) in moves.iter() {
261        rowop(&mut c, *a, *b);
262    }
263    let m = mult_f2(&c, matrix);
264    for i in 0..piece.len() {
265        assert_eq!(piece[i], m[i]);
266    }
267    (c, moves)
268}
269
270pub fn count_ones(matrix: &Matrix) -> usize {
271    matrix
272        .iter()
273        .map(|row| row.iter().filter(|x| **x).count())
274        .sum()
275}
276
277pub fn count_ones_except_diag(matrix: &Matrix) -> usize {
278    let mut count = 0;
279    for (i, row) in matrix.iter().enumerate() {
280        for (j, item) in row.iter().enumerate() {
281            if i != j && *item {
282                count += 1;
283            }
284        }
285    }
286    count
287}
288#[cfg(test)]
289mod tests {
290    use super::*;
291    use rand::Rng;
292
293    fn random_skinny(n: usize, m: usize) -> Matrix {
294        assert!(m <= n);
295        let mut rng = rand::thread_rng();
296        let mut matrix = vec![vec![false; m]; n];
297        for (i, row) in matrix.iter_mut().take(m).enumerate() {
298            row[i] = true;
299        }
300        for _ in 0..n * n {
301            let i = rng.gen_range(0..n);
302            let j = rng.gen_range(0..n);
303            if i != j {
304                rowop(&mut matrix, i, j);
305            }
306        }
307        matrix
308    }
309    #[test]
310    fn test_plu_facto() {
311        for _ in 0..10 {
312            let n = 40;
313            let m = 20;
314            let matrix = random_skinny(n, m);
315            let (p, l, u) = plu_facto(&matrix);
316
317            // Check if P is a permutation matrix
318            assert_eq!(
319                p.iter()
320                    .map(|r| r.iter().filter(|&&v| v).count())
321                    .collect::<Vec<_>>(),
322                vec![1; matrix.len()]
323            );
324            let mut perm = p
325                .iter()
326                .map(|r| r.iter().position(|&v| v).unwrap())
327                .collect::<Vec<_>>();
328            perm.sort();
329            assert_eq!(perm, (0..n).collect::<Vec<_>>());
330
331            for (i, row) in l.iter().enumerate() {
332                if i < row.len() - 1 {
333                    assert!(
334                        row[i + 1..].iter().all(|&v| !v),
335                        "L should be lower triangular"
336                    );
337                }
338            }
339
340            for (i, row) in u.iter().enumerate() {
341                assert!(row[..i].iter().all(|&v| !v), "U should be lower triangular");
342            }
343
344            // Check if PA = LU
345            let pa = mult_f2(&p, &matrix);
346            let lu = mult_f2(&l, &u);
347            assert_eq!(pa, lu, "PA should be equal to LU");
348        }
349    }
350
351    #[test]
352    fn test_lu_facto() {
353        for _ in 0..10 {
354            let n = 40;
355            let m = 22;
356            let matrix = random_skinny(n, m);
357            let (l, u, c, _) = lu_facto(&matrix);
358
359            for (i, row) in l.iter().enumerate() {
360                if i < row.len() - 1 {
361                    assert!(
362                        row[i + 1..].iter().all(|&v| !v),
363                        "L should be lower triangular"
364                    );
365                }
366            }
367
368            for (i, row) in u.iter().enumerate() {
369                assert!(row[..i].iter().all(|&v| !v), "U should be lower triangular");
370            }
371            let lu = mult_f2(&l, &u);
372            let clu = mult_f2(&inverse_f2(&c), &lu);
373            assert_eq!(clu, matrix, "C^-1 LU should be equal to A");
374        }
375    }
376
377    fn random_invertible(n: usize) -> Matrix {
378        let mut rng = rand::thread_rng();
379        let mut matrix = vec![vec![false; n]; n];
380        for (i, row) in matrix.iter_mut().enumerate() {
381            row[i] = true;
382        }
383        for _ in 0..n * n {
384            let i = rng.gen_range(0..n);
385            let j = rng.gen_range(0..n);
386            if i != j {
387                rowop(&mut matrix, i, j);
388            }
389        }
390        matrix
391    }
392
393    #[test]
394    fn test_leading_full_rank() {
395        for _ in 0..40 {
396            // Test data
397            let n = 40;
398            let mut test_matrix = random_invertible(n);
399            let (_, moves) = non_zero_leading_principal_minors(&test_matrix);
400            for (a, b) in moves {
401                rowop(&mut test_matrix, a, b);
402            }
403            for i in 1..n {
404                let piece: Matrix = test_matrix.clone().into_iter().take(i).collect();
405                assert_eq!(f2_rank_square(&piece), i);
406            }
407        }
408    }
409}
rustiq_core/routines/f2_linalg.rs

rustiq_core/routines/
f2_linalg.rs
