Line data Source code
1 : /* spv.c - "small prime vector" functions for arithmetic on vectors of
2 : residues modulo a single small prime
3 :
4 : Copyright 2005, 2006, 2007, 2008, 2009 Dave Newman, Jason Papadopoulos,
5 : Brian Gladman, Alexander Kruppa, Paul Zimmermann.
6 :
7 : The SP Library is free software; you can redistribute it and/or modify
8 : it under the terms of the GNU Lesser General Public License as published by
9 : the Free Software Foundation; either version 3 of the License, or (at your
10 : option) any later version.
11 :
12 : The SP Library is distributed in the hope that it will be useful, but
13 : WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14 : or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15 : License for more details.
16 :
17 : You should have received a copy of the GNU Lesser General Public License
18 : along with the SP Library; see the file COPYING.LIB. If not, write to
19 : the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
20 : MA 02110-1301, USA. */
21 :
22 : #include "config.h"
23 : #include <string.h> /* for memset */
24 : #include "sp.h"
25 : #ifdef USE_VALGRIND
26 : #include <valgrind/memcheck.h>
27 : #endif
28 : #include "ecm-impl.h"
29 :
30 : /* Routines for vectors of integers modulo r common small prime
31 : *
32 : * These are low-overhead routines that don't do memory allocation,
33 : * other than for temporary variables. Unless otherwise specified, any
34 : * of the input pointers can be equal. */
35 :
36 : #ifdef WANT_ASSERT
37 : int
38 : spv_verify (spv_t x, spv_size_t len, const sp_t m)
39 : {
40 : spv_size_t i;
41 : #ifdef USE_VALGRIND
42 : if (len > 0)
43 : {
44 : if (VALGRIND_CHECK_MEM_IS_DEFINED(x, len * sizeof(sp_t)) != 0)
45 : return 0;
46 : }
47 : #endif
48 : for (i = 0; i < len; i++)
49 : {
50 : if (x[i] >= m)
51 : return 0;
52 : }
53 :
54 : return 1;
55 : }
56 : #endif
57 :
58 : /* r = x */
59 : void
60 5895776 : spv_set (spv_t r, spv_t x, spv_size_t len)
61 : {
62 : #ifdef HAVE_MEMMOVE
63 : /* memmove doesn't rely on the assertion below */
64 5895776 : memmove (r, x, len * sizeof (sp_t));
65 : #else
66 : spv_size_t i;
67 :
68 : ASSERT (r >= x + len || x >= r);
69 :
70 : for (i = 0; i < len; i++)
71 : r[i] = x[i];
72 : #endif
73 5895776 : }
74 :
75 : /* r[0 ... len - 1] = x[len - 1 ... 0] */
76 : void
77 7695 : spv_rev (spv_t r, spv_t x, spv_size_t len)
78 : {
79 : spv_size_t i;
80 :
81 : ASSERT (r >= x + len || x >= r + len);
82 :
83 23749383 : for (i = 0; i < len; i++)
84 23741688 : r[i] = x[len - 1 - i];
85 7695 : }
86 :
87 : /* r = [y, y, ... ] */
88 : void
89 549593 : spv_set_sp (spv_t r, sp_t y, spv_size_t len)
90 : {
91 : spv_size_t i;
92 :
93 94867710 : for (i = 0; i < len; i++)
94 94318117 : r[i] = y;
95 549593 : }
96 :
97 : void
98 1447032 : spv_set_zero (spv_t r, spv_size_t len)
99 : {
100 1447032 : memset (r, 0, len * sizeof (sp_t));
101 1447032 : }
102 :
103 : /* r = x + y */
104 : void
105 2623 : spv_add (spv_t r, spv_t x, spv_t y, spv_size_t len, sp_t m)
106 : {
107 : spv_size_t i;
108 :
109 : ASSERT (r >= x + len || x >= r);
110 : ASSERT (r >= y + len || y >= r);
111 :
112 49719103 : for (i = 0; i < len; i++)
113 49716480 : r[i] = sp_add (x[i], y[i], m);
114 2623 : }
115 :
116 : /* r = [-x[0], -x[1], ... ] */
117 : void
118 6515 : spv_neg (spv_t r, spv_t x, spv_size_t len, sp_t m)
119 : {
120 : spv_size_t i;
121 :
122 16026739 : for (i = 0; i < len; i++)
123 16020224 : r[i] = sp_sub (0, x[i], m);
124 6515 : }
125 :
126 : /* Pointwise multiplication
127 : * r = [x[0] * y[0], x[1] * y[1], ... ] */
128 : void
129 1453040 : spv_pwmul (spv_t r, spv_t x, spv_t y, spv_size_t len, sp_t m, sp_t d)
130 : {
131 1453040 : spv_size_t i = 0;
132 :
133 : ASSERT (r >= x + len || x >= r);
134 : ASSERT (r >= y + len || y >= r);
135 :
136 : #if (defined(__GNUC__) || defined(__ICL)) && \
137 : defined(__i386__) && defined(HAVE_SSE2)
138 :
139 : asm volatile (
140 : "movd %6, %%xmm6 \n\t"
141 : "pshufd $0x44, %%xmm6, %%xmm5 \n\t"
142 : "pshufd $0, %%xmm6, %%xmm6 \n\t"
143 : "movd %7, %%xmm7 \n\t"
144 : "pshufd $0, %%xmm7, %%xmm7 \n\t"
145 :
146 : "0: \n\t"
147 : "movdqa (%1,%4,4), %%xmm0 \n\t"
148 : "movdqa (%2,%4,4), %%xmm2 \n\t"
149 : "pshufd $0x31, %%xmm0, %%xmm1\n\t"
150 : "pshufd $0x31, %%xmm2, %%xmm3\n\t"
151 : "pmuludq %%xmm2, %%xmm0 \n\t"
152 : "pmuludq %%xmm3, %%xmm1 \n\t"
153 :
154 : "movdqa %%xmm0, %%xmm2 \n\t"
155 : "movdqa %%xmm1, %%xmm3 \n\t"
156 : "psrlq $" STRING((2*SP_NUMB_BITS - W_TYPE_SIZE)) ", %%xmm2 \n\t"
157 : "pmuludq %%xmm7, %%xmm2 \n\t"
158 : "psrlq $" STRING((2*SP_NUMB_BITS - W_TYPE_SIZE)) ", %%xmm3 \n\t"
159 : "pmuludq %%xmm7, %%xmm3 \n\t"
160 :
161 : #if SP_NUMB_BITS < W_TYPE_SIZE - 1
162 : "psrlq $33, %%xmm2 \n\t"
163 : "pmuludq %%xmm6, %%xmm2 \n\t"
164 : "psrlq $33, %%xmm3 \n\t"
165 : "pmuludq %%xmm6, %%xmm3 \n\t"
166 : "psubq %%xmm2, %%xmm0 \n\t"
167 : "psubq %%xmm3, %%xmm1 \n\t"
168 : #else
169 : "pshufd $0xf5, %%xmm2, %%xmm2 \n\t"
170 : "pmuludq %%xmm6, %%xmm2 \n\t"
171 : "pshufd $0xf5, %%xmm3, %%xmm3 \n\t"
172 : "pmuludq %%xmm6, %%xmm3 \n\t"
173 : "psubq %%xmm2, %%xmm0 \n\t"
174 : "psubq %%xmm3, %%xmm1 \n\t"
175 :
176 : "psubq %%xmm5, %%xmm0 \n\t"
177 : "psubq %%xmm5, %%xmm1 \n\t"
178 : "pshufd $0xf5, %%xmm0, %%xmm2 \n\t"
179 : "pshufd $0xf5, %%xmm1, %%xmm3 \n\t"
180 : "pand %%xmm5, %%xmm2 \n\t"
181 : "pand %%xmm5, %%xmm3 \n\t"
182 : "paddq %%xmm2, %%xmm0 \n\t"
183 : "paddq %%xmm3, %%xmm1 \n\t"
184 : #endif
185 : "pshufd $0x8, %%xmm0, %%xmm0 \n\t"
186 : "pshufd $0x8, %%xmm1, %%xmm1 \n\t"
187 : "punpckldq %%xmm1, %%xmm0 \n\t"
188 : "psubd %%xmm6, %%xmm0 \n\t"
189 :
190 : "pxor %%xmm1, %%xmm1 \n\t"
191 : "pcmpgtd %%xmm0, %%xmm1 \n\t"
192 : "pand %%xmm6, %%xmm1 \n\t"
193 : "paddd %%xmm1, %%xmm0 \n\t"
194 : "movdqa %%xmm0, (%3,%4,4) \n\t"
195 :
196 : "addl $4, %4 \n\t" /* INC */
197 : "cmpl %5, %4 \n\t"
198 : "jne 0b \n\t"
199 :
200 : :"=r"(i)
201 : :"r"(x), "r"(y), "r"(r), "0"(i),
202 : "g"(len & (spv_size_t)(~3)), "g"(m), "g"(d)
203 : :"%xmm0", "%xmm1", "%xmm2", "%xmm3",
204 : "%xmm5", "%xmm6", "%xmm7", "cc", "memory");
205 : #elif defined( _MSC_VER ) && defined( SSE2)
206 : __asm
207 : { push esi
208 : push edi
209 : mov edi, x
210 : mov esi, y
211 : mov edx, r
212 : xor ecx, ecx
213 : mov eax, len
214 : and eax, ~3
215 :
216 : movd xmm6, m
217 : pshufd xmm5, xmm6, 0x44
218 : pshufd xmm6, xmm6, 0
219 : movd xmm7, d
220 : pshufd xmm7, xmm7, 0
221 :
222 : L0: movdqa xmm0, [edi+ecx*4]
223 : movdqa xmm2, [esi+ecx*4]
224 : pshufd xmm1, xmm0, 0x31
225 : pshufd xmm3, xmm2, 0x31
226 : pmuludq xmm0, xmm2
227 : pmuludq xmm1, xmm3
228 :
229 : movdqa xmm2, xmm0
230 : movdqa xmm3, xmm1
231 : psrlq xmm2, 2*SP_NUMB_BITS - W_TYPE_SIZE
232 : pmuludq xmm2, xmm7
233 : psrlq xmm3, 2*SP_NUMB_BITS - W_TYPE_SIZE
234 : pmuludq xmm3, xmm7
235 :
236 : #if SP_NUMB_BITS < W_TYPE_SIZE - 1
237 : psrlq xmm2, 33
238 : pmuludq xmm2, xmm6
239 : psrlq xmm3, 33
240 : pmuludq xmm3, xmm6
241 : psubq xmm0, xmm2
242 : psubq xmm1, xmm3
243 : #else
244 : pshufd xmm2, xmm2, 0xf5
245 : pmuludq xmm2, xmm6
246 : pshufd xmm3, xmm3, 0xf5
247 : pmuludq xmm3, xmm6
248 : psubq xmm0, xmm2
249 : psubq xmm1, xmm3
250 :
251 : psubq xmm0, xmm5
252 : psubq xmm1, xmm5
253 : pshufd xmm2, xmm0, 0xf5
254 : pshufd xmm3, xmm1, 0xf5
255 : pand xmm2, xmm5
256 : pand xmm3, xmm5
257 : paddq xmm0, xmm2
258 : paddq xmm1, xmm3
259 : #endif
260 : pshufd xmm0, xmm0, 0x8
261 : pshufd xmm1, xmm1, 0x8
262 : punpckldq xmm0, xmm1
263 : psubd xmm0, xmm6
264 :
265 : pxor xmm1, xmm1
266 : pcmpgtd xmm1, xmm0
267 : pand xmm1, xmm6
268 : paddd xmm0, xmm1
269 : movdqa [edx+ecx*4], xmm0
270 :
271 : add ecx, 4
272 : cmp eax, ecx
273 : jne L0
274 : mov i, ecx
275 : pop edi
276 : pop esi
277 : }
278 : #endif
279 :
280 837068320 : for (; i < len; i++)
281 835615280 : r[i] = sp_mul (x[i], y[i], m, d);
282 1453040 : }
283 :
284 : /* dst = src * y */
285 : void
286 22893072 : spv_mul_sp (spv_t r, spv_t x, sp_t c, spv_size_t len, sp_t m, sp_t d)
287 : {
288 22893072 : spv_size_t i = 0;
289 :
290 : ASSERT (r >= x + len || x >= r);
291 :
292 : #if (defined(__GNUC__) || defined(__ICL)) && \
293 : defined(__i386__) && defined(HAVE_SSE2)
294 :
295 : asm volatile (
296 : "movd %2, %%xmm4 \n\t"
297 : "pshufd $0, %%xmm4, %%xmm4 \n\t"
298 : "movd %6, %%xmm6 \n\t"
299 : "pshufd $0x44, %%xmm6, %%xmm5 \n\t"
300 : "pshufd $0, %%xmm6, %%xmm6 \n\t"
301 : "movd %7, %%xmm7 \n\t"
302 : "pshufd $0, %%xmm7, %%xmm7 \n\t"
303 :
304 : "0: \n\t"
305 : "movdqa (%1,%4,4), %%xmm0 \n\t"
306 : "pshufd $0x31, %%xmm0, %%xmm1\n\t"
307 : "pshufd $0x31, %%xmm4, %%xmm3\n\t"
308 : "pmuludq %%xmm4, %%xmm0 \n\t"
309 : "pmuludq %%xmm3, %%xmm1 \n\t"
310 :
311 : "movdqa %%xmm0, %%xmm2 \n\t"
312 : "movdqa %%xmm1, %%xmm3 \n\t"
313 : "psrlq $" STRING((2*SP_NUMB_BITS - W_TYPE_SIZE)) ", %%xmm2 \n\t"
314 : "pmuludq %%xmm7, %%xmm2 \n\t"
315 : "psrlq $" STRING((2*SP_NUMB_BITS - W_TYPE_SIZE)) ", %%xmm3 \n\t"
316 : "pmuludq %%xmm7, %%xmm3 \n\t"
317 :
318 : #if SP_NUMB_BITS < W_TYPE_SIZE - 1
319 : "psrlq $33, %%xmm2 \n\t"
320 : "pmuludq %%xmm6, %%xmm2 \n\t"
321 : "psrlq $33, %%xmm3 \n\t"
322 : "pmuludq %%xmm6, %%xmm3 \n\t"
323 : "psubq %%xmm2, %%xmm0 \n\t"
324 : "psubq %%xmm3, %%xmm1 \n\t"
325 : #else
326 : "pshufd $0xf5, %%xmm2, %%xmm2 \n\t"
327 : "pmuludq %%xmm6, %%xmm2 \n\t"
328 : "pshufd $0xf5, %%xmm3, %%xmm3 \n\t"
329 : "pmuludq %%xmm6, %%xmm3 \n\t"
330 : "psubq %%xmm2, %%xmm0 \n\t"
331 : "psubq %%xmm3, %%xmm1 \n\t"
332 :
333 : "psubq %%xmm5, %%xmm0 \n\t"
334 : "psubq %%xmm5, %%xmm1 \n\t"
335 : "pshufd $0xf5, %%xmm0, %%xmm2 \n\t"
336 : "pshufd $0xf5, %%xmm1, %%xmm3 \n\t"
337 : "pand %%xmm5, %%xmm2 \n\t"
338 : "pand %%xmm5, %%xmm3 \n\t"
339 : "paddq %%xmm2, %%xmm0 \n\t"
340 : "paddq %%xmm3, %%xmm1 \n\t"
341 : #endif
342 : "pshufd $0x8, %%xmm0, %%xmm0 \n\t"
343 : "pshufd $0x8, %%xmm1, %%xmm1 \n\t"
344 : "punpckldq %%xmm1, %%xmm0 \n\t"
345 : "psubd %%xmm6, %%xmm0 \n\t"
346 :
347 : "pxor %%xmm1, %%xmm1 \n\t"
348 : "pcmpgtd %%xmm0, %%xmm1 \n\t"
349 : "pand %%xmm6, %%xmm1 \n\t"
350 : "paddd %%xmm1, %%xmm0 \n\t"
351 : "movdqa %%xmm0, (%3,%4,4) \n\t"
352 :
353 : "addl $4, %4 \n\t" /* INC */
354 : "cmpl %5, %4 \n\t"
355 : "jne 0b \n\t"
356 :
357 : :"=r"(i)
358 : :"r"(x), "g"(c), "r"(r), "0"(i),
359 : "g"(len & (spv_size_t)(~3)), "g"(m), "g"(d)
360 : :"%xmm0", "%xmm1", "%xmm2", "%xmm3",
361 : "%xmm4", "%xmm5", "%xmm6", "%xmm7", "cc", "memory");
362 : #elif defined( _MSC_VER ) && defined( SSE2)
363 : __asm
364 : { push esi
365 : push edi
366 : xor ecx, ecx
367 : mov edi, x
368 : mov esi, c
369 : mov edx, r
370 : mov eax, len
371 : and eax, ~3
372 :
373 : movd xmm4, esi
374 : pshufd xmm4, xmm4, 0
375 : movd xmm6, m
376 : pshufd xmm5, xmm6, 0x44
377 : pshufd xmm6, xmm6, 0
378 : movd xmm7, d
379 : pshufd xmm7, xmm7, 0
380 :
381 : L0: movdqa xmm0, [edi+ecx*4]
382 : pshufd xmm1, xmm0, 0x31
383 : pshufd xmm3, xmm4, 0x31
384 : pmuludq xmm0, xmm4
385 : pmuludq xmm1, xmm3
386 :
387 : movdqa xmm2, xmm0
388 : movdqa xmm3, xmm1
389 : psrlq xmm2, 2*SP_NUMB_BITS - W_TYPE_SIZE
390 : pmuludq xmm2, xmm7
391 : psrlq xmm3, 2*SP_NUMB_BITS - W_TYPE_SIZE
392 : pmuludq xmm3, xmm7
393 :
394 : #if SP_NUMB_BITS < W_TYPE_SIZE - 1
395 : psrlq xmm2, 33
396 : pmuludq xmm2, xmm6
397 : psrlq xmm3, 33
398 : pmuludq xmm3, xmm6
399 : psubq xmm0, xmm2
400 : psubq xmm1, xmm3
401 : #else
402 : pshufd xmm2, xmm2, 0xf5
403 : pmuludq xmm2, xmm6
404 : pshufd xmm3, xmm3, 0xf5
405 : pmuludq xmm3, xmm6
406 : psubq xmm0, xmm2
407 : psubq xmm1, xmm3
408 :
409 : psubq xmm0, xmm5
410 : psubq xmm1, xmm5
411 : pshufd xmm2, xmm0, 0xf5
412 : pshufd xmm3, xmm1, 0xf5
413 : pand xmm2, xmm5
414 : pand xmm3, xmm5
415 : paddq xmm0, xmm2
416 : paddq xmm1, xmm3
417 : #endif
418 : pshufd xmm0, xmm0, 0x8
419 : pshufd xmm1, xmm1, 0x8
420 : punpckldq xmm0, xmm1
421 : psubd xmm0, xmm6
422 :
423 : pxor xmm1, xmm1
424 : pcmpgtd xmm1, xmm0
425 : pand xmm1, xmm6
426 : paddd xmm0, xmm1
427 : movdqa [edx+ecx*4], xmm0
428 :
429 : add ecx, 4
430 : cmp eax, ecx
431 : jne L0
432 : mov i, ecx
433 : pop edi
434 : pop esi
435 : }
436 : #endif
437 :
438 3634381864 : for (; i < len; i++)
439 3611488792 : r[i] = sp_mul (x[i], c, m, d);
440 22893072 : }
441 :
442 : void
443 24 : spv_random (spv_t x, spv_size_t len, sp_t m)
444 : {
445 : spv_size_t i;
446 24 : mpn_random ((mp_ptr) x, len);
447 :
448 6291480 : for (i = 0; i < len; i++)
449 15725866 : while (x[i] >= m)
450 9434410 : x[i] -= m;
451 24 : }
|
