adlo/
lib.rs

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

adlo/
lib.rs
