feanor_math/algorithms/
lll.rs

1use crate::algorithms::matmul::MatmulAlgorithm;
2use crate::algorithms::matmul::STANDARD_MATMUL;
3use crate::field::*;
4use crate::integer::*;
5use crate::homomorphism::*;
6use crate::matrix::*;
7use crate::matrix::transform::TransformTarget;
8use crate::primitive_int::*;
9use crate::rings::float_real::Real64;
10use crate::rings::float_real::Real64Base;
11use crate::rings::rational::*;
12use crate::ordered::{OrderedRing, OrderedRingStore};
13use crate::ring::*;
14
15use std::alloc::Allocator;
16use std::cmp::max;
17use std::marker::PhantomData;
18
19///
20/// Trait for (possibly approximations to) the real numbers that are used to keep
21/// an estimate of the size and orthogonality of vectors during executions of LLL.
22/// 
23/// Errors caused by approximation might reduce the quality of the LLL-reduced basis,
24/// but won't give any other correctness errors.
25/// 
26#[stability::unstable(feature = "enable")]
27pub trait LLLRealField<I>: OrderedRing + Field
28    where I: ?Sized + IntegerRing
29{
30    fn from_integer(&self, x: I::Element, ZZ: &I) -> Self::Element;
31    fn round_to_integer(&self, x: &Self::Element, ZZ: &I) -> I::Element;
32}
33
34impl<I, J> LLLRealField<J> for RationalFieldBase<I>
35    where I: IntegerRingStore,
36        I::Type: IntegerRing,
37        J: ?Sized + IntegerRing
38{
39    fn from_integer(&self, x: J::Element, ZZ: &J) -> Self::Element {
40        RingRef::new(self).inclusion().map(int_cast(x, self.base_ring(), RingRef::new(ZZ)))
41    }
42
43    fn round_to_integer(&self, x: &Self::Element, ZZ: &J) -> J::Element {
44        int_cast(self.base_ring().rounded_div(self.base_ring().clone_el(self.num(x)), self.den(x)), RingRef::new(ZZ), self.base_ring())
45    }
46}
47
48impl<J> LLLRealField<J> for Real64Base
49    where J: ?Sized + IntegerRing
50{
51    fn from_integer(&self, x: J::Element, ZZ: &J) -> Self::Element {
52        ZZ.to_float_approx(&x)
53    }
54
55    fn round_to_integer(&self, x: &Self::Element, ZZ: &J) -> J::Element {
56        int_cast(x.round() as i64, RingRef::new(ZZ), StaticRing::<i64>::RING)
57    }
58}
59
60fn size_reduce<R, I, V, T>(ring: R, int_ring: I, mut target: SubmatrixMut<V, El<R>>, target_j: usize, matrix: Submatrix<V, El<R>>, col_ops: &mut T)
61    where R: RingStore,
62        R::Type: LLLRealField<I::Type>,
63        I: IntegerRingStore + Copy,
64        I::Type: IntegerRing,
65        V: AsPointerToSlice<El<R>>,
66        T: TransformTarget<I::Type>
67{
68    for j in (0..matrix.col_count()).rev() {
69        let factor = ring.get_ring().round_to_integer(target.as_const().at(j, 0), int_ring.get_ring());
70        col_ops.subtract(int_ring, j, target_j, &factor);
71        let factor = ring.get_ring().from_integer(factor, int_ring.get_ring());
72        ring.sub_assign_ref(target.at_mut(j, 0), &factor);
73        for k in 0..j {
74            ring.sub_assign(target.at_mut(k, 0), ring.mul_ref(matrix.at(k, j), &factor));
75        }
76    }
77}
78
79///
80/// Updates the "gso"-matrix resulting from swapping cols i and i + 1.
81/// 
82/// This uses explicit formulas, in particular if `b_i` and `b_(i + 1)`
83/// represent the vectors, `b_i*` and `b_(i + 1)*` represent their GS-orthogonalizations
84/// and the corresponding `b'` values are the ones after swapping, then we find
85/// then this gives
86/// ```text
87///   b'_i = b_(i + 1)
88///   b'_(i + 1) = b_i
89/// 
90///   b'_i* = b_(i + 1)* + mu b_i*
91///   b'_(i + 1) = (1 - gamma^2 mu^2) b_i* - mu * gamma^2 b_(i + 1)*
92///     where gamma^2 = |b_i*|^2 / |b'_i*|^2
93///   mu' = gamma^2 mu
94/// ```
95/// 
96fn swap_gso_cols<R, V>(ring: R, mut gso: SubmatrixMut<V, El<R>>, i: usize, j: usize)
97    where R: RingStore,
98        R::Type: OrderedRing + Field,
99        V: AsPointerToSlice<El<R>>
100{
101    assert!(j == i + 1);
102
103    let col_count = gso.col_count();
104
105    // swap the columns
106    let (mut col_i, mut col_i1) = gso.reborrow().restrict_cols(i..(i + 2)).split_cols(0..1, 1..2);
107    for k in 0..i {
108        std::mem::swap(col_i.at_mut(k, 0), col_i1.at_mut(k, 0));
109    }
110
111    // re-orthogonalize the triangle `i..(i + 2) x i..(i + 2)`
112
113    // | b_i* |^2
114    let bi_star_norm_sqr = ring.clone_el(gso.at(i, i));
115    // | b_(i + 1)* |^2
116    let bi1_star_norm_sqr = ring.clone_el(gso.at(i + 1, i + 1));
117    // mu_(i + 1)i = <b_(i + 1), bi*> / <bi*, bi*>
118    let mu = ring.clone_el(gso.at(i, i + 1));
119    let mu_sqr = ring.pow(ring.clone_el(&mu), 2);
120
121    let new_bi_star_norm_sqr = ring.add_ref_fst(&bi1_star_norm_sqr, ring.mul_ref(&mu_sqr, &bi_star_norm_sqr));
122    // `|b_i*|^2 / |bnew_i*|^2`
123    let gamma_sqr = ring.div(&bi_star_norm_sqr, &new_bi_star_norm_sqr);
124    let new_bi1_star_norm_sqr = ring.mul_ref(&gamma_sqr, &bi1_star_norm_sqr);
125    let new_mu = ring.mul_ref(&gamma_sqr, &mu);
126
127    // we now update the `mu_ki` resp. `mu_k(i + 1)` by a linear transform
128    let lin_transform_muki = [ring.mul_ref(&gamma_sqr, &mu), ring.sub(ring.one(), ring.mul_ref(&gamma_sqr, &mu_sqr))];
129    let (mut row_i, mut row_i1) = gso.reborrow().restrict_rows(i..(i + 2)).split_rows(0..1, 1..2);
130    for k in (i + 2)..col_count {
131        let mu_ki = ring.clone_el(row_i.at(0, k));
132        std::mem::swap(row_i.at_mut(0, k), row_i1.at_mut(0, k));
133        ring.sub_assign(row_i1.at_mut(0, k), ring.mul_ref(&mu, row_i.at(0, k)));
134        ring.mul_assign_ref(row_i.at_mut(0, k), &lin_transform_muki[1]);
135        ring.add_assign(row_i.at_mut(0, k), ring.mul_ref_fst(&lin_transform_muki[0], mu_ki));
136    }
137
138    *gso.at_mut(i, i) = new_bi_star_norm_sqr;
139    *gso.at_mut(i, i + 1) = new_mu;
140    *gso.at_mut(i + 1, i + 1) = new_bi1_star_norm_sqr;
141}
142
143///
144/// gso contains on the diagonal the squared lengths of the GS-orthogonalized basis vectors `|bi*|^2`,
145/// and above it the GS-coefficients `mu_ij = <bi, bj*> / <bj*, bj*>`.
146/// 
147fn lll_base<R, I, V, T>(ring: R, int_ring: I, mut gso: SubmatrixMut<V, El<R>>, mut col_ops: T, delta: &El<R>)
148    where R: RingStore + Copy,
149        R::Type: LLLRealField<I::Type>,
150        I: IntegerRingStore + Copy,
151        I::Type: IntegerRing,
152        V: AsPointerToSlice<El<R>>,
153        T: TransformTarget<I::Type>
154{
155    let mut i = 0;
156    while i + 1 < gso.col_count() {
157        let (target, matrix) = gso.reborrow().split_cols((i + 1)..(i + 2), 0..(i + 1));
158        size_reduce(ring, int_ring, target, i + 1, matrix.as_const(), &mut col_ops);
159        if ring.is_gt(
160            &ring.mul_ref_snd(
161                ring.sub_ref_fst(delta, ring.mul_ref(gso.as_const().at(i, i + 1), gso.as_const().at(i, i + 1))),
162                gso.as_const().at(i, i)
163            ),
164            gso.as_const().at(i + 1, i + 1)
165        ) {
166            col_ops.swap(int_ring, i, i + 1);
167            swap_gso_cols(ring, gso.reborrow(), i, i + 1);
168            i = max(i, 1) - 1;
169        } else {
170            i += 1;
171        }
172    }
173}
174
175///
176/// Computes the LDL-decomposition of the given matrix, i.e. writes it as
177/// a product `L * D * L^T`, where `D` is diagonal and `L` is lower triangle.
178/// 
179/// Currently this requires that the input matrix is invertible (or in the 
180/// floating point case that no eigenvalues are very small in absolute value).
181/// 
182/// `D` is returned on the diagonal of the matrix, and `L^T` is returned in
183/// the upper triangle of the matrix.
184/// 
185fn ldl<R, V>(ring: R, mut matrix: SubmatrixMut<V, El<R>>)
186    where R: RingStore,
187        R::Type: Field, 
188        V: AsPointerToSlice<El<R>>
189{
190    // only the upper triangle part of matrix is used
191    assert_eq!(matrix.row_count(), matrix.col_count());
192    let n = matrix.row_count();
193    for i in 0..n {
194        let pivot = ring.clone_el(matrix.at(i, i));
195        let pivot_inv = ring.div(&ring.one(), matrix.at(i, i));
196        for j in (i + 1)..n {
197            ring.mul_assign_ref(matrix.at_mut(i, j), &pivot_inv);
198        }
199        for k in (i + 1)..n {
200            for l in k..n {
201                let subtract = ring.mul_ref_snd(ring.mul_ref(matrix.as_const().at(i, k), matrix.as_const().at(i, l)), &pivot);
202                ring.sub_assign(matrix.at_mut(k, l), subtract);
203            }
204        }
205    }
206}
207
208///
209/// LLL-reduces the given matrix, i.e. transforms `B` into `B'` such that
210/// `B' = BU` with `U in GL(Z)` and the columns of `B'` are "short".
211/// 
212/// The exact restrictions imposed on `B'` are that its columns `b1, ..., bn`
213/// are "LLL-reduced". This means
214///  - (size-reduced) `|<bi,bj*>| < <bj*,bj*> / 2` whenever `i > j`
215///  - (Lovasz-condition) `delta |b(k - 1)|^2 - <bk, b(k - 1)*>^2 / |b(k - 1)|^2 <= |bk*|^2`
216/// 
217/// Here the `bi*` refer to the Gram-Schmidt orthogonalization of the `bi`, and
218/// `<.,.>` is the inner product induced by `quadratic_form`.
219/// 
220/// # Internal computations with floating point numbers
221/// 
222/// For efficiency reasons, this function performs computations with floating point
223/// number internally. It is ensured that all operations are unimodular, so the result
224/// `B'` will always satisfy `B' = BU`. However (in particular if the conversion between
225/// `I` and `f64` is not implemented with high quality), the vectors in `B'` might be
226/// somewhat longer than explained above.
227/// 
228#[stability::unstable(feature = "enable")]
229pub fn lll_float<I, V1, V2, A>(ring: I, quadratic_form: Submatrix<V1, El<Real64>>, matrix: SubmatrixMut<V2, El<I>>, delta: f64, allocator: A)
230    where I: IntegerRingStore,
231        I::Type: IntegerRing,
232        V1: AsPointerToSlice<El<Real64>>,
233        V2: AsPointerToSlice<El<I>>,
234        A: Allocator
235{
236    assert!(delta < 1.);
237    assert!(delta > 0.25);
238    assert_eq!(quadratic_form.row_count(), quadratic_form.col_count());
239    assert_eq!(matrix.col_count(), quadratic_form.col_count());
240
241    let lll_reals = Real64::RING;
242    let hom = lll_reals.can_hom(&ring).unwrap();
243
244    let n = matrix.col_count();
245    let mut gso: OwnedMatrix<f64, &A> = OwnedMatrix::zero_in(n, n, &lll_reals, &allocator);
246    let mut tmp: OwnedMatrix<f64, &A> = OwnedMatrix::zero_in(n, n, &lll_reals, &allocator);
247    let matrix_RR: OwnedMatrix<f64, &A> = OwnedMatrix::from_fn_in(matrix.row_count(), matrix.col_count(), |i, j| hom.map_ref(matrix.at(i, j)), &allocator);
248    STANDARD_MATMUL.matmul(TransposableSubmatrix::from(matrix_RR.data()).transpose(), TransposableSubmatrix::from(quadratic_form), TransposableSubmatrixMut::from(tmp.data_mut()), lll_reals);
249    STANDARD_MATMUL.matmul(TransposableSubmatrix::from(tmp.data()), TransposableSubmatrix::from(matrix_RR.data()), TransposableSubmatrixMut::from(gso.data_mut()), lll_reals);
250
251    ldl(&lll_reals, gso.data_mut());
252    lll_base::<_, _, _, TransformLatticeBasis<I::Type, I::Type, _, _>>(&lll_reals, &ring, gso.data_mut(), TransformLatticeBasis { basis: matrix, hom: ring.identity(), int_ring: PhantomData }, &delta);
253}
254
255///
256/// LLL-reduces the given matrix, i.e. transforms `B` into `B'` such that
257/// `B' = BU` with `U in GL(Z)` and the columns of `B'` are "short".
258/// 
259/// The exact restrictions imposed on `B'` are that its columns `b1, ..., bn`
260/// are "LLL-reduced". This means
261///  - (size-reduced) `|<bi,bj*>| < <bj*,bj*> / 2` whenever `i > j`
262///  - (Lovasz-condition) `delta |b(k - 1)|^2 - <bk, b(k - 1)*>^2 / |b(k - 1)|^2 <= |bk*|^2`
263/// Here the `bi*` refer to the Gram-Schmidt orthogonalization of the `bi`.
264/// 
265/// # Internal computations with floating point numbers
266/// 
267/// The LLL algorithm has to handle fractional values internally, which can get huge
268/// denominators. Since only their size is relevant, it is usually much better to work
269/// with floating point numbers instead. This is done as a precomputation in this function,
270/// and as the only version in [`lll_float()`].
271/// 
272/// Despite precomputing the floating-point LLL in this function, we still perform the
273/// computation with rationals here, hence the runtime is likely to be vastly longer than
274/// the one of [`lll_float()`]. Since it is very unlikely that you need the exact guarantees
275/// of [`lll_exact()`], but just want a "short" matrix, prefer using [`lll_float()`] instead.
276/// 
277#[stability::unstable(feature = "enable")]
278pub fn lll_exact<I, V, A>(ring: I, mut matrix: SubmatrixMut<V, El<I>>, delta: f64, allocator: A)
279    where I: IntegerRingStore,
280        I::Type: IntegerRing,
281        V: AsPointerToSlice<El<I>>,
282        A: Allocator
283{
284    assert!(delta < 1.);
285    assert!(delta > 0.25);
286    lll_float(&ring, OwnedMatrix::<f64>::identity(matrix.col_count(), matrix.col_count(), Real64::RING).data(), matrix.reborrow(), delta, &allocator);
287
288    let lll_reals = RationalField::new(BigIntRing::RING);
289    let delta_int = ring.from_float_approx(delta * 2f64.powi(20)).unwrap();
290    let delta = lll_reals.div(&lll_reals.get_ring().from_integer(delta_int, ring.get_ring()), &lll_reals.get_ring().from_integer(ring.power_of_two(20), ring.get_ring()));
291
292    let n = matrix.col_count();
293    let mut gso = OwnedMatrix::zero_in(n, n, &lll_reals, &allocator);
294    let hom = lll_reals.inclusion().compose(BigIntRing::RING.can_hom(&ring).unwrap());
295    let matrix_RR: OwnedMatrix<_, &A> = OwnedMatrix::from_fn_in(matrix.row_count(), matrix.col_count(), |i, j| hom.map_ref(matrix.at(i, j)), &allocator);
296    STANDARD_MATMUL.matmul(TransposableSubmatrix::from(matrix_RR.data()).transpose(), TransposableSubmatrix::from(matrix_RR.data()), TransposableSubmatrixMut::from(gso.data_mut()), &lll_reals);
297
298    ldl(&lll_reals, gso.data_mut());
299    lll_base::<_, _, _, TransformLatticeBasis<I::Type, I::Type, _, _>>(&lll_reals, &ring, gso.data_mut(), TransformLatticeBasis { basis: matrix, hom: ring.identity(), int_ring: PhantomData }, &delta);
300}
301
302struct TransformLatticeBasis<'a, R, I, V, H>
303    where R: ?Sized + RingBase,
304        I: ?Sized + IntegerRing,
305        H: Homomorphism<I, R>,
306        V: AsPointerToSlice<R::Element>
307{
308    basis: SubmatrixMut<'a, V, R::Element>,
309    int_ring: PhantomData<I>,
310    hom: H
311}
312
313impl<'a, R, I, V, H> TransformTarget<I> for TransformLatticeBasis<'a, R, I, V, H>
314    where R: ?Sized + RingBase,
315        I: ?Sized + IntegerRing,
316        H: Homomorphism<I, R>,
317        V: AsPointerToSlice<R::Element>
318{
319    fn transform<S: Copy + RingStore<Type = I>>(&mut self, ring: S, i: usize, j: usize, transform: &[I::Element; 4]) {
320        assert!(ring.get_ring() == self.hom.domain().get_ring());
321        assert!(i != j);
322        let ring = self.hom.codomain();
323        for k in 0..self.basis.row_count() {
324            let a = ring.clone_el(self.basis.at(k, i));
325            let b = ring.clone_el(self.basis.at(k, j));
326            *self.basis.at_mut(k, i) = ring.add(self.hom.mul_ref_map(&a, &transform[0]), self.hom.mul_ref_map(&b, &transform[1]));
327            *self.basis.at_mut(k, i) = ring.add(self.hom.mul_ref_snd_map(a, &transform[2]), self.hom.mul_ref_snd_map(b, &transform[3]));
328        }
329    }
330
331    fn subtract<S: Copy + RingStore<Type = I>>(&mut self, ring: S, src: usize, dst: usize, factor: &I::Element) {
332        assert!(ring.get_ring() == self.hom.domain().get_ring());
333        assert!(src != dst);
334        let ring = self.hom.codomain();
335        for k in 0..self.basis.row_count() {
336            let subtract = self.hom.mul_ref_map(self.basis.at(k, src), factor);
337            ring.sub_assign(self.basis.at_mut(k, dst), subtract);
338        }
339    }
340
341    fn swap<S: Copy + RingStore<Type = I>>(&mut self, ring: S, i: usize, j: usize) {
342        assert!(ring.get_ring() == self.hom.domain().get_ring());
343        if i == j {
344            return;
345        }
346        let col_count = self.basis.col_count();
347        let (mut col_i, mut col_j) = self.basis.reborrow().split_cols(i..(i + 1), j..(j + 1));
348        for k in 0..col_count {
349            std::mem::swap(col_i.at_mut(k, 0), col_j.at_mut(k, 0));
350        }
351    }
352}
353
354#[cfg(test)]
355use crate::seq::*;
356#[cfg(test)]
357use test::Bencher;
358#[cfg(test)]
359use std::alloc::Global;
360#[cfg(test)]
361use crate::assert_matrix_eq;
362
363#[cfg(test)]
364const QQ: RationalField<StaticRing<i64>> = RationalField::new(StaticRing::<i64>::RING);
365
366#[cfg(test)]
367macro_rules! in_QQ {
368    ($hom:expr; $num:literal) => {
369        ($hom).map($num)
370    };
371    ($hom:expr; $num:literal, $den:literal) => {
372        ($hom).codomain().div(&($hom).map($num), &($hom).map($den))
373    };
374    ($([$($num:literal $(/ $den:literal)?),*]),*) => {
375        {
376            let ZZ_to_QQ = QQ.inclusion();
377            [
378                $([$(
379                    in_QQ!(ZZ_to_QQ; $num $(, $den)?)
380                ),*]),*
381            ]
382        }
383    };
384    ($(DerefArray::from([$($num:literal $(/ $den:literal)?),*])),*) => {
385        {
386            let ZZ_to_QQ = QQ.inclusion();
387            [
388                $(DerefArray::from([$(
389                    in_QQ!(ZZ_to_QQ; $num $(, $den)?)
390                ),*])),*
391            ]
392        }
393    };
394}
395
396#[test]
397fn test_ldl() {
398    let mut data = in_QQ![
399        DerefArray::from([1, 2, 1]),
400        DerefArray::from([2, 5, 0]),
401        DerefArray::from([1, 0, 7])
402    ];
403    let mut matrix = SubmatrixMut::<DerefArray<_, 3>, _>::from_2d(&mut data);
404    let mut expected = in_QQ![
405        [1, 2, 1],
406        [0, 1, -2],
407        [0, 0, 2]
408    ];
409    ldl(QQ, matrix.reborrow());
410
411    // only the upper triangle is filled
412    expected[1][0] = *matrix.at(1, 0);
413    expected[2][0] = *matrix.at(2, 0);
414    expected[2][1] = *matrix.at(2, 1);
415
416    assert_matrix_eq!(&QQ, &expected, &matrix);
417}
418
419#[test]
420fn test_swap_gso_cols() {
421    let mut matrix = in_QQ![
422        DerefArray::from([2, 1/2, 2/5]),
423        DerefArray::from([0, 3/2, 1/4]),
424        DerefArray::from([0,   0,   1])
425    ];
426    let expected = in_QQ![
427        [2, 1/2, 31/80],
428        [0, 3/2, 11/40],
429        [0,   0,     1]
430    ];
431    let matrix_view = SubmatrixMut::<DerefArray<_, 3>, _>::from_2d(&mut matrix);
432
433    swap_gso_cols(&QQ, matrix_view, 0, 1);
434
435    assert_matrix_eq!(&QQ, &expected, &matrix);
436}
437
438#[cfg(test)]
439fn norm_squared<V>(col: &Column<V, i64>) -> i64
440    where V: AsPointerToSlice<i64>
441{
442    StaticRing::<i64>::RING.sum((0..col.len()).map(|i| col.at(i) * col.at(i)))
443}
444
445#[cfg(test)]
446fn assert_lattice_isomorphic<V, const N: usize, const M: usize>(lhs: &[DerefArray<i64, M>; N], rhs: &Submatrix<V, i64>)
447    where V: AsPointerToSlice<i64>
448{
449    use crate::algorithms::linsolve::smith;
450
451    assert_eq!(rhs.row_count(), N);
452    assert_eq!(rhs.col_count(), M);
453    let ZZbig = BigIntRing::RING;
454    let mut A: OwnedMatrix<_> = OwnedMatrix::zero(N, M, ZZbig);
455    let mut B: OwnedMatrix<_> = OwnedMatrix::zero(N, M, ZZbig);
456    let int_to_ZZbig: CanHom<&RingValue<StaticRingBase<i64>>, &BigIntRing> = ZZbig.can_hom(&StaticRing::<i64>::RING).unwrap();
457    for i in 0..N {
458        for j in 0..M {
459            *A.at_mut(i, j) = int_to_ZZbig.map(lhs[i][j]);
460            *B.at_mut(i, j) = int_to_ZZbig.map(*rhs.at(i, j));
461        }
462    }
463    let mut U: OwnedMatrix<_> = OwnedMatrix::zero(N, M, ZZbig);
464    assert!(smith::solve_right_using_pre_smith(&ZZbig, A.clone_matrix(&ZZbig).data_mut(), B.clone_matrix(&ZZbig).data_mut(), U.data_mut(), Global).is_solved());
465    assert!(smith::solve_right_using_pre_smith(&ZZbig, B.clone_matrix(&ZZbig).data_mut(), A.clone_matrix(&ZZbig).data_mut(), U.data_mut(), Global).is_solved());
466}
467
468#[test]
469fn test_lll_float_2d() {
470    let ZZ = StaticRing::<i64>::RING;
471    let original = [
472        DerefArray::from([5,   9]),
473        DerefArray::from([11, 20])
474    ];
475    let mut reduced = original;
476    let mut reduced_matrix = SubmatrixMut::<DerefArray<_, 2>, _>::from_2d(&mut reduced);
477    lll_float(&ZZ, OwnedMatrix::<_>::identity(2, 2, Real64::RING).data(), reduced_matrix.reborrow(), 0.9, Global);
478
479    assert_lattice_isomorphic(&original, &reduced_matrix.as_const());
480    assert_eq!(1, norm_squared(&reduced_matrix.as_const().col_at(0)));
481    assert_eq!(1, norm_squared(&reduced_matrix.as_const().col_at(1)));
482
483    let original = [
484        DerefArray::from([10, 8]),
485        DerefArray::from([27, 22])
486    ];
487    let mut reduced = original;
488    let mut reduced_matrix = SubmatrixMut::<DerefArray<_, 2>, _>::from_2d(&mut reduced);
489    lll_float(&ZZ, OwnedMatrix::<_>::identity(2, 2, Real64::RING).data(), reduced_matrix.reborrow(), 0.9, Global);
490
491    assert_lattice_isomorphic(&original, &reduced_matrix.as_const());
492    assert_eq!(4, norm_squared(&reduced_matrix.as_const().col_at(0)));
493    assert_eq!(5, norm_squared(&reduced_matrix.as_const().col_at(1)));
494}
495
496#[test]
497fn test_lll_float_3d() {
498    let ZZ = StaticRing::<i64>::RING;
499    // in this case, the shortest vector is shorter than half the second successive minimum,
500    // so LLL will find it (for delta = 0.9 > 0.75)
501    let original = [
502        DerefArray::from([72, 0, 0]),
503        DerefArray::from([0,  9, 0]),
504        DerefArray::from([8432, 7344, 16864])
505    ];
506    let _expected = [
507        [144, 72, 72],
508        [0, 279, -72],
509        [0,   0, 272]
510    ];
511
512    let mut reduced = original;
513    let mut reduced_matrix = SubmatrixMut::<DerefArray<_, 3>, _>::from_2d(&mut reduced);
514    lll_float(&ZZ, OwnedMatrix::<_>::identity(3, 3, Real64::RING).data(), reduced_matrix.reborrow(), 0.999, Global);
515
516    assert_lattice_isomorphic(&original, &reduced_matrix.as_const());
517    assert_eq!(144 * 144, norm_squared(&reduced_matrix.as_const().col_at(0)));
518    assert_eq!(72 * 72 + 279 * 279, norm_squared(&reduced_matrix.as_const().col_at(1)));
519    assert_eq!(72 * 72 * 2 + 272 * 272, norm_squared(&reduced_matrix.as_const().col_at(2)));
520}
521
522#[bench]
523fn bench_lll_float_10d(bencher: &mut Bencher) {
524    let ZZ = StaticRing::<i64>::RING;
525
526    let _expected = [
527        [  2,   0,   0,  -2,  -6,  -2,  -3,   1,  -1,  -1],
528        [  0,   0,   1,  -2,  -1,   2,  -7,  -8,   8,   1],
529        [ -1,   1,   0,   4,  -1,   1,  -1,  -5,   1, -11],
530        [  3,   1,  -2,   0,   2,   1,  -2,   1,   5, -11],
531        [ -1,   5,   3,  -1,  -1,  -2,  -3,   1,  -3,   5],
532        [  1,  -1,   3,   1,   1,   2,  -1,   0,  -6,   2],
533        [  1,   1,   0,   3,   0,  -2,   1,  -1,   4,   6],
534        [  1,   1,   2,  -1,   0,   2,   7,   1,   2,   2],
535        [  1,   0,  -4,   2,   2,   4,  -1,   3,  -3,   8],
536        [ -1,  -2,   1,   1,   0,   3,   0,   7,   5,  -2]
537    ];
538    bencher.iter(|| {
539        let original = [
540            DerefArray::from([       1,        0,        0,        0,        0,        0,        0,        0,        0,        0]),
541            DerefArray::from([       0,        1,        0,        0,        0,        0,        0,        0,        0,        0]),
542            DerefArray::from([       0,        0,        1,        0,        0,        0,        0,        0,        0,        0]),
543            DerefArray::from([       0,        0,        0,        1,        0,        0,        0,        0,        0,        0]),
544            DerefArray::from([       0,        0,        0,        0,        1,        0,        0,        0,        0,        0]),
545            DerefArray::from([       0,        0,        0,        0,        0,        1,        0,        0,        0,        0]),
546            DerefArray::from([       0,        0,        0,        0,        0,        0,        1,        0,        0,        0]),
547            DerefArray::from([       2,        2,        2,        2,        0,        0,        1,        4,        0,        0]),
548            DerefArray::from([       4,        3,        3,        3,        1,        2,        1,        0,        5,        0]),
549            DerefArray::from([ 3433883, 14315221, 24549008,  6570781, 32725387, 33674813, 27390657, 15726308, 43003827, 43364304])
550        ];
551        let mut reduced = original;
552        let mut reduced_matrix = SubmatrixMut::<DerefArray<_, 10>, _>::from_2d(&mut reduced);
553        lll_float(&ZZ, OwnedMatrix::<_>::identity(10, 10, Real64::RING).data(), reduced_matrix.reborrow(), 0.9, Global);
554
555        assert_lattice_isomorphic(&original, &reduced_matrix.as_const());
556        assert!(16 * 16 > norm_squared(&reduced_matrix.as_const().col_at(0)));
557    });
558}
559
560#[test]
561fn test_lll_exact_2d() {
562    let ZZ = StaticRing::<i64>::RING;
563    let original = [
564        DerefArray::from([5,   9]),
565        DerefArray::from([11, 20])
566    ];
567    let mut reduced = original;
568    let mut reduced_matrix = SubmatrixMut::<DerefArray<_, 2>, _>::from_2d(&mut reduced);
569    lll_exact(&ZZ, reduced_matrix.reborrow(), 0.9, Global);
570
571    assert_lattice_isomorphic(&original, &reduced_matrix.as_const());
572    assert_eq!(1, norm_squared(&reduced_matrix.as_const().col_at(0)));
573    assert_eq!(1, norm_squared(&reduced_matrix.as_const().col_at(1)));
574
575    let original = [
576        DerefArray::from([10, 8]),
577        DerefArray::from([27, 22])
578    ];
579    let mut reduced = original;
580    let mut reduced_matrix = SubmatrixMut::<DerefArray<_, 2>, _>::from_2d(&mut reduced);
581    lll_exact(&ZZ, reduced_matrix.reborrow(), 0.9, Global);
582
583    assert_lattice_isomorphic(&original, &reduced_matrix.as_const());
584    assert_eq!(4, norm_squared(&reduced_matrix.as_const().col_at(0)));
585    assert_eq!(5, norm_squared(&reduced_matrix.as_const().col_at(1)));
586}
587
588#[test]
589fn test_lll_exact_3d() {
590    let ZZ = StaticRing::<i64>::RING;
591    // in this case, the shortest vector is shorter than half the second successive minimum,
592    // so LLL will find it (for delta = 0.9 > 0.75)
593    let original = [
594        DerefArray::from([72, 0, 0]),
595        DerefArray::from([0,  9, 0]),
596        DerefArray::from([8432, 7344, 16864])
597    ];
598    let _expected = [
599        [144, 72, 72],
600        [0, 279, -72],
601        [0,   0, 272]
602    ];
603
604    let mut reduced = original;
605    let mut reduced_matrix = SubmatrixMut::<DerefArray<_, 3>, _>::from_2d(&mut reduced);
606    lll_exact(&ZZ, reduced_matrix.reborrow(), 0.999, Global);
607
608    assert_lattice_isomorphic(&original, &reduced_matrix.as_const());
609    assert_eq!(144 * 144, norm_squared(&reduced_matrix.as_const().col_at(0)));
610    assert_eq!(72 * 72 + 279 * 279, norm_squared(&reduced_matrix.as_const().col_at(1)));
611    assert_eq!(72 * 72 * 2 + 272 * 272, norm_squared(&reduced_matrix.as_const().col_at(2)));
612}
613
614#[bench]
615fn bench_lll_exact_10d(bencher: &mut Bencher) {
616    let ZZ = StaticRing::<i64>::RING;
617
618    let _expected = [
619        [  2,   0,   0,  -2,  -6,  -2,  -3,   1,  -1,  -1],
620        [  0,   0,   1,  -2,  -1,   2,  -7,  -8,   8,   1],
621        [ -1,   1,   0,   4,  -1,   1,  -1,  -5,   1, -11],
622        [  3,   1,  -2,   0,   2,   1,  -2,   1,   5, -11],
623        [ -1,   5,   3,  -1,  -1,  -2,  -3,   1,  -3,   5],
624        [  1,  -1,   3,   1,   1,   2,  -1,   0,  -6,   2],
625        [  1,   1,   0,   3,   0,  -2,   1,  -1,   4,   6],
626        [  1,   1,   2,  -1,   0,   2,   7,   1,   2,   2],
627        [  1,   0,  -4,   2,   2,   4,  -1,   3,  -3,   8],
628        [ -1,  -2,   1,   1,   0,   3,   0,   7,   5,  -2]
629    ];
630    bencher.iter(|| {
631        let original = [
632            DerefArray::from([       1,        0,        0,        0,        0,        0,        0,        0,        0,        0]),
633            DerefArray::from([       0,        1,        0,        0,        0,        0,        0,        0,        0,        0]),
634            DerefArray::from([       0,        0,        1,        0,        0,        0,        0,        0,        0,        0]),
635            DerefArray::from([       0,        0,        0,        1,        0,        0,        0,        0,        0,        0]),
636            DerefArray::from([       0,        0,        0,        0,        1,        0,        0,        0,        0,        0]),
637            DerefArray::from([       0,        0,        0,        0,        0,        1,        0,        0,        0,        0]),
638            DerefArray::from([       0,        0,        0,        0,        0,        0,        1,        0,        0,        0]),
639            DerefArray::from([       2,        2,        2,        2,        0,        0,        1,        4,        0,        0]),
640            DerefArray::from([       4,        3,        3,        3,        1,        2,        1,        0,        5,        0]),
641            DerefArray::from([ 3433883, 14315221, 24549008,  6570781, 32725387, 33674813, 27390657, 15726308, 43003827, 43364304])
642        ];
643        let mut reduced = original;
644        let mut reduced_matrix = SubmatrixMut::<DerefArray<_, 10>, _>::from_2d(&mut reduced);
645        lll_exact(&ZZ, reduced_matrix.reborrow(), 0.9, Global);
646
647        assert_lattice_isomorphic(&original, &reduced_matrix.as_const());
648        assert!(16 * 16 > norm_squared(&reduced_matrix.as_const().col_at(0)));
649    });
650}
feanor_math/algorithms/lll.rs

feanor_math/algorithms/
lll.rs
