poulpy/spqlios/lib/spqlios/coeffs/coeffs_arithmetic.c
spqlios/lib/spqlios/coeffs/coeffs_arithmetic.cRaw
1#include "coeffs_arithmetic.h"
2
3#include <memory.h>
4#include <assert.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(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(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
133// 0 < p < 2nn
134EXPORT void rnx_automorphism_f64(uint64_t nn, int64_t p, double* res, const double* in) {
135 res[0] = in[0];
136 uint64_t a = 0;
137 uint64_t _2mn = 2 * nn - 1;
138 for (uint64_t i = 1; i < nn; i++) {
139 a = (a + p) & _2mn; // i*p mod 2n
140 if (a < nn) {
141 res[a] = in[i]; // res[ip mod 2n] = res[i]
142 } else {
143 res[a - nn] = -in[i];
144 }
145 }
146}
147
148EXPORT void znx_automorphism_i64(uint64_t nn, int64_t p, int64_t* res, const int64_t* in) {
149 res[0] = in[0];
150 uint64_t a = 0;
151 uint64_t _2mn = 2 * nn - 1;
152 for (uint64_t i = 1; i < nn; i++) {
153 a = (a + p) & _2mn;
154 if (a < nn) {
155 res[a] = in[i]; // res[ip mod 2n] = res[i]
156 } else {
157 res[a - nn] = -in[i];
158 }
159 }
160}
161
162EXPORT void rnx_rotate_inplace_f64(uint64_t nn, int64_t p, double* res) {
163 const uint64_t _2mn = 2 * nn - 1;
164 const uint64_t _mn = nn - 1;
165 uint64_t nb_modif = 0;
166 uint64_t j_start = 0;
167 while (nb_modif < nn) {
168 // follow the cycle that start with j_start
169 uint64_t j = j_start;
170 double tmp1 = res[j];
171 do {
172 // find where the value should go, and with which sign
173 uint64_t new_j = (j + p) & _2mn; // mod 2n to get the position and sign
174 uint64_t new_j_n = new_j & _mn; // mod n to get just the position
175 // exchange this position with tmp1 (and take care of the sign)
176 double tmp2 = res[new_j_n];
177 res[new_j_n] = (new_j < nn) ? tmp1 : -tmp1;
178 tmp1 = tmp2;
179 // move to the new location, and store the number of items modified
180 ++nb_modif;
181 j = new_j_n;
182 } while (j != j_start);
183 // move to the start of the next cycle:
184 // we need to find an index that has not been touched yet, and pick it as next j_start.
185 // in practice, it is enough to do +1, because the group of rotations is cyclic and 1 is a generator.
186 ++j_start;
187 }
188}
189
190EXPORT void znx_rotate_inplace_i64(uint64_t nn, int64_t p, int64_t* 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 int64_t 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 int64_t 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 rnx_mul_xp_minus_one_inplace(uint64_t nn, int64_t p, double* 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 double 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 double tmp2 = res[new_j_n];
233 res[new_j_n] = ((new_j < nn) ? tmp1 : -tmp1) - res[new_j_n];
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
246__always_inline int64_t get_base_k_digit(const int64_t x, const uint64_t base_k) {
247 return (x << (64 - base_k)) >> (64 - base_k);
248}
249
250__always_inline int64_t get_base_k_carry(const int64_t x, const int64_t digit, const uint64_t base_k) {
251 return (x - digit) >> base_k;
252}
253
254EXPORT void znx_normalize(uint64_t nn, uint64_t base_k, int64_t* out, int64_t* carry_out, const int64_t* in,
255 const int64_t* carry_in) {
256 assert(in);
257 if (out != 0) {
258 if (carry_in != 0x0 && carry_out != 0x0) {
259 // with carry in and carry out is computed
260 for (uint64_t i = 0; i < nn; ++i) {
261 const int64_t x = in[i];
262 const int64_t cin = carry_in[i];
263
264 int64_t digit = get_base_k_digit(x, base_k);
265 int64_t carry = get_base_k_carry(x, digit, base_k);
266 int64_t digit_plus_cin = digit + cin;
267 int64_t y = get_base_k_digit(digit_plus_cin, base_k);
268 int64_t cout = carry + get_base_k_carry(digit_plus_cin, y, base_k);
269
270 out[i] = y;
271 carry_out[i] = cout;
272 }
273 } else if (carry_in != 0) {
274 // with carry in and carry out is dropped
275 for (uint64_t i = 0; i < nn; ++i) {
276 const int64_t x = in[i];
277 const int64_t cin = carry_in[i];
278
279 int64_t digit = get_base_k_digit(x, base_k);
280 int64_t digit_plus_cin = digit + cin;
281 int64_t y = get_base_k_digit(digit_plus_cin, base_k);
282
283 out[i] = y;
284 }
285
286 } else if (carry_out != 0) {
287 // no carry in and carry out is computed
288 for (uint64_t i = 0; i < nn; ++i) {
289 const int64_t x = in[i];
290
291 int64_t y = get_base_k_digit(x, base_k);
292 int64_t cout = get_base_k_carry(x, y, base_k);
293
294 out[i] = y;
295 carry_out[i] = cout;
296 }
297
298 } else {
299 // no carry in and carry out is dropped
300 for (uint64_t i = 0; i < nn; ++i) {
301 out[i] = get_base_k_digit(in[i], base_k);
302 }
303 }
304 } else {
305 assert(carry_out);
306 if (carry_in != 0x0) {
307 // with carry in and carry out is computed
308 for (uint64_t i = 0; i < nn; ++i) {
309 const int64_t x = in[i];
310 const int64_t cin = carry_in[i];
311
312 int64_t digit = get_base_k_digit(x, base_k);
313 int64_t carry = get_base_k_carry(x, digit, base_k);
314 int64_t digit_plus_cin = digit + cin;
315 int64_t y = get_base_k_digit(digit_plus_cin, base_k);
316 int64_t cout = carry + get_base_k_carry(digit_plus_cin, y, base_k);
317
318 carry_out[i] = cout;
319 }
320 } else {
321 // no carry in and carry out is computed
322 for (uint64_t i = 0; i < nn; ++i) {
323 const int64_t x = in[i];
324
325 int64_t y = get_base_k_digit(x, base_k);
326 int64_t cout = get_base_k_carry(x, y, base_k);
327
328 carry_out[i] = cout;
329 }
330 }
331 }
332}
333
334void znx_automorphism_inplace_i64(uint64_t nn, int64_t p, int64_t* res) {
335 const uint64_t _2mn = 2 * nn - 1;
336 const uint64_t _mn = nn - 1;
337 const uint64_t m = nn >> 1;
338 // reduce p mod 2n
339 p &= _2mn;
340 // uint64_t vp = p & _2mn;
341 /// uint64_t target_modifs = m >> 1;
342 // we proceed by increasing binary valuation
343 for (uint64_t binval = 1, vp = p & _2mn, orb_size = m; binval < nn;
344 binval <<= 1, vp = (vp << 1) & _2mn, orb_size >>= 1) {
345 // In this loop, we are going to treat the orbit of indexes = binval mod 2.binval.
346 // At the beginning of this loop we have:
347 // vp = binval * p mod 2n
348 // target_modif = m / binval (i.e. order of the orbit binval % 2.binval)
349
350 // first, handle the orders 1 and 2.
351 // if p*binval == binval % 2n: we're done!
352 if (vp == binval) return;
353 // if p*binval == -binval % 2n: nega-mirror the orbit and all the sub-orbits and exit!
354 if (((vp + binval) & _2mn) == 0) {
355 for (uint64_t j = binval; j < m; j += binval) {
356 int64_t tmp = res[j];
357 res[j] = -res[nn - j];
358 res[nn - j] = -tmp;
359 }
360 res[m] = -res[m];
361 return;
362 }
363 // if p*binval == binval + n % 2n: negate the orbit and exit
364 if (((vp - binval) & _mn) == 0) {
365 for (uint64_t j = binval; j < nn; j += 2 * binval) {
366 res[j] = -res[j];
367 }
368 return;
369 }
370 // if p*binval == n - binval % 2n: mirror the orbit and continue!
371 if (((vp + binval) & _mn) == 0) {
372 for (uint64_t j = binval; j < m; j += 2 * binval) {
373 int64_t tmp = res[j];
374 res[j] = res[nn - j];
375 res[nn - j] = tmp;
376 }
377 continue;
378 }
379 // otherwise we will follow the orbit cycles,
380 // starting from binval and -binval in parallel
381 uint64_t j_start = binval;
382 uint64_t nb_modif = 0;
383 while (nb_modif < orb_size) {
384 // follow the cycle that start with j_start
385 uint64_t j = j_start;
386 int64_t tmp1 = res[j];
387 int64_t tmp2 = res[nn - j];
388 do {
389 // find where the value should go, and with which sign
390 uint64_t new_j = (j * p) & _2mn; // mod 2n to get the position and sign
391 uint64_t new_j_n = new_j & _mn; // mod n to get just the position
392 // exchange this position with tmp1 (and take care of the sign)
393 int64_t tmp1a = res[new_j_n];
394 int64_t tmp2a = res[nn - new_j_n];
395 if (new_j < nn) {
396 res[new_j_n] = tmp1;
397 res[nn - new_j_n] = tmp2;
398 } else {
399 res[new_j_n] = -tmp1;
400 res[nn - new_j_n] = -tmp2;
401 }
402 tmp1 = tmp1a;
403 tmp2 = tmp2a;
404 // move to the new location, and store the number of items modified
405 nb_modif += 2;
406 j = new_j_n;
407 } while (j != j_start);
408 // move to the start of the next cycle:
409 // we need to find an index that has not been touched yet, and pick it as next j_start.
410 // in practice, it is enough to do *5, because 5 is a generator.
411 j_start = (5 * j_start) & _mn;
412 }
413 }
414}
415
416void rnx_automorphism_inplace_f64(uint64_t nn, int64_t p, double* res) {
417 const uint64_t _2mn = 2 * nn - 1;
418 const uint64_t _mn = nn - 1;
419 const uint64_t m = nn >> 1;
420 // reduce p mod 2n
421 p &= _2mn;
422 // uint64_t vp = p & _2mn;
423 /// uint64_t target_modifs = m >> 1;
424 // we proceed by increasing binary valuation
425 for (uint64_t binval = 1, vp = p & _2mn, orb_size = m; binval < nn;
426 binval <<= 1, vp = (vp << 1) & _2mn, orb_size >>= 1) {
427 // In this loop, we are going to treat the orbit of indexes = binval mod 2.binval.
428 // At the beginning of this loop we have:
429 // vp = binval * p mod 2n
430 // target_modif = m / binval (i.e. order of the orbit binval % 2.binval)
431
432 // first, handle the orders 1 and 2.
433 // if p*binval == binval % 2n: we're done!
434 if (vp == binval) return;
435 // if p*binval == -binval % 2n: nega-mirror the orbit and all the sub-orbits and exit!
436 if (((vp + binval) & _2mn) == 0) {
437 for (uint64_t j = binval; j < m; j += binval) {
438 double tmp = res[j];
439 res[j] = -res[nn - j];
440 res[nn - j] = -tmp;
441 }
442 res[m] = -res[m];
443 return;
444 }
445 // if p*binval == binval + n % 2n: negate the orbit and exit
446 if (((vp - binval) & _mn) == 0) {
447 for (uint64_t j = binval; j < nn; j += 2 * binval) {
448 res[j] = -res[j];
449 }
450 return;
451 }
452 // if p*binval == n - binval % 2n: mirror the orbit and continue!
453 if (((vp + binval) & _mn) == 0) {
454 for (uint64_t j = binval; j < m; j += 2 * binval) {
455 double tmp = res[j];
456 res[j] = res[nn - j];
457 res[nn - j] = tmp;
458 }
459 continue;
460 }
461 // otherwise we will follow the orbit cycles,
462 // starting from binval and -binval in parallel
463 uint64_t j_start = binval;
464 uint64_t nb_modif = 0;
465 while (nb_modif < orb_size) {
466 // follow the cycle that start with j_start
467 uint64_t j = j_start;
468 double tmp1 = res[j];
469 double tmp2 = res[nn - j];
470 do {
471 // find where the value should go, and with which sign
472 uint64_t new_j = (j * p) & _2mn; // mod 2n to get the position and sign
473 uint64_t new_j_n = new_j & _mn; // mod n to get just the position
474 // exchange this position with tmp1 (and take care of the sign)
475 double tmp1a = res[new_j_n];
476 double tmp2a = res[nn - new_j_n];
477 if (new_j < nn) {
478 res[new_j_n] = tmp1;
479 res[nn - new_j_n] = tmp2;
480 } else {
481 res[new_j_n] = -tmp1;
482 res[nn - new_j_n] = -tmp2;
483 }
484 tmp1 = tmp1a;
485 tmp2 = tmp2a;
486 // move to the new location, and store the number of items modified
487 nb_modif += 2;
488 j = new_j_n;
489 } while (j != j_start);
490 // move to the start of the next cycle:
491 // we need to find an index that has not been touched yet, and pick it as next j_start.
492 // in practice, it is enough to do *5, because 5 is a generator.
493 j_start = (5 * j_start) & _mn;
494 }
495 }
496}