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 1016165 : ZX_factmod_init(GEN *F, GEN p)
32 : {
33 1016165 : if (lgefint(p) == 3)
34 : {
35 1014677 : ulong pp = p[2];
36 1014677 : if (pp == 2) { *F = ZX_to_F2x(*F); return 0; }
37 847810 : *F = ZX_to_Flx(*F, pp);
38 847820 : if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
39 847817 : return 1;
40 : }
41 1488 : *F = FpX_red(*F, p);
42 1489 : if (lg(*F) > 3) *F = FpX_normalize(*F, p);
43 1489 : return 2;
44 : }
45 : static void
46 86137 : ZX_rootmod_init(GEN *F, GEN p)
47 : {
48 86137 : if (lgefint(p) == 3)
49 : {
50 78166 : ulong pp = p[2];
51 78166 : *F = ZX_to_Flx(*F, pp);
52 78166 : if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
53 : }
54 : else
55 : {
56 7971 : *F = FpX_red(*F, p);
57 7971 : if (lg(*F) > 3) *F = FpX_normalize(*F, p);
58 : }
59 86137 : }
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 84 : for (i = 1; i < p; i++) r[i] = i;
70 : } else {
71 14 : r = cgetg(p+1, t_VECSMALL);
72 84 : 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 203 : Flx_Xnm1(long sv, long n, ulong p)
80 : {
81 203 : GEN t = cgetg(n+3, t_VECSMALL);
82 : long i;
83 203 : t[1] = sv;
84 203 : t[2] = p - 1;
85 581 : for (i = 3; i <= n+1; i++) t[i] = 0;
86 203 : t[i] = 1; return t;
87 : }
88 : /* X^n + 1 */
89 : static GEN
90 364 : Flx_Xn1(long sv, long n, ulong p)
91 : {
92 364 : GEN t = cgetg(n+3, t_VECSMALL);
93 : long i;
94 : (void) p;
95 364 : t[1] = sv;
96 364 : t[2] = 1;
97 1036 : for (i = 3; i <= n+1; i++) t[i] = 0;
98 364 : t[i] = 1; return t;
99 : }
100 :
101 : static GEN
102 11505 : Flx_root_mod_2(GEN f)
103 : {
104 11505 : int z1, z0 = !(f[2] & 1);
105 : long i,n;
106 : GEN y;
107 :
108 44019 : for (i=2, n=1; i < lg(f); i++) n += f[i];
109 11505 : z1 = n & 1;
110 11505 : y = cgetg(z0+z1+1, t_VECSMALL); i = 1;
111 11505 : if (z0) y[i++] = 0;
112 11505 : if (z1) y[i ] = 1;
113 11505 : return y;
114 : }
115 : static ulong
116 70 : Flx_oneroot_mod_2(GEN f)
117 : {
118 : long i,n;
119 70 : if (!(f[2] & 1)) return 0;
120 280 : for (i=2, n=1; i < lg(f); i++) n += f[i];
121 70 : if (n & 1) return 1;
122 28 : return 2;
123 : }
124 :
125 : static GEN FpX_roots_i(GEN f, GEN p);
126 : static GEN Flx_roots_i(GEN f, ulong p);
127 :
128 : static int
129 4394266 : cmpGuGu(GEN a, GEN b) { return (ulong)a < (ulong)b? -1: (a == b? 0: 1); }
130 :
131 : /* Generic driver to computes the roots of f modulo pp, using 'Roots' when
132 : * pp is a small prime.
133 : * if (gpwrap), check types thoroughly and return t_INTMODs, otherwise
134 : * assume that f is an FpX, pp a prime and return t_INTs */
135 : static GEN
136 73481 : rootmod_aux(GEN f, GEN pp)
137 : {
138 : GEN y;
139 73481 : switch(lg(f))
140 : {
141 14 : case 2: pari_err_ROOTS0("rootmod");
142 49 : case 3: return cgetg(1,t_COL);
143 : }
144 73418 : if (typ(f) == t_VECSMALL)
145 : {
146 70308 : ulong p = pp[2];
147 70308 : if (p == 2)
148 11505 : y = Flx_root_mod_2(f);
149 : else
150 : {
151 58803 : if (!odd(p)) pari_err_PRIME("rootmod",utoi(p));
152 58803 : y = Flx_roots_i(f, p);
153 : }
154 70301 : y = Flc_to_ZC(y);
155 : }
156 : else
157 3110 : y = FpX_roots_i(f, pp);
158 73404 : return y;
159 : }
160 : /* assume that f is a ZX and p a prime */
161 : GEN
162 73481 : FpX_roots(GEN f, GEN p)
163 : {
164 73481 : pari_sp av = avma;
165 73481 : GEN y; ZX_rootmod_init(&f, p); y = rootmod_aux(f, p);
166 73453 : return gerepileupto(av, y);
167 : }
168 :
169 : /* assume x reduced mod p > 2, monic. */
170 : static int
171 21 : FpX_quad_factortype(GEN x, GEN p)
172 : {
173 21 : GEN b = gel(x,3), c = gel(x,2);
174 21 : GEN D = subii(sqri(b), shifti(c,2));
175 21 : return kronecker(D,p);
176 : }
177 : /* assume x reduced mod p, monic. Return one root, or NULL if irreducible */
178 : static GEN
179 7566 : FpX_quad_root(GEN x, GEN p, int unknown)
180 : {
181 7566 : GEN s, D, b = gel(x,3), c = gel(x,2);
182 :
183 7566 : if (absequaliu(p, 2)) {
184 0 : if (!signe(b)) return c;
185 0 : return signe(c)? NULL: gen_1;
186 : }
187 7566 : D = subii(sqri(b), shifti(c,2));
188 7566 : D = remii(D,p);
189 7566 : if (unknown && kronecker(D,p) == -1) return NULL;
190 :
191 7025 : s = Fp_sqrt(D,p);
192 : /* p is not prime, go on and give e.g. maxord a chance to recover */
193 7025 : if (!s) return NULL;
194 7017 : return Fp_halve(Fp_sub(s,b, p), p);
195 : }
196 : static GEN
197 3152 : FpX_otherroot(GEN x, GEN r, GEN p)
198 3152 : { return Fp_neg(Fp_add(gel(x,3), r, p), p); }
199 :
200 : /* disc(x^2+bx+c) = b^2 - 4c */
201 : static ulong
202 21584767 : Fl_disc_bc(ulong b, ulong c, ulong p)
203 21584767 : { return Fl_sub(Fl_sqr(b,p), Fl_double(Fl_double(c,p),p), p); }
204 : /* p > 2 */
205 : static ulong
206 20456246 : Flx_quad_root(GEN x, ulong p, int unknown)
207 : {
208 20456246 : ulong s, b = x[3], c = x[2];
209 20456246 : ulong D = Fl_disc_bc(b, c, p);
210 20439092 : if (unknown && krouu(D,p) == -1) return p;
211 13649368 : s = Fl_sqrt(D,p);
212 13711376 : if (s==~0UL) return p;
213 13711363 : return Fl_halve(Fl_sub(s,b, p), p);
214 : }
215 : static ulong
216 12274073 : Flx_otherroot(GEN x, ulong r, ulong p)
217 12274073 : { return Fl_neg(Fl_add(x[3], r, p), p); }
218 :
219 : /* 'todo' contains the list of factors to be split.
220 : * 'done' the list of finished factors, no longer touched */
221 : struct split_t { GEN todo, done; };
222 : static void
223 4931171 : split_init(struct split_t *S, long max)
224 : {
225 4931171 : S->todo = vectrunc_init(max);
226 4930802 : S->done = vectrunc_init(max);
227 4930438 : }
228 : #if 0
229 : /* move todo[i] to done */
230 : static void
231 : split_convert(struct split_t *S, long i)
232 : {
233 : long n = lg(S->todo)-1;
234 : vectrunc_append(S->done, gel(S->todo,i));
235 : if (n) gel(S->todo,i) = gel(S->todo, n);
236 : setlg(S->todo, n);
237 : }
238 : #endif
239 : /* append t to todo */
240 : static void
241 5174123 : split_add(struct split_t *S, GEN t) { vectrunc_append(S->todo, t); }
242 : /* delete todo[i], add t to done */
243 : static void
244 5175374 : split_moveto_done(struct split_t *S, long i, GEN t)
245 : {
246 5175374 : long n = lg(S->todo)-1;
247 5175374 : vectrunc_append(S->done, t);
248 5175603 : if (n) gel(S->todo,i) = gel(S->todo, n);
249 5175603 : setlg(S->todo, n);
250 :
251 5175590 : }
252 : /* append t to done */
253 : static void
254 395376 : split_add_done(struct split_t *S, GEN t)
255 395376 : { vectrunc_append(S->done, t); }
256 : /* split todo[i] into a and b */
257 : static void
258 338598 : split_todo(struct split_t *S, long i, GEN a, GEN b)
259 : {
260 338598 : gel(S->todo, i) = a;
261 338598 : split_add(S, b);
262 338598 : }
263 : /* split todo[i] into a and b, moved to done */
264 : static void
265 374843 : split_done(struct split_t *S, long i, GEN a, GEN b)
266 : {
267 374843 : split_moveto_done(S, i, a);
268 374846 : split_add_done(S, b);
269 374845 : }
270 :
271 : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
272 : static GEN
273 3110 : FpX_roots_i(GEN f, GEN p)
274 : {
275 : GEN pol, pol0, a, q;
276 : struct split_t S;
277 :
278 3110 : split_init(&S, lg(f)-1);
279 3110 : settyp(S.done, t_COL);
280 3110 : if (ZX_valrem(f, &f)) split_add_done(&S, gen_0);
281 3110 : switch(degpol(f))
282 : {
283 7 : case 0: return ZC_copy(S.done);
284 14 : case 1: split_add_done(&S, subii(p, gel(f,2))); return ZC_copy(S.done);
285 1744 : case 2: {
286 1744 : GEN s, r = FpX_quad_root(f, p, 1);
287 1744 : if (r) {
288 1744 : split_add_done(&S, r);
289 1744 : s = FpX_otherroot(f,r, p);
290 : /* f not known to be square free yet */
291 1744 : if (!equalii(r, s)) split_add_done(&S, s);
292 : }
293 1744 : return sort(S.done);
294 : }
295 : }
296 :
297 1345 : a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
298 1345 : if (lg(a) < 3) pari_err_PRIME("rootmod",p);
299 1345 : a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
300 1345 : a = FpX_gcd(f,a, p);
301 1345 : if (!degpol(a)) return ZC_copy(S.done);
302 1345 : split_add(&S, FpX_normalize(a,p));
303 :
304 1345 : q = shifti(p,-1);
305 1345 : pol0 = icopy(gen_1); /* constant term, will vary in place */
306 1345 : pol = deg1pol_shallow(gen_1, pol0, varn(f));
307 1345 : for (pol0[2] = 1;; pol0[2]++)
308 1475 : {
309 2820 : long j, l = lg(S.todo);
310 2820 : if (l == 1) return sort(S.done);
311 1482 : if (pol0[2] == 100 && !BPSW_psp(p)) pari_err_PRIME("polrootsmod",p);
312 3080 : for (j = 1; j < l; j++)
313 : {
314 1605 : GEN b, r, s, c = gel(S.todo,j);
315 1605 : switch(degpol(c))
316 : { /* convert linear and quadratics to roots, try to split the rest */
317 196 : case 1:
318 196 : split_moveto_done(&S, j, subii(p, gel(c,2)));
319 196 : j--; l--; break;
320 1272 : case 2:
321 1272 : r = FpX_quad_root(c, p, 0);
322 1272 : if (!r) pari_err_PRIME("polrootsmod",p);
323 1265 : s = FpX_otherroot(c,r, p);
324 1265 : split_done(&S, j, r, s);
325 1265 : j--; l--; break;
326 137 : default:
327 137 : b = FpXQ_pow(pol,q, c,p);
328 137 : if (degpol(b) <= 0) continue;
329 123 : b = FpX_gcd(c,FpX_Fp_sub_shallow(b,gen_1,p), p);
330 123 : if (!degpol(b)) continue;
331 123 : b = FpX_normalize(b, p);
332 123 : c = FpX_div(c,b, p);
333 123 : split_todo(&S, j, b, c);
334 : }
335 : }
336 : }
337 : }
338 :
339 : /* Assume f is normalized */
340 : static ulong
341 174417 : Flx_cubic_root(GEN ff, ulong p)
342 : {
343 174417 : GEN f = Flx_normalize(ff,p);
344 174417 : ulong pi = get_Fl_red(p);
345 174417 : ulong a = f[4], b=f[3], c=f[2], p3 = p%3==1 ? (2*p+1)/3 :(p+1)/3;
346 174417 : ulong t = Fl_mul_pre(a, p3, p, pi), t2 = Fl_sqr_pre(t, p, pi);
347 174418 : ulong A = Fl_sub(b, Fl_triple(t2, p), p);
348 174418 : ulong B = Fl_addmul_pre(c, t, Fl_sub(Fl_double(t2, p), b, p), p, pi);
349 174418 : ulong A3 = Fl_mul_pre(A, p3, p, pi);
350 174418 : ulong A32 = Fl_sqr_pre(A3, p, pi), A33 = Fl_mul_pre(A3, A32, p, pi);
351 174418 : ulong S = Fl_neg(B,p), P = Fl_neg(A3,p);
352 174418 : ulong D = Fl_add(Fl_sqr_pre(S, p, pi), Fl_double(Fl_double(A33, p), p), p);
353 174416 : ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
354 174417 : if (s!=~0UL)
355 : {
356 109962 : ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
357 109962 : if (p%3==2) /* 1 solutions */
358 19446 : vS1 = Fl_powu_pre(S1, (2*p-1)/3, p, pi);
359 : else
360 : {
361 90516 : vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
362 90516 : if (vS1==~0UL) return p; /*0 solutions*/
363 : /*3 solutions*/
364 : }
365 79453 : vS2 = P? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): 0;
366 79453 : return Fl_sub(Fl_add(vS1,vS2, p), t, p);
367 : }
368 : else
369 : {
370 64455 : pari_sp av = avma;
371 64455 : GEN S1 = mkvecsmall2(Fl_halve(S, p), Fl_halve(1UL, p));
372 64455 : GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
373 : ulong Sa;
374 64456 : if (!vS1) return p; /*0 solutions, p%3==2*/
375 64456 : Sa = vS1[1];
376 64456 : if (p%3==1) /*1 solutions*/
377 : {
378 23992 : ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
379 23992 : if (Fa!=P)
380 15973 : Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
381 : }
382 64456 : set_avma(av);
383 64456 : return Fl_sub(Fl_double(Sa,p),t,p);
384 : }
385 : }
386 :
387 : /* assume p > 2 prime */
388 : static ulong
389 3140287 : Flx_oneroot_i(GEN f, ulong p, long fl)
390 : {
391 : GEN pol, a;
392 : ulong q;
393 : long da;
394 :
395 3140287 : if (Flx_val(f)) return 0;
396 3139428 : switch(degpol(f))
397 : {
398 12370 : case 1: return Fl_neg(f[2], p);
399 2689245 : case 2: return Flx_quad_root(f, p, 1);
400 158712 : case 3: if (p>3) return Flx_cubic_root(f, p); /*FALL THROUGH*/
401 : }
402 :
403 279108 : if (!fl)
404 : {
405 247689 : a = Flxq_powu(polx_Flx(f[1]), p - 1, f,p);
406 247689 : if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
407 247689 : a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
408 247689 : a = Flx_gcd(f,a, p);
409 31419 : } else a = f;
410 279108 : da = degpol(a);
411 279108 : if (!da) return p;
412 199426 : a = Flx_normalize(a,p);
413 :
414 199426 : q = p >> 1;
415 199426 : pol = polx_Flx(f[1]);
416 199424 : for(pol[2] = 1;; pol[2]++)
417 : {
418 299810 : if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
419 299811 : switch(da)
420 : {
421 124534 : case 1: return Fl_neg(a[2], p);
422 59189 : case 2: return Flx_quad_root(a, p, 0);
423 15712 : case 3: if (p>3) return Flx_cubic_root(a, p); /*FALL THROUGH*/
424 : default: {
425 100376 : GEN b = Flxq_powu(pol,q, a,p);
426 : long db;
427 100371 : if (degpol(b) <= 0) continue;
428 95330 : b = Flx_gcd(a,Flx_Fl_add(b,p-1,p), p);
429 95329 : db = degpol(b); if (!db) continue;
430 95323 : b = Flx_normalize(b, p);
431 95325 : if (db <= (da >> 1)) {
432 58163 : a = b;
433 58163 : da = db;
434 : } else {
435 37162 : a = Flx_div(a,b, p);
436 37175 : da -= db;
437 : }
438 : }
439 : }
440 : }
441 : }
442 :
443 : /* assume p > 2 prime */
444 : static GEN
445 4847 : FpX_oneroot_i(GEN f, GEN p)
446 : {
447 : GEN pol, pol0, a, q;
448 : long da;
449 :
450 4847 : if (ZX_val(f)) return gen_0;
451 4581 : switch(degpol(f))
452 : {
453 674 : case 1: return subii(p, gel(f,2));
454 3837 : case 2: return FpX_quad_root(f, p, 1);
455 : }
456 :
457 70 : a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
458 70 : if (lg(a) < 3) pari_err_PRIME("rootmod",p);
459 70 : a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
460 70 : a = FpX_gcd(f,a, p);
461 70 : da = degpol(a);
462 70 : if (!da) return NULL;
463 70 : a = FpX_normalize(a,p);
464 :
465 70 : q = shifti(p,-1);
466 70 : pol0 = icopy(gen_1); /* constant term, will vary in place */
467 70 : pol = deg1pol_shallow(gen_1, pol0, varn(f));
468 70 : for (pol0[2]=1; ; pol0[2]++)
469 : {
470 224 : if (pol0[2] == 1000 && !BPSW_psp(p)) pari_err_PRIME("FpX_oneroot",p);
471 224 : switch(da)
472 : {
473 42 : case 1: return subii(p, gel(a,2));
474 28 : case 2: return FpX_quad_root(a, p, 0);
475 154 : default: {
476 154 : GEN b = FpXQ_pow(pol,q, a,p);
477 : long db;
478 154 : if (degpol(b) <= 0) continue;
479 147 : b = FpX_gcd(a,FpX_Fp_sub_shallow(b,gen_1,p), p);
480 147 : db = degpol(b); if (!db) continue;
481 147 : b = FpX_normalize(b, p);
482 147 : if (db <= (da >> 1)) {
483 105 : a = b;
484 105 : da = db;
485 : } else {
486 42 : a = FpX_div(a,b, p);
487 42 : da -= db;
488 : }
489 : }
490 : }
491 : }
492 : }
493 :
494 : ulong
495 3100272 : Flx_oneroot(GEN f, ulong p)
496 : {
497 3100272 : pari_sp av = avma;
498 : ulong r;
499 3100272 : switch(lg(f))
500 : {
501 0 : case 2: return 0;
502 0 : case 3: return p;
503 : }
504 3100272 : if (p == 2) return Flx_oneroot_mod_2(f);
505 3100272 : r = Flx_oneroot_i(Flx_normalize(f, p), p, 0);
506 3100272 : return gc_ulong(av,r);
507 : }
508 :
509 : ulong
510 32277 : Flx_oneroot_split(GEN f, ulong p)
511 : {
512 32277 : pari_sp av = avma;
513 : ulong r;
514 32277 : switch(lg(f))
515 : {
516 0 : case 2: return 0;
517 0 : case 3: return p;
518 : }
519 32277 : if (p == 2) return Flx_oneroot_mod_2(f);
520 32277 : r = Flx_oneroot_i(Flx_normalize(f, p), p, 1);
521 32285 : return gc_ulong(av,r);
522 : }
523 :
524 : /* assume that p is prime */
525 : GEN
526 12656 : FpX_oneroot(GEN f, GEN pp) {
527 12656 : pari_sp av = avma;
528 12656 : ZX_rootmod_init(&f, pp);
529 12656 : switch(lg(f))
530 : {
531 0 : case 2: set_avma(av); return gen_0;
532 0 : case 3: return gc_NULL(av);
533 : }
534 12656 : if (typ(f) == t_VECSMALL)
535 : {
536 7809 : ulong r, p = pp[2];
537 7809 : if (p == 2)
538 70 : r = Flx_oneroot_mod_2(f);
539 : else
540 7739 : r = Flx_oneroot_i(f, p, 0);
541 7809 : set_avma(av);
542 7809 : return (r == p)? NULL: utoi(r);
543 : }
544 4847 : f = FpX_oneroot_i(f, pp);
545 4847 : if (!f) return gc_NULL(av);
546 4847 : return gerepileuptoint(av, f);
547 : }
548 :
549 : /* returns a root of unity in F_p that is suitable for finding a factor */
550 : /* of degree deg_factor of a polynomial of degree deg; the order is */
551 : /* returned in n */
552 : /* A good choice seems to be n close to deg/deg_factor; we choose n */
553 : /* twice as big and decrement until it divides p-1. */
554 : static GEN
555 651 : good_root_of_unity(GEN p, long deg, long deg_factor, long *pt_n)
556 : {
557 651 : pari_sp ltop = avma;
558 : GEN pm, factn, power, base, zeta;
559 : long n;
560 :
561 651 : pm = subis (p, 1ul);
562 917 : for (n = deg / 2 / deg_factor + 1; !dvdiu (pm, n); n--);
563 651 : factn = Z_factor(stoi(n));
564 651 : power = diviuexact (pm, n);
565 651 : base = gen_1;
566 : do {
567 938 : base = addis (base, 1l);
568 938 : zeta = Fp_pow (base, power, p);
569 : }
570 938 : while (!equaliu (Fp_order (zeta, factn, p), n));
571 651 : *pt_n = n;
572 651 : return gerepileuptoint (ltop, zeta);
573 : }
574 :
575 : GEN
576 924 : FpX_oneroot_split(GEN fact, GEN p)
577 : {
578 924 : pari_sp av = avma;
579 : long n, deg_f, i, dmin;
580 : GEN prim, expo, minfactor, xplusa, zeta, xpow;
581 924 : fact = FpX_normalize(fact, p);
582 924 : deg_f = degpol(fact);
583 924 : if (deg_f<=2) return FpX_oneroot(fact, p);
584 217 : minfactor = fact; /* factor of minimal degree found so far */
585 217 : dmin = degpol(minfactor);
586 217 : xplusa = pol_x(varn(fact));
587 868 : while (dmin != 1)
588 : {
589 : /* split minfactor by computing its gcd with (X+a)^exp-zeta, where */
590 : /* zeta varies over the roots of unity in F_p */
591 651 : fact = minfactor; deg_f = dmin;
592 651 : zeta = gen_1;
593 651 : prim = good_root_of_unity(p, deg_f, 1, &n);
594 651 : expo = diviuexact(subiu(p, 1), n);
595 : /* update X+a, avoid a=0 */
596 651 : gel (xplusa, 2) = addis (gel (xplusa, 2), 1);
597 651 : xpow = FpXQ_pow (xplusa, expo, fact, p);
598 1491 : for (i = 0; i < n; i++)
599 : {
600 1099 : GEN tmp = FpX_gcd(FpX_Fp_sub(xpow, zeta, p), fact, p);
601 1099 : long dtmp = degpol(tmp);
602 1099 : if (dtmp > 0 && dtmp < deg_f)
603 : {
604 455 : fact = FpX_div(fact, tmp, p); deg_f = degpol(fact);
605 455 : if (dtmp < dmin)
606 : {
607 448 : minfactor = FpX_normalize (tmp, p);
608 448 : dmin = dtmp;
609 448 : if (dmin == 1 || dmin <= (2 * deg_f) / n - 1)
610 : /* stop early to avoid too many gcds */
611 : break;
612 : }
613 : }
614 840 : zeta = Fp_mul (zeta, prim, p);
615 : }
616 : }
617 217 : return gerepileuptoint(av, Fp_neg(gel(minfactor,2), p));
618 : }
619 :
620 : /*******************************************************************/
621 : /* */
622 : /* FACTORISATION MODULO p */
623 : /* */
624 : /*******************************************************************/
625 :
626 : /* F / E a vector of vectors of factors / exponents of virtual length l
627 : * (their real lg may be larger). Set their lg to j, concat and return [F,E] */
628 : static GEN
629 583833 : FE_concat(GEN F, GEN E, long l)
630 : {
631 583833 : setlg(E,l); E = shallowconcat1(E);
632 583832 : setlg(F,l); F = shallowconcat1(F); return mkvec2(F,E);
633 : }
634 :
635 : static GEN
636 14 : ddf_to_ddf2_i(GEN V, long fl)
637 : {
638 : GEN F, D;
639 14 : long i, j, l = lg(V);
640 14 : F = cgetg(l, t_VEC);
641 14 : D = cgetg(l, t_VECSMALL);
642 112 : for (i = j = 1; i < l; i++)
643 : {
644 98 : GEN Vi = gel(V,i);
645 98 : if ((fl==2 && F2x_degree(Vi) == 0)
646 98 : ||(fl==0 && degpol(Vi) == 0)) continue;
647 35 : gel(F,j) = Vi;
648 35 : uel(D,j) = i; j++;
649 : }
650 14 : setlg(F,j);
651 14 : setlg(D,j); return mkvec2(F,D);
652 : }
653 :
654 : GEN
655 7 : ddf_to_ddf2(GEN V)
656 7 : { return ddf_to_ddf2_i(V, 0); }
657 :
658 : static GEN
659 7 : F2x_ddf_to_ddf2(GEN V)
660 7 : { return ddf_to_ddf2_i(V, 2); }
661 :
662 : GEN
663 493010 : vddf_to_simplefact(GEN V, long d)
664 : {
665 : GEN E, F;
666 493010 : long i, j, c, l = lg(V);
667 493010 : F = cgetg(d+1, t_VECSMALL);
668 493010 : E = cgetg(d+1, t_VECSMALL);
669 990507 : for (i = c = 1; i < l; i++)
670 : {
671 497497 : GEN Vi = gel(V,i);
672 497497 : long l = lg(Vi);
673 3580031 : for (j = 1; j < l; j++)
674 : {
675 3082534 : long k, n = degpol(gel(Vi,j)) / j;
676 4434010 : for (k = 1; k <= n; k++) { uel(F,c) = j; uel(E,c) = i; c++; }
677 : }
678 : }
679 493010 : setlg(F,c);
680 493010 : setlg(E,c);
681 493010 : return sort_factor(mkvec2(F,E), (void*)&cmpGuGu, cmp_nodata);
682 : }
683 :
684 : /* product of terms of degree 1 in factorization of f */
685 : GEN
686 147095 : FpX_split_part(GEN f, GEN p)
687 : {
688 147095 : long n = degpol(f);
689 147095 : GEN z, X = pol_x(varn(f));
690 147095 : if (n <= 1) return f;
691 145598 : f = FpX_red(f, p);
692 145598 : z = FpX_sub(FpX_Frobenius(f, p), X, p);
693 145598 : return FpX_gcd(z,f,p);
694 : }
695 :
696 : /* Compute the number of roots in Fp without counting multiplicity
697 : * return -1 for 0 polynomial. lc(f) must be prime to p. */
698 : long
699 99261 : FpX_nbroots(GEN f, GEN p)
700 : {
701 99261 : pari_sp av = avma;
702 99261 : GEN z = FpX_split_part(f, p);
703 99261 : return gc_long(av, degpol(z));
704 : }
705 :
706 : /* 1 < deg(f) <= p */
707 : static int
708 0 : Flx_is_totally_split_i(GEN f, ulong p)
709 : {
710 0 : GEN F = Flx_Frobenius(f, p);
711 0 : return degpol(F)==1 && uel(F,2)==0UL && uel(F,3)==1UL;
712 : }
713 : int
714 0 : Flx_is_totally_split(GEN f, ulong p)
715 : {
716 0 : pari_sp av = avma;
717 0 : ulong n = degpol(f);
718 0 : if (n <= 1) return 1;
719 0 : if (n > p) return 0; /* includes n < 0 */
720 0 : return gc_bool(av, Flx_is_totally_split_i(f,p));
721 : }
722 : int
723 0 : FpX_is_totally_split(GEN f, GEN p)
724 : {
725 0 : pari_sp av = avma;
726 0 : ulong n = degpol(f);
727 : int u;
728 0 : if (n <= 1) return 1;
729 0 : if (abscmpui(n, p) > 0) return 0; /* includes n < 0 */
730 0 : if (lgefint(p) != 3)
731 0 : u = gequalX(FpX_Frobenius(FpX_red(f,p), p));
732 : else
733 : {
734 0 : ulong pp = (ulong)p[2];
735 0 : u = Flx_is_totally_split_i(ZX_to_Flx(f,pp), pp);
736 : }
737 0 : return gc_bool(av, u);
738 : }
739 :
740 : long
741 2532564 : Flx_nbroots(GEN f, ulong p)
742 : {
743 2532564 : long n = degpol(f);
744 2532563 : pari_sp av = avma;
745 : GEN z;
746 2532563 : if (n <= 1) return n;
747 2522084 : if (n == 2)
748 : {
749 : ulong D;
750 19012 : if (p==2) return (f[2]==0) + (f[2]!=f[3]);
751 17859 : D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
752 17859 : return 1 + krouu(D,p);
753 : }
754 2503072 : z = Flx_sub(Flx_Frobenius(f, p), polx_Flx(f[1]), p);
755 2503068 : z = Flx_gcd(z, f, p);
756 2503083 : return gc_long(av, degpol(z));
757 : }
758 :
759 : long
760 4403 : FpX_ddf_degree(GEN T, GEN XP, GEN p)
761 : {
762 4403 : pari_sp av = avma;
763 : GEN X, b, g, xq;
764 : long i, j, n, v, B, l, m;
765 : pari_timer ti;
766 : hashtable h;
767 :
768 4403 : n = get_FpX_degree(T); v = get_FpX_var(T);
769 4403 : X = pol_x(v);
770 4403 : if (ZX_equal(X,XP)) return 1;
771 4403 : B = n/2;
772 4403 : l = usqrt(B);
773 4403 : m = (B+l-1)/l;
774 4403 : T = FpX_get_red(T, p);
775 4403 : hash_init_GEN(&h, l+2, ZX_equal, 1);
776 4403 : hash_insert_long(&h, X, 0);
777 4403 : hash_insert_long(&h, XP, 1);
778 4403 : if (DEBUGLEVEL>=7) timer_start(&ti);
779 4403 : b = XP;
780 4403 : xq = FpXQ_powers(b, brent_kung_optpow(n, l-1, 1), T, p);
781 4403 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq baby");
782 10402 : for (i = 3; i <= l+1; i++)
783 : {
784 6790 : b = FpX_FpXQV_eval(b, xq, T, p);
785 6790 : if (gequalX(b)) return gc_long(av,i-1);
786 5999 : hash_insert_long(&h, b, i-1);
787 : }
788 3612 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: baby");
789 3612 : g = b;
790 3612 : xq = FpXQ_powers(g, brent_kung_optpow(n, m, 1), T, p);
791 3612 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq giant");
792 12341 : for(i = 2; i <= m+1; i++)
793 : {
794 10745 : g = FpX_FpXQV_eval(g, xq, T, p);
795 10745 : if (hash_haskey_long(&h, g, &j)) return gc_long(av, l*i-j);
796 : }
797 1596 : return gc_long(av,n);
798 : }
799 :
800 : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
801 : static GEN
802 703 : FpX_ddf_Shoup(GEN T, GEN XP, GEN p)
803 : {
804 : GEN b, g, h, F, f, Tr, xq;
805 : long i, j, n, v, B, l, m;
806 : pari_timer ti;
807 :
808 703 : n = get_FpX_degree(T); v = get_FpX_var(T);
809 703 : if (n == 0) return cgetg(1, t_VEC);
810 703 : if (n == 1) return mkvec(get_FpX_mod(T));
811 657 : B = n/2;
812 657 : l = usqrt(B);
813 657 : m = (B+l-1)/l;
814 657 : T = FpX_get_red(T, p);
815 657 : b = cgetg(l+2, t_VEC);
816 657 : gel(b, 1) = pol_x(v);
817 657 : gel(b, 2) = XP;
818 657 : if (DEBUGLEVEL>=7) timer_start(&ti);
819 657 : xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1), T, p);
820 657 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq baby");
821 849 : for (i = 3; i <= l+1; i++)
822 192 : gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
823 657 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: baby");
824 657 : xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1), T, p);
825 657 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq giant");
826 657 : g = cgetg(m+1, t_VEC);
827 657 : gel(g, 1) = gel(xq, 2);
828 1394 : for(i = 2; i <= m; i++) gel(g, i) = FpX_FpXQV_eval(gel(g, i-1), xq, T, p);
829 657 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: giant");
830 657 : h = cgetg(m+1, t_VEC);
831 2051 : for (j = 1; j <= m; j++)
832 : {
833 1394 : pari_sp av = avma;
834 1394 : GEN gj = gel(g,j), e = FpX_sub(gj, gel(b,1), p);
835 2503 : for (i = 2; i <= l; i++) e = FpXQ_mul(e, FpX_sub(gj, gel(b,i), p), T, p);
836 1394 : gel(h,j) = gerepileupto(av, e);
837 : }
838 657 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: diff");
839 657 : Tr = get_FpX_mod(T);
840 657 : F = cgetg(m+1, t_VEC);
841 2051 : for (j = 1; j <= m; j++)
842 : {
843 1394 : GEN u = FpX_gcd(Tr, gel(h,j), p);
844 1394 : if (degpol(u))
845 : {
846 253 : u = FpX_normalize(u, p);
847 253 : Tr = FpX_div(Tr, u, p);
848 : }
849 1394 : gel(F,j) = u;
850 : }
851 657 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: F");
852 657 : f = const_vec(n, pol_1(v));
853 2051 : for (j = 1; j <= m; j++)
854 : {
855 1394 : GEN e = gel(F, j);
856 1430 : for (i=l-1; i >= 0; i--)
857 : {
858 1430 : GEN u = FpX_gcd(e, FpX_sub(gel(g, j), gel(b, i+1), p), p);
859 1430 : if (degpol(u))
860 : {
861 261 : u = FpX_normalize(u, p);
862 261 : gel(f, l*j-i) = u;
863 261 : e = FpX_div(e, u, p);
864 : }
865 1430 : if (!degpol(e)) break;
866 : }
867 : }
868 657 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: f");
869 657 : if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
870 657 : return f;
871 : }
872 :
873 : static void
874 0 : FpX_edf_simple(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
875 : {
876 0 : long n = degpol(Tp), r = n/d, ct = 0;
877 : GEN T, f, ff, p2;
878 0 : if (r==1) { gel(V, idx) = Tp; return; }
879 0 : p2 = shifti(p,-1);
880 0 : T = FpX_get_red(Tp, p);
881 0 : XP = FpX_rem(XP, T, p);
882 : while (1)
883 0 : {
884 0 : pari_sp btop = avma;
885 : long i;
886 0 : GEN g = random_FpX(n, varn(Tp), p);
887 0 : GEN t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
888 0 : if (signe(t) == 0) continue;
889 0 : for(i=1; i<=10; i++)
890 : {
891 0 : pari_sp btop2 = avma;
892 0 : GEN R = FpXQ_pow(FpX_Fp_add(t, randomi(p), p), p2, T, p);
893 0 : f = FpX_gcd(FpX_Fp_sub(R, gen_1, p), Tp, p);
894 0 : if (degpol(f) > 0 && degpol(f) < n) break;
895 0 : set_avma(btop2);
896 : }
897 0 : if (degpol(f) > 0 && degpol(f) < n) break;
898 0 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_simple",p);
899 0 : set_avma(btop);
900 : }
901 0 : f = FpX_normalize(f, p);
902 0 : ff = FpX_div(Tp, f ,p);
903 0 : FpX_edf_simple(f, XP, d, p, V, idx);
904 0 : FpX_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
905 : }
906 :
907 : static void
908 552 : FpX_edf_rec(GEN T, GEN hp, GEN t, long d, GEN p2, GEN p, GEN V, long idx)
909 : {
910 : pari_sp av;
911 552 : GEN Tp = get_FpX_mod(T);
912 552 : long n = degpol(hp), vT = varn(Tp), ct = 0;
913 : GEN u1, u2, f1, f2, R, h;
914 552 : h = FpX_get_red(hp, p);
915 552 : t = FpX_rem(t, T, p);
916 552 : av = avma;
917 : do
918 : {
919 892 : set_avma(av);
920 892 : R = FpXQ_pow(deg1pol(gen_1, randomi(p), vT), p2, h, p);
921 892 : u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
922 892 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_rec",p);
923 892 : } while (degpol(u1)==0 || degpol(u1)==n);
924 552 : f1 = FpX_gcd(FpX_FpXQ_eval(u1, t, T, p), Tp, p);
925 552 : f1 = FpX_normalize(f1, p);
926 552 : u2 = FpX_div(hp, u1, p);
927 552 : f2 = FpX_div(Tp, f1, p);
928 552 : if (degpol(u1)==1)
929 395 : gel(V, idx) = f1;
930 : else
931 157 : FpX_edf_rec(FpX_get_red(f1, p), u1, t, d, p2, p, V, idx);
932 552 : idx += degpol(f1)/d;
933 552 : if (degpol(u2)==1)
934 409 : gel(V, idx) = f2;
935 : else
936 143 : FpX_edf_rec(FpX_get_red(f2, p), u2, t, d, p2, p, V, idx);
937 552 : }
938 :
939 : /* assume Tp a squarefree product of r > 1 irred. factors of degree d */
940 : static void
941 252 : FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
942 : {
943 252 : long n = degpol(Tp), r = n/d, vT = varn(Tp), ct = 0;
944 : GEN T, h, t;
945 : pari_timer ti;
946 :
947 252 : T = FpX_get_red(Tp, p);
948 252 : XP = FpX_rem(XP, T, p);
949 252 : if (DEBUGLEVEL>=7) timer_start(&ti);
950 : do
951 : {
952 252 : GEN g = random_FpX(n, vT, p);
953 252 : t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
954 252 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
955 252 : h = FpXQ_minpoly(t, T, p);
956 252 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
957 252 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf",p);
958 252 : } while (degpol(h) != r);
959 252 : FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
960 252 : }
961 :
962 : static GEN
963 689 : FpX_factor_Shoup(GEN T, GEN p)
964 : {
965 689 : long i, n, s = 0;
966 : GEN XP, D, V;
967 689 : long e = expi(p);
968 : pari_timer ti;
969 689 : n = get_FpX_degree(T);
970 689 : T = FpX_get_red(T, p);
971 689 : if (DEBUGLEVEL>=6) timer_start(&ti);
972 689 : XP = FpX_Frobenius(T, p);
973 689 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
974 689 : D = FpX_ddf_Shoup(T, XP, p);
975 689 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
976 689 : s = ddf_to_nbfact(D);
977 689 : V = cgetg(s+1, t_COL);
978 5712 : for (i = 1, s = 1; i <= n; i++)
979 : {
980 5023 : GEN Di = gel(D,i);
981 5023 : long ni = degpol(Di), ri = ni/i;
982 5023 : if (ni == 0) continue;
983 699 : Di = FpX_normalize(Di, p);
984 699 : if (ni == i) { gel(V, s++) = Di; continue; }
985 252 : if (ri <= e*expu(e))
986 252 : FpX_edf(Di, XP, i, p, V, s);
987 : else
988 0 : FpX_edf_simple(Di, XP, i, p, V, s);
989 252 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
990 252 : s += ri;
991 : }
992 689 : return V;
993 : }
994 :
995 : long
996 668663 : ddf_to_nbfact(GEN D)
997 : {
998 668663 : long l = lg(D), i, s = 0;
999 4579284 : for(i = 1; i < l; i++) s += degpol(gel(D,i))/i;
1000 668662 : return s;
1001 : }
1002 :
1003 : /* Yun algorithm: Assume p > degpol(T) */
1004 : static GEN
1005 706 : FpX_factor_Yun(GEN T, GEN p)
1006 : {
1007 706 : long n = degpol(T), i = 1;
1008 706 : GEN a, b, c, d = FpX_deriv(T, p);
1009 706 : GEN V = cgetg(n+1,t_VEC);
1010 706 : a = FpX_gcd(T, d, p);
1011 706 : if (degpol(a) == 0) return mkvec(T);
1012 488 : b = FpX_div(T, a, p);
1013 : do
1014 : {
1015 2432 : c = FpX_div(d, a, p);
1016 2432 : d = FpX_sub(c, FpX_deriv(b, p), p);
1017 2432 : a = FpX_normalize(FpX_gcd(b, d, p), p);
1018 2432 : gel(V, i++) = a;
1019 2432 : b = FpX_div(b, a, p);
1020 2432 : } while (degpol(b));
1021 488 : setlg(V, i); return V;
1022 : }
1023 : GEN
1024 777 : FpX_factor_squarefree(GEN T, GEN p)
1025 : {
1026 777 : if (lgefint(p)==3)
1027 : {
1028 777 : ulong pp = (ulong)p[2];
1029 777 : GEN u = Flx_factor_squarefree(ZX_to_Flx(T,pp), pp);
1030 777 : return FlxV_to_ZXV(u);
1031 : }
1032 0 : return FpX_factor_Yun(T, p);
1033 : }
1034 :
1035 : long
1036 168 : FpX_ispower(GEN f, ulong k, GEN p, GEN *pt_r)
1037 : {
1038 168 : pari_sp av = avma;
1039 : GEN lc, F;
1040 168 : long i, l, n = degpol(f), v = varn(f);
1041 168 : if (n % k) return 0;
1042 168 : if (lgefint(p)==3)
1043 : {
1044 126 : ulong pp = p[2];
1045 126 : GEN fp = ZX_to_Flx(f, pp);
1046 126 : if (!Flx_ispower(fp, k, pp, pt_r)) return gc_long(av,0);
1047 105 : if (pt_r) *pt_r = gerepileupto(av, Flx_to_ZX(*pt_r)); else set_avma(av);
1048 105 : return 1;
1049 : }
1050 42 : lc = Fp_sqrtn(leading_coeff(f), stoi(k), p, NULL);
1051 42 : if (!lc) { av = avma; return 0; }
1052 42 : F = FpX_factor_Yun(f, p); l = lg(F)-1;
1053 1491 : for(i=1; i <= l; i++)
1054 1456 : if (i%k && degpol(gel(F,i))) return gc_long(av,0);
1055 35 : if (pt_r)
1056 : {
1057 35 : GEN r = scalarpol(lc, v), s = pol_1(v);
1058 1484 : for (i=l; i>=1; i--)
1059 : {
1060 1449 : if (i%k) continue;
1061 294 : s = FpX_mul(s, gel(F,i), p);
1062 294 : r = FpX_mul(r, s, p);
1063 : }
1064 35 : *pt_r = gerepileupto(av, r);
1065 0 : } else av = avma;
1066 35 : return 1;
1067 : }
1068 :
1069 : static GEN
1070 650 : FpX_factor_Cantor(GEN T, GEN p)
1071 : {
1072 650 : GEN E, F, V = FpX_factor_Yun(T, p);
1073 650 : long i, j, l = lg(V);
1074 650 : F = cgetg(l, t_VEC);
1075 650 : E = cgetg(l, t_VEC);
1076 1816 : for (i=1, j=1; i < l; i++)
1077 1166 : if (degpol(gel(V,i)))
1078 : {
1079 689 : GEN Fj = FpX_factor_Shoup(gel(V,i), p);
1080 689 : gel(F, j) = Fj;
1081 689 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1082 689 : j++;
1083 : }
1084 650 : return sort_factor_pol(FE_concat(F,E,j), cmpii);
1085 : }
1086 :
1087 : static GEN
1088 0 : FpX_ddf_i(GEN T, GEN p)
1089 : {
1090 : GEN XP;
1091 0 : T = FpX_get_red(T, p);
1092 0 : XP = FpX_Frobenius(T, p);
1093 0 : return ddf_to_ddf2(FpX_ddf_Shoup(T, XP, p));
1094 : }
1095 :
1096 : GEN
1097 7 : FpX_ddf(GEN f, GEN p)
1098 : {
1099 7 : pari_sp av = avma;
1100 : GEN F;
1101 7 : switch(ZX_factmod_init(&f, p))
1102 : {
1103 7 : case 0: F = F2x_ddf(f);
1104 7 : F2xV_to_ZXV_inplace(gel(F,1)); break;
1105 0 : case 1: F = Flx_ddf(f,p[2]);
1106 0 : FlxV_to_ZXV_inplace(gel(F,1)); break;
1107 0 : default: F = FpX_ddf_i(f,p); break;
1108 : }
1109 7 : return gerepilecopy(av, F);
1110 : }
1111 :
1112 : static GEN Flx_simplefact_Cantor(GEN T, ulong p);
1113 : static GEN
1114 14 : FpX_simplefact_Cantor(GEN T, GEN p)
1115 : {
1116 : GEN V;
1117 : long i, l;
1118 14 : if (lgefint(p) == 3)
1119 : {
1120 0 : ulong pp = p[2];
1121 0 : return Flx_simplefact_Cantor(ZX_to_Flx(T,pp), pp);
1122 : }
1123 14 : T = FpX_get_red(T, p);
1124 14 : V = FpX_factor_Yun(get_FpX_mod(T), p); l = lg(V);
1125 28 : for (i=1; i < l; i++)
1126 14 : gel(V,i) = FpX_ddf_Shoup(gel(V,i), FpX_Frobenius(gel(V,i), p), p);
1127 14 : return vddf_to_simplefact(V, get_FpX_degree(T));
1128 : }
1129 :
1130 : static int
1131 0 : FpX_isirred_Cantor(GEN Tp, GEN p)
1132 : {
1133 0 : pari_sp av = avma;
1134 : pari_timer ti;
1135 : long n;
1136 0 : GEN T = get_FpX_mod(Tp);
1137 0 : GEN dT = FpX_deriv(T, p);
1138 : GEN XP, D;
1139 0 : if (degpol(FpX_gcd(T, dT, p)) != 0) return gc_bool(av,0);
1140 0 : n = get_FpX_degree(T);
1141 0 : T = FpX_get_red(Tp, p);
1142 0 : if (DEBUGLEVEL>=6) timer_start(&ti);
1143 0 : XP = FpX_Frobenius(T, p);
1144 0 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
1145 0 : D = FpX_ddf_Shoup(T, XP, p);
1146 0 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
1147 0 : return gc_bool(av, degpol(gel(D,n)) == n);
1148 : }
1149 :
1150 : static GEN FpX_factor_deg2(GEN f, GEN p, long d, long flag);
1151 :
1152 : /*Assume that p is large and odd*/
1153 : static GEN
1154 1489 : FpX_factor_i(GEN f, GEN pp, long flag)
1155 : {
1156 1489 : long d = degpol(f);
1157 1489 : if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
1158 664 : switch(flag)
1159 : {
1160 650 : default: return FpX_factor_Cantor(f, pp);
1161 14 : case 1: return FpX_simplefact_Cantor(f, pp);
1162 0 : case 2: return FpX_isirred_Cantor(f, pp)? gen_1: NULL;
1163 : }
1164 : }
1165 :
1166 : long
1167 0 : FpX_nbfact_Frobenius(GEN T, GEN XP, GEN p)
1168 : {
1169 0 : pari_sp av = avma;
1170 0 : long s = ddf_to_nbfact(FpX_ddf_Shoup(T, XP, p));
1171 0 : return gc_long(av,s);
1172 : }
1173 :
1174 : long
1175 0 : FpX_nbfact(GEN T, GEN p)
1176 : {
1177 0 : pari_sp av = avma;
1178 0 : GEN XP = FpX_Frobenius(T, p);
1179 0 : long n = FpX_nbfact_Frobenius(T, XP, p);
1180 0 : return gc_long(av,n);
1181 : }
1182 :
1183 : /* p > 2 */
1184 : static GEN
1185 7 : FpX_is_irred_2(GEN f, GEN p, long d)
1186 : {
1187 7 : switch(d)
1188 : {
1189 0 : case -1:
1190 0 : case 0: return NULL;
1191 0 : case 1: return gen_1;
1192 : }
1193 7 : return FpX_quad_factortype(f, p) == -1? gen_1: NULL;
1194 : }
1195 : /* p > 2 */
1196 : static GEN
1197 14 : FpX_degfact_2(GEN f, GEN p, long d)
1198 : {
1199 14 : switch(d)
1200 : {
1201 0 : case -1:retmkvec2(mkvecsmall(-1),mkvecsmall(1));
1202 0 : case 0: return trivial_fact();
1203 0 : case 1: retmkvec2(mkvecsmall(1), mkvecsmall(1));
1204 : }
1205 14 : switch(FpX_quad_factortype(f, p)) {
1206 7 : case 1: retmkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1207 7 : case -1: retmkvec2(mkvecsmall(2), mkvecsmall(1));
1208 0 : default: retmkvec2(mkvecsmall(1), mkvecsmall(2));
1209 : }
1210 : }
1211 :
1212 : GEN
1213 70 : prime_fact(GEN x) { retmkmat2(mkcolcopy(x), mkcol(gen_1)); }
1214 : GEN
1215 61459 : trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
1216 :
1217 : /* not gerepile safe */
1218 : static GEN
1219 804 : FpX_factor_2(GEN f, GEN p, long d)
1220 : {
1221 : GEN r, s, R, S;
1222 : long v;
1223 : int sgn;
1224 804 : switch(d)
1225 : {
1226 7 : case -1: retmkvec2(mkcol(pol_0(varn(f))), mkvecsmall(1));
1227 30 : case 0: retmkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1228 82 : case 1: retmkvec2(mkcol(f), mkvecsmall(1));
1229 : }
1230 685 : r = FpX_quad_root(f, p, 1);
1231 685 : if (!r) return mkvec2(mkcol(f), mkvecsmall(1));
1232 143 : v = varn(f);
1233 143 : s = FpX_otherroot(f, r, p);
1234 143 : if (signe(r)) r = subii(p, r);
1235 143 : if (signe(s)) s = subii(p, s);
1236 143 : sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
1237 143 : R = deg1pol_shallow(gen_1, r, v);
1238 143 : if (!sgn) return mkvec2(mkcol(R), mkvecsmall(2));
1239 49 : S = deg1pol_shallow(gen_1, s, v);
1240 49 : return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
1241 : }
1242 : static GEN
1243 825 : FpX_factor_deg2(GEN f, GEN p, long d, long flag)
1244 : {
1245 825 : switch(flag) {
1246 7 : case 2: return FpX_is_irred_2(f, p, d);
1247 14 : case 1: return FpX_degfact_2(f, p, d);
1248 804 : default: return FpX_factor_2(f, p, d);
1249 : }
1250 : }
1251 :
1252 : static int
1253 76972 : F2x_quad_factortype(GEN x)
1254 76972 : { return x[2] == 7 ? -1: x[2] == 6 ? 1 :0; }
1255 :
1256 : static GEN
1257 7 : F2x_is_irred_2(GEN f, long d)
1258 7 : { return d == 1 || (d==2 && F2x_quad_factortype(f) == -1)? gen_1: NULL; }
1259 :
1260 : static GEN
1261 10808 : F2x_degfact_2(GEN f, long d)
1262 : {
1263 10808 : if (!d) return trivial_fact();
1264 10808 : if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1265 10598 : switch(F2x_quad_factortype(f)) {
1266 3332 : case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1267 3269 : case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1268 3997 : default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1269 : }
1270 : }
1271 :
1272 : static GEN
1273 122900 : F2x_factor_2(GEN f, long d)
1274 : {
1275 122900 : long v = f[1];
1276 122900 : if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1277 113302 : if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1278 64857 : switch(F2x_quad_factortype(f))
1279 : {
1280 15226 : case -1: return mkvec2(mkcol(f), mkvecsmall(1));
1281 29072 : case 0: return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
1282 20559 : default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
1283 : }
1284 : }
1285 : static GEN
1286 133715 : F2x_factor_deg2(GEN f, long d, long flag)
1287 : {
1288 133715 : switch(flag) {
1289 7 : case 2: return F2x_is_irred_2(f, d);
1290 10808 : case 1: return F2x_degfact_2(f, d);
1291 122900 : default: return F2x_factor_2(f, d);
1292 : }
1293 : }
1294 :
1295 : /* xt = NULL or x^(p-1)/2 mod g */
1296 : static void
1297 5668 : split_squares(struct split_t *S, GEN g, ulong p, GEN xt)
1298 : {
1299 5668 : ulong q = p >> 1;
1300 5668 : GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
1301 5668 : long d = degpol(a);
1302 5668 : if (d < 0)
1303 : {
1304 : ulong i;
1305 203 : split_add_done(S, (GEN)1);
1306 560 : for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr(i,p));
1307 : } else {
1308 5465 : if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
1309 5465 : if (d)
1310 : {
1311 5465 : if (xt) xt = Flx_Fl_add(xt, p-1, p); else xt = Flx_Xnm1(g[1], q, p);
1312 5465 : a = Flx_gcd(a, xt, p);
1313 5465 : if (degpol(a)) split_add(S, Flx_normalize(a, p));
1314 : }
1315 : }
1316 5668 : }
1317 : static void
1318 5668 : split_nonsquares(struct split_t *S, GEN g, ulong p, GEN xt)
1319 : {
1320 5668 : ulong q = p >> 1;
1321 5668 : GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
1322 5668 : long d = degpol(a);
1323 5668 : if (d < 0)
1324 : {
1325 28 : ulong i, z = nonsquare_Fl(p);
1326 28 : split_add_done(S, (GEN)z);
1327 63 : for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_mul(z, Fl_sqr(i,p), p));
1328 : } else {
1329 5640 : if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
1330 5640 : if (d)
1331 : {
1332 5626 : if (xt) xt = Flx_Fl_add(xt, 1, p); else xt = Flx_Xn1(g[1], q, p);
1333 5626 : a = Flx_gcd(a, xt, p);
1334 5626 : if (degpol(a)) split_add(S, Flx_normalize(a, p));
1335 : }
1336 : }
1337 5668 : }
1338 : /* p > 2. f monic Flx, f(0) != 0. Add to split_t structs coprime factors
1339 : * of g = \prod_{f(a) = 0} (X - a). Return 0 when f(x) = 0 for all x in Fp* */
1340 : static int
1341 4927311 : split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p)
1342 : {
1343 4927311 : GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
1344 4927033 : long d = degpol(g);
1345 4926883 : if (d < 0) return 0;
1346 4926841 : if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
1347 4926851 : if (!d) return 1;
1348 4910565 : if ((p >> 4) <= (ulong)d)
1349 : { /* small p; split directly using x^((p-1)/2) +/- 1 */
1350 5262 : GEN xt = ((ulong)d < (p>>1))? Flx_rem(monomial_Flx(1, p>>1, g[1]), g, p)
1351 5668 : : NULL;
1352 5668 : split_squares(S, g, p, xt);
1353 5668 : split_nonsquares(S, g, p, xt);
1354 : } else { /* large p; use x^(p-1) - 1 directly */
1355 4904897 : a = Flxq_powu(polx_Flx(f[1]), p-1, g,p);
1356 4892087 : if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
1357 4892087 : a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
1358 4893182 : g = Flx_gcd(g,a, p);
1359 4897594 : if (degpol(g)) split_add(S, Flx_normalize(g,p));
1360 : }
1361 4907848 : return 1;
1362 : }
1363 :
1364 : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
1365 : static GEN
1366 21933490 : Flx_roots_i(GEN f, ulong p)
1367 : {
1368 : GEN pol, g;
1369 21933490 : long v = Flx_valrem(f, &g);
1370 : ulong q;
1371 : struct split_t S;
1372 :
1373 : /* optimization: test for small degree first */
1374 21940723 : switch(degpol(g))
1375 : {
1376 27861 : case 1: {
1377 27861 : ulong r = p - g[2];
1378 27861 : return v? mkvecsmall2(0, r): mkvecsmall(r);
1379 : }
1380 16975186 : case 2: {
1381 16975186 : ulong r = Flx_quad_root(g, p, 1), s;
1382 16974579 : if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
1383 11583135 : s = Flx_otherroot(g,r, p);
1384 11654639 : if (r < s)
1385 2919353 : return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
1386 8735286 : else if (r > s)
1387 8762386 : return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
1388 : else
1389 138 : return v? mkvecsmall2(0, s): mkvecsmall(s);
1390 : }
1391 : }
1392 4928278 : q = p >> 1;
1393 4928278 : split_init(&S, lg(f)-1);
1394 4927262 : settyp(S.done, t_VECSMALL);
1395 4927262 : if (v) split_add_done(&S, (GEN)0);
1396 4927262 : if (! split_Flx_cut_out_roots(&S, g, p))
1397 42 : return all_roots_mod_p(p, lg(S.done) == 1);
1398 4924060 : pol = polx_Flx(f[1]);
1399 4925506 : for (pol[2]=1; ; pol[2]++)
1400 5253826 : {
1401 10179332 : long j, l = lg(S.todo);
1402 10179332 : if (l == 1) { vecsmall_sort(S.done); return S.done; }
1403 5253589 : if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
1404 10866793 : for (j = 1; j < l; j++)
1405 : {
1406 5612967 : GEN b, c = gel(S.todo,j);
1407 : ulong r, s;
1408 5612967 : switch(degpol(c))
1409 : {
1410 4800197 : case 1:
1411 4800197 : split_moveto_done(&S, j, (GEN)(p - c[2]));
1412 4800535 : j--; l--; break;
1413 373576 : case 2:
1414 373576 : r = Flx_quad_root(c, p, 0);
1415 373584 : if (r == p) pari_err_PRIME("polrootsmod",utoipos(p));
1416 373577 : s = Flx_otherroot(c,r, p);
1417 373581 : split_done(&S, j, (GEN)r, (GEN)s);
1418 373580 : j--; l--; break;
1419 439124 : default:
1420 439124 : b = Flxq_powu(pol,q, c,p); /* pol^(p-1)/2 */
1421 439217 : if (degpol(b) <= 0) continue;
1422 338589 : b = Flx_gcd(c,Flx_Fl_add(b,p-1,p), p);
1423 338600 : if (!degpol(b)) continue;
1424 338466 : b = Flx_normalize(b, p);
1425 338476 : c = Flx_div(c,b, p);
1426 338472 : split_todo(&S, j, b, c);
1427 : }
1428 : }
1429 : }
1430 : }
1431 :
1432 : GEN
1433 21880320 : Flx_roots(GEN f, ulong p)
1434 : {
1435 21880320 : pari_sp av = avma;
1436 21880320 : switch(lg(f))
1437 : {
1438 0 : case 2: pari_err_ROOTS0("Flx_roots");
1439 0 : case 3: set_avma(av); return cgetg(1, t_VECSMALL);
1440 : }
1441 21881610 : if (p == 2) return Flx_root_mod_2(f);
1442 21881610 : return gerepileuptoleaf(av, Flx_roots_i(Flx_normalize(f, p), p));
1443 : }
1444 :
1445 : /* assume x reduced mod p, monic. */
1446 : static int
1447 1134546 : Flx_quad_factortype(GEN x, ulong p)
1448 : {
1449 1134546 : ulong b = x[3], c = x[2];
1450 1134546 : return krouu(Fl_disc_bc(b, c, p), p);
1451 : }
1452 : static GEN
1453 56 : Flx_is_irred_2(GEN f, ulong p, long d)
1454 : {
1455 56 : if (!d) return NULL;
1456 56 : if (d == 1) return gen_1;
1457 56 : return Flx_quad_factortype(f, p) == -1? gen_1: NULL;
1458 : }
1459 : static GEN
1460 1161244 : Flx_degfact_2(GEN f, ulong p, long d)
1461 : {
1462 1161244 : if (!d) return trivial_fact();
1463 1161244 : if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1464 1134490 : switch(Flx_quad_factortype(f, p)) {
1465 541310 : case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1466 579859 : case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1467 13321 : default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1468 : }
1469 : }
1470 : /* p > 2 */
1471 : static GEN
1472 490773 : Flx_factor_2(GEN f, ulong p, long d)
1473 : {
1474 : ulong r, s;
1475 : GEN R,S;
1476 490773 : long v = f[1];
1477 490773 : if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1478 474132 : if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1479 368386 : r = Flx_quad_root(f, p, 1);
1480 368387 : if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
1481 247260 : s = Flx_otherroot(f, r, p);
1482 247270 : r = Fl_neg(r, p);
1483 247270 : s = Fl_neg(s, p);
1484 247270 : if (s < r) lswap(s,r);
1485 247270 : R = mkvecsmall3(v,r,1);
1486 247267 : if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
1487 212477 : S = mkvecsmall3(v,s,1);
1488 212477 : return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
1489 : }
1490 : static GEN
1491 1652074 : Flx_factor_deg2(GEN f, ulong p, long d, long flag)
1492 : {
1493 1652074 : switch(flag) {
1494 56 : case 2: return Flx_is_irred_2(f, p, d);
1495 1161244 : case 1: return Flx_degfact_2(f, p, d);
1496 490774 : default: return Flx_factor_2(f, p, d);
1497 : }
1498 : }
1499 :
1500 : static GEN
1501 10846 : F2x_Berlekamp_ker(GEN u)
1502 : {
1503 10846 : pari_sp ltop=avma;
1504 10846 : long j,N = F2x_degree(u);
1505 : GEN Q;
1506 : pari_timer T;
1507 10845 : timer_start(&T);
1508 10846 : Q = F2x_matFrobenius(u);
1509 251786 : for (j=1; j<=N; j++)
1510 240940 : F2m_flip(Q,j,j);
1511 10846 : if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
1512 10846 : Q = F2m_ker_sp(Q,0);
1513 10846 : if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
1514 10846 : return gerepileupto(ltop,Q);
1515 : }
1516 : #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
1517 : static long
1518 15199 : F2x_split_Berlekamp(GEN *t)
1519 : {
1520 15199 : GEN u = *t, a, b, vker;
1521 15199 : long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
1522 :
1523 15200 : if (du == 1) return 1;
1524 11613 : if (du == 2)
1525 : {
1526 767 : if (F2x_quad_factortype(u) == 1) /* 0 is a root: shouldn't occur */
1527 : {
1528 0 : t[0] = mkvecsmall2(sv, 2);
1529 0 : t[1] = mkvecsmall2(sv, 3);
1530 0 : return 2;
1531 : }
1532 767 : return 1;
1533 : }
1534 :
1535 10846 : vker = F2x_Berlekamp_ker(u);
1536 10846 : lb = lgcols(vker);
1537 10855 : d = lg(vker)-1;
1538 10855 : ir = 0;
1539 : /* t[i] irreducible for i < ir, still to be treated for i < L */
1540 42752 : for (L=1; L<d; )
1541 : {
1542 : GEN pol;
1543 31908 : if (d == 2)
1544 1855 : pol = F2v_to_F2x(gel(vker,2), sv);
1545 : else
1546 : {
1547 30053 : GEN v = zero_zv(lb);
1548 30057 : v[1] = du;
1549 30057 : v[2] = random_Fl(2); /*Assume vker[1]=1*/
1550 112266 : for (i=2; i<=d; i++)
1551 82208 : if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
1552 30058 : pol = F2v_to_F2x(v, sv);
1553 : }
1554 98873 : for (i=ir; i<L && L<d; i++)
1555 : {
1556 66976 : a = t[i]; la = F2x_degree(a);
1557 66967 : if (la == 1) { set_irred(i); }
1558 66785 : else if (la == 2)
1559 : {
1560 743 : if (F2x_quad_factortype(a) == 1) /* 0 is a root: shouldn't occur */
1561 : {
1562 0 : t[i] = mkvecsmall2(sv, 2);
1563 0 : t[L] = mkvecsmall2(sv, 3); L++;
1564 : }
1565 743 : set_irred(i);
1566 : }
1567 : else
1568 : {
1569 66042 : pari_sp av = avma;
1570 : long lb;
1571 66042 : b = F2x_rem(pol, a);
1572 66047 : if (F2x_degree(b) <= 0) { set_avma(av); continue; }
1573 21894 : b = F2x_gcd(a,b); lb = F2x_degree(b);
1574 21894 : if (lb && lb < la)
1575 : {
1576 21894 : t[L] = F2x_div(a,b);
1577 21895 : t[i]= b; L++;
1578 : }
1579 0 : else set_avma(av);
1580 : }
1581 : }
1582 : }
1583 10844 : return d;
1584 : }
1585 : /* assume deg f > 2 */
1586 : static GEN
1587 12383 : F2x_Berlekamp_i(GEN f, long flag)
1588 : {
1589 12383 : long lfact, val, d = F2x_degree(f), j, k, lV;
1590 : GEN y, E, t, V;
1591 :
1592 12384 : val = F2x_valrem(f, &f);
1593 12384 : if (flag == 2 && val) return NULL;
1594 12370 : V = F2x_factor_squarefree(f); lV = lg(V);
1595 12370 : if (flag == 2 && lV > 2) return NULL;
1596 :
1597 : /* to hold factors and exponents */
1598 12300 : t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
1599 12300 : E = cgetg(d+1,t_VECSMALL);
1600 12300 : lfact = 1;
1601 12300 : if (val) {
1602 3532 : if (flag == 1) t[1] = 1; else gel(t,1) = polx_F2x(f[1]);
1603 3532 : E[1] = val; lfact++;
1604 : }
1605 :
1606 59270 : for (k=1; k<lV; k++)
1607 : {
1608 47126 : if (F2x_degree(gel(V, k))==0) continue;
1609 15199 : gel(t,lfact) = gel(V, k);
1610 15199 : d = F2x_split_Berlekamp(&gel(t,lfact));
1611 15198 : if (flag == 2 && d != 1) return NULL;
1612 15044 : if (flag == 1)
1613 4095 : for (j=0; j<d; j++) t[lfact+j] = F2x_degree(gel(t,lfact+j));
1614 51569 : for (j=0; j<d; j++) E[lfact+j] = k;
1615 15044 : lfact += d;
1616 : }
1617 12144 : if (flag == 2) return gen_1; /* irreducible */
1618 12130 : setlg(t, lfact);
1619 12132 : setlg(E, lfact); y = mkvec2(t,E);
1620 2051 : return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
1621 14183 : : sort_factor_pol(y, cmpGuGu);
1622 : }
1623 :
1624 : /* Adapted from Shoup NTL */
1625 : GEN
1626 82433 : F2x_factor_squarefree(GEN f)
1627 : {
1628 : GEN r, t, v, tv;
1629 82433 : long i, q, n = F2x_degree(f);
1630 82431 : GEN u = const_vec(n+1, pol1_F2x(f[1]));
1631 82434 : for(q = 1;;q *= 2)
1632 : {
1633 147773 : r = F2x_gcd(f, F2x_deriv(f));
1634 147766 : if (F2x_degree(r) == 0)
1635 : {
1636 67464 : gel(u, q) = f;
1637 67464 : break;
1638 : }
1639 80304 : t = F2x_div(f, r);
1640 80304 : if (F2x_degree(t) > 0)
1641 : {
1642 : long j;
1643 33463 : for(j = 1;;j++)
1644 : {
1645 74135 : v = F2x_gcd(r, t);
1646 74135 : tv = F2x_div(t, v);
1647 74135 : if (F2x_degree(tv) > 0)
1648 34964 : gel(u, j*q) = tv;
1649 74133 : if (F2x_degree(v) <= 0) break;
1650 40671 : r = F2x_div(r, v);
1651 40672 : t = v;
1652 : }
1653 33463 : if (F2x_degree(r) == 0) break;
1654 : }
1655 65339 : f = F2x_sqrt(r);
1656 : }
1657 585018 : for (i = n; i; i--)
1658 584734 : if (F2x_degree(gel(u,i))) break;
1659 82413 : setlg(u,i+1); return u;
1660 : }
1661 :
1662 : static GEN
1663 86792 : F2x_ddf_simple(GEN T, GEN XP)
1664 : {
1665 86792 : pari_sp av = avma, av2;
1666 : GEN f, z, Tr, X;
1667 86792 : long j, n = F2x_degree(T), v = T[1], B = n/2;
1668 86792 : if (n == 0) return cgetg(1, t_VEC);
1669 86796 : if (n == 1) return mkvec(T);
1670 41281 : z = XP; Tr = T; X = polx_F2x(v);
1671 41285 : f = const_vec(n, pol1_F2x(v));
1672 41285 : av2 = avma;
1673 142256 : for (j = 1; j <= B; j++)
1674 : {
1675 106210 : GEN u = F2x_gcd(Tr, F2x_add(z, X));
1676 106212 : if (F2x_degree(u))
1677 : {
1678 23561 : gel(f, j) = u;
1679 23561 : Tr = F2x_div(Tr, u);
1680 23557 : av2 = avma;
1681 82666 : } else z = gerepileuptoleaf(av2, z);
1682 106269 : if (!F2x_degree(Tr)) break;
1683 100968 : z = F2xq_sqr(z, Tr);
1684 : }
1685 41282 : if (F2x_degree(Tr)) gel(f, F2x_degree(Tr)) = Tr;
1686 41285 : return gerepilecopy(av, f);
1687 : }
1688 :
1689 : GEN
1690 7 : F2x_ddf(GEN T)
1691 : {
1692 : GEN XP;
1693 7 : T = F2x_get_red(T);
1694 7 : XP = F2x_Frobenius(T);
1695 7 : return F2x_ddf_to_ddf2(F2x_ddf_simple(T, XP));
1696 : }
1697 :
1698 : static GEN
1699 7915 : F2xq_frobtrace(GEN a, long d, GEN T)
1700 : {
1701 7915 : pari_sp av = avma;
1702 : long i;
1703 7915 : GEN x = a;
1704 31878 : for(i=1; i<d; i++)
1705 : {
1706 23963 : x = F2x_add(a, F2xq_sqr(x,T));
1707 23963 : if (gc_needed(av, 2))
1708 0 : x = gerepileuptoleaf(av, x);
1709 : }
1710 7915 : return x;
1711 : }
1712 :
1713 : static void
1714 11933 : F2x_edf_simple(GEN Tp, GEN XP, long d, GEN V, long idx)
1715 : {
1716 11933 : long n = F2x_degree(Tp), r = n/d;
1717 : GEN T, f, ff;
1718 11933 : if (r==1) { gel(V, idx) = Tp; return; }
1719 4055 : T = Tp;
1720 4055 : XP = F2x_rem(XP, T);
1721 : while (1)
1722 3860 : {
1723 7915 : pari_sp btop = avma;
1724 : long df;
1725 7915 : GEN g = random_F2x(n, Tp[1]);
1726 7915 : GEN t = F2xq_frobtrace(g, d, T);
1727 7915 : if (lgpol(t) == 0) continue;
1728 5916 : f = F2x_gcd(t, Tp); df = F2x_degree(f);
1729 5916 : if (df > 0 && df < n) break;
1730 1861 : set_avma(btop);
1731 : }
1732 4055 : ff = F2x_div(Tp, f);
1733 4055 : F2x_edf_simple(f, XP, d, V, idx);
1734 4055 : F2x_edf_simple(ff, XP, d, V, idx+F2x_degree(f)/d);
1735 : }
1736 :
1737 : static GEN
1738 86792 : F2x_factor_Shoup(GEN T)
1739 : {
1740 86792 : long i, n, s = 0;
1741 : GEN XP, D, V;
1742 : pari_timer ti;
1743 86792 : n = F2x_degree(T);
1744 86789 : if (DEBUGLEVEL>=6) timer_start(&ti);
1745 86789 : XP = F2x_Frobenius(T);
1746 86785 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1747 86785 : D = F2x_ddf_simple(T, XP);
1748 86793 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1749 370549 : for (i = 1; i <= n; i++)
1750 283765 : s += F2x_degree(gel(D,i))/i;
1751 86784 : V = cgetg(s+1, t_COL);
1752 370533 : for (i = 1, s = 1; i <= n; i++)
1753 : {
1754 283686 : GEN Di = gel(D,i);
1755 283686 : long ni = F2x_degree(Di), ri = ni/i;
1756 283663 : if (ni == 0) continue;
1757 105026 : if (ni == i) { gel(V, s++) = Di; continue; }
1758 3746 : F2x_edf_simple(Di, XP, i, V, s);
1759 3823 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_edf(%ld)",i);
1760 3823 : s += ri;
1761 : }
1762 86847 : return V;
1763 : }
1764 :
1765 : static GEN
1766 70063 : F2x_factor_Cantor(GEN T)
1767 : {
1768 70063 : GEN E, F, V = F2x_factor_squarefree(T);
1769 70063 : long i, j, l = lg(V);
1770 70063 : E = cgetg(l, t_VEC);
1771 70063 : F = cgetg(l, t_VEC);
1772 278268 : for (i=1, j=1; i < l; i++)
1773 208204 : if (F2x_degree(gel(V,i)))
1774 : {
1775 86792 : GEN Fj = F2x_factor_Shoup(gel(V,i));
1776 86789 : gel(F, j) = Fj;
1777 86789 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1778 86792 : j++;
1779 : }
1780 70064 : return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
1781 : }
1782 :
1783 : #if 0
1784 : static GEN
1785 : F2x_simplefact_Shoup(GEN T)
1786 : {
1787 : long i, n, s = 0, j = 1, k;
1788 : GEN XP, D, V;
1789 : pari_timer ti;
1790 : n = F2x_degree(T);
1791 : if (DEBUGLEVEL>=6) timer_start(&ti);
1792 : XP = F2x_Frobenius(T);
1793 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1794 : D = F2x_ddf_simple(T, XP);
1795 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1796 : for (i = 1; i <= n; i++)
1797 : s += F2x_degree(gel(D,i))/i;
1798 : V = cgetg(s+1, t_VECSMALL);
1799 : for (i = 1; i <= n; i++)
1800 : {
1801 : long ni = F2x_degree(gel(D,i)), ri = ni/i;
1802 : if (ni == 0) continue;
1803 : for (k = 1; k <= ri; k++)
1804 : V[j++] = i;
1805 : }
1806 : return V;
1807 : }
1808 : static GEN
1809 : F2x_simplefact_Cantor(GEN T)
1810 : {
1811 : GEN E, F, V = F2x_factor_squarefree(T);
1812 : long i, j, l = lg(V);
1813 : F = cgetg(l, t_VEC);
1814 : E = cgetg(l, t_VEC);
1815 : for (i=1, j=1; i < l; i++)
1816 : if (F2x_degree(gel(V,i)))
1817 : {
1818 : GEN Fj = F2x_simplefact_Shoup(gel(V,i));
1819 : gel(F, j) = Fj;
1820 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1821 : j++;
1822 : }
1823 : return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
1824 : }
1825 : static int
1826 : F2x_isirred_Cantor(GEN T)
1827 : {
1828 : pari_sp av = avma;
1829 : pari_timer ti;
1830 : long n;
1831 : GEN dT = F2x_deriv(T);
1832 : GEN XP, D;
1833 : if (F2x_degree(F2x_gcd(T, dT)) != 0) return gc_bool(av,0);
1834 : n = F2x_degree(T);
1835 : if (DEBUGLEVEL>=6) timer_start(&ti);
1836 : XP = F2x_Frobenius(T);
1837 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1838 : D = F2x_ddf_simple(T, XP);
1839 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1840 : return gc_bool(av, F2x_degree(gel(D,n)) == n);
1841 : }
1842 : #endif
1843 :
1844 : /* driver for Cantor factorization, assume deg f > 2; not competitive for
1845 : * flag != 0, or as deg f increases */
1846 : static GEN
1847 70063 : F2x_Cantor_i(GEN f, long flag)
1848 : {
1849 : switch(flag)
1850 : {
1851 70063 : default: return F2x_factor_Cantor(f);
1852 : #if 0
1853 : case 1: return F2x_simplefact_Cantor(f);
1854 : case 2: return F2x_isirred_Cantor(f)? gen_1: NULL;
1855 : #endif
1856 : }
1857 : }
1858 : static GEN
1859 216162 : F2x_factor_i(GEN f, long flag)
1860 : {
1861 216162 : long d = F2x_degree(f);
1862 216161 : if (d <= 2) return F2x_factor_deg2(f,d,flag);
1863 80143 : return (flag == 0 && d <= 20)? F2x_Cantor_i(f, flag)
1864 162589 : : F2x_Berlekamp_i(f, flag);
1865 : }
1866 :
1867 : GEN
1868 0 : F2x_degfact(GEN f)
1869 : {
1870 0 : pari_sp av = avma;
1871 0 : GEN z = F2x_factor_i(f, 1);
1872 0 : return gerepilecopy(av, z);
1873 : }
1874 :
1875 : int
1876 238 : F2x_is_irred(GEN f) { return !!F2x_factor_i(f, 2); }
1877 :
1878 : /* Adapted from Shoup NTL */
1879 : GEN
1880 1009587 : Flx_factor_squarefree(GEN f, ulong p)
1881 : {
1882 1009587 : long i, q, n = degpol(f);
1883 1009587 : GEN u = const_vec(n+1, pol1_Flx(f[1]));
1884 1009584 : for(q = 1;;q *= p)
1885 54888 : {
1886 1064472 : GEN t, v, tv, r = Flx_gcd(f, Flx_deriv(f, p), p);
1887 1064446 : if (degpol(r) == 0) { gel(u, q) = f; break; }
1888 109866 : t = Flx_div(f, r, p);
1889 109866 : if (degpol(t) > 0)
1890 : {
1891 : long j;
1892 55818 : for(j = 1;;j++)
1893 : {
1894 190231 : v = Flx_gcd(r, t, p);
1895 190232 : tv = Flx_div(t, v, p);
1896 190232 : if (degpol(tv) > 0)
1897 82451 : gel(u, j*q) = Flx_normalize(tv, p);
1898 190231 : if (degpol(v) <= 0) break;
1899 134414 : r = Flx_div(r, v, p);
1900 134413 : t = v;
1901 : }
1902 55817 : if (degpol(r) == 0) break;
1903 : }
1904 54888 : f = Flx_normalize(Flx_deflate(r, p), p);
1905 : }
1906 5974775 : for (i = n; i; i--)
1907 5974765 : if (degpol(gel(u,i))) break;
1908 1009592 : setlg(u,i+1); return u;
1909 : }
1910 :
1911 : long
1912 2821 : Flx_ispower(GEN f, ulong k, ulong p, GEN *pt_r)
1913 : {
1914 2821 : pari_sp av = avma;
1915 : ulong lc;
1916 : GEN F;
1917 2821 : long i, n = degpol(f), v = f[1], l;
1918 2821 : if (n % k) return 0;
1919 2821 : lc = Fl_sqrtn(Flx_lead(f), k, p, NULL);
1920 2821 : if (lc == ULONG_MAX) { av = avma; return 0; }
1921 2821 : F = Flx_factor_squarefree(f, p); l = lg(F)-1;
1922 33299 : for (i = 1; i <= l; i++)
1923 30499 : if (i%k && degpol(gel(F,i))) return gc_long(av,0);
1924 2800 : if (pt_r)
1925 : {
1926 2800 : GEN r = Fl_to_Flx(lc, v), s = pol1_Flx(v);
1927 33278 : for(i = l; i >= 1; i--)
1928 : {
1929 30478 : if (i%k) continue;
1930 10808 : s = Flx_mul(s, gel(F,i), p);
1931 10808 : r = Flx_mul(r, s, p);
1932 : }
1933 2800 : *pt_r = gerepileuptoleaf(av, r);
1934 0 : } else set_avma(av);
1935 2800 : return 1;
1936 : }
1937 :
1938 : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
1939 : static GEN
1940 1451389 : Flx_ddf_Shoup(GEN T, GEN XP, ulong p)
1941 : {
1942 1451389 : pari_sp av = avma;
1943 : GEN b, g, h, F, f, Tr, xq;
1944 : long i, j, n, v, bo, ro;
1945 : long B, l, m;
1946 : pari_timer ti;
1947 1451389 : n = get_Flx_degree(T); v = get_Flx_var(T);
1948 1451398 : if (n == 0) return cgetg(1, t_VEC);
1949 1448185 : if (n == 1) return mkvec(get_Flx_mod(T));
1950 1344906 : B = n/2;
1951 1344906 : l = usqrt(B);
1952 1344907 : m = (B+l-1)/l;
1953 1344907 : T = Flx_get_red(T, p);
1954 1344901 : b = cgetg(l+2, t_VEC);
1955 1344900 : gel(b, 1) = polx_Flx(v);
1956 1344902 : gel(b, 2) = XP;
1957 1344902 : bo = brent_kung_optpow(n, l-1, 1);
1958 1344911 : ro = l<=1 ? 0:(bo-1)/(l-1) + ((n-1)/bo);
1959 1344911 : if (DEBUGLEVEL>=7) timer_start(&ti);
1960 1344911 : if (expu(p) <= ro)
1961 287165 : for (i = 3; i <= l+1; i++)
1962 161793 : gel(b, i) = Flxq_powu(gel(b, i-1), p, T, p);
1963 : else
1964 : {
1965 1219537 : xq = Flxq_powers(gel(b, 2), bo, T, p);
1966 1219532 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq baby");
1967 1453747 : for (i = 3; i <= l+1; i++)
1968 234221 : gel(b, i) = Flx_FlxqV_eval(gel(b, i-1), xq, T, p);
1969 : }
1970 1344898 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: baby");
1971 1344898 : xq = Flxq_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1), T, p);
1972 1344903 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq giant");
1973 1344903 : g = cgetg(m+1, t_VEC);
1974 1344908 : gel(g, 1) = gel(xq, 2);
1975 2600092 : for(i = 2; i <= m; i++)
1976 1255202 : gel(g, i) = Flx_FlxqV_eval(gel(g, i-1), xq, T, p);
1977 1344890 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: giant");
1978 1344890 : h = cgetg(m+1, t_VEC);
1979 3944981 : for (j = 1; j <= m; j++)
1980 : {
1981 2600111 : pari_sp av = avma;
1982 2600111 : GEN gj = gel(g, j);
1983 2600111 : GEN e = Flx_sub(gj, gel(b, 1), p);
1984 4035010 : for (i = 2; i <= l; i++)
1985 1435014 : e = Flxq_mul(e, Flx_sub(gj, gel(b, i), p), T, p);
1986 2599996 : gel(h, j) = gerepileupto(av, e);
1987 : }
1988 1344870 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: diff");
1989 1344870 : Tr = get_Flx_mod(T);
1990 1344911 : F = cgetg(m+1, t_VEC);
1991 3944967 : for (j = 1; j <= m; j++)
1992 : {
1993 2600067 : GEN u = Flx_gcd(Tr, gel(h, j), p);
1994 2600045 : if (degpol(u))
1995 : {
1996 1071837 : u = Flx_normalize(u, p);
1997 1071839 : Tr = Flx_div(Tr, u, p);
1998 : }
1999 2600011 : gel(F, j) = u;
2000 : }
2001 1344900 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: F");
2002 1344900 : f = const_vec(n, pol1_Flx(v));
2003 3944963 : for (j = 1; j <= m; j++)
2004 : {
2005 2600063 : GEN e = gel(F, j);
2006 2881689 : for (i=l-1; i >= 0; i--)
2007 : {
2008 2881724 : GEN u = Flx_gcd(e, Flx_sub(gel(g, j), gel(b, i+1), p), p);
2009 2881657 : if (degpol(u))
2010 : {
2011 1099618 : gel(f, l*j-i) = u;
2012 1099618 : e = Flx_div(e, u, p);
2013 : }
2014 2881644 : if (!degpol(e)) break;
2015 : }
2016 : }
2017 1344900 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: f");
2018 1344900 : if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
2019 1344902 : return gerepilecopy(av, f);
2020 : }
2021 :
2022 : static void
2023 152997 : Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
2024 : {
2025 152997 : long n = degpol(Tp), r = n/d;
2026 : GEN T, f, ff;
2027 : ulong p2;
2028 152998 : if (r==1) { gel(V, idx) = Tp; return; }
2029 66833 : p2 = p>>1;
2030 66833 : T = Flx_get_red(Tp, p);
2031 66833 : XP = Flx_rem(XP, T, p);
2032 : while (1)
2033 7383 : {
2034 74216 : pari_sp btop = avma;
2035 : long i;
2036 74216 : GEN g = random_Flx(n, Tp[1], p);
2037 74216 : GEN t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
2038 74216 : if (lgpol(t) == 0) continue;
2039 160285 : for(i=1; i<=10; i++)
2040 : {
2041 154731 : pari_sp btop2 = avma;
2042 154731 : GEN R = Flxq_powu(Flx_Fl_add(t, random_Fl(p), p), p2, T, p);
2043 154729 : f = Flx_gcd(Flx_Fl_add(R, p-1, p), Tp, p);
2044 154730 : if (degpol(f) > 0 && degpol(f) < n) break;
2045 87897 : set_avma(btop2);
2046 : }
2047 72387 : if (degpol(f) > 0 && degpol(f) < n) break;
2048 5554 : set_avma(btop);
2049 : }
2050 66833 : f = Flx_normalize(f, p);
2051 66833 : ff = Flx_div(Tp, f ,p);
2052 66833 : Flx_edf_simple(f, XP, d, p, V, idx);
2053 66833 : Flx_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
2054 : }
2055 : static void
2056 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx);
2057 :
2058 : static void
2059 514226 : Flx_edf_rec(GEN T, GEN XP, GEN hp, GEN t, long d, ulong p, GEN V, long idx)
2060 : {
2061 : pari_sp av;
2062 514226 : GEN Tp = get_Flx_mod(T);
2063 514229 : long n = degpol(hp), vT = Tp[1];
2064 : GEN u1, u2, f1, f2;
2065 514223 : ulong p2 = p>>1;
2066 : GEN R, h;
2067 514223 : h = Flx_get_red(hp, p);
2068 514210 : t = Flx_rem(t, T, p);
2069 514082 : av = avma;
2070 : do
2071 : {
2072 848209 : set_avma(av);
2073 848207 : R = Flxq_powu(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p);
2074 848071 : u1 = Flx_gcd(Flx_Fl_add(R, p-1, p), hp, p);
2075 848237 : } while (degpol(u1)==0 || degpol(u1)==n);
2076 514164 : f1 = Flx_gcd(Flx_Flxq_eval(u1, t, T, p), Tp, p);
2077 514153 : f1 = Flx_normalize(f1, p);
2078 514156 : u2 = Flx_div(hp, u1, p);
2079 514187 : f2 = Flx_div(Tp, f1, p);
2080 514182 : if (degpol(u1)==1)
2081 : {
2082 373989 : if (degpol(f1)==d)
2083 368377 : gel(V, idx) = f1;
2084 : else
2085 5605 : Flx_edf(f1, XP, d, p, V, idx);
2086 : }
2087 : else
2088 140249 : Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, V, idx);
2089 514223 : idx += degpol(f1)/d;
2090 514205 : if (degpol(u2)==1)
2091 : {
2092 371407 : if (degpol(f2)==d)
2093 366405 : gel(V, idx) = f2;
2094 : else
2095 5002 : Flx_edf(f2, XP, d, p, V, idx);
2096 : }
2097 : else
2098 142831 : Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, V, idx);
2099 514239 : }
2100 :
2101 : static void
2102 231190 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
2103 : {
2104 231190 : long n = degpol(Tp), r = n/d, vT = Tp[1];
2105 : GEN T, h, t;
2106 : pari_timer ti;
2107 231190 : if (r==1) { gel(V, idx) = Tp; return; }
2108 231190 : T = Flx_get_red(Tp, p);
2109 231187 : XP = Flx_rem(XP, T, p);
2110 231169 : if (DEBUGLEVEL>=7) timer_start(&ti);
2111 : do
2112 : {
2113 236774 : GEN g = random_Flx(n, vT, p);
2114 236790 : t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
2115 236807 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
2116 236807 : h = Flxq_minpoly(t, T, p);
2117 236796 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
2118 236796 : } while (degpol(h) <= 1);
2119 231190 : Flx_edf_rec(T, XP, h, t, d, p, V, idx);
2120 : }
2121 :
2122 : static GEN
2123 538842 : Flx_factor_Shoup(GEN T, ulong p)
2124 : {
2125 538842 : long i, n, s = 0;
2126 : GEN XP, D, V;
2127 538842 : long e = expu(p);
2128 : pari_timer ti;
2129 538841 : n = get_Flx_degree(T);
2130 538839 : T = Flx_get_red(T, p);
2131 538830 : if (DEBUGLEVEL>=6) timer_start(&ti);
2132 538830 : XP = Flx_Frobenius(T, p);
2133 538823 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2134 538823 : D = Flx_ddf_Shoup(T, XP, p);
2135 538861 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2136 538861 : s = ddf_to_nbfact(D);
2137 538853 : V = cgetg(s+1, t_COL);
2138 3335370 : for (i = 1, s = 1; i <= n; i++)
2139 : {
2140 2796520 : GEN Di = gel(D,i);
2141 2796520 : long ni = degpol(Di), ri = ni/i;
2142 2796533 : if (ni == 0) continue;
2143 676146 : Di = Flx_normalize(Di, p);
2144 676172 : if (ni == i) { gel(V, s++) = Di; continue; }
2145 239911 : if (ri <= e*expu(e))
2146 220579 : Flx_edf(Di, XP, i, p, V, s);
2147 : else
2148 19332 : Flx_edf_simple(Di, XP, i, p, V, s);
2149 239910 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
2150 239910 : s += ri;
2151 : }
2152 538850 : return V;
2153 : }
2154 :
2155 : static GEN
2156 513120 : Flx_factor_Cantor(GEN T, ulong p)
2157 : {
2158 513120 : GEN E, F, V = Flx_factor_squarefree(get_Flx_mod(T), p);
2159 513111 : long i, j, l = lg(V);
2160 513111 : F = cgetg(l, t_VEC);
2161 513112 : E = cgetg(l, t_VEC);
2162 1256207 : for (i=1, j=1; i < l; i++)
2163 743087 : if (degpol(gel(V,i)))
2164 : {
2165 538843 : GEN Fj = Flx_factor_Shoup(gel(V,i), p);
2166 538850 : gel(F, j) = Fj;
2167 538850 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
2168 538850 : j++;
2169 : }
2170 513120 : return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
2171 : }
2172 :
2173 : GEN
2174 0 : Flx_ddf(GEN T, ulong p)
2175 : {
2176 : GEN XP;
2177 0 : T = Flx_get_red(T, p);
2178 0 : XP = Flx_Frobenius(T, p);
2179 0 : return ddf_to_ddf2(Flx_ddf_Shoup(T, XP, p));
2180 : }
2181 :
2182 : static GEN
2183 492870 : Flx_simplefact_Cantor(GEN T, ulong p)
2184 : {
2185 : GEN V;
2186 : long i, l;
2187 492870 : T = Flx_get_red(T, p);
2188 492870 : V = Flx_factor_squarefree(get_Flx_mod(T), p); l = lg(V);
2189 990178 : for (i=1; i < l; i++)
2190 497308 : gel(V,i) = Flx_ddf_Shoup(gel(V,i), Flx_Frobenius(gel(V,i), p), p);
2191 492870 : return vddf_to_simplefact(V, get_Flx_degree(T));
2192 : }
2193 :
2194 : static int
2195 1078 : Flx_isirred_Cantor(GEN Tp, ulong p)
2196 : {
2197 1078 : pari_sp av = avma;
2198 : pari_timer ti;
2199 : long n;
2200 1078 : GEN T = get_Flx_mod(Tp), dT = Flx_deriv(T, p), XP, D;
2201 1078 : if (degpol(Flx_gcd(T, dT, p)) != 0) return gc_bool(av,0);
2202 791 : n = get_Flx_degree(T);
2203 791 : T = Flx_get_red(Tp, p);
2204 791 : if (DEBUGLEVEL>=6) timer_start(&ti);
2205 791 : XP = Flx_Frobenius(T, p);
2206 791 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2207 791 : D = Flx_ddf_Shoup(T, XP, p);
2208 791 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2209 791 : return gc_bool(av, degpol(gel(D,n)) == n);
2210 : }
2211 :
2212 : /* f monic */
2213 : static GEN
2214 2692987 : Flx_factor_i(GEN f, ulong pp, long flag)
2215 : {
2216 : long d;
2217 2692987 : if (pp==2) { /*We need to handle 2 specially */
2218 33846 : GEN F = F2x_factor_i(Flx_to_F2x(f),flag);
2219 33846 : if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
2220 33846 : return F;
2221 : }
2222 2659141 : d = degpol(f);
2223 2659140 : if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
2224 1007067 : switch(flag)
2225 : {
2226 513119 : default: return Flx_factor_Cantor(f, pp);
2227 492870 : case 1: return Flx_simplefact_Cantor(f, pp);
2228 1078 : case 2: return Flx_isirred_Cantor(f, pp)? gen_1: NULL;
2229 : }
2230 : }
2231 :
2232 : GEN
2233 1665139 : Flx_degfact(GEN f, ulong p)
2234 : {
2235 1665139 : pari_sp av = avma;
2236 1665139 : GEN z = Flx_factor_i(Flx_normalize(f,p),p,1);
2237 1665139 : return gerepilecopy(av, z);
2238 : }
2239 :
2240 : /* T must be squarefree mod p*/
2241 : GEN
2242 285936 : Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
2243 : {
2244 : GEN XP, D;
2245 : pari_timer ti;
2246 285936 : long i, s, n = get_Flx_degree(T);
2247 285936 : GEN V = const_vecsmall(n, 0);
2248 285936 : pari_sp av = avma;
2249 285936 : T = Flx_get_red(T, p);
2250 285936 : if (DEBUGLEVEL>=6) timer_start(&ti);
2251 285936 : XP = Flx_Frobenius(T, p);
2252 285936 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2253 285936 : D = Flx_ddf_Shoup(T, XP, p);
2254 285936 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2255 1678782 : for (i = 1, s = 0; i <= n; i++)
2256 : {
2257 1392846 : V[i] = degpol(gel(D,i))/i;
2258 1392846 : s += V[i];
2259 : }
2260 285936 : *nb = s;
2261 285936 : set_avma(av); return V;
2262 : }
2263 :
2264 : long
2265 128525 : Flx_nbfact_Frobenius(GEN T, GEN XP, ulong p)
2266 : {
2267 128525 : pari_sp av = avma;
2268 128525 : long s = ddf_to_nbfact(Flx_ddf_Shoup(T, XP, p));
2269 128526 : return gc_long(av,s);
2270 : }
2271 :
2272 : /* T must be squarefree mod p*/
2273 : long
2274 128521 : Flx_nbfact(GEN T, ulong p)
2275 : {
2276 128521 : pari_sp av = avma;
2277 128521 : GEN XP = Flx_Frobenius(T, p);
2278 128525 : long n = Flx_nbfact_Frobenius(T, XP, p);
2279 128526 : return gc_long(av,n);
2280 : }
2281 :
2282 : int
2283 1057 : Flx_is_irred(GEN f, ulong p)
2284 : {
2285 1057 : pari_sp av = avma;
2286 1057 : f = Flx_normalize(f,p);
2287 1057 : return gc_bool(av, !!Flx_factor_i(f,p,2));
2288 : }
2289 :
2290 : /* Use this function when you think f is reducible, and that there are lots of
2291 : * factors. If you believe f has few factors, use FpX_nbfact(f,p)==1 instead */
2292 : int
2293 105 : FpX_is_irred(GEN f, GEN p)
2294 : {
2295 105 : pari_sp av = avma;
2296 : int z;
2297 105 : switch(ZX_factmod_init(&f,p))
2298 : {
2299 21 : case 0: z = !!F2x_factor_i(f,2); break;
2300 77 : case 1: z = !!Flx_factor_i(f,p[2],2); break;
2301 7 : default: z = !!FpX_factor_i(f,p,2); break;
2302 : }
2303 105 : return gc_bool(av,z);
2304 : }
2305 : GEN
2306 1862 : FpX_degfact(GEN f, GEN p) {
2307 1862 : pari_sp av = avma;
2308 : GEN F;
2309 1862 : switch(ZX_factmod_init(&f,p))
2310 : {
2311 7 : case 0: F = F2x_factor_i(f,1); break;
2312 1827 : case 1: F = Flx_factor_i(f,p[2],1); break;
2313 28 : default: F = FpX_factor_i(f,p,1); break;
2314 : }
2315 1862 : return gerepilecopy(av, F);
2316 : }
2317 :
2318 : #if 0
2319 : /* set x <-- x + c*y mod p */
2320 : /* x is not required to be normalized.*/
2321 : static void
2322 : Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
2323 : {
2324 : long i, lx, ly;
2325 : ulong *x=(ulong *)gx;
2326 : ulong *y=(ulong *)gy;
2327 : if (!c) return;
2328 : lx = lg(gx);
2329 : ly = lg(gy);
2330 : if (lx<ly) pari_err_BUG("lx<ly in Flx_addmul_inplace");
2331 : if (SMALL_ULONG(p))
2332 : for (i=2; i<ly; i++) x[i] = (x[i] + c*y[i]) % p;
2333 : else
2334 : for (i=2; i<ly; i++) x[i] = Fl_add(x[i], Fl_mul(c,y[i],p),p);
2335 : }
2336 : #endif
2337 :
2338 : GEN
2339 1014192 : FpX_factor(GEN f, GEN p)
2340 : {
2341 1014192 : pari_sp av = avma;
2342 : GEN F;
2343 1014192 : switch(ZX_factmod_init(&f, p))
2344 : {
2345 166832 : case 0: F = F2x_factor_i(f,0);
2346 166832 : F2xV_to_ZXV_inplace(gel(F,1)); break;
2347 845912 : case 1: F = Flx_factor_i(f,p[2],0);
2348 845922 : FlxV_to_ZXV_inplace(gel(F,1)); break;
2349 1453 : default: F = FpX_factor_i(f,p,0); break;
2350 : }
2351 1014199 : return gerepilecopy(av, F);
2352 : }
2353 :
2354 : GEN
2355 178976 : Flx_factor(GEN f, ulong p)
2356 : {
2357 178976 : pari_sp av = avma;
2358 178976 : return gerepilecopy(av, Flx_factor_i(Flx_normalize(f,p),p,0));
2359 : }
2360 : GEN
2361 15217 : F2x_factor(GEN f)
2362 : {
2363 15217 : pari_sp av = avma;
2364 15217 : return gerepilecopy(av, F2x_factor_i(f,0));
2365 : }
|
