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