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;
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#[stability::unstable(feature = "enable")]
87pub 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
88    where R: RingStore + Copy,
89        R::Type: PrincipalIdealRing,
90        V1: AsPointerToSlice<El<R>>, V2: AsPointerToSlice<El<R>>, V3: AsPointerToSlice<El<R>>,
91        A: Allocator
92{
93    assert_eq!(lhs.row_count(), rhs.row_count());
94    assert_eq!(lhs.col_count(), out.row_count());
95    assert_eq!(rhs.col_count(), out.col_count());
96
97    let mut R = TransformList::new(lhs.col_count());
98    pre_smith(ring, &mut TransformRows(rhs.reborrow(), ring.get_ring()), &mut R, lhs.reborrow());
99
100    for i in out.row_count()..rhs.row_count() {
101        for j in 0..rhs.col_count() {
102            if !ring.is_zero(rhs.at(i, j)) {
103                return SolveResult::NoSolution;
104            }
105        }
106    }
107    // the value of out[lhs.row_count().., ..] is irrelevant, since lhs
108    // is zero in these places anyway. Thus we just leave it unchanged
109
110    for i in 0..min(lhs.row_count(), lhs.col_count()) {
111        let pivot = lhs.at(i, i);
112        for j in 0..rhs.col_count() {
113            if let Some(quo) = ring.checked_left_div(rhs.at(i, j), pivot) {
114                *out.at_mut(i, j) = quo;
115            } else {
116                return SolveResult::NoSolution;
117            }
118        }
119    }
120
121    R.replay_transposed(ring, TransformRows(out, ring.get_ring()));
122    return SolveResult::FoundSomeSolution;
123}
124
125#[stability::unstable(feature = "enable")]
126pub fn determinant_using_pre_smith<R, V, A>(ring: R, mut matrix: SubmatrixMut<V, El<R>>, _allocator: A) -> El<R>
127    where R: RingStore + Copy,
128        R::Type: PrincipalIdealRing,
129        V:AsPointerToSlice<El<R>>,
130        A: Allocator
131{
132    assert_eq!(matrix.row_count(), matrix.col_count());
133    let mut unit_part_rows = ring.one();
134    let mut unit_part_cols = ring.one();
135    pre_smith(ring, &mut DetUnit { current_unit: &mut unit_part_rows }, &mut DetUnit { current_unit: &mut unit_part_cols }, matrix.reborrow());
136    return ring.checked_div(
137        &ring.prod((0..matrix.row_count()).map(|i| ring.clone_el(matrix.at(i, i)))),
138        &ring.prod([unit_part_rows, unit_part_cols])
139    ).unwrap();
140}
141
142struct DetUnit<'a, R: ?Sized + RingBase> {
143    current_unit: &'a mut R::Element
144}
145
146impl<'a, R> TransformTarget<R> for DetUnit<'a, R>
147    where R: ?Sized + RingBase
148{
149    fn subtract<S: Copy + RingStore<Type = R>>(&mut self, _ring: S, _src: usize, _dst: usize, _factor: &<R as RingBase>::Element) {
150        // determinant does not change
151    }
152
153    fn swap<S: Copy + RingStore<Type = R>>(&mut self, ring: S, _i: usize, _j: usize) {
154        ring.negate_inplace(&mut self.current_unit)
155    }
156
157    fn transform<S: Copy + RingStore<Type = R>>(&mut self, ring: S, _i: usize, _j: usize, transform: &[<R as RingBase>::Element; 4]) {
158        let unit = ring.sub(ring.mul_ref(&transform[0], &transform[3]), ring.mul_ref(&transform[1], &transform[2]));
159        ring.mul_assign(&mut self.current_unit, unit);
160    }
161}
162
163#[cfg(test)]
164use crate::primitive_int::StaticRing;
165#[cfg(test)]
166use crate::rings::extension::extension_impl::FreeAlgebraImpl;
167#[cfg(test)]
168use crate::rings::extension::galois_field::GaloisField;
169#[cfg(test)]
170use crate::rings::extension::FreeAlgebraStore;
171#[cfg(test)]
172use crate::rings::zn::zn_64::Zn;
173#[cfg(test)]
174use crate::rings::zn::zn_static;
175#[cfg(test)]
176use crate::assert_matrix_eq;
177#[cfg(test)]
178use crate::rings::zn::ZnRingStore;
179#[cfg(test)]
180use crate::seq::VectorView;
181#[cfg(test)]
182use crate::algorithms::matmul::ComputeInnerProduct;
183#[cfg(test)]
184use std::alloc::Global;
185#[cfg(test)]
186use test::Bencher;
187#[cfg(test)]
188use crate::homomorphism::Homomorphism;
189#[cfg(test)]
190use crate::algorithms::linsolve::LinSolveRing;
191#[cfg(test)]
192use std::ptr::Alignment;
193#[cfg(test)]
194use std::rc::Rc;
195#[cfg(test)]
196use std::time::Instant;
197#[cfg(test)]
198use crate::algorithms::convolution::STANDARD_CONVOLUTION;
199#[cfg(test)]
200use crate::algorithms::linsolve::extension::solve_right_over_extension;
201#[cfg(test)]
202use crate::algorithms::matmul::*;
203
204#[cfg(test)]
205fn multiply<'a, R: RingStore, V: AsPointerToSlice<El<R>>, I: IntoIterator<Item = Submatrix<'a, V, El<R>>>>(matrices: I, ring: R) -> OwnedMatrix<El<R>>
206    where R::Type: 'a,
207        V: 'a
208{
209    let mut it = matrices.into_iter();
210    let fst = it.next().unwrap();
211    let snd = it.next().unwrap();
212    let mut new_result = OwnedMatrix::zero(fst.row_count(), snd.col_count(), &ring);
213    STANDARD_MATMUL.matmul(TransposableSubmatrix::from(fst), TransposableSubmatrix::from(snd), TransposableSubmatrixMut::from(new_result.data_mut()), &ring);
214    let mut result = new_result;
215
216    for m in it {
217        let mut new_result = OwnedMatrix::zero(result.row_count(), m.col_count(), &ring);
218        STANDARD_MATMUL.matmul(TransposableSubmatrix::from(result.data()), TransposableSubmatrix::from(m), TransposableSubmatrixMut::from(new_result.data_mut()), &ring);
219        result = new_result;
220    }
221    return result;
222}
223
224#[test]
225fn test_smith_integers() {
226    let ring = StaticRing::<i64>::RING;
227    let mut A = OwnedMatrix::new(
228        vec![ 1, 2, 3, 4, 
229                    2, 3, 4, 5,
230                    3, 4, 5, 6 ], 
231        4
232    );
233    let original_A = A.clone_matrix(&ring);
234    let mut L: OwnedMatrix<i64> = OwnedMatrix::identity(3, 3, StaticRing::<i64>::RING);
235    let mut R: OwnedMatrix<i64> = OwnedMatrix::identity(4, 4, StaticRing::<i64>::RING);
236    pre_smith(ring, &mut TransformRows(L.data_mut(), ring.get_ring()), &mut TransformCols(R.data_mut(), ring.get_ring()), A.data_mut());
237    
238    assert_matrix_eq!(&ring, &[
239        [1, 0, 0, 0],
240        [0,-1, 0, 0],
241        [0, 0, 0, 0]], &A);
242
243    assert_matrix_eq!(&ring, &multiply([L.data(), original_A.data(), R.data()], ring), &A);
244}
245
246#[test]
247fn test_smith_zn() {
248    let ring = zn_static::Zn::<45>::RING;
249    let mut A = OwnedMatrix::new(
250        vec![ 8, 3, 5, 8,
251                    0, 9, 0, 9,
252                    5, 9, 5, 14,
253                    8, 3, 5, 23,
254                    3,39, 0, 39 ],
255        4
256    );
257    let original_A = A.clone_matrix(&ring);
258    let mut L: OwnedMatrix<u64> = OwnedMatrix::identity(5, 5, ring);
259    let mut R: OwnedMatrix<u64> = OwnedMatrix::identity(4, 4, ring);
260    pre_smith(ring, &mut TransformRows(L.data_mut(), ring.get_ring()), &mut TransformCols(R.data_mut(), ring.get_ring()), A.data_mut());
261
262    assert_matrix_eq!(&ring, &[
263        [8, 0, 0, 0],
264        [0, 3, 0, 0],
265        [0, 0, 0, 0], 
266        [0, 0, 0, 15],
267        [0, 0, 0, 0]], &A);
268        
269    assert_matrix_eq!(&ring, &multiply([L.data(), original_A.data(), R.data()], ring), &A);
270}
271
272#[test]
273fn test_solve_zn() {
274    let ring = zn_static::Zn::<45>::RING;
275    let A = OwnedMatrix::new(
276        vec![ 8, 3, 5, 8,
277                    0, 9, 0, 9,
278                    5, 9, 5, 14,
279                    8, 3, 5, 23,
280                    3,39, 0, 39 ],
281        4
282    );
283    let B = OwnedMatrix::new(
284        vec![11, 43, 10, 22,
285                   18,  9, 27, 27,
286                    8, 34,  7, 22,
287                   41, 13, 40, 37,
288                    3,  9,  3,  0],
289        4
290    );
291    let mut solution: OwnedMatrix<_> = OwnedMatrix::zero(4, 4, ring);
292    ring.get_ring().solve_right(A.clone_matrix(ring).data_mut(), B.clone_matrix(ring).data_mut(), solution.data_mut(), Global).assert_solved();
293
294    assert_matrix_eq!(&ring, &multiply([A.data(), solution.data()], ring), &B);
295}
296
297#[test]
298fn test_solve_int() {
299    let ring = StaticRing::<i64>::RING;
300    let A = OwnedMatrix::new(
301        vec![3, 6, 2, 0, 4, 7,
302                   5, 5, 4, 5, 5, 5],
303        6
304    );
305    let B: OwnedMatrix<i64> = OwnedMatrix::identity(2, 2, ring);
306    let mut solution: OwnedMatrix<i64> = OwnedMatrix::zero(6, 2, ring);
307    ring.get_ring().solve_right(A.clone_matrix(ring).data_mut(), B.clone_matrix(ring).data_mut(), solution.data_mut(), Global).assert_solved();
308
309    assert_matrix_eq!(&ring, &multiply([A.data(), solution.data()], &ring), &B);
310}
311
312#[test]
313fn test_large() {
314    let ring = zn_static::Zn::<16>::RING;
315    let data_A = [
316        [0, 0, 0, 0, 0, 0, 0, 0,11, 0, 0],
317        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
318        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0,10],
319        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
320        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 8],
321        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
322    ];
323    let mut A: OwnedMatrix<u64> = OwnedMatrix::zero(6, 11, &ring);
324    for i in 0..6 {
325        for j in 0..11 {
326            *A.at_mut(i, j) = data_A[i][j];
327        }
328    }
329    let mut solution: OwnedMatrix<_> = OwnedMatrix::zero(11, 11, ring);
330    assert!(ring.get_ring().solve_right(A.clone_matrix(&ring).data_mut(), A.data_mut(), solution.data_mut(), Global).is_solved());
331}
332
333#[test]
334fn test_determinant() {
335    let ring = StaticRing::<i64>::RING;
336    let A = OwnedMatrix::new(
337        vec![1, 0, 3, 
338                   2, 1, 0, 
339                   9, 8, 7],
340        3
341    );
342    assert_el_eq!(ring, (7 + 48 - 27), determinant_using_pre_smith(ring, A.clone_matrix(&ring).data_mut(), Global));
343    
344    // we need a ring that has units of order > 2 to test whether an inversion is necessary for
345    // the accumulated determinant units
346    #[derive(PartialEq, Clone, Copy)]
347    struct TestRing;
348    use crate::delegate::DelegateRing;
349    impl DelegateRing for TestRing {
350        type Base = zn_static::ZnBase<45, false>;
351        type Element = u64;
352
353        fn get_delegate(&self) -> &Self::Base { zn_static::Zn::RING.get_ring() }
354        fn delegate(&self, el: Self::Element) -> <Self::Base as RingBase>::Element { el }
355        fn delegate_mut<'a>(&self, el: &'a mut Self::Element) -> &'a mut <Self::Base as RingBase>::Element { el }
356        fn delegate_ref<'a>(&self, el: &'a Self::Element) -> &'a <Self::Base as RingBase>::Element { el }
357        fn rev_delegate(&self, el: <Self::Base as RingBase>::Element) -> Self::Element { el }
358    }
359    impl PrincipalIdealRing for TestRing {
360
361        fn extended_ideal_gen(&self, lhs: &Self::Element, rhs: &Self::Element) -> (Self::Element, Self::Element, Self::Element) {
362            self.get_delegate().extended_ideal_gen(lhs, rhs)
363        }
364
365        fn checked_div_min(&self, lhs: &Self::Element, rhs: &Self::Element) -> Option<Self::Element> {
366            self.get_delegate().checked_div_min(lhs, rhs)
367        }
368
369        fn create_elimination_matrix(&self, a: &Self::Element, b: &Self::Element) -> ([Self::Element; 4], Self::Element) {
370            assert_eq!(9, *a);
371            assert_eq!(15, *b);
372            assert_eq!(3, self.add(self.mul(42, *a), self.mul(2, *b)));
373            assert_eq!(0, self.add(self.mul(5, *a), self.mul(3, *b)));
374            return ([42, 2, 10, 6], 3);
375        }
376    }
377
378    let ring = RingValue::from(TestRing);
379    let A = OwnedMatrix::new(
380        vec![ 9, 0, 
381                   15, 3],
382        2
383    );
384    assert_el_eq!(ring, 27, determinant_using_pre_smith(ring, A.clone_matrix(&ring).data_mut(), Global));
385}
386
387#[test]
388#[ignore]
389fn time_solve_right_using_pre_smith_galois_field() {
390    let n = 100;
391    let base_field = Zn::new(257).as_field().ok().unwrap();
392    let allocator = feanor_mempool::AllocRc(Rc::new(feanor_mempool::dynsize::DynLayoutMempool::new_global(Alignment::of::<u64>())));
393    let field = GaloisField::new_with_convolution(base_field, 21, allocator, STANDARD_CONVOLUTION);
394    let matrix = OwnedMatrix::from_fn(n, n, |i, j| field.pow(field.int_hom().mul_map(field.canonical_gen(), i as i32 + 1), j));
395    
396    let mut inv = OwnedMatrix::zero(n, n, &field);
397    let mut copy = matrix.clone_matrix(&field);
398    let start = Instant::now();
399    solve_right_using_pre_smith(&field, copy.data_mut(), OwnedMatrix::identity(n, n, &field).data_mut(), inv.data_mut(), Global).assert_solved();
400    let end = Instant::now();
401    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())));
402    
403    println!("total: {} us", (end - start).as_micros());
404}
405
406#[test]
407#[ignore]
408fn time_solve_right_using_extension() {
409    let n = 126;
410    let base_field = Zn::new(257).as_field().ok().unwrap();
411    let allocator = feanor_mempool::AllocRc(Rc::new(feanor_mempool::dynsize::DynLayoutMempool::new_global(Alignment::of::<u64>())));
412    let field = GaloisField::new_with_convolution(base_field, 21, allocator, STANDARD_CONVOLUTION);
413    let matrix = OwnedMatrix::from_fn(n, n, |i, j| field.pow(field.int_hom().mul_map(field.canonical_gen(), i as i32 + 1), j));
414    
415    let mut inv = OwnedMatrix::zero(n, n, &field);
416    let mut copy = matrix.clone_matrix(&field);
417    let start = Instant::now();
418    solve_right_over_extension(&field, copy.data_mut(), OwnedMatrix::identity(n, n, &field).data_mut(), inv.data_mut(), Global).assert_solved();
419    let end = Instant::now();
420    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())));
421
422    println!("total: {} us", (end - start).as_micros());
423}
424
425#[bench]
426fn bench_solve_right_using_pre_smith_galois_field(bencher: &mut Bencher) {
427    let base_field = Zn::new(257).as_field().ok().unwrap();
428    let allocator = feanor_mempool::AllocRc(Rc::new(feanor_mempool::dynsize::DynLayoutMempool::new_global(Alignment::of::<u64>())));
429    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());
430    let matrix = OwnedMatrix::from_fn(10, 10, |i, j| field.pow(field.int_hom().mul_map(field.canonical_gen(), i as i32 + 1), j));
431    bencher.iter(|| {
432        let mut inv = OwnedMatrix::zero(10, 10, &field);
433        let mut copy = matrix.clone_matrix(&field);
434        solve_right_using_pre_smith(&field, copy.data_mut(), OwnedMatrix::identity(10, 10, &field).data_mut(), inv.data_mut(), Global).assert_solved();
435        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())));
436    });
437}
feanor_math/algorithms/linsolve/smith.rs

feanor_math/algorithms/linsolve/
smith.rs
