adlo/
lib.rs

1//! ## Example
2//!
3//! ```rust
4//! use adlo::*;
5//! use nalgebra::{DVector, DMatrix};
6//!
7//! fn basis_2d_test() {
8//!        let mut basis_2d = vec![
9//!            AdloVector::new(DVector::from_vec(vec![5, 3])),
10//!            AdloVector::new(DVector::from_vec(vec![2, 2]))
11//!         ];
12//!        adaptive_lll(&mut basis_2d);
13//!        let expected = AdloVector::new(DVector::from_vec(vec![1, -1]));
14//!        assert_eq!(expected, basis_2d[basis_2d.len()-1]);
15//!    }
16//!```
17
18use nalgebra::{DMatrix, DVector};
19
20// Used for calculating adaptive lovaz factor;
21pub const PSI: f64 = 0.618;
22const S_PSI: f64 = 1.0 - PSI;
23
24/// The Vector struct represents a distinction in n-dimensional space.  
25#[derive(Debug, Clone, Default, PartialEq, PartialOrd)]
26pub struct AdloVector {
27    pub coords: DVector<i64>,
28    pub f64_coords: DVector<f64>,
29    pub norm_sq: i128,
30}
31
32impl AdloVector {
33    /// Create a new Adaptive LLL Vector
34    pub fn new(coords: DVector<i64>) -> Self {
35        let f64_coords = coords.clone().cast::<f64>();
36        let norm_sq = Self::calculate_norm_sq(&coords);
37        AdloVector {
38            coords,
39            f64_coords,
40            norm_sq,
41        }
42    }
43    fn calculate_norm_sq(coords: &DVector<i64>) -> i128 {
44        let mut sum: i128 = 0;
45        for &coord in coords.iter() {
46            let coord_i128 = coord as i128;
47            sum += coord_i128 * coord_i128;
48        }
49        sum
50    }
51    fn norm_squared(&self) -> i128 {
52        self.norm_sq
53    }
54    pub fn norm(&self) -> f64 {
55        (self.norm_squared() as f64).sqrt()
56    }
57    pub fn set_f64_coords(&mut self) {
58        self.f64_coords = self.coords.clone().cast::<f64>();
59    }
60}
61
62fn calculate_local_density(basis: &[AdloVector], k: usize) -> f64 {
63    let n = basis.len();
64    let mut count = 0;
65    // TODO: Experiment with search radius
66    let radius = basis[k].norm() * 1.5;
67
68    for i in 0..n {
69        if i != k {
70            let dist = basis[k].norm() - basis[i].norm();
71            if dist <= radius {
72                count += 1;
73            }
74        }
75    }
76    // A simple density measure: number of nearby vectors
77    (count as f64) / (n as f64 - 1.0)
78}
79
80fn calculate_lovasz_factor(density: f64) -> f64 {
81    // Higher density -> factor closer to 1
82    // TODO: Experiment with different functions to map density to the factor
83    PSI + S_PSI * density.min(1.0) // Linear interpolation (example)
84}
85
86/// The LLL (Lenstra–Lenstra–Lovász) algorithm is a lattice basis reduction algorithm that finds a basis with short,
87///
88/// nearly orthogonal vectors for any given lattice. It is particularly useful for solving problems in number theory and cryptography,
89///
90/// such as finding integer relations and breaking certain types of encryption schemes.
91///
92/// The Gram-Schmidt process is a method used in linear algebra to transform a set of linearly independent vectors into
93///
94/// an orthogonal set of vectors that span the same space. This process can further be used to create an orthonormal set
95///
96/// of vectors by normalizing each vector in the orthogonal set to have a unit length
97///
98/// The lovasz_factor is calculated via `PSI + 1 - PSI * density.min(1.0)` where density is the
99///
100/// measure of nearby vectors.
101pub fn adaptive_lll(basis: &mut [AdloVector]) {
102    let n = basis.len();
103    let mut b_star: Vec<AdloVector> = Vec::with_capacity(n);
104    let mut mu: DMatrix<f64> = DMatrix::from_element(n, n, 0.0);
105
106    // Gram-Schmidt
107    for i in 0..n {
108        b_star.push(basis[i].clone());
109        for j in 0..i {
110            mu[(i, j)] = b_star[i].f64_coords.dot(&b_star[j].f64_coords) / b_star[j].norm().powi(2);
111            // Calculate new_coords_f64 using i128 for intermediate values:
112            let mut new_coords_f64 = DVector::<f64>::zeros(n);
113            for dim in 0..n {
114                let mu_val = mu[(i, j)];
115                let b_star_j_coord = b_star[j].f64_coords[dim];
116                let intermediate_val = (mu_val * b_star_j_coord) as f64; // Cast to f64 to prevent overflow
117                new_coords_f64[dim] = intermediate_val;
118            }
119            b_star[i].f64_coords -= new_coords_f64.clone();
120            b_star[i].coords = new_coords_f64.map(|x| x.round() as i64);
121        }
122    }
123
124    let mut k = 1;
125    while k < n {
126        // Size reduction
127        for i in (0..k).rev() {
128            let r_f64 = mu[(k, i)];
129            let r_i128 = r_f64.round() as i128; // Use i128 here
130            let mut change_basis = DVector::<i64>::zeros(n);
131
132            for dim in 0..n {
133                change_basis[dim] = (r_i128 * basis[i].coords[dim] as i128).try_into().unwrap(); // Convert i128 back to i64
134            }
135
136            basis[k].coords -= change_basis;
137            basis[k].set_f64_coords();
138
139            for l in 0..=i {
140                mu[(k, l)] -= r_f64 * mu[(i, l)];
141            }
142        }
143        // Adaptive Lovász condition
144        let local_density = calculate_local_density(basis, k);
145        let lovasz_factor = calculate_lovasz_factor(local_density);
146        // Lovász condition
147        if b_star[k].norm().powi(2) + (mu[(k, k - 1)].round() as f64 * b_star[k - 1].norm().powi(2))
148            < lovasz_factor * b_star[k - 1].norm().powi(2)
149        {
150            basis.swap(k - 1, k);
151            // Recompute mu and b_star after swap
152            b_star.clear();
153            mu = DMatrix::zeros(n, n);
154            for i in 0..n {
155                b_star.push(basis[i].clone());
156                for j in 0..i {
157                    mu[(i, j)] =
158                        b_star[i].f64_coords.dot(&b_star[j].f64_coords) / b_star[j].norm().powi(2);
159                    let new_coords_f64 = mu[(i, j)] * &b_star[j].f64_coords;
160                    b_star[i].f64_coords -= new_coords_f64.clone();
161                    b_star[i].coords = new_coords_f64.map(|x| x.round() as i64);
162                }
163            }
164            k = 1;
165        } else {
166            k += 1;
167        }
168    }
169}
170
171pub fn create_rotation_matrix(n: usize, i: usize, j: usize, theta: f64) -> DMatrix<f64> {
172    let mut matrix = DMatrix::identity(n, n);
173    let cos_theta = theta.cos();
174    let sin_theta = theta.sin();
175
176    matrix[(i, i)] = cos_theta;
177    matrix[(j, j)] = cos_theta;
178    matrix[(i, j)] = -sin_theta;
179    matrix[(j, i)] = sin_theta;
180
181    matrix
182}
183
184pub fn rotate_vector(v: &DVector<f64>, i: usize, j: usize, theta: f64) -> DVector<f64> {
185    let n = v.len();
186    let rotation_matrix = create_rotation_matrix(n, i, j, theta);
187    rotation_matrix * v
188}
189
190/// Multi-frame search rotates the basis by 45 degrees and
191///
192/// calculates best basis via `adaptive_lll` on each frame.
193pub fn multi_frame_search(basis: &mut Vec<AdloVector>) {
194    let mut best_basis = basis.to_vec();
195    let mut best_norm = best_basis[0].norm();
196    let n = basis.len();
197    // Generate multiple frames
198    for angle in (0..360).step_by(45) {
199        let angle_rad = (angle as f64).to_radians();
200        // Iterate over all pairs of dimensions
201        for i in 0..n {
202            for j in i + 1..n {
203                // Avoid rotating a dimension with itself and avoid duplicates
204                let mut rotated_basis = basis.to_vec();
205                for vec in &mut rotated_basis {
206                    // Rotate in the i-j plane
207                    let rotated_f64_coords = rotate_vector(&vec.f64_coords, i, j, angle_rad);
208                    let rotated_i64_coords = rotated_f64_coords.map(|x| x.round() as i64);
209                    vec.f64_coords = rotated_f64_coords;
210                    vec.coords = rotated_i64_coords;
211                }
212                adaptive_lll(&mut rotated_basis);
213                if rotated_basis[0].norm() < best_norm {
214                    best_norm = rotated_basis[0].norm();
215                    best_basis = rotated_basis;
216                }
217            }
218        }
219    }
220    *basis = best_basis;
221}
222
223// Tests
224//-------------------------------------------------------------------------------
225#[cfg(test)]
226mod tests {
227
228    use super::*;
229
230    #[test]
231    fn basis_2d_test() {
232        let mut basis_2d = vec![
233            AdloVector::new(DVector::from_vec(vec![5, 3])),
234            AdloVector::new(DVector::from_vec(vec![2, 2])),
235        ];
236        multi_frame_search(&mut basis_2d);
237        let mut expected = AdloVector::new(DVector::from_vec(vec![1, -1]));
238        expected.norm_sq = 34;
239        assert_eq!(expected, basis_2d[basis_2d.len() - 1]);
240    }
241
242    #[test]
243    fn basis_5d_test() {
244        let mut basis_5d = vec![
245            AdloVector::new(DVector::from_vec(vec![1, 1, 1])),
246            AdloVector::new(DVector::from_vec(vec![-1, 0, 2])),
247            AdloVector::new(DVector::from_vec(vec![3, 5, 6])),
248        ];
249        multi_frame_search(&mut basis_5d);
250        let mut expected = AdloVector::new(DVector::from_vec(vec![-1, 0, 2]));
251        expected.norm_sq = 5;
252        assert_eq!(expected, basis_5d[basis_5d.len() - 2]);
253    }
254}
adlo/lib.rs

adlo/
lib.rs
