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 : /* Check types and replace F by a monic normalized FpX having the same roots
29 : * Don't bother to make constant polynomials monic */
30 : static void
31 350 : factmod_init(GEN *F, GEN p)
32 : {
33 350 : if (typ(p)!=t_INT) pari_err_TYPE("factmod",p);
34 350 : if (signe(p) < 0) pari_err_PRIME("factmod",p);
35 350 : if (typ(*F)!=t_POL) pari_err_TYPE("factmod",*F);
36 350 : if (lgefint(p) == 3)
37 : {
38 224 : ulong pp = p[2];
39 224 : if (pp < 2) pari_err_PRIME("factmod", p);
40 224 : *F = RgX_to_Flx(*F, pp);
41 224 : if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
42 : }
43 : else
44 : {
45 126 : *F = RgX_to_FpX(*F, p);
46 126 : if (lg(*F) > 3) *F = FpX_normalize(*F, p);
47 : }
48 350 : }
49 : /* as above, assume p prime and *F a ZX */
50 : static void
51 938680 : ZX_factmod_init(GEN *F, GEN p)
52 : {
53 938680 : if (lgefint(p) == 3)
54 : {
55 931592 : ulong pp = p[2];
56 931592 : *F = ZX_to_Flx(*F, pp);
57 931592 : if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
58 : }
59 : else
60 : {
61 7088 : *F = FpX_red(*F, p);
62 7088 : if (lg(*F) > 3) *F = FpX_normalize(*F, p);
63 : }
64 938680 : }
65 :
66 : /* return 1,...,p-1 [not_0 = 1] or 0,...,p [not_0 = 0] */
67 : static GEN
68 42 : all_roots_mod_p(ulong p, int not_0)
69 : {
70 : GEN r;
71 : ulong i;
72 42 : if (not_0) {
73 28 : r = cgetg(p, t_VECSMALL);
74 28 : for (i = 1; i < p; i++) r[i] = i;
75 : } else {
76 14 : r = cgetg(p+1, t_VECSMALL);
77 14 : for (i = 0; i < p; i++) r[i+1] = i;
78 : }
79 42 : return r;
80 : }
81 :
82 : /* X^n - 1 */
83 : static GEN
84 13 : Flx_Xnm1(long sv, long n, ulong p)
85 : {
86 13 : GEN t = cgetg(n+3, t_VECSMALL);
87 : long i;
88 13 : t[1] = sv;
89 13 : t[2] = p - 1;
90 13 : for (i = 3; i <= n+1; i++) t[i] = 0;
91 13 : t[i] = 1; return t;
92 : }
93 : /* X^n + 1 */
94 : static GEN
95 6 : Flx_Xn1(long sv, long n, ulong p)
96 : {
97 6 : GEN t = cgetg(n+3, t_VECSMALL);
98 : long i;
99 : (void) p;
100 6 : t[1] = sv;
101 6 : t[2] = 1;
102 6 : for (i = 3; i <= n+1; i++) t[i] = 0;
103 6 : t[i] = 1; return t;
104 : }
105 :
106 : static ulong
107 6 : Fl_nonsquare(ulong p)
108 : {
109 6 : long k = 2;
110 6 : for (;; k++)
111 : {
112 12 : long i = krouu(k, p);
113 12 : if (!i) pari_err_PRIME("Fl_nonsquare",utoipos(p));
114 18 : if (i < 0) return k;
115 6 : }
116 : }
117 :
118 : /* f monic Flx, f(0) != 0. Return a monic squarefree g with the same
119 : * roots as f */
120 : static GEN
121 14 : Flx_cut_out_roots(GEN f, ulong p)
122 : {
123 14 : GEN g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
124 14 : if (g != f && degpol(g) >= 0) {
125 7 : (void)Flx_valrem(g, &g); /* reduction may introduce 0 root */
126 7 : g = Flx_gcd(g, Flx_Xnm1(g[1], p-1, p), p);
127 7 : g = Flx_normalize(g, p);
128 : }
129 14 : return g;
130 : }
131 :
132 : /* by checking f(0..p-1) */
133 : GEN
134 14 : Flx_roots_naive(GEN f, ulong p)
135 : {
136 14 : long d, n = 0;
137 14 : ulong s = 1UL, r;
138 14 : GEN q, y = cgetg(degpol(f) + 1, t_VECSMALL);
139 14 : pari_sp av2, av = avma;
140 :
141 14 : if (Flx_valrem(f, &f)) y[++n] = 0;
142 14 : f = Flx_cut_out_roots(f, p);
143 14 : d = degpol(f);
144 14 : if (d < 0) return all_roots_mod_p(p, n == 0);
145 14 : av2 = avma;
146 308 : while (d > 1) /* d = current degree of f */
147 : {
148 287 : q = Flx_div_by_X_x(f, s, p, &r); /* TODO: FFT-type multi-evaluation */
149 287 : if (r) avma = av2; else { y[++n] = s; d--; f = q; av2 = avma; }
150 287 : if (++s == p) break;
151 : }
152 14 : if (d == 1)
153 : { /* -f[2]/f[3], root of deg 1 polynomial */
154 7 : r = Fl_mul(p - Fl_inv(f[3], p), f[2], p);
155 7 : if (r >= s) y[++n] = r; /* otherwise double root */
156 : }
157 14 : avma = av; fixlg(y, n+1); return y;
158 : }
159 : static GEN
160 8750 : Flx_root_mod_2(GEN f)
161 : {
162 8750 : int z1, z0 = !(f[2] & 1);
163 : long i,n;
164 : GEN y;
165 :
166 8750 : for (i=2, n=1; i < lg(f); i++) n += f[i];
167 8750 : z1 = n & 1;
168 8750 : y = cgetg(z0+z1+1, t_VECSMALL); i = 1;
169 8750 : if (z0) y[i++] = 0;
170 8750 : if (z1) y[i ] = 1;
171 8750 : return y;
172 : }
173 : static ulong
174 35 : Flx_oneroot_mod_2(GEN f)
175 : {
176 : long i,n;
177 35 : if (!(f[2] & 1)) return 0;
178 35 : for (i=2, n=1; i < lg(f); i++) n += f[i];
179 35 : if (n & 1) return 1;
180 21 : return 2;
181 : }
182 :
183 : static GEN FpX_roots_i(GEN f, GEN p);
184 : static GEN Flx_roots_i(GEN f, ulong p);
185 : static GEN FpX_Berlekamp_i(GEN f, GEN pp, long flag);
186 :
187 : static int
188 2883278 : cmpGuGu(GEN a, GEN b) { return (ulong)a < (ulong)b? -1: (a == b? 0: 1); }
189 :
190 : /* Generic driver to computes the roots of f modulo pp, using 'Roots' when
191 : * pp is a small prime.
192 : * if (gpwrap), check types thoroughly and return t_INTMODs, otherwise
193 : * assume that f is an FpX, pp a prime and return t_INTs */
194 : static GEN
195 70512 : rootmod_aux(GEN f, GEN pp, GEN (*Roots)(GEN,ulong), int gpwrap)
196 : {
197 70512 : pari_sp av = avma;
198 : GEN y;
199 70512 : if (gpwrap)
200 84 : factmod_init(&f, pp);
201 : else
202 70428 : ZX_factmod_init(&f, pp);
203 70512 : switch(lg(f))
204 : {
205 14 : case 2: pari_err_ROOTS0("rootmod");
206 49 : case 3: avma = av; return cgetg(1,t_COL);
207 : }
208 70449 : if (typ(f) == t_VECSMALL)
209 : {
210 67652 : ulong p = pp[2];
211 67652 : if (p == 2)
212 8750 : y = Flx_root_mod_2(f);
213 : else
214 : {
215 58902 : if (!odd(p)) pari_err_PRIME("rootmod",utoi(p));
216 58902 : y = Roots(f, p);
217 : }
218 67645 : y = Flc_to_ZC(y);
219 : }
220 : else
221 2797 : y = FpX_roots_i(f, pp);
222 70435 : if (gpwrap) y = FpC_to_mod(y, pp);
223 70435 : return gerepileupto(av, y);
224 : }
225 : /* assume that f is a ZX an pp a prime */
226 : GEN
227 70428 : FpX_roots(GEN f, GEN pp)
228 70428 : { return rootmod_aux(f, pp, Flx_roots_i, 0); }
229 : /* no assumptions on f and pp */
230 : GEN
231 14 : rootmod2(GEN f, GEN pp) { return rootmod_aux(f, pp, &Flx_roots_naive, 1); }
232 : GEN
233 70 : rootmod(GEN f, GEN pp) { return rootmod_aux(f, pp, &Flx_roots_i, 1); }
234 : GEN
235 84 : rootmod0(GEN f, GEN p, long flag)
236 : {
237 84 : switch(flag)
238 : {
239 70 : case 0: return rootmod(f,p);
240 14 : case 1: return rootmod2(f,p);
241 0 : default: pari_err_FLAG("polrootsmod");
242 : }
243 : return NULL; /* LCOV_EXCL_LINE */
244 : }
245 :
246 : /* assume x reduced mod p > 2, monic. */
247 : static int
248 21 : FpX_quad_factortype(GEN x, GEN p)
249 : {
250 21 : GEN b = gel(x,3), c = gel(x,2);
251 21 : GEN D = subii(sqri(b), shifti(c,2));
252 21 : return kronecker(D,p);
253 : }
254 : /* assume x reduced mod p, monic. Return one root, or NULL if irreducible */
255 : GEN
256 6724 : FpX_quad_root(GEN x, GEN p, int unknown)
257 : {
258 6724 : GEN s, D, b = gel(x,3), c = gel(x,2);
259 :
260 6724 : if (absequaliu(p, 2)) {
261 0 : if (!signe(b)) return c;
262 0 : return signe(c)? NULL: gen_1;
263 : }
264 6724 : D = subii(sqri(b), shifti(c,2));
265 6724 : D = remii(D,p);
266 6724 : if (unknown && kronecker(D,p) == -1) return NULL;
267 :
268 6694 : s = Fp_sqrt(D,p);
269 : /* p is not prime, go on and give e.g. maxord a chance to recover */
270 6694 : if (!s) return NULL;
271 6684 : return Fp_halve(Fp_sub(s,b, p), p);
272 : }
273 : static GEN
274 2867 : FpX_otherroot(GEN x, GEN r, GEN p)
275 2867 : { return Fp_neg(Fp_add(gel(x,3), r, p), p); }
276 :
277 : /* disc(x^2+bx+c) = b^2 - 4c */
278 : static ulong
279 19022430 : Fl_disc_bc(ulong b, ulong c, ulong p)
280 19022430 : { return Fl_sub(Fl_sqr(b,p), Fl_double(Fl_double(c,p),p), p); }
281 : /* p > 2 */
282 : static ulong
283 18282974 : Flx_quad_root(GEN x, ulong p, int unknown)
284 : {
285 18282974 : ulong s, b = x[3], c = x[2];
286 18282974 : ulong D = Fl_disc_bc(b, c, p);
287 18278335 : if (unknown && krouu(D,p) == -1) return p;
288 12644877 : s = Fl_sqrt(D,p);
289 12634742 : if (s==~0UL) return p;
290 12634717 : return Fl_halve(Fl_sub(s,b, p), p);
291 : }
292 : static ulong
293 12151958 : Flx_otherroot(GEN x, ulong r, ulong p)
294 12151958 : { return Fl_neg(Fl_add(x[3], r, p), p); }
295 :
296 :
297 : /* 'todo' contains the list of factors to be split.
298 : * 'done' the list of finished factors, no longer touched */
299 : struct split_t { GEN todo, done; };
300 : static void
301 4996045 : split_init(struct split_t *S, long max)
302 : {
303 4996045 : S->todo = vectrunc_init(max);
304 4996238 : S->done = vectrunc_init(max);
305 4996127 : }
306 : #if 0
307 : /* move todo[i] to done */
308 : static void
309 : split_convert(struct split_t *S, long i)
310 : {
311 : long n = lg(S->todo)-1;
312 : vectrunc_append(S->done, gel(S->todo,i));
313 : if (n) gel(S->todo,i) = gel(S->todo, n);
314 : setlg(S->todo, n);
315 : }
316 : #endif
317 : /* append t to todo */
318 : static void
319 5221964 : split_add(struct split_t *S, GEN t) { vectrunc_append(S->todo, t); }
320 : /* delete todo[i], add t to done */
321 : static void
322 5222146 : split_moveto_done(struct split_t *S, long i, GEN t)
323 : {
324 5222146 : long n = lg(S->todo)-1;
325 5222146 : vectrunc_append(S->done, t);
326 5222736 : if (n) gel(S->todo,i) = gel(S->todo, n);
327 5222736 : setlg(S->todo, n);
328 :
329 5222845 : }
330 : /* append t to done */
331 : static void
332 370620 : split_add_done(struct split_t *S, GEN t)
333 370620 : { vectrunc_append(S->done, t); }
334 : /* split todo[i] into a and b */
335 : static void
336 318640 : split_todo(struct split_t *S, long i, GEN a, GEN b)
337 : {
338 318640 : gel(S->todo, i) = a;
339 318640 : split_add(S, b);
340 318647 : }
341 : /* split todo[i] into a and b, moved to done */
342 : static void
343 351643 : split_done(struct split_t *S, long i, GEN a, GEN b)
344 : {
345 351643 : split_moveto_done(S, i, a);
346 351644 : split_add_done(S, b);
347 351644 : }
348 :
349 : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
350 : static GEN
351 2797 : FpX_roots_i(GEN f, GEN p)
352 : {
353 : GEN pol, pol0, a, q;
354 : struct split_t S;
355 :
356 2797 : split_init(&S, lg(f)-1);
357 2797 : settyp(S.done, t_COL);
358 2797 : if (ZX_valrem(f, &f)) split_add_done(&S, gen_0);
359 2797 : switch(degpol(f))
360 : {
361 7 : case 0: return ZC_copy(S.done);
362 14 : case 1: split_add_done(&S, subii(p, gel(f,2))); return ZC_copy(S.done);
363 : case 2: {
364 1603 : GEN s, r = FpX_quad_root(f, p, 1);
365 1603 : if (r) {
366 1603 : split_add_done(&S, r);
367 1603 : s = FpX_otherroot(f,r, p);
368 : /* f not known to be square free yet */
369 1603 : if (!equalii(r, s)) split_add_done(&S, s);
370 : }
371 1603 : return sort(S.done);
372 : }
373 : }
374 :
375 1173 : a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
376 1173 : if (lg(a) < 3) pari_err_PRIME("rootmod",p);
377 1173 : a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
378 1173 : a = FpX_gcd(f,a, p);
379 1173 : if (!degpol(a)) return ZC_copy(S.done);
380 1173 : split_add(&S, FpX_normalize(a,p));
381 :
382 1173 : q = shifti(p,-1);
383 1173 : pol0 = icopy(gen_1); /* constant term, will vary in place */
384 1173 : pol = deg1pol_shallow(gen_1, pol0, varn(f));
385 2444 : for (pol0[2] = 1;; pol0[2]++)
386 : {
387 2444 : long j, l = lg(S.todo);
388 2444 : if (l == 1) return sort(S.done);
389 1278 : if (pol0[2] == 100 && !BPSW_psp(p)) pari_err_PRIME("polrootsmod",p);
390 2647 : for (j = 1; j < l; j++)
391 : {
392 1376 : GEN c = gel(S.todo,j);
393 1376 : switch(degpol(c))
394 : { /* convert linear and quadratics to roots, try to split the rest */
395 : case 1:
396 63 : split_moveto_done(&S, j, subii(p, gel(c,2)));
397 63 : j--; l--; break;
398 : case 2: {
399 1208 : GEN r = FpX_quad_root(c, p, 0), s;
400 1208 : if (!r) pari_err_PRIME("polrootsmod",p);
401 1201 : s = FpX_otherroot(c,r, p);
402 1201 : split_done(&S, j, r, s);
403 1201 : j--; l--; break;
404 : }
405 : default: {
406 : /* b = pol^(p-1)/2 - 1 */
407 105 : GEN b = FpX_Fp_sub_shallow(FpXQ_pow(pol,q, c,p), gen_1, p);
408 : long db;
409 105 : b = FpX_gcd(c,b, p); db = degpol(b);
410 105 : if (db && db < degpol(c))
411 : {
412 98 : b = FpX_normalize(b, p);
413 98 : c = FpX_div(c,b, p);
414 98 : split_todo(&S, j, b, c);
415 : }
416 : }
417 : }
418 : }
419 1271 : }
420 : }
421 :
422 : /* Assume f is normalized */
423 : static ulong
424 254902 : Flx_cubic_root(GEN ff, ulong p)
425 : {
426 254902 : GEN f = Flx_normalize(ff,p);
427 254901 : ulong pi = get_Fl_red(p);
428 254901 : ulong a = f[4], b=f[3], c=f[2], p3 = p%3==1 ? (2*p+1)/3 :(p+1)/3;
429 254901 : ulong t = Fl_mul_pre(a, p3, p, pi), t2 = Fl_sqr_pre(t, p, pi);
430 254900 : ulong A = Fl_sub(b, Fl_triple(t2, p), p);
431 254899 : ulong B = Fl_addmul_pre(c, t, Fl_sub(Fl_double(t2, p), b, p), p, pi);
432 254900 : ulong A3 = Fl_mul_pre(A, p3, p, pi);
433 254901 : ulong A32 = Fl_sqr_pre(A3, p, pi), A33 = Fl_mul_pre(A3, A32, p, pi);
434 254900 : ulong S = Fl_neg(B,p), P = Fl_neg(A3,p);
435 254900 : ulong D = Fl_add(Fl_sqr_pre(S, p, pi), Fl_double(Fl_double(A33, p), p), p);
436 254902 : ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
437 254901 : if (s!=~0UL)
438 : {
439 147237 : ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
440 147238 : if (p%3==2) /* 1 solutions */
441 31578 : vS1 = Fl_powu_pre(S1, (2*p-1)/3, p, pi);
442 : else
443 : {
444 115660 : vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
445 115660 : if (vS1==~0UL) return p; /*0 solutions*/
446 : /*3 solutions*/
447 : }
448 93404 : vS2 = P? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): 0;
449 93404 : return Fl_sub(Fl_add(vS1,vS2, p), t, p);
450 : }
451 : else
452 : {
453 107664 : pari_sp av = avma;
454 107664 : GEN S1 = mkvecsmall2(Fl_halve(S, p), Fl_halve(1UL, p));
455 107663 : GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
456 : ulong Sa;
457 107664 : if (!vS1) return p; /*0 solutions, p%3==2*/
458 104419 : Sa = vS1[1];
459 104419 : if (p%3==1) /*1 solutions*/
460 : {
461 64827 : ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
462 64827 : if (Fa!=P)
463 42715 : Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
464 : }
465 104419 : avma = av;
466 104419 : return Fl_sub(Fl_double(Sa,p),t,p);
467 : }
468 : }
469 :
470 : /* assume p > 2 prime */
471 : static ulong
472 1594314 : Flx_oneroot_i(GEN f, ulong p, long fl)
473 : {
474 : GEN pol, a;
475 : ulong q;
476 : long da;
477 :
478 1594314 : if (Flx_val(f)) return 0;
479 1594143 : switch(degpol(f))
480 : {
481 10790 : case 1: return Fl_neg(f[2], p);
482 402330 : case 2: return Flx_quad_root(f, p, 1);
483 234783 : case 3: if (p>3) return Flx_cubic_root(f, p); /*FALL THROUGH*/
484 : }
485 :
486 946174 : if (!fl)
487 : {
488 916249 : a = Flxq_powu(polx_Flx(f[1]), p - 1, f,p);
489 916249 : if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
490 916249 : a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
491 916249 : a = Flx_gcd(f,a, p);
492 29925 : } else a = f;
493 946174 : da = degpol(a);
494 946188 : if (!da) return p;
495 432657 : a = Flx_normalize(a,p);
496 :
497 432681 : q = p >> 1;
498 432681 : pol = polx_Flx(f[1]);
499 624545 : for(pol[2] = 1;; pol[2]++)
500 : {
501 624545 : if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
502 624475 : switch(da)
503 : {
504 155522 : case 1: return Fl_neg(a[2], p);
505 257147 : case 2: return Flx_quad_root(a, p, 0);
506 20118 : case 3: if (p>3) return Flx_cubic_root(a, p); /*FALL THROUGH*/
507 : default: {
508 191688 : GEN b = Flx_Fl_add(Flxq_powu(pol,q, a,p), p-1, p);
509 : long db;
510 191896 : b = Flx_gcd(a,b, p); db = degpol(b);
511 191894 : if (db && db < da)
512 : {
513 175594 : b = Flx_normalize(b, p);
514 175598 : if (db <= (da >> 1)) {
515 115309 : a = b;
516 115309 : da = db;
517 : } else {
518 60289 : a = Flx_div(a,b, p);
519 60288 : da -= db;
520 : }
521 : }
522 : }
523 : }
524 191897 : }
525 : }
526 :
527 : /* assume p > 2 prime */
528 : static GEN
529 4195 : FpX_oneroot_i(GEN f, GEN p)
530 : {
531 : GEN pol, pol0, a, q;
532 : long da;
533 :
534 4195 : if (ZX_val(f)) return gen_0;
535 4041 : switch(degpol(f))
536 : {
537 315 : case 1: return subii(p, gel(f,2));
538 3656 : case 2: return FpX_quad_root(f, p, 1);
539 : }
540 :
541 70 : a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
542 70 : if (lg(a) < 3) pari_err_PRIME("rootmod",p);
543 70 : a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
544 70 : a = FpX_gcd(f,a, p);
545 70 : da = degpol(a);
546 70 : if (!da) return NULL;
547 70 : a = FpX_normalize(a,p);
548 :
549 70 : q = shifti(p,-1);
550 70 : pol0 = icopy(gen_1); /* constant term, will vary in place */
551 70 : pol = deg1pol_shallow(gen_1, pol0, varn(f));
552 224 : for (pol0[2]=1; ; pol0[2]++)
553 : {
554 224 : if (pol0[2] == 1000 && !BPSW_psp(p)) pari_err_PRIME("FpX_oneroot",p);
555 224 : switch(da)
556 : {
557 42 : case 1: return subii(p, gel(a,2));
558 28 : case 2: return FpX_quad_root(a, p, 0);
559 : default: {
560 154 : GEN b = FpX_Fp_sub_shallow(FpXQ_pow(pol,q, a,p), gen_1, p);
561 : long db;
562 154 : b = FpX_gcd(a,b, p); db = degpol(b);
563 154 : if (db && db < da)
564 : {
565 147 : b = FpX_normalize(b, p);
566 147 : if (db <= (da >> 1)) {
567 105 : a = b;
568 105 : da = db;
569 : } else {
570 42 : a = FpX_div(a,b, p);
571 42 : da -= db;
572 : }
573 : }
574 : }
575 : }
576 154 : }
577 : }
578 :
579 : ulong
580 1556794 : Flx_oneroot(GEN f, ulong p)
581 : {
582 1556794 : pari_sp av = avma;
583 : ulong r;
584 1556794 : switch(lg(f))
585 : {
586 0 : case 2: return 0;
587 1 : case 3: avma = av; return p;
588 : }
589 1556793 : if (p == 2) return Flx_oneroot_mod_2(f);
590 1556793 : r = Flx_oneroot_i(Flx_normalize(f, p), p, 0);
591 1556793 : avma = av; return r;
592 : }
593 :
594 : ulong
595 30350 : Flx_oneroot_split(GEN f, ulong p)
596 : {
597 30350 : pari_sp av = avma;
598 : ulong r;
599 30350 : switch(lg(f))
600 : {
601 0 : case 2: return 0;
602 0 : case 3: avma = av; return p;
603 : }
604 30350 : if (p == 2) return Flx_oneroot_mod_2(f);
605 30350 : r = Flx_oneroot_i(Flx_normalize(f, p), p, 1);
606 30489 : avma = av; return r;
607 : }
608 :
609 : /* assume that p is prime */
610 : GEN
611 11389 : FpX_oneroot(GEN f, GEN pp) {
612 11389 : pari_sp av = avma;
613 11389 : ZX_factmod_init(&f, pp);
614 11389 : switch(lg(f))
615 : {
616 0 : case 2: avma = av; return gen_0;
617 0 : case 3: avma = av; return NULL;
618 : }
619 11389 : if (typ(f) == t_VECSMALL)
620 : {
621 7194 : ulong r, p = pp[2];
622 7194 : if (p == 2)
623 35 : r = Flx_oneroot_mod_2(f);
624 : else
625 7159 : r = Flx_oneroot_i(f, p, 0);
626 7194 : avma = av;
627 7194 : return (r == p)? NULL: utoi(r);
628 : }
629 4195 : f = FpX_oneroot_i(f, pp);
630 4195 : if (!f) { avma = av; return NULL; }
631 4195 : return gerepileuptoint(av, f);
632 : }
633 :
634 : /* returns a root of unity in F_p that is suitable for finding a factor */
635 : /* of degree deg_factor of a polynomial of degree deg; the order is */
636 : /* returned in n */
637 : /* A good choice seems to be n close to deg/deg_factor; we choose n */
638 : /* twice as big and decrement until it divides p-1. */
639 : static GEN
640 238 : good_root_of_unity(GEN p, long deg, long deg_factor, long *pt_n)
641 : {
642 238 : pari_sp ltop = avma;
643 : GEN pm, factn, power, base, zeta;
644 : long n;
645 :
646 238 : pm = subis (p, 1ul);
647 238 : for (n = deg / 2 / deg_factor + 1; !dvdiu (pm, n); n--);
648 238 : factn = Z_factor(stoi(n));
649 238 : power = diviuexact (pm, n);
650 238 : base = gen_1;
651 : do {
652 392 : base = addis (base, 1l);
653 392 : zeta = Fp_pow (base, power, p);
654 : }
655 392 : while (!equaliu (Fp_order (zeta, factn, p), n));
656 238 : *pt_n = n;
657 238 : return gerepileuptoint (ltop, zeta);
658 : }
659 :
660 : GEN
661 665 : FpX_oneroot_split(GEN fact, GEN p)
662 : {
663 665 : pari_sp av = avma;
664 : long n, deg_f, i, dmin;
665 : GEN prim, expo, minfactor, xplusa, zeta, xpow;
666 665 : fact = FpX_normalize(fact, p);
667 665 : deg_f = degpol(fact);
668 665 : if (deg_f<=2) return FpX_oneroot(fact, p);
669 238 : minfactor = fact; /* factor of minimal degree found so far */
670 238 : dmin = degpol(minfactor);
671 238 : prim = good_root_of_unity(p, deg_f, 1, &n);
672 238 : expo = diviuexact(subiu(p, 1), n);
673 238 : xplusa = pol_x(varn(fact));
674 238 : zeta = gen_1;
675 945 : while (dmin != 1)
676 : {
677 : /* split minfactor by computing its gcd with (X+a)^exp-zeta, where */
678 : /* zeta varies over the roots of unity in F_p */
679 469 : fact = minfactor; deg_f = dmin;
680 : /* update X+a, avoid a=0 */
681 469 : gel (xplusa, 2) = addis (gel (xplusa, 2), 1);
682 469 : xpow = FpXQ_pow (xplusa, expo, fact, p);
683 833 : for (i = 0; i < n; i++)
684 : {
685 735 : GEN tmp = FpX_gcd(FpX_Fp_sub(xpow, zeta, p), fact, p);
686 735 : long dtmp = degpol(tmp);
687 735 : if (dtmp > 0 && dtmp < deg_f)
688 : {
689 427 : fact = FpX_div(fact, tmp, p); deg_f = degpol(fact);
690 427 : if (dtmp < dmin)
691 : {
692 427 : minfactor = FpX_normalize (tmp, p);
693 427 : dmin = dtmp;
694 427 : if (dmin == 1 || dmin <= deg_f / (n / 2) + 1)
695 : /* stop early to avoid too many gcds */
696 : break;
697 : }
698 : }
699 364 : zeta = Fp_mul (zeta, prim, p);
700 : }
701 : }
702 238 : return gerepileuptoint(av, Fp_neg(gel(minfactor,2), p));
703 : }
704 :
705 : /*******************************************************************/
706 : /* */
707 : /* FACTORISATION MODULO p */
708 : /* */
709 : /*******************************************************************/
710 :
711 : /* Functions giving information on the factorisation. */
712 :
713 : /* u in Z[X], return kernel of (Frob - Id) over Fp[X] / u */
714 : static GEN
715 60 : FpX_Berlekamp_ker(GEN u, GEN p)
716 : {
717 60 : pari_sp ltop=avma;
718 60 : long j,N = degpol(u);
719 60 : GEN Q = FpX_matFrobenius(u, p);
720 2164 : for (j=1; j<=N; j++)
721 2104 : gcoeff(Q,j,j) = Fp_sub(gcoeff(Q,j,j), gen_1, p);
722 60 : return gerepileupto(ltop, FpM_ker(Q,p));
723 : }
724 :
725 : static GEN
726 19441 : F2x_Berlekamp_ker(GEN u)
727 : {
728 19441 : pari_sp ltop=avma;
729 19441 : long j,N = F2x_degree(u);
730 : GEN Q;
731 : pari_timer T;
732 19441 : timer_start(&T);
733 19441 : Q = F2x_matFrobenius(u);
734 116425 : for (j=1; j<=N; j++)
735 96984 : F2m_flip(Q,j,j);
736 19441 : if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
737 19441 : Q = F2m_ker_sp(Q,0);
738 19441 : if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
739 19441 : return gerepileupto(ltop,Q);
740 : }
741 :
742 : static GEN
743 282216 : Flx_Berlekamp_ker(GEN u, ulong l)
744 : {
745 282216 : pari_sp ltop=avma;
746 282216 : long j,N = degpol(u);
747 : GEN Q;
748 : pari_timer T;
749 282216 : timer_start(&T);
750 282216 : Q = Flx_matFrobenius(u, l);
751 1651341 : for (j=1; j<=N; j++)
752 1369125 : coeff(Q,j,j) = Fl_sub(coeff(Q,j,j),1,l);
753 282216 : if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
754 282216 : Q = Flm_ker_sp(Q,l,0);
755 282216 : if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
756 282216 : return gerepileupto(ltop,Q);
757 : }
758 :
759 : /* product of terms of degree 1 in factorization of f */
760 : GEN
761 142913 : FpX_split_part(GEN f, GEN p)
762 : {
763 142913 : long n = degpol(f);
764 142913 : GEN z, X = pol_x(varn(f));
765 142913 : if (n <= 1) return f;
766 142612 : f = FpX_red(f, p);
767 142612 : z = FpX_sub(FpX_Frobenius(f, p), X, p);
768 142612 : return FpX_gcd(z,f,p);
769 : }
770 :
771 : /* Compute the number of roots in Fp without counting multiplicity
772 : * return -1 for 0 polynomial. lc(f) must be prime to p. */
773 : long
774 99120 : FpX_nbroots(GEN f, GEN p)
775 : {
776 99120 : pari_sp av = avma;
777 99120 : GEN z = FpX_split_part(f, p);
778 99120 : avma = av; return degpol(z);
779 : }
780 :
781 : int
782 0 : FpX_is_totally_split(GEN f, GEN p)
783 : {
784 0 : long n=degpol(f);
785 0 : pari_sp av = avma;
786 0 : if (n <= 1) return 1;
787 0 : if (abscmpui(n, p) > 0) return 0;
788 0 : f = FpX_red(f, p);
789 0 : avma = av; return gequalX(FpX_Frobenius(f, p));
790 : }
791 :
792 : long
793 2374417 : Flx_nbroots(GEN f, ulong p)
794 : {
795 2374417 : long n = degpol(f);
796 2374417 : pari_sp av = avma;
797 : GEN z;
798 2374417 : if (n <= 1) return n;
799 2372310 : if (n == 2)
800 : {
801 : ulong D;
802 11515 : if (p==2) return (f[2]==0) + (f[2]!=f[3]);
803 10577 : D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
804 10577 : return 1 + krouu(D,p);
805 : }
806 2360795 : z = Flx_sub(Flx_Frobenius(f, p), polx_Flx(f[1]), p);
807 2360795 : z = Flx_gcd(z, f, p);
808 2360795 : avma = av; return degpol(z);
809 : }
810 :
811 : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
812 : static GEN
813 4116 : FpX_ddf(GEN T, GEN XP, GEN p)
814 : {
815 4116 : pari_sp av = avma;
816 : GEN b, g, h, F, f, Tr, xq;
817 : long i, j, n, v;
818 : long B, l, m;
819 : pari_timer ti;
820 4116 : n = get_FpX_degree(T); v = get_FpX_var(T);
821 4116 : if (n == 0) return cgetg(1, t_VEC);
822 4116 : if (n == 1) return mkvec(get_FpX_mod(T));
823 4109 : B = n/2;
824 4109 : l = usqrt(B);
825 4109 : m = (B+l-1)/l;
826 4109 : T = FpX_get_red(T, p);
827 4109 : b = cgetg(l+2, t_VEC);
828 4109 : gel(b, 1) = pol_x(v);
829 4109 : gel(b, 2) = XP;
830 4109 : if (DEBUGLEVEL>=7) timer_start(&ti);
831 4109 : xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1), T, p);
832 4109 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: xq baby");
833 10675 : for (i = 3; i <= l+1; i++)
834 6566 : gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
835 4109 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: baby");
836 4109 : xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1), T, p);
837 4109 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: xq giant");
838 4109 : g = cgetg(m+1, t_VEC);
839 4109 : gel(g, 1) = gel(xq, 2);
840 14959 : for(i = 2; i <= m; i++)
841 10850 : gel(g, i) = FpX_FpXQV_eval(gel(g, i-1), xq, T, p);
842 4109 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: giant");
843 4109 : h = cgetg(m+1, t_VEC);
844 19068 : for (j = 1; j <= m; j++)
845 : {
846 14959 : pari_sp av = avma;
847 14959 : GEN gj = gel(g, j);
848 14959 : GEN e = FpX_sub(gj, gel(b, 1), p);
849 42966 : for (i = 2; i <= l; i++)
850 28007 : e = FpXQ_mul(e, FpX_sub(gj, gel(b, i), p), T, p);
851 14959 : gel(h, j) = gerepileupto(av, e);
852 : }
853 4109 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: diff");
854 4109 : Tr = get_FpX_mod(T);
855 4109 : F = cgetg(m+1, t_VEC);
856 19068 : for (j = 1; j <= m; j++)
857 : {
858 14959 : gel(F, j) = FpX_gcd(Tr, gel(h, j), p);
859 14959 : Tr = FpX_div(Tr, gel(F,j), p);
860 : }
861 4109 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: F");
862 4109 : f = const_vec(n, pol_1(v));
863 19068 : for (j = 1; j <= m; j++)
864 : {
865 14959 : GEN e = gel(F, j);
866 16947 : for (i=l-1; i >= 0; i--)
867 : {
868 16947 : GEN u = FpX_gcd(e, FpX_sub(gel(g, j), gel(b, i+1), p), p);
869 16947 : if (degpol(u))
870 : {
871 2716 : gel(f, l*j-i) = u;
872 2716 : e = FpX_div(e, u, p);
873 : }
874 16947 : if (!degpol(e)) break;
875 : }
876 : }
877 4109 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: f");
878 4109 : if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
879 4109 : return gerepilecopy(av, f);
880 : }
881 :
882 : static void
883 0 : FpX_edf_simple(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
884 : {
885 0 : long n = degpol(Tp), r = n/d;
886 : GEN T, f, ff;
887 : GEN p2;
888 0 : if (r==1) { gel(V, idx) = Tp; return; }
889 0 : p2 = shifti(p,-1);
890 0 : T = FpX_get_red(Tp, p);
891 0 : XP = FpX_rem(XP, T, p);
892 : while (1)
893 : {
894 0 : pari_sp btop = avma;
895 : long i;
896 0 : GEN g = random_FpX(n, varn(Tp), p);
897 0 : GEN t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
898 0 : if (signe(t) == 0) continue;
899 0 : for(i=1; i<=10; i++)
900 : {
901 0 : pari_sp btop2 = avma;
902 0 : GEN R = FpXQ_pow(FpX_Fp_add(t, randomi(p), p), p2, T, p);
903 0 : f = FpX_gcd(FpX_Fp_sub(R, gen_1, p), Tp, p);
904 0 : if (degpol(f) > 0 && degpol(f) < n) break;
905 0 : avma = btop2;
906 : }
907 0 : if (degpol(f) > 0 && degpol(f) < n) break;
908 0 : avma = btop;
909 0 : }
910 0 : f = FpX_normalize(f, p);
911 0 : ff = FpX_div(Tp, f ,p);
912 0 : FpX_edf_simple(f, XP, d, p, V, idx);
913 0 : FpX_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
914 : }
915 :
916 : static void
917 210 : FpX_edf_rec(GEN T, GEN hp, GEN t, long d, GEN p2, GEN p, GEN V, long idx)
918 : {
919 : pari_sp av;
920 210 : GEN Tp = get_FpX_mod(T);
921 210 : long n = degpol(hp), vT = varn(Tp);
922 : GEN u1, u2, f1, f2;
923 : GEN R, h;
924 210 : h = FpX_get_red(hp, p);
925 210 : t = FpX_rem(t, T, p);
926 210 : av = avma;
927 : do
928 : {
929 343 : avma = av;
930 343 : R = FpXQ_pow(deg1pol(gen_1, randomi(p), vT), p2, h, p);
931 343 : u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
932 343 : } while (degpol(u1)==0 || degpol(u1)==n);
933 210 : f1 = FpX_gcd(FpX_FpXQ_eval(u1, t, T, p), Tp, p);
934 210 : f1 = FpX_normalize(f1, p);
935 210 : u2 = FpX_div(hp, u1, p);
936 210 : f2 = FpX_div(Tp, f1, p);
937 210 : if (degpol(u1)==1)
938 105 : gel(V, idx) = f1;
939 : else
940 105 : FpX_edf_rec(FpX_get_red(f1, p), u1, t, d, p2, p, V, idx);
941 210 : idx += degpol(f1)/d;
942 210 : if (degpol(u2)==1)
943 119 : gel(V, idx) = f2;
944 : else
945 91 : FpX_edf_rec(FpX_get_red(f2, p), u2, t, d, p2, p, V, idx);
946 210 : }
947 :
948 : static void
949 14 : FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
950 : {
951 14 : long n = degpol(Tp), r = n/d, vT = varn(Tp);
952 : GEN T, h, t;
953 : pari_timer ti;
954 28 : if (r==1) { gel(V, idx) = Tp; return; }
955 14 : T = FpX_get_red(Tp, p);
956 14 : XP = FpX_rem(XP, T, p);
957 14 : if (DEBUGLEVEL>=7) timer_start(&ti);
958 : do
959 : {
960 14 : GEN g = random_FpX(n, vT, p);
961 14 : t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
962 14 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
963 14 : h = FpXQ_minpoly(t, T, p);
964 14 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
965 14 : } while (degpol(h) != r);
966 14 : FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
967 : }
968 :
969 : static GEN
970 14 : FpX_factor_Shoup(GEN T, GEN p)
971 : {
972 14 : long i, n, s = 0;
973 : GEN XP, D, V;
974 14 : long e = expi(p);
975 : pari_timer ti;
976 14 : n = get_FpX_degree(T);
977 14 : T = FpX_get_red(T, p);
978 14 : if (DEBUGLEVEL>=6) timer_start(&ti);
979 14 : XP = FpX_Frobenius(T, p);
980 14 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
981 14 : D = FpX_ddf(T, XP, p);
982 14 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf");
983 301 : for (i = 1; i <= n; i++)
984 287 : s += degpol(gel(D,i))/i;
985 14 : V = cgetg(s+1, t_COL);
986 301 : for (i = 1, s = 1; i <= n; i++)
987 : {
988 287 : GEN Di = gel(D,i);
989 287 : long ni = degpol(Di), ri = ni/i;
990 287 : if (ni == 0) continue;
991 21 : Di = FpX_normalize(Di, p);
992 21 : if (ni == i) { gel(V, s++) = Di; continue; }
993 14 : if (ri <= e*expu(e))
994 14 : FpX_edf(Di, XP, i, p, V, s);
995 : else
996 0 : FpX_edf_simple(Di, XP, i, p, V, s);
997 14 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
998 14 : s += ri;
999 : }
1000 14 : return V;
1001 : }
1002 :
1003 : GEN
1004 352130 : ddf_to_simplefact(GEN D, long n)
1005 : {
1006 : GEN V;
1007 352130 : long i, s = 0, j = 1, k;
1008 2341266 : for (i = 1; i <= n; i++)
1009 1989136 : s += degpol(gel(D,i))/i;
1010 352130 : V = cgetg(s+1, t_VECSMALL);
1011 2341266 : for (i = 1; i <= n; i++)
1012 : {
1013 1989136 : long ni = degpol(gel(D,i)), ri = ni/i;
1014 1989136 : if (ni == 0) continue;
1015 1352907 : for (k = 1; k <= ri; k++)
1016 905845 : uel(V, j++) = i;
1017 : }
1018 352130 : return V;
1019 : }
1020 :
1021 : static GEN
1022 14 : FpX_simplefact_Shoup(GEN T, GEN p)
1023 : {
1024 : GEN XP, D;
1025 14 : long n = get_FpX_degree(T);
1026 : pari_timer ti;
1027 14 : T = FpX_get_red(T, p);
1028 14 : if (DEBUGLEVEL>=6) timer_start(&ti);
1029 14 : XP = FpX_Frobenius(T, p);
1030 14 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
1031 14 : D = FpX_ddf(T, XP, p);
1032 14 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf");
1033 14 : return ddf_to_simplefact(D, n);
1034 : }
1035 :
1036 : /* Yun algorithm: Assume p > degpol(T) */
1037 : static GEN
1038 293 : FpX_factor_Yun(GEN T, GEN p)
1039 : {
1040 293 : long n = degpol(T);
1041 293 : long i = 1;
1042 293 : GEN d = FpX_deriv(T, p);
1043 : GEN a, b, c;
1044 293 : GEN V = cgetg(n+1,t_VEC);
1045 293 : a = FpX_gcd(T, d, p);
1046 293 : if (degpol(a) == 0) return mkvec(T);
1047 159 : b = FpX_div(T, a, p);
1048 : do
1049 : {
1050 346 : c = FpX_div(d, a, p);
1051 346 : d = FpX_sub(c, FpX_deriv(b, p), p);
1052 346 : a = FpX_normalize(FpX_gcd(b, d, p), p);
1053 346 : gel(V, i++) = a;
1054 346 : b = FpX_div(b, a, p);
1055 346 : } while (degpol(b));
1056 159 : setlg(V, i); return V;
1057 : }
1058 : GEN
1059 728 : FpX_factor_squarefree(GEN T, GEN p)
1060 : {
1061 728 : if (abscmpiu(p, degpol(T)) <= 0)
1062 : {
1063 518 : ulong pp = (ulong)p[2];
1064 518 : GEN u = Flx_factor_squarefree(ZX_to_Flx(T,pp), pp);
1065 518 : return FlxV_to_ZXV(u);
1066 : }
1067 210 : return FpX_factor_Yun(T, p);
1068 : }
1069 :
1070 : /* F / E a vector of factors / exponents of virtual length l
1071 : * (their real lg may be larger). Set their lg to j and return [F,E] */
1072 : static GEN
1073 419636 : FE_setlg(GEN F, GEN E, long l)
1074 : {
1075 419636 : setlg(E,l);
1076 419636 : setlg(F,l); return mkvec2(F,E);
1077 : }
1078 : /* F / E a vector of vectors of factors / exponents of virtual length l
1079 : * (their real lg may be larger). Set their lg to j, concat and return [F,E] */
1080 : static GEN
1081 400235 : FE_concat(GEN F, GEN E, long l)
1082 : {
1083 400235 : setlg(E,l); E = shallowconcat1(E);
1084 400235 : setlg(F,l); F = shallowconcat1(F); return mkvec2(F,E);
1085 : }
1086 :
1087 : static GEN
1088 7 : FpX_factor_Cantor(GEN T, GEN p)
1089 : {
1090 7 : GEN E, F, V = FpX_factor_Yun(T, p);
1091 7 : long i, j, l = lg(V);
1092 7 : F = cgetg(l, t_VEC);
1093 7 : E = cgetg(l, t_VEC);
1094 21 : for (i=1, j=1; i < l; i++)
1095 14 : if (degpol(gel(V,i)))
1096 : {
1097 14 : GEN Fj = FpX_factor_Shoup(gel(V,i), p);
1098 14 : gel(F, j) = Fj;
1099 14 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1100 14 : j++;
1101 : }
1102 7 : return sort_factor_pol(FE_concat(F,E,j), cmpii);
1103 : }
1104 :
1105 : static GEN
1106 14 : FpX_simplefact_Cantor(GEN T, GEN p)
1107 : {
1108 : GEN E, F, V;
1109 : long i, j, l;
1110 14 : V = FpX_factor_Yun(get_FpX_mod(T), p);
1111 14 : l = lg(V);
1112 14 : E = cgetg(l, t_VEC);
1113 14 : F = cgetg(l, t_VEC);
1114 28 : for (i=1, j=1; i < l; i++)
1115 14 : if (degpol(gel(V,i)))
1116 : {
1117 14 : GEN Fj = FpX_simplefact_Shoup(gel(V,i), p);
1118 14 : gel(F, j) = Fj;
1119 14 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1120 14 : j++;
1121 : }
1122 14 : return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
1123 : }
1124 :
1125 : static int
1126 0 : FpX_isirred_Cantor(GEN Tp, GEN p)
1127 : {
1128 0 : pari_sp av = avma;
1129 : pari_timer ti;
1130 : long n, d;
1131 0 : GEN T = get_FpX_mod(Tp);
1132 0 : GEN dT = FpX_deriv(T, p);
1133 : GEN XP, D;
1134 0 : if (degpol(FpX_gcd(T, dT, p)) != 0) { avma = av; return 0; }
1135 0 : n = get_FpX_degree(T);
1136 0 : T = FpX_get_red(Tp, p);
1137 0 : if (DEBUGLEVEL>=6) timer_start(&ti);
1138 0 : XP = FpX_Frobenius(T, p);
1139 0 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
1140 0 : D = FpX_ddf(T, XP, p);
1141 0 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf");
1142 0 : d = degpol(gel(D, n));
1143 0 : avma = av; return d==n;
1144 : }
1145 :
1146 : static GEN FpX_factor_deg2(GEN f, GEN p, long d, long flag);
1147 :
1148 : /*Assume that p is large and odd*/
1149 : static GEN
1150 49 : FpX_factcantor_i(GEN f, GEN pp, long flag)
1151 : {
1152 49 : long d = degpol(f);
1153 49 : if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
1154 21 : switch(flag)
1155 : {
1156 7 : default: return FpX_factor_Cantor(f, pp);
1157 14 : case 1: return FpX_simplefact_Cantor(f, pp);
1158 0 : case 2: return FpX_isirred_Cantor(f, pp)? gen_1: NULL;
1159 : }
1160 : }
1161 :
1162 : long
1163 86520 : ddf_to_nbfact(GEN D)
1164 : {
1165 86520 : long l = lg(D), i, s = 0;
1166 793954 : for(i = 1; i < l; i++)
1167 707434 : s += degpol(gel(D,i))/i;
1168 86520 : return s;
1169 : }
1170 :
1171 : long
1172 4088 : FpX_nbfact_Frobenius(GEN T, GEN XP, GEN p)
1173 : {
1174 4088 : pari_sp av = avma;
1175 4088 : long s = ddf_to_nbfact(FpX_ddf(T, XP, p));
1176 4088 : avma = av; return s;
1177 : }
1178 :
1179 : long
1180 28 : FpX_nbfact(GEN T, GEN p)
1181 : {
1182 28 : pari_sp av = avma;
1183 28 : GEN XP = FpX_Frobenius(T, p);
1184 28 : long n = FpX_nbfact_Frobenius(T, XP, p);
1185 28 : avma = av; return n;
1186 : }
1187 :
1188 : /* p > 2 */
1189 : static GEN
1190 7 : FpX_is_irred_2(GEN f, GEN p, long d)
1191 : {
1192 7 : switch(d)
1193 : {
1194 : case -1:
1195 0 : case 0: return NULL;
1196 0 : case 1: return gen_1;
1197 : }
1198 7 : return FpX_quad_factortype(f, p) == -1? gen_1: NULL;
1199 : }
1200 : /* p > 2 */
1201 : static GEN
1202 14 : FpX_degfact_2(GEN f, GEN p, long d)
1203 : {
1204 14 : switch(d)
1205 : {
1206 0 : case -1:retmkvec2(mkvecsmall(-1),mkvecsmall(1));
1207 0 : case 0: return trivial_fact();
1208 0 : case 1: retmkvec2(mkvecsmall(1), mkvecsmall(1));
1209 : }
1210 14 : switch(FpX_quad_factortype(f, p)) {
1211 7 : case 1: retmkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1212 7 : case -1: retmkvec2(mkvecsmall(2), mkvecsmall(1));
1213 0 : default: retmkvec2(mkvecsmall(1), mkvecsmall(2));
1214 : }
1215 : }
1216 :
1217 : GEN
1218 70 : prime_fact(GEN x) { retmkmat2(mkcolcopy(x), mkcol(gen_1)); }
1219 : GEN
1220 945 : trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
1221 :
1222 : /* Mod(0,p) * x, where x is f's main variable */
1223 : static GEN
1224 14 : Mod0pX(GEN f, GEN p)
1225 14 : { return scalarpol(mkintmod(gen_0, p), varn(f)); }
1226 : static GEN
1227 14 : zero_fact_intmod(GEN f, GEN p) { return prime_fact(Mod0pX(f,p)); }
1228 :
1229 : /* not gerepile safe */
1230 : static GEN
1231 90 : FpX_factor_2(GEN f, GEN p, long d)
1232 : {
1233 : GEN r, s, R, S;
1234 : long v;
1235 : int sgn;
1236 90 : switch(d)
1237 : {
1238 0 : case -1: retmkvec2(mkcol(pol_0(varn(f))), mkvecsmall(1));
1239 2 : case 0: retmkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1240 3 : case 1: retmkvec2(mkcol(f), mkvecsmall(1));
1241 : }
1242 85 : r = FpX_quad_root(f, p, 1);
1243 85 : if (!r) return mkvec2(mkcol(f), mkvecsmall(1));
1244 54 : v = varn(f);
1245 54 : s = FpX_otherroot(f, r, p);
1246 54 : if (signe(r)) r = subii(p, r);
1247 54 : if (signe(s)) s = subii(p, s);
1248 54 : sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
1249 54 : R = deg1pol_shallow(gen_1, r, v);
1250 54 : if (!sgn) return mkvec2(mkcol(R), mkvecsmall(2));
1251 47 : S = deg1pol_shallow(gen_1, s, v);
1252 47 : return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
1253 : }
1254 : static GEN
1255 111 : FpX_factor_deg2(GEN f, GEN p, long d, long flag)
1256 : {
1257 111 : switch(flag) {
1258 7 : case 2: return FpX_is_irred_2(f, p, d);
1259 14 : case 1: return FpX_degfact_2(f, p, d);
1260 90 : default: return FpX_factor_2(f, p, d);
1261 : }
1262 : }
1263 :
1264 : static int
1265 58765 : F2x_quad_factortype(GEN x)
1266 58765 : { return x[2] == 7 ? -1: x[2] == 6 ? 1 :0; }
1267 :
1268 : static GEN
1269 7 : F2x_is_irred_2(GEN f, long d)
1270 7 : { return d == 1 || (d==2 && F2x_quad_factortype(f) == -1)? gen_1: NULL; }
1271 :
1272 : static GEN
1273 7252 : F2x_degfact_2(GEN f, long d)
1274 : {
1275 7252 : if (!d) return trivial_fact();
1276 7252 : if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1277 7049 : switch(F2x_quad_factortype(f)) {
1278 2345 : case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1279 2233 : case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1280 2471 : default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1281 : }
1282 : }
1283 :
1284 : static GEN
1285 99680 : F2x_factor_2(GEN f, long d)
1286 : {
1287 99680 : long v = f[1];
1288 99680 : if (d < 0) pari_err(e_ROOTS0,"Flx_factor_2");
1289 99680 : if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1290 92022 : if (d == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1291 45122 : switch(F2x_quad_factortype(f))
1292 : {
1293 9324 : case -1: return mkvec2(mkcol(f), mkvecsmall(1));
1294 21413 : case 0: return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
1295 14385 : default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
1296 : }
1297 : }
1298 : static GEN
1299 106939 : F2x_factor_deg2(GEN f, long d, long flag)
1300 : {
1301 106939 : switch(flag) {
1302 7 : case 2: return F2x_is_irred_2(f, d);
1303 7252 : case 1: return F2x_degfact_2(f, d);
1304 99680 : default: return F2x_factor_2(f, d);
1305 : }
1306 : }
1307 :
1308 : static void
1309 12 : split_squares(struct split_t *S, GEN g, ulong p)
1310 : {
1311 12 : ulong q = p >> 1;
1312 12 : GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
1313 12 : if (degpol(a) < 0)
1314 : {
1315 : ulong i;
1316 6 : split_add_done(S, (GEN)1);
1317 6 : for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr(i,p));
1318 : } else {
1319 6 : (void)Flx_valrem(a, &a);
1320 6 : if (degpol(a))
1321 : {
1322 6 : a = Flx_gcd(a, Flx_Xnm1(g[1], q, p), p);
1323 6 : if (degpol(a)) split_add(S, Flx_normalize(a, p));
1324 : }
1325 : }
1326 12 : }
1327 : static void
1328 12 : split_nonsquares(struct split_t *S, GEN g, ulong p)
1329 : {
1330 12 : ulong q = p >> 1;
1331 12 : GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
1332 12 : if (degpol(a) < 0)
1333 : {
1334 6 : ulong i, z = Fl_nonsquare(p);
1335 6 : split_add_done(S, (GEN)z);
1336 6 : for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_mul(z, Fl_sqr(i,p), p));
1337 : } else {
1338 6 : (void)Flx_valrem(a, &a);
1339 6 : if (degpol(a))
1340 : {
1341 6 : a = Flx_gcd(a, Flx_Xn1(g[1], q, p), p);
1342 6 : if (degpol(a)) split_add(S, Flx_normalize(a, p));
1343 : }
1344 : }
1345 12 : }
1346 : /* p > 2. f monic Flx, f(0) != 0. Add to split_t structs coprime factors
1347 : * of g = \prod_{f(a) = 0} (X - a). Return 0 when f(x) = 0 for all x in Fp* */
1348 : static int
1349 4993079 : split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p)
1350 : {
1351 4993079 : GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
1352 4993177 : long d = degpol(g);
1353 4993384 : if (d < 0) return 0;
1354 4993132 : if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
1355 4993432 : if (!d) return 1;
1356 4977836 : if (p <= 1.4 * (ulong)d) {
1357 : /* small p; split further using x^((p-1)/2) +/- 1.
1358 : * 30% degree drop makes the extra gcd worth it. */
1359 12 : split_squares(S, g, p);
1360 12 : split_nonsquares(S, g, p);
1361 : } else { /* large p; use x^(p-1) - 1 directly */
1362 4977824 : a = Flxq_powu(polx_Flx(f[1]), p-1, g,p);
1363 4976974 : if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
1364 4976974 : a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
1365 4976489 : g = Flx_gcd(g,a, p);
1366 4976749 : if (degpol(g)) split_add(S, Flx_normalize(g,p));
1367 : }
1368 4977318 : return 1;
1369 : }
1370 :
1371 : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
1372 : static GEN
1373 21851716 : Flx_roots_i(GEN f, ulong p)
1374 : {
1375 : GEN pol, g;
1376 21851716 : long v = Flx_valrem(f, &g);
1377 : ulong q;
1378 : struct split_t S;
1379 :
1380 : /* optimization: test for small degree first */
1381 21853043 : switch(degpol(g))
1382 : {
1383 : case 1: {
1384 25529 : ulong r = p - g[2];
1385 25529 : return v? mkvecsmall2(0, r): mkvecsmall(r);
1386 : }
1387 : case 2: {
1388 16833867 : ulong r = Flx_quad_root(g, p, 1), s;
1389 16841539 : if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
1390 11524111 : s = Flx_otherroot(g,r, p);
1391 11534235 : if (r < s)
1392 2874948 : return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
1393 8659287 : else if (r > s)
1394 8659126 : return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
1395 : else
1396 161 : return v? mkvecsmall2(0, s): mkvecsmall(s);
1397 : }
1398 : }
1399 4993355 : q = p >> 1;
1400 4993355 : split_init(&S, lg(f)-1);
1401 4993315 : settyp(S.done, t_VECSMALL);
1402 4993315 : if (v) split_add_done(&S, (GEN)0);
1403 4993315 : if (! split_Flx_cut_out_roots(&S, g, p))
1404 42 : return all_roots_mod_p(p, lg(S.done) == 1);
1405 4992608 : pol = polx_Flx(f[1]);
1406 10294016 : for (pol[2]=1; ; pol[2]++)
1407 : {
1408 10294016 : long j, l = lg(S.todo);
1409 10294016 : if (l == 1) { vecsmall_sort(S.done); return S.done; }
1410 5300855 : if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
1411 10934266 : for (j = 1; j < l; j++)
1412 : {
1413 5633834 : GEN c = gel(S.todo,j);
1414 5633834 : switch(degpol(c))
1415 : {
1416 : case 1:
1417 4870474 : split_moveto_done(&S, j, (GEN)(p - c[2]));
1418 4870857 : j--; l--; break;
1419 : case 2: {
1420 350447 : ulong r = Flx_quad_root(c, p, 0), s;
1421 350449 : if (r == p) pari_err_PRIME("polrootsmod",utoipos(p));
1422 350442 : s = Flx_otherroot(c,r, p);
1423 350442 : split_done(&S, j, (GEN)r, (GEN)s);
1424 350443 : j--; l--; break;
1425 : }
1426 : default: {
1427 412813 : GEN b = Flx_Fl_add(Flxq_powu(pol,q, c,p), p-1, p); /* pol^(p-1)/2 */
1428 : long db;
1429 412803 : b = Flx_gcd(c,b, p); db = degpol(b);
1430 412811 : if (db && db < degpol(c))
1431 : {
1432 318549 : b = Flx_normalize(b, p);
1433 318551 : c = Flx_div(c,b, p);
1434 318545 : split_todo(&S, j, b, c);
1435 : }
1436 : }
1437 : }
1438 : }
1439 5300432 : }
1440 : }
1441 :
1442 : GEN
1443 21792363 : Flx_roots(GEN f, ulong p)
1444 : {
1445 21792363 : pari_sp av = avma;
1446 21792363 : switch(lg(f))
1447 : {
1448 0 : case 2: pari_err_ROOTS0("Flx_roots");
1449 0 : case 3: avma = av; return cgetg(1, t_VECSMALL);
1450 : }
1451 21799465 : if (p == 2) return Flx_root_mod_2(f);
1452 21799465 : return gerepileuptoleaf(av, Flx_roots_i(Flx_normalize(f, p), p));
1453 : }
1454 :
1455 : /* assume x reduced mod p, monic. */
1456 : static int
1457 741216 : Flx_quad_factortype(GEN x, ulong p)
1458 : {
1459 741216 : ulong b = x[3], c = x[2];
1460 741216 : return krouu(Fl_disc_bc(b, c, p), p);
1461 : }
1462 : static GEN
1463 0 : Flx_is_irred_2(GEN f, ulong p, long d)
1464 : {
1465 0 : if (!d) return NULL;
1466 0 : if (d == 1) return gen_1;
1467 0 : return Flx_quad_factortype(f, p) == -1? gen_1: NULL;
1468 : }
1469 : static GEN
1470 758359 : Flx_degfact_2(GEN f, ulong p, long d)
1471 : {
1472 758359 : if (!d) return trivial_fact();
1473 758359 : if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1474 741216 : switch(Flx_quad_factortype(f, p)) {
1475 352366 : case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1476 380093 : case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1477 8757 : default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1478 : }
1479 : }
1480 : /* p > 2 */
1481 : static GEN
1482 413860 : Flx_factor_2(GEN f, ulong p, long d)
1483 : {
1484 : ulong r, s;
1485 : GEN R,S;
1486 413860 : long v = f[1];
1487 413860 : if (d < 0) pari_err(e_ROOTS0,"Flx_factor_2");
1488 413860 : if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1489 399940 : if (d == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1490 301950 : r = Flx_quad_root(f, p, 1);
1491 301950 : if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
1492 199106 : s = Flx_otherroot(f, r, p);
1493 199106 : r = Fl_neg(r, p);
1494 199106 : s = Fl_neg(s, p);
1495 199106 : if (s < r) lswap(s,r);
1496 199106 : R = mkvecsmall3(v,r,1);
1497 199106 : if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
1498 171953 : S = mkvecsmall3(v,s,1);
1499 171953 : return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
1500 : }
1501 : static GEN
1502 1172219 : Flx_factor_deg2(GEN f, ulong p, long d, long flag)
1503 : {
1504 1172219 : switch(flag) {
1505 0 : case 2: return Flx_is_irred_2(f, p, d);
1506 758359 : case 1: return Flx_degfact_2(f, p, d);
1507 413860 : default: return Flx_factor_2(f, p, d);
1508 : }
1509 : }
1510 :
1511 : void
1512 18214 : F2xV_to_FlxV_inplace(GEN v)
1513 : {
1514 : long i;
1515 18214 : for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_Flx(gel(v,i));
1516 18214 : }
1517 : void
1518 718160 : FlxV_to_ZXV_inplace(GEN v)
1519 : {
1520 : long i;
1521 718160 : for(i=1;i<lg(v);i++) gel(v,i)= Flx_to_ZX(gel(v,i));
1522 718160 : }
1523 : void
1524 141106 : F2xV_to_ZXV_inplace(GEN v)
1525 : {
1526 : long i;
1527 141106 : for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_ZX(gel(v,i));
1528 141106 : }
1529 :
1530 : /* Adapted from Shoup NTL */
1531 : GEN
1532 76563 : F2x_factor_squarefree(GEN f)
1533 : {
1534 : GEN r, t, v, tv;
1535 76563 : long i, q, n = F2x_degree(f);
1536 76563 : GEN u = const_vec(n+1, pol1_F2x(f[1]));
1537 121743 : for(q = 1;;q *= 2)
1538 : {
1539 121743 : r = F2x_gcd(f, F2x_deriv(f));
1540 121738 : if (F2x_degree(r) == 0)
1541 : {
1542 63997 : gel(u, q) = f;
1543 63997 : break;
1544 : }
1545 57746 : t = F2x_div(f, r);
1546 57746 : if (F2x_degree(t) > 0)
1547 : {
1548 : long j;
1549 53608 : for(j = 1;;j++)
1550 : {
1551 53608 : v = F2x_gcd(r, t);
1552 53608 : tv = F2x_div(t, v);
1553 53608 : if (F2x_degree(tv) > 0)
1554 20352 : gel(u, j*q) = tv;
1555 53608 : if (F2x_degree(v) <= 0) break;
1556 34654 : r = F2x_div(r, v);
1557 34654 : t = v;
1558 34654 : }
1559 18954 : if (F2x_degree(r) == 0) break;
1560 : }
1561 45179 : f = F2x_sqrt(r);
1562 45179 : }
1563 521225 : for (i = n; i; i--)
1564 515802 : if (F2x_degree(gel(u,i))) break;
1565 76582 : setlg(u,i+1); return u;
1566 : }
1567 :
1568 : static GEN
1569 26232 : F2x_ddf_simple(GEN T, GEN XP)
1570 : {
1571 26232 : pari_sp av = avma, av2;
1572 : GEN f, z, Tr, X;
1573 26232 : long j, n = F2x_degree(T), v = T[1], B = n/2;
1574 26234 : if (n == 0) return cgetg(1, t_VEC);
1575 26234 : if (n == 1) return mkvec(T);
1576 19654 : z = XP; Tr = T; X = polx_F2x(v);
1577 19655 : f = const_vec(n, pol1_F2x(v));
1578 19655 : av2 = avma;
1579 158330 : for (j = 1; j <= B; j++)
1580 : {
1581 148606 : GEN u = F2x_gcd(Tr, F2x_add(z, X));
1582 148603 : if (F2x_degree(u))
1583 : {
1584 42302 : gel(f, j) = u;
1585 42302 : Tr = F2x_div(Tr, u);
1586 42299 : av2 = avma;
1587 106321 : } else z = gerepileuptoleaf(av2, z);
1588 148634 : if (!F2x_degree(Tr)) break;
1589 138679 : z = F2xq_sqr(z, Tr);
1590 : }
1591 19653 : if (F2x_degree(Tr)) gel(f, F2x_degree(Tr)) = Tr;
1592 19655 : return gerepilecopy(av, f);
1593 : }
1594 :
1595 : static GEN
1596 8324 : F2xq_frobtrace(GEN a, long d, GEN T)
1597 : {
1598 8324 : pari_sp av = avma;
1599 : long i;
1600 8324 : GEN x = a;
1601 54543 : for(i=1; i<d; i++)
1602 : {
1603 46218 : x = F2x_add(a, F2xq_sqr(x,T));
1604 46219 : if (gc_needed(av, 2))
1605 0 : x = gerepileuptoleaf(av, x);
1606 : }
1607 8325 : return x;
1608 : }
1609 :
1610 : static void
1611 12812 : F2x_edf_simple(GEN Tp, GEN XP, long d, GEN V, long idx)
1612 : {
1613 12812 : long n = F2x_degree(Tp), r = n/d;
1614 : GEN T, f, ff;
1615 25624 : if (r==1) { gel(V, idx) = Tp; return; }
1616 4392 : T = Tp;
1617 4392 : XP = F2x_rem(XP, T);
1618 : while (1)
1619 : {
1620 8324 : pari_sp btop = avma;
1621 : long df;
1622 8324 : GEN g = random_F2x(n, Tp[1]);
1623 8325 : GEN t = F2xq_frobtrace(g, d, T);
1624 8325 : if (lgpol(t) == 0) continue;
1625 6447 : f = F2x_gcd(t, Tp); df = F2x_degree(f);
1626 6447 : if (df > 0 && df < n) break;
1627 2055 : avma = btop;
1628 3932 : }
1629 4392 : ff = F2x_div(Tp, f);
1630 4392 : F2x_edf_simple(f, XP, d, V, idx);
1631 4392 : F2x_edf_simple(ff, XP, d, V, idx+F2x_degree(f)/d);
1632 : }
1633 :
1634 : static GEN
1635 24288 : F2x_factor_Shoup(GEN T)
1636 : {
1637 24288 : long i, n, s = 0;
1638 : GEN XP, D, V;
1639 : pari_timer ti;
1640 24288 : n = F2x_degree(T);
1641 24288 : if (DEBUGLEVEL>=6) timer_start(&ti);
1642 24288 : XP = F2x_Frobenius(T);
1643 24286 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1644 24286 : D = F2x_ddf_simple(T, XP);
1645 24289 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf");
1646 360882 : for (i = 1; i <= n; i++)
1647 336595 : s += F2x_degree(gel(D,i))/i;
1648 24287 : V = cgetg(s+1, t_COL);
1649 360892 : for (i = 1, s = 1; i <= n; i++)
1650 : {
1651 336603 : GEN Di = gel(D,i);
1652 336603 : long ni = F2x_degree(Di), ri = ni/i;
1653 336591 : if (ni == 0) continue;
1654 55955 : if (ni == i) { gel(V, s++) = Di; continue; }
1655 4028 : F2x_edf_simple(Di, XP, i, V, s);
1656 4028 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_edf(%ld)",i);
1657 4028 : s += ri;
1658 : }
1659 24289 : return V;
1660 : }
1661 :
1662 : static GEN
1663 1778 : F2x_simplefact_Shoup(GEN T)
1664 : {
1665 1778 : long i, n, s = 0, j = 1, k;
1666 : GEN XP, D, V;
1667 : pari_timer ti;
1668 1778 : n = F2x_degree(T);
1669 1778 : if (DEBUGLEVEL>=6) timer_start(&ti);
1670 1778 : XP = F2x_Frobenius(T);
1671 1778 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1672 1778 : D = F2x_ddf_simple(T, XP);
1673 1778 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf");
1674 9065 : for (i = 1; i <= n; i++)
1675 7287 : s += F2x_degree(gel(D,i))/i;
1676 1778 : V = cgetg(s+1, t_VECSMALL);
1677 9065 : for (i = 1; i <= n; i++)
1678 : {
1679 7287 : long ni = F2x_degree(gel(D,i)), ri = ni/i;
1680 7287 : if (ni == 0) continue;
1681 4333 : for (k = 1; k <= ri; k++)
1682 2219 : V[j++] = i;
1683 : }
1684 1778 : return V;
1685 : }
1686 :
1687 : static GEN
1688 17616 : F2x_factor_Cantor(GEN T)
1689 : {
1690 17616 : GEN E, F, V = F2x_factor_squarefree(T);
1691 17616 : long i, j, l = lg(V);
1692 17616 : E = cgetg(l, t_VEC);
1693 17617 : F = cgetg(l, t_VEC);
1694 46398 : for (i=1, j=1; i < l; i++)
1695 28781 : if (F2x_degree(gel(V,i)))
1696 : {
1697 24288 : GEN Fj = F2x_factor_Shoup(gel(V,i));
1698 24288 : gel(F, j) = Fj;
1699 24288 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1700 24288 : j++;
1701 : }
1702 17617 : return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
1703 : }
1704 :
1705 : static GEN
1706 1736 : F2x_simplefact_Cantor(GEN T)
1707 : {
1708 1736 : GEN E, F, V = F2x_factor_squarefree(T);
1709 1736 : long i, j, l = lg(V);
1710 1736 : F = cgetg(l, t_VEC);
1711 1736 : E = cgetg(l, t_VEC);
1712 4690 : for (i=1, j=1; i < l; i++)
1713 2954 : if (F2x_degree(gel(V,i)))
1714 : {
1715 1778 : GEN Fj = F2x_simplefact_Shoup(gel(V,i));
1716 1778 : gel(F, j) = Fj;
1717 1778 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1718 1778 : j++;
1719 : }
1720 1736 : return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
1721 : }
1722 :
1723 : static int
1724 238 : F2x_isirred_Cantor(GEN T)
1725 : {
1726 238 : pari_sp av = avma;
1727 : pari_timer ti;
1728 : long n, d;
1729 238 : GEN dT = F2x_deriv(T);
1730 : GEN XP, D;
1731 238 : if (F2x_degree(F2x_gcd(T, dT)) != 0) { avma = av; return 0; }
1732 168 : n = F2x_degree(T);
1733 168 : if (DEBUGLEVEL>=6) timer_start(&ti);
1734 168 : XP = F2x_Frobenius(T);
1735 168 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1736 168 : D = F2x_ddf_simple(T, XP);
1737 168 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf");
1738 168 : d = F2x_degree(gel(D, n));
1739 168 : avma = av; return d==n;
1740 : }
1741 :
1742 : static GEN
1743 30777 : F2x_factcantor_i(GEN f, long flag)
1744 : {
1745 30777 : long d = F2x_degree(f);
1746 30776 : if (d <= 2) return F2x_factor_deg2(f,d,flag);
1747 19590 : switch(flag)
1748 : {
1749 17616 : default: return F2x_factor_Cantor(f);
1750 1736 : case 1: return F2x_simplefact_Cantor(f);
1751 238 : case 2: return F2x_isirred_Cantor(f)? gen_1: NULL;
1752 : }
1753 : }
1754 :
1755 : GEN
1756 15174 : F2x_factcantor(GEN f, long flag)
1757 : {
1758 15174 : pari_sp av = avma;
1759 15174 : GEN z = F2x_factcantor_i(f, flag);
1760 15173 : if (flag == 2) { avma = av; return z; }
1761 15173 : return gerepilecopy(av, z);
1762 : }
1763 :
1764 : GEN
1765 0 : F2x_degfact(GEN f)
1766 : {
1767 0 : pari_sp av = avma;
1768 0 : GEN z = F2x_factcantor_i(f, 1);
1769 0 : return gerepilecopy(av, z);
1770 : }
1771 :
1772 : int
1773 238 : F2x_is_irred(GEN f) { return !!F2x_factcantor_i(f, 2); }
1774 :
1775 : /* Adapted from Shoup NTL */
1776 : GEN
1777 743742 : Flx_factor_squarefree(GEN f, ulong p)
1778 : {
1779 743742 : long i, q, n = degpol(f);
1780 743742 : GEN u = const_vec(n+1, pol1_Flx(f[1]));
1781 796906 : for(q = 1;;q *= p)
1782 : {
1783 796906 : GEN t, v, tv, r = Flx_gcd(f, Flx_deriv(f, p), p);
1784 796906 : if (degpol(r) == 0) { gel(u, q) = f; break; }
1785 89046 : t = Flx_div(f, r, p);
1786 89046 : if (degpol(t) > 0)
1787 : {
1788 : long j;
1789 112707 : for(j = 1;;j++)
1790 : {
1791 112707 : v = Flx_gcd(r, t, p);
1792 112707 : tv = Flx_div(t, v, p);
1793 112707 : if (degpol(tv) > 0) gel(u, j*q) = tv;
1794 112707 : if (degpol(v) <= 0) break;
1795 76358 : r = Flx_div(r, v, p);
1796 76358 : t = v;
1797 76358 : }
1798 36349 : if (degpol(r) == 0) break;
1799 : }
1800 53164 : f = Flx_deflate(r, p);
1801 53164 : }
1802 3820345 : for (i = n; i; i--)
1803 3817685 : if (degpol(gel(u,i))) break;
1804 743742 : setlg(u,i+1); return u;
1805 : }
1806 :
1807 : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
1808 : static GEN
1809 659712 : Flx_ddf(GEN T, GEN XP, ulong p)
1810 : {
1811 659712 : pari_sp av = avma;
1812 : GEN b, g, h, F, f, Tr, xq;
1813 : long i, j, n, v, bo, ro;
1814 : long B, l, m;
1815 : pari_timer ti;
1816 659712 : n = get_Flx_degree(T); v = get_Flx_var(T);
1817 659712 : if (n == 0) return cgetg(1, t_VEC);
1818 659712 : if (n == 1) return mkvec(get_Flx_mod(T));
1819 648148 : B = n/2;
1820 648148 : l = usqrt(B);
1821 648148 : m = (B+l-1)/l;
1822 648148 : T = Flx_get_red(T, p);
1823 648148 : b = cgetg(l+2, t_VEC);
1824 648148 : gel(b, 1) = polx_Flx(v);
1825 648148 : gel(b, 2) = XP;
1826 648148 : bo = brent_kung_optpow(n, l-1, 1);
1827 648148 : ro = l<=1 ? 0:(bo-1)/(l-1) + ((n-1)/bo);
1828 648148 : if (DEBUGLEVEL>=7) timer_start(&ti);
1829 648148 : if (expu(p) <= ro)
1830 151322 : for (i = 3; i <= l+1; i++)
1831 85398 : gel(b, i) = Flxq_powu(gel(b, i-1), p, T, p);
1832 : else
1833 : {
1834 582224 : xq = Flxq_powers(gel(b, 2), bo, T, p);
1835 582224 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: xq baby");
1836 681959 : for (i = 3; i <= l+1; i++)
1837 99735 : gel(b, i) = Flx_FlxqV_eval(gel(b, i-1), xq, T, p);
1838 : }
1839 648148 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: baby");
1840 648148 : xq = Flxq_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1), T, p);
1841 648148 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: xq giant");
1842 648148 : g = cgetg(m+1, t_VEC);
1843 648148 : gel(g, 1) = gel(xq, 2);
1844 1258203 : for(i = 2; i <= m; i++)
1845 610055 : gel(g, i) = Flx_FlxqV_eval(gel(g, i-1), xq, T, p);
1846 648148 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: giant");
1847 648148 : h = cgetg(m+1, t_VEC);
1848 1906351 : for (j = 1; j <= m; j++)
1849 : {
1850 1258203 : pari_sp av = avma;
1851 1258203 : GEN gj = gel(g, j);
1852 1258203 : GEN e = Flx_sub(gj, gel(b, 1), p);
1853 1929446 : for (i = 2; i <= l; i++)
1854 671245 : e = Flxq_mul(e, Flx_sub(gj, gel(b, i), p), T, p);
1855 1258201 : gel(h, j) = gerepileupto(av, e);
1856 : }
1857 648148 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: diff");
1858 648148 : Tr = get_Flx_mod(T);
1859 648148 : F = cgetg(m+1, t_VEC);
1860 1906350 : for (j = 1; j <= m; j++)
1861 : {
1862 1258203 : gel(F, j) = Flx_gcd(Tr, gel(h, j), p);
1863 1258203 : Tr = Flx_div(Tr, gel(F,j), p);
1864 : }
1865 648147 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: F");
1866 648147 : f = const_vec(n, pol1_Flx(v));
1867 1906350 : for (j = 1; j <= m; j++)
1868 : {
1869 1258203 : GEN e = gel(F, j);
1870 1396385 : for (i=l-1; i >= 0; i--)
1871 : {
1872 1396385 : GEN u = Flx_gcd(e, Flx_sub(gel(g, j), gel(b, i+1), p), p);
1873 1396384 : if (degpol(u))
1874 : {
1875 528076 : gel(f, l*j-i) = u;
1876 528076 : e = Flx_div(e, u, p);
1877 : }
1878 1396384 : if (!degpol(e)) break;
1879 : }
1880 : }
1881 648147 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: f");
1882 648147 : if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
1883 648147 : return gerepilecopy(av, f);
1884 : }
1885 :
1886 : static void
1887 45736 : Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
1888 : {
1889 45736 : long n = degpol(Tp), r = n/d;
1890 : GEN T, f, ff;
1891 : ulong p2;
1892 91472 : if (r==1) { gel(V, idx) = Tp; return; }
1893 20646 : p2 = p>>1;
1894 20646 : T = Flx_get_red(Tp, p);
1895 20646 : XP = Flx_rem(XP, T, p);
1896 : while (1)
1897 : {
1898 22517 : pari_sp btop = avma;
1899 : long i;
1900 22517 : GEN g = random_Flx(n, Tp[1], p);
1901 22517 : GEN t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
1902 22517 : if (lgpol(t) == 0) continue;
1903 46031 : for(i=1; i<=10; i++)
1904 : {
1905 44616 : pari_sp btop2 = avma;
1906 44616 : GEN R = Flxq_powu(Flx_Fl_add(t, random_Fl(p), p), p2, T, p);
1907 44616 : f = Flx_gcd(Flx_Fl_add(R, p-1, p), Tp, p);
1908 44616 : if (degpol(f) > 0 && degpol(f) < n) break;
1909 23970 : avma = btop2;
1910 : }
1911 22061 : if (degpol(f) > 0 && degpol(f) < n) break;
1912 1415 : avma = btop;
1913 1871 : }
1914 20646 : f = Flx_normalize(f, p);
1915 20646 : ff = Flx_div(Tp, f ,p);
1916 20646 : Flx_edf_simple(f, XP, d, p, V, idx);
1917 20646 : Flx_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
1918 : }
1919 : static void
1920 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx);
1921 :
1922 : static void
1923 32747 : Flx_edf_rec(GEN T, GEN XP, GEN hp, GEN t, long d, ulong p, GEN V, long idx)
1924 : {
1925 : pari_sp av;
1926 32747 : GEN Tp = get_Flx_mod(T);
1927 32747 : long n = degpol(hp), vT = Tp[1];
1928 : GEN u1, u2, f1, f2;
1929 32747 : ulong p2 = p>>1;
1930 : GEN R, h;
1931 32747 : h = Flx_get_red(hp, p);
1932 32747 : t = Flx_rem(t, T, p);
1933 32747 : av = avma;
1934 : do
1935 : {
1936 56331 : avma = av;
1937 56331 : R = Flxq_powu(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p);
1938 56331 : u1 = Flx_gcd(Flx_Fl_add(R, p-1, p), hp, p);
1939 56331 : } while (degpol(u1)==0 || degpol(u1)==n);
1940 32747 : f1 = Flx_gcd(Flx_Flxq_eval(u1, t, T, p), Tp, p);
1941 32747 : f1 = Flx_normalize(f1, p);
1942 32747 : u2 = Flx_div(hp, u1, p);
1943 32747 : f2 = Flx_div(Tp, f1, p);
1944 32747 : if (degpol(u1)==1)
1945 : {
1946 25748 : if (degpol(f1)==d)
1947 25390 : gel(V, idx) = f1;
1948 : else
1949 358 : Flx_edf(f1, XP, d, p, V, idx);
1950 : }
1951 : else
1952 6999 : Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, V, idx);
1953 32747 : idx += degpol(f1)/d;
1954 32747 : if (degpol(u2)==1)
1955 : {
1956 25909 : if (degpol(f2)==d)
1957 25469 : gel(V, idx) = f2;
1958 : else
1959 440 : Flx_edf(f2, XP, d, p, V, idx);
1960 : }
1961 : else
1962 6838 : Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, V, idx);
1963 32747 : }
1964 :
1965 : static void
1966 18910 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
1967 : {
1968 18910 : long n = degpol(Tp), r = n/d, vT = Tp[1];
1969 : GEN T, h, t;
1970 : pari_timer ti;
1971 37820 : if (r==1) { gel(V, idx) = Tp; return; }
1972 18910 : T = Flx_get_red(Tp, p);
1973 18910 : XP = Flx_rem(XP, T, p);
1974 18910 : if (DEBUGLEVEL>=7) timer_start(&ti);
1975 : do
1976 : {
1977 20914 : GEN g = random_Flx(n, vT, p);
1978 20914 : t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
1979 20914 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
1980 20914 : h = Flxq_minpoly(t, T, p);
1981 20914 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
1982 20914 : } while (degpol(h) <= 1);
1983 18910 : Flx_edf_rec(T, XP, h, t, d, p, V, idx);
1984 : }
1985 :
1986 : static GEN
1987 32834 : Flx_factor_Shoup(GEN T, ulong p)
1988 : {
1989 32834 : long i, n, s = 0;
1990 : GEN XP, D, V;
1991 32834 : long e = expu(p);
1992 : pari_timer ti;
1993 32834 : n = get_Flx_degree(T);
1994 32834 : T = Flx_get_red(T, p);
1995 32834 : if (DEBUGLEVEL>=6) timer_start(&ti);
1996 32834 : XP = Flx_Frobenius(T, p);
1997 32834 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
1998 32834 : D = Flx_ddf(T, XP, p);
1999 32834 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
2000 336772 : for (i = 1; i <= n; i++)
2001 303938 : s += degpol(gel(D,i))/i;
2002 32834 : V = cgetg(s+1, t_COL);
2003 336770 : for (i = 1, s = 1; i <= n; i++)
2004 : {
2005 303936 : GEN Di = gel(D,i);
2006 303936 : long ni = degpol(Di), ri = ni/i;
2007 303936 : if (ni == 0) continue;
2008 49489 : Di = Flx_normalize(Di, p);
2009 49488 : if (ni == i) { gel(V, s++) = Di; continue; }
2010 22556 : if (ri <= e*expu(e))
2011 18112 : Flx_edf(Di, XP, i, p, V, s);
2012 : else
2013 4444 : Flx_edf_simple(Di, XP, i, p, V, s);
2014 22556 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
2015 22556 : s += ri;
2016 : }
2017 32834 : return V;
2018 : }
2019 :
2020 : static GEN
2021 352116 : Flx_simplefact_Shoup(GEN T, ulong p)
2022 : {
2023 : GEN XP, D;
2024 : pari_timer ti;
2025 352116 : long n = get_Flx_degree(T);
2026 352116 : T = Flx_get_red(T, p);
2027 352116 : if (DEBUGLEVEL>=6) timer_start(&ti);
2028 352116 : XP = Flx_Frobenius(T, p);
2029 352116 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2030 352116 : D = Flx_ddf(T, XP, p);
2031 352116 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
2032 352116 : return ddf_to_simplefact(D, n);
2033 : }
2034 :
2035 : static GEN
2036 29727 : Flx_factor_Cantor(GEN T, ulong p)
2037 : {
2038 29727 : GEN E, F, V = Flx_factor_squarefree(get_Flx_mod(T), p);
2039 29727 : long i, j, l = lg(V);
2040 29727 : F = cgetg(l, t_VEC);
2041 29727 : E = cgetg(l, t_VEC);
2042 65491 : for (i=1, j=1; i < l; i++)
2043 35764 : if (degpol(gel(V,i)))
2044 : {
2045 32834 : GEN Fj = Flx_factor_Shoup(gel(V,i), p);
2046 32834 : gel(F, j) = Fj;
2047 32834 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
2048 32834 : j++;
2049 : }
2050 29727 : return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
2051 : }
2052 :
2053 : static GEN
2054 351134 : Flx_simplefact_Cantor(GEN T, ulong p)
2055 : {
2056 351134 : GEN E, F, V = Flx_factor_squarefree(get_Flx_mod(T), p);
2057 351134 : long i, j, l = lg(V);
2058 351134 : F = cgetg(l, t_VEC);
2059 351134 : E = cgetg(l, t_VEC);
2060 705253 : for (i=1, j=1; i < l; i++)
2061 354119 : if (degpol(gel(V,i)))
2062 : {
2063 352116 : GEN Fj = Flx_simplefact_Shoup(gel(V,i), p);
2064 352116 : gel(F, j) = Fj;
2065 352116 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
2066 352116 : j++;
2067 : }
2068 351134 : return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
2069 : }
2070 :
2071 : static int
2072 581 : Flx_isirred_Cantor(GEN Tp, ulong p)
2073 : {
2074 581 : pari_sp av = avma;
2075 : pari_timer ti;
2076 : long n, d;
2077 581 : GEN T = get_Flx_mod(Tp);
2078 581 : GEN dT = Flx_deriv(T, p);
2079 : GEN XP, D;
2080 581 : if (degpol(Flx_gcd(T, dT, p)) != 0) { avma = av; return 0; }
2081 441 : n = get_Flx_degree(T);
2082 441 : T = Flx_get_red(Tp, p);
2083 441 : if (DEBUGLEVEL>=6) timer_start(&ti);
2084 441 : XP = Flx_Frobenius(T, p);
2085 441 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2086 441 : D = Flx_ddf(T, XP, p);
2087 441 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
2088 441 : d = degpol(gel(D, n));
2089 441 : avma = av; return d==n;
2090 : }
2091 :
2092 : static GEN
2093 1159595 : Flx_factcantor_i(GEN f, ulong pp, long flag)
2094 : {
2095 : long d;
2096 1159595 : if (pp==2) { /*We need to handle 2 specially */
2097 15330 : GEN F = F2x_factcantor_i(Flx_to_F2x(f),flag);
2098 15330 : if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
2099 15330 : return F;
2100 : }
2101 1144265 : d = degpol(f);
2102 1144265 : if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
2103 381442 : switch(flag)
2104 : {
2105 29727 : default: return Flx_factor_Cantor(f, pp);
2106 351134 : case 1: return Flx_simplefact_Cantor(f, pp);
2107 581 : case 2: return Flx_isirred_Cantor(f, pp)? gen_1: NULL;
2108 : }
2109 : }
2110 :
2111 : GEN
2112 7610 : Flx_factcantor(GEN f, ulong p, long flag)
2113 : {
2114 7610 : pari_sp av = avma;
2115 7610 : GEN z = Flx_factcantor_i(Flx_normalize(f,p),p,flag);
2116 7610 : if (flag == 2) { avma = av; return z; }
2117 7610 : return gerepilecopy(av, z);
2118 : }
2119 :
2120 : GEN
2121 1118187 : Flx_degfact(GEN f, ulong p)
2122 : {
2123 1118187 : pari_sp av = avma;
2124 1118187 : GEN z = Flx_factcantor_i(Flx_normalize(f,p),p,1);
2125 1118187 : return gerepilecopy(av, z);
2126 : }
2127 :
2128 : /* T must be squarefree mod p*/
2129 : GEN
2130 192071 : Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
2131 : {
2132 : GEN XP, D;
2133 : pari_timer ti;
2134 192071 : long i, s, n = get_Flx_degree(T);
2135 192071 : GEN V = const_vecsmall(n, 0);
2136 192071 : pari_sp av = avma;
2137 192071 : T = Flx_get_red(T, p);
2138 192071 : if (DEBUGLEVEL>=6) timer_start(&ti);
2139 192071 : XP = Flx_Frobenius(T, p);
2140 192071 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2141 192071 : D = Flx_ddf(T, XP, p);
2142 192071 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
2143 1220094 : for (i = 1, s = 0; i <= n; i++)
2144 : {
2145 1028023 : V[i] = degpol(gel(D,i))/i;
2146 1028023 : s += V[i];
2147 : }
2148 192071 : *nb = s;
2149 192071 : avma = av; return V;
2150 : }
2151 :
2152 : long
2153 82250 : Flx_nbfact_Frobenius(GEN T, GEN XP, ulong p)
2154 : {
2155 82250 : pari_sp av = avma;
2156 82250 : long s = ddf_to_nbfact(Flx_ddf(T, XP, p));
2157 82250 : avma = av; return s;
2158 : }
2159 :
2160 : /* T must be squarefree mod p*/
2161 : long
2162 82250 : Flx_nbfact(GEN T, ulong p)
2163 : {
2164 82250 : pari_sp av = avma;
2165 82250 : GEN XP = Flx_Frobenius(T, p);
2166 82250 : long n = Flx_nbfact_Frobenius(T, XP, p);
2167 82250 : avma = av; return n;
2168 : }
2169 :
2170 : int
2171 581 : Flx_is_irred(GEN f, ulong p) { return !!Flx_factcantor_i(f,p,2); }
2172 :
2173 : /* factor f (FpX or Flx) mod pp.
2174 : * flag = 1: return the degrees, not the factors
2175 : * flag = 2: return NULL if f is not irreducible.
2176 : * Not gerepile-safe */
2177 : static GEN
2178 112 : factcantor_i(GEN f, GEN pp, long flag)
2179 : {
2180 112 : if (typ(f) == t_VECSMALL)
2181 : { /* lgefint(pp) = 3 */
2182 : GEN F;
2183 63 : ulong p = pp[2];
2184 63 : if (p==2) {
2185 35 : F = F2x_factcantor_i(Flx_to_F2x(f),flag);
2186 35 : if (flag==0) F2xV_to_ZXV_inplace(gel(F,1));
2187 : } else {
2188 28 : F = Flx_factcantor_i(f,p,flag);
2189 28 : if (flag==0) FlxV_to_ZXV_inplace(gel(F,1));
2190 : }
2191 63 : return F;
2192 : }
2193 49 : return FpX_factcantor_i(f, pp, flag);
2194 : }
2195 : GEN
2196 0 : FpX_factcantor(GEN f, GEN pp, long flag)
2197 : {
2198 0 : pari_sp av = avma;
2199 : GEN z;
2200 0 : ZX_factmod_init(&f,pp);
2201 0 : z = factcantor_i(f,pp,flag);
2202 0 : if (flag == 2) { avma = av; return z; }
2203 0 : return gerepilecopy(av, z);
2204 : }
2205 :
2206 : static GEN
2207 266 : factmod_aux(GEN f, GEN p, GEN (*Factor)(GEN,GEN,long), long flag)
2208 : {
2209 266 : pari_sp av = avma;
2210 : long j, lfact;
2211 : GEN z, t, E, y, u, v;
2212 :
2213 266 : factmod_init(&f, p);
2214 266 : switch(lg(f))
2215 : {
2216 49 : case 3: avma = av; return trivial_fact();
2217 14 : case 2: return gerepileupto(av, zero_fact_intmod(f, p));
2218 : }
2219 203 : z = Factor(f,p,flag); t = gel(z,1); E = gel(z,2);
2220 203 : lfact = lg(t); y = cgetg(3, t_MAT);
2221 203 : gel(y,1) = u = cgetg(lfact,t_COL);
2222 203 : gel(y,2) = v = cgetg(lfact,t_COL);
2223 203 : if (flag)
2224 119 : for (j=1; j<lfact; j++)
2225 : {
2226 77 : gel(u,j) = utoi(uel(t,j));
2227 77 : gel(v,j) = utoi(uel(E,j));
2228 : }
2229 : else
2230 8078 : for (j=1; j<lfact; j++)
2231 : {
2232 7917 : gel(u,j) = FpX_to_mod(gel(t,j), p);
2233 7917 : gel(v,j) = utoi(uel(E,j));
2234 : }
2235 203 : return gerepileupto(av, y);
2236 : }
2237 : GEN
2238 112 : factcantor0(GEN f, GEN p, long flag)
2239 112 : { return factmod_aux(f, p, &factcantor_i, flag); }
2240 : GEN
2241 154 : factmod(GEN f, GEN p)
2242 154 : { return factmod_aux(f, p, &FpX_Berlekamp_i, 0); }
2243 :
2244 : /* Use this function when you think f is reducible, and that there are lots of
2245 : * factors. If you believe f has few factors, use FpX_nbfact(f,p)==1 instead */
2246 : int
2247 14 : FpX_is_irred(GEN f, GEN p) {
2248 14 : ZX_factmod_init(&f,p);
2249 14 : return !!factcantor_i(f,p,2);
2250 : }
2251 : GEN
2252 0 : FpX_degfact(GEN f, GEN p) {
2253 0 : pari_sp av = avma;
2254 : GEN z;
2255 0 : ZX_factmod_init(&f,p);
2256 0 : z = factcantor_i(f,p,1);
2257 0 : return gerepilecopy(av, z);
2258 : }
2259 : GEN
2260 70 : factcantor(GEN f, GEN p) { return factcantor0(f,p,0); }
2261 : GEN
2262 42 : simplefactmod(GEN f, GEN p) { return factcantor0(f,p,1); }
2263 :
2264 : /* set x <-- x + c*y mod p */
2265 : /* x is not required to be normalized.*/
2266 : static void
2267 1254947 : Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
2268 : {
2269 : long i, lx, ly;
2270 1254947 : ulong *x=(ulong *)gx;
2271 1254947 : ulong *y=(ulong *)gy;
2272 2509894 : if (!c) return;
2273 1221181 : lx = lg(gx);
2274 1221181 : ly = lg(gy);
2275 1221181 : if (lx<ly) pari_err_BUG("lx<ly in Flx_addmul_inplace");
2276 1221181 : if (SMALL_ULONG(p))
2277 1218680 : for (i=2; i<ly; i++) x[i] = (x[i] + c*y[i]) % p;
2278 : else
2279 2501 : for (i=2; i<ly; i++) x[i] = Fl_add(x[i], Fl_mul(c,y[i],p),p);
2280 : }
2281 :
2282 : #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
2283 : /* assume x1 != 0 */
2284 : static GEN
2285 154010 : deg1_Flx(ulong x1, ulong x0, ulong sv)
2286 : {
2287 154010 : return mkvecsmall3(sv, x0, x1);
2288 : }
2289 :
2290 : static long
2291 52878 : F2x_split_Berlekamp(GEN *t)
2292 : {
2293 52878 : GEN u = *t, a, b, vker;
2294 52878 : long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
2295 :
2296 52878 : if (du == 1) return 1;
2297 26014 : if (du == 2)
2298 : {
2299 6573 : if (F2x_quad_factortype(u) == 1) /* 0 is a root: shouldn't occur */
2300 : {
2301 0 : t[0] = mkvecsmall2(sv, 2);
2302 0 : t[1] = mkvecsmall2(sv, 3);
2303 0 : return 2;
2304 : }
2305 6573 : return 1;
2306 : }
2307 :
2308 19441 : vker = F2x_Berlekamp_ker(u);
2309 19441 : lb = lgcols(vker);
2310 19441 : d = lg(vker)-1;
2311 19441 : ir = 0;
2312 : /* t[i] irreducible for i < ir, still to be treated for i < L */
2313 46092 : for (L=1; L<d; )
2314 : {
2315 : GEN pol;
2316 7210 : if (d == 2)
2317 6139 : pol = F2v_to_F2x(gel(vker,2), sv);
2318 : else
2319 : {
2320 1071 : GEN v = zero_zv(lb);
2321 1071 : v[1] = du;
2322 1071 : v[2] = random_Fl(2); /*Assume vker[1]=1*/
2323 3465 : for (i=2; i<=d; i++)
2324 2394 : if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
2325 1071 : pol = F2v_to_F2x(v, sv);
2326 : }
2327 15435 : for (i=ir; i<L && L<d; i++)
2328 : {
2329 8225 : a = t[i]; la = F2x_degree(a);
2330 8225 : if (la == 1) { set_irred(i); }
2331 8204 : else if (la == 2)
2332 : {
2333 14 : if (F2x_quad_factortype(a) == 1) /* 0 is a root: shouldn't occur */
2334 : {
2335 0 : t[i] = mkvecsmall2(sv, 2);
2336 0 : t[L] = mkvecsmall2(sv, 3); L++;
2337 : }
2338 14 : set_irred(i);
2339 : }
2340 : else
2341 : {
2342 8190 : pari_sp av = avma;
2343 : long lb;
2344 8190 : b = F2x_rem(pol, a);
2345 8190 : if (F2x_degree(b) <= 0) { avma=av; continue; }
2346 6748 : b = F2x_gcd(a,b); lb = F2x_degree(b);
2347 6748 : if (lb && lb < la)
2348 : {
2349 6748 : t[L] = F2x_div(a,b);
2350 6748 : t[i]= b; L++;
2351 : }
2352 0 : else avma = av;
2353 : }
2354 : }
2355 : }
2356 19441 : return d;
2357 : }
2358 :
2359 : /* p != 2 */
2360 : static long
2361 372098 : Flx_split_Berlekamp(GEN *t, ulong p)
2362 : {
2363 372098 : GEN u = *t, a,b,vker;
2364 372098 : long d, i, ir, L, la, lb, sv = u[1];
2365 372098 : long l = lg(u);
2366 : ulong po2;
2367 :
2368 372098 : if (p == 2)
2369 : {
2370 0 : *t = Flx_to_F2x(*t);
2371 0 : d = F2x_split_Berlekamp(t);
2372 0 : for (i = 1; i <= d; i++) t[i] = F2x_to_Flx(t[i]);
2373 0 : return d;
2374 : }
2375 372098 : la = degpol(u);
2376 372098 : if (la == 1) return 1;
2377 297214 : if (la == 2)
2378 : {
2379 14998 : ulong r = Flx_quad_root(u,p,1);
2380 14998 : if (r != p)
2381 : {
2382 7151 : t[0] = deg1_Flx(1, p - r, sv); r = Flx_otherroot(u,r,p);
2383 7151 : t[1] = deg1_Flx(1, p - r, sv);
2384 7151 : return 2;
2385 : }
2386 7847 : return 1;
2387 : }
2388 :
2389 282216 : vker = Flx_Berlekamp_ker(u,p);
2390 282216 : vker = Flm_to_FlxV(vker, sv);
2391 282216 : d = lg(vker)-1;
2392 282216 : po2 = p >> 1; /* (p-1) / 2 */
2393 282216 : ir = 0;
2394 : /* t[i] irreducible for i < ir, still to be treated for i < L */
2395 1063230 : for (L=1; L<d; )
2396 : {
2397 498798 : GEN pol = zero_zv(l-2);
2398 498798 : pol[1] = sv;
2399 498798 : pol[2] = random_Fl(p); /*Assume vker[1]=1*/
2400 1753745 : for (i=2; i<=d; i++)
2401 1254947 : Flx_addmul_inplace(pol, gel(vker,i), random_Fl(p), p);
2402 498798 : (void)Flx_renormalize(pol,l-1);
2403 :
2404 1301667 : for (i=ir; i<L && L<d; i++)
2405 : {
2406 802869 : a = t[i]; la = degpol(a);
2407 802869 : if (la == 1) { set_irred(i); }
2408 744154 : else if (la == 2)
2409 : {
2410 130198 : ulong r = Flx_quad_root(a,p,1);
2411 130198 : if (r != p)
2412 : {
2413 69854 : t[i] = deg1_Flx(1, p - r, sv); r = Flx_otherroot(a,r,p);
2414 69854 : t[L] = deg1_Flx(1, p - r, sv); L++;
2415 : }
2416 130198 : set_irred(i);
2417 : }
2418 : else
2419 : {
2420 613956 : pari_sp av = avma;
2421 613956 : b = Flx_rem(pol, a, p);
2422 613956 : if (degpol(b) <= 0) { avma=av; continue; }
2423 563582 : b = Flx_Fl_add(Flxq_powu(b,po2, a,p), p-1, p);
2424 563582 : b = Flx_gcd(a,b, p); lb = degpol(b);
2425 563582 : if (lb && lb < la)
2426 : {
2427 313731 : b = Flx_normalize(b, p);
2428 313731 : t[L] = Flx_div(a,b,p);
2429 313731 : t[i]= b; L++;
2430 : }
2431 249851 : else avma = av;
2432 : }
2433 : }
2434 : }
2435 282216 : return d;
2436 : }
2437 :
2438 : static long
2439 66 : FpX_split_Berlekamp(GEN *t, GEN p)
2440 : {
2441 66 : GEN u = *t, a,b,po2,vker;
2442 66 : long d, i, ir, L, la, lb, vu = varn(u);
2443 66 : if (lgefint(p) == 3)
2444 : {
2445 0 : ulong up = p[2];
2446 0 : if (up == 2)
2447 : {
2448 0 : *t = ZX_to_F2x(*t);
2449 0 : d = F2x_split_Berlekamp(t);
2450 0 : for (i = 0; i < d; i++) t[i] = F2x_to_ZX(t[i]);
2451 : }
2452 : else
2453 : {
2454 0 : *t = ZX_to_Flx(*t, up);
2455 0 : d = Flx_split_Berlekamp(t, up);
2456 0 : for (i = 0; i < d; i++) t[i] = Flx_to_ZX(t[i]);
2457 : }
2458 0 : return d;
2459 : }
2460 66 : la = degpol(u);
2461 66 : if (la == 1) return 1;
2462 62 : if (la == 2)
2463 : {
2464 2 : GEN r = FpX_quad_root(u,p,1);
2465 2 : if (r)
2466 : {
2467 0 : t[0] = deg1pol_shallow(gen_1, subii(p,r), vu); r = FpX_otherroot(u,r,p);
2468 0 : t[1] = deg1pol_shallow(gen_1, subii(p,r), vu);
2469 0 : return 2;
2470 : }
2471 2 : return 1;
2472 : }
2473 60 : vker = FpX_Berlekamp_ker(u,p);
2474 60 : vker = RgM_to_RgXV(vker,vu);
2475 60 : d = lg(vker)-1;
2476 60 : po2 = shifti(p, -1); /* (p-1) / 2 */
2477 60 : ir = 0;
2478 : /* t[i] irreducible for i < ir, still to be treated for i < L */
2479 164 : for (L=1; L<d; )
2480 : {
2481 44 : GEN pol = scalar_ZX_shallow(randomi(p), vu);
2482 176 : for (i=2; i<=d; i++)
2483 132 : pol = ZX_add(pol, ZX_Z_mul(gel(vker,i), randomi(p)));
2484 44 : pol = FpX_red(pol,p);
2485 124 : for (i=ir; i<L && L<d; i++)
2486 : {
2487 80 : a = t[i]; la = degpol(a);
2488 80 : if (la == 1) { set_irred(i); }
2489 71 : else if (la == 2)
2490 : {
2491 9 : GEN r = FpX_quad_root(a,p,1);
2492 9 : if (r)
2493 : {
2494 9 : t[i] = deg1pol_shallow(gen_1, subii(p,r), vu); r = FpX_otherroot(a,r,p);
2495 9 : t[L] = deg1pol_shallow(gen_1, subii(p,r), vu); L++;
2496 : }
2497 9 : set_irred(i);
2498 : }
2499 : else
2500 : {
2501 62 : pari_sp av = avma;
2502 62 : b = FpX_rem(pol, a, p);
2503 62 : if (degpol(b) <= 0) { avma=av; continue; }
2504 62 : b = FpX_Fp_sub_shallow(FpXQ_pow(b,po2, a,p), gen_1, p);
2505 62 : b = FpX_gcd(a,b, p); lb = degpol(b);
2506 62 : if (lb && lb < la)
2507 : {
2508 49 : b = FpX_normalize(b, p);
2509 49 : t[L] = FpX_div(a,b,p);
2510 49 : t[i]= b; L++;
2511 : }
2512 13 : else avma = av;
2513 : }
2514 : }
2515 : }
2516 60 : return d;
2517 : }
2518 :
2519 : static GEN
2520 152964 : F2x_Berlekamp_i(GEN f, long flag)
2521 : {
2522 152964 : long lfact, val, d = F2x_degree(f), j, k, lV;
2523 : GEN y, E, t, V;
2524 :
2525 152964 : if (d <= 2) return F2x_factor_deg2(f, d, flag);
2526 :
2527 57211 : val = F2x_valrem(f, &f);
2528 57211 : if (flag == 2 && val > 1) return NULL;
2529 57211 : V = F2x_factor_squarefree(f); lV = lg(V);
2530 57211 : if (flag == 2 && lV > 2) return NULL;
2531 :
2532 : /* to hold factors and exponents */
2533 57211 : t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
2534 57211 : E = cgetg(d+1,t_VECSMALL);
2535 57211 : lfact = 1;
2536 57211 : if (val) {
2537 18639 : if (flag == 1)
2538 0 : t[1] = 1;
2539 : else
2540 18639 : gel(t,1) = polx_F2x(f[1]);
2541 18639 : E[1] = val; lfact++;
2542 : }
2543 :
2544 221822 : for (k=1; k<lV; k++)
2545 : {
2546 164611 : if (F2x_degree(gel(V, k))==0) continue;
2547 52878 : gel(t,lfact) = gel(V, k);
2548 52878 : d = F2x_split_Berlekamp(&gel(t,lfact));
2549 52878 : if (flag == 2 && d != 1) return NULL;
2550 52878 : if (flag == 1)
2551 0 : for (j=0; j<d; j++) t[lfact+j] = F2x_degree(gel(t,lfact+j));
2552 52878 : for (j=0; j<d; j++) E[lfact+j] = k;
2553 52878 : lfact += d;
2554 : }
2555 57211 : if (flag == 2) return gen_1; /* irreducible */
2556 57211 : y = FE_setlg(t,E, lfact);
2557 57211 : return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
2558 57211 : : sort_factor_pol(y, cmpGuGu);
2559 : }
2560 :
2561 : static GEN
2562 783624 : Flx_Berlekamp_i(GEN f, ulong p, long flag)
2563 : {
2564 783624 : long lfact, val, d = degpol(f), j, k, lV;
2565 : GEN y, E, t, V;
2566 :
2567 783624 : if (p == 2)
2568 : {
2569 11865 : GEN F = F2x_Berlekamp_i(Flx_to_F2x(f),flag);
2570 11865 : if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
2571 11865 : return F;
2572 : }
2573 771759 : if (d <= 2) return Flx_factor_deg2(f,p,d,flag);
2574 362363 : val = Flx_valrem(f, &f);
2575 362363 : if (flag == 2 && val > 1) return NULL;
2576 362363 : V = Flx_factor_squarefree(f, p); lV = lg(V);
2577 362363 : if (flag == 2 && lV > 2) return NULL;
2578 :
2579 : /* to hold factors and exponents */
2580 362363 : t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
2581 362363 : E = cgetg(d+1,t_VECSMALL);
2582 362363 : lfact = 1;
2583 362363 : if (val) {
2584 14713 : if (flag == 1)
2585 0 : t[1] = 1;
2586 : else
2587 14713 : gel(t,1) = polx_Flx(f[1]);
2588 14713 : E[1] = val; lfact++;
2589 : }
2590 :
2591 908546 : for (k=1; k<lV; k++)
2592 : {
2593 546183 : if (degpol(gel(V, k))==0) continue;
2594 372098 : gel(t,lfact) = Flx_normalize(gel(V, k), p);
2595 372098 : d = Flx_split_Berlekamp(&gel(t,lfact), p);
2596 372098 : if (flag == 2 && d != 1) return NULL;
2597 372098 : if (flag == 1)
2598 0 : for (j=0; j<d; j++) t[lfact+j] = degpol(gel(t,lfact+j));
2599 372098 : for (j=0; j<d; j++) E[lfact+j] = k;
2600 372098 : lfact += d;
2601 : }
2602 362363 : if (flag == 2) return gen_1; /* irreducible */
2603 362363 : y = FE_setlg(t,E,lfact);
2604 362363 : return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
2605 362363 : : sort_factor_pol(y, cmpGuGu);
2606 : }
2607 :
2608 : /* f an FpX or an Flx */
2609 : static GEN
2610 856954 : FpX_Berlekamp_i(GEN f, GEN p, long flag)
2611 : {
2612 856954 : long lfact, val, d = degpol(f), j, k, lV;
2613 : GEN y, E, t ,V;
2614 :
2615 856954 : if (typ(f) == t_VECSMALL)
2616 : {/* lgefint(p) == 3 */
2617 856809 : ulong pp = p[2];
2618 : GEN F;
2619 856809 : if (pp == 2) {
2620 141085 : F = F2x_Berlekamp_i(Flx_to_F2x(f), flag);
2621 141085 : if (flag==0) F2xV_to_ZXV_inplace(gel(F,1));
2622 : } else {
2623 715724 : F = Flx_Berlekamp_i(f, pp, flag);
2624 715724 : if (flag==0) FlxV_to_ZXV_inplace(gel(F,1));
2625 : }
2626 856809 : return F;
2627 : }
2628 : /* p is large (and odd) */
2629 145 : if (d <= 2) return FpX_factor_deg2(f, p, d, flag);
2630 62 : val = ZX_valrem(f, &f);
2631 62 : if (flag == 2 && val > 1) return NULL;
2632 62 : V = FpX_factor_Yun(f, p); lV = lg(V);
2633 62 : if (flag == 2 && lg(V) > 2) return NULL;
2634 :
2635 : /* to hold factors and exponents */
2636 62 : t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
2637 62 : E = cgetg(d+1,t_VECSMALL);
2638 62 : lfact = 1;
2639 62 : if (val) {
2640 2 : if (flag == 1)
2641 0 : t[1] = 1;
2642 : else
2643 2 : gel(t,1) = pol_x(varn(f));
2644 2 : E[1] = val; lfact++;
2645 : }
2646 :
2647 129 : for (k=1; k<lV; k++)
2648 : {
2649 67 : if (degpol(gel(V,k))==0) continue;
2650 66 : gel(t,lfact) = FpX_normalize(gel(V, k), p);
2651 66 : d = FpX_split_Berlekamp(&gel(t,lfact), p);
2652 66 : if (flag == 2 && d != 1) return NULL;
2653 66 : if (flag == 1)
2654 0 : for (j=0; j<d; j++) t[lfact+j] = degpol(gel(t,lfact+j));
2655 66 : for (j=0; j<d; j++) E[lfact+j] = k;
2656 66 : lfact += d;
2657 : }
2658 62 : if (flag == 2) return gen_1; /* irreducible */
2659 62 : y = FE_setlg(t,E, lfact);
2660 62 : return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
2661 62 : : sort_factor_pol(y, cmpii);
2662 : }
2663 : GEN
2664 856849 : FpX_factor(GEN f, GEN p)
2665 : {
2666 856849 : pari_sp av = avma;
2667 856849 : ZX_factmod_init(&f, p);
2668 856849 : return gerepilecopy(av, FpX_Berlekamp_i(f, p, 0));
2669 : }
2670 : GEN
2671 101089 : Flx_factor(GEN f, ulong p)
2672 : {
2673 101089 : pari_sp av = avma;
2674 101089 : GEN F = (degpol(f)>log2(p))? Flx_factcantor_i(f,p,0): Flx_Berlekamp_i(f,p,0);
2675 101089 : return gerepilecopy(av, F);
2676 : }
2677 : GEN
2678 14 : F2x_factor(GEN f)
2679 : {
2680 14 : pari_sp av = avma;
2681 14 : return gerepilecopy(av, F2x_Berlekamp_i(f, 0));
2682 : }
2683 :
2684 : GEN
2685 126 : factormod0(GEN f, GEN p, long flag)
2686 : {
2687 126 : switch(flag)
2688 : {
2689 84 : case 0: return factmod(f,p);
2690 42 : case 1: return simplefactmod(f,p);
2691 0 : default: pari_err_FLAG("factormod");
2692 : }
2693 : return NULL; /* LCOV_EXCL_LINE */
2694 : }
|
