Line data Source code
1 : /* ntt_gfp.c - low-level radix-2 dif/dit ntt routines over GF(p)
2 :
3 : Copyright 2005, 2006, 2007, 2008, 2009 Dave Newman, Jason Papadopoulos,
4 : Brian Gladman, Alexander Kruppa, Paul Zimmermann.
5 :
6 : The SP Library is free software; you can redistribute it and/or modify
7 : it under the terms of the GNU Lesser General Public License as published by
8 : the Free Software Foundation; either version 3 of the License, or (at your
9 : option) any later version.
10 :
11 : The SP Library is distributed in the hope that it will be useful, but
12 : WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
13 : or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
14 : License for more details.
15 :
16 : You should have received a copy of the GNU Lesser General Public License
17 : along with the SP Library; see the file COPYING.LIB. If not, write to
18 : the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
19 : MA 02110-1301, USA. */
20 :
21 : #include "sp.h"
22 : #include "ecm-impl.h"
23 :
24 : /*--------------------------- FORWARD NTT --------------------------------*/
25 223078692 : static void bfly_dif(spv_t x0, spv_t x1, spv_t w,
26 : spv_size_t len, sp_t p, sp_t d)
27 : {
28 223078692 : spv_size_t i = 0;
29 :
30 : #if (defined(__GNUC__) || defined(__ICL)) && \
31 : defined(__i386__) && defined(HAVE_SSE2)
32 :
33 : asm volatile (
34 : "movd %6, %%xmm6 \n\t"
35 : "pshufd $0x44, %%xmm6, %%xmm5 \n\t"
36 : "pshufd $0, %%xmm6, %%xmm6 \n\t"
37 : "movd %7, %%xmm7 \n\t"
38 : "pshufd $0, %%xmm7, %%xmm7 \n\t"
39 :
40 : "0: \n\t"
41 : "movdqa (%1,%4,4), %%xmm0 \n\t"
42 : "movdqa (%2,%4,4), %%xmm1 \n\t"
43 : "movdqa %%xmm1, %%xmm2 \n\t"
44 : "paddd %%xmm0, %%xmm1 \n\t"
45 : "psubd %%xmm2, %%xmm0 \n\t"
46 : "psubd %%xmm6, %%xmm1 \n\t"
47 :
48 : "pxor %%xmm2, %%xmm2 \n\t"
49 : "pcmpgtd %%xmm1, %%xmm2 \n\t"
50 : "pand %%xmm6, %%xmm2 \n\t"
51 : "paddd %%xmm2, %%xmm1 \n\t"
52 : "movdqa %%xmm1, (%1,%4,4) \n\t"
53 :
54 : "pxor %%xmm2, %%xmm2 \n\t"
55 : "pcmpgtd %%xmm0, %%xmm2 \n\t"
56 : "pand %%xmm6, %%xmm2 \n\t"
57 : "paddd %%xmm2, %%xmm0 \n\t"
58 :
59 : "movdqa (%3,%4,4), %%xmm2 \n\t"
60 : "addl $4, %4 \n\t" /* INC */
61 : "pshufd $0x31, %%xmm0, %%xmm1\n\t"
62 : "pshufd $0x31, %%xmm2, %%xmm3\n\t"
63 : "pmuludq %%xmm2, %%xmm0 \n\t"
64 : "pmuludq %%xmm3, %%xmm1 \n\t"
65 :
66 : "movdqa %%xmm0, %%xmm2 \n\t"
67 : "movdqa %%xmm1, %%xmm3 \n\t"
68 : "psrlq $" STRING((2*SP_NUMB_BITS - W_TYPE_SIZE)) ", %%xmm2 \n\t"
69 : "pmuludq %%xmm7, %%xmm2 \n\t"
70 : "psrlq $" STRING((2*SP_NUMB_BITS - W_TYPE_SIZE)) ", %%xmm3 \n\t"
71 : "pmuludq %%xmm7, %%xmm3 \n\t"
72 :
73 : #if SP_NUMB_BITS < W_TYPE_SIZE - 1
74 : "psrlq $33, %%xmm2 \n\t"
75 : "pmuludq %%xmm6, %%xmm2 \n\t"
76 : "psrlq $33, %%xmm3 \n\t"
77 : "pmuludq %%xmm6, %%xmm3 \n\t"
78 : "psubq %%xmm2, %%xmm0 \n\t"
79 : "psubq %%xmm3, %%xmm1 \n\t"
80 : #else
81 : "pshufd $0xf5, %%xmm2, %%xmm2 \n\t"
82 : "pmuludq %%xmm6, %%xmm2 \n\t"
83 : "pshufd $0xf5, %%xmm3, %%xmm3 \n\t"
84 : "pmuludq %%xmm6, %%xmm3 \n\t"
85 : "psubq %%xmm2, %%xmm0 \n\t"
86 : "psubq %%xmm3, %%xmm1 \n\t"
87 :
88 : "psubq %%xmm5, %%xmm0 \n\t"
89 : "psubq %%xmm5, %%xmm1 \n\t"
90 : "pshufd $0xf5, %%xmm0, %%xmm2 \n\t"
91 : "pshufd $0xf5, %%xmm1, %%xmm3 \n\t"
92 : "pand %%xmm5, %%xmm2 \n\t"
93 : "pand %%xmm5, %%xmm3 \n\t"
94 : "paddq %%xmm2, %%xmm0 \n\t"
95 : "paddq %%xmm3, %%xmm1 \n\t"
96 : #endif
97 : "pshufd $0x8, %%xmm0, %%xmm0 \n\t"
98 : "pshufd $0x8, %%xmm1, %%xmm1 \n\t"
99 : "punpckldq %%xmm1, %%xmm0 \n\t"
100 : "psubd %%xmm6, %%xmm0 \n\t"
101 : "pxor %%xmm1, %%xmm1 \n\t"
102 : "pcmpgtd %%xmm0, %%xmm1 \n\t"
103 : "pand %%xmm6, %%xmm1 \n\t"
104 : "paddd %%xmm1, %%xmm0 \n\t"
105 : "movdqa %%xmm0, -16(%2,%4,4) \n\t"
106 :
107 : "cmpl %5, %4 \n\t"
108 : "jne 0b \n\t"
109 :
110 : :"=r"(i)
111 : :"r"(x0), "r"(x1), "r"(w), "0"(i), "g"(len), "g"(p), "g"(d)
112 : :"%xmm0", "%xmm1", "%xmm2", "%xmm3",
113 : "%xmm5", "%xmm6", "%xmm7", "cc", "memory");
114 : #elif defined( _MSC_VER ) && defined( SSE2)
115 : __asm
116 : { push esi
117 : push edi
118 : mov edi, x0
119 : mov esi, x1
120 : mov edx, w
121 : xor ecx, ecx
122 : mov eax, len
123 : movd xmm6, p
124 : pshufd xmm5, xmm6, 0x44
125 : pshufd xmm6, xmm6, 0
126 : movd xmm7, d
127 : pshufd xmm7, xmm7, 0
128 :
129 : L0: movdqa xmm0, [edi+ecx*4]
130 : movdqa xmm1, [esi+ecx*4]
131 : movdqa xmm2, xmm1
132 : paddd xmm1, xmm0
133 : psubd xmm0, xmm2
134 : psubd xmm1, xmm6
135 :
136 : pxor xmm2, xmm2
137 : pcmpgtd xmm2, xmm1
138 : pand xmm2, xmm6
139 : paddd xmm1, xmm2
140 : movdqa [edi+ecx*4], xmm1
141 :
142 : pxor xmm2, xmm2
143 : pcmpgtd xmm2, xmm0
144 : pand xmm2, xmm6
145 : paddd xmm0, xmm2
146 :
147 : movdqa xmm2, [edx+ecx*4]
148 : add ecx, 4
149 : pshufd xmm1, xmm0, 0x31
150 : pshufd xmm3, xmm2, 0x31
151 : pmuludq xmm0, xmm2
152 : pmuludq xmm1, xmm3
153 :
154 : movdqa xmm2, xmm0
155 : movdqa xmm3, xmm1
156 : psrlq xmm2, 2*SP_NUMB_BITS - W_TYPE_SIZE
157 : pmuludq xmm2, xmm7
158 : psrlq xmm3, 2*SP_NUMB_BITS - W_TYPE_SIZE
159 : pmuludq xmm3, xmm7
160 :
161 : #if SP_NUMB_BITS < W_TYPE_SIZE - 1
162 : psrlq xmm2, 33
163 : pmuludq xmm2, xmm6
164 : psrlq xmm3, 33
165 : pmuludq xmm3, xmm6
166 : psubq xmm0, xmm2
167 : psubq xmm1, xmm3
168 : #else
169 : pshufd xmm2, xmm2, 0xf5
170 : pmuludq xmm2, xmm6
171 : pshufd xmm3, xmm3, 0xf5
172 : pmuludq xmm3, xmm6
173 : psubq xmm0, xmm2
174 : psubq xmm1, xmm3
175 :
176 : psubq xmm0, xmm5
177 : psubq xmm1, xmm5
178 : pshufd xmm2, xmm0, 0xf5
179 : pshufd xmm3, xmm1, 0xf5
180 : pand xmm2, xmm5
181 : pand xmm3, xmm5
182 : paddq xmm0, xmm2
183 : paddq xmm1, xmm3
184 : #endif
185 : pshufd xmm0, xmm0, 0x8
186 : pshufd xmm1, xmm1, 0x8
187 : punpckldq xmm0, xmm1
188 : psubd xmm0, xmm6
189 : pxor xmm1, xmm1
190 : pcmpgtd xmm1, xmm0
191 : pand xmm1, xmm6
192 : paddd xmm0, xmm1
193 : movdqa [esi+ecx*4-16], xmm0
194 : cmp eax, ecx
195 : jne L0
196 : pop edi
197 : pop esi
198 : }
199 : #else
200 9611137932 : for (i = 0; i < len; i++)
201 : {
202 9388059240 : sp_t w0 = w[i];
203 9388059240 : sp_t t0 = x0[i];
204 9388059240 : sp_t t1 = x1[i];
205 : sp_t t2, t3;
206 9388059240 : t2 = sp_add (t0, t1, p);
207 9388059240 : t3 = sp_sub (t0, t1, p);
208 9388059240 : t3 = sp_mul (t3, w0, p, d);
209 9388059240 : x0[i] = t2;
210 9388059240 : x1[i] = t3;
211 : }
212 : #endif
213 223078692 : }
214 :
215 : static void
216 409907895 : spv_ntt_dif_core (spv_t x, spv_t w,
217 : spv_size_t log2_len, sp_t p, sp_t d)
218 : {
219 : spv_size_t len;
220 : spv_t x0, x1;
221 :
222 : /* handle small transforms immediately */
223 409907895 : switch (log2_len) {
224 0 : case 0:
225 0 : return;
226 270 : case 1:
227 : {
228 270 : sp_t t0 = x[0];
229 270 : sp_t t1 = x[1];
230 270 : x[0] = sp_add (t0, t1, p);
231 270 : x[1] = sp_sub (t0, t1, p);
232 270 : return;
233 : }
234 207390 : case 2:
235 : {
236 207390 : sp_t t0 = x[0];
237 207390 : sp_t t1 = x[1];
238 207390 : sp_t t2 = x[2];
239 207390 : sp_t t3 = x[3];
240 : sp_t t4, t5, t6, t7;
241 207390 : t4 = sp_add (t0, t2, p);
242 207390 : t6 = sp_sub (t0, t2, p);
243 207390 : t5 = sp_add (t1, t3, p);
244 207390 : t7 = sp_sub (t1, t3, p);
245 207390 : x[0] = sp_add (t4, t5, p);
246 207390 : x[1] = sp_sub (t4, t5, p);
247 207390 : t7 = sp_mul (t7, w[1], p, d);
248 207390 : x[2] = sp_add (t6, t7, p);
249 207390 : x[3] = sp_sub (t6, t7, p);
250 207390 : return;
251 : }
252 206186375 : case 3:
253 : {
254 206186375 : sp_t t0 = x[0];
255 206186375 : sp_t t1 = x[1];
256 206186375 : sp_t t2 = x[2];
257 206186375 : sp_t t3 = x[3];
258 206186375 : sp_t t4 = x[4];
259 206186375 : sp_t t5 = x[5];
260 206186375 : sp_t t6 = x[6];
261 206186375 : sp_t t7 = x[7];
262 : sp_t t8, t9, t10, t11, t12, t13, t14, t15;
263 :
264 206186375 : t8 = sp_add (t0, t4, p);
265 206186375 : t12 = sp_sub (t0, t4, p);
266 206186375 : t9 = sp_add (t1, t5, p);
267 206186375 : t13 = sp_sub (t1, t5, p);
268 206186375 : t13 = sp_mul (t13, w[1], p, d);
269 206186375 : t10 = sp_add (t2, t6, p);
270 206186375 : t14 = sp_sub (t2, t6, p);
271 206186375 : t14 = sp_mul (t14, w[2], p, d);
272 206186375 : t11 = sp_add (t3, t7, p);
273 206186375 : t15 = sp_sub (t3, t7, p);
274 206186375 : t15 = sp_mul (t15, w[3], p, d);
275 :
276 206186375 : t0 = sp_add (t8, t10, p);
277 206186375 : t2 = sp_sub (t8, t10, p);
278 206186375 : t1 = sp_add (t9, t11, p);
279 206186375 : t3 = sp_sub (t9, t11, p);
280 206186375 : t3 = sp_mul (t3, w[2], p, d);
281 206186375 : x[0] = sp_add (t0, t1, p);
282 206186375 : x[1] = sp_sub (t0, t1, p);
283 206186375 : x[2] = sp_add (t2, t3, p);
284 206186375 : x[3] = sp_sub (t2, t3, p);
285 :
286 206186375 : t0 = sp_add (t12, t14, p);
287 206186375 : t2 = sp_sub (t12, t14, p);
288 206186375 : t1 = sp_add (t13, t15, p);
289 206186375 : t3 = sp_sub (t13, t15, p);
290 206186375 : t3 = sp_mul (t3, w[2], p, d);
291 206186375 : x[4] = sp_add (t0, t1, p);
292 206186375 : x[5] = sp_sub (t0, t1, p);
293 206186375 : x[6] = sp_add (t2, t3, p);
294 206186375 : x[7] = sp_sub (t2, t3, p);
295 206186375 : return;
296 : }
297 : }
298 :
299 203513860 : len = 1 << (log2_len - 1);
300 203513860 : x0 = x;
301 203513860 : x1 = x + len;
302 203513860 : bfly_dif (x0, x1, w, len, p, d);
303 203513860 : spv_ntt_dif_core (x0, w + len, log2_len - 1, p, d);
304 203513860 : spv_ntt_dif_core (x1, w + len, log2_len - 1, p, d);
305 : }
306 :
307 : void
308 3551797 : spv_ntt_gfp_dif (spv_t x, spv_size_t log2_len, spm_t data)
309 : {
310 3551797 : sp_t p = data->sp;
311 3551797 : sp_t d = data->mul_c;
312 :
313 3551797 : if (log2_len <= NTT_GFP_TWIDDLE_DIF_BREAKOVER)
314 : {
315 2880175 : spv_t w = data->nttdata->twiddle +
316 2880175 : data->nttdata->twiddle_size - (1 << log2_len);
317 2880175 : spv_ntt_dif_core (x, w, log2_len, p, d);
318 : }
319 : else
320 : {
321 : /* recursive version for data that
322 : doesn't fit in the L1 cache */
323 671622 : spv_size_t len = 1 << (log2_len - 1);
324 671622 : spv_t x0 = x;
325 671622 : spv_t x1 = x + len;
326 671622 : spv_t roots = data->nttdata->ntt_roots;
327 :
328 : {
329 : spv_size_t i;
330 671622 : spv_size_t block_size = MIN(len, MAX_NTT_BLOCK_SIZE);
331 671622 : sp_t root = roots[log2_len];
332 671622 : spv_t w = data->scratch;
333 :
334 671622 : w[0] = 1;
335 47432448 : for (i = 1; i < block_size; i++)
336 46760826 : w[i] = sp_mul (w[i-1], root, p, d);
337 :
338 671622 : root = sp_pow (root, block_size, p, d);
339 :
340 20236454 : for (i = 0; i < len; i += block_size)
341 : {
342 19564832 : if (i)
343 18893210 : spv_mul_sp (w, w, root, block_size, p, d);
344 :
345 19564832 : bfly_dif (x0 + i, x1 + i, w, block_size, p, d);
346 : }
347 : }
348 :
349 671622 : spv_ntt_gfp_dif (x0, log2_len - 1, data);
350 671622 : spv_ntt_gfp_dif (x1, log2_len - 1, data);
351 : }
352 3551797 : }
353 :
354 : /*--------------------------- INVERSE NTT --------------------------------*/
355 118942961 : static inline void bfly_dit(spv_t x0, spv_t x1, spv_t w,
356 : spv_size_t len, sp_t p, sp_t d)
357 : {
358 118942961 : spv_size_t i = 0;
359 :
360 : #if (defined(__GNUC__) || defined(__ICL)) && \
361 : defined(__i386__) && defined(HAVE_SSE2)
362 :
363 : asm volatile (
364 : "movd %6, %%xmm6 \n\t"
365 : "pshufd $0x44, %%xmm6, %%xmm5 \n\t"
366 : "pshufd $0, %%xmm6, %%xmm6 \n\t"
367 : "movd %7, %%xmm7 \n\t"
368 : "pshufd $0, %%xmm7, %%xmm7 \n\t"
369 :
370 : "0: \n\t"
371 : "movdqa (%2,%4,4), %%xmm0 \n\t"
372 : "movdqa (%3,%4,4), %%xmm2 \n\t"
373 : "pshufd $0x31, %%xmm0, %%xmm1\n\t"
374 : "pshufd $0x31, %%xmm2, %%xmm3\n\t"
375 : "pmuludq %%xmm2, %%xmm0 \n\t"
376 : "pmuludq %%xmm3, %%xmm1 \n\t"
377 :
378 : "movdqa %%xmm0, %%xmm2 \n\t"
379 : "movdqa %%xmm1, %%xmm3 \n\t"
380 : "psrlq $" STRING((2*SP_NUMB_BITS - W_TYPE_SIZE)) ", %%xmm2 \n\t"
381 : "pmuludq %%xmm7, %%xmm2 \n\t"
382 : "psrlq $" STRING((2*SP_NUMB_BITS - W_TYPE_SIZE)) ", %%xmm3 \n\t"
383 : "pmuludq %%xmm7, %%xmm3 \n\t"
384 :
385 : #if SP_NUMB_BITS < W_TYPE_SIZE - 1
386 : "psrlq $33, %%xmm2 \n\t"
387 : "pmuludq %%xmm6, %%xmm2 \n\t"
388 : "psrlq $33, %%xmm3 \n\t"
389 : "pmuludq %%xmm6, %%xmm3 \n\t"
390 : "psubq %%xmm2, %%xmm0 \n\t"
391 : "psubq %%xmm3, %%xmm1 \n\t"
392 : #else
393 : "pshufd $0xf5, %%xmm2, %%xmm2 \n\t"
394 : "pmuludq %%xmm6, %%xmm2 \n\t"
395 : "pshufd $0xf5, %%xmm3, %%xmm3 \n\t"
396 : "pmuludq %%xmm6, %%xmm3 \n\t"
397 : "psubq %%xmm2, %%xmm0 \n\t"
398 : "psubq %%xmm3, %%xmm1 \n\t"
399 :
400 : "psubq %%xmm5, %%xmm0 \n\t"
401 : "psubq %%xmm5, %%xmm1 \n\t"
402 : "pshufd $0xf5, %%xmm0, %%xmm2 \n\t"
403 : "pshufd $0xf5, %%xmm1, %%xmm3 \n\t"
404 : "pand %%xmm5, %%xmm2 \n\t"
405 : "pand %%xmm5, %%xmm3 \n\t"
406 : "paddq %%xmm2, %%xmm0 \n\t"
407 : "paddq %%xmm3, %%xmm1 \n\t"
408 : #endif
409 : "pshufd $0x8, %%xmm0, %%xmm0 \n\t"
410 : "pshufd $0x8, %%xmm1, %%xmm1 \n\t"
411 : "punpckldq %%xmm1, %%xmm0 \n\t"
412 : "psubd %%xmm6, %%xmm0 \n\t"
413 : "pxor %%xmm1, %%xmm1 \n\t"
414 : "pcmpgtd %%xmm0, %%xmm1 \n\t"
415 : "pand %%xmm6, %%xmm1 \n\t"
416 : "paddd %%xmm0, %%xmm1 \n\t"
417 :
418 : "movdqa (%1,%4,4), %%xmm0 \n\t"
419 : "movdqa %%xmm1, %%xmm2 \n\t"
420 : "paddd %%xmm0, %%xmm1 \n\t"
421 : "psubd %%xmm2, %%xmm0 \n\t"
422 : "psubd %%xmm6, %%xmm1 \n\t"
423 :
424 : "pxor %%xmm2, %%xmm2 \n\t"
425 : "pcmpgtd %%xmm1, %%xmm2 \n\t"
426 : "pand %%xmm6, %%xmm2 \n\t"
427 : "paddd %%xmm2, %%xmm1 \n\t"
428 : "movdqa %%xmm1, (%1,%4,4) \n\t"
429 :
430 : "pxor %%xmm2, %%xmm2 \n\t"
431 : "pcmpgtd %%xmm0, %%xmm2 \n\t"
432 : "pand %%xmm6, %%xmm2 \n\t"
433 : "paddd %%xmm2, %%xmm0 \n\t"
434 : "movdqa %%xmm0, (%2,%4,4) \n\t"
435 :
436 : "addl $4, %4 \n\t" /* INC */
437 : "cmpl %5, %4 \n\t"
438 : "jne 0b \n\t"
439 :
440 : :"=r"(i)
441 : :"r"(x0), "r"(x1), "r"(w), "0"(i), "g"(len), "g"(p), "g"(d)
442 : :"%xmm0", "%xmm1", "%xmm2", "%xmm3",
443 : "%xmm5", "%xmm6", "%xmm7", "cc", "memory");
444 : #elif defined( _MSC_VER ) && defined( SSE2)
445 : __asm
446 : { push esi
447 : push edi
448 : mov edi, x0
449 : mov esi, x1
450 : mov edx, w
451 : xor ecx, ecx
452 : mov eax, len
453 : movd xmm6, p
454 : pshufd xmm5, xmm6, 0x44
455 : pshufd xmm6, xmm6, 0
456 : movd xmm7, d
457 : pshufd xmm7, xmm7, 0
458 :
459 : L0: movdqa xmm0, [esi+ecx*4]
460 : movdqa xmm2, [edx+ecx*4]
461 : pshufd xmm1, xmm0, 0x31
462 : pshufd xmm3, xmm2, 0x31
463 : pmuludq xmm0, xmm2
464 : pmuludq xmm1, xmm3
465 :
466 : movdqa xmm2, xmm0
467 : movdqa xmm3, xmm1
468 : psrlq xmm2, 2*SP_NUMB_BITS - W_TYPE_SIZE
469 : pmuludq xmm2, xmm7
470 : psrlq xmm3, 2*SP_NUMB_BITS - W_TYPE_SIZE
471 : pmuludq xmm3, xmm7
472 :
473 : #if SP_NUMB_BITS < W_TYPE_SIZE - 1
474 : psrlq xmm2, 33
475 : pmuludq xmm2, xmm6
476 : psrlq xmm3, 33
477 : pmuludq xmm3, xmm6
478 : psubq xmm0, xmm2
479 : psubq xmm1, xmm3
480 : #else
481 : pshufd xmm2, xmm2, 0xf5
482 : pmuludq xmm2, xmm6
483 : pshufd xmm3, xmm3, 0xf5
484 : pmuludq xmm3, xmm6
485 : psubq xmm0, xmm2
486 : psubq xmm1, xmm3
487 :
488 : psubq xmm0, xmm5
489 : psubq xmm1, xmm5
490 : pshufd xmm2, xmm0, 0xf5
491 : pshufd xmm3, xmm1, 0xf5
492 : pand xmm2, xmm5
493 : pand xmm3, xmm5
494 : paddq xmm0, xmm2
495 : paddq xmm1, xmm3
496 : #endif
497 : pshufd xmm0, xmm0, 0x8
498 : pshufd xmm1, xmm1, 0x8
499 : punpckldq xmm0, xmm1
500 : psubd xmm0, xmm6
501 : pxor xmm1, xmm1
502 : pcmpgtd xmm1, xmm0
503 : pand xmm1, xmm6
504 : paddd xmm1, xmm0
505 :
506 : movdqa xmm0, [edi+ecx*4]
507 : movdqa xmm2, xmm1
508 : paddd xmm1, xmm0
509 : psubd xmm0, xmm2
510 : psubd xmm1, xmm6
511 :
512 : pxor xmm2, xmm2
513 : pcmpgtd xmm2, xmm1
514 : pand xmm2, xmm6
515 : paddd xmm1, xmm2
516 : movdqa [edi+ecx*4], xmm1
517 :
518 : pxor xmm2, xmm2
519 : pcmpgtd xmm2, xmm0
520 : pand xmm2, xmm6
521 : paddd xmm0, xmm2
522 : movdqa [esi+ecx*4], xmm0
523 : add ecx, 4
524 : cmp eax, ecx
525 : jne L0
526 : pop edi
527 : pop esi
528 : }
529 : #else
530 5325555969 : for (i = 0; i < len; i++)
531 : {
532 5206613008 : sp_t w0 = w[i];
533 5206613008 : sp_t t0 = x0[i];
534 5206613008 : sp_t t1 = x1[i];
535 5206613008 : t1 = sp_mul (t1, w0, p, d);
536 5206613008 : x0[i] = sp_add (t0, t1, p);
537 5206613008 : x1[i] = sp_sub (t0, t1, p);
538 : }
539 : #endif
540 118942961 : }
541 :
542 : static void
543 233715685 : spv_ntt_dit_core (spv_t x, spv_t w,
544 : spv_size_t log2_len, sp_t p, sp_t d)
545 : {
546 : spv_size_t len;
547 : spv_t x0, x1;
548 :
549 : /* handle small transforms immediately */
550 233715685 : switch (log2_len) {
551 0 : case 0:
552 0 : return;
553 0 : case 1:
554 : {
555 0 : sp_t t0 = x[0];
556 0 : sp_t t1 = x[1];
557 0 : x[0] = sp_add (t0, t1, p);
558 0 : x[1] = sp_sub (t0, t1, p);
559 0 : return;
560 : }
561 206732 : case 2:
562 : {
563 206732 : sp_t t0 = x[0];
564 206732 : sp_t t1 = x[1];
565 206732 : sp_t t2 = x[2];
566 206732 : sp_t t3 = x[3];
567 : sp_t t4, t5, t6, t7;
568 206732 : t4 = sp_add (t0, t1, p);
569 206732 : t5 = sp_sub (t0, t1, p);
570 206732 : t6 = sp_add (t2, t3, p);
571 206732 : t7 = sp_sub (t2, t3, p);
572 206732 : x[0] = sp_add (t4, t6, p);
573 206732 : x[2] = sp_sub (t4, t6, p);
574 206732 : t7 = sp_mul (t7, w[1], p, d);
575 206732 : x[1] = sp_add (t5, t7, p);
576 206732 : x[3] = sp_sub (t5, t7, p);
577 206732 : return;
578 : }
579 117617640 : case 3:
580 : {
581 117617640 : sp_t t0 = x[0];
582 117617640 : sp_t t1 = x[1];
583 117617640 : sp_t t2 = x[2];
584 117617640 : sp_t t3 = x[3];
585 117617640 : sp_t t4 = x[4];
586 117617640 : sp_t t5 = x[5];
587 117617640 : sp_t t6 = x[6];
588 117617640 : sp_t t7 = x[7];
589 : sp_t t8, t9, t10, t11;
590 :
591 117617640 : t8 = sp_add(t0, t1, p);
592 117617640 : t9 = sp_sub(t0, t1, p);
593 117617640 : t10 = sp_add(t2, t3, p);
594 117617640 : t11 = sp_sub(t2, t3, p);
595 117617640 : t0 = sp_add(t8, t10, p);
596 117617640 : t2 = sp_sub(t8, t10, p);
597 117617640 : t11 = sp_mul (t11, w[2], p, d);
598 117617640 : t1 = sp_add(t9, t11, p);
599 117617640 : t3 = sp_sub(t9, t11, p);
600 :
601 117617640 : t8 = sp_add(t4, t5, p);
602 117617640 : t9 = sp_sub(t4, t5, p);
603 117617640 : t10 = sp_add(t6, t7, p);
604 117617640 : t11 = sp_sub(t6, t7, p);
605 117617640 : t4 = sp_add(t8, t10, p);
606 117617640 : t6 = sp_sub(t8, t10, p);
607 117617640 : t11 = sp_mul (t11, w[2], p, d);
608 117617640 : t5 = sp_add(t9, t11, p);
609 117617640 : t7 = sp_sub(t9, t11, p);
610 :
611 117617640 : x[0] = sp_add(t0, t4, p);
612 117617640 : x[4] = sp_sub(t0, t4, p);
613 117617640 : t5 = sp_mul (t5, w[1], p, d);
614 117617640 : x[1] = sp_add(t1, t5, p);
615 117617640 : x[5] = sp_sub(t1, t5, p);
616 117617640 : t6 = sp_mul (t6, w[2], p, d);
617 117617640 : x[2] = sp_add(t2, t6, p);
618 117617640 : x[6] = sp_sub(t2, t6, p);
619 117617640 : t7 = sp_mul (t7, w[3], p, d);
620 117617640 : x[3] = sp_add(t3, t7, p);
621 117617640 : x[7] = sp_sub(t3, t7, p);
622 117617640 : return;
623 : }
624 : }
625 :
626 115891313 : len = 1 << (log2_len - 1);
627 115891313 : x0 = x;
628 115891313 : x1 = x + len;
629 115891313 : spv_ntt_dit_core (x0, w + len, log2_len - 1, p, d);
630 115891313 : spv_ntt_dit_core (x1, w + len, log2_len - 1, p, d);
631 115891313 : bfly_dit (x0, x1, w, len, p, d);
632 : }
633 :
634 : void
635 2402287 : spv_ntt_gfp_dit (spv_t x, spv_size_t log2_len, spm_t data)
636 : {
637 2402287 : sp_t p = data->sp;
638 2402287 : sp_t d = data->mul_c;
639 :
640 2402287 : if (log2_len <= NTT_GFP_TWIDDLE_DIT_BREAKOVER)
641 : {
642 1933059 : spv_t w = data->inttdata->twiddle +
643 1933059 : data->inttdata->twiddle_size - (1 << log2_len);
644 1933059 : spv_ntt_dit_core (x, w, log2_len, p, d);
645 : }
646 : else
647 : {
648 469228 : spv_size_t len = 1 << (log2_len - 1);
649 469228 : spv_t x0 = x;
650 469228 : spv_t x1 = x + len;
651 469228 : spv_t roots = data->inttdata->ntt_roots;
652 :
653 469228 : spv_ntt_gfp_dit (x0, log2_len - 1, data);
654 469228 : spv_ntt_gfp_dit (x1, log2_len - 1, data);
655 :
656 : {
657 : spv_size_t i;
658 469228 : spv_size_t block_size = MIN(len, MAX_NTT_BLOCK_SIZE);
659 469228 : sp_t root = roots[log2_len];
660 469228 : spv_t w = data->scratch;
661 :
662 469228 : w[0] = 1;
663 21526016 : for (i = 1; i < block_size; i++)
664 21056788 : w[i] = sp_mul (w[i-1], root, p, d);
665 :
666 469228 : root = sp_pow (root, block_size, p, d);
667 :
668 3520876 : for (i = 0; i < len; i += block_size)
669 : {
670 3051648 : if (i)
671 2582420 : spv_mul_sp (w, w, root, block_size, p, d);
672 :
673 3051648 : bfly_dit (x0 + i, x1 + i, w, block_size, p, d);
674 : }
675 : }
676 : }
677 2402287 : }
|
