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