algebraeon_rings/matrix/smith_normal_form.rs
1use super::*;
2
3impl<RS: BezoutDomainSignature, RSB: BorrowedStructure<RS>> MatrixStructure<RS, RSB> {
4 //return (u, s, v, k) such that self = usv and s is in smith normal form (with diagonal entries their favorite associates) and u, v are invertible and k is the number of non-zero elements in the diagonal of s
5 pub fn smith_algorithm(
6 &self,
7 mut m: Matrix<RS::Set>,
8 ) -> (Matrix<RS::Set>, Matrix<RS::Set>, Matrix<RS::Set>, usize) {
9 let mut u = self.ident(m.rows());
10 let mut v = self.ident(m.cols());
11
12 let mut n = 0;
13 'inductive_loop: while n < m.rows() && n < m.cols() {
14 //search for a non-zero element to make the new starting point for (n, n)
15 //having a non-zero element is necessary later in the algorithm
16 'search_for_nonzero_element: {
17 if self.ring().equal(m.at(n, n).unwrap(), &self.ring().zero()) {
18 //search the first row to start with
19 for c in n + 1..m.cols() {
20 if !self.ring().equal(m.at(n, c).unwrap(), &self.ring().zero()) {
21 //swap column n and column c
22 let col_opp = ElementaryOpp::new_col_opp(
23 self.ring().clone(),
24 ElementaryOppType::Swap(n, c),
25 );
26 col_opp.apply(&mut m);
27 col_opp.apply(&mut v);
28
29 break 'search_for_nonzero_element;
30 }
31 }
32 //search all the rows below row n
33 for r in n + 1..m.rows() {
34 for c in n..m.cols() {
35 if !self.ring().equal(m.at(r, c).unwrap(), &self.ring().zero()) {
36 //swap column n and column c
37 let col_opp = ElementaryOpp::new_col_opp(
38 self.ring().clone(),
39 ElementaryOppType::Swap(n, c),
40 );
41 col_opp.apply(&mut m);
42 col_opp.apply(&mut v);
43
44 //swap row n and row r
45 let row_opp = ElementaryOpp::new_col_opp(
46 self.ring().clone(),
47 ElementaryOppType::Swap(n, r),
48 );
49 row_opp.apply(&mut m);
50 row_opp.apply(&mut u);
51
52 break 'search_for_nonzero_element;
53 }
54 }
55 }
56 //everything remaining is zero. The algorithm can terminate here
57 break 'inductive_loop;
58 }
59 }
60
61 //turn (n, n) into its favorite associate
62 let (unit, _assoc) = self.ring().factor_fav_assoc(m.at(n, n).unwrap());
63 let row_opp = ElementaryOpp::new_row_opp(
64 self.ring().clone(),
65 ElementaryOppType::UnitMul {
66 row: n,
67 unit: self.ring().try_reciprocal(&unit).unwrap(),
68 },
69 );
70 row_opp.apply(&mut m);
71 row_opp.apply(&mut u);
72
73 let mut first = true;
74 let mut all_divideisible;
75 'zero_first_row_and_column_loop: loop {
76 //replace the first row (a0, a1, ..., ak) with (gcd, 0, ..., 0). Might mess up the first column in the process
77 all_divideisible = true;
78 for c in n + 1..m.cols() {
79 let a = m.at(n, n).unwrap();
80 let b = m.at(n, c).unwrap();
81
82 if self.ring().is_zero(a) {
83 //swap a and b
84 //a=0 so this does have the effect of (a, b) -> (gcd(a, b), 0)
85 let col_opp = ElementaryOpp::new_col_opp(
86 self.ring().clone(),
87 ElementaryOppType::Swap(n, c),
88 );
89 col_opp.apply(&mut m);
90 col_opp.apply(&mut v);
91 } else {
92 match self.ring().try_divide(b, a) {
93 Some(q) => {
94 //b is a multiple of a
95 //replace (a, b) with (a, 0) by subtracting a multiple of a from b
96 let col_opp = ElementaryOpp::new_col_opp(
97 self.ring().clone(),
98 ElementaryOppType::AddRowMul {
99 i: c,
100 j: n,
101 x: self.ring().neg(&q),
102 },
103 );
104 col_opp.apply(&mut m);
105 col_opp.apply(&mut v);
106 }
107 None => {
108 all_divideisible = false;
109 //b is not a multiple of a
110 //replace (a, b) with (gcd, 0)
111 let (d, x, y) = self.ring().xgcd(a, b);
112 debug_assert!(
113 self.ring().equal(
114 &self
115 .ring()
116 .add(&self.ring().mul(&x, a), &self.ring().mul(&y, b)),
117 &d
118 )
119 );
120 let col_opp = ElementaryOpp::new_col_opp(
121 self.ring().clone(),
122 ElementaryOppType::TwoInv {
123 i: n,
124 j: c,
125 a: x,
126 b: y,
127 c: self.ring().neg(&self.ring().try_divide(b, &d).unwrap()),
128 d: self.ring().try_divide(a, &d).unwrap(),
129 },
130 );
131 col_opp.apply(&mut m);
132 col_opp.apply(&mut v);
133 }
134 }
135 }
136 }
137 if all_divideisible && !first {
138 break 'zero_first_row_and_column_loop;
139 }
140 first = false;
141
142 //replace the first column (a0, a1, ..., ak) with (gcd, 0, ..., 0). Might mess up the first row in the process
143 all_divideisible = true;
144 for r in n + 1..m.rows() {
145 let a = m.at(n, n).unwrap();
146 let b = m.at(r, n).unwrap();
147 if self.ring().is_zero(a) {
148 //swap a and b
149 //a=0 so this does have the effect of (a, b) -> (gcd(a, b), 0)
150 let col_opp = ElementaryOpp::new_row_opp(
151 self.ring().clone(),
152 ElementaryOppType::Swap(n, r),
153 );
154 col_opp.apply(&mut m);
155 col_opp.apply(&mut u);
156 } else {
157 match self.ring().try_divide(b, a) {
158 Some(q) => {
159 //b is a multiple of a
160 //replace (a, b) with (a, 0) by subtracting a multiple of a from b
161 let col_opp = ElementaryOpp::new_row_opp(
162 self.ring().clone(),
163 ElementaryOppType::AddRowMul {
164 i: r,
165 j: n,
166 x: self.ring().neg(&q),
167 },
168 );
169 col_opp.apply(&mut m);
170 col_opp.apply(&mut u);
171 }
172 None => {
173 all_divideisible = false;
174 //b is not a multiple of a
175 //replace (a, b) with (gcd, 0)
176 let (d, x, y) = self.ring().xgcd(a, b);
177 debug_assert!(
178 self.ring().equal(
179 &self
180 .ring()
181 .add(&self.ring().mul(&x, a), &self.ring().mul(&y, b)),
182 &d
183 )
184 );
185 let row_opp = ElementaryOpp::new_row_opp(
186 self.ring().clone(),
187 ElementaryOppType::TwoInv {
188 i: n,
189 j: r,
190 a: x,
191 b: y,
192 c: self.ring().neg(&self.ring().try_divide(b, &d).unwrap()),
193 d: self.ring().try_divide(a, &d).unwrap(),
194 },
195 );
196 row_opp.apply(&mut m);
197 row_opp.apply(&mut u);
198 }
199 }
200 }
201 }
202 if all_divideisible {
203 break 'zero_first_row_and_column_loop;
204 }
205 }
206 //now the first row and the first column are all zero except the top left element at (n, n) which is non-zero
207 debug_assert!(!self.ring().equal(m.at(n, n).unwrap(), &self.ring().zero()));
208 //some more fiddling is needed now to make sure the top left element divides everything else
209 for r in n + 1..m.rows() {
210 //row(n) = row(n) + row(r)
211 let row_opp = ElementaryOpp::new_row_opp(
212 self.ring().clone(),
213 ElementaryOppType::AddRowMul {
214 i: n,
215 j: r,
216 x: self.ring().one(),
217 },
218 );
219 row_opp.apply(&mut m);
220 row_opp.apply(&mut u);
221
222 //row(n) goes from (g, a1, a2, ..., an) to (gcd, 0, 0, ..., 0) by applying col operations
223 for c in n + 1..m.cols() {
224 let a = m.at(n, n).unwrap();
225 let b = m.at(n, c).unwrap();
226 // println!("a={:?} b={:?}", a, b);
227 //if a=0 and b=0 then there is nothing to do and the following step would fail when dividing by g=0
228 if !self.ring().is_zero(a) || !self.ring().is_zero(b) {
229 //b might not be a multiple of a
230 //replace (a, b) with (gcd, 0) to fix this
231 let (g, x, y) = self.ring().xgcd(a, b);
232 debug_assert!(
233 self.ring().equal(
234 &self
235 .ring()
236 .add(&self.ring().mul(&x, a), &self.ring().mul(&y, b)),
237 &g
238 )
239 );
240 let col_opp = ElementaryOpp::new_col_opp(
241 self.ring().clone(),
242 ElementaryOppType::TwoInv {
243 i: n,
244 j: c,
245 a: x,
246 b: y,
247 c: self.ring().neg(&self.ring().try_divide(b, &g).unwrap()),
248 d: self.ring().try_divide(a, &g).unwrap(),
249 },
250 );
251 col_opp.apply(&mut m);
252 col_opp.apply(&mut v);
253 }
254 }
255 }
256
257 //fix the first column
258 for fix_r in n + 1..m.rows() {
259 let a = m.at(n, n).unwrap();
260 let b = m.at(fix_r, n).unwrap();
261 let q = self.ring().try_divide(b, a).unwrap();
262 let col_opp = ElementaryOpp::new_row_opp(
263 self.ring().clone(),
264 ElementaryOppType::AddRowMul {
265 i: fix_r,
266 j: n,
267 x: self.ring().neg(&q),
268 },
269 );
270 col_opp.apply(&mut m);
271 col_opp.apply(&mut u);
272 }
273
274 if self.ring().equal(m.at(n, n).unwrap(), &self.ring().zero()) {
275 //the bottom right submatrix is all zero
276 break 'inductive_loop;
277 }
278 n += 1;
279 }
280
281 (u, m, v, n)
282 }
283}
284
285impl<R: MetaType> Matrix<R>
286where
287 R::Signature: BezoutDomainSignature,
288{
289 pub fn smith_algorithm(&self) -> (Self, Self, Self, usize) {
290 Self::structure().smith_algorithm(self.clone())
291 }
292}
293
294#[cfg(test)]
295mod tests {
296 use super::*;
297
298 #[test]
299 fn test_smith_algorithm() {
300 {
301 let a = Matrix::<Integer>::from_rows(vec![
302 vec![Integer::from(2), Integer::from(4), Integer::from(4)],
303 vec![Integer::from(-6), Integer::from(6), Integer::from(12)],
304 vec![Integer::from(10), Integer::from(4), Integer::from(16)],
305 ]);
306 let (u, s, v, k) = a.clone().smith_algorithm();
307 assert_eq!(s, Matrix::mul(&Matrix::mul(&u, &a).unwrap(), &v).unwrap());
308 assert_eq!(k, 3);
309 assert_eq!(
310 s,
311 Matrix::from_rows(vec![
312 vec![Integer::from(2), Integer::from(0), Integer::from(0)],
313 vec![Integer::from(0), Integer::from(2), Integer::from(0)],
314 vec![Integer::from(0), Integer::from(0), Integer::from(156)]
315 ])
316 );
317
318 let a = Matrix::<Integer>::from_rows(vec![
319 vec![
320 Integer::from(-6),
321 Integer::from(111),
322 Integer::from(-36),
323 Integer::from(6),
324 ],
325 vec![
326 Integer::from(5),
327 Integer::from(-672),
328 Integer::from(210),
329 Integer::from(74),
330 ],
331 vec![
332 Integer::from(0),
333 Integer::from(-255),
334 Integer::from(81),
335 Integer::from(24),
336 ],
337 vec![
338 Integer::from(-7),
339 Integer::from(255),
340 Integer::from(-81),
341 Integer::from(-10),
342 ],
343 ]);
344 let (u, s, v, k) = a.clone().smith_algorithm();
345 assert_eq!(s, Matrix::mul(&Matrix::mul(&u, &a).unwrap(), &v).unwrap());
346 assert_eq!(k, 3);
347 assert_eq!(
348 s,
349 Matrix::from_rows(vec![
350 vec![
351 Integer::from(1),
352 Integer::from(0),
353 Integer::from(0),
354 Integer::from(0)
355 ],
356 vec![
357 Integer::from(0),
358 Integer::from(3),
359 Integer::from(0),
360 Integer::from(0)
361 ],
362 vec![
363 Integer::from(0),
364 Integer::from(0),
365 Integer::from(21),
366 Integer::from(0)
367 ],
368 vec![
369 Integer::from(0),
370 Integer::from(0),
371 Integer::from(0),
372 Integer::from(0)
373 ]
374 ])
375 );
376 }
377
378 {
379 //used to cause a divide by zero
380 let a = Matrix::<Integer>::from_rows(vec![
381 vec![
382 Integer::from(0),
383 Integer::from(0),
384 Integer::from(0),
385 Integer::from(0),
386 ],
387 vec![
388 Integer::from(0),
389 Integer::from(0),
390 Integer::from(0),
391 Integer::from(0),
392 ],
393 vec![
394 Integer::from(0),
395 Integer::from(0),
396 Integer::from(0),
397 Integer::from(0),
398 ],
399 vec![
400 Integer::from(0),
401 Integer::from(0),
402 Integer::from(0),
403 Integer::from(1),
404 ],
405 ]);
406 let (_u, _s, _v, k) = a.clone().smith_algorithm();
407 assert_eq!(k, 1);
408 }
409 }
410}