feanor_math/algorithms/linsolve/
smith.rs

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

feanor_math/algorithms/linsolve/
smith.rs
