1use anyhow::Result;
  2
  3use arith::{Matrix, Ring, Rq, C, R};
  4
  5#[derive(Clone, Debug)]
  6pub struct SecretKey(Rq);
  7
  8#[derive(Clone, Debug)]
  9pub struct PublicKey(Rq, Rq);
 10
 11pub struct Encoder {
 12    n: usize,
 13    scale_factor: C<f64>, // Δ (delta)
 14    primitive: C<f64>,
 15    basis: Matrix<C<f64>>,
 16    basis_t: Matrix<C<f64>>, // transposed basis
 17}
 18
 19/// returns the mitive root of unity
 20fn primitive_root_of_unity(m: usize) -> C<f64> {
 21    let pi = C::<f64>::from(std::f64::consts::PI);
 22    ((C::<f64>::from(2f64) * pi * C::<f64>::i()) / C::<f64>::new(m as f64, 0f64)).exp()
 23}
 24
 25/// where 'w' is 'omega', the primitive root of unity
 26fn vandermonde(n: usize, w: C<f64>) -> Matrix<C<f64>> {
 27    let mut v: Vec<Vec<C<f64>>> = vec![];
 28    for i in 0..n {
 29        let root = w.pow(2 * i as u32 + 1);
 30        let mut row: Vec<C<f64>> = vec![];
 31        for j in 0..n {
 32            row.push(root.pow(j as u32));
 33        }
 34        v.push(row);
 35    }
 36    Matrix::<C<f64>>(v)
 37}
 38impl Encoder {
 39    pub fn new(n: usize, scale_factor: C<f64>) -> Self {
 40        let primitive: C<f64> = primitive_root_of_unity(2 * n);
 41        let basis = vandermonde(n, primitive);
 42        let basis_t = basis.transpose();
 43
 44        Self {
 45            n,
 46            scale_factor,
 47            primitive,
 48            basis,
 49            basis_t,
 50        }
 51    }
 52
 53    /// encode as described in the CKKS paper.
 54    /// from $\mathbb{C}^{N/2} \longrightarrow \mathbb{Z_q}[X]/(X^N +1) = R$
 55    // TODO use alg.1 from 2018-1043,
 56    //      or as in 2018-1073: $f(x) = 1N (U^T.conj() m + U^T m.conj())$
 57    pub fn encode(&self, z: &[C<f64>]) -> Result<R> {
 58        // $pi^{-1}: \mathbb{C}^{N/2} \longrightarrow \mathbb{H}$
 59        let expanded = self.pi_inv(z);
 60
 61        // scale the values
 62        let scaled: Vec<C<f64>> = expanded.iter().map(|e| *e * self.scale_factor).collect();
 63
 64        // but $\mathbb{H} \neq \sigma(R)$, since $\sigma(R) \subseteq \mathbb{H}$, so we need to
 65        // discretize $\pi^{-1}(z)$ into an element of $\sigma(R)$.
 66
 67        // discretize \pi^-1(z_projected) to \sigma(R)
 68        // project 'scaled' into \sigma(R):
 69        // get the orthogonal basis (note: that would be doing Gram-Schmidt, which is not this, but
 70        // we're fine since the basis=Vandermonde matrix which is orthogonal, so we project z to it):
 71        // $z = \sum z_i * b_i, with z_i = <z,b_i>/||b_i||^2$
 72        let z_projected = self
 73            .basis_t
 74            .0
 75            .iter()
 76            .map(|b_i| {
 77                // TODO: the b_j.conj() can be precomputed at initialization (of the basis)
 78                let num: C<f64> = scaled
 79                    .iter()
 80                    .zip(b_i.iter())
 81                    .map(|(z_j, b_j)| *z_j * b_j.conj())
 82                    .sum::<C<f64>>();
 83                let den: C<f64> = b_i.iter().map(|b_j| *b_j * b_j.conj()).sum::<C<f64>>();
 84                let mut z_i = num / den;
 85                z_i.im = 0.0; // get only the real component
 86                z_i
 87            })
 88            .collect::<Vec<C<f64>>>();
 89
 90        // V * z_projected (V: Vandermonde matrix)
 91        let discretized = self.basis.mul_vec(&z_projected)?;
 92
 93        // sigma_inv
 94        let r = self.sigma_inv(&discretized)?;
 95
 96        // TMP: naive round, maybe do gaussian
 97        let coeffs = r.iter().map(|e| e.re.round() as i64).collect::<Vec<i64>>();
 98        Ok(R::from_vec(self.n, coeffs))
 99    }
100
101    pub fn decode(&self, p: &R) -> Result<Vec<C<f64>>> {
102        let p: Vec<C<f64>> = p
103            .coeffs()
104            .iter()
105            .map(|&e| C::<f64>::new(e as f64, 0_f64)) // TODO review u64 to f64 conversion overflow
106            .collect();
107        let in_sigma = self.sigma(&p)?;
108
109        let deescalated: Vec<C<f64>> = in_sigma.iter().map(|e| *e / self.scale_factor).collect();
110        Ok(self.pi(&deescalated))
111    }
112
113    /// pi: \mathbb{H} \longrightarrow \mathbb{C}^{N/2}
114    fn pi(&self, z: &[C<f64>]) -> Vec<C<f64>> {
115        z[..self.n / 2].to_vec()
116    }
117    /// pi^{-1}: \mathbb{C}^{N/2} \longrightarrow \mathbb{H}
118    fn pi_inv(&self, z: &[C<f64>]) -> Vec<C<f64>> {
119        z.iter()
120            .cloned()
121            .chain(z.iter().rev().map(|z_i| z_i.conj()))
122            .collect()
123    }
124
125    fn sigma(&self, p: &[C<f64>]) -> Result<Vec<C<f64>>> {
126        // the roots of unity are already calculated in the 2nd row of the transpose of the
127        // Vandermonde matrix used as the basis (ie. the 2nd column of the Vandermonde matrix).
128        // let roots = &self.basis_t[1];
129        // // Approach 1: evaluate p at the roots of unity
130        // let mut z = vec![];
131        // for root_i in roots.iter() {
132        //     z.push(eval(p, root_i));
133        // }
134
135        // Approach 2: Vandermonde * p
136        let z: Vec<C<f64>> = self.basis.mul_vec(&p.to_vec())?;
137
138        // TODO check using NTT-ish (2018-1043) for the encode/decode
139
140        Ok(z)
141    }
142
143    fn sigma_inv(&self, z: &Vec<C<f64>>) -> Result<Vec<C<f64>>> {
144        // $\alpha = A^{-1} * z$
145        let a = self.basis.solve(z)?;
146        Ok(a.to_vec())
147    }
148}
149
150#[cfg(test)]
151mod tests {
152    use super::*;
153    use rand::Rng;
154
155    #[test]
156    fn test_encode_decode() -> Result<()> {
157        const Q: u64 = 1024;
158        const N: usize = 32;
159        let n: usize = 32;
160
161        let T = 128; // WIP
162        let mut rng = rand::thread_rng();
163
164        for _ in 0..100 {
165            let z: Vec<C<f64>> = std::iter::repeat_with(|| {
166                C::<f64>::new(rng.gen_range(0..T) as f64, rng.gen_range(0..T) as f64)
167            })
168            .take(N / 2)
169            .collect();
170
171            let delta = C::<f64>::new(64.0, 0.0); // delta = scaling factor
172            let encoder = Encoder::new(n, delta);
173
174            let m: R = encoder.encode(&z)?; // polynomial (encoded vec) \in R
175
176            let z_decoded = encoder.decode(&m)?;
177
178            // round it to compare it to the initial value
179            let rounded_z_decoded: Vec<C<f64>> = z_decoded
180                .iter()
181                .map(|c| C::<f64>::new(c.re.round(), c.im.round()))
182                .collect();
183            assert_eq!(rounded_z_decoded, z);
184        }
185
186        Ok(())
187    }
188}