poulpy/poulpy-backend/src/cpu_spqlios/spqlios-arithmetic/spqlios/coeffs/coeffs_arithmetic.c
poulpy-backend/src/cpu_spqlios/spqlios-arithmetic/spqlios/coeffs/coeffs_arithmetic.cRaw
1#include "coeffs_arithmetic.h"
2
3#include <assert.h>
4#include <memory.h>
5
6/** res = a + b */
7EXPORT void znx_add_i64_ref(uint64_t nn, int64_t* res, const int64_t* a, const int64_t* b) {
8 for (uint64_t i = 0; i < nn; ++i) {
9 res[i] = a[i] + b[i];
10 }
11}
12/** res = a - b */
13EXPORT void znx_sub_i64_ref(uint64_t nn, int64_t* res, const int64_t* a, const int64_t* b) {
14 for (uint64_t i = 0; i < nn; ++i) {
15 res[i] = a[i] - b[i];
16 }
17}
18
19EXPORT void znx_negate_i64_ref(uint64_t nn, int64_t* res, const int64_t* a) {
20 for (uint64_t i = 0; i < nn; ++i) {
21 res[i] = -a[i];
22 }
23}
24EXPORT void znx_copy_i64_ref(uint64_t nn, int64_t* res, const int64_t* a) { memcpy(res, a, nn * sizeof(int64_t)); }
25
26EXPORT void znx_zero_i64_ref(uint64_t nn, int64_t* res) { memset(res, 0, nn * sizeof(int64_t)); }
27
28EXPORT void rnx_divide_by_m_ref(uint64_t n, double m, double* res, const double* a) {
29 const double invm = 1. / m;
30 for (uint64_t i = 0; i < n; ++i) {
31 res[i] = a[i] * invm;
32 }
33}
34
35EXPORT void rnx_rotate_f64(uint64_t nn, int64_t p, double* res, const double* in) {
36 uint64_t a = (-p) & (2 * nn - 1); // a= (-p) (pos)mod (2*nn)
37
38 if (a < nn) { // rotate to the left
39 uint64_t nma = nn - a;
40 // rotate first half
41 for (uint64_t j = 0; j < nma; j++) {
42 res[j] = in[j + a];
43 }
44 for (uint64_t j = nma; j < nn; j++) {
45 res[j] = -in[j - nma];
46 }
47 } else {
48 a -= nn;
49 uint64_t nma = nn - a;
50 for (uint64_t j = 0; j < nma; j++) {
51 res[j] = -in[j + a];
52 }
53 for (uint64_t j = nma; j < nn; j++) {
54 // rotate first half
55 res[j] = in[j - nma];
56 }
57 }
58}
59
60EXPORT void znx_rotate_i64(uint64_t nn, int64_t p, int64_t* res, const int64_t* in) {
61 uint64_t a = (-p) & (2 * nn - 1); // a= (-p) (pos)mod (2*nn)
62
63 if (a < nn) { // rotate to the left
64 uint64_t nma = nn - a;
65 // rotate first half
66 for (uint64_t j = 0; j < nma; j++) {
67 res[j] = in[j + a];
68 }
69 for (uint64_t j = nma; j < nn; j++) {
70 res[j] = -in[j - nma];
71 }
72 } else {
73 a -= nn;
74 uint64_t nma = nn - a;
75 for (uint64_t j = 0; j < nma; j++) {
76 res[j] = -in[j + a];
77 }
78 for (uint64_t j = nma; j < nn; j++) {
79 // rotate first half
80 res[j] = in[j - nma];
81 }
82 }
83}
84
85EXPORT void rnx_mul_xp_minus_one_f64(uint64_t nn, int64_t p, double* res, const double* in) {
86 uint64_t a = (-p) & (2 * nn - 1); // a= (-p) (pos)mod (2*nn)
87 if (a < nn) { // rotate to the left
88 uint64_t nma = nn - a;
89 // rotate first half
90 for (uint64_t j = 0; j < nma; j++) {
91 res[j] = in[j + a] - in[j];
92 }
93 for (uint64_t j = nma; j < nn; j++) {
94 res[j] = -in[j - nma] - in[j];
95 }
96 } else {
97 a -= nn;
98 uint64_t nma = nn - a;
99 for (uint64_t j = 0; j < nma; j++) {
100 res[j] = -in[j + a] - in[j];
101 }
102 for (uint64_t j = nma; j < nn; j++) {
103 // rotate first half
104 res[j] = in[j - nma] - in[j];
105 }
106 }
107}
108
109EXPORT void znx_mul_xp_minus_one_i64(uint64_t nn, int64_t p, int64_t* res, const int64_t* in) {
110 uint64_t a = (-p) & (2 * nn - 1); // a= (-p) (pos)mod (2*nn)
111 if (a < nn) { // rotate to the left
112 uint64_t nma = nn - a;
113 // rotate first half
114 for (uint64_t j = 0; j < nma; j++) {
115 res[j] = in[j + a] - in[j];
116 }
117 for (uint64_t j = nma; j < nn; j++) {
118 res[j] = -in[j - nma] - in[j];
119 }
120 } else {
121 a -= nn;
122 uint64_t nma = nn - a;
123 for (uint64_t j = 0; j < nma; j++) {
124 res[j] = -in[j + a] - in[j];
125 }
126 for (uint64_t j = nma; j < nn; j++) {
127 // rotate first half
128 res[j] = in[j - nma] - in[j];
129 }
130 }
131}
132
133EXPORT void znx_mul_xp_minus_one_inplace_i64(uint64_t nn, int64_t p, int64_t* res) {
134 const uint64_t _2mn = 2 * nn - 1;
135 const uint64_t _mn = nn - 1;
136 uint64_t nb_modif = 0;
137 uint64_t j_start = 0;
138 while (nb_modif < nn) {
139 // follow the cycle that start with j_start
140 uint64_t j = j_start;
141 int64_t tmp1 = res[j];
142 do {
143 // find where the value should go, and with which sign
144 uint64_t new_j = (j + p) & _2mn; // mod 2n to get the position and sign
145 uint64_t new_j_n = new_j & _mn; // mod n to get just the position
146 // exchange this position with tmp1 (and take care of the sign)
147 int64_t tmp2 = res[new_j_n];
148 res[new_j_n] = ((new_j < nn) ? tmp1 : -tmp1) - res[new_j_n];
149 tmp1 = tmp2;
150 // move to the new location, and store the number of items modified
151 ++nb_modif;
152 j = new_j_n;
153 } while (j != j_start);
154 // move to the start of the next cycle:
155 // we need to find an index that has not been touched yet, and pick it as next j_start.
156 // in practice, it is enough to do +1, because the group of rotations is cyclic and 1 is a generator.
157 ++j_start;
158 }
159}
160
161// 0 < p < 2nn
162EXPORT void rnx_automorphism_f64(uint64_t nn, int64_t p, double* res, const double* in) {
163 res[0] = in[0];
164 uint64_t a = 0;
165 uint64_t _2mn = 2 * nn - 1;
166 for (uint64_t i = 1; i < nn; i++) {
167 a = (a + p) & _2mn; // i*p mod 2n
168 if (a < nn) {
169 res[a] = in[i]; // res[ip mod 2n] = res[i]
170 } else {
171 res[a - nn] = -in[i];
172 }
173 }
174}
175
176EXPORT void znx_automorphism_i64(uint64_t nn, int64_t p, int64_t* res, const int64_t* in) {
177 res[0] = in[0];
178 uint64_t a = 0;
179 uint64_t _2mn = 2 * nn - 1;
180 for (uint64_t i = 1; i < nn; i++) {
181 a = (a + p) & _2mn;
182 if (a < nn) {
183 res[a] = in[i]; // res[ip mod 2n] = res[i]
184 } else {
185 res[a - nn] = -in[i];
186 }
187 }
188}
189
190EXPORT void rnx_rotate_inplace_f64(uint64_t nn, int64_t p, double* res) {
191 const uint64_t _2mn = 2 * nn - 1;
192 const uint64_t _mn = nn - 1;
193 uint64_t nb_modif = 0;
194 uint64_t j_start = 0;
195 while (nb_modif < nn) {
196 // follow the cycle that start with j_start
197 uint64_t j = j_start;
198 double tmp1 = res[j];
199 do {
200 // find where the value should go, and with which sign
201 uint64_t new_j = (j + p) & _2mn; // mod 2n to get the position and sign
202 uint64_t new_j_n = new_j & _mn; // mod n to get just the position
203 // exchange this position with tmp1 (and take care of the sign)
204 double tmp2 = res[new_j_n];
205 res[new_j_n] = (new_j < nn) ? tmp1 : -tmp1;
206 tmp1 = tmp2;
207 // move to the new location, and store the number of items modified
208 ++nb_modif;
209 j = new_j_n;
210 } while (j != j_start);
211 // move to the start of the next cycle:
212 // we need to find an index that has not been touched yet, and pick it as next j_start.
213 // in practice, it is enough to do +1, because the group of rotations is cyclic and 1 is a generator.
214 ++j_start;
215 }
216}
217
218EXPORT void znx_rotate_inplace_i64(uint64_t nn, int64_t p, int64_t* res) {
219 const uint64_t _2mn = 2 * nn - 1;
220 const uint64_t _mn = nn - 1;
221 uint64_t nb_modif = 0;
222 uint64_t j_start = 0;
223 while (nb_modif < nn) {
224 // follow the cycle that start with j_start
225 uint64_t j = j_start;
226 int64_t tmp1 = res[j];
227 do {
228 // find where the value should go, and with which sign
229 uint64_t new_j = (j + p) & _2mn; // mod 2n to get the position and sign
230 uint64_t new_j_n = new_j & _mn; // mod n to get just the position
231 // exchange this position with tmp1 (and take care of the sign)
232 int64_t tmp2 = res[new_j_n];
233 res[new_j_n] = (new_j < nn) ? tmp1 : -tmp1;
234 tmp1 = tmp2;
235 // move to the new location, and store the number of items modified
236 ++nb_modif;
237 j = new_j_n;
238 } while (j != j_start);
239 // move to the start of the next cycle:
240 // we need to find an index that has not been touched yet, and pick it as next j_start.
241 // in practice, it is enough to do +1, because the group of rotations is cyclic and 1 is a generator.
242 ++j_start;
243 }
244}
245
246EXPORT void rnx_mul_xp_minus_one_inplace_f64(uint64_t nn, int64_t p, double* res) {
247 const uint64_t _2mn = 2 * nn - 1;
248 const uint64_t _mn = nn - 1;
249 uint64_t nb_modif = 0;
250 uint64_t j_start = 0;
251 while (nb_modif < nn) {
252 // follow the cycle that start with j_start
253 uint64_t j = j_start;
254 double tmp1 = res[j];
255 do {
256 // find where the value should go, and with which sign
257 uint64_t new_j = (j + p) & _2mn; // mod 2n to get the position and sign
258 uint64_t new_j_n = new_j & _mn; // mod n to get just the position
259 // exchange this position with tmp1 (and take care of the sign)
260 double tmp2 = res[new_j_n];
261 res[new_j_n] = ((new_j < nn) ? tmp1 : -tmp1) - res[new_j_n];
262 tmp1 = tmp2;
263 // move to the new location, and store the number of items modified
264 ++nb_modif;
265 j = new_j_n;
266 } while (j != j_start);
267 // move to the start of the next cycle:
268 // we need to find an index that has not been touched yet, and pick it as next j_start.
269 // in practice, it is enough to do +1, because the group of rotations is cyclic and 1 is a generator.
270 ++j_start;
271 }
272}
273
274__always_inline int64_t get_base_k_digit(const int64_t x, const uint64_t base_k) {
275 return (x << (64 - base_k)) >> (64 - base_k);
276}
277
278__always_inline int64_t get_base_k_carry(const int64_t x, const int64_t digit, const uint64_t base_k) {
279 return (x - digit) >> base_k;
280}
281
282EXPORT void znx_normalize(uint64_t nn, uint64_t base_k, int64_t* out, int64_t* carry_out, const int64_t* in,
283 const int64_t* carry_in) {
284 assert(in);
285 if (out != 0) {
286 if (carry_in != 0x0 && carry_out != 0x0) {
287 // with carry in and carry out is computed
288 for (uint64_t i = 0; i < nn; ++i) {
289 const int64_t x = in[i];
290 const int64_t cin = carry_in[i];
291
292 int64_t digit = get_base_k_digit(x, base_k);
293 int64_t carry = get_base_k_carry(x, digit, base_k);
294 int64_t digit_plus_cin = digit + cin;
295 int64_t y = get_base_k_digit(digit_plus_cin, base_k);
296 int64_t cout = carry + get_base_k_carry(digit_plus_cin, y, base_k);
297
298 out[i] = y;
299 carry_out[i] = cout;
300 }
301 } else if (carry_in != 0) {
302 // with carry in and carry out is dropped
303 for (uint64_t i = 0; i < nn; ++i) {
304 const int64_t x = in[i];
305 const int64_t cin = carry_in[i];
306
307 int64_t digit = get_base_k_digit(x, base_k);
308 int64_t digit_plus_cin = digit + cin;
309 int64_t y = get_base_k_digit(digit_plus_cin, base_k);
310
311 out[i] = y;
312 }
313
314 } else if (carry_out != 0) {
315 // no carry in and carry out is computed
316 for (uint64_t i = 0; i < nn; ++i) {
317 const int64_t x = in[i];
318
319 int64_t y = get_base_k_digit(x, base_k);
320 int64_t cout = get_base_k_carry(x, y, base_k);
321
322 out[i] = y;
323 carry_out[i] = cout;
324 }
325
326 } else {
327 // no carry in and carry out is dropped
328 for (uint64_t i = 0; i < nn; ++i) {
329 out[i] = get_base_k_digit(in[i], base_k);
330 }
331 }
332 } else {
333 assert(carry_out);
334 if (carry_in != 0x0) {
335 // with carry in and carry out is computed
336 for (uint64_t i = 0; i < nn; ++i) {
337 const int64_t x = in[i];
338 const int64_t cin = carry_in[i];
339
340 int64_t digit = get_base_k_digit(x, base_k);
341 int64_t carry = get_base_k_carry(x, digit, base_k);
342 int64_t digit_plus_cin = digit + cin;
343 int64_t y = get_base_k_digit(digit_plus_cin, base_k);
344 int64_t cout = carry + get_base_k_carry(digit_plus_cin, y, base_k);
345
346 carry_out[i] = cout;
347 }
348 } else {
349 // no carry in and carry out is computed
350 for (uint64_t i = 0; i < nn; ++i) {
351 const int64_t x = in[i];
352
353 int64_t y = get_base_k_digit(x, base_k);
354 int64_t cout = get_base_k_carry(x, y, base_k);
355
356 carry_out[i] = cout;
357 }
358 }
359 }
360}
361
362void znx_automorphism_inplace_i64(uint64_t nn, int64_t p, int64_t* res) {
363 const uint64_t _2mn = 2 * nn - 1;
364 const uint64_t _mn = nn - 1;
365 const uint64_t m = nn >> 1;
366 // reduce p mod 2n
367 p &= _2mn;
368 // uint64_t vp = p & _2mn;
369 /// uint64_t target_modifs = m >> 1;
370 // we proceed by increasing binary valuation
371 for (uint64_t binval = 1, vp = p & _2mn, orb_size = m; binval < nn;
372 binval <<= 1, vp = (vp << 1) & _2mn, orb_size >>= 1) {
373 // In this loop, we are going to treat the orbit of indexes = binval mod 2.binval.
374 // At the beginning of this loop we have:
375 // vp = binval * p mod 2n
376 // target_modif = m / binval (i.e. order of the orbit binval % 2.binval)
377
378 // first, handle the orders 1 and 2.
379 // if p*binval == binval % 2n: we're done!
380 if (vp == binval) return;
381 // if p*binval == -binval % 2n: nega-mirror the orbit and all the sub-orbits and exit!
382 if (((vp + binval) & _2mn) == 0) {
383 for (uint64_t j = binval; j < m; j += binval) {
384 int64_t tmp = res[j];
385 res[j] = -res[nn - j];
386 res[nn - j] = -tmp;
387 }
388 res[m] = -res[m];
389 return;
390 }
391 // if p*binval == binval + n % 2n: negate the orbit and exit
392 if (((vp - binval) & _mn) == 0) {
393 for (uint64_t j = binval; j < nn; j += 2 * binval) {
394 res[j] = -res[j];
395 }
396 return;
397 }
398 // if p*binval == n - binval % 2n: mirror the orbit and continue!
399 if (((vp + binval) & _mn) == 0) {
400 for (uint64_t j = binval; j < m; j += 2 * binval) {
401 int64_t tmp = res[j];
402 res[j] = res[nn - j];
403 res[nn - j] = tmp;
404 }
405 continue;
406 }
407 // otherwise we will follow the orbit cycles,
408 // starting from binval and -binval in parallel
409 uint64_t j_start = binval;
410 uint64_t nb_modif = 0;
411 while (nb_modif < orb_size) {
412 // follow the cycle that start with j_start
413 uint64_t j = j_start;
414 int64_t tmp1 = res[j];
415 int64_t tmp2 = res[nn - j];
416 do {
417 // find where the value should go, and with which sign
418 uint64_t new_j = (j * p) & _2mn; // mod 2n to get the position and sign
419 uint64_t new_j_n = new_j & _mn; // mod n to get just the position
420 // exchange this position with tmp1 (and take care of the sign)
421 int64_t tmp1a = res[new_j_n];
422 int64_t tmp2a = res[nn - new_j_n];
423 if (new_j < nn) {
424 res[new_j_n] = tmp1;
425 res[nn - new_j_n] = tmp2;
426 } else {
427 res[new_j_n] = -tmp1;
428 res[nn - new_j_n] = -tmp2;
429 }
430 tmp1 = tmp1a;
431 tmp2 = tmp2a;
432 // move to the new location, and store the number of items modified
433 nb_modif += 2;
434 j = new_j_n;
435 } while (j != j_start);
436 // move to the start of the next cycle:
437 // we need to find an index that has not been touched yet, and pick it as next j_start.
438 // in practice, it is enough to do *5, because 5 is a generator.
439 j_start = (5 * j_start) & _mn;
440 }
441 }
442}
443
444void rnx_automorphism_inplace_f64(uint64_t nn, int64_t p, double* res) {
445 const uint64_t _2mn = 2 * nn - 1;
446 const uint64_t _mn = nn - 1;
447 const uint64_t m = nn >> 1;
448 // reduce p mod 2n
449 p &= _2mn;
450 // uint64_t vp = p & _2mn;
451 /// uint64_t target_modifs = m >> 1;
452 // we proceed by increasing binary valuation
453 for (uint64_t binval = 1, vp = p & _2mn, orb_size = m; binval < nn;
454 binval <<= 1, vp = (vp << 1) & _2mn, orb_size >>= 1) {
455 // In this loop, we are going to treat the orbit of indexes = binval mod 2.binval.
456 // At the beginning of this loop we have:
457 // vp = binval * p mod 2n
458 // target_modif = m / binval (i.e. order of the orbit binval % 2.binval)
459
460 // first, handle the orders 1 and 2.
461 // if p*binval == binval % 2n: we're done!
462 if (vp == binval) return;
463 // if p*binval == -binval % 2n: nega-mirror the orbit and all the sub-orbits and exit!
464 if (((vp + binval) & _2mn) == 0) {
465 for (uint64_t j = binval; j < m; j += binval) {
466 double tmp = res[j];
467 res[j] = -res[nn - j];
468 res[nn - j] = -tmp;
469 }
470 res[m] = -res[m];
471 return;
472 }
473 // if p*binval == binval + n % 2n: negate the orbit and exit
474 if (((vp - binval) & _mn) == 0) {
475 for (uint64_t j = binval; j < nn; j += 2 * binval) {
476 res[j] = -res[j];
477 }
478 return;
479 }
480 // if p*binval == n - binval % 2n: mirror the orbit and continue!
481 if (((vp + binval) & _mn) == 0) {
482 for (uint64_t j = binval; j < m; j += 2 * binval) {
483 double tmp = res[j];
484 res[j] = res[nn - j];
485 res[nn - j] = tmp;
486 }
487 continue;
488 }
489 // otherwise we will follow the orbit cycles,
490 // starting from binval and -binval in parallel
491 uint64_t j_start = binval;
492 uint64_t nb_modif = 0;
493 while (nb_modif < orb_size) {
494 // follow the cycle that start with j_start
495 uint64_t j = j_start;
496 double tmp1 = res[j];
497 double tmp2 = res[nn - j];
498 do {
499 // find where the value should go, and with which sign
500 uint64_t new_j = (j * p) & _2mn; // mod 2n to get the position and sign
501 uint64_t new_j_n = new_j & _mn; // mod n to get just the position
502 // exchange this position with tmp1 (and take care of the sign)
503 double tmp1a = res[new_j_n];
504 double tmp2a = res[nn - new_j_n];
505 if (new_j < nn) {
506 res[new_j_n] = tmp1;
507 res[nn - new_j_n] = tmp2;
508 } else {
509 res[new_j_n] = -tmp1;
510 res[nn - new_j_n] = -tmp2;
511 }
512 tmp1 = tmp1a;
513 tmp2 = tmp2a;
514 // move to the new location, and store the number of items modified
515 nb_modif += 2;
516 j = new_j_n;
517 } while (j != j_start);
518 // move to the start of the next cycle:
519 // we need to find an index that has not been touched yet, and pick it as next j_start.
520 // in practice, it is enough to do *5, because 5 is a generator.
521 j_start = (5 * j_start) & _mn;
522 }
523 }
524}