poulpy/poulpy-backend/src/cpu_spqlios/spqlios-arithmetic/spqlios/q120/q120_arithmetic_avx2.c
poulpy-backend/src/cpu_spqlios/spqlios-arithmetic/spqlios/q120/q120_arithmetic_avx2.cRaw
1#include <immintrin.h>
2
3#include "q120_arithmetic.h"
4#include "q120_arithmetic_private.h"
5
6EXPORT void q120_vec_mat1col_product_baa_avx2(q120_mat1col_product_baa_precomp* precomp, const uint64_t ell,
7 q120b* const res, const q120a* const x, const q120a* const y) {
8 /**
9 * Algorithm:
10 * - res = acc1 + acc2 . ((2^H) % Q)
11 * - acc1 is the sum of H LSB of products x[i].y[i]
12 * - acc2 is the sum of 64-H MSB of products x[i]].y[i]
13 * - for l < 10k acc1 will have H + log2(10000) and acc2 64 - H + log2(10000) bits
14 * - final sum has max(H, 64 - H + bit_size((2^H) % Q)) + log2(10000) + 1 bits
15 */
16
17 const uint64_t H = precomp->h;
18 const __m256i MASK = _mm256_set1_epi64x((UINT64_C(1) << H) - 1);
19
20 __m256i acc1 = _mm256_setzero_si256();
21 __m256i acc2 = _mm256_setzero_si256();
22
23 const __m256i* x_ptr = (__m256i*)x;
24 const __m256i* y_ptr = (__m256i*)y;
25
26 for (uint64_t i = 0; i < ell; ++i) {
27 __m256i a = _mm256_loadu_si256(x_ptr);
28 __m256i b = _mm256_loadu_si256(y_ptr);
29 __m256i t = _mm256_mul_epu32(a, b);
30
31 acc1 = _mm256_add_epi64(acc1, _mm256_and_si256(t, MASK));
32 acc2 = _mm256_add_epi64(acc2, _mm256_srli_epi64(t, H));
33
34 x_ptr++;
35 y_ptr++;
36 }
37
38 const __m256i H_POW_RED = _mm256_loadu_si256((__m256i*)precomp->h_pow_red);
39
40 __m256i t = _mm256_add_epi64(acc1, _mm256_mul_epu32(acc2, H_POW_RED));
41 _mm256_storeu_si256((__m256i*)res, t);
42}
43
44EXPORT void q120_vec_mat1col_product_bbb_avx2(q120_mat1col_product_bbb_precomp* precomp, const uint64_t ell,
45 q120b* const res, const q120b* const x, const q120b* const y) {
46 /**
47 * Algorithm:
48 * 1. Split x_i and y_i in 2 32-bit parts and compute the cross-products:
49 * - x_i = xl_i + xh_i . 2^32
50 * - y_i = yl_i + yh_i . 2^32
51 * - A_i = xl_i . yl_i
52 * - B_i = xl_i . yh_i
53 * - C_i = xh_i . yl_i
54 * - D_i = xh_i . yh_i
55 * - we have x_i . y_i == A_i + (B_i + C_i) . 2^32 + D_i . 2^64
56 * 2. Split A_i, B_i, C_i and D_i into 2 32-bit parts
57 * - A_i = Al_i + Ah_i . 2^32
58 * - B_i = Bl_i + Bh_i . 2^32
59 * - C_i = Cl_i + Ch_i . 2^32
60 * - D_i = Dl_i + Dh_i . 2^32
61 * 3. Compute the sums:
62 * - S1 = \sum Al_i
63 * - S2 = \sum (Ah_i + Bl_i + Cl_i)
64 * - S3 = \sum (Bh_i + Ch_i + Dl_i)
65 * - S4 = \sum Dh_i
66 * - here S1, S4 have 32 + log2(ell) bits and S2, S3 have 32 + log2(ell) +
67 * log2(3) bits
68 * - for ell == 10000 S2, S3 have < 47 bits
69 * 4. Split S1, S2, S3 and S4 in 2 24-bit parts (24 = ceil(47/2))
70 * - S1 = S1l + S1h . 2^24
71 * - S2 = S2l + S2h . 2^24
72 * - S3 = S3l + S3h . 2^24
73 * - S4 = S4l + S4h . 2^24
74 * 5. Compute final result as:
75 * - \sum x_i . y_i = S1l + S1h . 2^24
76 * + S2l . 2^32 + S2h . 2^(32+24)
77 * + S3l . 2^64 + S3h . 2^(64 + 24)
78 * + S4l . 2^96 + S4l . 2^(96+24)
79 * - here the powers of 2 are reduced modulo the primes Q before
80 * multiplications
81 * - the result will be on 24 + 3 + bit size of primes Q
82 */
83 const uint64_t H1 = 32;
84 const __m256i MASK1 = _mm256_set1_epi64x((UINT64_C(1) << H1) - 1);
85
86 __m256i s1 = _mm256_setzero_si256();
87 __m256i s2 = _mm256_setzero_si256();
88 __m256i s3 = _mm256_setzero_si256();
89 __m256i s4 = _mm256_setzero_si256();
90
91 const __m256i* x_ptr = (__m256i*)x;
92 const __m256i* y_ptr = (__m256i*)y;
93
94 for (uint64_t i = 0; i < ell; ++i) {
95 __m256i x = _mm256_loadu_si256(x_ptr);
96 __m256i xl = _mm256_and_si256(x, MASK1);
97 __m256i xh = _mm256_srli_epi64(x, H1);
98
99 __m256i y = _mm256_loadu_si256(y_ptr);
100 __m256i yl = _mm256_and_si256(y, MASK1);
101 __m256i yh = _mm256_srli_epi64(y, H1);
102
103 __m256i a = _mm256_mul_epu32(xl, yl);
104 __m256i b = _mm256_mul_epu32(xl, yh);
105 __m256i c = _mm256_mul_epu32(xh, yl);
106 __m256i d = _mm256_mul_epu32(xh, yh);
107
108 s1 = _mm256_add_epi64(s1, _mm256_and_si256(a, MASK1));
109
110 s2 = _mm256_add_epi64(s2, _mm256_srli_epi64(a, H1));
111 s2 = _mm256_add_epi64(s2, _mm256_and_si256(b, MASK1));
112 s2 = _mm256_add_epi64(s2, _mm256_and_si256(c, MASK1));
113
114 s3 = _mm256_add_epi64(s3, _mm256_srli_epi64(b, H1));
115 s3 = _mm256_add_epi64(s3, _mm256_srli_epi64(c, H1));
116 s3 = _mm256_add_epi64(s3, _mm256_and_si256(d, MASK1));
117
118 s4 = _mm256_add_epi64(s4, _mm256_srli_epi64(d, H1));
119
120 x_ptr++;
121 y_ptr++;
122 }
123
124 const uint64_t H2 = precomp->h;
125 const __m256i MASK2 = _mm256_set1_epi64x((UINT64_C(1) << H2) - 1);
126
127 const __m256i S1H_POW_RED = _mm256_loadu_si256((__m256i*)precomp->s1h_pow_red);
128 __m256i s1l = _mm256_and_si256(s1, MASK2);
129 __m256i s1h = _mm256_srli_epi64(s1, H2);
130 __m256i t = _mm256_add_epi64(s1l, _mm256_mul_epu32(s1h, S1H_POW_RED));
131
132 const __m256i S2L_POW_RED = _mm256_loadu_si256((__m256i*)precomp->s2l_pow_red);
133 const __m256i S2H_POW_RED = _mm256_loadu_si256((__m256i*)precomp->s2h_pow_red);
134 __m256i s2l = _mm256_and_si256(s2, MASK2);
135 __m256i s2h = _mm256_srli_epi64(s2, H2);
136 t = _mm256_add_epi64(t, _mm256_mul_epu32(s2l, S2L_POW_RED));
137 t = _mm256_add_epi64(t, _mm256_mul_epu32(s2h, S2H_POW_RED));
138
139 const __m256i S3L_POW_RED = _mm256_loadu_si256((__m256i*)precomp->s3l_pow_red);
140 const __m256i S3H_POW_RED = _mm256_loadu_si256((__m256i*)precomp->s3h_pow_red);
141 __m256i s3l = _mm256_and_si256(s3, MASK2);
142 __m256i s3h = _mm256_srli_epi64(s3, H2);
143 t = _mm256_add_epi64(t, _mm256_mul_epu32(s3l, S3L_POW_RED));
144 t = _mm256_add_epi64(t, _mm256_mul_epu32(s3h, S3H_POW_RED));
145
146 const __m256i S4L_POW_RED = _mm256_loadu_si256((__m256i*)precomp->s4l_pow_red);
147 const __m256i S4H_POW_RED = _mm256_loadu_si256((__m256i*)precomp->s4h_pow_red);
148 __m256i s4l = _mm256_and_si256(s4, MASK2);
149 __m256i s4h = _mm256_srli_epi64(s4, H2);
150 t = _mm256_add_epi64(t, _mm256_mul_epu32(s4l, S4L_POW_RED));
151 t = _mm256_add_epi64(t, _mm256_mul_epu32(s4h, S4H_POW_RED));
152
153 _mm256_storeu_si256((__m256i*)res, t);
154}
155
156EXPORT void q120_vec_mat1col_product_bbc_avx2(q120_mat1col_product_bbc_precomp* precomp, const uint64_t ell,
157 q120b* const res, const q120b* const x, const q120c* const y) {
158 /**
159 * Algorithm:
160 * 0. We have
161 * - y0_i == y_i % Q and y1_i == (y_i . 2^32) % Q
162 * 1. Split x_i in 2 32-bit parts and compute the cross-products:
163 * - x_i = xl_i + xh_i . 2^32
164 * - A_i = xl_i . y1_i
165 * - B_i = xh_i . y2_i
166 * - we have x_i . y_i == A_i + B_i
167 * 2. Split A_i and B_i into 2 32-bit parts
168 * - A_i = Al_i + Ah_i . 2^32
169 * - B_i = Bl_i + Bh_i . 2^32
170 * 3. Compute the sums:
171 * - S1 = \sum Al_i + Bl_i
172 * - S2 = \sum Ah_i + Bh_i
173 * - here S1 and S2 have 32 + log2(ell) bits
174 * - for ell == 10000 S1, S2 have < 46 bits
175 * 4. Split S2 in 27-bit and 19-bit parts (27+19 == 46)
176 * - S2 = S2l + S2h . 2^27
177 * 5. Compute final result as:
178 * - \sum x_i . y_i = S1 + S2l . 2^32 + S2h . 2^(32+27)
179 * - here the powers of 2 are reduced modulo the primes Q before
180 * multiplications
181 * - the result will be on < 52 bits
182 */
183
184 const uint64_t H1 = 32;
185 const __m256i MASK1 = _mm256_set1_epi64x((UINT64_C(1) << H1) - 1);
186
187 __m256i s1 = _mm256_setzero_si256();
188 __m256i s2 = _mm256_setzero_si256();
189
190 const __m256i* x_ptr = (__m256i*)x;
191 const __m256i* y_ptr = (__m256i*)y;
192
193 for (uint64_t i = 0; i < ell; ++i) {
194 __m256i x = _mm256_loadu_si256(x_ptr);
195 __m256i xl = _mm256_and_si256(x, MASK1);
196 __m256i xh = _mm256_srli_epi64(x, H1);
197
198 __m256i y = _mm256_loadu_si256(y_ptr);
199 __m256i y0 = _mm256_and_si256(y, MASK1);
200 __m256i y1 = _mm256_srli_epi64(y, H1);
201
202 __m256i a = _mm256_mul_epu32(xl, y0);
203 __m256i b = _mm256_mul_epu32(xh, y1);
204
205 s1 = _mm256_add_epi64(s1, _mm256_and_si256(a, MASK1));
206 s1 = _mm256_add_epi64(s1, _mm256_and_si256(b, MASK1));
207
208 s2 = _mm256_add_epi64(s2, _mm256_srli_epi64(a, H1));
209 s2 = _mm256_add_epi64(s2, _mm256_srli_epi64(b, H1));
210
211 x_ptr++;
212 y_ptr++;
213 }
214
215 const uint64_t H2 = precomp->h;
216 const __m256i MASK2 = _mm256_set1_epi64x((UINT64_C(1) << H2) - 1);
217
218 __m256i t = s1;
219
220 const __m256i S2L_POW_RED = _mm256_loadu_si256((__m256i*)precomp->s2l_pow_red);
221 const __m256i S2H_POW_RED = _mm256_loadu_si256((__m256i*)precomp->s2h_pow_red);
222 __m256i s2l = _mm256_and_si256(s2, MASK2);
223 __m256i s2h = _mm256_srli_epi64(s2, H2);
224 t = _mm256_add_epi64(t, _mm256_mul_epu32(s2l, S2L_POW_RED));
225 t = _mm256_add_epi64(t, _mm256_mul_epu32(s2h, S2H_POW_RED));
226
227 _mm256_storeu_si256((__m256i*)res, t);
228}
229
230/**
231 * @deprecated keeping this one for history only.
232 * There is a slight register starvation condition on the q120x2_vec_mat2cols
233 * strategy below sounds better.
234 */
235EXPORT void q120x2_vec_mat2cols_product_bbc_avx2_old(q120_mat1col_product_bbc_precomp* precomp, const uint64_t ell,
236 q120b* const res, const q120b* const x, const q120c* const y) {
237 __m256i s0 = _mm256_setzero_si256(); // col 1a
238 __m256i s1 = _mm256_setzero_si256();
239 __m256i s2 = _mm256_setzero_si256(); // col 1b
240 __m256i s3 = _mm256_setzero_si256();
241 __m256i s4 = _mm256_setzero_si256(); // col 2a
242 __m256i s5 = _mm256_setzero_si256();
243 __m256i s6 = _mm256_setzero_si256(); // col 2b
244 __m256i s7 = _mm256_setzero_si256();
245 __m256i s8, s9, s10, s11;
246 __m256i s12, s13, s14, s15;
247
248 const __m256i* x_ptr = (__m256i*)x;
249 const __m256i* y_ptr = (__m256i*)y;
250 __m256i* res_ptr = (__m256i*)res;
251 for (uint64_t i = 0; i < ell; ++i) {
252 s8 = _mm256_loadu_si256(x_ptr);
253 s9 = _mm256_loadu_si256(x_ptr + 1);
254 s10 = _mm256_srli_epi64(s8, 32);
255 s11 = _mm256_srli_epi64(s9, 32);
256
257 s12 = _mm256_loadu_si256(y_ptr);
258 s13 = _mm256_loadu_si256(y_ptr + 1);
259 s14 = _mm256_srli_epi64(s12, 32);
260 s15 = _mm256_srli_epi64(s13, 32);
261
262 s12 = _mm256_mul_epu32(s8, s12); // -> s0,s1
263 s13 = _mm256_mul_epu32(s9, s13); // -> s2,s3
264 s14 = _mm256_mul_epu32(s10, s14); // -> s0,s1
265 s15 = _mm256_mul_epu32(s11, s15); // -> s2,s3
266
267 s10 = _mm256_slli_epi64(s12, 32); // -> s0
268 s11 = _mm256_slli_epi64(s13, 32); // -> s2
269 s12 = _mm256_srli_epi64(s12, 32); // -> s1
270 s13 = _mm256_srli_epi64(s13, 32); // -> s3
271 s10 = _mm256_srli_epi64(s10, 32); // -> s0
272 s11 = _mm256_srli_epi64(s11, 32); // -> s2
273
274 s0 = _mm256_add_epi64(s0, s10);
275 s1 = _mm256_add_epi64(s1, s12);
276 s2 = _mm256_add_epi64(s2, s11);
277 s3 = _mm256_add_epi64(s3, s13);
278
279 s10 = _mm256_slli_epi64(s14, 32); // -> s0
280 s11 = _mm256_slli_epi64(s15, 32); // -> s2
281 s14 = _mm256_srli_epi64(s14, 32); // -> s1
282 s15 = _mm256_srli_epi64(s15, 32); // -> s3
283 s10 = _mm256_srli_epi64(s10, 32); // -> s0
284 s11 = _mm256_srli_epi64(s11, 32); // -> s2
285
286 s0 = _mm256_add_epi64(s0, s10);
287 s1 = _mm256_add_epi64(s1, s14);
288 s2 = _mm256_add_epi64(s2, s11);
289 s3 = _mm256_add_epi64(s3, s15);
290
291 // deal with the second column
292 // s8,s9 are still in place!
293 s10 = _mm256_srli_epi64(s8, 32);
294 s11 = _mm256_srli_epi64(s9, 32);
295
296 s12 = _mm256_loadu_si256(y_ptr + 2);
297 s13 = _mm256_loadu_si256(y_ptr + 3);
298 s14 = _mm256_srli_epi64(s12, 32);
299 s15 = _mm256_srli_epi64(s13, 32);
300
301 s12 = _mm256_mul_epu32(s8, s12); // -> s4,s5
302 s13 = _mm256_mul_epu32(s9, s13); // -> s6,s7
303 s14 = _mm256_mul_epu32(s10, s14); // -> s4,s5
304 s15 = _mm256_mul_epu32(s11, s15); // -> s6,s7
305
306 s10 = _mm256_slli_epi64(s12, 32); // -> s4
307 s11 = _mm256_slli_epi64(s13, 32); // -> s6
308 s12 = _mm256_srli_epi64(s12, 32); // -> s5
309 s13 = _mm256_srli_epi64(s13, 32); // -> s7
310 s10 = _mm256_srli_epi64(s10, 32); // -> s4
311 s11 = _mm256_srli_epi64(s11, 32); // -> s6
312
313 s4 = _mm256_add_epi64(s4, s10);
314 s5 = _mm256_add_epi64(s5, s12);
315 s6 = _mm256_add_epi64(s6, s11);
316 s7 = _mm256_add_epi64(s7, s13);
317
318 s10 = _mm256_slli_epi64(s14, 32); // -> s4
319 s11 = _mm256_slli_epi64(s15, 32); // -> s6
320 s14 = _mm256_srli_epi64(s14, 32); // -> s5
321 s15 = _mm256_srli_epi64(s15, 32); // -> s7
322 s10 = _mm256_srli_epi64(s10, 32); // -> s4
323 s11 = _mm256_srli_epi64(s11, 32); // -> s6
324
325 s4 = _mm256_add_epi64(s4, s10);
326 s5 = _mm256_add_epi64(s5, s14);
327 s6 = _mm256_add_epi64(s6, s11);
328 s7 = _mm256_add_epi64(s7, s15);
329
330 x_ptr += 2;
331 y_ptr += 4;
332 }
333 // final reduction
334 const uint64_t H2 = precomp->h;
335 s8 = _mm256_set1_epi64x((UINT64_C(1) << H2) - 1); // MASK2
336 s9 = _mm256_loadu_si256((__m256i*)precomp->s2l_pow_red); // S2L_POW_RED
337 s10 = _mm256_loadu_si256((__m256i*)precomp->s2h_pow_red); // S2H_POW_RED
338 //--- s0,s1
339 s11 = _mm256_and_si256(s1, s8);
340 s12 = _mm256_srli_epi64(s1, H2);
341 s13 = _mm256_mul_epu32(s11, s9);
342 s14 = _mm256_mul_epu32(s12, s10);
343 s0 = _mm256_add_epi64(s0, s13);
344 s0 = _mm256_add_epi64(s0, s14);
345 _mm256_storeu_si256(res_ptr + 0, s0);
346 //--- s2,s3
347 s11 = _mm256_and_si256(s3, s8);
348 s12 = _mm256_srli_epi64(s3, H2);
349 s13 = _mm256_mul_epu32(s11, s9);
350 s14 = _mm256_mul_epu32(s12, s10);
351 s2 = _mm256_add_epi64(s2, s13);
352 s2 = _mm256_add_epi64(s2, s14);
353 _mm256_storeu_si256(res_ptr + 1, s2);
354 //--- s4,s5
355 s11 = _mm256_and_si256(s5, s8);
356 s12 = _mm256_srli_epi64(s5, H2);
357 s13 = _mm256_mul_epu32(s11, s9);
358 s14 = _mm256_mul_epu32(s12, s10);
359 s4 = _mm256_add_epi64(s4, s13);
360 s4 = _mm256_add_epi64(s4, s14);
361 _mm256_storeu_si256(res_ptr + 2, s4);
362 //--- s6,s7
363 s11 = _mm256_and_si256(s7, s8);
364 s12 = _mm256_srli_epi64(s7, H2);
365 s13 = _mm256_mul_epu32(s11, s9);
366 s14 = _mm256_mul_epu32(s12, s10);
367 s6 = _mm256_add_epi64(s6, s13);
368 s6 = _mm256_add_epi64(s6, s14);
369 _mm256_storeu_si256(res_ptr + 3, s6);
370}
371
372EXPORT void q120x2_vec_mat2cols_product_bbc_avx2(q120_mat1col_product_bbc_precomp* precomp, const uint64_t ell,
373 q120b* const res, const q120b* const x, const q120c* const y) {
374 __m256i s0 = _mm256_setzero_si256(); // col 1a
375 __m256i s1 = _mm256_setzero_si256();
376 __m256i s2 = _mm256_setzero_si256(); // col 1b
377 __m256i s3 = _mm256_setzero_si256();
378 __m256i s4 = _mm256_setzero_si256(); // col 2a
379 __m256i s5 = _mm256_setzero_si256();
380 __m256i s6 = _mm256_setzero_si256(); // col 2b
381 __m256i s7 = _mm256_setzero_si256();
382 __m256i s8, s9, s10, s11;
383 __m256i s12, s13, s14, s15;
384
385 s11 = _mm256_set1_epi64x(0xFFFFFFFFUL);
386 const __m256i* x_ptr = (__m256i*)x;
387 const __m256i* y_ptr = (__m256i*)y;
388 __m256i* res_ptr = (__m256i*)res;
389 for (uint64_t i = 0; i < ell; ++i) {
390 // treat item a
391 s8 = _mm256_loadu_si256(x_ptr);
392 s9 = _mm256_srli_epi64(s8, 32);
393
394 s12 = _mm256_loadu_si256(y_ptr);
395 s13 = _mm256_loadu_si256(y_ptr + 2);
396 s14 = _mm256_srli_epi64(s12, 32);
397 s15 = _mm256_srli_epi64(s13, 32);
398
399 s12 = _mm256_mul_epu32(s8, s12); // c1a -> s0,s1
400 s13 = _mm256_mul_epu32(s8, s13); // c2a -> s4,s5
401 s14 = _mm256_mul_epu32(s9, s14); // c1a -> s0,s1
402 s15 = _mm256_mul_epu32(s9, s15); // c2a -> s4,s5
403
404 s8 = _mm256_and_si256(s12, s11); // -> s0
405 s9 = _mm256_and_si256(s13, s11); // -> s4
406 s12 = _mm256_srli_epi64(s12, 32); // -> s1
407 s13 = _mm256_srli_epi64(s13, 32); // -> s5
408 s0 = _mm256_add_epi64(s0, s8);
409 s1 = _mm256_add_epi64(s1, s12);
410 s4 = _mm256_add_epi64(s4, s9);
411 s5 = _mm256_add_epi64(s5, s13);
412
413 s8 = _mm256_and_si256(s14, s11); // -> s0
414 s9 = _mm256_and_si256(s15, s11); // -> s4
415 s14 = _mm256_srli_epi64(s14, 32); // -> s1
416 s15 = _mm256_srli_epi64(s15, 32); // -> s5
417 s0 = _mm256_add_epi64(s0, s8);
418 s1 = _mm256_add_epi64(s1, s14);
419 s4 = _mm256_add_epi64(s4, s9);
420 s5 = _mm256_add_epi64(s5, s15);
421
422 // treat item b
423 s8 = _mm256_loadu_si256(x_ptr + 1);
424 s9 = _mm256_srli_epi64(s8, 32);
425
426 s12 = _mm256_loadu_si256(y_ptr + 1);
427 s13 = _mm256_loadu_si256(y_ptr + 3);
428 s14 = _mm256_srli_epi64(s12, 32);
429 s15 = _mm256_srli_epi64(s13, 32);
430
431 s12 = _mm256_mul_epu32(s8, s12); // c1b -> s2,s3
432 s13 = _mm256_mul_epu32(s8, s13); // c2b -> s6,s7
433 s14 = _mm256_mul_epu32(s9, s14); // c1b -> s2,s3
434 s15 = _mm256_mul_epu32(s9, s15); // c2b -> s6,s7
435
436 s8 = _mm256_and_si256(s12, s11); // -> s2
437 s9 = _mm256_and_si256(s13, s11); // -> s6
438 s12 = _mm256_srli_epi64(s12, 32); // -> s3
439 s13 = _mm256_srli_epi64(s13, 32); // -> s7
440 s2 = _mm256_add_epi64(s2, s8);
441 s3 = _mm256_add_epi64(s3, s12);
442 s6 = _mm256_add_epi64(s6, s9);
443 s7 = _mm256_add_epi64(s7, s13);
444
445 s8 = _mm256_and_si256(s14, s11); // -> s2
446 s9 = _mm256_and_si256(s15, s11); // -> s6
447 s14 = _mm256_srli_epi64(s14, 32); // -> s3
448 s15 = _mm256_srli_epi64(s15, 32); // -> s7
449 s2 = _mm256_add_epi64(s2, s8);
450 s3 = _mm256_add_epi64(s3, s14);
451 s6 = _mm256_add_epi64(s6, s9);
452 s7 = _mm256_add_epi64(s7, s15);
453
454 x_ptr += 2;
455 y_ptr += 4;
456 }
457 // final reduction
458 const uint64_t H2 = precomp->h;
459 s8 = _mm256_set1_epi64x((UINT64_C(1) << H2) - 1); // MASK2
460 s9 = _mm256_loadu_si256((__m256i*)precomp->s2l_pow_red); // S2L_POW_RED
461 s10 = _mm256_loadu_si256((__m256i*)precomp->s2h_pow_red); // S2H_POW_RED
462 //--- s0,s1
463 s11 = _mm256_and_si256(s1, s8);
464 s12 = _mm256_srli_epi64(s1, H2);
465 s13 = _mm256_mul_epu32(s11, s9);
466 s14 = _mm256_mul_epu32(s12, s10);
467 s0 = _mm256_add_epi64(s0, s13);
468 s0 = _mm256_add_epi64(s0, s14);
469 _mm256_storeu_si256(res_ptr + 0, s0);
470 //--- s2,s3
471 s11 = _mm256_and_si256(s3, s8);
472 s12 = _mm256_srli_epi64(s3, H2);
473 s13 = _mm256_mul_epu32(s11, s9);
474 s14 = _mm256_mul_epu32(s12, s10);
475 s2 = _mm256_add_epi64(s2, s13);
476 s2 = _mm256_add_epi64(s2, s14);
477 _mm256_storeu_si256(res_ptr + 1, s2);
478 //--- s4,s5
479 s11 = _mm256_and_si256(s5, s8);
480 s12 = _mm256_srli_epi64(s5, H2);
481 s13 = _mm256_mul_epu32(s11, s9);
482 s14 = _mm256_mul_epu32(s12, s10);
483 s4 = _mm256_add_epi64(s4, s13);
484 s4 = _mm256_add_epi64(s4, s14);
485 _mm256_storeu_si256(res_ptr + 2, s4);
486 //--- s6,s7
487 s11 = _mm256_and_si256(s7, s8);
488 s12 = _mm256_srli_epi64(s7, H2);
489 s13 = _mm256_mul_epu32(s11, s9);
490 s14 = _mm256_mul_epu32(s12, s10);
491 s6 = _mm256_add_epi64(s6, s13);
492 s6 = _mm256_add_epi64(s6, s14);
493 _mm256_storeu_si256(res_ptr + 3, s6);
494}
495
496EXPORT void q120x2_vec_mat1col_product_bbc_avx2(q120_mat1col_product_bbc_precomp* precomp, const uint64_t ell,
497 q120b* const res, const q120b* const x, const q120c* const y) {
498 __m256i s0 = _mm256_setzero_si256(); // col 1a
499 __m256i s1 = _mm256_setzero_si256();
500 __m256i s2 = _mm256_setzero_si256(); // col 1b
501 __m256i s3 = _mm256_setzero_si256();
502 __m256i s4 = _mm256_set1_epi64x(0xFFFFFFFFUL);
503 __m256i s8, s9, s10, s11;
504 __m256i s12, s13, s14, s15;
505
506 const __m256i* x_ptr = (__m256i*)x;
507 const __m256i* y_ptr = (__m256i*)y;
508 __m256i* res_ptr = (__m256i*)res;
509 for (uint64_t i = 0; i < ell; ++i) {
510 s8 = _mm256_loadu_si256(x_ptr);
511 s9 = _mm256_loadu_si256(x_ptr + 1);
512 s10 = _mm256_srli_epi64(s8, 32);
513 s11 = _mm256_srli_epi64(s9, 32);
514
515 s12 = _mm256_loadu_si256(y_ptr);
516 s13 = _mm256_loadu_si256(y_ptr + 1);
517 s14 = _mm256_srli_epi64(s12, 32);
518 s15 = _mm256_srli_epi64(s13, 32);
519
520 s12 = _mm256_mul_epu32(s8, s12); // -> s0,s1
521 s13 = _mm256_mul_epu32(s9, s13); // -> s2,s3
522 s14 = _mm256_mul_epu32(s10, s14); // -> s0,s1
523 s15 = _mm256_mul_epu32(s11, s15); // -> s2,s3
524
525 s8 = _mm256_and_si256(s12, s4); // -> s0
526 s9 = _mm256_and_si256(s13, s4); // -> s2
527 s10 = _mm256_and_si256(s14, s4); // -> s0
528 s11 = _mm256_and_si256(s15, s4); // -> s2
529 s12 = _mm256_srli_epi64(s12, 32); // -> s1
530 s13 = _mm256_srli_epi64(s13, 32); // -> s3
531 s14 = _mm256_srli_epi64(s14, 32); // -> s1
532 s15 = _mm256_srli_epi64(s15, 32); // -> s3
533
534 s0 = _mm256_add_epi64(s0, s8);
535 s1 = _mm256_add_epi64(s1, s12);
536 s2 = _mm256_add_epi64(s2, s9);
537 s3 = _mm256_add_epi64(s3, s13);
538 s0 = _mm256_add_epi64(s0, s10);
539 s1 = _mm256_add_epi64(s1, s14);
540 s2 = _mm256_add_epi64(s2, s11);
541 s3 = _mm256_add_epi64(s3, s15);
542
543 x_ptr += 2;
544 y_ptr += 2;
545 }
546 // final reduction
547 const uint64_t H2 = precomp->h;
548 s8 = _mm256_set1_epi64x((UINT64_C(1) << H2) - 1); // MASK2
549 s9 = _mm256_loadu_si256((__m256i*)precomp->s2l_pow_red); // S2L_POW_RED
550 s10 = _mm256_loadu_si256((__m256i*)precomp->s2h_pow_red); // S2H_POW_RED
551 //--- s0,s1
552 s11 = _mm256_and_si256(s1, s8);
553 s12 = _mm256_srli_epi64(s1, H2);
554 s13 = _mm256_mul_epu32(s11, s9);
555 s14 = _mm256_mul_epu32(s12, s10);
556 s0 = _mm256_add_epi64(s0, s13);
557 s0 = _mm256_add_epi64(s0, s14);
558 _mm256_storeu_si256(res_ptr + 0, s0);
559 //--- s2,s3
560 s11 = _mm256_and_si256(s3, s8);
561 s12 = _mm256_srli_epi64(s3, H2);
562 s13 = _mm256_mul_epu32(s11, s9);
563 s14 = _mm256_mul_epu32(s12, s10);
564 s2 = _mm256_add_epi64(s2, s13);
565 s2 = _mm256_add_epi64(s2, s14);
566 _mm256_storeu_si256(res_ptr + 1, s2);
567}