feanor_math/algorithms/linsolve/
smith.rs

1use std::alloc::Allocator;
2use std::cmp::min;
3
4use transform::TransformList;
5
6use crate::algorithms::linsolve::SolveResult;
7use crate::divisibility::*;
8use crate::matrix::transform::{TransformCols, TransformRows, TransformTarget};
9use crate::matrix::*;
10use crate::pid::{PrincipalIdealRing, PrincipalIdealRingStore};
11use crate::ring::*;
12
13/// Transforms `A` into `A'` via transformations `L, R` such that
14/// `L A R = A'` and `A'` is diagonal.
15///
16/// # (Non-)Uniqueness of the solution
17///
18/// Note that this is not the complete Smith normal form,
19/// as that requires the entries on the diagonal to divide
20/// each other. However, computing this pre-smith form is much
21/// faster, and can still be used for solving equations (the main
22/// use case). However, it is not unique.
23///
24/// # Warning on infinite rings (in particular Z)
25///
26/// For infinite principal ideal rings, this function is correct,
27/// but in some situations, the performance can be terrible. The
28/// reason is that no care is taken which of the many possible results
29/// is returned - and in fact, this algorithm can sometimes choose
30/// one that has exponential size in the input. Hence, in these
31/// cases it is recommended to use another algorithm, e.g. based on
32/// LLL to perform intermediate lattice reductions (not yet implemented
33/// in feanor_math).
34#[stability::unstable(feature = "enable")]
35pub fn pre_smith<R, TL, TR, V>(ring: R, L: &mut TL, R: &mut TR, mut A: SubmatrixMut<V, El<R>>)
36where
37    R: RingStore + Copy,
38    R::Type: PrincipalIdealRing,
39    TL: TransformTarget<R::Type>,
40    TR: TransformTarget<R::Type>,
41    V: AsPointerToSlice<El<R>>,
42{
43    // otherwise we might not terminate...
44    assert!(ring.is_noetherian());
45    assert!(ring.is_commutative());
46
47    for k in 0..min(A.row_count(), A.col_count()) {
48        let mut changed_row = true;
49        while changed_row {
50            changed_row = false;
51
52            // eliminate the column
53            for i in (k + 1)..A.row_count() {
54                if ring.is_zero(A.at(i, k)) {
55                    continue;
56                } else if let Some(quo) = ring.checked_div(A.at(i, k), A.at(k, k)) {
57                    TransformRows(A.reborrow(), ring.get_ring()).subtract(ring, k, i, &quo);
58                    L.subtract(ring, k, i, &quo);
59                } else {
60                    let (transform, _) = ring.get_ring().create_elimination_matrix(A.at(k, k), A.at(i, k));
61                    TransformRows(A.reborrow(), ring.get_ring()).transform(ring, k, i, &transform);
62                    L.transform(ring, k, i, &transform);
63                }
64            }
65
66            // now eliminate the row
67            for j in (k + 1)..A.col_count() {
68                if ring.is_zero(A.at(k, j)) {
69                    continue;
70                } else if let Some(quo) = ring.checked_div(A.at(k, j), A.at(k, k)) {
71                    changed_row = true;
72                    TransformCols(A.reborrow(), ring.get_ring()).subtract(ring, k, j, &quo);
73                    R.subtract(ring, k, j, &quo);
74                } else {
75                    changed_row = true;
76                    let (transform, _) = ring.get_ring().create_elimination_matrix(A.at(k, k), A.at(k, j));
77                    TransformCols(A.reborrow(), ring.get_ring()).transform(ring, k, j, &transform);
78                    R.transform(ring, k, j, &transform);
79                }
80            }
81        }
82    }
83}
84
85/// Computes a matrix `X` such that `lhs * X = rhs`, if it exists.
86///
87/// This function will change the value of `A`, more concretely overwrite it
88/// with its "pre-smith" form as specified by [`pre_smith()`].
89#[stability::unstable(feature = "enable")]
90pub fn solve_right_using_pre_smith<R, V1, V2, V3, A>(
91    ring: R,
92    mut lhs: SubmatrixMut<V1, El<R>>,
93    mut rhs: SubmatrixMut<V2, El<R>>,
94    mut out: SubmatrixMut<V3, El<R>>,
95    _allocator: A,
96) -> SolveResult
97where
98    R: RingStore + Copy,
99    R::Type: PrincipalIdealRing,
100    V1: AsPointerToSlice<El<R>>,
101    V2: AsPointerToSlice<El<R>>,
102    V3: AsPointerToSlice<El<R>>,
103    A: Allocator,
104{
105    assert_eq!(lhs.row_count(), rhs.row_count());
106    assert_eq!(lhs.col_count(), out.row_count());
107    assert_eq!(rhs.col_count(), out.col_count());
108
109    let mut R = TransformList::new(lhs.col_count());
110    pre_smith(
111        ring,
112        &mut TransformRows(rhs.reborrow(), ring.get_ring()),
113        &mut R,
114        lhs.reborrow(),
115    );
116
117    for i in out.row_count()..rhs.row_count() {
118        for j in 0..rhs.col_count() {
119            if !ring.is_zero(rhs.at(i, j)) {
120                return SolveResult::NoSolution;
121            }
122        }
123    }
124    // the value of out[lhs.row_count().., ..] is irrelevant, since lhs
125    // is zero in these places anyway. Thus we just leave it unchanged
126
127    let mut solution_unique = true;
128
129    for i in 0..min(lhs.row_count(), lhs.col_count()) {
130        let pivot = lhs.at(i, i);
131        for j in 0..rhs.col_count() {
132            if let Some(quo) = ring.checked_left_div(rhs.at(i, j), pivot) {
133                solution_unique &= ring.is_zero(&ring.annihilator(pivot));
134                *out.at_mut(i, j) = quo;
135            } else {
136                return SolveResult::NoSolution;
137            }
138        }
139    }
140
141    R.replay_transposed(ring, TransformRows(out, ring.get_ring()));
142    return if solution_unique {
143        SolveResult::FoundUniqueSolution
144    } else {
145        SolveResult::FoundSomeSolution
146    };
147}
148
149/// Computes a basis of the right-kernel of the matrix `A`, i.e. a full-rank
150/// matrix `B` such that `AB = 0` and every `x` with `Ax = 0` is of the form
151/// `x = By` for some `y`.
152///
153/// This function will change the value of `A`, more concretely overwrite it
154/// with its "pre-smith" form as specified by [`pre_smith()`].
155///
156/// **Note on rings with zero-divisors** If the given ring has zero-divisors,
157/// the notion of being "full-rank" is not well-defined. Instead, `B` will have
158/// the property that every minor is nonzero. Note that it is not required that
159/// every minor is invertible - this is usually not possible. For example, in
160/// `Z/6Z`, the right-kernel of `A = (2)` is `{0, 3}`, and thus the only generating
161/// set is `B = (3)`, which doesn't have a unit minor.
162#[stability::unstable(feature = "enable")]
163pub fn kernel_basis_using_pre_smith<R, V, A>(
164    ring: R,
165    mut A: SubmatrixMut<V, El<R>>,
166    allocator: A,
167) -> OwnedMatrix<El<R>, A>
168where
169    R: RingStore + Copy,
170    R::Type: PrincipalIdealRing,
171    V: AsPointerToSlice<El<R>>,
172    A: Allocator,
173{
174    let mut R = TransformList::new(A.col_count());
175    pre_smith(ring, &mut (), &mut R, A.reborrow());
176
177    let annihilators = (0..A.col_count())
178        .map(|i| {
179            if i < A.row_count() {
180                ring.annihilator(A.at(i, i))
181            } else {
182                ring.one()
183            }
184        })
185        .collect::<Vec<_>>();
186    let annihilators = annihilators
187        .iter()
188        .enumerate()
189        .filter(|(_, a)| !ring.is_zero(a))
190        .enumerate();
191    let mut B = OwnedMatrix::zero_in(A.col_count(), annihilators.clone().count(), ring, allocator);
192    for (i, (j, a)) in annihilators {
193        *B.at_mut(i, j) = ring.clone_el(a);
194    }
195    R.replay_transposed(ring, &mut TransformRows(B.data_mut(), ring.get_ring()));
196    return B;
197}
198
199#[stability::unstable(feature = "enable")]
200pub fn determinant_using_pre_smith<R, V, A>(ring: R, mut matrix: SubmatrixMut<V, El<R>>, _allocator: A) -> El<R>
201where
202    R: RingStore + Copy,
203    R::Type: PrincipalIdealRing,
204    V: AsPointerToSlice<El<R>>,
205    A: Allocator,
206{
207    assert_eq!(matrix.row_count(), matrix.col_count());
208    let mut unit_part_rows = ring.one();
209    let mut unit_part_cols = ring.one();
210    pre_smith(
211        ring,
212        &mut DetUnit {
213            current_unit: &mut unit_part_rows,
214        },
215        &mut DetUnit {
216            current_unit: &mut unit_part_cols,
217        },
218        matrix.reborrow(),
219    );
220    return ring
221        .checked_div(
222            &ring.prod((0..matrix.row_count()).map(|i| ring.clone_el(matrix.at(i, i)))),
223            &ring.prod([unit_part_rows, unit_part_cols]),
224        )
225        .unwrap();
226}
227
228struct DetUnit<'a, R: ?Sized + RingBase> {
229    current_unit: &'a mut R::Element,
230}
231
232impl<'a, R> TransformTarget<R> for DetUnit<'a, R>
233where
234    R: ?Sized + RingBase,
235{
236    fn subtract<S: Copy + RingStore<Type = R>>(
237        &mut self,
238        _ring: S,
239        _src: usize,
240        _dst: usize,
241        _factor: &<R as RingBase>::Element,
242    ) {
243        // determinant does not change
244    }
245
246    fn swap<S: Copy + RingStore<Type = R>>(&mut self, ring: S, _i: usize, _j: usize) {
247        ring.negate_inplace(self.current_unit)
248    }
249
250    fn transform<S: Copy + RingStore<Type = R>>(
251        &mut self,
252        ring: S,
253        _i: usize,
254        _j: usize,
255        transform: &[<R as RingBase>::Element; 4],
256    ) {
257        let unit = ring.sub(
258            ring.mul_ref(&transform[0], &transform[3]),
259            ring.mul_ref(&transform[1], &transform[2]),
260        );
261        ring.mul_assign(self.current_unit, unit);
262    }
263}
264
265#[cfg(test)]
266use std::alloc::Global;
267#[cfg(test)]
268use std::ptr::Alignment;
269#[cfg(test)]
270use std::rc::Rc;
271#[cfg(test)]
272use std::time::Instant;
273
274#[cfg(test)]
275use test::Bencher;
276
277#[cfg(test)]
278use crate::algorithms::convolution::STANDARD_CONVOLUTION;
279#[cfg(test)]
280use crate::algorithms::linsolve::LinSolveRing;
281#[cfg(test)]
282use crate::algorithms::linsolve::extension::solve_right_over_extension;
283#[cfg(test)]
284use crate::algorithms::matmul::ComputeInnerProduct;
285#[cfg(test)]
286use crate::algorithms::matmul::*;
287#[cfg(test)]
288use crate::assert_matrix_eq;
289#[cfg(test)]
290use crate::homomorphism::Homomorphism;
291#[cfg(test)]
292use crate::primitive_int::StaticRing;
293#[cfg(test)]
294use crate::rings::extension::FreeAlgebraStore;
295#[cfg(test)]
296use crate::rings::extension::extension_impl::FreeAlgebraImpl;
297#[cfg(test)]
298use crate::rings::extension::galois_field::GaloisField;
299#[cfg(test)]
300use crate::rings::zn::ZnRingStore;
301#[cfg(test)]
302use crate::rings::zn::zn_64::Zn;
303#[cfg(test)]
304use crate::rings::zn::zn_static;
305#[cfg(test)]
306use crate::seq::VectorView;
307
308#[cfg(test)]
309fn multiply<'a, R: RingStore, V: AsPointerToSlice<El<R>>, I: IntoIterator<Item = Submatrix<'a, V, El<R>>>>(
310    matrices: I,
311    ring: R,
312) -> OwnedMatrix<El<R>>
313where
314    R::Type: 'a,
315    V: 'a,
316{
317    let mut it = matrices.into_iter();
318    let fst = it.next().unwrap();
319    let snd = it.next().unwrap();
320    let mut new_result = OwnedMatrix::zero(fst.row_count(), snd.col_count(), &ring);
321    STANDARD_MATMUL.matmul(
322        TransposableSubmatrix::from(fst),
323        TransposableSubmatrix::from(snd),
324        TransposableSubmatrixMut::from(new_result.data_mut()),
325        &ring,
326    );
327    let mut result = new_result;
328
329    for m in it {
330        let mut new_result = OwnedMatrix::zero(result.row_count(), m.col_count(), &ring);
331        STANDARD_MATMUL.matmul(
332            TransposableSubmatrix::from(result.data()),
333            TransposableSubmatrix::from(m),
334            TransposableSubmatrixMut::from(new_result.data_mut()),
335            &ring,
336        );
337        result = new_result;
338    }
339    return result;
340}
341
342#[test]
343fn test_smith_integers() {
344    let ring = StaticRing::<i64>::RING;
345    let mut A = OwnedMatrix::new(vec![1, 2, 3, 4, 2, 3, 4, 5, 3, 4, 5, 6], 4);
346    let original_A = A.clone_matrix(&ring);
347    let mut L: OwnedMatrix<i64> = OwnedMatrix::identity(3, 3, StaticRing::<i64>::RING);
348    let mut R: OwnedMatrix<i64> = OwnedMatrix::identity(4, 4, StaticRing::<i64>::RING);
349    pre_smith(
350        ring,
351        &mut TransformRows(L.data_mut(), ring.get_ring()),
352        &mut TransformCols(R.data_mut(), ring.get_ring()),
353        A.data_mut(),
354    );
355
356    assert_matrix_eq!(&ring, &[[1, 0, 0, 0], [0, -1, 0, 0], [0, 0, 0, 0]], &A);
357
358    assert_matrix_eq!(&ring, &multiply([L.data(), original_A.data(), R.data()], ring), &A);
359}
360
361#[test]
362fn test_smith_zn() {
363    let ring = zn_static::Zn::<45>::RING;
364    let mut A = OwnedMatrix::new(vec![8, 3, 5, 8, 0, 9, 0, 9, 5, 9, 5, 14, 8, 3, 5, 23, 3, 39, 0, 39], 4);
365    let original_A = A.clone_matrix(&ring);
366    let mut L: OwnedMatrix<u64> = OwnedMatrix::identity(5, 5, ring);
367    let mut R: OwnedMatrix<u64> = OwnedMatrix::identity(4, 4, ring);
368    pre_smith(
369        ring,
370        &mut TransformRows(L.data_mut(), ring.get_ring()),
371        &mut TransformCols(R.data_mut(), ring.get_ring()),
372        A.data_mut(),
373    );
374
375    assert_matrix_eq!(
376        &ring,
377        &[[8, 0, 0, 0], [0, 3, 0, 0], [0, 0, 0, 0], [0, 0, 0, 15], [0, 0, 0, 0]],
378        &A
379    );
380
381    assert_matrix_eq!(&ring, &multiply([L.data(), original_A.data(), R.data()], ring), &A);
382}
383
384#[test]
385fn test_solve_zn() {
386    let ring = zn_static::Zn::<45>::RING;
387    let A = OwnedMatrix::new(vec![8, 3, 5, 8, 0, 9, 0, 9, 5, 9, 5, 14, 8, 3, 5, 23, 3, 39, 0, 39], 4);
388    let B = OwnedMatrix::new(
389        vec![11, 43, 10, 22, 18, 9, 27, 27, 8, 34, 7, 22, 41, 13, 40, 37, 3, 9, 3, 0],
390        4,
391    );
392    let mut solution: OwnedMatrix<_> = OwnedMatrix::zero(4, 4, ring);
393    ring.get_ring()
394        .solve_right(
395            A.clone_matrix(ring).data_mut(),
396            B.clone_matrix(ring).data_mut(),
397            solution.data_mut(),
398            Global,
399        )
400        .assert_solved();
401
402    assert_matrix_eq!(&ring, &multiply([A.data(), solution.data()], ring), &B);
403}
404
405#[test]
406fn test_unique_solution_correct() {
407    let ring = zn_static::Zn::<45>::RING;
408    let A = OwnedMatrix::new(vec![1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1], 4);
409    let B = OwnedMatrix::new(vec![1, 1, 0, 0], 1);
410    let mut solution: OwnedMatrix<_> = OwnedMatrix::zero(4, 1, ring);
411    assert_eq!(
412        SolveResult::FoundUniqueSolution,
413        ring.get_ring().solve_right(
414            A.clone_matrix(ring).data_mut(),
415            B.clone_matrix(ring).data_mut(),
416            solution.data_mut(),
417            Global
418        )
419    );
420    assert_matrix_eq!(&ring, &multiply([A.data(), solution.data()], ring), &B);
421
422    let A = OwnedMatrix::new(vec![1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], 4);
423    let B = OwnedMatrix::new(vec![1, 1, 0, 0], 1);
424    let mut solution: OwnedMatrix<_> = OwnedMatrix::zero(4, 1, ring);
425    assert_eq!(
426        SolveResult::FoundSomeSolution,
427        ring.get_ring().solve_right(
428            A.clone_matrix(ring).data_mut(),
429            B.clone_matrix(ring).data_mut(),
430            solution.data_mut(),
431            Global
432        )
433    );
434    assert_matrix_eq!(&ring, &multiply([A.data(), solution.data()], ring), &B);
435
436    let A = OwnedMatrix::new(vec![1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 3], 4);
437    let B = OwnedMatrix::new(vec![1, 1, 0, 0], 1);
438    let mut solution: OwnedMatrix<_> = OwnedMatrix::zero(4, 1, ring);
439    assert_eq!(
440        SolveResult::FoundSomeSolution,
441        ring.get_ring().solve_right(
442            A.clone_matrix(ring).data_mut(),
443            B.clone_matrix(ring).data_mut(),
444            solution.data_mut(),
445            Global
446        )
447    );
448    assert_matrix_eq!(&ring, &multiply([A.data(), solution.data()], ring), &B);
449}
450
451#[test]
452fn test_solve_int() {
453    let ring = StaticRing::<i64>::RING;
454    let A = OwnedMatrix::new(vec![3, 6, 2, 0, 4, 7, 5, 5, 4, 5, 5, 5], 6);
455    let B: OwnedMatrix<i64> = OwnedMatrix::identity(2, 2, ring);
456    let mut solution: OwnedMatrix<i64> = OwnedMatrix::zero(6, 2, ring);
457    ring.get_ring()
458        .solve_right(
459            A.clone_matrix(ring).data_mut(),
460            B.clone_matrix(ring).data_mut(),
461            solution.data_mut(),
462            Global,
463        )
464        .assert_solved();
465
466    assert_matrix_eq!(&ring, &multiply([A.data(), solution.data()], &ring), &B);
467}
468
469#[test]
470fn test_large() {
471    let ring = zn_static::Zn::<16>::RING;
472    let data_A = [
473        [0, 0, 0, 0, 0, 0, 0, 0, 11, 0, 0],
474        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
475        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 10],
476        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
477        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 8],
478        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
479    ];
480    let mut A: OwnedMatrix<u64> = OwnedMatrix::zero(6, 11, &ring);
481    for i in 0..6 {
482        for j in 0..11 {
483            *A.at_mut(i, j) = data_A[i][j];
484        }
485    }
486    let mut solution: OwnedMatrix<_> = OwnedMatrix::zero(11, 11, ring);
487    assert!(
488        ring.get_ring()
489            .solve_right(
490                A.clone_matrix(&ring).data_mut(),
491                A.data_mut(),
492                solution.data_mut(),
493                Global
494            )
495            .is_solved()
496    );
497}
498
499#[test]
500fn test_determinant() {
501    let ring = StaticRing::<i64>::RING;
502    let A = OwnedMatrix::new(vec![1, 0, 3, 2, 1, 0, 9, 8, 7], 3);
503    assert_el_eq!(
504        ring,
505        (7 + 48 - 27),
506        determinant_using_pre_smith(ring, A.clone_matrix(&ring).data_mut(), Global)
507    );
508
509    // we need a ring that has units of order > 2 to test whether an inversion is necessary for
510    // the accumulated determinant units
511    #[derive(PartialEq, Clone, Copy)]
512    struct TestRing;
513    use crate::delegate::DelegateRing;
514    impl DelegateRing for TestRing {
515        type Base = zn_static::ZnBase<45, false>;
516        type Element = u64;
517
518        fn get_delegate(&self) -> &Self::Base { zn_static::Zn::RING.get_ring() }
519        fn delegate(&self, el: Self::Element) -> <Self::Base as RingBase>::Element { el }
520        fn delegate_mut<'a>(&self, el: &'a mut Self::Element) -> &'a mut <Self::Base as RingBase>::Element { el }
521        fn delegate_ref<'a>(&self, el: &'a Self::Element) -> &'a <Self::Base as RingBase>::Element { el }
522        fn rev_delegate(&self, el: <Self::Base as RingBase>::Element) -> Self::Element { el }
523    }
524    impl PrincipalIdealRing for TestRing {
525        fn extended_ideal_gen(
526            &self,
527            lhs: &Self::Element,
528            rhs: &Self::Element,
529        ) -> (Self::Element, Self::Element, Self::Element) {
530            self.get_delegate().extended_ideal_gen(lhs, rhs)
531        }
532
533        fn checked_div_min(&self, lhs: &Self::Element, rhs: &Self::Element) -> Option<Self::Element> {
534            self.get_delegate().checked_div_min(lhs, rhs)
535        }
536
537        fn create_elimination_matrix(
538            &self,
539            a: &Self::Element,
540            b: &Self::Element,
541        ) -> ([Self::Element; 4], Self::Element) {
542            assert_eq!(9, *a);
543            assert_eq!(15, *b);
544            assert_eq!(3, self.add(self.mul(42, *a), self.mul(2, *b)));
545            assert_eq!(0, self.add(self.mul(5, *a), self.mul(3, *b)));
546            return ([42, 2, 10, 6], 3);
547        }
548    }
549
550    let ring = RingValue::from(TestRing);
551    let A = OwnedMatrix::new(vec![9, 0, 15, 3], 2);
552    assert_el_eq!(
553        ring,
554        27,
555        determinant_using_pre_smith(ring, A.clone_matrix(&ring).data_mut(), Global)
556    );
557}
558
559#[test]
560fn test_kernel_basis() {
561    let ring = StaticRing::<i64>::RING;
562    let mut A = OwnedMatrix::identity(2, 2, ring);
563    assert_matrix_eq!(ring, [[], []], kernel_basis_using_pre_smith(ring, A.data_mut(), Global));
564
565    let mut A = OwnedMatrix::zero(2, 2, ring);
566    assert_matrix_eq!(
567        ring,
568        [[1, 0], [0, 1]],
569        kernel_basis_using_pre_smith(ring, A.data_mut(), Global)
570    );
571
572    let A = OwnedMatrix::new(vec![1, 1, 2, 3, 2, 1], 3);
573    let B = kernel_basis_using_pre_smith(ring, A.clone_matrix(ring).data_mut(), Global);
574    assert_eq!(1, B.col_count());
575    assert!(!ring.is_zero(B.at(0, 0)));
576    let mut product = OwnedMatrix::zero(2, 1, ring);
577    STANDARD_MATMUL.matmul(
578        TransposableSubmatrix::from(A.data()),
579        TransposableSubmatrix::from(B.data()),
580        TransposableSubmatrixMut::from(product.data_mut()),
581        ring,
582    );
583    assert_matrix_eq!(ring, [[0], [0]], product);
584
585    let A = OwnedMatrix::new(vec![1, 1, 1, 1, 1, 1], 2);
586    let B = kernel_basis_using_pre_smith(ring, A.clone_matrix(ring).data_mut(), Global);
587    assert_eq!(1, B.col_count());
588    assert!(!ring.is_zero(B.at(0, 0)));
589    let mut product = OwnedMatrix::zero(3, 1, ring);
590    STANDARD_MATMUL.matmul(
591        TransposableSubmatrix::from(A.data()),
592        TransposableSubmatrix::from(B.data()),
593        TransposableSubmatrixMut::from(product.data_mut()),
594        ring,
595    );
596    assert_matrix_eq!(ring, [[0], [0], [0]], product);
597
598    let ring = Zn::new(6);
599    let A = OwnedMatrix::new(vec![ring.int_hom().map(2)], 1);
600    let B = kernel_basis_using_pre_smith(ring, A.clone_matrix(ring).data_mut(), Global);
601    assert_eq!(1, B.col_count());
602    assert_matrix_eq!(ring, [[ring.int_hom().map(3)]], B);
603}
604
605#[test]
606#[ignore]
607fn time_solve_right_using_pre_smith_galois_field() {
608    let n = 100;
609    let base_field = Zn::new(257).as_field().ok().unwrap();
610    let allocator = feanor_mempool::AllocRc(Rc::new(feanor_mempool::dynsize::DynLayoutMempool::new_global(
611        Alignment::of::<u64>(),
612    )));
613    let field = GaloisField::new_with_convolution(base_field, 21, allocator, STANDARD_CONVOLUTION);
614    let matrix = OwnedMatrix::from_fn(n, n, |i, j| {
615        field.pow(field.int_hom().mul_map(field.canonical_gen(), i as i32 + 1), j)
616    });
617
618    let mut inv = OwnedMatrix::zero(n, n, &field);
619    let mut copy = matrix.clone_matrix(&field);
620    let start = Instant::now();
621    solve_right_using_pre_smith(
622        &field,
623        copy.data_mut(),
624        OwnedMatrix::identity(n, n, &field).data_mut(),
625        inv.data_mut(),
626        Global,
627    )
628    .assert_solved();
629    let end = Instant::now();
630    assert_el_eq!(
631        &field,
632        field.one(),
633        <_ as ComputeInnerProduct>::inner_product_ref(
634            field.get_ring(),
635            inv.data().col_at(4).as_iter().zip(matrix.data().row_at(4).as_iter())
636        )
637    );
638
639    println!("total: {} us", (end - start).as_micros());
640}
641
642#[test]
643#[ignore]
644fn time_solve_right_using_extension() {
645    let n = 126;
646    let base_field = Zn::new(257).as_field().ok().unwrap();
647    let allocator = feanor_mempool::AllocRc(Rc::new(feanor_mempool::dynsize::DynLayoutMempool::new_global(
648        Alignment::of::<u64>(),
649    )));
650    let field = GaloisField::new_with_convolution(base_field, 21, allocator, STANDARD_CONVOLUTION);
651    let matrix = OwnedMatrix::from_fn(n, n, |i, j| {
652        field.pow(field.int_hom().mul_map(field.canonical_gen(), i as i32 + 1), j)
653    });
654
655    let mut inv = OwnedMatrix::zero(n, n, &field);
656    let mut copy = matrix.clone_matrix(&field);
657    let start = Instant::now();
658    solve_right_over_extension(
659        &field,
660        copy.data_mut(),
661        OwnedMatrix::identity(n, n, &field).data_mut(),
662        inv.data_mut(),
663        Global,
664    )
665    .assert_solved();
666    let end = Instant::now();
667    assert_el_eq!(
668        &field,
669        field.one(),
670        <_ as ComputeInnerProduct>::inner_product_ref(
671            field.get_ring(),
672            inv.data().col_at(4).as_iter().zip(matrix.data().row_at(4).as_iter())
673        )
674    );
675
676    println!("total: {} us", (end - start).as_micros());
677}
678
679#[bench]
680fn bench_solve_right_using_pre_smith_galois_field(bencher: &mut Bencher) {
681    let base_field = Zn::new(257).as_field().ok().unwrap();
682    let allocator = feanor_mempool::AllocRc(Rc::new(feanor_mempool::dynsize::DynLayoutMempool::new_global(
683        Alignment::of::<u64>(),
684    )));
685    let field = GaloisField::create(
686        FreeAlgebraImpl::new_with_convolution(
687            base_field,
688            5,
689            [base_field.int_hom().map(3), base_field.int_hom().map(-4)],
690            "x",
691            allocator,
692            STANDARD_CONVOLUTION,
693        )
694        .as_field()
695        .ok()
696        .unwrap(),
697    );
698    let matrix = OwnedMatrix::from_fn(10, 10, |i, j| {
699        field.pow(field.int_hom().mul_map(field.canonical_gen(), i as i32 + 1), j)
700    });
701    bencher.iter(|| {
702        let mut inv = OwnedMatrix::zero(10, 10, &field);
703        let mut copy = matrix.clone_matrix(&field);
704        solve_right_using_pre_smith(
705            &field,
706            copy.data_mut(),
707            OwnedMatrix::identity(10, 10, &field).data_mut(),
708            inv.data_mut(),
709            Global,
710        )
711        .assert_solved();
712        assert_el_eq!(
713            &field,
714            field.one(),
715            <_ as ComputeInnerProduct>::inner_product_ref(
716                field.get_ring(),
717                inv.data().col_at(4).as_iter().zip(matrix.data().row_at(4).as_iter())
718            )
719        );
720    });
721}
feanor_math/algorithms/linsolve/smith.rs

feanor_math/algorithms/linsolve/
smith.rs
