svp/algebra/
vector.rs

1use rug::{Float, Integer};
2
3/**
4
5Generic n-vectors used to represent a basis
6
7# Examples
8
9```rust
10use svp::{Vector, nvec};
11
12// Build an integer vector
13let u = Vector{ vec: vec![0i64; 3], norm: None};
14let v = Vector{ vec: vec![0, 1, 0], norm: None};
15
16// Shorthand notation
17let u = nvec![0i64; 3];
18let v = nvec![0, 1, 0];
19```
20
21Vector multiplication corresponds to the scalar product
22
23```rust
24use svp::{Vector, nvec};
25
26let u = nvec![1, 0, 0];
27let v = nvec![0, 1, 0];
28
29assert_eq!(&u * &v, 0);
30```
31
32It can be useful to keep track of the squared norm
33
34```rust
35use svp::{Vector, nvec};
36
37// Build an integer vector
38let mut u = nvec![1, 0, 0];
39
40// Compute the squared norm (inner product)
41u.norm = Some(&u * &u);
42
43assert_eq!(u.norm.unwrap(), 1);
44```
45
46As a special case, big integers can be used for arbitrary precision
47
48```rust
49use rug::Integer;
50use svp::{Vector, nvec};
51
52// Big Integers with gmp
53let p = "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f";
54let u = nvec![Integer::from_str_radix(p, 16).unwrap(); 3];
55
56let v = nvec![Integer::new(), Integer::from(1), Integer::new()];
57assert_eq!(&v * &v, Integer::from(1));
58```
59**/
60
61#[derive(Debug, Clone)]
62pub struct Vector<T> {
63    pub vec: Vec<T>,     // n-vector
64    pub norm: Option<T>, // squared norm
65}
66
67/// `GaussReduce` with respect to v
68pub trait GaussReduce<T> {
69    fn reduce(&mut self, v: &Vector<T>) -> bool;
70}
71
72impl std::ops::Mul for &Vector<i64> {
73    /// The resulting scalar type of the inner product
74    type Output = i64;
75    /// Compute the inner product of two n-vectors
76    #[inline]
77    fn mul(self, _rhs: &Vector<i64>) -> i64 {
78        assert!(!self.vec.is_empty() && self.vec.len() == _rhs.vec.len());
79        let mut res = self.vec[0] * _rhs.vec[0];
80        for i in 1..self.vec.len() {
81            res += self.vec[i] * _rhs.vec[i];
82        }
83        res
84    }
85}
86
87impl std::ops::Mul for &Vector<f64> {
88    /// The resulting scalar type of the inner product
89    type Output = f64;
90    /// Compute the inner product of two n-vectors
91    #[inline]
92    fn mul(self, _rhs: &Vector<f64>) -> f64 {
93        assert!(!self.vec.is_empty() && self.vec.len() == _rhs.vec.len());
94        let mut res: f64 = self.vec[0] * _rhs.vec[0];
95        for i in 1..self.vec.len() {
96            res += self.vec[i] * _rhs.vec[i];
97        }
98        res
99    }
100}
101
102impl std::ops::Mul<&Vector<f64>> for &Vector<i64> {
103    /// The resulting scalar type of the inner product
104    type Output = f64;
105    /// Compute the (truncated) inner product of two n-vectors
106    #[inline]
107    fn mul(self, _rhs: &Vector<f64>) -> f64 {
108        assert!(!self.vec.is_empty() && self.vec.len() == _rhs.vec.len());
109        let mut res: f64 = self.vec[0] as f64 * _rhs.vec[0];
110        for i in 1..self.vec.len() {
111            res += self.vec[i] as f64 * _rhs.vec[i];
112        }
113        res
114    }
115}
116
117impl std::ops::Mul for &Vector<Integer> {
118    /// The resulting scalar type of the inner product
119    type Output = Integer;
120    /// Compute the inner product of two arbitrary precision n-vectors
121    #[inline]
122    fn mul(self, _rhs: &Vector<Integer>) -> Integer {
123        assert!(!self.vec.is_empty() && self.vec.len() == _rhs.vec.len());
124        let mut res: Integer = Integer::from(&self.vec[0] * &_rhs.vec[0]);
125        for i in 1..self.vec.len() {
126            res += &self.vec[i] * &_rhs.vec[i];
127        }
128        res
129    }
130}
131
132impl std::ops::Mul for &Vector<Float> {
133    /// The resulting scalar type of the inner product
134    type Output = Float;
135    /// Compute the inner product of two arbitrary precision n-vectors
136    #[inline]
137    fn mul(self, _rhs: &Vector<Float>) -> Float {
138        assert!(!self.vec.is_empty() && self.vec.len() == _rhs.vec.len());
139        let mut res: Float = Float::with_val(self.vec[0].prec(), &self.vec[0] * &_rhs.vec[0]);
140        for i in 1..self.vec.len() {
141            res += &self.vec[i] * &_rhs.vec[i];
142        }
143        res
144    }
145}
146
147impl std::ops::Mul<&Vector<Float>> for &Vector<Integer> {
148    /// The resulting scalar type of the inner product
149    type Output = Float;
150    /// Compute the inner product of two arbitrary precision n-vectors
151    #[inline]
152    fn mul(self, _rhs: &Vector<Float>) -> Float {
153        assert!(!self.vec.is_empty() && self.vec.len() == _rhs.vec.len());
154        let prec = _rhs.vec[0].prec();
155        let mut res: Float = Float::with_val(prec, &_rhs.vec[0] * &self.vec[0]);
156        for i in 1..self.vec.len() {
157            res += Float::with_val(prec, &self.vec[i] * &_rhs.vec[i]);
158        }
159        res
160    }
161}
162
163impl GaussReduce<i64> for Vector<i64> {
164    /// `GaussReduce` with respect to v
165    fn reduce(&mut self, v: &Vector<i64>) -> bool {
166        let ip = &*self * v;
167        if v.norm.unwrap() < (ip << 1).abs() {
168            let q = (ip as f64 / v.norm.unwrap() as f64).round() as i64;
169            for i in 0..self.vec.len() {
170                self.vec[i] -= q * v.vec[i];
171            }
172            self.norm = Some(&*self * &*self);
173            return true;
174        }
175        false
176    }
177}
178
179impl GaussReduce<Integer> for Vector<Integer> {
180    /// `GaussReduce` with respect to v
181    fn reduce(&mut self, v: &Vector<Integer>) -> bool {
182        let ip = &*self * v;
183        let ip2: Integer = ip.clone() * 2;
184        if v.norm.as_ref().unwrap() < &ip2.abs() {
185            let (q, _) = ip.div_rem_round(v.norm.clone().unwrap());
186            for i in 0..self.vec.len() {
187                self.vec[i] -= &q * &v.vec[i];
188            }
189            self.norm = Some(&*self * &*self);
190            return true;
191        }
192        false
193    }
194}
195
196/**
197Shorthand notation for declaring n-vectors
198
199# Examples
200
201```rust
202use svp::{Vector, nvec};
203
204// Build an integer vector
205let u = Vector{ vec: vec![0i64; 3], norm: None};
206let v = Vector{ vec: vec![0, 1, 0], norm: None};
207
208// Shorthand notation
209let u = nvec![0i64; 3];
210let v = nvec![0, 1, 0];
211```
212**/
213#[macro_export]
214macro_rules! nvec {
215    ($elem:expr; $n:expr) => (
216        Vector {
217            vec: vec![$elem; $n],
218            norm: None,
219        }
220    );
221    ($($x:expr),*) => (
222        Vector {
223            vec: <[_]>::into_vec(Box::new([$($x),*])),
224            norm: None,
225        }
226    );
227    ($($x:expr,)*) => (nvec![$($x),*])
228}
229
230#[cfg(test)]
231mod tests {
232    use super::*;
233
234    #[test]
235    fn test_prim() {
236        let mut v1 = nvec![1; 5];
237        v1.norm = Some(&v1 * &v1);
238        assert_eq!(v1.norm.unwrap(), 5);
239
240        let mut e0 = nvec![1, 0, 0];
241        let mut e1 = nvec![0, 1, 0];
242        let mut e2 = nvec![0, 0, 1];
243
244        e0.norm = Some(&e0 * &e0);
245        e1.norm = Some(&e0 * &e0);
246        e2.norm = Some(&e0 * &e0);
247        assert!(e0.norm == e1.norm && e1.norm == e2.norm);
248    }
249
250    #[test]
251    fn test_mp() {
252        let mut v1 = nvec![Integer::from(1); 5];
253        v1.norm = Some(&v1 * &v1);
254        assert_eq!(v1.norm.unwrap(), 5);
255
256        let mut e0 = nvec![Integer::from(1), Integer::new(), Integer::new()];
257        let mut e1 = nvec![Integer::new(), Integer::from(1), Integer::new()];
258        let mut e2 = nvec![Integer::new(), Integer::new(), Integer::from(1)];
259
260        e0.norm = Some(&e0 * &e0);
261        e1.norm = Some(&e1 * &e1);
262        e2.norm = Some(&e2 * &e2);
263        assert!(e0.norm == e1.norm && e1.norm == e2.norm);
264    }
265}
svp/algebra/vector.rs

svp/algebra/
vector.rs
