feanor_math/algorithms/linsolve/
gauss.rs

1use crate::divisibility::{DivisibilityRing, DivisibilityRingStore};
2use crate::local::{PrincipalLocalRing, PrincipalLocalRingStore};
3use crate::matrix::{AsPointerToSlice, SubmatrixMut, TransposableSubmatrixMut};
4use crate::ring::*;
5
6///
7/// Computes the largest square submatrix of `A` that has nonzero determinant. In particular,
8/// if the ring is a field, this computes the rank of `A`.
9/// 
10/// The result are two sorted lists `rows_idxs` and `col_idxs` both of same length `k` such 
11/// that the submatrix is of size `k x k` and has entries `A[row_idxs[i], col_idxs[j]]` for `i, j < k`.
12/// 
13/// This function changes the content of `A` in an unspecified way.
14/// 
15/// While this functionality seems to be well-defined for arbitrary rings, the current algorithm
16/// (based on gaussian elimination) only works for (possibly non-integral) valuation rings, since it
17/// relies on divisibility inducing a total order.
18/// 
19/// # Example
20/// ```rust
21/// # use feanor_math::assert_el_eq;
22/// # use feanor_math::ring::*;
23/// # use feanor_math::primitive_int::*;
24/// # use feanor_math::rings::rational::*;
25/// # use feanor_math::homomorphism::Homomorphism;
26/// # use feanor_math::rings::local::*;
27/// # use feanor_math::algorithms;
28/// # use feanor_math::matrix::SubmatrixMut;
29/// let ZZ = StaticRing::<i64>::RING;
30/// let QQ = RationalField::new(ZZ);
31/// // due to a suboptimal interface of `RationalField`, it currently doesn't implement `PrincipalLocalRing`;
32/// // Thus we must use a wrapper
33/// let QQ = AsLocalPIR::from_field(QQ);
34/// let hom = QQ.can_hom(&ZZ).unwrap();
35/// let (row_idxs, col_idxs) = algorithms::linsolve::gauss::largest_nonzero_minor(SubmatrixMut::<Vec<_>, _>::from_2d(&mut [vec![hom.map(0), hom.map(1)], vec![hom.map(0), hom.map(0)]]), QQ);
36/// let rank = row_idxs.len();
37/// assert_eq!(rank, col_idxs.len());
38/// assert_eq!(vec![0], row_idxs);
39/// assert_eq!(vec![1], col_idxs);
40/// ```
41/// 
42#[stability::unstable(feature = "enable")]
43pub fn largest_nonzero_minor<R, V>(A: SubmatrixMut<V, El<R>>, ring: R) -> (Vec<usize>, Vec<usize>)
44    where R: RingStore + Copy,
45        R::Type: PrincipalLocalRing,
46        V: AsPointerToSlice<El<R>>
47{
48    assert!(ring.is_noetherian());
49    assert!(ring.is_commutative());
50    let n = A.row_count();
51    let m = A.col_count();
52
53    fn largest_nonzero_minor_impl<R, V, const T: bool>(mut A: TransposableSubmatrixMut<V, El<R>, T>, ring: R) -> (Vec<usize>, Vec<usize>)
54        where R: RingStore + Copy,
55            R::Type: PrincipalLocalRing,
56            V: AsPointerToSlice<El<R>>
57    {
58        let n = A.row_count();
59        let m = A.col_count();
60        assert!(n >= m);
61        let mut row_perm: Vec<usize> = (0..n).collect();
62        let mut col_perm: Vec<usize> = (0..m).collect();
63
64        fn swap_rows<V, E, const TRANSPOSED: bool>(mut A: TransposableSubmatrixMut<V, E, TRANSPOSED>, i: usize, j: usize)
65            where V: AsPointerToSlice<E>
66        {
67            if i == j {
68                return;
69            }
70            let m = A.col_count();
71            let (mut fst, mut snd) = A.reborrow().split_rows(i..(i + 1), j..(j + 1));
72            for l in 0..m {
73                std::mem::swap(fst.at_mut(0, l), snd.at_mut(0, l));
74            }
75        }
76
77        fn swap_cols<V, E, const TRANSPOSED: bool>(mut A: TransposableSubmatrixMut<V, E, TRANSPOSED>, i: usize, j: usize)
78            where V: AsPointerToSlice<E>
79        {
80            if i == j {
81                return;
82            }
83            let n = A.row_count();
84            let (mut fst, mut snd) = A.reborrow().split_cols(i..(i + 1), j..(j + 1));
85            for k in 0..n {
86                std::mem::swap(fst.at_mut(k, 0), snd.at_mut(k, 0));
87            }
88        }
89
90        fn elim_row<V, R, const TRANSPOSED: bool>(mut A: TransposableSubmatrixMut<V, El<R>, TRANSPOSED>, pivot: usize, src: usize, dst: usize, ring: R)
91            where R: RingStore,
92                R::Type: DivisibilityRing,
93                V: AsPointerToSlice<El<R>>
94        {
95            let m = A.col_count();
96            let (src, mut dst) = A.reborrow().split_rows(src..(src + 1), dst..(dst + 1));
97            let factor = ring.checked_div(dst.at(0, pivot), src.at(0, pivot)).unwrap();
98            for l in 0..m {
99                ring.sub_assign(dst.at_mut(0, l), ring.mul_ref(&factor, src.at(0, l)));
100            }
101        }
102
103        let mut i = 0;
104        while i < m {
105            let (pivot_i, pivot_j) = (i..n).flat_map(|k| (i..m).map(move |l| (k, l))).min_by_key(|(k, l)| ring.valuation(A.at(*k, *l)).unwrap_or(usize::MAX)).unwrap();
106            if ring.valuation(A.at(pivot_i, pivot_j)).is_none() {
107                break;
108            }
109            row_perm.swap(i, pivot_i);
110            col_perm.swap(i, pivot_j);
111            swap_rows(A.reborrow(), i, pivot_i);
112            swap_cols(A.reborrow(), i, pivot_j);
113            for k in (i + 1)..n {
114                elim_row(A.reborrow(), i, i, k, &ring);
115            }
116            i += 1;
117        }
118
119        let mut k = 0;
120        let mut current = ring.clone_el(A.at(0, 0));
121        while (k + 1) < m && !ring.is_zero(&current) {
122            k += 1;
123            ring.mul_assign_ref(&mut current, A.at(k, k));
124        }
125        if !ring.is_zero(&current) {
126            k += 1;
127        }
128        let mut row_result = row_perm;
129        _ = row_result.drain(k..);
130        let mut col_result = col_perm;
131        _ = col_result.drain(k..);
132        row_result.sort_unstable();
133        col_result.sort_unstable();
134        return (row_result, col_result);
135    }
136
137    if n >= m {
138        return largest_nonzero_minor_impl(TransposableSubmatrixMut::from(A), ring);
139    } else {
140        let (col_res, row_res) = largest_nonzero_minor_impl(TransposableSubmatrixMut::from(A).transpose(), ring);
141        return (row_res, col_res);
142    }
143}
144
145#[cfg(test)]
146use crate::rings::local::AsLocalPIR;
147#[cfg(test)]
148use crate::rings::zn::zn_64::Zn;
149#[cfg(test)]
150use crate::rings::zn::zn_static::Fp;
151#[cfg(test)]
152use crate::homomorphism::Homomorphism;
153
154#[test]
155fn test_largest_nonzero_minor_field() {
156    let field = Fp::<17>::RING;
157
158    let mut matrix = [vec![1, 0], vec![1, 0]];
159    let (rows, cols) = largest_nonzero_minor(SubmatrixMut::<Vec<_>, _>::from_2d(&mut matrix), field);
160    assert_eq!(1, rows.len());
161    assert_eq!(vec![0], cols);
162
163    let mut matrix = [vec![0, 0], vec![0, 1]];
164    let (rows, cols) = largest_nonzero_minor(SubmatrixMut::<Vec<_>, _>::from_2d(&mut matrix), field);
165    assert_eq!(vec![1], rows);
166    assert_eq!(vec![1], cols);
167
168    let mut matrix = [vec![1, 2, 3], vec![1, 2, 3], vec![2, 3, 4]];
169    let (rows, cols) = largest_nonzero_minor(SubmatrixMut::<Vec<_>, _>::from_2d(&mut matrix), field);
170    assert!(rows == vec![0, 2] || rows == vec![1, 2]);
171    assert_eq!(2, cols.len());
172
173    let mut matrix = [
174        vec![15,  3,  9, 15,  9,],
175        vec![10,  6,  7,  3,  9,],
176        vec![ 2, 14, 14,  8,  6,],
177        vec![12, 16,  8,  6, 16,],
178        vec![15,  4, 14,  1, 11,]
179    ];
180    let (rows, cols) = largest_nonzero_minor(SubmatrixMut::<Vec<_>, _>::from_2d(&mut matrix), field);
181    assert_eq!(3, rows.len());
182    assert_eq!(3, cols.len());
183}
184
185#[test]
186fn test_largest_nonzero_minor_localpir() {
187    let ring = AsLocalPIR::from_zn(Zn::new(8)).unwrap();
188    let i = |x: i32| ring.int_hom().map(x);
189
190    let mut matrix = [vec![i(0), i(0)], vec![i(0), i(1)]];
191    let (rows, cols) = largest_nonzero_minor(SubmatrixMut::<Vec<_>, _>::from_2d(&mut matrix), ring);
192    assert_eq!(vec![1], rows);
193    assert_eq!(vec![1], cols);
194
195    let mut matrix = [vec![i(4), i(0), i(0)], vec![i(0), i(0), i(2)], vec![i(0), i(1), i(0)]];
196    let (rows, cols) = largest_nonzero_minor(SubmatrixMut::<Vec<_>, _>::from_2d(&mut matrix), ring);
197    assert!((&vec![0, 2], &vec![0, 1]) == (&rows, &cols) || (&vec![1, 2], &vec![1, 2]) == (&rows, &cols));
198}
feanor_math/algorithms/linsolve/gauss.rs

feanor_math/algorithms/linsolve/
gauss.rs
