polycore 0.1.6

Fields, prime fields, sparse and dense polynomials, and sparse echelon forms.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100

use polycore::modp_echelon;
use polycore::sample::Rng;
use polycore::{Echelon, Fp, Lead, Modular};

#[test]
fn agrees_with_generic_elimination() {
    for p in [
        2,
        101,
        4_611_686_018_427_387_847,
        18_446_744_073_709_551_557,
    ] {
        for lead in [Lead::Low, Lead::High] {
            let mut rng = Rng::new(19);
            let mut raw = modp_echelon::Echelon::recording(lead, p);
            let mut generic = Echelon::recording(lead);
            for _ in 0..28 {
                // Deliberately include duplicates, zeros, and unordered columns.
                let row: Vec<_> = (0..8)
                    .map(|_| ((rng.next_u64() % 32) as usize, rng.next_u64() % p))
                    .collect();
                let fp: Vec<_> = row.iter().map(|&(j, v)| (j, Fp::new(v, p))).collect();
                assert_eq!(raw.insert(&row), generic.insert(&fp));
            }
            for i in 0..35 {
                let (got, used) = raw.solve(i);
                let (want, dependencies) = generic.solve(i);
                assert_eq!(
                    got,
                    want.iter()
                        .map(|&(j, v)| (j, v.residue_mod(p)))
                        .collect::<Vec<_>>()
                );
                assert_eq!(used, dependencies);
                assert_eq!(raw.uses(i), generic.uses(i));
            }
            let got = raw.nullspace(35);
            let want: Vec<Vec<_>> = generic
                .nullspace(35)
                .iter()
                .map(|r| r.iter().map(|&(j, v)| (j, v.residue_mod(p))).collect())
                .collect();
            assert_eq!(got, want);
        }
    }
}

#[test]
fn rectangular_solve() {
    use polycore::dense::SolveError;
    assert_eq!(
        modp_echelon::solve(&[vec![1, 2], vec![3, 4], vec![4, 6]], &[5, 11, 16], 101),
        Ok(vec![1, 2])
    );
    assert_eq!(
        modp_echelon::solve(&[vec![1, 2], vec![2, 4]], &[5, 11], 101),
        Err(SolveError::Inconsistent)
    );
    assert_eq!(
        modp_echelon::solve(&[vec![1, 2]], &[5], 101),
        Err(SolveError::Deficient)
    );
}

#[test]
fn inversion_rejects_singular_and_multiplies_to_identity() {
    use polycore::modp::{add, mul};
    for p in [
        2,
        101,
        4_611_686_018_427_387_847,
        18_446_744_073_709_551_557,
    ] {
        assert!(modp_echelon::invert(&[vec![1, 1], vec![1, 1]], p).is_none());
        assert!(modp_echelon::invert(&[vec![1, 0], vec![1]], p).is_none());
        let mut rng = Rng::new(13);
        for n in 1..8 {
            let a: Vec<Vec<_>> = (0..n)
                .map(|i| {
                    (0..n)
                        .map(|j| match i.cmp(&j) {
                            std::cmp::Ordering::Equal => 1,
                            std::cmp::Ordering::Less => rng.next_u64() % p,
                            std::cmp::Ordering::Greater => 0,
                        })
                        .collect()
                })
                .collect();
            let b = modp_echelon::invert(&a, p).unwrap();
            for (i, row) in a.iter().enumerate() {
                for (j, _) in b.iter().enumerate() {
                    assert_eq!(
                        (0..n).fold(0, |s, k| add(s, mul(row[k], b[k][j], p), p)),
                        u64::from(i == j)
                    );
                }
            }
        }
    }
}