algebraeon_rings/matrix/
hermite_reduction.rs

1use crate::module::{
2    finitely_free_affine::FinitelyFreeSubmoduleAffineSubset,
3    finitely_free_module::FinitelyFreeModuleStructure,
4    finitely_free_submodule::FinitelyFreeSubmodule,
5};
6
7use super::*;
8
9/// Rings for which hermite normal forms can be computed
10pub trait HermiteAlgorithmSignature: BezoutDomainSignature {}
11impl<Ring: BezoutDomainSignature> HermiteAlgorithmSignature for Ring {}
12
13/// Rings for which reduced hermite normal forms can be computed
14pub trait ReducedHermiteAlgorithmSignature:
15    HermiteAlgorithmSignature + EuclideanDomainSignature + FavoriteAssociateSignature
16{
17}
18impl<Ring: HermiteAlgorithmSignature + EuclideanDomainSignature + FavoriteAssociateSignature>
19    ReducedHermiteAlgorithmSignature for Ring
20{
21}
22
23/// Rings for which reduced hermite normal forms can be computed and are unique
24pub trait UniqueReducedHermiteAlgorithmSignature: ReducedHermiteAlgorithmSignature {}
25impl UniqueReducedHermiteAlgorithmSignature for IntegerCanonicalStructure {}
26impl<Field: FieldSignature> UniqueReducedHermiteAlgorithmSignature for Field {}
27
28impl<Ring: HermiteAlgorithmSignature, RingB: BorrowedStructure<Ring>> MatrixStructure<Ring, RingB> {
29    /// Return (H, U, u_det, pivots) such that
30    /// - H is in row hermite normal form, meaning
31    /// - U is invertible
32    /// - UM=H
33    /// - u_det is the determinant of u
34    /// - pivots[r] is the column of the rth pivot
35    pub fn row_hermite_algorithm(
36        &self,
37        mut m: Matrix<Ring::Set>,
38    ) -> (Matrix<Ring::Set>, Matrix<Ring::Set>, Ring::Set, Vec<usize>) {
39        //build up U by applying row operations to the identity as we go
40        let mut u = self.ident(m.rows());
41        let mut u_det = self.ring().one();
42        let mut pivs = vec![];
43
44        let (mut pr, mut pc) = (0, 0);
45        'pivot_loop: while pr < m.rows() {
46            //find the next pivot row
47            //the next pivot row is the next row with a non-zero element below the previous pivot
48            'next_pivot_loop: loop {
49                if pc == m.cols() {
50                    break 'pivot_loop;
51                }
52
53                for r in pr..m.rows() {
54                    if !self.ring().equal(m.at(r, pc).unwrap(), &self.ring().zero()) {
55                        break 'next_pivot_loop;
56                    }
57                }
58
59                pc += 1;
60            }
61            pivs.push(pc);
62
63            if pr + 1 < m.rows() {
64                //reduce everything below the pivot to zero
65                for r in pr + 1..m.rows() {
66                    let a = m.at(pr, pc).unwrap();
67                    let b = m.at(r, pc).unwrap();
68                    //if a=0 and b=0 there is nothing to do. The reduction step would fail because d=0 and we divide by d, so just skip it in this case
69                    if !self.ring().equal(a, &self.ring().zero())
70                        || !self.ring().equal(b, &self.ring().zero())
71                    {
72                        let (d, x, y) = self.ring().xgcd(a, b);
73                        debug_assert!(
74                            self.ring().equal(
75                                &self
76                                    .ring()
77                                    .add(&self.ring().mul(&x, a), &self.ring().mul(&y, b)),
78                                &d
79                            )
80                        );
81                        // perform the following row operations on self
82                        // / x  -b/d \
83                        // \ y   a/d /
84                        let row_opp = ElementaryOpp::new_row_opp(
85                            self.ring().clone(),
86                            ElementaryOppType::TwoInv {
87                                i: pr,
88                                j: r,
89                                a: x,
90                                b: y,
91                                //TODO: compute b/d and a/d at the same time d is computed?
92                                c: self.ring().neg(&self.ring().try_divide(b, &d).unwrap()),
93                                d: self.ring().try_divide(a, &d).unwrap(),
94                            },
95                        );
96                        //this will implicitly put the pivot into fav assoc form because that is what gcd does
97                        row_opp.apply(&mut m);
98                        row_opp.apply(&mut u);
99                        self.ring().mul_mut(&mut u_det, &row_opp.det());
100                    }
101                }
102            } else {
103                //explicitly put the pivot into fav assoc form
104                let (unit, _assoc) = self.ring().factor_fav_assoc(m.at(pr, pc).unwrap());
105                let row_opp = ElementaryOpp::new_row_opp(
106                    self.ring().clone(),
107                    ElementaryOppType::UnitMul {
108                        row: pr,
109                        unit: self.ring().try_reciprocal(&unit).unwrap(),
110                    },
111                );
112                //this will implicitly put the pivot into fav assoc form because that is what the gcd returns
113                row_opp.apply(&mut m);
114                row_opp.apply(&mut u);
115                self.ring().mul_mut(&mut u_det, &row_opp.det());
116            }
117
118            //should have eliminated everything below the pivot
119            for r in pr + 1..m.rows() {
120                debug_assert!(self.ring().equal(m.at(r, pc).unwrap(), &self.ring().zero()));
121            }
122            pr += 1;
123        }
124
125        if m.rows() <= 4 {
126            debug_assert!(self.ring().equal(&self.det_naive(&u).unwrap(), &u_det));
127        }
128
129        (m, u, u_det, pivs)
130    }
131
132    pub fn col_hermite_algorithm(
133        &self,
134        a: Matrix<Ring::Set>,
135    ) -> (Matrix<Ring::Set>, Matrix<Ring::Set>, Ring::Set, Vec<usize>) {
136        let (rh, ru, u_det, pivs) = self.row_hermite_algorithm(a.transpose());
137        (rh.transpose(), ru.transpose(), u_det, pivs)
138    }
139
140    fn det_hermite(&self, a: Matrix<Ring::Set>) -> Ring::Set {
141        let n = a.rows();
142        debug_assert_eq!(n, a.cols());
143        let (h, _u, u_det, _pivs) = self.row_hermite_algorithm(a);
144        //h = u * self, we know det(u), and h is upper triangular
145        let mut h_det = self.ring().one();
146        for i in 0..n {
147            self.ring().mul_mut(&mut h_det, h.at(i, i).unwrap());
148        }
149        self.ring().try_divide(&h_det, &u_det).unwrap()
150    }
151
152    pub fn det(&self, a: Matrix<Ring::Set>) -> Result<Ring::Set, MatOppErr> {
153        let n = a.rows();
154        if n != a.cols() {
155            Err(MatOppErr::NotSquare)
156        } else if n <= 3 {
157            //for speed
158            Ok(self.det_naive(&a).unwrap())
159        } else {
160            Ok(self.det_hermite(a))
161        }
162    }
163
164    pub fn rank(&self, a: Matrix<Ring::Set>) -> usize {
165        let (_h, _u, _u_det, pivs) = self.row_hermite_algorithm(a);
166        pivs.len()
167    }
168}
169
170impl<Ring: ReducedHermiteAlgorithmSignature, RingB: BorrowedStructure<Ring>>
171    MatrixStructure<Ring, RingB>
172{
173    /// Returns (H, U, u_det, pivots) such that
174    /// - H is in row reduced hermite normal form, meaning entries above pivots have euclidean norm strictly less than the pivot
175    /// - U is invertible
176    /// - UM=H
177    /// - pivots[r] is the column of the rth pivot and pivots.len() == rank(A)
178    pub fn row_reduced_hermite_algorithm(
179        &self,
180        m: Matrix<Ring::Set>,
181    ) -> (Matrix<Ring::Set>, Matrix<Ring::Set>, Ring::Set, Vec<usize>) {
182        let (mut h, mut u, u_det, pivs) = self.row_hermite_algorithm(m);
183
184        for (pr, pc) in pivs.iter().enumerate() {
185            for r in 0..pr {
186                //reduce h[r, pc] so that it has norm less than h[pr, pc]
187                let a = h.at(r, *pc).unwrap();
188                let b = h.at(pr, *pc).unwrap();
189                //a = b*q + r
190                let q = self.ring().quo(a, b).unwrap();
191                let row_opp = ElementaryOpp::new_row_opp(
192                    self.ring().clone(),
193                    ElementaryOppType::AddRowMul {
194                        i: r,
195                        j: pr,
196                        x: self.ring().neg(&q),
197                    },
198                );
199                row_opp.apply(&mut h);
200                row_opp.apply(&mut u);
201                // adding a multiple of a row does not change the determinant so no need to update u_det
202            }
203        }
204
205        (h, u, u_det, pivs)
206    }
207
208    pub fn row_reduced_hermite_normal_form(&self, m: Matrix<Ring::Set>) -> Matrix<Ring::Set> {
209        self.row_reduced_hermite_algorithm(m).0
210    }
211
212    pub fn col_reduced_hermite_algorithm(
213        &self,
214        m: Matrix<Ring::Set>,
215    ) -> (Matrix<Ring::Set>, Matrix<Ring::Set>, Ring::Set, Vec<usize>) {
216        let (rh, ru, ru_det, pivs) = self.row_reduced_hermite_algorithm(m.transpose());
217        (rh.transpose(), ru.transpose(), ru_det, pivs)
218    }
219
220    pub fn col_reduced_hermite_normal_form(&self, m: Matrix<Ring::Set>) -> Matrix<Ring::Set> {
221        self.col_reduced_hermite_algorithm(m).0
222    }
223
224    pub fn inv(&self, a: Matrix<Ring::Set>) -> Result<Matrix<Ring::Set>, MatOppErr> {
225        let n = a.rows();
226        if n == a.cols() {
227            let (h, u, _u_det, _pivs) = self.row_reduced_hermite_algorithm(a);
228            //h = u*a
229            if self.equal(&h, &self.ident(n)) {
230                Ok(u)
231            } else {
232                Err(MatOppErr::Singular)
233            }
234        } else {
235            Err(MatOppErr::NotSquare)
236        }
237    }
238
239    pub fn row_span(&self, matrix: Matrix<Ring::Set>) -> FinitelyFreeSubmodule<Ring::Set> {
240        FinitelyFreeModuleStructure::<Ring, _>::new(self.ring(), matrix.cols())
241            .into_submodules()
242            .matrix_row_span(matrix)
243    }
244
245    pub fn col_span(&self, matrix: Matrix<Ring::Set>) -> FinitelyFreeSubmodule<Ring::Set> {
246        FinitelyFreeModuleStructure::<Ring, _>::new(self.ring(), matrix.rows())
247            .into_submodules()
248            .matrix_col_span(matrix)
249    }
250
251    pub fn row_kernel(&self, matrix: Matrix<Ring::Set>) -> FinitelyFreeSubmodule<Ring::Set> {
252        FinitelyFreeModuleStructure::<Ring, _>::new(self.ring(), matrix.rows())
253            .into_submodules()
254            .matrix_row_kernel(matrix)
255    }
256
257    pub fn col_kernel(&self, matrix: Matrix<Ring::Set>) -> FinitelyFreeSubmodule<Ring::Set> {
258        FinitelyFreeModuleStructure::<Ring, _>::new(self.ring(), matrix.cols())
259            .into_submodules()
260            .matrix_col_kernel(matrix)
261    }
262
263    pub fn row_preimage(
264        &self,
265        matrix: &Matrix<Ring::Set>,
266        space: &FinitelyFreeSubmodule<Ring::Set>,
267    ) -> FinitelyFreeSubmodule<Ring::Set> {
268        FinitelyFreeModuleStructure::<Ring, _>::new(self.ring(), matrix.rows())
269            .into_submodules()
270            .matrix_row_preimage(matrix, space)
271    }
272
273    pub fn col_preimage(
274        &self,
275        matrix: &Matrix<Ring::Set>,
276        space: &FinitelyFreeSubmodule<Ring::Set>,
277    ) -> FinitelyFreeSubmodule<Ring::Set> {
278        FinitelyFreeModuleStructure::<Ring, _>::new(self.ring(), matrix.rows())
279            .into_submodules()
280            .matrix_col_preimage(matrix, space)
281    }
282
283    pub fn row_affine_span(
284        &self,
285        matrix: Matrix<Ring::Set>,
286    ) -> FinitelyFreeSubmoduleAffineSubset<Ring::Set> {
287        let span = (0..matrix.rows())
288            .map(|r| matrix.get_row(r))
289            .collect::<Vec<_>>();
290        FinitelyFreeModuleStructure::<Ring, _>::new(self.ring(), matrix.cols())
291            .affine_subsets()
292            .from_affine_span(span.iter().collect())
293    }
294
295    pub fn col_affine_span(
296        &self,
297        matrix: Matrix<Ring::Set>,
298    ) -> FinitelyFreeSubmoduleAffineSubset<Ring::Set> {
299        self.row_affine_span(matrix.transpose())
300    }
301
302    pub fn row_solve(
303        &self,
304        matrix: Matrix<Ring::Set>,
305        y: &Vec<Ring::Set>,
306    ) -> Option<Vec<Ring::Set>> {
307        let submodules = FinitelyFreeModuleStructure::<Ring, _>::new(self.ring(), matrix.cols())
308            .into_submodules();
309        let (row_span_submodule, basis_in_terms_of_matrix_rows) =
310            submodules.matrix_row_span_and_basis(matrix);
311        let (offset, y_reduced) = submodules.reduce_element(&row_span_submodule, y);
312        if y_reduced.iter().all(|v| self.ring().is_zero(v)) {
313            Some(
314                self.mul(
315                    &Matrix::from_rows(vec![offset]),
316                    &basis_in_terms_of_matrix_rows,
317                )
318                .unwrap()
319                .get_row(0),
320            )
321        } else {
322            None
323        }
324    }
325
326    pub fn col_solve(
327        &self,
328        matrix: Matrix<Ring::Set>,
329        y: &Vec<Ring::Set>,
330    ) -> Option<Vec<Ring::Set>> {
331        self.row_solve(matrix.transpose(), y)
332    }
333
334    pub fn row_solution_set(
335        &self,
336        matrix: Matrix<Ring::Set>,
337        y: &Vec<Ring::Set>,
338    ) -> FinitelyFreeSubmoduleAffineSubset<Ring::Set> {
339        let module = FinitelyFreeModuleStructure::<Ring, _>::new(self.ring(), matrix.rows());
340        match self.row_solve(matrix.clone(), y) {
341            Some(offset) => FinitelyFreeSubmoduleAffineSubset::NonEmpty(
342                module
343                    .cosets()
344                    .from_offset_and_submodule(&offset, self.row_kernel(matrix)),
345            ),
346            None => FinitelyFreeSubmoduleAffineSubset::Empty,
347        }
348    }
349
350    pub fn col_solution_set(
351        &self,
352        matrix: Matrix<Ring::Set>,
353        y: &Vec<Ring::Set>,
354    ) -> FinitelyFreeSubmoduleAffineSubset<Ring::Set> {
355        self.row_solution_set(matrix.transpose(), y)
356    }
357}
358
359impl<R: MetaType> Matrix<R>
360where
361    R::Signature: BezoutDomainSignature,
362{
363    pub fn row_hermite_algorithm(&self) -> (Self, Self, R, Vec<usize>) {
364        Self::structure().row_hermite_algorithm(self.clone())
365    }
366
367    pub fn col_hermite_algorithm(&self) -> (Self, Self, R, Vec<usize>) {
368        Self::structure().col_hermite_algorithm(self.clone())
369    }
370
371    pub fn det(&self) -> Result<R, MatOppErr> {
372        Self::structure().det(self.clone())
373    }
374
375    pub fn rank(&self) -> usize {
376        Self::structure().rank(self.clone())
377    }
378}
379
380impl<R: MetaType> Matrix<R>
381where
382    R::Signature: EuclideanDomainSignature + BezoutDomainSignature + FavoriteAssociateSignature,
383{
384    pub fn row_reduced_hermite_algorithm(&self) -> (Self, Self, R, Vec<usize>) {
385        Self::structure().row_reduced_hermite_algorithm(self.clone())
386    }
387
388    pub fn row_reduced_hermite_normal_form(&self) -> Self {
389        Self::structure().row_reduced_hermite_normal_form(self.clone())
390    }
391
392    pub fn col_reduced_hermite_algorithm(&self) -> (Self, Self, R, Vec<usize>) {
393        Self::structure().col_reduced_hermite_algorithm(self.clone())
394    }
395
396    pub fn col_reduced_hermite_normal_form(&self) -> Self {
397        Self::structure().col_reduced_hermite_normal_form(self.clone())
398    }
399
400    pub fn inv(&self) -> Result<Matrix<R>, MatOppErr> {
401        Self::structure().inv(self.clone())
402    }
403
404    pub fn row_span(self) -> FinitelyFreeSubmodule<R> {
405        Self::structure().row_span(self)
406    }
407
408    pub fn col_span(self) -> FinitelyFreeSubmodule<R> {
409        Self::structure().col_span(self)
410    }
411
412    pub fn row_kernel(self) -> FinitelyFreeSubmodule<R> {
413        Self::structure().row_kernel(self)
414    }
415
416    pub fn col_kernel(self) -> FinitelyFreeSubmodule<R> {
417        Self::structure().col_kernel(self)
418    }
419
420    pub fn row_preimage(&self, space: &FinitelyFreeSubmodule<R>) -> FinitelyFreeSubmodule<R> {
421        Self::structure().row_preimage(self, space)
422    }
423
424    pub fn col_preimage(&self, space: &FinitelyFreeSubmodule<R>) -> FinitelyFreeSubmodule<R> {
425        Self::structure().col_preimage(self, space)
426    }
427
428    pub fn row_affine_span(self) -> FinitelyFreeSubmoduleAffineSubset<R> {
429        Self::structure().row_affine_span(self)
430    }
431
432    pub fn col_affine_span(self) -> FinitelyFreeSubmoduleAffineSubset<R> {
433        Self::structure().col_affine_span(self)
434    }
435
436    pub fn row_solve(self, y: &Vec<R>) -> Option<Vec<R>> {
437        Self::structure().row_solve(self, y)
438    }
439
440    pub fn col_solve(self, y: &Vec<R>) -> Option<Vec<R>> {
441        Self::structure().col_solve(self, y)
442    }
443
444    pub fn row_solution_set(self, y: &Vec<R>) -> FinitelyFreeSubmoduleAffineSubset<R> {
445        Self::structure().row_solution_set(self, y)
446    }
447
448    pub fn col_solution_set(self, y: &Vec<R>) -> FinitelyFreeSubmoduleAffineSubset<R> {
449        Self::structure().col_solution_set(self, y)
450    }
451}
452
453#[cfg(test)]
454mod tests {
455    use crate::module::finitely_free_module::RingToFinitelyFreeModuleSignature;
456
457    use super::*;
458
459    #[test]
460    fn hermite_algorithm() {
461        for a in [
462            Matrix::from_rows(vec![
463                vec![
464                    Integer::from(2),
465                    Integer::from(3),
466                    Integer::from(6),
467                    Integer::from(2),
468                ],
469                vec![
470                    Integer::from(5),
471                    Integer::from(6),
472                    Integer::from(1),
473                    Integer::from(6),
474                ],
475                vec![
476                    Integer::from(8),
477                    Integer::from(3),
478                    Integer::from(1),
479                    Integer::from(1),
480                ],
481            ]),
482            Matrix::from_rows(vec![
483                vec![
484                    Integer::from(2),
485                    Integer::from(3),
486                    Integer::from(6),
487                    Integer::from(2),
488                ],
489                vec![
490                    Integer::from(5),
491                    Integer::from(6),
492                    Integer::from(1),
493                    Integer::from(6),
494                ],
495                vec![
496                    Integer::from(8),
497                    Integer::from(3),
498                    Integer::from(1),
499                    Integer::from(1),
500                ],
501            ])
502            .transpose(),
503            Matrix::from_rows(vec![
504                vec![
505                    Integer::from(2),
506                    Integer::from(3),
507                    Integer::from(5),
508                    Integer::from(2),
509                ],
510                vec![
511                    Integer::from(5),
512                    Integer::from(6),
513                    Integer::from(11),
514                    Integer::from(6),
515                ],
516                vec![
517                    Integer::from(8),
518                    Integer::from(3),
519                    Integer::from(11),
520                    Integer::from(1),
521                ],
522            ]),
523            Matrix::from_rows(vec![
524                vec![Integer::from(0), Integer::from(4), Integer::from(4)],
525                vec![Integer::from(1), Integer::from(6), Integer::from(12)],
526                vec![Integer::from(1), Integer::from(4), Integer::from(16)],
527            ]),
528        ] {
529            println!();
530            println!("hermite reduced row algorithm for");
531            a.pprint();
532            let (h, u, _u_det, pivs) = a.clone().row_reduced_hermite_algorithm();
533            println!("H =");
534            h.pprint();
535            println!("pivs = {:?}", pivs);
536            println!("U =");
537            u.pprint();
538            assert_eq!(h, Matrix::mul(&u, &a).unwrap());
539
540            //trace the boundary of zeros and check that everything under is zero
541            let mut rz = 0;
542            for cz in 0..h.cols() {
543                if let Some(cp) = pivs.get(rz)
544                    && cp == &cz
545                {
546                    rz += 1;
547                }
548                for r in rz..h.rows() {
549                    assert_eq!(h.at(r, cz).unwrap(), &Integer::zero());
550                }
551            }
552
553            //check pivot columns
554            for (pr, pc) in pivs.iter().enumerate() {
555                assert!(h.at(pr, *pc).unwrap() != &Integer::zero());
556                for r in 0..h.rows() {
557                    #[allow(clippy::comparison_chain)]
558                    if r > pr {
559                        assert_eq!(h.at(r, *pc).unwrap(), &Integer::zero());
560                    } else if r == pr {
561                        let (_unit, assoc) = Integer::factor_fav_assoc(h.at(r, *pc).unwrap());
562                        assert_eq!(&assoc, h.at(r, *pc).unwrap());
563                    } else {
564                        assert!(
565                            Integer::norm(h.at(r, *pc).unwrap())
566                                < Integer::norm(h.at(pr, *pc).unwrap())
567                        );
568                    }
569                }
570            }
571
572            println!();
573            println!("hermite reduced col algorithm for");
574            a.pprint();
575            let (h, u, _u_det, pivs) = a.clone().col_reduced_hermite_algorithm();
576            println!("H =");
577            h.pprint();
578            println!("pivs = {:?}", pivs);
579            println!("U =");
580            u.pprint();
581
582            //trace the boundary of zeros and check that everything to the right is zero
583            let mut cz = 0;
584            for rz in 0..h.rows() {
585                if let Some(rp) = pivs.get(cz)
586                    && rp == &rz
587                {
588                    cz += 1;
589                }
590                for c in cz..h.cols() {
591                    assert_eq!(h.at(rz, c).unwrap(), &Integer::zero());
592                }
593            }
594
595            //check the pivot rows
596            assert_eq!(h, Matrix::mul(&a, &u).unwrap());
597            for (pc, pr) in pivs.iter().enumerate() {
598                assert!(h.at(*pr, pc).unwrap() != &Integer::zero());
599                for c in 0..h.cols() {
600                    #[allow(clippy::comparison_chain)]
601                    if c > pc {
602                        assert_eq!(h.at(*pr, c).unwrap(), &Integer::zero());
603                    } else if c == pc {
604                        let (_unit, assoc) = Integer::factor_fav_assoc(h.at(*pr, c).unwrap());
605                        assert_eq!(&assoc, h.at(*pr, c).unwrap());
606                    } else {
607                        assert!(
608                            Integer::norm(h.at(*pr, c).unwrap())
609                                < Integer::norm(h.at(*pr, pc).unwrap())
610                        );
611                    }
612                }
613            }
614        }
615
616        {
617            //integer reduced hermite normal form is unique, so we can fully check an example
618            let a = Matrix::<Integer>::from_rows(vec![
619                vec![
620                    Integer::from(2),
621                    Integer::from(3),
622                    Integer::from(6),
623                    Integer::from(2),
624                ],
625                vec![
626                    Integer::from(5),
627                    Integer::from(6),
628                    Integer::from(1),
629                    Integer::from(6),
630                ],
631                vec![
632                    Integer::from(8),
633                    Integer::from(3),
634                    Integer::from(1),
635                    Integer::from(1),
636                ],
637            ]);
638
639            let expected_h = Matrix::from_rows(vec![
640                vec![
641                    Integer::from(1),
642                    Integer::from(0),
643                    Integer::from(50),
644                    Integer::from(-11),
645                ],
646                vec![
647                    Integer::from(0),
648                    Integer::from(3),
649                    Integer::from(28),
650                    Integer::from(-2),
651                ],
652                vec![
653                    Integer::from(0),
654                    Integer::from(0),
655                    Integer::from(61),
656                    Integer::from(-13),
657                ],
658            ]);
659
660            let expected_u = Matrix::from_rows(vec![
661                vec![Integer::from(9), Integer::from(-5), Integer::from(1)],
662                vec![Integer::from(5), Integer::from(-2), Integer::from(0)],
663                vec![Integer::from(11), Integer::from(-6), Integer::from(1)],
664            ]);
665
666            let (h, u, _u_det, pivs) = a.clone().row_reduced_hermite_algorithm();
667
668            assert_eq!(h, expected_h);
669            assert_eq!(u, expected_u);
670            assert_eq!(pivs, vec![0, 1, 2]);
671        }
672
673        {
674            //this one used to cause a dividion by zero error when replacing (a, b) with (gcd, 0) when (a, b) = (0, 0)
675            let a = Matrix::<Integer>::from_rows(vec![
676                vec![Integer::from(1), Integer::from(0), Integer::from(0)],
677                vec![Integer::from(0), Integer::from(1), Integer::from(0)],
678                vec![Integer::from(0), Integer::from(0), Integer::from(0)],
679                vec![Integer::from(0), Integer::from(0), Integer::from(0)],
680                vec![Integer::from(0), Integer::from(0), Integer::from(1)],
681            ]);
682
683            let expected_h = Matrix::from_rows(vec![
684                vec![Integer::from(1), Integer::from(0), Integer::from(0)],
685                vec![Integer::from(0), Integer::from(1), Integer::from(0)],
686                vec![Integer::from(0), Integer::from(0), Integer::from(1)],
687                vec![Integer::from(0), Integer::from(0), Integer::from(0)],
688                vec![Integer::from(0), Integer::from(0), Integer::from(0)],
689            ]);
690
691            let expected_u = Matrix::from_rows(vec![
692                vec![
693                    Integer::from(1),
694                    Integer::from(0),
695                    Integer::from(0),
696                    Integer::from(0),
697                    Integer::from(0),
698                ],
699                vec![
700                    Integer::from(0),
701                    Integer::from(1),
702                    Integer::from(0),
703                    Integer::from(0),
704                    Integer::from(0),
705                ],
706                vec![
707                    Integer::from(0),
708                    Integer::from(0),
709                    Integer::from(0),
710                    Integer::from(0),
711                    Integer::from(1),
712                ],
713                vec![
714                    Integer::from(0),
715                    Integer::from(0),
716                    Integer::from(0),
717                    Integer::from(1),
718                    Integer::from(0),
719                ],
720                vec![
721                    Integer::from(0),
722                    Integer::from(0),
723                    Integer::from(-1),
724                    Integer::from(0),
725                    Integer::from(0),
726                ],
727            ]);
728
729            let (h, u, _u_det, pivs) = a.clone().row_reduced_hermite_algorithm();
730
731            assert_eq!(h, expected_h);
732            assert_eq!(u, expected_u);
733            assert_eq!(pivs, vec![0, 1, 2]);
734        }
735    }
736
737    #[test]
738    fn invert() {
739        let a = Matrix::<Rational>::from_rows(vec![
740            vec![Rational::from(2), Rational::from(4), Rational::from(4)],
741            vec![Rational::from(-6), Rational::from(6), Rational::from(12)],
742            vec![Rational::from(10), Rational::from(7), Rational::from(17)],
743        ]);
744        a.pprint();
745        println!("{:?}", a.rank());
746        let b = Matrix::inv(&a).unwrap();
747        b.pprint();
748        debug_assert_eq!(Matrix::mul(&a, &b).unwrap(), Matrix::ident(3));
749        debug_assert_eq!(Matrix::mul(&b, &a).unwrap(), Matrix::ident(3));
750    }
751
752    #[test]
753    fn span_and_kernel_rank() {
754        let mat = Matrix::<Integer>::from_rows(vec![
755            vec![Integer::from(1), Integer::from(1), Integer::from(1)],
756            vec![Integer::from(1), Integer::from(2), Integer::from(1)],
757            vec![Integer::from(1), Integer::from(1), Integer::from(1)],
758            vec![Integer::from(1), Integer::from(1), Integer::from(1)],
759        ]);
760
761        assert_eq!(mat.clone().row_span().rank(), 2);
762        assert_eq!(mat.clone().col_span().rank(), 2);
763        assert_eq!(mat.clone().row_kernel().rank(), 2);
764        assert_eq!(mat.clone().col_kernel().rank(), 1);
765    }
766
767    #[test]
768    fn affine_span() {
769        {
770            let module = Integer::structure().into_free_module(2);
771
772            //row affine span
773            let lat1 = Matrix::<Integer>::from_rows(vec![
774                vec![Integer::from(1), Integer::from(1)],
775                vec![Integer::from(3), Integer::from(1)],
776                vec![Integer::from(2), Integer::from(3)],
777            ])
778            .row_affine_span();
779
780            let lat2 = FinitelyFreeSubmoduleAffineSubset::NonEmpty(
781                module.cosets().from_offset_and_submodule(
782                    &vec![Integer::from(2), Integer::from(3)],
783                    Matrix::<Integer>::from_rows(vec![vec![1, 2], vec![-1, 2]]).row_span(),
784                ),
785            );
786
787            println!("lat1 = {:?}", lat1);
788            println!("lat2 = {:?}", lat2);
789
790            assert!(module.affine_subsets().equal(&lat1, &lat2));
791            assert!(module.affine_subsets().equal_slow(&lat1, &lat2));
792        }
793
794        {
795            let module = Integer::structure().into_free_module(2);
796
797            //column affine span
798            let lat1 = Matrix::<Integer>::from_rows(vec![
799                vec![Integer::from(1), Integer::from(3), Integer::from(2)],
800                vec![Integer::from(1), Integer::from(1), Integer::from(3)],
801            ])
802            .col_affine_span();
803
804            let lat2 = FinitelyFreeSubmoduleAffineSubset::NonEmpty(
805                module.cosets().from_offset_and_submodule(
806                    &vec![Integer::from(2), Integer::from(3)],
807                    Matrix::<Integer>::from_rows(vec![vec![1, -1], vec![2, 2]]).col_span(),
808                ),
809            );
810
811            println!("lat1 = {:?}", lat1);
812            println!("lat2 = {:?}", lat2);
813
814            assert!(module.affine_subsets().equal(&lat1, &lat2));
815            assert!(module.affine_subsets().equal_slow(&lat1, &lat2));
816        }
817    }
818
819    #[test]
820    fn span_and_kernel_points() {
821        let module = Integer::structure().into_free_module(4);
822
823        let mat = Matrix::<Integer>::from_rows(vec![
824            vec![
825                Integer::from(1),
826                Integer::from(1),
827                Integer::from(0),
828                Integer::from(0),
829            ],
830            vec![
831                Integer::from(1),
832                Integer::from(1),
833                Integer::from(0),
834                Integer::from(0),
835            ],
836            vec![
837                Integer::from(0),
838                Integer::from(0),
839                Integer::from(3),
840                Integer::from(5),
841            ],
842        ]);
843
844        println!("matrix");
845        mat.pprint();
846
847        let k = mat.col_kernel();
848
849        assert!(module.submodules().contains_element(
850            &k,
851            &vec![
852                Integer::from(-1),
853                Integer::from(1),
854                Integer::from(5),
855                Integer::from(-3)
856            ]
857        ));
858    }
859
860    #[test]
861    fn test_row_solve() {
862        let module = Integer::structure().into_free_module(3);
863
864        let matrix =
865            Matrix::<Integer>::from_rows(vec![vec![1, 0, 0], vec![1, 0, 1], vec![1, 1, 1]]);
866        let x =
867            matrix
868                .clone()
869                .row_solve(&vec![Integer::from(4), Integer::from(4), Integer::from(7)]);
870        assert_eq!(
871            x.unwrap(),
872            vec![Integer::from(-3), Integer::from(3), Integer::from(4),]
873        );
874
875        let a = Matrix::<Integer>::from_rows(vec![vec![1, 0, 0], vec![1, 0, 1]])
876            .row_solution_set(&vec![Integer::from(2), Integer::from(3), Integer::from(2)]);
877        for s in module.affine_subsets().affine_basis(&a) {
878            println!("{:?}", s);
879        }
880    }
881}
algebraeon_rings/matrix/hermite_reduction.rs

algebraeon_rings/matrix/
hermite_reduction.rs
