poulpy/core/src/blind_rotation/ccgi.rs
core/src/blind_rotation/ccgi.rsRaw
1use backend::{
2 FFT64, MatZnxDftOps, MatZnxDftScratch, Module, ScalarZnx, ScalarZnxAlloc, ScalarZnxDft, ScalarZnxDftAlloc, ScalarZnxDftOps,
3 Scratch, VecZnxAlloc, VecZnxBigAlloc, VecZnxBigOps, VecZnxBigScratch, VecZnxDftAlloc, VecZnxDftOps, VecZnxOps, ZnxView,
4 ZnxViewMut, ZnxZero,
5};
6use itertools::izip;
7
8use crate::{
9 GLWECiphertext, GLWECiphertextToMut, GLWEOps, Infos, LWECiphertext, ScratchCore,
10 blind_rotation::{key::BlindRotationKeyCGGI, lut::LookUpTable},
11 dist::Distribution,
12 lwe::ciphertext::LWECiphertextToRef,
13};
14
15pub fn cggi_blind_rotate_scratch_space(
16 module: &Module<FFT64>,
17 block_size: usize,
18 extension_factor: usize,
19 basek: usize,
20 k_res: usize,
21 k_brk: usize,
22 rows: usize,
23 rank: usize,
24) -> usize {
25 let brk_size: usize = k_brk.div_ceil(basek);
26
27 if block_size > 1 {
28 let cols: usize = rank + 1;
29 let acc_dft: usize = module.bytes_of_vec_znx_dft(cols, rows) * extension_factor;
30 let acc_big: usize = module.bytes_of_vec_znx_big(1, brk_size);
31 let vmp_res: usize = module.bytes_of_vec_znx_dft(cols, brk_size) * extension_factor;
32 let acc_dft_add: usize = vmp_res;
33 let xai_plus_y: usize = module.bytes_of_scalar_znx(1);
34 let xai_plus_y_dft: usize = module.bytes_of_scalar_znx_dft(1);
35 let vmp: usize = module.vmp_apply_tmp_bytes(brk_size, rows, rows, 2, 2, brk_size); // GGSW product: (1 x 2) x (2 x 2)
36
37 let acc: usize;
38 if extension_factor > 1 {
39 acc = module.bytes_of_vec_znx(cols, k_res.div_ceil(basek)) * extension_factor;
40 } else {
41 acc = 0;
42 }
43
44 return acc
45 + acc_dft
46 + acc_dft_add
47 + vmp_res
48 + xai_plus_y
49 + xai_plus_y_dft
50 + (vmp | (acc_big + (module.vec_znx_big_normalize_tmp_bytes() | module.vec_znx_idft_tmp_bytes())));
51 } else {
52 2 * GLWECiphertext::bytes_of(module, basek, k_res, rank)
53 + GLWECiphertext::external_product_scratch_space(module, basek, k_res, k_res, k_brk, 1, rank)
54 }
55}
56
57pub fn cggi_blind_rotate<DataRes, DataIn>(
58 module: &Module<FFT64>,
59 res: &mut GLWECiphertext<DataRes>,
60 lwe: &LWECiphertext<DataIn>,
61 lut: &LookUpTable,
62 brk: &BlindRotationKeyCGGI<FFT64>,
63 scratch: &mut Scratch,
64) where
65 DataRes: AsRef<[u8]> + AsMut<[u8]>,
66 DataIn: AsRef<[u8]>,
67{
68 match brk.dist {
69 Distribution::BinaryBlock(_) | Distribution::BinaryFixed(_) | Distribution::BinaryProb(_) | Distribution::ZERO => {
70 if lut.extension_factor() > 1 {
71 cggi_blind_rotate_block_binary_extended(module, res, lwe, lut, brk, scratch);
72 } else if brk.block_size() > 1 {
73 cggi_blind_rotate_block_binary(module, res, lwe, lut, brk, scratch);
74 } else {
75 cggi_blind_rotate_binary_standard(module, res, lwe, lut, brk, scratch);
76 }
77 }
78 // TODO: ternary distribution ?
79 _ => panic!(
80 "invalid BlindRotationKeyCGGI distribution: must be BinaryBlock, BinaryFixed or BinaryProb (or ZERO for debugging)"
81 ),
82 }
83}
84
85// TODO: ENSURE DOMAIN EXTENSION AS
86pub(crate) fn cggi_blind_rotate_block_binary_extended<DataRes, DataIn>(
87 module: &Module<FFT64>,
88 res: &mut GLWECiphertext<DataRes>,
89 lwe: &LWECiphertext<DataIn>,
90 lut: &LookUpTable,
91 brk: &BlindRotationKeyCGGI<FFT64>,
92 scratch: &mut Scratch,
93) where
94 DataRes: AsRef<[u8]> + AsMut<[u8]>,
95 DataIn: AsRef<[u8]>,
96{
97 let extension_factor: usize = lut.extension_factor();
98 let basek: usize = res.basek();
99 let rows: usize = brk.rows();
100 let cols: usize = res.rank() + 1;
101
102 let (mut acc, scratch1) = scratch.tmp_slice_vec_znx(extension_factor, module, cols, res.size());
103 let (mut acc_dft, scratch2) = scratch1.tmp_slice_vec_znx_dft(extension_factor, module, cols, rows);
104 let (mut vmp_res, scratch3) = scratch2.tmp_slice_vec_znx_dft(extension_factor, module, cols, brk.size());
105 let (mut acc_add_dft, scratch4) = scratch3.tmp_slice_vec_znx_dft(extension_factor, module, cols, brk.size());
106 let (mut xai_plus_y, scratch5) = scratch4.tmp_scalar_znx(module, 1);
107 let (mut xai_plus_y_dft, scratch6) = scratch5.tmp_scalar_znx_dft(module, 1);
108
109 (0..extension_factor).for_each(|i| {
110 acc[i].zero();
111 });
112
113 let mut lwe_2n: Vec<i64> = vec![0i64; lwe.n() + 1]; // TODO: from scratch space
114 let lwe_ref: LWECiphertext<&[u8]> = lwe.to_ref();
115
116 let two_n_ext: usize = 2 * lut.domain_size();
117
118 negate_and_mod_switch_2n(two_n_ext, &mut lwe_2n, &lwe_ref);
119
120 let a: &[i64] = &lwe_2n[1..];
121 let b_pos: usize = ((lwe_2n[0] + two_n_ext as i64) % two_n_ext as i64) as usize;
122
123 let b_hi: usize = b_pos / extension_factor;
124 let b_lo: usize = b_pos % extension_factor;
125
126 for (i, j) in (0..b_lo).zip(extension_factor - b_lo..extension_factor) {
127 module.vec_znx_rotate(b_hi as i64 + 1, &mut acc[i], 0, &lut.data[j], 0);
128 }
129 for (i, j) in (b_lo..extension_factor).zip(0..extension_factor - b_lo) {
130 module.vec_znx_rotate(b_hi as i64, &mut acc[i], 0, &lut.data[j], 0);
131 }
132
133 let block_size: usize = brk.block_size();
134
135 izip!(
136 a.chunks_exact(block_size),
137 brk.data.chunks_exact(block_size)
138 )
139 .for_each(|(ai, ski)| {
140 (0..extension_factor).for_each(|i| {
141 (0..cols).for_each(|j| {
142 module.vec_znx_dft(1, 0, &mut acc_dft[i], j, &acc[i], j);
143 });
144 acc_add_dft[i].zero();
145 });
146
147 // TODO: first & last iterations can be optimized
148 izip!(ai.iter(), ski.iter()).for_each(|(aii, skii)| {
149 let ai_pos: usize = ((aii + two_n_ext as i64) % two_n_ext as i64) as usize;
150 let ai_hi: usize = ai_pos / extension_factor;
151 let ai_lo: usize = ai_pos % extension_factor;
152
153 // vmp_res = DFT(acc) * BRK[i]
154 (0..extension_factor).for_each(|i| {
155 module.vmp_apply(&mut vmp_res[i], &acc_dft[i], &skii.data, scratch6);
156 });
157
158 // Trivial case: no rotation between polynomials, we can directly multiply with (X^{-ai} - 1)
159 if ai_lo == 0 {
160 // DFT X^{-ai}
161 set_xai_plus_y(module, ai_hi, -1, &mut xai_plus_y_dft, &mut xai_plus_y);
162
163 // Sets acc_add_dft[i] = (acc[i] * sk) * (X^{-ai} - 1)
164 (0..extension_factor).for_each(|j| {
165 (0..cols).for_each(|i| {
166 module.svp_apply_inplace(&mut vmp_res[j], i, &xai_plus_y_dft, 0);
167 module.vec_znx_dft_add_inplace(&mut acc_add_dft[j], i, &vmp_res[j], i);
168 });
169 });
170 // Non trivial case: rotation between polynomials
171 // In this case we can't directly multiply with (X^{-ai} - 1) because of the
172 // ring homomorphism R^{N} -> prod R^{N/extension_factor}, so we split the
173 // computation in two steps: acc_add_dft = (acc * sk) * (-1) + (acc * sk) * X^{-ai}
174 } else {
175 // Sets acc_add_dft[i] = acc[i] * sk
176 (0..extension_factor).for_each(|i| {
177 (0..cols).for_each(|k| {
178 module.vec_znx_dft_sub_ab_inplace(&mut acc_add_dft[i], k, &vmp_res[i], k);
179 })
180 });
181
182 // DFT X^{-ai}
183 set_xai_plus_y(module, ai_hi + 1, 0, &mut xai_plus_y_dft, &mut xai_plus_y);
184
185 // Sets acc_add_dft[0..ai_lo] += (acc[extension_factor - ai_lo..extension_factor] * sk) * X^{-ai+1}
186 for (i, j) in (0..ai_lo).zip(extension_factor - ai_lo..extension_factor) {
187 (0..cols).for_each(|k| {
188 module.svp_apply_inplace(&mut vmp_res[j], k, &xai_plus_y_dft, 0);
189 module.vec_znx_dft_add_inplace(&mut acc_add_dft[i], k, &vmp_res[j], k);
190 });
191 }
192
193 // DFT X^{-ai}
194 set_xai_plus_y(module, ai_hi, 0, &mut xai_plus_y_dft, &mut xai_plus_y);
195
196 // Sets acc_add_dft[ai_lo..extension_factor] += (acc[0..extension_factor - ai_lo] * sk) * X^{-ai}
197 for (i, j) in (ai_lo..extension_factor).zip(0..extension_factor - ai_lo) {
198 (0..cols).for_each(|k| {
199 module.svp_apply_inplace(&mut vmp_res[j], k, &xai_plus_y_dft, 0);
200 module.vec_znx_dft_add_inplace(&mut acc_add_dft[i], k, &vmp_res[j], k);
201 });
202 }
203 }
204 });
205
206 {
207 let (mut acc_add_big, scratch7) = scratch6.tmp_vec_znx_big(module, 1, brk.size());
208
209 (0..extension_factor).for_each(|j| {
210 (0..cols).for_each(|i| {
211 module.vec_znx_idft(&mut acc_add_big, 0, &acc_add_dft[j], i, scratch7);
212 module.vec_znx_big_add_small_inplace(&mut acc_add_big, 0, &acc[j], i);
213 module.vec_znx_big_normalize(basek, &mut acc[j], i, &acc_add_big, 0, scratch7);
214 });
215 });
216 }
217 });
218
219 (0..cols).for_each(|i| {
220 module.vec_znx_copy(&mut res.data, i, &acc[0], i);
221 });
222}
223
224fn set_xai_plus_y(
225 module: &Module<FFT64>,
226 ai: usize,
227 y: i64,
228 res: &mut ScalarZnxDft<&mut [u8], FFT64>,
229 buf: &mut ScalarZnx<&mut [u8]>,
230) {
231 let n: usize = module.n();
232
233 {
234 let raw: &mut [i64] = buf.at_mut(0, 0);
235 if ai < n {
236 raw[ai] = 1;
237 } else {
238 raw[(ai - n) & (n - 1)] = -1;
239 }
240 raw[0] += y;
241 }
242
243 module.svp_prepare(res, 0, buf, 0);
244
245 {
246 let raw: &mut [i64] = buf.at_mut(0, 0);
247
248 if ai < n {
249 raw[ai] = 0;
250 } else {
251 raw[(ai - n) & (n - 1)] = 0;
252 }
253 raw[0] = 0;
254 }
255}
256
257pub(crate) fn cggi_blind_rotate_block_binary<DataRes, DataIn>(
258 module: &Module<FFT64>,
259 res: &mut GLWECiphertext<DataRes>,
260 lwe: &LWECiphertext<DataIn>,
261 lut: &LookUpTable,
262 brk: &BlindRotationKeyCGGI<FFT64>,
263 scratch: &mut Scratch,
264) where
265 DataRes: AsRef<[u8]> + AsMut<[u8]>,
266 DataIn: AsRef<[u8]>,
267{
268 let mut lwe_2n: Vec<i64> = vec![0i64; lwe.n() + 1]; // TODO: from scratch space
269 let mut out_mut: GLWECiphertext<&mut [u8]> = res.to_mut();
270 let lwe_ref: LWECiphertext<&[u8]> = lwe.to_ref();
271 let two_n: usize = module.n() << 1;
272 let basek: usize = brk.basek();
273 let rows = brk.rows();
274
275 let cols: usize = out_mut.rank() + 1;
276
277 negate_and_mod_switch_2n(2 * lut.domain_size(), &mut lwe_2n, &lwe_ref);
278
279 let a: &[i64] = &lwe_2n[1..];
280 let b: i64 = lwe_2n[0];
281
282 out_mut.data.zero();
283
284 // Initialize out to X^{b} * LUT(X)
285 module.vec_znx_rotate(b, &mut out_mut.data, 0, &lut.data[0], 0);
286
287 let block_size: usize = brk.block_size();
288
289 // ACC + [sum DFT(X^ai -1) * (DFT(ACC) x BRKi)]
290
291 let (mut acc_dft, scratch1) = scratch.tmp_vec_znx_dft(module, cols, rows);
292 let (mut vmp_res, scratch2) = scratch1.tmp_vec_znx_dft(module, cols, brk.size());
293 let (mut acc_add_dft, scratch3) = scratch2.tmp_vec_znx_dft(module, cols, brk.size());
294 let (mut xai_plus_y, scratch4) = scratch3.tmp_scalar_znx(module, 1);
295 let (mut xai_plus_y_dft, scratch5) = scratch4.tmp_scalar_znx_dft(module, 1);
296
297 izip!(
298 a.chunks_exact(block_size),
299 brk.data.chunks_exact(block_size)
300 )
301 .for_each(|(ai, ski)| {
302 (0..cols).for_each(|j| {
303 module.vec_znx_dft(1, 0, &mut acc_dft, j, &out_mut.data, j);
304 });
305
306 acc_add_dft.zero();
307
308 izip!(ai.iter(), ski.iter()).for_each(|(aii, skii)| {
309 let ai_pos: usize = ((aii + two_n as i64) % two_n as i64) as usize;
310
311 // vmp_res = DFT(acc) * BRK[i]
312 module.vmp_apply(&mut vmp_res, &acc_dft, &skii.data, scratch5);
313
314 // DFT(X^ai -1)
315 set_xai_plus_y(module, ai_pos, -1, &mut xai_plus_y_dft, &mut xai_plus_y);
316
317 // DFT(X^ai -1) * (DFT(acc) * BRK[i])
318 (0..cols).for_each(|i| {
319 module.svp_apply_inplace(&mut vmp_res, i, &xai_plus_y_dft, 0);
320 module.vec_znx_dft_add_inplace(&mut acc_add_dft, i, &vmp_res, i);
321 });
322 });
323
324 (0..cols).for_each(|i| {
325 module.vec_znx_dft_add_inplace(&mut acc_dft, i, &acc_add_dft, i);
326 });
327
328 {
329 let (mut acc_add_big, scratch6) = scratch5.tmp_vec_znx_big(module, 1, brk.size());
330
331 (0..cols).for_each(|i| {
332 module.vec_znx_idft(&mut acc_add_big, 0, &acc_add_dft, i, scratch6);
333 module.vec_znx_big_add_small_inplace(&mut acc_add_big, 0, &out_mut.data, i);
334 module.vec_znx_big_normalize(basek, &mut out_mut.data, i, &acc_add_big, 0, scratch6);
335 });
336 }
337 });
338}
339
340pub(crate) fn cggi_blind_rotate_binary_standard<DataRes, DataIn>(
341 module: &Module<FFT64>,
342 res: &mut GLWECiphertext<DataRes>,
343 lwe: &LWECiphertext<DataIn>,
344 lut: &LookUpTable,
345 brk: &BlindRotationKeyCGGI<FFT64>,
346 scratch: &mut Scratch,
347) where
348 DataRes: AsRef<[u8]> + AsMut<[u8]>,
349 DataIn: AsRef<[u8]>,
350{
351 #[cfg(debug_assertions)]
352 {
353 assert_eq!(
354 res.n(),
355 module.n(),
356 "res.n(): {} != brk.n(): {}",
357 res.n(),
358 module.n()
359 );
360 assert_eq!(
361 lut.domain_size(),
362 module.n(),
363 "lut.n(): {} != brk.n(): {}",
364 lut.domain_size(),
365 module.n()
366 );
367 assert_eq!(
368 brk.n(),
369 module.n(),
370 "brk.n(): {} != brk.n(): {}",
371 brk.n(),
372 module.n()
373 );
374 assert_eq!(
375 res.rank(),
376 brk.rank(),
377 "res.rank(): {} != brk.rank(): {}",
378 res.rank(),
379 brk.rank()
380 );
381 assert_eq!(
382 lwe.n(),
383 brk.data.len(),
384 "lwe.n(): {} != brk.data.len(): {}",
385 lwe.n(),
386 brk.data.len()
387 );
388 }
389
390 let mut lwe_2n: Vec<i64> = vec![0i64; lwe.n() + 1]; // TODO: from scratch space
391 let mut out_mut: GLWECiphertext<&mut [u8]> = res.to_mut();
392 let lwe_ref: LWECiphertext<&[u8]> = lwe.to_ref();
393 let basek: usize = brk.basek();
394
395 negate_and_mod_switch_2n(2 * lut.domain_size(), &mut lwe_2n, &lwe_ref);
396
397 let a: &[i64] = &lwe_2n[1..];
398 let b: i64 = lwe_2n[0];
399
400 out_mut.data.zero();
401
402 // Initialize out to X^{b} * LUT(X)
403 module.vec_znx_rotate(b, &mut out_mut.data, 0, &lut.data[0], 0);
404
405 // ACC + [sum DFT(X^ai -1) * (DFT(ACC) x BRKi)]
406 let (mut acc_tmp, scratch1) = scratch.tmp_glwe_ct(module, basek, out_mut.k(), out_mut.rank());
407 let (mut acc_tmp_rot, scratch2) = scratch1.tmp_glwe_ct(module, basek, out_mut.k(), out_mut.rank());
408
409 // TODO: see if faster by skipping normalization in external product and keeping acc in big coeffs
410 // TODO: first iteration can be optimized to be a gglwe product
411 izip!(a.iter(), brk.data.iter()).for_each(|(ai, ski)| {
412 // acc_tmp = sk[i] * acc
413 acc_tmp.external_product(module, &out_mut, ski, scratch2);
414
415 // acc_tmp = (sk[i] * acc) * X^{ai}
416 acc_tmp_rot.rotate(module, *ai, &acc_tmp);
417
418 // acc = acc + (sk[i] * acc) * X^{ai}
419 out_mut.add_inplace(module, &acc_tmp_rot);
420
421 // acc = acc + (sk[i] * acc) * X^{ai} - (sk[i] * acc) = acc + (sk[i] * acc) * (X^{ai} - 1)
422 out_mut.sub_inplace_ab(module, &acc_tmp);
423 });
424
425 // We can normalize only at the end because we add normalized values in [-2^{basek-1}, 2^{basek-1}]
426 // on top of each others, thus ~ 2^{63-basek} additions are supported before overflow.
427 out_mut.normalize_inplace(module, scratch2);
428}
429
430pub(crate) fn negate_and_mod_switch_2n(n: usize, res: &mut [i64], lwe: &LWECiphertext<&[u8]>) {
431 let basek: usize = lwe.basek();
432
433 let log2n: usize = usize::BITS as usize - (n - 1).leading_zeros() as usize + 1;
434
435 res.copy_from_slice(&lwe.data.at(0, 0));
436 res.iter_mut().for_each(|x| *x = -*x);
437
438 if basek > log2n {
439 let diff: usize = basek - log2n;
440 res.iter_mut().for_each(|x| {
441 *x = div_ceil_signed_by_pow2(x, diff);
442 })
443 } else {
444 let rem: usize = basek - (log2n % basek);
445 let size: usize = log2n.div_ceil(basek);
446 (1..size).for_each(|i| {
447 if i == size - 1 && rem != basek {
448 let k_rem: usize = basek - rem;
449 izip!(lwe.data.at(0, i).iter(), res.iter_mut()).for_each(|(x, y)| {
450 *y = (*y << k_rem) + (x >> rem);
451 });
452 } else {
453 izip!(lwe.data.at(0, i).iter(), res.iter_mut()).for_each(|(x, y)| {
454 *y = (*y << basek) + x;
455 });
456 }
457 })
458 }
459}
460
461#[inline(always)]
462fn div_round_by_pow2(x: &i64, k: usize) -> i64 {
463 if x >= &0 {
464 (x + (1 << (k - 1))) >> k
465 } else {
466 (x + (-1 << (k - 1))) >> k
467 }
468}
469
470// #[inline(always)]
471// fn div_floor_signed_by_pow2(x: &i64, k: usize) -> i64{
472// let bias: i64 = (1 << k) - 1;
473// (x + ((x >> 63) & bias)) >> k
474// }
475
476#[inline(always)]
477fn div_ceil_signed_by_pow2(x: &i64, k: usize) -> i64 {
478 let bias: i64 = (1 << k) - 1;
479 (x + ((x >> 63) & bias)) >> k
480}