Line data Source code
1 : /* Copyright (C) 2012 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
8 : ANY WARRANTY WHATSOEVER.
9 :
10 : Check the License for details. You should have received a copy of it, along
11 : with the package; see the file 'COPYING'. If not, write to the Free Software
12 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
13 :
14 : #include "pari.h"
15 : #include "paripriv.h"
16 :
17 : /***********************************************************************/
18 : /** **/
19 : /** Factorisation over finite field **/
20 : /** **/
21 : /***********************************************************************/
22 :
23 : /*******************************************************************/
24 : /* */
25 : /* ROOTS MODULO a prime p (no multiplicities) */
26 : /* */
27 : /*******************************************************************/
28 : /* Replace F by a monic normalized FpX having the same factors;
29 : * assume p prime and *F a ZX */
30 : static int
31 863571 : ZX_factmod_init(GEN *F, GEN p)
32 : {
33 863571 : if (lgefint(p) == 3)
34 : {
35 862187 : ulong pp = p[2];
36 862187 : if (pp == 2) { *F = ZX_to_F2x(*F); return 0; }
37 718367 : *F = ZX_to_Flx(*F, pp);
38 718367 : if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
39 718367 : return 1;
40 : }
41 1384 : *F = FpX_red(*F, p);
42 1384 : if (lg(*F) > 3) *F = FpX_normalize(*F, p);
43 1384 : return 2;
44 : }
45 : static void
46 82678 : ZX_rootmod_init(GEN *F, GEN p)
47 : {
48 82678 : if (lgefint(p) == 3)
49 : {
50 74821 : ulong pp = p[2];
51 74821 : *F = ZX_to_Flx(*F, pp);
52 74821 : if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
53 : }
54 : else
55 : {
56 7857 : *F = FpX_red(*F, p);
57 7857 : if (lg(*F) > 3) *F = FpX_normalize(*F, p);
58 : }
59 82678 : }
60 :
61 : /* return 1,...,p-1 [not_0 = 1] or 0,...,p [not_0 = 0] */
62 : static GEN
63 42 : all_roots_mod_p(ulong p, int not_0)
64 : {
65 : GEN r;
66 : ulong i;
67 42 : if (not_0) {
68 28 : r = cgetg(p, t_VECSMALL);
69 28 : for (i = 1; i < p; i++) r[i] = i;
70 : } else {
71 14 : r = cgetg(p+1, t_VECSMALL);
72 14 : for (i = 0; i < p; i++) r[i+1] = i;
73 : }
74 42 : return r;
75 : }
76 :
77 : /* X^n - 1 */
78 : static GEN
79 126 : Flx_Xnm1(long sv, long n, ulong p)
80 : {
81 126 : GEN t = cgetg(n+3, t_VECSMALL);
82 : long i;
83 126 : t[1] = sv;
84 126 : t[2] = p - 1;
85 126 : for (i = 3; i <= n+1; i++) t[i] = 0;
86 126 : t[i] = 1; return t;
87 : }
88 : /* X^n + 1 */
89 : static GEN
90 140 : Flx_Xn1(long sv, long n, ulong p)
91 : {
92 140 : GEN t = cgetg(n+3, t_VECSMALL);
93 : long i;
94 : (void) p;
95 140 : t[1] = sv;
96 140 : t[2] = 1;
97 140 : for (i = 3; i <= n+1; i++) t[i] = 0;
98 140 : t[i] = 1; return t;
99 : }
100 :
101 : static ulong
102 28 : Fl_nonsquare(ulong p)
103 : {
104 28 : long k = 2;
105 7 : for (;; k++)
106 : {
107 35 : long i = krouu(k, p);
108 35 : if (!i) pari_err_PRIME("Fl_nonsquare",utoipos(p));
109 63 : if (i < 0) return k;
110 7 : }
111 : }
112 :
113 : static GEN
114 8839 : Flx_root_mod_2(GEN f)
115 : {
116 8839 : int z1, z0 = !(f[2] & 1);
117 : long i,n;
118 : GEN y;
119 :
120 8839 : for (i=2, n=1; i < lg(f); i++) n += f[i];
121 8839 : z1 = n & 1;
122 8839 : y = cgetg(z0+z1+1, t_VECSMALL); i = 1;
123 8839 : if (z0) y[i++] = 0;
124 8839 : if (z1) y[i ] = 1;
125 8839 : return y;
126 : }
127 : static ulong
128 56 : Flx_oneroot_mod_2(GEN f)
129 : {
130 : long i,n;
131 56 : if (!(f[2] & 1)) return 0;
132 56 : for (i=2, n=1; i < lg(f); i++) n += f[i];
133 56 : if (n & 1) return 1;
134 28 : return 2;
135 : }
136 :
137 : static GEN FpX_roots_i(GEN f, GEN p);
138 : static GEN Flx_roots_i(GEN f, ulong p);
139 :
140 : static int
141 3185845 : cmpGuGu(GEN a, GEN b) { return (ulong)a < (ulong)b? -1: (a == b? 0: 1); }
142 :
143 : /* Generic driver to computes the roots of f modulo pp, using 'Roots' when
144 : * pp is a small prime.
145 : * if (gpwrap), check types thoroughly and return t_INTMODs, otherwise
146 : * assume that f is an FpX, pp a prime and return t_INTs */
147 : static GEN
148 70526 : rootmod_aux(GEN f, GEN pp)
149 : {
150 : GEN y;
151 70526 : switch(lg(f))
152 : {
153 14 : case 2: pari_err_ROOTS0("rootmod");
154 49 : case 3: return cgetg(1,t_COL);
155 : }
156 70463 : if (typ(f) == t_VECSMALL)
157 : {
158 67468 : ulong p = pp[2];
159 67468 : if (p == 2)
160 8839 : y = Flx_root_mod_2(f);
161 : else
162 : {
163 58629 : if (!odd(p)) pari_err_PRIME("rootmod",utoi(p));
164 58629 : y = Flx_roots_i(f, p);
165 : }
166 67461 : y = Flc_to_ZC(y);
167 : }
168 : else
169 2995 : y = FpX_roots_i(f, pp);
170 70449 : return y;
171 : }
172 : /* assume that f is a ZX and p a prime */
173 : GEN
174 70526 : FpX_roots(GEN f, GEN p)
175 : {
176 70526 : pari_sp av = avma;
177 70526 : GEN y; ZX_rootmod_init(&f, p); y = rootmod_aux(f, p);
178 70498 : return gerepileupto(av, y);
179 : }
180 :
181 : /* assume x reduced mod p > 2, monic. */
182 : static int
183 21 : FpX_quad_factortype(GEN x, GEN p)
184 : {
185 21 : GEN b = gel(x,3), c = gel(x,2);
186 21 : GEN D = subii(sqri(b), shifti(c,2));
187 21 : return kronecker(D,p);
188 : }
189 : /* assume x reduced mod p, monic. Return one root, or NULL if irreducible */
190 : GEN
191 7659 : FpX_quad_root(GEN x, GEN p, int unknown)
192 : {
193 7659 : GEN s, D, b = gel(x,3), c = gel(x,2);
194 :
195 7659 : if (absequaliu(p, 2)) {
196 0 : if (!signe(b)) return c;
197 0 : return signe(c)? NULL: gen_1;
198 : }
199 7659 : D = subii(sqri(b), shifti(c,2));
200 7659 : D = remii(D,p);
201 7659 : if (unknown && kronecker(D,p) == -1) return NULL;
202 :
203 7146 : s = Fp_sqrt(D,p);
204 : /* p is not prime, go on and give e.g. maxord a chance to recover */
205 7146 : if (!s) return NULL;
206 7138 : return Fp_halve(Fp_sub(s,b, p), p);
207 : }
208 : static GEN
209 3140 : FpX_otherroot(GEN x, GEN r, GEN p)
210 3140 : { return Fp_neg(Fp_add(gel(x,3), r, p), p); }
211 :
212 : /* disc(x^2+bx+c) = b^2 - 4c */
213 : static ulong
214 21232860 : Fl_disc_bc(ulong b, ulong c, ulong p)
215 21232860 : { return Fl_sub(Fl_sqr(b,p), Fl_double(Fl_double(c,p),p), p); }
216 : /* p > 2 */
217 : static ulong
218 20508449 : Flx_quad_root(GEN x, ulong p, int unknown)
219 : {
220 20508449 : ulong s, b = x[3], c = x[2];
221 20508449 : ulong D = Fl_disc_bc(b, c, p);
222 20513417 : if (unknown && krouu(D,p) == -1) return p;
223 13708380 : s = Fl_sqrt(D,p);
224 13693563 : if (s==~0UL) return p;
225 13693550 : return Fl_halve(Fl_sub(s,b, p), p);
226 : }
227 : static ulong
228 12273793 : Flx_otherroot(GEN x, ulong r, ulong p)
229 12273793 : { return Fl_neg(Fl_add(x[3], r, p), p); }
230 :
231 :
232 : /* 'todo' contains the list of factors to be split.
233 : * 'done' the list of finished factors, no longer touched */
234 : struct split_t { GEN todo, done; };
235 : static void
236 4917904 : split_init(struct split_t *S, long max)
237 : {
238 4917904 : S->todo = vectrunc_init(max);
239 4917995 : S->done = vectrunc_init(max);
240 4917759 : }
241 : #if 0
242 : /* move todo[i] to done */
243 : static void
244 : split_convert(struct split_t *S, long i)
245 : {
246 : long n = lg(S->todo)-1;
247 : vectrunc_append(S->done, gel(S->todo,i));
248 : if (n) gel(S->todo,i) = gel(S->todo, n);
249 : setlg(S->todo, n);
250 : }
251 : #endif
252 : /* append t to todo */
253 : static void
254 5160812 : split_add(struct split_t *S, GEN t) { vectrunc_append(S->todo, t); }
255 : /* delete todo[i], add t to done */
256 : static void
257 5160872 : split_moveto_done(struct split_t *S, long i, GEN t)
258 : {
259 5160872 : long n = lg(S->todo)-1;
260 5160872 : vectrunc_append(S->done, t);
261 5161409 : if (n) gel(S->todo,i) = gel(S->todo, n);
262 5161409 : setlg(S->todo, n);
263 :
264 5161492 : }
265 : /* append t to done */
266 : static void
267 387967 : split_add_done(struct split_t *S, GEN t)
268 387967 : { vectrunc_append(S->done, t); }
269 : /* split todo[i] into a and b */
270 : static void
271 333341 : split_todo(struct split_t *S, long i, GEN a, GEN b)
272 : {
273 333341 : gel(S->todo, i) = a;
274 333341 : split_add(S, b);
275 333345 : }
276 : /* split todo[i] into a and b, moved to done */
277 : static void
278 368608 : split_done(struct split_t *S, long i, GEN a, GEN b)
279 : {
280 368608 : split_moveto_done(S, i, a);
281 368611 : split_add_done(S, b);
282 368609 : }
283 :
284 : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
285 : static GEN
286 2995 : FpX_roots_i(GEN f, GEN p)
287 : {
288 : GEN pol, pol0, a, q;
289 : struct split_t S;
290 :
291 2995 : split_init(&S, lg(f)-1);
292 2995 : settyp(S.done, t_COL);
293 2995 : if (ZX_valrem(f, &f)) split_add_done(&S, gen_0);
294 2995 : switch(degpol(f))
295 : {
296 7 : case 0: return ZC_copy(S.done);
297 14 : case 1: split_add_done(&S, subii(p, gel(f,2))); return ZC_copy(S.done);
298 : case 2: {
299 1743 : GEN s, r = FpX_quad_root(f, p, 1);
300 1743 : if (r) {
301 1743 : split_add_done(&S, r);
302 1743 : s = FpX_otherroot(f,r, p);
303 : /* f not known to be square free yet */
304 1743 : if (!equalii(r, s)) split_add_done(&S, s);
305 : }
306 1743 : return sort(S.done);
307 : }
308 : }
309 :
310 1231 : a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
311 1231 : if (lg(a) < 3) pari_err_PRIME("rootmod",p);
312 1231 : a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
313 1231 : a = FpX_gcd(f,a, p);
314 1231 : if (!degpol(a)) return ZC_copy(S.done);
315 1231 : split_add(&S, FpX_normalize(a,p));
316 :
317 1231 : q = shifti(p,-1);
318 1231 : pol0 = icopy(gen_1); /* constant term, will vary in place */
319 1231 : pol = deg1pol_shallow(gen_1, pol0, varn(f));
320 2574 : for (pol0[2] = 1;; pol0[2]++)
321 : {
322 2574 : long j, l = lg(S.todo);
323 2574 : if (l == 1) return sort(S.done);
324 1350 : if (pol0[2] == 100 && !BPSW_psp(p)) pari_err_PRIME("polrootsmod",p);
325 2805 : for (j = 1; j < l; j++)
326 : {
327 1462 : GEN b, r, s, c = gel(S.todo,j);
328 1462 : switch(degpol(c))
329 : { /* convert linear and quadratics to roots, try to split the rest */
330 : case 1:
331 77 : split_moveto_done(&S, j, subii(p, gel(c,2)));
332 77 : j--; l--; break;
333 : case 2:
334 1266 : r = FpX_quad_root(c, p, 0);
335 1266 : if (!r) pari_err_PRIME("polrootsmod",p);
336 1259 : s = FpX_otherroot(c,r, p);
337 1259 : split_done(&S, j, r, s);
338 1259 : j--; l--; break;
339 : default:
340 119 : b = FpXQ_pow(pol,q, c,p);
341 119 : if (degpol(b) <= 0) continue;
342 112 : b = FpX_gcd(c,FpX_Fp_sub_shallow(b,gen_1,p), p);
343 112 : if (!degpol(b)) continue;
344 112 : b = FpX_normalize(b, p);
345 112 : c = FpX_div(c,b, p);
346 112 : split_todo(&S, j, b, c);
347 : }
348 : }
349 1343 : }
350 : }
351 :
352 : /* Assume f is normalized */
353 : static ulong
354 173413 : Flx_cubic_root(GEN ff, ulong p)
355 : {
356 173413 : GEN f = Flx_normalize(ff,p);
357 173413 : ulong pi = get_Fl_red(p);
358 173414 : ulong a = f[4], b=f[3], c=f[2], p3 = p%3==1 ? (2*p+1)/3 :(p+1)/3;
359 173414 : ulong t = Fl_mul_pre(a, p3, p, pi), t2 = Fl_sqr_pre(t, p, pi);
360 173415 : ulong A = Fl_sub(b, Fl_triple(t2, p), p);
361 173412 : ulong B = Fl_addmul_pre(c, t, Fl_sub(Fl_double(t2, p), b, p), p, pi);
362 173416 : ulong A3 = Fl_mul_pre(A, p3, p, pi);
363 173416 : ulong A32 = Fl_sqr_pre(A3, p, pi), A33 = Fl_mul_pre(A3, A32, p, pi);
364 173414 : ulong S = Fl_neg(B,p), P = Fl_neg(A3,p);
365 173413 : ulong D = Fl_add(Fl_sqr_pre(S, p, pi), Fl_double(Fl_double(A33, p), p), p);
366 173415 : ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
367 173414 : if (s!=~0UL)
368 : {
369 109293 : ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
370 109293 : if (p%3==2) /* 1 solutions */
371 19240 : vS1 = Fl_powu_pre(S1, (2*p-1)/3, p, pi);
372 : else
373 : {
374 90053 : vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
375 90053 : if (vS1==~0UL) return p; /*0 solutions*/
376 : /*3 solutions*/
377 : }
378 78976 : vS2 = P? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): 0;
379 78975 : return Fl_sub(Fl_add(vS1,vS2, p), t, p);
380 : }
381 : else
382 : {
383 64121 : pari_sp av = avma;
384 64121 : GEN S1 = mkvecsmall2(Fl_halve(S, p), Fl_halve(1UL, p));
385 64123 : GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
386 : ulong Sa;
387 64122 : if (!vS1) return p; /*0 solutions, p%3==2*/
388 64122 : Sa = vS1[1];
389 64122 : if (p%3==1) /*1 solutions*/
390 : {
391 23740 : ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
392 23740 : if (Fa!=P)
393 15785 : Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
394 : }
395 64122 : avma = av;
396 64122 : return Fl_sub(Fl_double(Sa,p),t,p);
397 : }
398 : }
399 :
400 : /* assume p > 2 prime */
401 : static ulong
402 3112832 : Flx_oneroot_i(GEN f, ulong p, long fl)
403 : {
404 : GEN pol, a;
405 : ulong q;
406 : long da;
407 :
408 3112832 : if (Flx_val(f)) return 0;
409 3112129 : switch(degpol(f))
410 : {
411 11754 : case 1: return Fl_neg(f[2], p);
412 2664919 : case 2: return Flx_quad_root(f, p, 1);
413 157651 : case 3: if (p>3) return Flx_cubic_root(f, p); /*FALL THROUGH*/
414 : }
415 :
416 277773 : if (!fl)
417 : {
418 246542 : a = Flxq_powu(polx_Flx(f[1]), p - 1, f,p);
419 246542 : if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
420 246542 : a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
421 246542 : a = Flx_gcd(f,a, p);
422 31231 : } else a = f;
423 277773 : da = degpol(a);
424 277770 : if (!da) return p;
425 198796 : a = Flx_normalize(a,p);
426 :
427 198808 : q = p >> 1;
428 198808 : pol = polx_Flx(f[1]);
429 298646 : for(pol[2] = 1;; pol[2]++)
430 : {
431 298646 : if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
432 298630 : switch(da)
433 : {
434 124102 : case 1: return Fl_neg(a[2], p);
435 59003 : case 2: return Flx_quad_root(a, p, 0);
436 15762 : case 3: if (p>3) return Flx_cubic_root(a, p); /*FALL THROUGH*/
437 : default: {
438 99763 : GEN b = Flxq_powu(pol,q, a,p);
439 : long db;
440 99852 : if (degpol(b) <= 0) continue;
441 94822 : b = Flx_gcd(a,Flx_Fl_add(b,p-1,p), p);
442 94822 : db = degpol(b); if (!db) continue;
443 94821 : b = Flx_normalize(b, p);
444 94832 : if (db <= (da >> 1)) {
445 57708 : a = b;
446 57708 : da = db;
447 : } else {
448 37124 : a = Flx_div(a,b, p);
449 37117 : da -= db;
450 : }
451 : }
452 : }
453 99858 : }
454 : }
455 :
456 : /* assume p > 2 prime */
457 : static GEN
458 4848 : FpX_oneroot_i(GEN f, GEN p)
459 : {
460 : GEN pol, pol0, a, q;
461 : long da;
462 :
463 4848 : if (ZX_val(f)) return gen_0;
464 4582 : switch(degpol(f))
465 : {
466 675 : case 1: return subii(p, gel(f,2));
467 3837 : case 2: return FpX_quad_root(f, p, 1);
468 : }
469 :
470 70 : a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
471 70 : if (lg(a) < 3) pari_err_PRIME("rootmod",p);
472 70 : a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
473 70 : a = FpX_gcd(f,a, p);
474 70 : da = degpol(a);
475 70 : if (!da) return NULL;
476 70 : a = FpX_normalize(a,p);
477 :
478 70 : q = shifti(p,-1);
479 70 : pol0 = icopy(gen_1); /* constant term, will vary in place */
480 70 : pol = deg1pol_shallow(gen_1, pol0, varn(f));
481 224 : for (pol0[2]=1; ; pol0[2]++)
482 : {
483 224 : if (pol0[2] == 1000 && !BPSW_psp(p)) pari_err_PRIME("FpX_oneroot",p);
484 224 : switch(da)
485 : {
486 42 : case 1: return subii(p, gel(a,2));
487 28 : case 2: return FpX_quad_root(a, p, 0);
488 : default: {
489 154 : GEN b = FpXQ_pow(pol,q, a,p);
490 : long db;
491 154 : if (degpol(b) <= 0) continue;
492 147 : b = FpX_gcd(a,FpX_Fp_sub_shallow(b,gen_1,p), p);
493 147 : db = degpol(b); if (!db) continue;
494 147 : b = FpX_normalize(b, p);
495 147 : if (db <= (da >> 1)) {
496 105 : a = b;
497 105 : da = db;
498 : } else {
499 42 : a = FpX_div(a,b, p);
500 42 : da -= db;
501 : }
502 : }
503 : }
504 154 : }
505 : }
506 :
507 : ulong
508 3073491 : Flx_oneroot(GEN f, ulong p)
509 : {
510 3073491 : pari_sp av = avma;
511 : ulong r;
512 3073491 : switch(lg(f))
513 : {
514 0 : case 2: return 0;
515 0 : case 3: avma = av; return p;
516 : }
517 3073491 : if (p == 2) return Flx_oneroot_mod_2(f);
518 3073491 : r = Flx_oneroot_i(Flx_normalize(f, p), p, 0);
519 3073491 : avma = av; return r;
520 : }
521 :
522 : ulong
523 32090 : Flx_oneroot_split(GEN f, ulong p)
524 : {
525 32090 : pari_sp av = avma;
526 : ulong r;
527 32090 : switch(lg(f))
528 : {
529 0 : case 2: return 0;
530 0 : case 3: avma = av; return p;
531 : }
532 32090 : if (p == 2) return Flx_oneroot_mod_2(f);
533 32090 : r = Flx_oneroot_i(Flx_normalize(f, p), p, 1);
534 32162 : avma = av; return r;
535 : }
536 :
537 : /* assume that p is prime */
538 : GEN
539 12152 : FpX_oneroot(GEN f, GEN pp) {
540 12152 : pari_sp av = avma;
541 12152 : ZX_rootmod_init(&f, pp);
542 12152 : switch(lg(f))
543 : {
544 0 : case 2: avma = av; return gen_0;
545 0 : case 3: avma = av; return NULL;
546 : }
547 12152 : if (typ(f) == t_VECSMALL)
548 : {
549 7304 : ulong r, p = pp[2];
550 7304 : if (p == 2)
551 56 : r = Flx_oneroot_mod_2(f);
552 : else
553 7248 : r = Flx_oneroot_i(f, p, 0);
554 7304 : avma = av;
555 7304 : return (r == p)? NULL: utoi(r);
556 : }
557 4848 : f = FpX_oneroot_i(f, pp);
558 4848 : if (!f) { avma = av; return NULL; }
559 4848 : return gerepileuptoint(av, f);
560 : }
561 :
562 : /* returns a root of unity in F_p that is suitable for finding a factor */
563 : /* of degree deg_factor of a polynomial of degree deg; the order is */
564 : /* returned in n */
565 : /* A good choice seems to be n close to deg/deg_factor; we choose n */
566 : /* twice as big and decrement until it divides p-1. */
567 : static GEN
568 217 : good_root_of_unity(GEN p, long deg, long deg_factor, long *pt_n)
569 : {
570 217 : pari_sp ltop = avma;
571 : GEN pm, factn, power, base, zeta;
572 : long n;
573 :
574 217 : pm = subis (p, 1ul);
575 217 : for (n = deg / 2 / deg_factor + 1; !dvdiu (pm, n); n--);
576 217 : factn = Z_factor(stoi(n));
577 217 : power = diviuexact (pm, n);
578 217 : base = gen_1;
579 : do {
580 378 : base = addis (base, 1l);
581 378 : zeta = Fp_pow (base, power, p);
582 : }
583 378 : while (!equaliu (Fp_order (zeta, factn, p), n));
584 217 : *pt_n = n;
585 217 : return gerepileuptoint (ltop, zeta);
586 : }
587 :
588 : GEN
589 924 : FpX_oneroot_split(GEN fact, GEN p)
590 : {
591 924 : pari_sp av = avma;
592 : long n, deg_f, i, dmin;
593 : GEN prim, expo, minfactor, xplusa, zeta, xpow;
594 924 : fact = FpX_normalize(fact, p);
595 924 : deg_f = degpol(fact);
596 924 : if (deg_f<=2) return FpX_oneroot(fact, p);
597 217 : minfactor = fact; /* factor of minimal degree found so far */
598 217 : dmin = degpol(minfactor);
599 217 : prim = good_root_of_unity(p, deg_f, 1, &n);
600 217 : expo = diviuexact(subiu(p, 1), n);
601 217 : xplusa = pol_x(varn(fact));
602 217 : zeta = gen_1;
603 1071 : while (dmin != 1)
604 : {
605 : /* split minfactor by computing its gcd with (X+a)^exp-zeta, where */
606 : /* zeta varies over the roots of unity in F_p */
607 637 : fact = minfactor; deg_f = dmin;
608 : /* update X+a, avoid a=0 */
609 637 : gel (xplusa, 2) = addis (gel (xplusa, 2), 1);
610 637 : xpow = FpXQ_pow (xplusa, expo, fact, p);
611 1260 : for (i = 0; i < n; i++)
612 : {
613 980 : GEN tmp = FpX_gcd(FpX_Fp_sub(xpow, zeta, p), fact, p);
614 980 : long dtmp = degpol(tmp);
615 980 : if (dtmp > 0 && dtmp < deg_f)
616 : {
617 427 : fact = FpX_div(fact, tmp, p); deg_f = degpol(fact);
618 427 : if (dtmp < dmin)
619 : {
620 427 : minfactor = FpX_normalize (tmp, p);
621 427 : dmin = dtmp;
622 427 : if (dmin == 1 || dmin <= deg_f / (n / 2) + 1)
623 : /* stop early to avoid too many gcds */
624 : break;
625 : }
626 : }
627 623 : zeta = Fp_mul (zeta, prim, p);
628 : }
629 : }
630 217 : return gerepileuptoint(av, Fp_neg(gel(minfactor,2), p));
631 : }
632 :
633 : /*******************************************************************/
634 : /* */
635 : /* FACTORISATION MODULO p */
636 : /* */
637 : /*******************************************************************/
638 :
639 : /* F / E a vector of vectors of factors / exponents of virtual length l
640 : * (their real lg may be larger). Set their lg to j, concat and return [F,E] */
641 : static GEN
642 476805 : FE_concat(GEN F, GEN E, long l)
643 : {
644 476805 : setlg(E,l); E = shallowconcat1(E);
645 476805 : setlg(F,l); F = shallowconcat1(F); return mkvec2(F,E);
646 : }
647 :
648 : static GEN
649 14 : ddf_to_ddf2_i(GEN V, long fl)
650 : {
651 : GEN F, D;
652 14 : long i, j, l = lg(V);
653 14 : F = cgetg(l, t_VEC);
654 14 : D = cgetg(l, t_VECSMALL);
655 112 : for (i = j = 1; i < l; i++)
656 : {
657 98 : GEN Vi = gel(V,i);
658 98 : if ((fl==2 && F2x_degree(Vi) == 0)
659 84 : ||(fl==0 && degpol(Vi) == 0)) continue;
660 35 : gel(F,j) = Vi;
661 35 : uel(D,j) = i; j++;
662 : }
663 14 : setlg(F,j);
664 14 : setlg(D,j); return mkvec2(F,D);
665 : }
666 :
667 : GEN
668 7 : ddf_to_ddf2(GEN V)
669 7 : { return ddf_to_ddf2_i(V, 0); }
670 :
671 : static GEN
672 7 : F2x_ddf_to_ddf2(GEN V)
673 7 : { return ddf_to_ddf2_i(V, 2); }
674 :
675 : GEN
676 354779 : vddf_to_simplefact(GEN V, long d)
677 : {
678 : GEN E, F;
679 354779 : long i, j, c, l = lg(V);
680 354779 : F = cgetg(d+1, t_VECSMALL);
681 354779 : E = cgetg(d+1, t_VECSMALL);
682 712611 : for (i = c = 1; i < l; i++)
683 : {
684 357832 : GEN Vi = gel(V,i);
685 357832 : long l = lg(Vi);
686 2463663 : for (j = 1; j < l; j++)
687 : {
688 2105831 : long k, n = degpol(gel(Vi,j)) / j;
689 2105831 : for (k = 1; k <= n; k++) { uel(F,c) = j; uel(E,c) = i; c++; }
690 : }
691 : }
692 354779 : setlg(F,c);
693 354779 : setlg(E,c);
694 354779 : return sort_factor(mkvec2(F,E), (void*)&cmpGuGu, cmp_nodata);
695 : }
696 :
697 : /* product of terms of degree 1 in factorization of f */
698 : GEN
699 142606 : FpX_split_part(GEN f, GEN p)
700 : {
701 142606 : long n = degpol(f);
702 142606 : GEN z, X = pol_x(varn(f));
703 142606 : if (n <= 1) return f;
704 142291 : f = FpX_red(f, p);
705 142291 : z = FpX_sub(FpX_Frobenius(f, p), X, p);
706 142291 : return FpX_gcd(z,f,p);
707 : }
708 :
709 : /* Compute the number of roots in Fp without counting multiplicity
710 : * return -1 for 0 polynomial. lc(f) must be prime to p. */
711 : long
712 99148 : FpX_nbroots(GEN f, GEN p)
713 : {
714 99148 : pari_sp av = avma;
715 99148 : GEN z = FpX_split_part(f, p);
716 99148 : avma = av; return degpol(z);
717 : }
718 :
719 : int
720 0 : FpX_is_totally_split(GEN f, GEN p)
721 : {
722 0 : long n=degpol(f);
723 0 : pari_sp av = avma;
724 0 : if (n <= 1) return 1;
725 0 : if (abscmpui(n, p) > 0) return 0;
726 0 : f = FpX_red(f, p);
727 0 : avma = av; return gequalX(FpX_Frobenius(f, p));
728 : }
729 :
730 : long
731 2403613 : Flx_nbroots(GEN f, ulong p)
732 : {
733 2403613 : long n = degpol(f);
734 2403613 : pari_sp av = avma;
735 : GEN z;
736 2403613 : if (n <= 1) return n;
737 2401408 : if (n == 2)
738 : {
739 : ulong D;
740 11319 : if (p==2) return (f[2]==0) + (f[2]!=f[3]);
741 10437 : D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
742 10437 : return 1 + krouu(D,p);
743 : }
744 2390089 : z = Flx_sub(Flx_Frobenius(f, p), polx_Flx(f[1]), p);
745 2390088 : z = Flx_gcd(z, f, p);
746 2390089 : avma = av; return degpol(z);
747 : }
748 :
749 : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
750 : static GEN
751 5078 : FpX_ddf_Shoup(GEN T, GEN XP, GEN p)
752 : {
753 : GEN b, g, h, F, f, Tr, xq;
754 : long i, j, n, v, B, l, m;
755 : pari_timer ti;
756 :
757 5078 : n = get_FpX_degree(T); v = get_FpX_var(T);
758 5078 : if (n == 0) return cgetg(1, t_VEC);
759 5078 : if (n == 1) return mkvec(get_FpX_mod(T));
760 5032 : B = n/2;
761 5032 : l = usqrt(B);
762 5032 : m = (B+l-1)/l;
763 5032 : T = FpX_get_red(T, p);
764 5032 : b = cgetg(l+2, t_VEC);
765 5032 : gel(b, 1) = pol_x(v);
766 5032 : gel(b, 2) = XP;
767 5032 : if (DEBUGLEVEL>=7) timer_start(&ti);
768 5032 : xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1), T, p);
769 5032 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq baby");
770 12267 : for (i = 3; i <= l+1; i++)
771 7235 : gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
772 5032 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: baby");
773 5032 : xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1), T, p);
774 5032 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq giant");
775 5032 : g = cgetg(m+1, t_VEC);
776 5032 : gel(g, 1) = gel(xq, 2);
777 5032 : for(i = 2; i <= m; i++) gel(g, i) = FpX_FpXQV_eval(gel(g, i-1), xq, T, p);
778 5032 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: giant");
779 5032 : h = cgetg(m+1, t_VEC);
780 22411 : for (j = 1; j <= m; j++)
781 : {
782 17379 : pari_sp av = avma;
783 17379 : GEN gj = gel(g,j), e = FpX_sub(gj, gel(b,1), p);
784 17379 : for (i = 2; i <= l; i++) e = FpXQ_mul(e, FpX_sub(gj, gel(b,i), p), T, p);
785 17379 : gel(h,j) = gerepileupto(av, e);
786 : }
787 5032 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: diff");
788 5032 : Tr = get_FpX_mod(T);
789 5032 : F = cgetg(m+1, t_VEC);
790 22411 : for (j = 1; j <= m; j++)
791 : {
792 17379 : GEN u = FpX_gcd(Tr, gel(h,j), p);
793 17379 : if (degpol(u))
794 : {
795 3031 : u = FpX_normalize(u, p);
796 3031 : Tr = FpX_div(Tr, u, p);
797 : }
798 17379 : gel(F,j) = u;
799 : }
800 5032 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: F");
801 5032 : f = const_vec(n, pol_1(v));
802 22411 : for (j = 1; j <= m; j++)
803 : {
804 17379 : GEN e = gel(F, j);
805 19507 : for (i=l-1; i >= 0; i--)
806 : {
807 19507 : GEN u = FpX_gcd(e, FpX_sub(gel(g, j), gel(b, i+1), p), p);
808 19507 : if (degpol(u))
809 : {
810 3038 : u = FpX_normalize(u, p);
811 3038 : gel(f, l*j-i) = u;
812 3038 : e = FpX_div(e, u, p);
813 : }
814 19507 : if (!degpol(e)) break;
815 : }
816 : }
817 5032 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: f");
818 5032 : if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
819 5032 : return f;
820 : }
821 :
822 : static void
823 0 : FpX_edf_simple(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
824 : {
825 0 : long n = degpol(Tp), r = n/d, ct = 0;
826 : GEN T, f, ff, p2;
827 0 : if (r==1) { gel(V, idx) = Tp; return; }
828 0 : p2 = shifti(p,-1);
829 0 : T = FpX_get_red(Tp, p);
830 0 : XP = FpX_rem(XP, T, p);
831 : while (1)
832 : {
833 0 : pari_sp btop = avma;
834 : long i;
835 0 : GEN g = random_FpX(n, varn(Tp), p);
836 0 : GEN t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
837 0 : if (signe(t) == 0) continue;
838 0 : for(i=1; i<=10; i++)
839 : {
840 0 : pari_sp btop2 = avma;
841 0 : GEN R = FpXQ_pow(FpX_Fp_add(t, randomi(p), p), p2, T, p);
842 0 : f = FpX_gcd(FpX_Fp_sub(R, gen_1, p), Tp, p);
843 0 : if (degpol(f) > 0 && degpol(f) < n) break;
844 0 : avma = btop2;
845 : }
846 0 : if (degpol(f) > 0 && degpol(f) < n) break;
847 0 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_simple",p);
848 0 : avma = btop;
849 0 : }
850 0 : f = FpX_normalize(f, p);
851 0 : ff = FpX_div(Tp, f ,p);
852 0 : FpX_edf_simple(f, XP, d, p, V, idx);
853 0 : FpX_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
854 : }
855 :
856 : static void
857 448 : FpX_edf_rec(GEN T, GEN hp, GEN t, long d, GEN p2, GEN p, GEN V, long idx)
858 : {
859 : pari_sp av;
860 448 : GEN Tp = get_FpX_mod(T);
861 448 : long n = degpol(hp), vT = varn(Tp), ct = 0;
862 : GEN u1, u2, f1, f2, R, h;
863 448 : h = FpX_get_red(hp, p);
864 448 : t = FpX_rem(t, T, p);
865 448 : av = avma;
866 : do
867 : {
868 736 : avma = av;
869 736 : R = FpXQ_pow(deg1pol(gen_1, randomi(p), vT), p2, h, p);
870 736 : u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
871 736 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_rec",p);
872 736 : } while (degpol(u1)==0 || degpol(u1)==n);
873 448 : f1 = FpX_gcd(FpX_FpXQ_eval(u1, t, T, p), Tp, p);
874 448 : f1 = FpX_normalize(f1, p);
875 448 : u2 = FpX_div(hp, u1, p);
876 448 : f2 = FpX_div(Tp, f1, p);
877 448 : if (degpol(u1)==1)
878 329 : gel(V, idx) = f1;
879 : else
880 119 : FpX_edf_rec(FpX_get_red(f1, p), u1, t, d, p2, p, V, idx);
881 448 : idx += degpol(f1)/d;
882 448 : if (degpol(u2)==1)
883 341 : gel(V, idx) = f2;
884 : else
885 107 : FpX_edf_rec(FpX_get_red(f2, p), u2, t, d, p2, p, V, idx);
886 448 : }
887 :
888 : /* assume Tp a squarefree product of r > 1 irred. factors of degree d */
889 : static void
890 222 : FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
891 : {
892 222 : long n = degpol(Tp), r = n/d, vT = varn(Tp), ct = 0;
893 : GEN T, h, t;
894 : pari_timer ti;
895 :
896 222 : T = FpX_get_red(Tp, p);
897 222 : XP = FpX_rem(XP, T, p);
898 222 : if (DEBUGLEVEL>=7) timer_start(&ti);
899 : do
900 : {
901 222 : GEN g = random_FpX(n, vT, p);
902 222 : t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
903 222 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
904 222 : h = FpXQ_minpoly(t, T, p);
905 222 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
906 222 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf",p);
907 222 : } while (degpol(h) != r);
908 222 : FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
909 222 : }
910 :
911 : static GEN
912 654 : FpX_factor_Shoup(GEN T, GEN p)
913 : {
914 654 : long i, n, s = 0;
915 : GEN XP, D, V;
916 654 : long e = expi(p);
917 : pari_timer ti;
918 654 : n = get_FpX_degree(T);
919 654 : T = FpX_get_red(T, p);
920 654 : if (DEBUGLEVEL>=6) timer_start(&ti);
921 654 : XP = FpX_Frobenius(T, p);
922 654 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
923 654 : D = FpX_ddf_Shoup(T, XP, p);
924 654 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
925 654 : s = ddf_to_nbfact(D);
926 654 : V = cgetg(s+1, t_COL);
927 5496 : for (i = 1, s = 1; i <= n; i++)
928 : {
929 4842 : GEN Di = gel(D,i);
930 4842 : long ni = degpol(Di), ri = ni/i;
931 4842 : if (ni == 0) continue;
932 663 : Di = FpX_normalize(Di, p);
933 663 : if (ni == i) { gel(V, s++) = Di; continue; }
934 222 : if (ri <= e*expu(e))
935 222 : FpX_edf(Di, XP, i, p, V, s);
936 : else
937 0 : FpX_edf_simple(Di, XP, i, p, V, s);
938 222 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
939 222 : s += ri;
940 : }
941 654 : return V;
942 : }
943 :
944 : long
945 529985 : ddf_to_nbfact(GEN D)
946 : {
947 529985 : long l = lg(D), i, s = 0;
948 529985 : for(i = 1; i < l; i++) s += degpol(gel(D,i))/i;
949 529985 : return s;
950 : }
951 :
952 : /* Yun algorithm: Assume p > degpol(T) */
953 : static GEN
954 671 : FpX_factor_Yun(GEN T, GEN p)
955 : {
956 671 : long n = degpol(T), i = 1;
957 671 : GEN a, b, c, d = FpX_deriv(T, p);
958 671 : GEN V = cgetg(n+1,t_VEC);
959 671 : a = FpX_gcd(T, d, p);
960 671 : if (degpol(a) == 0) return mkvec(T);
961 488 : b = FpX_div(T, a, p);
962 : do
963 : {
964 2432 : c = FpX_div(d, a, p);
965 2432 : d = FpX_sub(c, FpX_deriv(b, p), p);
966 2432 : a = FpX_normalize(FpX_gcd(b, d, p), p);
967 2432 : gel(V, i++) = a;
968 2432 : b = FpX_div(b, a, p);
969 2432 : } while (degpol(b));
970 488 : setlg(V, i); return V;
971 : }
972 : GEN
973 749 : FpX_factor_squarefree(GEN T, GEN p)
974 : {
975 749 : if (lgefint(p)==3)
976 : {
977 749 : ulong pp = (ulong)p[2];
978 749 : GEN u = Flx_factor_squarefree(ZX_to_Flx(T,pp), pp);
979 749 : return FlxV_to_ZXV(u);
980 : }
981 0 : return FpX_factor_Yun(T, p);
982 : }
983 :
984 : long
985 168 : FpX_ispower(GEN f, ulong k, GEN p, GEN *pt_r)
986 : {
987 168 : pari_sp av = avma;
988 : GEN lc, F;
989 168 : long i, l, n = degpol(f), v = varn(f);
990 168 : if (n % k) return 0;
991 168 : if (lgefint(p)==3)
992 : {
993 126 : ulong pp = p[2];
994 126 : GEN fp = ZX_to_Flx(f, pp);
995 126 : if (!Flx_ispower(fp, k, pp, pt_r)) { avma = av; return 0; }
996 105 : if (pt_r) *pt_r = gerepileupto(av, Flx_to_ZX(*pt_r)); else avma = av;
997 105 : return 1;
998 : }
999 42 : lc = Fp_sqrtn(leading_coeff(f), stoi(k), p, NULL);
1000 42 : if (!lc) { av = avma; return 0; }
1001 42 : F = FpX_factor_Yun(f, p); l = lg(F)-1;
1002 1491 : for(i=1; i <= l; i++)
1003 1456 : if (i%k && degpol(gel(F,i))) { avma = av; return 0; }
1004 35 : if (pt_r)
1005 : {
1006 35 : GEN r = scalarpol(lc, v), s = pol_1(v);
1007 1484 : for (i=l; i>=1; i--)
1008 : {
1009 1449 : if (i%k) continue;
1010 294 : s = FpX_mul(s, gel(F,i), p);
1011 294 : r = FpX_mul(r, s, p);
1012 : }
1013 35 : *pt_r = gerepileupto(av, r);
1014 0 : } else av = avma;
1015 35 : return 1;
1016 : }
1017 :
1018 : static GEN
1019 615 : FpX_factor_Cantor(GEN T, GEN p)
1020 : {
1021 615 : GEN E, F, V = FpX_factor_Yun(T, p);
1022 615 : long i, j, l = lg(V);
1023 615 : F = cgetg(l, t_VEC);
1024 615 : E = cgetg(l, t_VEC);
1025 1746 : for (i=1, j=1; i < l; i++)
1026 1131 : if (degpol(gel(V,i)))
1027 : {
1028 654 : GEN Fj = FpX_factor_Shoup(gel(V,i), p);
1029 654 : gel(F, j) = Fj;
1030 654 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1031 654 : j++;
1032 : }
1033 615 : return sort_factor_pol(FE_concat(F,E,j), cmpii);
1034 : }
1035 :
1036 : static GEN
1037 0 : FpX_ddf_i(GEN T, GEN p)
1038 : {
1039 : GEN XP;
1040 0 : T = FpX_get_red(T, p);
1041 0 : XP = FpX_Frobenius(T, p);
1042 0 : return ddf_to_ddf2(FpX_ddf_Shoup(T, XP, p));
1043 : }
1044 :
1045 : GEN
1046 7 : FpX_ddf(GEN f, GEN p)
1047 : {
1048 7 : pari_sp av = avma;
1049 : GEN F;
1050 7 : switch(ZX_factmod_init(&f, p))
1051 : {
1052 7 : case 0: F = F2x_ddf(f);
1053 7 : F2xV_to_ZXV_inplace(gel(F,1)); break;
1054 0 : case 1: F = Flx_ddf(f,p[2]);
1055 0 : FlxV_to_ZXV_inplace(gel(F,1)); break;
1056 0 : default: F = FpX_ddf_i(f,p); break;
1057 : }
1058 7 : return gerepilecopy(av, F);
1059 : }
1060 :
1061 : static GEN Flx_simplefact_Cantor(GEN T, ulong p);
1062 : static GEN
1063 14 : FpX_simplefact_Cantor(GEN T, GEN p)
1064 : {
1065 : GEN V, XP;
1066 : long i, l;
1067 14 : if (lgefint(p) == 3)
1068 : {
1069 0 : ulong pp = p[2];
1070 0 : return Flx_simplefact_Cantor(ZX_to_Flx(T,pp), pp);
1071 : }
1072 14 : T = FpX_get_red(T, p);
1073 14 : XP = FpX_Frobenius(T, p);
1074 14 : V = FpX_factor_Yun(get_FpX_mod(T), p); l = lg(V);
1075 14 : for (i=1; i < l; i++) gel(V,i) = FpX_ddf_Shoup(gel(V,i), XP, p);
1076 14 : return vddf_to_simplefact(V, get_FpX_degree(T));
1077 : }
1078 :
1079 : static int
1080 0 : FpX_isirred_Cantor(GEN Tp, GEN p)
1081 : {
1082 0 : pari_sp av = avma;
1083 : pari_timer ti;
1084 : long n, d;
1085 0 : GEN T = get_FpX_mod(Tp);
1086 0 : GEN dT = FpX_deriv(T, p);
1087 : GEN XP, D;
1088 0 : if (degpol(FpX_gcd(T, dT, p)) != 0) { avma = av; return 0; }
1089 0 : n = get_FpX_degree(T);
1090 0 : T = FpX_get_red(Tp, p);
1091 0 : if (DEBUGLEVEL>=6) timer_start(&ti);
1092 0 : XP = FpX_Frobenius(T, p);
1093 0 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
1094 0 : D = FpX_ddf_Shoup(T, XP, p);
1095 0 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
1096 0 : d = degpol(gel(D, n));
1097 0 : avma = av; return d==n;
1098 : }
1099 :
1100 : static GEN FpX_factor_deg2(GEN f, GEN p, long d, long flag);
1101 :
1102 : /*Assume that p is large and odd*/
1103 : static GEN
1104 1384 : FpX_factor_i(GEN f, GEN pp, long flag)
1105 : {
1106 1384 : long d = degpol(f);
1107 1384 : if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
1108 629 : switch(flag)
1109 : {
1110 615 : default: return FpX_factor_Cantor(f, pp);
1111 14 : case 1: return FpX_simplefact_Cantor(f, pp);
1112 0 : case 2: return FpX_isirred_Cantor(f, pp)? gen_1: NULL;
1113 : }
1114 : }
1115 :
1116 : long
1117 4410 : FpX_nbfact_Frobenius(GEN T, GEN XP, GEN p)
1118 : {
1119 4410 : pari_sp av = avma;
1120 4410 : long s = ddf_to_nbfact(FpX_ddf_Shoup(T, XP, p));
1121 4410 : avma = av; return s;
1122 : }
1123 :
1124 : long
1125 0 : FpX_nbfact(GEN T, GEN p)
1126 : {
1127 0 : pari_sp av = avma;
1128 0 : GEN XP = FpX_Frobenius(T, p);
1129 0 : long n = FpX_nbfact_Frobenius(T, XP, p);
1130 0 : avma = av; return n;
1131 : }
1132 :
1133 : /* p > 2 */
1134 : static GEN
1135 7 : FpX_is_irred_2(GEN f, GEN p, long d)
1136 : {
1137 7 : switch(d)
1138 : {
1139 : case -1:
1140 0 : case 0: return NULL;
1141 0 : case 1: return gen_1;
1142 : }
1143 7 : return FpX_quad_factortype(f, p) == -1? gen_1: NULL;
1144 : }
1145 : /* p > 2 */
1146 : static GEN
1147 14 : FpX_degfact_2(GEN f, GEN p, long d)
1148 : {
1149 14 : switch(d)
1150 : {
1151 0 : case -1:retmkvec2(mkvecsmall(-1),mkvecsmall(1));
1152 0 : case 0: return trivial_fact();
1153 0 : case 1: retmkvec2(mkvecsmall(1), mkvecsmall(1));
1154 : }
1155 14 : switch(FpX_quad_factortype(f, p)) {
1156 7 : case 1: retmkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1157 7 : case -1: retmkvec2(mkvecsmall(2), mkvecsmall(1));
1158 0 : default: retmkvec2(mkvecsmall(1), mkvecsmall(2));
1159 : }
1160 : }
1161 :
1162 : GEN
1163 63 : prime_fact(GEN x) { retmkmat2(mkcolcopy(x), mkcol(gen_1)); }
1164 : GEN
1165 1029 : trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
1166 :
1167 : /* not gerepile safe */
1168 : static GEN
1169 734 : FpX_factor_2(GEN f, GEN p, long d)
1170 : {
1171 : GEN r, s, R, S;
1172 : long v;
1173 : int sgn;
1174 734 : switch(d)
1175 : {
1176 7 : case -1: retmkvec2(mkcol(pol_0(varn(f))), mkvecsmall(1));
1177 30 : case 0: retmkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1178 45 : case 1: retmkvec2(mkcol(f), mkvecsmall(1));
1179 : }
1180 652 : r = FpX_quad_root(f, p, 1);
1181 652 : if (!r) return mkvec2(mkcol(f), mkvecsmall(1));
1182 138 : v = varn(f);
1183 138 : s = FpX_otherroot(f, r, p);
1184 138 : if (signe(r)) r = subii(p, r);
1185 138 : if (signe(s)) s = subii(p, s);
1186 138 : sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
1187 138 : R = deg1pol_shallow(gen_1, r, v);
1188 138 : if (!sgn) return mkvec2(mkcol(R), mkvecsmall(2));
1189 47 : S = deg1pol_shallow(gen_1, s, v);
1190 47 : return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
1191 : }
1192 : static GEN
1193 755 : FpX_factor_deg2(GEN f, GEN p, long d, long flag)
1194 : {
1195 755 : switch(flag) {
1196 7 : case 2: return FpX_is_irred_2(f, p, d);
1197 14 : case 1: return FpX_degfact_2(f, p, d);
1198 734 : default: return FpX_factor_2(f, p, d);
1199 : }
1200 : }
1201 :
1202 : static int
1203 57555 : F2x_quad_factortype(GEN x)
1204 57555 : { return x[2] == 7 ? -1: x[2] == 6 ? 1 :0; }
1205 :
1206 : static GEN
1207 7 : F2x_is_irred_2(GEN f, long d)
1208 7 : { return d == 1 || (d==2 && F2x_quad_factortype(f) == -1)? gen_1: NULL; }
1209 :
1210 : static GEN
1211 7105 : F2x_degfact_2(GEN f, long d)
1212 : {
1213 7105 : if (!d) return trivial_fact();
1214 7105 : if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1215 6902 : switch(F2x_quad_factortype(f)) {
1216 2261 : case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1217 2219 : case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1218 2422 : default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1219 : }
1220 : }
1221 :
1222 : static GEN
1223 104454 : F2x_factor_2(GEN f, long d)
1224 : {
1225 104454 : long v = f[1];
1226 104454 : if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1227 96607 : if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1228 49336 : switch(F2x_quad_factortype(f))
1229 : {
1230 12467 : case -1: return mkvec2(mkcol(f), mkvecsmall(1));
1231 22498 : case 0: return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
1232 14371 : default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
1233 : }
1234 : }
1235 : static GEN
1236 111566 : F2x_factor_deg2(GEN f, long d, long flag)
1237 : {
1238 111566 : switch(flag) {
1239 7 : case 2: return F2x_is_irred_2(f, d);
1240 7105 : case 1: return F2x_degfact_2(f, d);
1241 104454 : default: return F2x_factor_2(f, d);
1242 : }
1243 : }
1244 :
1245 : /* xt = NULL or x^(p-1)/2 mod g */
1246 : static void
1247 4816 : split_squares(struct split_t *S, GEN g, ulong p, GEN xt)
1248 : {
1249 4816 : ulong q = p >> 1;
1250 4816 : GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
1251 4816 : long d = degpol(a);
1252 4816 : if (d < 0)
1253 : {
1254 : ulong i;
1255 42 : split_add_done(S, (GEN)1);
1256 42 : for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr(i,p));
1257 : } else {
1258 4774 : if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
1259 4774 : if (d)
1260 : {
1261 4774 : if (xt) xt = Flx_Fl_add(xt, p-1, p); else xt = Flx_Xnm1(g[1], q, p);
1262 4774 : a = Flx_gcd(a, xt, p);
1263 4774 : if (degpol(a)) split_add(S, Flx_normalize(a, p));
1264 : }
1265 : }
1266 4816 : }
1267 : static void
1268 4816 : split_nonsquares(struct split_t *S, GEN g, ulong p, GEN xt)
1269 : {
1270 4816 : ulong q = p >> 1;
1271 4816 : GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
1272 4816 : long d = degpol(a);
1273 4816 : if (d < 0)
1274 : {
1275 28 : ulong i, z = Fl_nonsquare(p);
1276 28 : split_add_done(S, (GEN)z);
1277 28 : for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_mul(z, Fl_sqr(i,p), p));
1278 : } else {
1279 4788 : if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
1280 4788 : if (d)
1281 : {
1282 4788 : if (xt) xt = Flx_Fl_add(xt, 1, p); else xt = Flx_Xn1(g[1], q, p);
1283 4788 : a = Flx_gcd(a, xt, p);
1284 4788 : if (degpol(a)) split_add(S, Flx_normalize(a, p));
1285 : }
1286 : }
1287 4816 : }
1288 : /* p > 2. f monic Flx, f(0) != 0. Add to split_t structs coprime factors
1289 : * of g = \prod_{f(a) = 0} (X - a). Return 0 when f(x) = 0 for all x in Fp* */
1290 : static int
1291 4914629 : split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p)
1292 : {
1293 4914629 : GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
1294 4914530 : long d = degpol(g);
1295 4914464 : if (d < 0) return 0;
1296 4914483 : if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
1297 4914656 : if (!d) return 1;
1298 4899067 : if ((p >> 4) <= (ulong)d)
1299 : { /* small p; split directly using x^((p-1)/2) +/- 1 */
1300 14280 : GEN xt = ((ulong)d < (p>>1))? Flx_rem(monomial_Flx(1, p>>1, g[1]), g, p)
1301 9464 : : NULL;
1302 4816 : split_squares(S, g, p, xt);
1303 4816 : split_nonsquares(S, g, p, xt);
1304 : } else { /* large p; use x^(p-1) - 1 directly */
1305 4894251 : a = Flxq_powu(polx_Flx(f[1]), p-1, g,p);
1306 4894153 : if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
1307 4894153 : a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
1308 4893645 : g = Flx_gcd(g,a, p);
1309 4893929 : if (degpol(g)) split_add(S, Flx_normalize(g,p));
1310 : }
1311 4899124 : return 1;
1312 : }
1313 :
1314 : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
1315 : static GEN
1316 22057661 : Flx_roots_i(GEN f, ulong p)
1317 : {
1318 : GEN pol, g;
1319 22057661 : long v = Flx_valrem(f, &g);
1320 : ulong q;
1321 : struct split_t S;
1322 :
1323 : /* optimization: test for small degree first */
1324 22055431 : switch(degpol(g))
1325 : {
1326 : case 1: {
1327 25368 : ulong r = p - g[2];
1328 25368 : return v? mkvecsmall2(0, r): mkvecsmall(r);
1329 : }
1330 : case 2: {
1331 17115915 : ulong r = Flx_quad_root(g, p, 1), s;
1332 17123391 : if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
1333 11706450 : s = Flx_otherroot(g,r, p);
1334 11716617 : if (r < s)
1335 2924256 : return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
1336 8792361 : else if (r > s)
1337 8792200 : return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
1338 : else
1339 161 : return v? mkvecsmall2(0, s): mkvecsmall(s);
1340 : }
1341 : }
1342 4914927 : q = p >> 1;
1343 4914927 : split_init(&S, lg(f)-1);
1344 4914685 : settyp(S.done, t_VECSMALL);
1345 4914685 : if (v) split_add_done(&S, (GEN)0);
1346 4914685 : if (! split_Flx_cut_out_roots(&S, g, p))
1347 42 : return all_roots_mod_p(p, lg(S.done) == 1);
1348 4914482 : pol = polx_Flx(f[1]);
1349 10155351 : for (pol[2]=1; ; pol[2]++)
1350 : {
1351 10155351 : long j, l = lg(S.todo);
1352 10155351 : if (l == 1) { vecsmall_sort(S.done); return S.done; }
1353 5240803 : if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
1354 10832857 : for (j = 1; j < l; j++)
1355 : {
1356 5592506 : GEN b, c = gel(S.todo,j);
1357 : ulong r, s;
1358 5592506 : switch(degpol(c))
1359 : {
1360 : case 1:
1361 4792180 : split_moveto_done(&S, j, (GEN)(p - c[2]));
1362 4792194 : j--; l--; break;
1363 : case 2:
1364 367358 : r = Flx_quad_root(c, p, 0);
1365 367358 : if (r == p) pari_err_PRIME("polrootsmod",utoipos(p));
1366 367351 : s = Flx_otherroot(c,r, p);
1367 367350 : split_done(&S, j, (GEN)r, (GEN)s);
1368 367350 : j--; l--; break;
1369 : default:
1370 432804 : b = Flxq_powu(pol,q, c,p); /* pol^(p-1)/2 */
1371 432797 : if (degpol(b) <= 0) continue;
1372 333364 : b = Flx_gcd(c,Flx_Fl_add(b,p-1,p), p);
1373 333369 : if (!degpol(b)) continue;
1374 333237 : b = Flx_normalize(b, p);
1375 333233 : c = Flx_div(c,b, p);
1376 333231 : split_todo(&S, j, b, c);
1377 : }
1378 : }
1379 5240351 : }
1380 : }
1381 :
1382 : GEN
1383 21994385 : Flx_roots(GEN f, ulong p)
1384 : {
1385 21994385 : pari_sp av = avma;
1386 21994385 : switch(lg(f))
1387 : {
1388 0 : case 2: pari_err_ROOTS0("Flx_roots");
1389 0 : case 3: avma = av; return cgetg(1, t_VECSMALL);
1390 : }
1391 22006386 : if (p == 2) return Flx_root_mod_2(f);
1392 22006386 : return gerepileuptoleaf(av, Flx_roots_i(Flx_normalize(f, p), p));
1393 : }
1394 :
1395 : /* assume x reduced mod p, monic. */
1396 : static int
1397 725592 : Flx_quad_factortype(GEN x, ulong p)
1398 : {
1399 725592 : ulong b = x[3], c = x[2];
1400 725592 : return krouu(Fl_disc_bc(b, c, p), p);
1401 : }
1402 : static GEN
1403 35 : Flx_is_irred_2(GEN f, ulong p, long d)
1404 : {
1405 35 : if (!d) return NULL;
1406 35 : if (d == 1) return gen_1;
1407 35 : return Flx_quad_factortype(f, p) == -1? gen_1: NULL;
1408 : }
1409 : static GEN
1410 742700 : Flx_degfact_2(GEN f, ulong p, long d)
1411 : {
1412 742700 : if (!d) return trivial_fact();
1413 742700 : if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1414 725557 : switch(Flx_quad_factortype(f, p)) {
1415 344687 : case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1416 372428 : case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1417 8442 : default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1418 : }
1419 : }
1420 : /* p > 2 */
1421 : static GEN
1422 422491 : Flx_factor_2(GEN f, ulong p, long d)
1423 : {
1424 : ulong r, s;
1425 : GEN R,S;
1426 422491 : long v = f[1];
1427 422491 : if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1428 408388 : if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1429 309726 : r = Flx_quad_root(f, p, 1);
1430 309726 : if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
1431 200831 : s = Flx_otherroot(f, r, p);
1432 200831 : r = Fl_neg(r, p);
1433 200831 : s = Fl_neg(s, p);
1434 200831 : if (s < r) lswap(s,r);
1435 200831 : R = mkvecsmall3(v,r,1);
1436 200831 : if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
1437 173307 : S = mkvecsmall3(v,s,1);
1438 173307 : return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
1439 : }
1440 : static GEN
1441 1165226 : Flx_factor_deg2(GEN f, ulong p, long d, long flag)
1442 : {
1443 1165226 : switch(flag) {
1444 35 : case 2: return Flx_is_irred_2(f, p, d);
1445 742700 : case 1: return Flx_degfact_2(f, p, d);
1446 422491 : default: return Flx_factor_2(f, p, d);
1447 : }
1448 : }
1449 :
1450 : void
1451 19978 : F2xV_to_FlxV_inplace(GEN v)
1452 : {
1453 : long i;
1454 19978 : for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_Flx(gel(v,i));
1455 19978 : }
1456 : void
1457 720748 : FlxV_to_ZXV_inplace(GEN v)
1458 : {
1459 : long i;
1460 720748 : for(i=1;i<lg(v);i++) gel(v,i)= Flx_to_ZX(gel(v,i));
1461 720748 : }
1462 : void
1463 143806 : F2xV_to_ZXV_inplace(GEN v)
1464 : {
1465 : long i;
1466 143806 : for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_ZX(gel(v,i));
1467 143806 : }
1468 :
1469 : static GEN
1470 10368 : F2x_Berlekamp_ker(GEN u)
1471 : {
1472 10368 : pari_sp ltop=avma;
1473 10368 : long j,N = F2x_degree(u);
1474 : GEN Q;
1475 : pari_timer T;
1476 10368 : timer_start(&T);
1477 10368 : Q = F2x_matFrobenius(u);
1478 243637 : for (j=1; j<=N; j++)
1479 233269 : F2m_flip(Q,j,j);
1480 10368 : if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
1481 10368 : Q = F2m_ker_sp(Q,0);
1482 10368 : if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
1483 10368 : return gerepileupto(ltop,Q);
1484 : }
1485 : #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
1486 : static long
1487 14260 : F2x_split_Berlekamp(GEN *t)
1488 : {
1489 14260 : GEN u = *t, a, b, vker;
1490 14260 : long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
1491 :
1492 14260 : if (du == 1) return 1;
1493 10967 : if (du == 2)
1494 : {
1495 599 : if (F2x_quad_factortype(u) == 1) /* 0 is a root: shouldn't occur */
1496 : {
1497 0 : t[0] = mkvecsmall2(sv, 2);
1498 0 : t[1] = mkvecsmall2(sv, 3);
1499 0 : return 2;
1500 : }
1501 599 : return 1;
1502 : }
1503 :
1504 10368 : vker = F2x_Berlekamp_ker(u);
1505 10368 : lb = lgcols(vker);
1506 10371 : d = lg(vker)-1;
1507 10371 : ir = 0;
1508 : /* t[i] irreducible for i < ir, still to be treated for i < L */
1509 51638 : for (L=1; L<d; )
1510 : {
1511 : GEN pol;
1512 30899 : if (d == 2)
1513 1706 : pol = F2v_to_F2x(gel(vker,2), sv);
1514 : else
1515 : {
1516 29193 : GEN v = zero_zv(lb);
1517 29195 : v[1] = du;
1518 29195 : v[2] = random_Fl(2); /*Assume vker[1]=1*/
1519 108983 : for (i=2; i<=d; i++)
1520 79788 : if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
1521 29195 : pol = F2v_to_F2x(v, sv);
1522 : }
1523 95880 : for (i=ir; i<L && L<d; i++)
1524 : {
1525 64984 : a = t[i]; la = F2x_degree(a);
1526 64976 : if (la == 1) { set_irred(i); }
1527 64805 : else if (la == 2)
1528 : {
1529 711 : if (F2x_quad_factortype(a) == 1) /* 0 is a root: shouldn't occur */
1530 : {
1531 0 : t[i] = mkvecsmall2(sv, 2);
1532 0 : t[L] = mkvecsmall2(sv, 3); L++;
1533 : }
1534 711 : set_irred(i);
1535 : }
1536 : else
1537 : {
1538 64094 : pari_sp av = avma;
1539 : long lb;
1540 64094 : b = F2x_rem(pol, a);
1541 64091 : if (F2x_degree(b) <= 0) { avma=av; continue; }
1542 21149 : b = F2x_gcd(a,b); lb = F2x_degree(b);
1543 21149 : if (lb && lb < la)
1544 : {
1545 21149 : t[L] = F2x_div(a,b);
1546 21151 : t[i]= b; L++;
1547 : }
1548 0 : else avma = av;
1549 : }
1550 : }
1551 : }
1552 10368 : return d;
1553 : }
1554 : /* assume deg f > 2 */
1555 : static GEN
1556 11468 : F2x_Berlekamp_i(GEN f, long flag)
1557 : {
1558 11468 : long lfact, val, d = F2x_degree(f), j, k, lV;
1559 : GEN y, E, t, V;
1560 :
1561 11468 : val = F2x_valrem(f, &f);
1562 11468 : if (flag == 2 && val > 1) return NULL;
1563 11468 : V = F2x_factor_squarefree(f); lV = lg(V);
1564 11468 : if (flag == 2 && lV > 2) return NULL;
1565 :
1566 : /* to hold factors and exponents */
1567 11398 : t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
1568 11398 : E = cgetg(d+1,t_VECSMALL);
1569 11398 : lfact = 1;
1570 11398 : if (val) {
1571 3402 : if (flag == 1) t[1] = 1; else gel(t,1) = polx_F2x(f[1]);
1572 3402 : E[1] = val; lfact++;
1573 : }
1574 :
1575 48456 : for (k=1; k<lV; k++)
1576 : {
1577 37212 : if (F2x_degree(gel(V, k))==0) continue;
1578 14260 : gel(t,lfact) = gel(V, k);
1579 14260 : d = F2x_split_Berlekamp(&gel(t,lfact));
1580 14260 : if (flag == 2 && d != 1) return NULL;
1581 14106 : if (flag == 1)
1582 1533 : for (j=0; j<d; j++) t[lfact+j] = F2x_degree(gel(t,lfact+j));
1583 14106 : for (j=0; j<d; j++) E[lfact+j] = k;
1584 14106 : lfact += d;
1585 : }
1586 11244 : if (flag == 2) return gen_1; /* irreducible */
1587 11230 : setlg(t, lfact);
1588 11230 : setlg(E, lfact); y = mkvec2(t,E);
1589 11230 : return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
1590 11230 : : sort_factor_pol(y, cmpGuGu);
1591 : }
1592 :
1593 : /* Adapted from Shoup NTL */
1594 : GEN
1595 76273 : F2x_factor_squarefree(GEN f)
1596 : {
1597 : GEN r, t, v, tv;
1598 76273 : long i, q, n = F2x_degree(f);
1599 76271 : GEN u = const_vec(n+1, pol1_F2x(f[1]));
1600 134473 : for(q = 1;;q *= 2)
1601 : {
1602 134473 : r = F2x_gcd(f, F2x_deriv(f));
1603 134469 : if (F2x_degree(r) == 0)
1604 : {
1605 62564 : gel(u, q) = f;
1606 62564 : break;
1607 : }
1608 71908 : t = F2x_div(f, r);
1609 71909 : if (F2x_degree(t) > 0)
1610 : {
1611 : long j;
1612 68835 : for(j = 1;;j++)
1613 : {
1614 68835 : v = F2x_gcd(r, t);
1615 68836 : tv = F2x_div(t, v);
1616 68832 : if (F2x_degree(tv) > 0)
1617 32259 : gel(u, j*q) = tv;
1618 68834 : if (F2x_degree(v) <= 0) break;
1619 37840 : r = F2x_div(r, v);
1620 37841 : t = v;
1621 37841 : }
1622 30994 : if (F2x_degree(r) == 0) break;
1623 : }
1624 58199 : f = F2x_sqrt(r);
1625 58199 : }
1626 555947 : for (i = n; i; i--)
1627 555670 : if (F2x_degree(gel(u,i))) break;
1628 76290 : setlg(u,i+1); return u;
1629 : }
1630 :
1631 : static GEN
1632 80170 : F2x_ddf_simple(GEN T, GEN XP)
1633 : {
1634 80170 : pari_sp av = avma, av2;
1635 : GEN f, z, Tr, X;
1636 80170 : long j, n = F2x_degree(T), v = T[1], B = n/2;
1637 80170 : if (n == 0) return cgetg(1, t_VEC);
1638 80170 : if (n == 1) return mkvec(T);
1639 37937 : z = XP; Tr = T; X = polx_F2x(v);
1640 37938 : f = const_vec(n, pol1_F2x(v));
1641 37938 : av2 = avma;
1642 134346 : for (j = 1; j <= B; j++)
1643 : {
1644 100912 : GEN u = F2x_gcd(Tr, F2x_add(z, X));
1645 100913 : if (F2x_degree(u))
1646 : {
1647 21931 : gel(f, j) = u;
1648 21931 : Tr = F2x_div(Tr, u);
1649 21931 : av2 = avma;
1650 78896 : } else z = gerepileuptoleaf(av2, z);
1651 100931 : if (!F2x_degree(Tr)) break;
1652 96416 : z = F2xq_sqr(z, Tr);
1653 : }
1654 37937 : if (F2x_degree(Tr)) gel(f, F2x_degree(Tr)) = Tr;
1655 37937 : return gerepilecopy(av, f);
1656 : }
1657 :
1658 : GEN
1659 7 : F2x_ddf(GEN T)
1660 : {
1661 : GEN XP;
1662 7 : T = F2x_get_red(T);
1663 7 : XP = F2x_Frobenius(T);
1664 7 : return F2x_ddf_to_ddf2(F2x_ddf_simple(T, XP));
1665 : }
1666 :
1667 : static GEN
1668 6496 : F2xq_frobtrace(GEN a, long d, GEN T)
1669 : {
1670 6496 : pari_sp av = avma;
1671 : long i;
1672 6496 : GEN x = a;
1673 26617 : for(i=1; i<d; i++)
1674 : {
1675 20122 : x = F2x_add(a, F2xq_sqr(x,T));
1676 20121 : if (gc_needed(av, 2))
1677 0 : x = gerepileuptoleaf(av, x);
1678 : }
1679 6495 : return x;
1680 : }
1681 :
1682 : static void
1683 9719 : F2x_edf_simple(GEN Tp, GEN XP, long d, GEN V, long idx)
1684 : {
1685 9719 : long n = F2x_degree(Tp), r = n/d;
1686 : GEN T, f, ff;
1687 19438 : if (r==1) { gel(V, idx) = Tp; return; }
1688 3312 : T = Tp;
1689 3312 : XP = F2x_rem(XP, T);
1690 : while (1)
1691 : {
1692 6495 : pari_sp btop = avma;
1693 : long df;
1694 6495 : GEN g = random_F2x(n, Tp[1]);
1695 6495 : GEN t = F2xq_frobtrace(g, d, T);
1696 6495 : if (lgpol(t) == 0) continue;
1697 4888 : f = F2x_gcd(t, Tp); df = F2x_degree(f);
1698 4888 : if (df > 0 && df < n) break;
1699 1576 : avma = btop;
1700 3183 : }
1701 3312 : ff = F2x_div(Tp, f);
1702 3312 : F2x_edf_simple(f, XP, d, V, idx);
1703 3312 : F2x_edf_simple(ff, XP, d, V, idx+F2x_degree(f)/d);
1704 : }
1705 :
1706 : static GEN
1707 80163 : F2x_factor_Shoup(GEN T)
1708 : {
1709 80163 : long i, n, s = 0;
1710 : GEN XP, D, V;
1711 : pari_timer ti;
1712 80163 : n = F2x_degree(T);
1713 80163 : if (DEBUGLEVEL>=6) timer_start(&ti);
1714 80163 : XP = F2x_Frobenius(T);
1715 80163 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1716 80163 : D = F2x_ddf_simple(T, XP);
1717 80164 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1718 349207 : for (i = 1; i <= n; i++)
1719 269043 : s += F2x_degree(gel(D,i))/i;
1720 80164 : V = cgetg(s+1, t_COL);
1721 349208 : for (i = 1, s = 1; i <= n; i++)
1722 : {
1723 269045 : GEN Di = gel(D,i);
1724 269045 : long ni = F2x_degree(Di), ri = ni/i;
1725 269038 : if (ni == 0) continue;
1726 97586 : if (ni == i) { gel(V, s++) = Di; continue; }
1727 3095 : F2x_edf_simple(Di, XP, i, V, s);
1728 3095 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_edf(%ld)",i);
1729 3095 : s += ri;
1730 : }
1731 80163 : return V;
1732 : }
1733 :
1734 : static GEN
1735 64803 : F2x_factor_Cantor(GEN T)
1736 : {
1737 64803 : GEN E, F, V = F2x_factor_squarefree(T);
1738 64806 : long i, j, l = lg(V);
1739 64806 : E = cgetg(l, t_VEC);
1740 64806 : F = cgetg(l, t_VEC);
1741 255689 : for (i=1, j=1; i < l; i++)
1742 190883 : if (F2x_degree(gel(V,i)))
1743 : {
1744 80163 : GEN Fj = F2x_factor_Shoup(gel(V,i));
1745 80163 : gel(F, j) = Fj;
1746 80163 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1747 80163 : j++;
1748 : }
1749 64806 : return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
1750 : }
1751 :
1752 : #if 0
1753 : static GEN
1754 : F2x_simplefact_Shoup(GEN T)
1755 : {
1756 : long i, n, s = 0, j = 1, k;
1757 : GEN XP, D, V;
1758 : pari_timer ti;
1759 : n = F2x_degree(T);
1760 : if (DEBUGLEVEL>=6) timer_start(&ti);
1761 : XP = F2x_Frobenius(T);
1762 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1763 : D = F2x_ddf_simple(T, XP);
1764 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1765 : for (i = 1; i <= n; i++)
1766 : s += F2x_degree(gel(D,i))/i;
1767 : V = cgetg(s+1, t_VECSMALL);
1768 : for (i = 1; i <= n; i++)
1769 : {
1770 : long ni = F2x_degree(gel(D,i)), ri = ni/i;
1771 : if (ni == 0) continue;
1772 : for (k = 1; k <= ri; k++)
1773 : V[j++] = i;
1774 : }
1775 : return V;
1776 : }
1777 : static GEN
1778 : F2x_simplefact_Cantor(GEN T)
1779 : {
1780 : GEN E, F, V = F2x_factor_squarefree(T);
1781 : long i, j, l = lg(V);
1782 : F = cgetg(l, t_VEC);
1783 : E = cgetg(l, t_VEC);
1784 : for (i=1, j=1; i < l; i++)
1785 : if (F2x_degree(gel(V,i)))
1786 : {
1787 : GEN Fj = F2x_simplefact_Shoup(gel(V,i));
1788 : gel(F, j) = Fj;
1789 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1790 : j++;
1791 : }
1792 : return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
1793 : }
1794 : static int
1795 : F2x_isirred_Cantor(GEN T)
1796 : {
1797 : pari_sp av = avma;
1798 : pari_timer ti;
1799 : long n, d;
1800 : GEN dT = F2x_deriv(T);
1801 : GEN XP, D;
1802 : if (F2x_degree(F2x_gcd(T, dT)) != 0) { avma = av; return 0; }
1803 : n = F2x_degree(T);
1804 : if (DEBUGLEVEL>=6) timer_start(&ti);
1805 : XP = F2x_Frobenius(T);
1806 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1807 : D = F2x_ddf_simple(T, XP);
1808 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1809 : d = F2x_degree(gel(D, n));
1810 : avma = av; return d==n;
1811 : }
1812 : #endif
1813 :
1814 : /* driver for Cantor factorization, assume deg f > 2; not competitive for
1815 : * flag != 0, or as deg f increases */
1816 : static GEN
1817 64803 : F2x_Cantor_i(GEN f, long flag)
1818 : {
1819 : switch(flag)
1820 : {
1821 64803 : default: return F2x_factor_Cantor(f);
1822 : #if 0
1823 : case 1: return F2x_simplefact_Cantor(f);
1824 : case 2: return F2x_isirred_Cantor(f)? gen_1: NULL;
1825 : #endif
1826 : }
1827 : }
1828 : static GEN
1829 187840 : F2x_factor_i(GEN f, long flag)
1830 : {
1831 187840 : long d = F2x_degree(f);
1832 187837 : if (d <= 2) return F2x_factor_deg2(f,d,flag);
1833 150536 : return (flag == 0 && d <= 20)? F2x_Cantor_i(f, flag)
1834 141074 : : F2x_Berlekamp_i(f, flag);
1835 : }
1836 :
1837 : GEN
1838 0 : F2x_degfact(GEN f)
1839 : {
1840 0 : pari_sp av = avma;
1841 0 : GEN z = F2x_factor_i(f, 1);
1842 0 : return gerepilecopy(av, z);
1843 : }
1844 :
1845 : int
1846 238 : F2x_is_irred(GEN f) { return !!F2x_factor_i(f, 2); }
1847 :
1848 : /* Adapted from Shoup NTL */
1849 : GEN
1850 766898 : Flx_factor_squarefree(GEN f, ulong p)
1851 : {
1852 766898 : long i, q, n = degpol(f);
1853 766898 : GEN u = const_vec(n+1, pol1_Flx(f[1]));
1854 820958 : for(q = 1;;q *= p)
1855 : {
1856 820958 : GEN t, v, tv, r = Flx_gcd(f, Flx_deriv(f, p), p);
1857 820957 : if (degpol(r) == 0) { gel(u, q) = f; break; }
1858 98151 : t = Flx_div(f, r, p);
1859 98151 : if (degpol(t) > 0)
1860 : {
1861 : long j;
1862 143728 : for(j = 1;;j++)
1863 : {
1864 143728 : v = Flx_gcd(r, t, p);
1865 143728 : tv = Flx_div(t, v, p);
1866 143728 : if (degpol(tv) > 0)
1867 67666 : gel(u, j*q) = Flx_normalize(tv, p);
1868 143728 : if (degpol(v) <= 0) break;
1869 99071 : r = Flx_div(r, v, p);
1870 99071 : t = v;
1871 99071 : }
1872 44657 : if (degpol(r) == 0) break;
1873 : }
1874 54060 : f = Flx_normalize(Flx_deflate(r, p), p);
1875 54060 : }
1876 4157020 : for (i = n; i; i--)
1877 4157019 : if (degpol(gel(u,i))) break;
1878 766899 : setlg(u,i+1); return u;
1879 : }
1880 :
1881 : long
1882 126 : Flx_ispower(GEN f, ulong k, ulong p, GEN *pt_r)
1883 : {
1884 126 : pari_sp av = avma;
1885 : ulong lc;
1886 : GEN F;
1887 126 : long i, n = degpol(f), v = f[1], l;
1888 126 : if (n % k) return 0;
1889 126 : lc = Fl_sqrtn(Flx_lead(f), k, p, NULL);
1890 126 : if (lc == ULONG_MAX) { av = avma; return 0; }
1891 126 : F = Flx_factor_squarefree(f, p); l = lg(F)-1;
1892 6405 : for (i = 1; i <= l; i++)
1893 6300 : if (i%k && degpol(gel(F,i))) { avma = av; return 0; }
1894 105 : if (pt_r)
1895 : {
1896 105 : GEN r = Fl_to_Flx(lc, v), s = pol1_Flx(v);
1897 6384 : for(i = l; i >= 1; i--)
1898 : {
1899 6279 : if (i%k) continue;
1900 1274 : s = Flx_mul(s, gel(F,i), p);
1901 1274 : r = Flx_mul(r, s, p);
1902 : }
1903 105 : *pt_r = gerepileuptoleaf(av, r);
1904 0 : } else av = avma;
1905 105 : return 1;
1906 : }
1907 :
1908 : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
1909 : static GEN
1910 1070307 : Flx_ddf_Shoup(GEN T, GEN XP, ulong p)
1911 : {
1912 1070307 : pari_sp av = avma;
1913 : GEN b, g, h, F, f, Tr, xq;
1914 : long i, j, n, v, bo, ro;
1915 : long B, l, m;
1916 : pari_timer ti;
1917 1070307 : n = get_Flx_degree(T); v = get_Flx_var(T);
1918 1070307 : if (n == 0) return cgetg(1, t_VEC);
1919 1068292 : if (n == 1) return mkvec(get_Flx_mod(T));
1920 973144 : B = n/2;
1921 973144 : l = usqrt(B);
1922 973144 : m = (B+l-1)/l;
1923 973144 : T = Flx_get_red(T, p);
1924 973144 : b = cgetg(l+2, t_VEC);
1925 973144 : gel(b, 1) = polx_Flx(v);
1926 973144 : gel(b, 2) = XP;
1927 973144 : bo = brent_kung_optpow(n, l-1, 1);
1928 973144 : ro = l<=1 ? 0:(bo-1)/(l-1) + ((n-1)/bo);
1929 973144 : if (DEBUGLEVEL>=7) timer_start(&ti);
1930 973144 : if (expu(p) <= ro)
1931 227824 : for (i = 3; i <= l+1; i++)
1932 128661 : gel(b, i) = Flxq_powu(gel(b, i-1), p, T, p);
1933 : else
1934 : {
1935 873981 : xq = Flxq_powers(gel(b, 2), bo, T, p);
1936 873981 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq baby");
1937 1017822 : for (i = 3; i <= l+1; i++)
1938 143841 : gel(b, i) = Flx_FlxqV_eval(gel(b, i-1), xq, T, p);
1939 : }
1940 973144 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: baby");
1941 973144 : xq = Flxq_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1), T, p);
1942 973144 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq giant");
1943 973144 : g = cgetg(m+1, t_VEC);
1944 973144 : gel(g, 1) = gel(xq, 2);
1945 1816740 : for(i = 2; i <= m; i++)
1946 843596 : gel(g, i) = Flx_FlxqV_eval(gel(g, i-1), xq, T, p);
1947 973144 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: giant");
1948 973144 : h = cgetg(m+1, t_VEC);
1949 2789884 : for (j = 1; j <= m; j++)
1950 : {
1951 1816740 : pari_sp av = avma;
1952 1816740 : GEN gj = gel(g, j);
1953 1816740 : GEN e = Flx_sub(gj, gel(b, 1), p);
1954 2835041 : for (i = 2; i <= l; i++)
1955 1018301 : e = Flxq_mul(e, Flx_sub(gj, gel(b, i), p), T, p);
1956 1816740 : gel(h, j) = gerepileupto(av, e);
1957 : }
1958 973144 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: diff");
1959 973144 : Tr = get_Flx_mod(T);
1960 973144 : F = cgetg(m+1, t_VEC);
1961 2789884 : for (j = 1; j <= m; j++)
1962 : {
1963 1816740 : GEN u = Flx_gcd(Tr, gel(h, j), p);
1964 1816739 : if (degpol(u))
1965 : {
1966 745452 : u = Flx_normalize(u, p);
1967 745452 : Tr = Flx_div(Tr, u, p);
1968 : }
1969 1816739 : gel(F, j) = u;
1970 : }
1971 973144 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: F");
1972 973144 : f = const_vec(n, pol1_Flx(v));
1973 2789882 : for (j = 1; j <= m; j++)
1974 : {
1975 1816738 : GEN e = gel(F, j);
1976 2000589 : for (i=l-1; i >= 0; i--)
1977 : {
1978 2000589 : GEN u = Flx_gcd(e, Flx_sub(gel(g, j), gel(b, i+1), p), p);
1979 2000590 : if (degpol(u))
1980 : {
1981 764537 : gel(f, l*j-i) = u;
1982 764537 : e = Flx_div(e, u, p);
1983 : }
1984 2000591 : if (!degpol(e)) break;
1985 : }
1986 : }
1987 973144 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: f");
1988 973144 : if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
1989 973144 : return gerepilecopy(av, f);
1990 : }
1991 :
1992 : static void
1993 121141 : Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
1994 : {
1995 121141 : long n = degpol(Tp), r = n/d;
1996 : GEN T, f, ff;
1997 : ulong p2;
1998 242282 : if (r==1) { gel(V, idx) = Tp; return; }
1999 52967 : p2 = p>>1;
2000 52967 : T = Flx_get_red(Tp, p);
2001 52967 : XP = Flx_rem(XP, T, p);
2002 : while (1)
2003 : {
2004 58135 : pari_sp btop = avma;
2005 : long i;
2006 58135 : GEN g = random_Flx(n, Tp[1], p);
2007 58135 : GEN t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
2008 58135 : if (lgpol(t) == 0) continue;
2009 122751 : for(i=1; i<=10; i++)
2010 : {
2011 118853 : pari_sp btop2 = avma;
2012 118853 : GEN R = Flxq_powu(Flx_Fl_add(t, random_Fl(p), p), p2, T, p);
2013 118853 : f = Flx_gcd(Flx_Fl_add(R, p-1, p), Tp, p);
2014 118853 : if (degpol(f) > 0 && degpol(f) < n) break;
2015 65886 : avma = btop2;
2016 : }
2017 56865 : if (degpol(f) > 0 && degpol(f) < n) break;
2018 3898 : avma = btop;
2019 5168 : }
2020 52967 : f = Flx_normalize(f, p);
2021 52967 : ff = Flx_div(Tp, f ,p);
2022 52967 : Flx_edf_simple(f, XP, d, p, V, idx);
2023 52967 : Flx_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
2024 : }
2025 : static void
2026 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx);
2027 :
2028 : static void
2029 342534 : Flx_edf_rec(GEN T, GEN XP, GEN hp, GEN t, long d, ulong p, GEN V, long idx)
2030 : {
2031 : pari_sp av;
2032 342534 : GEN Tp = get_Flx_mod(T);
2033 342534 : long n = degpol(hp), vT = Tp[1];
2034 : GEN u1, u2, f1, f2;
2035 342534 : ulong p2 = p>>1;
2036 : GEN R, h;
2037 342534 : h = Flx_get_red(hp, p);
2038 342534 : t = Flx_rem(t, T, p);
2039 342534 : av = avma;
2040 : do
2041 : {
2042 560568 : avma = av;
2043 560568 : R = Flxq_powu(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p);
2044 560568 : u1 = Flx_gcd(Flx_Fl_add(R, p-1, p), hp, p);
2045 560568 : } while (degpol(u1)==0 || degpol(u1)==n);
2046 342534 : f1 = Flx_gcd(Flx_Flxq_eval(u1, t, T, p), Tp, p);
2047 342534 : f1 = Flx_normalize(f1, p);
2048 342534 : u2 = Flx_div(hp, u1, p);
2049 342534 : f2 = Flx_div(Tp, f1, p);
2050 342534 : if (degpol(u1)==1)
2051 : {
2052 248074 : if (degpol(f1)==d)
2053 244046 : gel(V, idx) = f1;
2054 : else
2055 4028 : Flx_edf(f1, XP, d, p, V, idx);
2056 : }
2057 : else
2058 94460 : Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, V, idx);
2059 342534 : idx += degpol(f1)/d;
2060 342534 : if (degpol(u2)==1)
2061 : {
2062 245709 : if (degpol(f2)==d)
2063 241994 : gel(V, idx) = f2;
2064 : else
2065 3715 : Flx_edf(f2, XP, d, p, V, idx);
2066 : }
2067 : else
2068 96825 : Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, V, idx);
2069 342534 : }
2070 :
2071 : static void
2072 151249 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
2073 : {
2074 151249 : long n = degpol(Tp), r = n/d, vT = Tp[1];
2075 : GEN T, h, t;
2076 : pari_timer ti;
2077 302498 : if (r==1) { gel(V, idx) = Tp; return; }
2078 151249 : T = Flx_get_red(Tp, p);
2079 151249 : XP = Flx_rem(XP, T, p);
2080 151249 : if (DEBUGLEVEL>=7) timer_start(&ti);
2081 : do
2082 : {
2083 156261 : GEN g = random_Flx(n, vT, p);
2084 156261 : t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
2085 156261 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
2086 156261 : h = Flxq_minpoly(t, T, p);
2087 156261 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
2088 156261 : } while (degpol(h) <= 1);
2089 151249 : Flx_edf_rec(T, XP, h, t, d, p, V, idx);
2090 : }
2091 :
2092 : static GEN
2093 433459 : Flx_factor_Shoup(GEN T, ulong p)
2094 : {
2095 433459 : long i, n, s = 0;
2096 : GEN XP, D, V;
2097 433459 : long e = expu(p);
2098 : pari_timer ti;
2099 433459 : n = get_Flx_degree(T);
2100 433459 : T = Flx_get_red(T, p);
2101 433459 : if (DEBUGLEVEL>=6) timer_start(&ti);
2102 433459 : XP = Flx_Frobenius(T, p);
2103 433459 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2104 433459 : D = Flx_ddf_Shoup(T, XP, p);
2105 433459 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2106 433459 : s = ddf_to_nbfact(D);
2107 433459 : V = cgetg(s+1, t_COL);
2108 2415966 : for (i = 1, s = 1; i <= n; i++)
2109 : {
2110 1982507 : GEN Di = gel(D,i);
2111 1982507 : long ni = degpol(Di), ri = ni/i;
2112 1982507 : if (ni == 0) continue;
2113 542154 : Di = Flx_normalize(Di, p);
2114 542154 : if (ni == i) { gel(V, s++) = Di; continue; }
2115 158713 : if (ri <= e*expu(e))
2116 143506 : Flx_edf(Di, XP, i, p, V, s);
2117 : else
2118 15207 : Flx_edf_simple(Di, XP, i, p, V, s);
2119 158713 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
2120 158713 : s += ri;
2121 : }
2122 433459 : return V;
2123 : }
2124 :
2125 : static GEN
2126 411384 : Flx_factor_Cantor(GEN T, ulong p)
2127 : {
2128 411384 : GEN E, F, V = Flx_factor_squarefree(get_Flx_mod(T), p);
2129 411384 : long i, j, l = lg(V);
2130 411384 : F = cgetg(l, t_VEC);
2131 411384 : E = cgetg(l, t_VEC);
2132 1037992 : for (i=1, j=1; i < l; i++)
2133 626608 : if (degpol(gel(V,i)))
2134 : {
2135 433459 : GEN Fj = Flx_factor_Shoup(gel(V,i), p);
2136 433459 : gel(F, j) = Fj;
2137 433459 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
2138 433459 : j++;
2139 : }
2140 411384 : return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
2141 : }
2142 :
2143 : GEN
2144 0 : Flx_ddf(GEN T, ulong p)
2145 : {
2146 : GEN XP;
2147 0 : T = Flx_get_red(T, p);
2148 0 : XP = Flx_Frobenius(T, p);
2149 0 : return ddf_to_ddf2(Flx_ddf_Shoup(T, XP, p));
2150 : }
2151 :
2152 : static GEN
2153 354639 : Flx_simplefact_Cantor(GEN T, ulong p)
2154 : {
2155 : GEN XP, V;
2156 : long i, l;
2157 354639 : T = Flx_get_red(T, p);
2158 354639 : XP = Flx_Frobenius(T, p);
2159 354639 : V = Flx_factor_squarefree(get_Flx_mod(T), p); l = lg(V);
2160 354639 : for (i=1; i < l; i++) gel(V,i) = Flx_ddf_Shoup(gel(V,i), XP, p);
2161 354639 : return vddf_to_simplefact(V, get_Flx_degree(T));
2162 : }
2163 :
2164 : static int
2165 581 : Flx_isirred_Cantor(GEN Tp, ulong p)
2166 : {
2167 581 : pari_sp av = avma;
2168 : pari_timer ti;
2169 : long n, d;
2170 581 : GEN T = get_Flx_mod(Tp);
2171 581 : GEN dT = Flx_deriv(T, p);
2172 : GEN XP, D;
2173 581 : if (degpol(Flx_gcd(T, dT, p)) != 0) { avma = av; return 0; }
2174 441 : n = get_Flx_degree(T);
2175 441 : T = Flx_get_red(Tp, p);
2176 441 : if (DEBUGLEVEL>=6) timer_start(&ti);
2177 441 : XP = Flx_Frobenius(T, p);
2178 441 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2179 441 : D = Flx_ddf_Shoup(T, XP, p);
2180 441 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2181 441 : d = degpol(gel(D, n));
2182 441 : avma = av; return d==n;
2183 : }
2184 :
2185 : /* f monic */
2186 : static GEN
2187 1960677 : Flx_factor_i(GEN f, ulong pp, long flag)
2188 : {
2189 : long d;
2190 1960677 : if (pp==2) { /*We need to handle 2 specially */
2191 28847 : GEN F = F2x_factor_i(Flx_to_F2x(f),flag);
2192 28847 : if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
2193 28847 : return F;
2194 : }
2195 1931830 : d = degpol(f);
2196 1931830 : if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
2197 766604 : switch(flag)
2198 : {
2199 411384 : default: return Flx_factor_Cantor(f, pp);
2200 354639 : case 1: return Flx_simplefact_Cantor(f, pp);
2201 581 : case 2: return Flx_isirred_Cantor(f, pp)? gen_1: NULL;
2202 : }
2203 : }
2204 :
2205 : GEN
2206 1106201 : Flx_degfact(GEN f, ulong p)
2207 : {
2208 1106201 : pari_sp av = avma;
2209 1106201 : GEN z = Flx_factor_i(Flx_normalize(f,p),p,1);
2210 1106201 : return gerepilecopy(av, z);
2211 : }
2212 :
2213 : /* T must be squarefree mod p*/
2214 : GEN
2215 188023 : Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
2216 : {
2217 : GEN XP, D;
2218 : pari_timer ti;
2219 188023 : long i, s, n = get_Flx_degree(T);
2220 188023 : GEN V = const_vecsmall(n, 0);
2221 188023 : pari_sp av = avma;
2222 188023 : T = Flx_get_red(T, p);
2223 188023 : if (DEBUGLEVEL>=6) timer_start(&ti);
2224 188023 : XP = Flx_Frobenius(T, p);
2225 188023 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2226 188023 : D = Flx_ddf_Shoup(T, XP, p);
2227 188023 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2228 1208547 : for (i = 1, s = 0; i <= n; i++)
2229 : {
2230 1020524 : V[i] = degpol(gel(D,i))/i;
2231 1020524 : s += V[i];
2232 : }
2233 188023 : *nb = s;
2234 188023 : avma = av; return V;
2235 : }
2236 :
2237 : long
2238 90741 : Flx_nbfact_Frobenius(GEN T, GEN XP, ulong p)
2239 : {
2240 90741 : pari_sp av = avma;
2241 90741 : long s = ddf_to_nbfact(Flx_ddf_Shoup(T, XP, p));
2242 90741 : avma = av; return s;
2243 : }
2244 :
2245 : /* T must be squarefree mod p*/
2246 : long
2247 90741 : Flx_nbfact(GEN T, ulong p)
2248 : {
2249 90741 : pari_sp av = avma;
2250 90741 : GEN XP = Flx_Frobenius(T, p);
2251 90741 : long n = Flx_nbfact_Frobenius(T, XP, p);
2252 90741 : avma = av; return n;
2253 : }
2254 :
2255 : int
2256 581 : Flx_is_irred(GEN f, ulong p)
2257 : {
2258 581 : pari_sp av = avma;
2259 581 : int z = !!Flx_factor_i(Flx_normalize(f,p),p,2);
2260 581 : avma = av; return z;
2261 : }
2262 :
2263 : /* Use this function when you think f is reducible, and that there are lots of
2264 : * factors. If you believe f has few factors, use FpX_nbfact(f,p)==1 instead */
2265 : int
2266 49 : FpX_is_irred(GEN f, GEN p)
2267 : {
2268 49 : pari_sp av = avma;
2269 : int z;
2270 49 : switch(ZX_factmod_init(&f,p))
2271 : {
2272 7 : case 0: z = !!F2x_factor_i(f,2); break;
2273 35 : case 1: z = !!Flx_factor_i(f,p[2],2); break;
2274 7 : default: z = !!FpX_factor_i(f,p,2); break;
2275 : }
2276 49 : avma = av; return z;
2277 : }
2278 : GEN
2279 42 : FpX_degfact(GEN f, GEN p) {
2280 42 : pari_sp av = avma;
2281 : GEN F;
2282 42 : switch(ZX_factmod_init(&f,p))
2283 : {
2284 7 : case 0: F = F2x_factor_i(f,1); break;
2285 7 : case 1: F = Flx_factor_i(f,p[2],1); break;
2286 28 : default: F = FpX_factor_i(f,p,1); break;
2287 : }
2288 42 : return gerepilecopy(av, F);
2289 : }
2290 :
2291 : #if 0
2292 : /* set x <-- x + c*y mod p */
2293 : /* x is not required to be normalized.*/
2294 : static void
2295 : Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
2296 : {
2297 : long i, lx, ly;
2298 : ulong *x=(ulong *)gx;
2299 : ulong *y=(ulong *)gy;
2300 : if (!c) return;
2301 : lx = lg(gx);
2302 : ly = lg(gy);
2303 : if (lx<ly) pari_err_BUG("lx<ly in Flx_addmul_inplace");
2304 : if (SMALL_ULONG(p))
2305 : for (i=2; i<ly; i++) x[i] = (x[i] + c*y[i]) % p;
2306 : else
2307 : for (i=2; i<ly; i++) x[i] = Fl_add(x[i], Fl_mul(c,y[i],p),p);
2308 : }
2309 : #endif
2310 :
2311 : GEN
2312 863473 : FpX_factor(GEN f, GEN p)
2313 : {
2314 863473 : pari_sp av = avma;
2315 : GEN F;
2316 863473 : switch(ZX_factmod_init(&f, p))
2317 : {
2318 143799 : case 0: F = F2x_factor_i(f,0);
2319 143799 : F2xV_to_ZXV_inplace(gel(F,1)); break;
2320 718325 : case 1: F = Flx_factor_i(f,p[2],0);
2321 718325 : FlxV_to_ZXV_inplace(gel(F,1)); break;
2322 1349 : default: F = FpX_factor_i(f,p,0); break;
2323 : }
2324 863473 : return gerepilecopy(av, F);
2325 : }
2326 :
2327 : GEN
2328 135528 : Flx_factor(GEN f, ulong p)
2329 : {
2330 135528 : pari_sp av = avma;
2331 135528 : return gerepilecopy(av, Flx_factor_i(f,p,0));
2332 : }
2333 : GEN
2334 14942 : F2x_factor(GEN f)
2335 : {
2336 14942 : pari_sp av = avma;
2337 14942 : return gerepilecopy(av, F2x_factor_i(f,0));
2338 : }
|
