1 | |
---|
2 | #include <NTL/ZZ_pXFactoring.h> |
---|
3 | #include <NTL/vec_ZZVec.h> |
---|
4 | #include <NTL/fileio.h> |
---|
5 | #include <NTL/FacVec.h> |
---|
6 | |
---|
7 | #include <stdio.h> |
---|
8 | |
---|
9 | #include <NTL/new.h> |
---|
10 | |
---|
11 | NTL_START_IMPL |
---|
12 | |
---|
13 | |
---|
14 | |
---|
15 | |
---|
16 | |
---|
17 | void SquareFreeDecomp(vec_pair_ZZ_pX_long& u, const ZZ_pX& ff) |
---|
18 | { |
---|
19 | ZZ_pX f = ff; |
---|
20 | |
---|
21 | if (!IsOne(LeadCoeff(f))) |
---|
22 | Error("SquareFreeDecomp: bad args"); |
---|
23 | |
---|
24 | ZZ_pX r, t, v, tmp1; |
---|
25 | long m, j, finished, done; |
---|
26 | |
---|
27 | u.SetLength(0); |
---|
28 | |
---|
29 | if (deg(f) == 0) |
---|
30 | return; |
---|
31 | |
---|
32 | m = 1; |
---|
33 | finished = 0; |
---|
34 | |
---|
35 | do { |
---|
36 | j = 1; |
---|
37 | diff(tmp1, f); |
---|
38 | GCD(r, f, tmp1); |
---|
39 | div(t, f, r); |
---|
40 | |
---|
41 | if (deg(t) > 0) { |
---|
42 | done = 0; |
---|
43 | do { |
---|
44 | GCD(v, r, t); |
---|
45 | div(tmp1, t, v); |
---|
46 | if (deg(tmp1) > 0) append(u, cons(tmp1, j*m)); |
---|
47 | if (deg(v) > 0) { |
---|
48 | div(r, r, v); |
---|
49 | t = v; |
---|
50 | j++; |
---|
51 | } |
---|
52 | else |
---|
53 | done = 1; |
---|
54 | } while (!done); |
---|
55 | if (deg(r) == 0) finished = 1; |
---|
56 | } |
---|
57 | |
---|
58 | if (!finished) { |
---|
59 | /* r is a p-th power */ |
---|
60 | long p, k, d; |
---|
61 | conv(p, ZZ_p::modulus()); |
---|
62 | d = deg(r)/p; |
---|
63 | f.rep.SetLength(d+1); |
---|
64 | for (k = 0; k <= d; k++) |
---|
65 | f.rep[k] = r.rep[k*p]; |
---|
66 | m = m*p; |
---|
67 | } |
---|
68 | } while (!finished); |
---|
69 | } |
---|
70 | |
---|
71 | |
---|
72 | |
---|
73 | static |
---|
74 | void NullSpace(long& r, vec_long& D, vec_ZZVec& M, long verbose) |
---|
75 | { |
---|
76 | long k, l, n; |
---|
77 | long i, j; |
---|
78 | long pos; |
---|
79 | ZZ t1, t2; |
---|
80 | ZZ *x, *y; |
---|
81 | |
---|
82 | const ZZ& p = ZZ_p::modulus(); |
---|
83 | |
---|
84 | n = M.length(); |
---|
85 | |
---|
86 | D.SetLength(n); |
---|
87 | for (j = 0; j < n; j++) D[j] = -1; |
---|
88 | |
---|
89 | r = 0; |
---|
90 | |
---|
91 | l = 0; |
---|
92 | for (k = 0; k < n; k++) { |
---|
93 | |
---|
94 | pos = -1; |
---|
95 | for (i = l; i < n; i++) { |
---|
96 | rem(t1, M[i][k], p); |
---|
97 | M[i][k] = t1; |
---|
98 | if (pos == -1 && !IsZero(t1)) |
---|
99 | pos = i; |
---|
100 | } |
---|
101 | |
---|
102 | if (pos != -1) { |
---|
103 | swap(M[pos], M[l]); |
---|
104 | |
---|
105 | // make M[l, k] == -1 mod p, and make row l reduced |
---|
106 | |
---|
107 | InvMod(t1, M[l][k], p); |
---|
108 | NegateMod(t1, t1, p); |
---|
109 | for (j = k+1; j < n; j++) { |
---|
110 | rem(t2, M[l][j], p); |
---|
111 | MulMod(M[l][j], t2, t1, p); |
---|
112 | } |
---|
113 | |
---|
114 | for (i = l+1; i < n; i++) { |
---|
115 | // M[i] = M[i] + M[l]*M[i,k] |
---|
116 | |
---|
117 | t1 = M[i][k]; // this is already reduced |
---|
118 | |
---|
119 | x = M[i].elts() + (k+1); |
---|
120 | y = M[l].elts() + (k+1); |
---|
121 | |
---|
122 | for (j = k+1; j < n; j++, x++, y++) { |
---|
123 | // *x = *x + (*y)*t1 |
---|
124 | |
---|
125 | mul(t2, *y, t1); |
---|
126 | add(*x, *x, t2); |
---|
127 | } |
---|
128 | } |
---|
129 | |
---|
130 | D[k] = l; // variable k is defined by row l |
---|
131 | l++; |
---|
132 | |
---|
133 | } |
---|
134 | else { |
---|
135 | r++; |
---|
136 | } |
---|
137 | } |
---|
138 | } |
---|
139 | |
---|
140 | |
---|
141 | |
---|
142 | static |
---|
143 | void BuildMatrix(vec_ZZVec& M, long n, const ZZ_pX& g, const ZZ_pXModulus& F, |
---|
144 | long verbose) |
---|
145 | { |
---|
146 | long i, j, m; |
---|
147 | ZZ_pXMultiplier G; |
---|
148 | ZZ_pX h; |
---|
149 | |
---|
150 | ZZ t; |
---|
151 | sqr(t, ZZ_p::modulus()); |
---|
152 | mul(t, t, n); |
---|
153 | |
---|
154 | long size = t.size(); |
---|
155 | |
---|
156 | M.SetLength(n); |
---|
157 | for (i = 0; i < n; i++) |
---|
158 | M[i].SetSize(n, size); |
---|
159 | |
---|
160 | build(G, g, F); |
---|
161 | |
---|
162 | set(h); |
---|
163 | for (j = 0; j < n; j++) { |
---|
164 | |
---|
165 | m = deg(h); |
---|
166 | for (i = 0; i < n; i++) { |
---|
167 | if (i <= m) |
---|
168 | M[i][j] = rep(h.rep[i]); |
---|
169 | else |
---|
170 | clear(M[i][j]); |
---|
171 | } |
---|
172 | |
---|
173 | if (j < n-1) |
---|
174 | MulMod(h, h, G, F); |
---|
175 | } |
---|
176 | |
---|
177 | for (i = 0; i < n; i++) |
---|
178 | AddMod(M[i][i], M[i][i], -1, ZZ_p::modulus()); |
---|
179 | |
---|
180 | } |
---|
181 | |
---|
182 | |
---|
183 | |
---|
184 | static |
---|
185 | void RecFindRoots(vec_ZZ_p& x, const ZZ_pX& f) |
---|
186 | { |
---|
187 | if (deg(f) == 0) return; |
---|
188 | |
---|
189 | if (deg(f) == 1) { |
---|
190 | long k = x.length(); |
---|
191 | x.SetLength(k+1); |
---|
192 | negate(x[k], ConstTerm(f)); |
---|
193 | return; |
---|
194 | } |
---|
195 | |
---|
196 | ZZ_pX h; |
---|
197 | |
---|
198 | ZZ_p r; |
---|
199 | ZZ p1; |
---|
200 | |
---|
201 | |
---|
202 | RightShift(p1, ZZ_p::modulus(), 1); |
---|
203 | |
---|
204 | { |
---|
205 | ZZ_pXModulus F; |
---|
206 | build(F, f); |
---|
207 | |
---|
208 | do { |
---|
209 | random(r); |
---|
210 | PowerXPlusAMod(h, r, p1, F); |
---|
211 | add(h, h, -1); |
---|
212 | GCD(h, h, f); |
---|
213 | } while (deg(h) <= 0 || deg(h) == deg(f)); |
---|
214 | } |
---|
215 | |
---|
216 | RecFindRoots(x, h); |
---|
217 | div(h, f, h); |
---|
218 | RecFindRoots(x, h); |
---|
219 | } |
---|
220 | |
---|
221 | void FindRoots(vec_ZZ_p& x, const ZZ_pX& ff) |
---|
222 | { |
---|
223 | ZZ_pX f = ff; |
---|
224 | |
---|
225 | if (!IsOne(LeadCoeff(f))) |
---|
226 | Error("FindRoots: bad args"); |
---|
227 | |
---|
228 | x.SetMaxLength(deg(f)); |
---|
229 | x.SetLength(0); |
---|
230 | RecFindRoots(x, f); |
---|
231 | } |
---|
232 | |
---|
233 | |
---|
234 | static |
---|
235 | void RandomBasisElt(ZZ_pX& g, const vec_long& D, const vec_ZZVec& M) |
---|
236 | { |
---|
237 | ZZ t1, t2; |
---|
238 | |
---|
239 | long n = D.length(); |
---|
240 | |
---|
241 | long i, j, s; |
---|
242 | |
---|
243 | g.rep.SetLength(n); |
---|
244 | |
---|
245 | vec_ZZ_p& v = g.rep; |
---|
246 | |
---|
247 | for (j = n-1; j >= 0; j--) { |
---|
248 | if (D[j] == -1) |
---|
249 | random(v[j]); |
---|
250 | else { |
---|
251 | i = D[j]; |
---|
252 | |
---|
253 | // v[j] = sum_{s=j+1}^{n-1} v[s]*M[i,s] |
---|
254 | |
---|
255 | clear(t1); |
---|
256 | |
---|
257 | for (s = j+1; s < n; s++) { |
---|
258 | mul(t2, rep(v[s]), M[i][s]); |
---|
259 | add(t1, t1, t2); |
---|
260 | } |
---|
261 | |
---|
262 | conv(v[j], t1); |
---|
263 | } |
---|
264 | } |
---|
265 | |
---|
266 | g.normalize(); |
---|
267 | } |
---|
268 | |
---|
269 | |
---|
270 | |
---|
271 | static |
---|
272 | void split(ZZ_pX& f1, ZZ_pX& g1, ZZ_pX& f2, ZZ_pX& g2, |
---|
273 | const ZZ_pX& f, const ZZ_pX& g, |
---|
274 | const vec_ZZ_p& roots, long lo, long mid) |
---|
275 | { |
---|
276 | long r = mid-lo+1; |
---|
277 | |
---|
278 | ZZ_pXModulus F; |
---|
279 | build(F, f); |
---|
280 | |
---|
281 | vec_ZZ_p lroots(INIT_SIZE, r); |
---|
282 | long i; |
---|
283 | |
---|
284 | for (i = 0; i < r; i++) |
---|
285 | lroots[i] = roots[lo+i]; |
---|
286 | |
---|
287 | |
---|
288 | ZZ_pX h, a, d; |
---|
289 | BuildFromRoots(h, lroots); |
---|
290 | CompMod(a, h, g, F); |
---|
291 | |
---|
292 | |
---|
293 | GCD(f1, a, f); |
---|
294 | |
---|
295 | div(f2, f, f1); |
---|
296 | |
---|
297 | rem(g1, g, f1); |
---|
298 | rem(g2, g, f2); |
---|
299 | } |
---|
300 | |
---|
301 | static |
---|
302 | void RecFindFactors(vec_ZZ_pX& factors, const ZZ_pX& f, const ZZ_pX& g, |
---|
303 | const vec_ZZ_p& roots, long lo, long hi) |
---|
304 | { |
---|
305 | long r = hi-lo+1; |
---|
306 | |
---|
307 | if (r == 0) return; |
---|
308 | |
---|
309 | if (r == 1) { |
---|
310 | append(factors, f); |
---|
311 | return; |
---|
312 | } |
---|
313 | |
---|
314 | ZZ_pX f1, g1, f2, g2; |
---|
315 | |
---|
316 | long mid = (lo+hi)/2; |
---|
317 | |
---|
318 | split(f1, g1, f2, g2, f, g, roots, lo, mid); |
---|
319 | |
---|
320 | RecFindFactors(factors, f1, g1, roots, lo, mid); |
---|
321 | RecFindFactors(factors, f2, g2, roots, mid+1, hi); |
---|
322 | } |
---|
323 | |
---|
324 | |
---|
325 | static |
---|
326 | void FindFactors(vec_ZZ_pX& factors, const ZZ_pX& f, const ZZ_pX& g, |
---|
327 | const vec_ZZ_p& roots) |
---|
328 | { |
---|
329 | long r = roots.length(); |
---|
330 | |
---|
331 | factors.SetMaxLength(r); |
---|
332 | factors.SetLength(0); |
---|
333 | |
---|
334 | RecFindFactors(factors, f, g, roots, 0, r-1); |
---|
335 | } |
---|
336 | |
---|
337 | #if 0 |
---|
338 | |
---|
339 | static |
---|
340 | void IterFindFactors(vec_ZZ_pX& factors, const ZZ_pX& f, |
---|
341 | const ZZ_pX& g, const vec_ZZ_p& roots) |
---|
342 | { |
---|
343 | long r = roots.length(); |
---|
344 | long i; |
---|
345 | ZZ_pX h; |
---|
346 | |
---|
347 | factors.SetLength(r); |
---|
348 | |
---|
349 | for (i = 0; i < r; i++) { |
---|
350 | sub(h, g, roots[i]); |
---|
351 | GCD(factors[i], f, h); |
---|
352 | } |
---|
353 | } |
---|
354 | |
---|
355 | #endif |
---|
356 | |
---|
357 | |
---|
358 | |
---|
359 | |
---|
360 | void SFBerlekamp(vec_ZZ_pX& factors, const ZZ_pX& ff, long verbose) |
---|
361 | { |
---|
362 | ZZ_pX f = ff; |
---|
363 | |
---|
364 | if (!IsOne(LeadCoeff(f))) |
---|
365 | Error("SFBerlekamp: bad args"); |
---|
366 | |
---|
367 | if (deg(f) == 0) { |
---|
368 | factors.SetLength(0); |
---|
369 | return; |
---|
370 | } |
---|
371 | |
---|
372 | if (deg(f) == 1) { |
---|
373 | factors.SetLength(1); |
---|
374 | factors[0] = f; |
---|
375 | return; |
---|
376 | } |
---|
377 | |
---|
378 | double t; |
---|
379 | |
---|
380 | const ZZ& p = ZZ_p::modulus(); |
---|
381 | |
---|
382 | long n = deg(f); |
---|
383 | |
---|
384 | ZZ_pXModulus F; |
---|
385 | |
---|
386 | build(F, f); |
---|
387 | |
---|
388 | ZZ_pX g, h; |
---|
389 | |
---|
390 | PowerXMod(g, p, F); |
---|
391 | |
---|
392 | vec_long D; |
---|
393 | long r; |
---|
394 | |
---|
395 | vec_ZZVec M; |
---|
396 | |
---|
397 | BuildMatrix(M, n, g, F, verbose); |
---|
398 | |
---|
399 | NullSpace(r, D, M, verbose); |
---|
400 | |
---|
401 | if (r == 1) { |
---|
402 | factors.SetLength(1); |
---|
403 | factors[0] = f; |
---|
404 | return; |
---|
405 | } |
---|
406 | |
---|
407 | vec_ZZ_p roots; |
---|
408 | |
---|
409 | RandomBasisElt(g, D, M); |
---|
410 | MinPolyMod(h, g, F, r); |
---|
411 | if (deg(h) == r) M.kill(); |
---|
412 | FindRoots(roots, h); |
---|
413 | FindFactors(factors, f, g, roots); |
---|
414 | |
---|
415 | ZZ_pX g1; |
---|
416 | vec_ZZ_pX S, S1; |
---|
417 | long i; |
---|
418 | |
---|
419 | while (factors.length() < r) { |
---|
420 | RandomBasisElt(g, D, M); |
---|
421 | S.kill(); |
---|
422 | for (i = 0; i < factors.length(); i++) { |
---|
423 | const ZZ_pX& f = factors[i]; |
---|
424 | if (deg(f) == 1) { |
---|
425 | append(S, f); |
---|
426 | continue; |
---|
427 | } |
---|
428 | build(F, f); |
---|
429 | rem(g1, g, F); |
---|
430 | if (deg(g1) <= 0) { |
---|
431 | append(S, f); |
---|
432 | continue; |
---|
433 | } |
---|
434 | MinPolyMod(h, g1, F, min(deg(f), r-factors.length()+1)); |
---|
435 | FindRoots(roots, h); |
---|
436 | S1.kill(); |
---|
437 | FindFactors(S1, f, g1, roots); |
---|
438 | append(S, S1); |
---|
439 | } |
---|
440 | swap(factors, S); |
---|
441 | } |
---|
442 | |
---|
443 | } |
---|
444 | |
---|
445 | |
---|
446 | void berlekamp(vec_pair_ZZ_pX_long& factors, const ZZ_pX& f, long verbose) |
---|
447 | { |
---|
448 | double t; |
---|
449 | vec_pair_ZZ_pX_long sfd; |
---|
450 | vec_ZZ_pX x; |
---|
451 | |
---|
452 | if (!IsOne(LeadCoeff(f))) |
---|
453 | Error("berlekamp: bad args"); |
---|
454 | |
---|
455 | |
---|
456 | SquareFreeDecomp(sfd, f); |
---|
457 | |
---|
458 | factors.SetLength(0); |
---|
459 | |
---|
460 | long i, j; |
---|
461 | |
---|
462 | for (i = 0; i < sfd.length(); i++) { |
---|
463 | |
---|
464 | SFBerlekamp(x, sfd[i].a, verbose); |
---|
465 | |
---|
466 | for (j = 0; j < x.length(); j++) |
---|
467 | append(factors, cons(x[j], sfd[i].b)); |
---|
468 | } |
---|
469 | } |
---|
470 | |
---|
471 | |
---|
472 | |
---|
473 | static |
---|
474 | void AddFactor(vec_pair_ZZ_pX_long& factors, const ZZ_pX& g, long d, long verbose) |
---|
475 | { |
---|
476 | append(factors, cons(g, d)); |
---|
477 | } |
---|
478 | |
---|
479 | static |
---|
480 | void ProcessTable(ZZ_pX& f, vec_pair_ZZ_pX_long& factors, |
---|
481 | const ZZ_pXModulus& F, long limit, const vec_ZZ_pX& tbl, |
---|
482 | long d, long verbose) |
---|
483 | |
---|
484 | { |
---|
485 | if (limit == 0) return; |
---|
486 | |
---|
487 | ZZ_pX t1; |
---|
488 | |
---|
489 | if (limit == 1) { |
---|
490 | GCD(t1, f, tbl[0]); |
---|
491 | if (deg(t1) > 0) { |
---|
492 | AddFactor(factors, t1, d, verbose); |
---|
493 | div(f, f, t1); |
---|
494 | } |
---|
495 | |
---|
496 | return; |
---|
497 | } |
---|
498 | |
---|
499 | long i; |
---|
500 | |
---|
501 | t1 = tbl[0]; |
---|
502 | for (i = 1; i < limit; i++) |
---|
503 | MulMod(t1, t1, tbl[i], F); |
---|
504 | |
---|
505 | GCD(t1, f, t1); |
---|
506 | |
---|
507 | if (deg(t1) == 0) return; |
---|
508 | |
---|
509 | div(f, f, t1); |
---|
510 | |
---|
511 | ZZ_pX t2; |
---|
512 | |
---|
513 | i = 0; |
---|
514 | d = d - limit + 1; |
---|
515 | |
---|
516 | while (2*d <= deg(t1)) { |
---|
517 | GCD(t2, tbl[i], t1); |
---|
518 | if (deg(t2) > 0) { |
---|
519 | AddFactor(factors, t2, d, verbose); |
---|
520 | div(t1, t1, t2); |
---|
521 | } |
---|
522 | |
---|
523 | i++; |
---|
524 | d++; |
---|
525 | } |
---|
526 | |
---|
527 | if (deg(t1) > 0) |
---|
528 | AddFactor(factors, t1, deg(t1), verbose); |
---|
529 | } |
---|
530 | |
---|
531 | |
---|
532 | void TraceMap(ZZ_pX& w, const ZZ_pX& a, long d, const ZZ_pXModulus& F, |
---|
533 | const ZZ_pX& b) |
---|
534 | |
---|
535 | { |
---|
536 | if (d < 0) Error("TraceMap: bad args"); |
---|
537 | |
---|
538 | ZZ_pX y, z, t; |
---|
539 | |
---|
540 | z = b; |
---|
541 | y = a; |
---|
542 | clear(w); |
---|
543 | |
---|
544 | while (d) { |
---|
545 | if (d == 1) { |
---|
546 | if (IsZero(w)) |
---|
547 | w = y; |
---|
548 | else { |
---|
549 | CompMod(w, w, z, F); |
---|
550 | add(w, w, y); |
---|
551 | } |
---|
552 | } |
---|
553 | else if ((d & 1) == 0) { |
---|
554 | Comp2Mod(z, t, z, y, z, F); |
---|
555 | add(y, t, y); |
---|
556 | } |
---|
557 | else if (IsZero(w)) { |
---|
558 | w = y; |
---|
559 | Comp2Mod(z, t, z, y, z, F); |
---|
560 | add(y, t, y); |
---|
561 | } |
---|
562 | else { |
---|
563 | Comp3Mod(z, t, w, z, y, w, z, F); |
---|
564 | add(w, w, y); |
---|
565 | add(y, t, y); |
---|
566 | } |
---|
567 | |
---|
568 | d = d >> 1; |
---|
569 | } |
---|
570 | } |
---|
571 | |
---|
572 | |
---|
573 | void PowerCompose(ZZ_pX& y, const ZZ_pX& h, long q, const ZZ_pXModulus& F) |
---|
574 | { |
---|
575 | if (q < 0) Error("PowerCompose: bad args"); |
---|
576 | |
---|
577 | ZZ_pX z(INIT_SIZE, F.n); |
---|
578 | long sw; |
---|
579 | |
---|
580 | z = h; |
---|
581 | SetX(y); |
---|
582 | |
---|
583 | while (q) { |
---|
584 | sw = 0; |
---|
585 | |
---|
586 | if (q > 1) sw = 2; |
---|
587 | if (q & 1) { |
---|
588 | if (IsX(y)) |
---|
589 | y = z; |
---|
590 | else |
---|
591 | sw = sw | 1; |
---|
592 | } |
---|
593 | |
---|
594 | switch (sw) { |
---|
595 | case 0: |
---|
596 | break; |
---|
597 | |
---|
598 | case 1: |
---|
599 | CompMod(y, y, z, F); |
---|
600 | break; |
---|
601 | |
---|
602 | case 2: |
---|
603 | CompMod(z, z, z, F); |
---|
604 | break; |
---|
605 | |
---|
606 | case 3: |
---|
607 | Comp2Mod(y, z, y, z, z, F); |
---|
608 | break; |
---|
609 | } |
---|
610 | |
---|
611 | q = q >> 1; |
---|
612 | } |
---|
613 | } |
---|
614 | |
---|
615 | |
---|
616 | long ProbIrredTest(const ZZ_pX& f, long iter) |
---|
617 | { |
---|
618 | long n = deg(f); |
---|
619 | |
---|
620 | if (n <= 0) return 0; |
---|
621 | if (n == 1) return 1; |
---|
622 | |
---|
623 | const ZZ& p = ZZ_p::modulus(); |
---|
624 | |
---|
625 | ZZ_pXModulus F; |
---|
626 | |
---|
627 | build(F, f); |
---|
628 | |
---|
629 | ZZ_pX b, r, s; |
---|
630 | |
---|
631 | PowerXMod(b, p, F); |
---|
632 | |
---|
633 | long i; |
---|
634 | |
---|
635 | for (i = 0; i < iter; i++) { |
---|
636 | random(r, n); |
---|
637 | TraceMap(s, r, n, F, b); |
---|
638 | |
---|
639 | if (deg(s) > 0) return 0; |
---|
640 | } |
---|
641 | |
---|
642 | if (p >= n) return 1; |
---|
643 | |
---|
644 | long pp; |
---|
645 | |
---|
646 | conv(pp, p); |
---|
647 | |
---|
648 | if (n % pp != 0) return 1; |
---|
649 | |
---|
650 | PowerCompose(s, b, n/pp, F); |
---|
651 | return !IsX(s); |
---|
652 | } |
---|
653 | |
---|
654 | long ZZ_pX_BlockingFactor = 10; |
---|
655 | |
---|
656 | void DDF(vec_pair_ZZ_pX_long& factors, const ZZ_pX& ff, const ZZ_pX& hh, |
---|
657 | long verbose) |
---|
658 | { |
---|
659 | ZZ_pX f = ff; |
---|
660 | ZZ_pX h = hh; |
---|
661 | |
---|
662 | if (!IsOne(LeadCoeff(f))) |
---|
663 | Error("DDF: bad args"); |
---|
664 | |
---|
665 | factors.SetLength(0); |
---|
666 | |
---|
667 | if (deg(f) == 0) |
---|
668 | return; |
---|
669 | |
---|
670 | if (deg(f) == 1) { |
---|
671 | AddFactor(factors, f, 1, verbose); |
---|
672 | return; |
---|
673 | } |
---|
674 | |
---|
675 | long CompTableSize = 2*SqrRoot(deg(f)); |
---|
676 | |
---|
677 | long GCDTableSize = ZZ_pX_BlockingFactor; |
---|
678 | |
---|
679 | ZZ_pXModulus F; |
---|
680 | build(F, f); |
---|
681 | |
---|
682 | ZZ_pXArgument H; |
---|
683 | |
---|
684 | build(H, h, F, min(CompTableSize, deg(f))); |
---|
685 | |
---|
686 | long i, d, limit, old_n; |
---|
687 | ZZ_pX g, X; |
---|
688 | |
---|
689 | |
---|
690 | vec_ZZ_pX tbl(INIT_SIZE, GCDTableSize); |
---|
691 | |
---|
692 | SetX(X); |
---|
693 | |
---|
694 | i = 0; |
---|
695 | g = h; |
---|
696 | d = 1; |
---|
697 | limit = GCDTableSize; |
---|
698 | |
---|
699 | |
---|
700 | while (2*d <= deg(f)) { |
---|
701 | |
---|
702 | old_n = deg(f); |
---|
703 | sub(tbl[i], g, X); |
---|
704 | i++; |
---|
705 | if (i == limit) { |
---|
706 | ProcessTable(f, factors, F, i, tbl, d, verbose); |
---|
707 | i = 0; |
---|
708 | } |
---|
709 | |
---|
710 | d = d + 1; |
---|
711 | if (2*d <= deg(f)) { |
---|
712 | // we need to go further |
---|
713 | |
---|
714 | if (deg(f) < old_n) { |
---|
715 | // f has changed |
---|
716 | |
---|
717 | build(F, f); |
---|
718 | rem(h, h, f); |
---|
719 | rem(g, g, f); |
---|
720 | build(H, h, F, min(CompTableSize, deg(f))); |
---|
721 | } |
---|
722 | |
---|
723 | CompMod(g, g, H, F); |
---|
724 | } |
---|
725 | } |
---|
726 | |
---|
727 | ProcessTable(f, factors, F, i, tbl, d-1, verbose); |
---|
728 | |
---|
729 | if (!IsOne(f)) AddFactor(factors, f, deg(f), verbose); |
---|
730 | } |
---|
731 | |
---|
732 | |
---|
733 | |
---|
734 | void RootEDF(vec_ZZ_pX& factors, const ZZ_pX& f, long verbose) |
---|
735 | { |
---|
736 | vec_ZZ_p roots; |
---|
737 | double t; |
---|
738 | |
---|
739 | FindRoots(roots, f); |
---|
740 | |
---|
741 | long r = roots.length(); |
---|
742 | factors.SetLength(r); |
---|
743 | for (long j = 0; j < r; j++) { |
---|
744 | SetX(factors[j]); |
---|
745 | sub(factors[j], factors[j], roots[j]); |
---|
746 | } |
---|
747 | } |
---|
748 | |
---|
749 | static |
---|
750 | void EDFSplit(vec_ZZ_pX& v, const ZZ_pX& f, const ZZ_pX& b, long d) |
---|
751 | { |
---|
752 | ZZ_pX a, g, h; |
---|
753 | ZZ_pXModulus F; |
---|
754 | vec_ZZ_p roots; |
---|
755 | |
---|
756 | build(F, f); |
---|
757 | long n = F.n; |
---|
758 | long r = n/d; |
---|
759 | random(a, n); |
---|
760 | TraceMap(g, a, d, F, b); |
---|
761 | MinPolyMod(h, g, F, r); |
---|
762 | FindRoots(roots, h); |
---|
763 | FindFactors(v, f, g, roots); |
---|
764 | } |
---|
765 | |
---|
766 | static |
---|
767 | void RecEDF(vec_ZZ_pX& factors, const ZZ_pX& f, const ZZ_pX& b, long d, |
---|
768 | long verbose) |
---|
769 | { |
---|
770 | vec_ZZ_pX v; |
---|
771 | long i; |
---|
772 | ZZ_pX bb; |
---|
773 | |
---|
774 | EDFSplit(v, f, b, d); |
---|
775 | for (i = 0; i < v.length(); i++) { |
---|
776 | if (deg(v[i]) == d) { |
---|
777 | append(factors, v[i]); |
---|
778 | } |
---|
779 | else { |
---|
780 | ZZ_pX bb; |
---|
781 | rem(bb, b, v[i]); |
---|
782 | RecEDF(factors, v[i], bb, d, verbose); |
---|
783 | } |
---|
784 | } |
---|
785 | } |
---|
786 | |
---|
787 | |
---|
788 | void EDF(vec_ZZ_pX& factors, const ZZ_pX& ff, const ZZ_pX& bb, |
---|
789 | long d, long verbose) |
---|
790 | |
---|
791 | { |
---|
792 | ZZ_pX f = ff; |
---|
793 | ZZ_pX b = bb; |
---|
794 | |
---|
795 | if (!IsOne(LeadCoeff(f))) |
---|
796 | Error("EDF: bad args"); |
---|
797 | |
---|
798 | long n = deg(f); |
---|
799 | long r = n/d; |
---|
800 | |
---|
801 | if (r == 0) { |
---|
802 | factors.SetLength(0); |
---|
803 | return; |
---|
804 | } |
---|
805 | |
---|
806 | if (r == 1) { |
---|
807 | factors.SetLength(1); |
---|
808 | factors[0] = f; |
---|
809 | return; |
---|
810 | } |
---|
811 | |
---|
812 | if (d == 1) { |
---|
813 | RootEDF(factors, f, verbose); |
---|
814 | return; |
---|
815 | } |
---|
816 | |
---|
817 | factors.SetLength(0); |
---|
818 | |
---|
819 | RecEDF(factors, f, b, d, verbose); |
---|
820 | |
---|
821 | } |
---|
822 | |
---|
823 | |
---|
824 | void SFCanZass(vec_ZZ_pX& factors, const ZZ_pX& ff, long verbose) |
---|
825 | { |
---|
826 | ZZ_pX f = ff; |
---|
827 | |
---|
828 | if (!IsOne(LeadCoeff(f))) |
---|
829 | Error("SFCanZass: bad args"); |
---|
830 | |
---|
831 | if (deg(f) == 0) { |
---|
832 | factors.SetLength(0); |
---|
833 | return; |
---|
834 | } |
---|
835 | |
---|
836 | if (deg(f) == 1) { |
---|
837 | factors.SetLength(1); |
---|
838 | factors[0] = f; |
---|
839 | return; |
---|
840 | } |
---|
841 | |
---|
842 | factors.SetLength(0); |
---|
843 | |
---|
844 | double t; |
---|
845 | |
---|
846 | const ZZ& p = ZZ_p::modulus(); |
---|
847 | |
---|
848 | |
---|
849 | ZZ_pXModulus F; |
---|
850 | build(F, f); |
---|
851 | |
---|
852 | ZZ_pX h; |
---|
853 | |
---|
854 | PowerXMod(h, p, F); |
---|
855 | |
---|
856 | vec_pair_ZZ_pX_long u; |
---|
857 | NewDDF(u, f, h, verbose); |
---|
858 | |
---|
859 | ZZ_pX hh; |
---|
860 | vec_ZZ_pX v; |
---|
861 | |
---|
862 | long i; |
---|
863 | for (i = 0; i < u.length(); i++) { |
---|
864 | const ZZ_pX& g = u[i].a; |
---|
865 | long d = u[i].b; |
---|
866 | long r = deg(g)/d; |
---|
867 | |
---|
868 | if (r == 1) { |
---|
869 | // g is already irreducible |
---|
870 | |
---|
871 | append(factors, g); |
---|
872 | } |
---|
873 | else { |
---|
874 | // must perform EDF |
---|
875 | |
---|
876 | if (d == 1) { |
---|
877 | // root finding |
---|
878 | RootEDF(v, g, verbose); |
---|
879 | append(factors, v); |
---|
880 | } |
---|
881 | else { |
---|
882 | // general case |
---|
883 | rem(hh, h, g); |
---|
884 | EDF(v, g, hh, d, verbose); |
---|
885 | append(factors, v); |
---|
886 | } |
---|
887 | } |
---|
888 | } |
---|
889 | } |
---|
890 | |
---|
891 | void CanZass(vec_pair_ZZ_pX_long& factors, const ZZ_pX& f, long verbose) |
---|
892 | { |
---|
893 | if (!IsOne(LeadCoeff(f))) |
---|
894 | Error("CanZass: bad args"); |
---|
895 | |
---|
896 | double t; |
---|
897 | vec_pair_ZZ_pX_long sfd; |
---|
898 | vec_ZZ_pX x; |
---|
899 | |
---|
900 | |
---|
901 | SquareFreeDecomp(sfd, f); |
---|
902 | |
---|
903 | factors.SetLength(0); |
---|
904 | |
---|
905 | long i, j; |
---|
906 | |
---|
907 | for (i = 0; i < sfd.length(); i++) { |
---|
908 | |
---|
909 | SFCanZass(x, sfd[i].a, verbose); |
---|
910 | |
---|
911 | for (j = 0; j < x.length(); j++) |
---|
912 | append(factors, cons(x[j], sfd[i].b)); |
---|
913 | } |
---|
914 | } |
---|
915 | |
---|
916 | void mul(ZZ_pX& f, const vec_pair_ZZ_pX_long& v) |
---|
917 | { |
---|
918 | long i, j, n; |
---|
919 | |
---|
920 | n = 0; |
---|
921 | for (i = 0; i < v.length(); i++) |
---|
922 | n += v[i].b*deg(v[i].a); |
---|
923 | |
---|
924 | ZZ_pX g(INIT_SIZE, n+1); |
---|
925 | |
---|
926 | set(g); |
---|
927 | for (i = 0; i < v.length(); i++) |
---|
928 | for (j = 0; j < v[i].b; j++) { |
---|
929 | mul(g, g, v[i].a); |
---|
930 | } |
---|
931 | |
---|
932 | f = g; |
---|
933 | } |
---|
934 | |
---|
935 | |
---|
936 | |
---|
937 | |
---|
938 | static |
---|
939 | long BaseCase(const ZZ_pX& h, long q, long a, const ZZ_pXModulus& F) |
---|
940 | { |
---|
941 | long b, e; |
---|
942 | ZZ_pX lh(INIT_SIZE, F.n); |
---|
943 | |
---|
944 | lh = h; |
---|
945 | b = 1; |
---|
946 | e = 0; |
---|
947 | while (e < a-1 && !IsX(lh)) { |
---|
948 | e++; |
---|
949 | b *= q; |
---|
950 | PowerCompose(lh, lh, q, F); |
---|
951 | } |
---|
952 | |
---|
953 | if (!IsX(lh)) b *= q; |
---|
954 | |
---|
955 | return b; |
---|
956 | } |
---|
957 | |
---|
958 | |
---|
959 | |
---|
960 | static |
---|
961 | void TandemPowerCompose(ZZ_pX& y1, ZZ_pX& y2, const ZZ_pX& h, |
---|
962 | long q1, long q2, const ZZ_pXModulus& F) |
---|
963 | { |
---|
964 | ZZ_pX z(INIT_SIZE, F.n); |
---|
965 | long sw; |
---|
966 | |
---|
967 | z = h; |
---|
968 | SetX(y1); |
---|
969 | SetX(y2); |
---|
970 | |
---|
971 | while (q1 || q2) { |
---|
972 | sw = 0; |
---|
973 | |
---|
974 | if (q1 > 1 || q2 > 1) sw = 4; |
---|
975 | |
---|
976 | if (q1 & 1) { |
---|
977 | if (IsX(y1)) |
---|
978 | y1 = z; |
---|
979 | else |
---|
980 | sw = sw | 2; |
---|
981 | } |
---|
982 | |
---|
983 | if (q2 & 1) { |
---|
984 | if (IsX(y2)) |
---|
985 | y2 = z; |
---|
986 | else |
---|
987 | sw = sw | 1; |
---|
988 | } |
---|
989 | |
---|
990 | switch (sw) { |
---|
991 | case 0: |
---|
992 | break; |
---|
993 | |
---|
994 | case 1: |
---|
995 | CompMod(y2, y2, z, F); |
---|
996 | break; |
---|
997 | |
---|
998 | case 2: |
---|
999 | CompMod(y1, y1, z, F); |
---|
1000 | break; |
---|
1001 | |
---|
1002 | case 3: |
---|
1003 | Comp2Mod(y1, y2, y1, y2, z, F); |
---|
1004 | break; |
---|
1005 | |
---|
1006 | case 4: |
---|
1007 | CompMod(z, z, z, F); |
---|
1008 | break; |
---|
1009 | |
---|
1010 | case 5: |
---|
1011 | Comp2Mod(z, y2, z, y2, z, F); |
---|
1012 | break; |
---|
1013 | |
---|
1014 | case 6: |
---|
1015 | Comp2Mod(z, y1, z, y1, z, F); |
---|
1016 | break; |
---|
1017 | |
---|
1018 | case 7: |
---|
1019 | Comp3Mod(z, y1, y2, z, y1, y2, z, F); |
---|
1020 | break; |
---|
1021 | } |
---|
1022 | |
---|
1023 | q1 = q1 >> 1; |
---|
1024 | q2 = q2 >> 1; |
---|
1025 | } |
---|
1026 | } |
---|
1027 | |
---|
1028 | |
---|
1029 | |
---|
1030 | static |
---|
1031 | long RecComputeDegree(long u, const ZZ_pX& h, const ZZ_pXModulus& F, |
---|
1032 | FacVec& fvec) |
---|
1033 | { |
---|
1034 | if (IsX(h)) return 1; |
---|
1035 | |
---|
1036 | if (fvec[u].link == -1) return BaseCase(h, fvec[u].q, fvec[u].a, F); |
---|
1037 | |
---|
1038 | ZZ_pX h1, h2; |
---|
1039 | long q1, q2, r1, r2; |
---|
1040 | |
---|
1041 | q1 = fvec[fvec[u].link].val; |
---|
1042 | q2 = fvec[fvec[u].link+1].val; |
---|
1043 | |
---|
1044 | TandemPowerCompose(h1, h2, h, q1, q2, F); |
---|
1045 | r1 = RecComputeDegree(fvec[u].link, h2, F, fvec); |
---|
1046 | r2 = RecComputeDegree(fvec[u].link+1, h1, F, fvec); |
---|
1047 | return r1*r2; |
---|
1048 | } |
---|
1049 | |
---|
1050 | |
---|
1051 | |
---|
1052 | |
---|
1053 | long ComputeDegree(const ZZ_pX& h, const ZZ_pXModulus& F) |
---|
1054 | // f = F.f is assumed to be an "equal degree" polynomial |
---|
1055 | // h = X^p mod f |
---|
1056 | // the common degree of the irreducible factors of f is computed |
---|
1057 | { |
---|
1058 | if (F.n == 1 || IsX(h)) return 1; |
---|
1059 | |
---|
1060 | FacVec fvec; |
---|
1061 | |
---|
1062 | FactorInt(fvec, F.n); |
---|
1063 | |
---|
1064 | return RecComputeDegree(fvec.length()-1, h, F, fvec); |
---|
1065 | } |
---|
1066 | |
---|
1067 | long ProbComputeDegree(const ZZ_pX& h, const ZZ_pXModulus& F) |
---|
1068 | { |
---|
1069 | if (F.n == 1 || IsX(h)) |
---|
1070 | return 1; |
---|
1071 | |
---|
1072 | long n = F.n; |
---|
1073 | |
---|
1074 | ZZ_pX P1, P2, P3; |
---|
1075 | |
---|
1076 | random(P1, n); |
---|
1077 | TraceMap(P2, P1, n, F, h); |
---|
1078 | ProbMinPolyMod(P3, P2, F, n/2); |
---|
1079 | |
---|
1080 | long r = deg(P3); |
---|
1081 | |
---|
1082 | if (r <= 0 || n % r != 0) |
---|
1083 | return 0; |
---|
1084 | else |
---|
1085 | return n/r; |
---|
1086 | } |
---|
1087 | |
---|
1088 | |
---|
1089 | |
---|
1090 | void FindRoot(ZZ_p& root, const ZZ_pX& ff) |
---|
1091 | // finds a root of ff. |
---|
1092 | // assumes that ff is monic and splits into distinct linear factors |
---|
1093 | |
---|
1094 | { |
---|
1095 | ZZ_pXModulus F; |
---|
1096 | ZZ_pX h, h1, f; |
---|
1097 | ZZ_p r; |
---|
1098 | ZZ p1; |
---|
1099 | |
---|
1100 | f = ff; |
---|
1101 | |
---|
1102 | if (!IsOne(LeadCoeff(f))) |
---|
1103 | Error("FindRoot: bad args"); |
---|
1104 | |
---|
1105 | if (deg(f) == 0) |
---|
1106 | Error("FindRoot: bad args"); |
---|
1107 | |
---|
1108 | RightShift(p1, ZZ_p::modulus(), 1); |
---|
1109 | h1 = 1; |
---|
1110 | |
---|
1111 | while (deg(f) > 1) { |
---|
1112 | build(F, f); |
---|
1113 | random(r); |
---|
1114 | PowerXPlusAMod(h, r, p1, F); |
---|
1115 | sub(h, h, h1); |
---|
1116 | GCD(h, h, f); |
---|
1117 | if (deg(h) > 0 && deg(h) < deg(f)) { |
---|
1118 | if (deg(h) > deg(f)/2) |
---|
1119 | div(f, f, h); |
---|
1120 | else |
---|
1121 | f = h; |
---|
1122 | } |
---|
1123 | } |
---|
1124 | |
---|
1125 | negate(root, ConstTerm(f)); |
---|
1126 | } |
---|
1127 | |
---|
1128 | |
---|
1129 | static |
---|
1130 | long power(long a, long e) |
---|
1131 | { |
---|
1132 | long i, res; |
---|
1133 | |
---|
1134 | res = 1; |
---|
1135 | for (i = 1; i <= e; i++) |
---|
1136 | res = res * a; |
---|
1137 | |
---|
1138 | return res; |
---|
1139 | } |
---|
1140 | |
---|
1141 | |
---|
1142 | static |
---|
1143 | long IrredBaseCase(const ZZ_pX& h, long q, long a, const ZZ_pXModulus& F) |
---|
1144 | { |
---|
1145 | long e; |
---|
1146 | ZZ_pX X, s, d; |
---|
1147 | |
---|
1148 | e = power(q, a-1); |
---|
1149 | PowerCompose(s, h, e, F); |
---|
1150 | SetX(X); |
---|
1151 | sub(s, s, X); |
---|
1152 | GCD(d, F.f, s); |
---|
1153 | return IsOne(d); |
---|
1154 | } |
---|
1155 | |
---|
1156 | |
---|
1157 | static |
---|
1158 | long RecIrredTest(long u, const ZZ_pX& h, const ZZ_pXModulus& F, |
---|
1159 | const FacVec& fvec) |
---|
1160 | { |
---|
1161 | long q1, q2; |
---|
1162 | ZZ_pX h1, h2; |
---|
1163 | |
---|
1164 | if (IsX(h)) return 0; |
---|
1165 | |
---|
1166 | if (fvec[u].link == -1) { |
---|
1167 | return IrredBaseCase(h, fvec[u].q, fvec[u].a, F); |
---|
1168 | } |
---|
1169 | |
---|
1170 | |
---|
1171 | q1 = fvec[fvec[u].link].val; |
---|
1172 | q2 = fvec[fvec[u].link+1].val; |
---|
1173 | |
---|
1174 | TandemPowerCompose(h1, h2, h, q1, q2, F); |
---|
1175 | return RecIrredTest(fvec[u].link, h2, F, fvec) |
---|
1176 | && RecIrredTest(fvec[u].link+1, h1, F, fvec); |
---|
1177 | } |
---|
1178 | |
---|
1179 | long DetIrredTest(const ZZ_pX& f) |
---|
1180 | { |
---|
1181 | if (deg(f) <= 0) return 0; |
---|
1182 | if (deg(f) == 1) return 1; |
---|
1183 | |
---|
1184 | ZZ_pXModulus F; |
---|
1185 | |
---|
1186 | build(F, f); |
---|
1187 | |
---|
1188 | ZZ_pX h; |
---|
1189 | |
---|
1190 | PowerXMod(h, ZZ_p::modulus(), F); |
---|
1191 | |
---|
1192 | ZZ_pX s; |
---|
1193 | PowerCompose(s, h, F.n, F); |
---|
1194 | if (!IsX(s)) return 0; |
---|
1195 | |
---|
1196 | FacVec fvec; |
---|
1197 | |
---|
1198 | FactorInt(fvec, F.n); |
---|
1199 | |
---|
1200 | return RecIrredTest(fvec.length()-1, h, F, fvec); |
---|
1201 | } |
---|
1202 | |
---|
1203 | |
---|
1204 | |
---|
1205 | long IterIrredTest(const ZZ_pX& f) |
---|
1206 | { |
---|
1207 | if (deg(f) <= 0) return 0; |
---|
1208 | if (deg(f) == 1) return 1; |
---|
1209 | |
---|
1210 | ZZ_pXModulus F; |
---|
1211 | |
---|
1212 | build(F, f); |
---|
1213 | |
---|
1214 | ZZ_pX h; |
---|
1215 | |
---|
1216 | PowerXMod(h, ZZ_p::modulus(), F); |
---|
1217 | |
---|
1218 | long CompTableSize = 2*SqrRoot(deg(f)); |
---|
1219 | |
---|
1220 | ZZ_pXArgument H; |
---|
1221 | |
---|
1222 | build(H, h, F, CompTableSize); |
---|
1223 | |
---|
1224 | long i, d, limit, limit_sqr; |
---|
1225 | ZZ_pX g, X, t, prod; |
---|
1226 | |
---|
1227 | |
---|
1228 | SetX(X); |
---|
1229 | |
---|
1230 | i = 0; |
---|
1231 | g = h; |
---|
1232 | d = 1; |
---|
1233 | limit = 2; |
---|
1234 | limit_sqr = limit*limit; |
---|
1235 | |
---|
1236 | set(prod); |
---|
1237 | |
---|
1238 | |
---|
1239 | while (2*d <= deg(f)) { |
---|
1240 | sub(t, g, X); |
---|
1241 | MulMod(prod, prod, t, F); |
---|
1242 | i++; |
---|
1243 | if (i == limit_sqr) { |
---|
1244 | GCD(t, f, prod); |
---|
1245 | if (!IsOne(t)) return 0; |
---|
1246 | |
---|
1247 | set(prod); |
---|
1248 | limit++; |
---|
1249 | limit_sqr = limit*limit; |
---|
1250 | i = 0; |
---|
1251 | } |
---|
1252 | |
---|
1253 | d = d + 1; |
---|
1254 | if (2*d <= deg(f)) { |
---|
1255 | CompMod(g, g, H, F); |
---|
1256 | } |
---|
1257 | } |
---|
1258 | |
---|
1259 | if (i > 0) { |
---|
1260 | GCD(t, f, prod); |
---|
1261 | if (!IsOne(t)) return 0; |
---|
1262 | } |
---|
1263 | |
---|
1264 | return 1; |
---|
1265 | } |
---|
1266 | |
---|
1267 | |
---|
1268 | static |
---|
1269 | void MulByXPlusY(vec_ZZ_pX& h, const ZZ_pX& f, const ZZ_pX& g) |
---|
1270 | // h represents the bivariate polynomial h[0] + h[1]*Y + ... + h[n-1]*Y^k, |
---|
1271 | // where the h[i]'s are polynomials in X, each of degree < deg(f), |
---|
1272 | // and k < deg(g). |
---|
1273 | // h is replaced by the bivariate polynomial h*(X+Y) (mod f(X), g(Y)). |
---|
1274 | |
---|
1275 | { |
---|
1276 | long n = deg(g); |
---|
1277 | long k = h.length()-1; |
---|
1278 | |
---|
1279 | if (k < 0) return; |
---|
1280 | |
---|
1281 | if (k < n-1) { |
---|
1282 | h.SetLength(k+2); |
---|
1283 | h[k+1] = h[k]; |
---|
1284 | for (long i = k; i >= 1; i--) { |
---|
1285 | MulByXMod(h[i], h[i], f); |
---|
1286 | add(h[i], h[i], h[i-1]); |
---|
1287 | } |
---|
1288 | MulByXMod(h[0], h[0], f); |
---|
1289 | } |
---|
1290 | else { |
---|
1291 | ZZ_pX b, t; |
---|
1292 | |
---|
1293 | b = h[n-1]; |
---|
1294 | for (long i = n-1; i >= 1; i--) { |
---|
1295 | mul(t, b, g.rep[i]); |
---|
1296 | MulByXMod(h[i], h[i], f); |
---|
1297 | add(h[i], h[i], h[i-1]); |
---|
1298 | sub(h[i], h[i], t); |
---|
1299 | } |
---|
1300 | mul(t, b, g.rep[0]); |
---|
1301 | MulByXMod(h[0], h[0], f); |
---|
1302 | sub(h[0], h[0], t); |
---|
1303 | } |
---|
1304 | |
---|
1305 | // normalize |
---|
1306 | |
---|
1307 | k = h.length()-1; |
---|
1308 | while (k >= 0 && IsZero(h[k])) k--; |
---|
1309 | h.SetLength(k+1); |
---|
1310 | } |
---|
1311 | |
---|
1312 | |
---|
1313 | static |
---|
1314 | void IrredCombine(ZZ_pX& x, const ZZ_pX& f, const ZZ_pX& g) |
---|
1315 | { |
---|
1316 | if (deg(f) < deg(g)) { |
---|
1317 | IrredCombine(x, g, f); |
---|
1318 | return; |
---|
1319 | } |
---|
1320 | |
---|
1321 | // deg(f) >= deg(g)...not necessary, but maybe a little more |
---|
1322 | // time & space efficient |
---|
1323 | |
---|
1324 | long df = deg(f); |
---|
1325 | long dg = deg(g); |
---|
1326 | long m = df*dg; |
---|
1327 | |
---|
1328 | vec_ZZ_pX h(INIT_SIZE, dg); |
---|
1329 | |
---|
1330 | long i; |
---|
1331 | for (i = 0; i < dg; i++) h[i].SetMaxLength(df); |
---|
1332 | |
---|
1333 | h.SetLength(1); |
---|
1334 | set(h[0]); |
---|
1335 | |
---|
1336 | vec_ZZ_p a; |
---|
1337 | |
---|
1338 | a.SetLength(2*m); |
---|
1339 | |
---|
1340 | for (i = 0; i < 2*m; i++) { |
---|
1341 | a[i] = ConstTerm(h[0]); |
---|
1342 | if (i < 2*m-1) |
---|
1343 | MulByXPlusY(h, f, g); |
---|
1344 | } |
---|
1345 | |
---|
1346 | MinPolySeq(x, a, m); |
---|
1347 | } |
---|
1348 | |
---|
1349 | |
---|
1350 | static |
---|
1351 | void BuildPrimePowerIrred(ZZ_pX& f, long q, long e) |
---|
1352 | { |
---|
1353 | long n = power(q, e); |
---|
1354 | |
---|
1355 | do { |
---|
1356 | random(f, n); |
---|
1357 | SetCoeff(f, n); |
---|
1358 | } while (!IterIrredTest(f)); |
---|
1359 | } |
---|
1360 | |
---|
1361 | static |
---|
1362 | void RecBuildIrred(ZZ_pX& f, long u, const FacVec& fvec) |
---|
1363 | { |
---|
1364 | if (fvec[u].link == -1) |
---|
1365 | BuildPrimePowerIrred(f, fvec[u].q, fvec[u].a); |
---|
1366 | else { |
---|
1367 | ZZ_pX g, h; |
---|
1368 | RecBuildIrred(g, fvec[u].link, fvec); |
---|
1369 | RecBuildIrred(h, fvec[u].link+1, fvec); |
---|
1370 | IrredCombine(f, g, h); |
---|
1371 | } |
---|
1372 | } |
---|
1373 | |
---|
1374 | |
---|
1375 | void BuildIrred(ZZ_pX& f, long n) |
---|
1376 | { |
---|
1377 | if (n <= 0) |
---|
1378 | Error("BuildIrred: n must be positive"); |
---|
1379 | |
---|
1380 | if (NTL_OVERFLOW(n, 1, 0)) Error("overflow in BuildIrred"); |
---|
1381 | |
---|
1382 | if (n == 1) { |
---|
1383 | SetX(f); |
---|
1384 | return; |
---|
1385 | } |
---|
1386 | |
---|
1387 | FacVec fvec; |
---|
1388 | |
---|
1389 | FactorInt(fvec, n); |
---|
1390 | |
---|
1391 | RecBuildIrred(f, fvec.length()-1, fvec); |
---|
1392 | } |
---|
1393 | |
---|
1394 | |
---|
1395 | |
---|
1396 | void BuildRandomIrred(ZZ_pX& f, const ZZ_pX& g) |
---|
1397 | { |
---|
1398 | ZZ_pXModulus G; |
---|
1399 | ZZ_pX h, ff; |
---|
1400 | |
---|
1401 | build(G, g); |
---|
1402 | do { |
---|
1403 | random(h, deg(g)); |
---|
1404 | IrredPolyMod(ff, h, G); |
---|
1405 | } while (deg(ff) < deg(g)); |
---|
1406 | |
---|
1407 | f = ff; |
---|
1408 | } |
---|
1409 | |
---|
1410 | |
---|
1411 | /************* NEW DDF ****************/ |
---|
1412 | |
---|
1413 | long ZZ_pX_GCDTableSize = 4; |
---|
1414 | char ZZ_pX_stem[256] = ""; |
---|
1415 | |
---|
1416 | double ZZ_pXFileThresh = 256; |
---|
1417 | |
---|
1418 | static vec_ZZ_pX BabyStepFile; |
---|
1419 | static vec_ZZ_pX GiantStepFile; |
---|
1420 | |
---|
1421 | |
---|
1422 | |
---|
1423 | static |
---|
1424 | double CalcTableSize(long n, long k) |
---|
1425 | { |
---|
1426 | double sz = ZZ_p::storage(); |
---|
1427 | sz = sz * n; |
---|
1428 | sz = sz + NTL_VECTOR_HEADER_SIZE + sizeof(vec_ZZ_p); |
---|
1429 | sz = sz * k; |
---|
1430 | sz = sz/1024; |
---|
1431 | return sz; |
---|
1432 | } |
---|
1433 | |
---|
1434 | |
---|
1435 | static |
---|
1436 | void GenerateBabySteps(ZZ_pX& h1, const ZZ_pX& f, const ZZ_pX& h, long k, |
---|
1437 | long verbose) |
---|
1438 | |
---|
1439 | { |
---|
1440 | double t; |
---|
1441 | |
---|
1442 | ZZ_pXModulus F; |
---|
1443 | build(F, f); |
---|
1444 | |
---|
1445 | ZZ_pXArgument H; |
---|
1446 | build(H, h, F, 2*SqrRoot(F.n)); |
---|
1447 | |
---|
1448 | |
---|
1449 | h1 = h; |
---|
1450 | |
---|
1451 | long i; |
---|
1452 | |
---|
1453 | BabyStepFile.kill(); |
---|
1454 | BabyStepFile.SetLength(k-1); |
---|
1455 | |
---|
1456 | for (i = 1; i <= k-1; i++) { |
---|
1457 | BabyStepFile(i) = h1; |
---|
1458 | |
---|
1459 | CompMod(h1, h1, H, F); |
---|
1460 | } |
---|
1461 | |
---|
1462 | } |
---|
1463 | |
---|
1464 | |
---|
1465 | static |
---|
1466 | void GenerateGiantSteps(const ZZ_pX& f, const ZZ_pX& h, long l, long verbose) |
---|
1467 | { |
---|
1468 | |
---|
1469 | double t; |
---|
1470 | |
---|
1471 | |
---|
1472 | ZZ_pXModulus F; |
---|
1473 | build(F, f); |
---|
1474 | |
---|
1475 | ZZ_pXArgument H; |
---|
1476 | build(H, h, F, 2*SqrRoot(F.n)); |
---|
1477 | |
---|
1478 | ZZ_pX h1; |
---|
1479 | |
---|
1480 | h1 = h; |
---|
1481 | |
---|
1482 | long i; |
---|
1483 | |
---|
1484 | GiantStepFile.kill(); |
---|
1485 | GiantStepFile.SetLength(l); |
---|
1486 | |
---|
1487 | for (i = 1; i <= l-1; i++) { |
---|
1488 | GiantStepFile(i) = h1; |
---|
1489 | |
---|
1490 | CompMod(h1, h1, H, F); |
---|
1491 | } |
---|
1492 | |
---|
1493 | GiantStepFile(i) = h1; |
---|
1494 | |
---|
1495 | } |
---|
1496 | |
---|
1497 | static |
---|
1498 | void FileCleanup(long k, long l) |
---|
1499 | { |
---|
1500 | BabyStepFile.kill(); |
---|
1501 | GiantStepFile.kill(); |
---|
1502 | } |
---|
1503 | |
---|
1504 | |
---|
1505 | static |
---|
1506 | void NewAddFactor(vec_pair_ZZ_pX_long& u, const ZZ_pX& g, long m, long verbose) |
---|
1507 | { |
---|
1508 | long len = u.length(); |
---|
1509 | |
---|
1510 | u.SetLength(len+1); |
---|
1511 | u[len].a = g; |
---|
1512 | u[len].b = m; |
---|
1513 | |
---|
1514 | } |
---|
1515 | |
---|
1516 | |
---|
1517 | |
---|
1518 | |
---|
1519 | static |
---|
1520 | void NewProcessTable(vec_pair_ZZ_pX_long& u, ZZ_pX& f, const ZZ_pXModulus& F, |
---|
1521 | vec_ZZ_pX& buf, long size, long StartInterval, |
---|
1522 | long IntervalLength, long verbose) |
---|
1523 | |
---|
1524 | { |
---|
1525 | if (size == 0) return; |
---|
1526 | |
---|
1527 | ZZ_pX& g = buf[size-1]; |
---|
1528 | |
---|
1529 | long i; |
---|
1530 | |
---|
1531 | for (i = 0; i < size-1; i++) |
---|
1532 | MulMod(g, g, buf[i], F); |
---|
1533 | |
---|
1534 | GCD(g, f, g); |
---|
1535 | |
---|
1536 | if (deg(g) == 0) return; |
---|
1537 | |
---|
1538 | div(f, f, g); |
---|
1539 | |
---|
1540 | long d = (StartInterval-1)*IntervalLength + 1; |
---|
1541 | i = 0; |
---|
1542 | long interval = StartInterval; |
---|
1543 | |
---|
1544 | while (i < size-1 && 2*d <= deg(g)) { |
---|
1545 | GCD(buf[i], buf[i], g); |
---|
1546 | if (deg(buf[i]) > 0) { |
---|
1547 | NewAddFactor(u, buf[i], interval, verbose); |
---|
1548 | div(g, g, buf[i]); |
---|
1549 | } |
---|
1550 | |
---|
1551 | i++; |
---|
1552 | interval++; |
---|
1553 | d += IntervalLength; |
---|
1554 | } |
---|
1555 | |
---|
1556 | if (deg(g) > 0) { |
---|
1557 | if (i == size-1) |
---|
1558 | NewAddFactor(u, g, interval, verbose); |
---|
1559 | else |
---|
1560 | NewAddFactor(u, g, (deg(g)+IntervalLength-1)/IntervalLength, verbose); |
---|
1561 | } |
---|
1562 | } |
---|
1563 | |
---|
1564 | |
---|
1565 | static |
---|
1566 | void FetchGiantStep(ZZ_pX& g, long gs, const ZZ_pXModulus& F) |
---|
1567 | { |
---|
1568 | g = GiantStepFile(gs); |
---|
1569 | |
---|
1570 | rem(g, g, F); |
---|
1571 | } |
---|
1572 | |
---|
1573 | |
---|
1574 | static |
---|
1575 | void FetchBabySteps(vec_ZZ_pX& v, long k) |
---|
1576 | { |
---|
1577 | v.SetLength(k); |
---|
1578 | |
---|
1579 | SetX(v[0]); |
---|
1580 | |
---|
1581 | long i; |
---|
1582 | for (i = 1; i <= k-1; i++) { |
---|
1583 | v[i] = BabyStepFile(i); |
---|
1584 | } |
---|
1585 | } |
---|
1586 | |
---|
1587 | |
---|
1588 | |
---|
1589 | static |
---|
1590 | void GiantRefine(vec_pair_ZZ_pX_long& u, const ZZ_pX& ff, long k, long l, |
---|
1591 | long verbose) |
---|
1592 | |
---|
1593 | { |
---|
1594 | u.SetLength(0); |
---|
1595 | |
---|
1596 | vec_ZZ_pX BabyStep; |
---|
1597 | |
---|
1598 | FetchBabySteps(BabyStep, k); |
---|
1599 | |
---|
1600 | vec_ZZ_pX buf(INIT_SIZE, ZZ_pX_GCDTableSize); |
---|
1601 | |
---|
1602 | ZZ_pX f; |
---|
1603 | f = ff; |
---|
1604 | |
---|
1605 | ZZ_pXModulus F; |
---|
1606 | build(F, f); |
---|
1607 | |
---|
1608 | ZZ_pX g; |
---|
1609 | ZZ_pX h; |
---|
1610 | |
---|
1611 | long size = 0; |
---|
1612 | |
---|
1613 | long first_gs; |
---|
1614 | |
---|
1615 | long d = 1; |
---|
1616 | |
---|
1617 | while (2*d <= deg(f)) { |
---|
1618 | |
---|
1619 | long old_n = deg(f); |
---|
1620 | |
---|
1621 | long gs = (d+k-1)/k; |
---|
1622 | long bs = gs*k - d; |
---|
1623 | |
---|
1624 | if (bs == k-1) { |
---|
1625 | size++; |
---|
1626 | if (size == 1) first_gs = gs; |
---|
1627 | FetchGiantStep(g, gs, F); |
---|
1628 | sub(buf[size-1], g, BabyStep[bs]); |
---|
1629 | } |
---|
1630 | else { |
---|
1631 | sub(h, g, BabyStep[bs]); |
---|
1632 | MulMod(buf[size-1], buf[size-1], h, F); |
---|
1633 | } |
---|
1634 | |
---|
1635 | if (size == ZZ_pX_GCDTableSize && bs == 0) { |
---|
1636 | NewProcessTable(u, f, F, buf, size, first_gs, k, verbose); |
---|
1637 | size = 0; |
---|
1638 | } |
---|
1639 | |
---|
1640 | d++; |
---|
1641 | |
---|
1642 | if (2*d <= deg(f) && deg(f) < old_n) { |
---|
1643 | build(F, f); |
---|
1644 | |
---|
1645 | long i; |
---|
1646 | for (i = 1; i <= k-1; i++) |
---|
1647 | rem(BabyStep[i], BabyStep[i], F); |
---|
1648 | } |
---|
1649 | } |
---|
1650 | |
---|
1651 | if (size > 0) { |
---|
1652 | NewProcessTable(u, f, F, buf, size, first_gs, k, verbose); |
---|
1653 | } |
---|
1654 | |
---|
1655 | if (deg(f) > 0) |
---|
1656 | NewAddFactor(u, f, 0, verbose); |
---|
1657 | |
---|
1658 | } |
---|
1659 | |
---|
1660 | |
---|
1661 | static |
---|
1662 | void IntervalRefine(vec_pair_ZZ_pX_long& factors, const ZZ_pX& ff, |
---|
1663 | long k, long gs, const vec_ZZ_pX& BabyStep, long verbose) |
---|
1664 | |
---|
1665 | { |
---|
1666 | vec_ZZ_pX buf(INIT_SIZE, ZZ_pX_GCDTableSize); |
---|
1667 | |
---|
1668 | ZZ_pX f; |
---|
1669 | f = ff; |
---|
1670 | |
---|
1671 | ZZ_pXModulus F; |
---|
1672 | build(F, f); |
---|
1673 | |
---|
1674 | ZZ_pX g; |
---|
1675 | |
---|
1676 | FetchGiantStep(g, gs, F); |
---|
1677 | |
---|
1678 | long size = 0; |
---|
1679 | |
---|
1680 | long first_d; |
---|
1681 | |
---|
1682 | long d = (gs-1)*k + 1; |
---|
1683 | long bs = k-1; |
---|
1684 | |
---|
1685 | while (bs >= 0 && 2*d <= deg(f)) { |
---|
1686 | |
---|
1687 | long old_n = deg(f); |
---|
1688 | |
---|
1689 | if (size == 0) first_d = d; |
---|
1690 | rem(buf[size], BabyStep[bs], F); |
---|
1691 | sub(buf[size], buf[size], g); |
---|
1692 | size++; |
---|
1693 | |
---|
1694 | if (size == ZZ_pX_GCDTableSize) { |
---|
1695 | NewProcessTable(factors, f, F, buf, size, first_d, 1, verbose); |
---|
1696 | size = 0; |
---|
1697 | } |
---|
1698 | |
---|
1699 | d++; |
---|
1700 | bs--; |
---|
1701 | |
---|
1702 | if (bs >= 0 && 2*d <= deg(f) && deg(f) < old_n) { |
---|
1703 | build(F, f); |
---|
1704 | rem(g, g, F); |
---|
1705 | } |
---|
1706 | } |
---|
1707 | |
---|
1708 | NewProcessTable(factors, f, F, buf, size, first_d, 1, verbose); |
---|
1709 | |
---|
1710 | if (deg(f) > 0) |
---|
1711 | NewAddFactor(factors, f, deg(f), verbose); |
---|
1712 | } |
---|
1713 | |
---|
1714 | |
---|
1715 | |
---|
1716 | |
---|
1717 | static |
---|
1718 | void BabyRefine(vec_pair_ZZ_pX_long& factors, const vec_pair_ZZ_pX_long& u, |
---|
1719 | long k, long l, long verbose) |
---|
1720 | |
---|
1721 | { |
---|
1722 | |
---|
1723 | factors.SetLength(0); |
---|
1724 | |
---|
1725 | vec_ZZ_pX BabyStep; |
---|
1726 | |
---|
1727 | long i; |
---|
1728 | for (i = 0; i < u.length(); i++) { |
---|
1729 | const ZZ_pX& g = u[i].a; |
---|
1730 | long gs = u[i].b; |
---|
1731 | |
---|
1732 | if (gs == 0 || 2*((gs-1)*k+1) > deg(g)) |
---|
1733 | NewAddFactor(factors, g, deg(g), verbose); |
---|
1734 | else { |
---|
1735 | if (BabyStep.length() == 0) |
---|
1736 | FetchBabySteps(BabyStep, k); |
---|
1737 | IntervalRefine(factors, g, k, gs, BabyStep, verbose); |
---|
1738 | } |
---|
1739 | } |
---|
1740 | |
---|
1741 | } |
---|
1742 | |
---|
1743 | |
---|
1744 | |
---|
1745 | |
---|
1746 | |
---|
1747 | |
---|
1748 | void NewDDF(vec_pair_ZZ_pX_long& factors, |
---|
1749 | const ZZ_pX& f, |
---|
1750 | const ZZ_pX& h, |
---|
1751 | long verbose) |
---|
1752 | |
---|
1753 | { |
---|
1754 | if (!IsOne(LeadCoeff(f))) |
---|
1755 | Error("NewDDF: bad args"); |
---|
1756 | |
---|
1757 | if (deg(f) == 0) { |
---|
1758 | factors.SetLength(0); |
---|
1759 | return; |
---|
1760 | } |
---|
1761 | |
---|
1762 | if (deg(f) == 1) { |
---|
1763 | factors.SetLength(0); |
---|
1764 | append(factors, cons(f, 1)); |
---|
1765 | return; |
---|
1766 | } |
---|
1767 | |
---|
1768 | if (!ZZ_pX_stem[0]) |
---|
1769 | sprintf(ZZ_pX_stem, "ddf-%ld", RandomBnd(10000)); |
---|
1770 | |
---|
1771 | long B = deg(f)/2; |
---|
1772 | long k = SqrRoot(B); |
---|
1773 | long l = (B+k-1)/k; |
---|
1774 | |
---|
1775 | ZZ_pX h1; |
---|
1776 | |
---|
1777 | |
---|
1778 | GenerateBabySteps(h1, f, h, k, verbose); |
---|
1779 | |
---|
1780 | GenerateGiantSteps(f, h1, l, verbose); |
---|
1781 | |
---|
1782 | vec_pair_ZZ_pX_long u; |
---|
1783 | GiantRefine(u, f, k, l, verbose); |
---|
1784 | BabyRefine(factors, u, k, l, verbose); |
---|
1785 | |
---|
1786 | FileCleanup(k, l); |
---|
1787 | } |
---|
1788 | |
---|
1789 | NTL_END_IMPL |
---|