Line data Source code
1 : /* Arithmetic on lists of residues modulo n.
2 :
3 : Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2012
4 : Paul Zimmermann and Alexander Kruppa.
5 :
6 : This file is part of the ECM Library.
7 :
8 : The ECM Library is free software; you can redistribute it and/or modify
9 : it under the terms of the GNU Lesser General Public License as published by
10 : the Free Software Foundation; either version 3 of the License, or (at your
11 : option) any later version.
12 :
13 : The ECM Library is distributed in the hope that it will be useful, but
14 : WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 : or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 : License for more details.
17 :
18 : You should have received a copy of the GNU Lesser General Public License
19 : along with the ECM Library; see the file COPYING.LIB. If not, see
20 : http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 : 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22 :
23 : #include <stdlib.h>
24 : #include "ecm-impl.h"
25 :
26 : #ifdef DEBUG
27 : #define ASSERTD(x) assert(x)
28 : #else
29 : #define ASSERTD(x)
30 : #endif
31 :
32 : extern unsigned int Fermat;
33 :
34 : /* returns a bound on the auxiliary memory needed by list_mult_n */
35 : int
36 78557 : list_mul_mem (unsigned int len)
37 : {
38 78557 : return 2 * len;
39 : }
40 :
41 : /* creates a list of n integers, return NULL if error */
42 : listz_t
43 7618538 : init_list (unsigned int n)
44 : {
45 : listz_t p;
46 : unsigned int i;
47 :
48 7618538 : p = (mpz_t*) malloc (n * sizeof (mpz_t));
49 7618538 : if (p == NULL)
50 0 : return NULL;
51 41933917 : for (i = 0; i < n; i++)
52 34315379 : mpz_init (p[i]);
53 7618538 : return p;
54 : }
55 :
56 : /* creates a list of n integers, return NULL if error. Allocates each
57 : mpz_t to the size of N bits */
58 : listz_t
59 88710 : init_list2 (unsigned int n, unsigned int N)
60 : {
61 : listz_t p;
62 : unsigned int i;
63 :
64 88710 : p = (mpz_t*) malloc (n * sizeof (mpz_t));
65 88710 : if (p == NULL)
66 0 : return NULL;
67 35332033 : for (i = 0; i < n; i++)
68 35243323 : mpz_init2 (p[i], N);
69 88710 : return p;
70 : }
71 :
72 : /* clears a list of n integers */
73 : void
74 7711388 : clear_list (listz_t p, unsigned int n)
75 : {
76 : unsigned int i;
77 :
78 7711388 : if (p == NULL)
79 2054 : return;
80 77323274 : for (i = 0; i < n; i++)
81 69613940 : mpz_clear (p[i]);
82 7709334 : free (p);
83 : }
84 :
85 : #ifdef DEBUG
86 : /* prints a list of n coefficients as a polynomial */
87 : void
88 : print_list (listz_t p, unsigned int n)
89 : {
90 : unsigned int i;
91 :
92 : for (i = 0; i < n; i++)
93 : {
94 : if (i > 0 && mpz_cmp_ui (p[i], 0) >= 0)
95 : fprintf (ECM_STDOUT, "+");
96 : mpz_out_str (ECM_STDOUT, 10, p[i]);
97 : fprintf (ECM_STDOUT, "*x^%u", i);
98 : }
99 : fprintf (ECM_STDOUT, "\n");
100 : }
101 :
102 : static int
103 : list_check (listz_t a, unsigned int l, mpz_t n)
104 : {
105 : unsigned int i;
106 :
107 : for (i = 0; i < l; i++)
108 : if (mpz_cmp_ui (a[i], 0) < 0 || mpz_cmp (n, a[i]) <= 0)
109 : {
110 : fprintf (ECM_STDOUT, "l=%u i=%u\n", l, i);
111 : mpz_out_str (ECM_STDOUT, 10, a[i]);
112 : fprintf (ECM_STDOUT, "\n");
113 : return 0;
114 : }
115 : return 1;
116 : }
117 : #endif /* DEBUG */
118 :
119 : /* Read all entries in list from stream.
120 : Return 0 on success, ECM_ERROR on error */
121 : int
122 113566 : list_inp_raw (listz_t a, FILE *f, unsigned int n)
123 : {
124 : unsigned int i;
125 :
126 688856 : for (i = 0; i < n; i++)
127 575290 : if (mpz_inp_raw (a[i], f) == 0)
128 0 : return ECM_ERROR;
129 :
130 113566 : return 0;
131 : }
132 :
133 : /* Write all entries in list to stream.
134 : Return 0 on success, ECM_ERROR on error */
135 : int
136 1816 : list_out_raw (FILE *f, listz_t a, unsigned int n)
137 : {
138 : unsigned int i;
139 :
140 599634 : for (i = 0; i < n; i++)
141 597818 : if (mpz_out_raw (f, a[i]) == 0)
142 0 : return ECM_ERROR;
143 :
144 1816 : return 0;
145 : }
146 :
147 : /* p <- q */
148 : void
149 444665 : list_set (listz_t p, listz_t q, unsigned int n)
150 : {
151 : unsigned int i;
152 :
153 1897612 : for (i = 0; i < n; i++)
154 1452947 : mpz_set (p[i], q[i]);
155 444665 : }
156 :
157 : /* p[0] <-> p[n-1], p[1] <-> p[n-2], ... */
158 : void
159 12041 : list_revert (listz_t p, unsigned int n)
160 : {
161 : unsigned int i;
162 :
163 3911767 : for (i = 0; i < n - 1 - i; i++)
164 3899726 : mpz_swap (p[i], p[n - 1 - i]);
165 12041 : }
166 :
167 : void
168 1037 : list_swap (listz_t p, listz_t q, unsigned int n)
169 : {
170 : unsigned int i;
171 :
172 1455760 : for (i = 0; i < n; i++)
173 1454723 : mpz_swap (p[i], q[i]);
174 1037 : }
175 :
176 : /* p <- -q, keeps residues normalized */
177 : void
178 7465 : list_neg (listz_t p, listz_t q, unsigned int l, mpz_t n)
179 : {
180 : unsigned int i;
181 :
182 614359 : for (i = 0; i < l; i++)
183 : {
184 606894 : if (mpz_sgn (q[i]))
185 606894 : mpz_sub (p[i], n, q[i]);
186 : else
187 0 : mpz_set_ui (p[i], 0);
188 : }
189 7465 : }
190 :
191 : /* p <- q modulo mod */
192 : void
193 4568183 : list_mod (listz_t p, listz_t q, unsigned int n, mpz_t mod)
194 : {
195 : unsigned int i;
196 :
197 68006515 : for (i = 0; i < n; i++)
198 63438332 : mpz_mod (p[i], q[i], mod);
199 4568183 : }
200 :
201 : /* p <- q + r */
202 : void
203 192863205 : list_add (listz_t p, listz_t q, listz_t r, unsigned int l)
204 : {
205 : unsigned int i;
206 :
207 684698744 : for (i = 0; i < l; i++)
208 491835539 : mpz_add (p[i], q[i], r[i]);
209 192863205 : }
210 :
211 : /* p <- q - r */
212 : void
213 113309501 : list_sub (listz_t p, listz_t q, listz_t r, unsigned int l)
214 : {
215 : unsigned int i;
216 :
217 468717370 : for (i = 0; i < l; i++)
218 355407869 : mpz_sub (p[i], q[i], r[i]);
219 113309501 : }
220 :
221 : /* Multiply up the integers in l, modulo n. Each entry becomes the
222 : product (mod n) of itself and all previous entries */
223 :
224 : void
225 1027 : list_mulup (listz_t l, unsigned int k, mpz_t n, mpz_t t)
226 : {
227 : unsigned int i;
228 :
229 1114755 : for (i = 1; i < k; i++)
230 : {
231 1113728 : mpz_mul (t, l[i - 1], l[i]);
232 1113728 : mpz_mod (l[i], t, n);
233 : }
234 1027 : }
235 :
236 : /* p <- 0 */
237 : void
238 0 : list_zero (listz_t p, unsigned int n)
239 : {
240 : unsigned int i;
241 :
242 0 : for (i = 0; i < n; i++)
243 0 : mpz_set_ui (p[i], 0);
244 0 : }
245 :
246 : /* puts in a[K-1]..a[2K-2] the K high terms of the product
247 : of b[0..K-1] and c[0..K-1].
248 : Assumes K >= 1, and a[0..2K-2] exist.
249 : Needs space for list_mul_mem(K) in t.
250 : */
251 : void
252 75965 : list_mul_high (listz_t a, listz_t b, listz_t c, unsigned int K)
253 : {
254 75965 : list_mult_n (a, b, c, K);
255 75965 : }
256 :
257 : /* multiplies b[0]+...+b[k-1]*x^(k-1)+x^k by c[0]+...+c[l-1]*x^(l-1)+x^l
258 : and puts the results in a[0]+...+a[k+l-1]*x^(k+l-1)
259 : [the leading monomial x^(k+l) is implicit].
260 : If monic is 0, don't consider x^k in b (and x^l in c).
261 : Assumes k = l or k = l+1.
262 : The auxiliary array t contains at least list_mul_mem(l) entries.
263 : a and t should not overlap.
264 : */
265 : void
266 4220832 : list_mul (listz_t a, listz_t b, unsigned int k,
267 : listz_t c, unsigned int l, int monic, listz_t t)
268 : {
269 : unsigned int i, po2;
270 :
271 : ASSERT(k == l || k == l + 1);
272 :
273 8543468 : for (po2 = l; (po2 & 1) == 0; po2 >>= 1);
274 4220832 : po2 = (po2 == 1);
275 :
276 : #ifdef DEBUG
277 : if (Fermat && !(po2 && l == k))
278 : fprintf (ECM_STDOUT, "list_mul: Fermat number, but poly lengths %d and %d\n", k, l);
279 : #endif
280 :
281 4220832 : if (po2 && Fermat)
282 : {
283 459163 : if (monic && l == k)
284 : {
285 459163 : F_mul (a, b, c, l, MONIC, Fermat, t);
286 459163 : monic = 0;
287 : }
288 : else
289 0 : F_mul (a, b, c, l, DEFAULT, Fermat, t);
290 : }
291 : else
292 3761669 : list_mult_n (a, b, c, l); /* set a[0]...a[2l-2] */
293 :
294 4220832 : if (k > l) /* multiply b[l]*x^l by c[0]+...+c[l-1]*x^(l-1) */
295 : {
296 298144 : for (i = 0; i < l - 1; i++)
297 219056 : mpz_addmul (a[l+i], b[l], c[i]);
298 79088 : mpz_mul (a[2*l-1], b[l], c[l-1]);
299 : }
300 :
301 : /* deal with x^k and x^l */
302 4220832 : if (monic)
303 : {
304 3754833 : mpz_set_ui (a[k + l - 1], 0);
305 :
306 : /* Single pass over a[] */
307 :
308 : /* a += b * x^l + c * x^k, so a[i] += b[i-l]; a[i] += c[i-k]
309 : if 0 <= i-l < k or 0 <= i-k < l, respectively */
310 3754833 : if (k > l) /* case k = l+1 */
311 79088 : mpz_add (a[l], a[l], b[0]);
312 20808363 : for (i = k; i < k + l; i++)
313 : {
314 17053530 : mpz_add (a[i], a[i], b[i-l]); /* i-l < k */
315 17053530 : mpz_add (a[i], a[i], c[i-k]); /* i-k < l */
316 : }
317 : }
318 4220832 : }
319 :
320 : /*
321 : Multiplies b[0..k-1] by c[0..k-1], stores the result in a[0..2k-2],
322 : and stores the reduced product in a2[0..2k-2].
323 : (Here, there is no implicit monic leading monomial.)
324 : Requires at least list_mul_mem(k) cells in t.
325 : */
326 : void
327 127 : list_mulmod (listz_t a2, listz_t a, listz_t b, listz_t c, unsigned int k,
328 : listz_t t, mpz_t n)
329 : {
330 : int i;
331 :
332 : /* keep the semicolon on a separate line to silence a warning with clang */
333 912 : for (i = k; (i & 1) == 0; i >>= 1)
334 : ;
335 :
336 : ASSERTD(list_check(b,k,n));
337 : ASSERTD(list_check(c,k,n));
338 127 : if (i == 1 && Fermat)
339 29 : F_mul (a, b, c, k, DEFAULT, Fermat, t);
340 : else
341 98 : list_mult_n (a, b, c, k); /* set a[0]...a[2l-2] */
342 :
343 127 : list_mod (a2, a, 2 * k - 1, n);
344 127 : }
345 :
346 : /* puts in G[0]..G[k-1] the coefficients from (x+a[0])...(x+a[k-1])
347 : Warning: doesn't fill the coefficient 1 of G[k], which is implicit.
348 : Needs k + list_mul_mem(k/2) cells in T.
349 : G == a is allowed. T must not overlap with anything else.
350 : */
351 : void
352 5606412 : PolyFromRoots (listz_t G, listz_t a, unsigned int k, listz_t T, mpz_t n)
353 : {
354 : unsigned int l, m;
355 :
356 : ASSERT (T != G && T != a);
357 : ASSERT (k >= 1);
358 :
359 5606412 : if (k == 1)
360 : {
361 : /* we consider x + a[0], which mean we consider negated roots */
362 2807378 : mpz_mod (G[0], a[0], n);
363 2807378 : return;
364 : }
365 :
366 2799034 : m = k / 2; /* m >= 1 */
367 2799034 : l = k - m; /* l >= 1 */
368 :
369 2799034 : PolyFromRoots (G, a, l, T, n);
370 2799034 : PolyFromRoots (G + l, a + l, m, T, n);
371 2799034 : list_mul (T, G, l, G + l, m, 1, T + k);
372 2799034 : list_mod (G, T, k, n);
373 : }
374 :
375 : /* puts in G[0]..G[k-1] the coefficients from (x+a[0])...(x+a[k-1])
376 : Warning: doesn't fill the coefficient 1 of G[k], which is implicit.
377 : Needs k + list_mul_mem(k/2) cells in T.
378 : The product tree is stored in:
379 : G[0..k-1] (degree k)
380 : Tree[0][0..k-1] (degree k/2)
381 : Tree[1][0..k-1] (degree k/4), ...,
382 : Tree[lgk-1][0..k-1] (degree 1)
383 : (then we should have initially Tree[lgk-1] = a).
384 :
385 : The parameter dolvl signals that only level 'dolvl' of
386 : the tree should be computed (dolvl < 0 means all levels).
387 :
388 : Either Tree <> NULL and TreeFile == NULL, and we write the tree to memory,
389 : or Tree == NULL and TreeFile <> NULL, and we write the tree to disk.
390 : */
391 : int
392 257544 : PolyFromRoots_Tree (listz_t G, listz_t a, unsigned int k, listz_t T,
393 : int dolvl, mpz_t n, listz_t *Tree, FILE *TreeFile,
394 : unsigned int sh)
395 : {
396 : unsigned int l, m;
397 : listz_t H1, *NextTree;
398 :
399 : ASSERT (k >= 1);
400 :
401 257544 : if (k == 1)
402 : {
403 : /* we consider x + a[0], which mean we consider negated roots */
404 128384 : mpz_mod (G[0], a[0], n);
405 128384 : return 0;
406 : }
407 :
408 129160 : if (Tree == NULL) /* -treefile case */
409 : {
410 1076 : H1 = G;
411 1076 : NextTree = NULL;
412 : }
413 : else
414 : {
415 128084 : H1 = Tree[0] + sh;
416 128084 : NextTree = Tree + 1;
417 : }
418 :
419 129160 : m = k / 2;
420 129160 : l = k - m;
421 :
422 129160 : if (dolvl != 0) /* either dolvl < 0 and we need to compute all levels,
423 : or dolvl > 0 and we need first to compute lower levels */
424 : {
425 128706 : PolyFromRoots_Tree (H1, a, l, T, dolvl - 1, n, NextTree, TreeFile, sh);
426 128706 : PolyFromRoots_Tree (H1 + l, a + l, m, T, dolvl - 1, n, NextTree,
427 : TreeFile, sh + l);
428 : }
429 129160 : if (dolvl <= 0)
430 : {
431 : /* Write this level to disk, if requested */
432 128538 : if (TreeFile != NULL)
433 : {
434 908 : if (list_out_raw (TreeFile, H1, l) == ECM_ERROR ||
435 454 : list_out_raw (TreeFile, H1 + l, m) == ECM_ERROR)
436 : {
437 0 : outputf (OUTPUT_ERROR, "Error writing product tree of F\n");
438 0 : return ECM_ERROR;
439 : }
440 : }
441 128538 : list_mul (T, H1, l, H1 + l, m, 1, T + k);
442 128538 : list_mod (G, T, k, n);
443 : }
444 :
445 129160 : return 0;
446 : }
447 :
448 : /* puts in q[0..K-1] the quotient of x^(2K-2) by B
449 : where B = b[0]+b[1]*x+...+b[K-1]*x^(K-1) with b[K-1]=1.
450 : */
451 : void
452 8545 : PolyInvert (listz_t q, listz_t b, unsigned int K, listz_t t, mpz_t n)
453 : {
454 8545 : if (K == 1)
455 : {
456 1080 : mpz_set_ui (q[0], 1);
457 1080 : return;
458 : }
459 : else
460 : {
461 7465 : int k, l, po2, use_middle_product = 0;
462 :
463 7465 : use_middle_product = 1;
464 :
465 7465 : k = K / 2;
466 7465 : l = K - k;
467 :
468 37485 : for (po2 = K; (po2 & 1) == 0; po2 >>= 1);
469 7465 : po2 = (po2 == 1 && Fermat != 0);
470 :
471 : /* first determine l most-significant coeffs of Q */
472 7465 : PolyInvert (q + k, b + k, l, t, n); /* Q1 = {q+k, l} */
473 :
474 : /* now Q1 * B = x^(2K-2) + O(x^(2K-2-l)) = x^(2K-2) + O(x^(K+k-2)).
475 : We need the coefficients of degree K-1 to K+k-2 of Q1*B */
476 :
477 : ASSERTD(list_check(q+k,l,n) && list_check(b,l,n));
478 7465 : if (po2 == 0 && use_middle_product)
479 : {
480 7426 : TMulKS (t, k - 1, q + k, l - 1, b, K - 1, n, 0);
481 7426 : list_neg (t, t, k, n);
482 : }
483 39 : else if (po2)
484 : {
485 39 : list_revert (q + k, l);
486 : /* This expects the leading monomials explicitly in q[2k-1] and b[k+l-1] */
487 39 : F_mul_trans (t, q + k, b, K / 2, K, Fermat, t + k);
488 39 : list_revert (q + k, l);
489 39 : list_neg (t, t, k, n);
490 : }
491 : else
492 : {
493 0 : list_mult_n (t, q + k, b, l); /* t[0..2l-1] = Q1 * B0 */
494 0 : list_neg (t, t + l - 1, k, n);
495 :
496 0 : if (k > 1)
497 : {
498 0 : list_mul (t + k, q + k, l - 1, b + l, k - 1, 1,
499 0 : t + k + K - 2); /* Q1 * B1 */
500 0 : list_sub (t + 1, t + 1, t + k, k - 1);
501 : }
502 : }
503 7465 : list_mod (t, t, k, n); /* high(1-B*Q1) */
504 :
505 : ASSERTD(list_check(t,k,n) && list_check(q+l,k,n));
506 7465 : if (po2)
507 39 : F_mul (t + k, t, q + l, k, DEFAULT, Fermat, t + 3 * k);
508 : else
509 7426 : list_mult_n (t + k, t, q + l, k);
510 7465 : list_mod (q, t + 2 * k - 1, k, n);
511 : }
512 : }
513 :
514 : /*
515 : Returns in a[0]+a[1]*x+...+a[K-1]*x^(K-1)
516 : the remainder of the division of
517 : A = a[0]+a[1]*x+...+a[2K-2]*x^(2K-2)
518 : by B = b[0]+b[1]*x+...+b[K-1]*x^(K-1)+b[K]*x^K with b[K]=1 *explicit*.
519 : (We have A = Q*B + R with deg(Q)=K-2 and deg(R)=K-1.)
520 : Assumes invb[0]+invb[1]*x+...+invb[K-2]*x^(K-2) equals Quo(x^(2K-2), B).
521 : Assumes K >= 2.
522 : Requires 2K-1 + list_mul_mem(K) cells in t.
523 :
524 : Notations: R = r[0..K-1], A = a[0..2K-2], low(A) = a[0..K-1],
525 : high(A) = a[K..2K-2], Q = t[0..K-2]
526 : Return non-zero iff an error occurred.
527 : */
528 : int
529 75877 : PrerevertDivision (listz_t a, listz_t b, listz_t invb,
530 : unsigned int K, listz_t t, mpz_t n)
531 : {
532 : int po2, wrap;
533 75877 : listz_t t2 = NULL;
534 :
535 75877 : wrap = ks_wrapmul_m (K + 1, K + 1, n) <= 2 * K - 1 + list_mul_mem (K);
536 :
537 : /* Q <- high(high(A) * INVB) with a short product */
538 251591 : for (po2 = K; (po2 & 1) == 0; po2 >>= 1);
539 75877 : po2 = (po2 == 1);
540 75877 : if (Fermat && po2)
541 : {
542 29 : mpz_set_ui (a[2 * K - 1], 0);
543 29 : if (K <= 4 * Fermat)
544 : {
545 26 : F_mul (t, a + K, invb, K, DEFAULT, Fermat, t + 2 * K);
546 : /* Put Q in T, as we still need high(A) later on */
547 26 : list_mod (t, t + K - 2, K, n);
548 : }
549 : else
550 : {
551 3 : F_mul (t, a + K, invb, K, DEFAULT, Fermat, t + 2 * K);
552 3 : list_mod (a + K, t + K - 2, K, n);
553 : }
554 : }
555 : else /* non-Fermat case */
556 : {
557 75848 : list_mul_high (t, a + K, invb, K - 1);
558 : /* the high part of A * INVB is now in {t+K-2, K-1} */
559 75848 : if (wrap)
560 : {
561 75829 : t2 = init_list2 (K - 1, mpz_sizeinbase (n, 2));
562 75829 : ASSERT_ALWAYS(t2 != NULL);
563 75829 : list_mod (t2, t + K - 2, K - 1, n);
564 : }
565 : else /* we can store in high(A) which is no longer needed */
566 19 : list_mod (a + K, t + K - 2, K - 1, n);
567 : }
568 :
569 : /* the quotient Q = trunc(A / B) has degree K-2, i.e. K-1 terms */
570 :
571 : /* T <- low(Q * B) with a short product */
572 75877 : mpz_set_ui (a[2 * K - 1], 0);
573 75877 : if (Fermat && po2)
574 : {
575 29 : if (K <= 4 * Fermat)
576 : {
577 : /* Multiply without zero padding, result is (mod x^K - 1) */
578 26 : F_mul (t + K, t, b, K, NOPAD, Fermat, t + 2 * K);
579 : /* Take the leading monomial x^K of B into account */
580 26 : list_add (t, t + K, t, K);
581 : /* Subtract high(A) */
582 26 : list_sub(t, t, a + K, K);
583 : }
584 : else
585 3 : F_mul (t, a + K, b, K, DEFAULT, Fermat, t + 2 * K);
586 : }
587 : else /* non-Fermat case */
588 : {
589 75848 : if (wrap)
590 : /* Q = {t2, K-1}, B = {b, K+1}
591 : We know that Q*B vanishes with the coefficients of degree
592 : K to 2K-2 of {A, 2K-1} */
593 : {
594 : unsigned int m;
595 75829 : m = ks_wrapmul (t, K + 1, b, K + 1, t2, K - 1, n);
596 75829 : clear_list (t2, K - 1);
597 : /* coefficients of degree m..2K-2 wrap around,
598 : i.e. were added to 0..2K-2-m */
599 75829 : if (m < 2 * K - 1) /* otherwise product is exact */
600 8878 : list_sub (t, t, a + m, 2 * K - 1 - m);
601 : }
602 : else
603 19 : list_mult_n (t, a + K, b, K);
604 : }
605 :
606 : /* now {t, K} contains the low K terms from Q*B */
607 75877 : list_sub (a, a, t, K);
608 75877 : list_mod (a, a, K, n);
609 :
610 75877 : return 0;
611 : }
|
