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 874607 : ZX_factmod_init(GEN *F, GEN p)
32 : {
33 874607 : if (lgefint(p) == 3)
34 : {
35 873222 : ulong pp = p[2];
36 873222 : if (pp == 2) { *F = ZX_to_F2x(*F); return 0; }
37 728520 : *F = ZX_to_Flx(*F, pp);
38 728520 : if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
39 728520 : return 1;
40 : }
41 1385 : *F = FpX_red(*F, p);
42 1385 : if (lg(*F) > 3) *F = FpX_normalize(*F, p);
43 1385 : return 2;
44 : }
45 : static void
46 83252 : ZX_rootmod_init(GEN *F, GEN p)
47 : {
48 83252 : if (lgefint(p) == 3)
49 : {
50 75280 : ulong pp = p[2];
51 75280 : *F = ZX_to_Flx(*F, pp);
52 75280 : if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
53 : }
54 : else
55 : {
56 7972 : *F = FpX_red(*F, p);
57 7972 : if (lg(*F) > 3) *F = FpX_normalize(*F, p);
58 : }
59 83252 : }
60 :
61 : /* return 1,...,p-1 [not_0 = 1] or 0,...,p [not_0 = 0] */
62 : static GEN
63 42 : all_roots_mod_p(ulong p, int not_0)
64 : {
65 : GEN r;
66 : ulong i;
67 42 : if (not_0) {
68 28 : r = cgetg(p, t_VECSMALL);
69 28 : for (i = 1; i < p; i++) r[i] = i;
70 : } else {
71 14 : r = cgetg(p+1, t_VECSMALL);
72 14 : for (i = 0; i < p; i++) r[i+1] = i;
73 : }
74 42 : return r;
75 : }
76 :
77 : /* X^n - 1 */
78 : static GEN
79 126 : Flx_Xnm1(long sv, long n, ulong p)
80 : {
81 126 : GEN t = cgetg(n+3, t_VECSMALL);
82 : long i;
83 126 : t[1] = sv;
84 126 : t[2] = p - 1;
85 126 : for (i = 3; i <= n+1; i++) t[i] = 0;
86 126 : t[i] = 1; return t;
87 : }
88 : /* X^n + 1 */
89 : static GEN
90 140 : Flx_Xn1(long sv, long n, ulong p)
91 : {
92 140 : GEN t = cgetg(n+3, t_VECSMALL);
93 : long i;
94 : (void) p;
95 140 : t[1] = sv;
96 140 : t[2] = 1;
97 140 : for (i = 3; i <= n+1; i++) t[i] = 0;
98 140 : t[i] = 1; return t;
99 : }
100 :
101 : static ulong
102 28 : Fl_nonsquare(ulong p)
103 : {
104 28 : long k = 2;
105 7 : for (;; k++)
106 7 : {
107 35 : long i = krouu(k, p);
108 35 : if (!i) pari_err_PRIME("Fl_nonsquare",utoipos(p));
109 63 : if (i < 0) return k;
110 : }
111 : }
112 :
113 : static GEN
114 8881 : Flx_root_mod_2(GEN f)
115 : {
116 8881 : int z1, z0 = !(f[2] & 1);
117 : long i,n;
118 : GEN y;
119 :
120 8881 : for (i=2, n=1; i < lg(f); i++) n += f[i];
121 8881 : z1 = n & 1;
122 8881 : y = cgetg(z0+z1+1, t_VECSMALL); i = 1;
123 8881 : if (z0) y[i++] = 0;
124 8881 : if (z1) y[i ] = 1;
125 8881 : return y;
126 : }
127 : static ulong
128 56 : Flx_oneroot_mod_2(GEN f)
129 : {
130 : long i,n;
131 56 : if (!(f[2] & 1)) return 0;
132 56 : for (i=2, n=1; i < lg(f); i++) n += f[i];
133 56 : if (n & 1) return 1;
134 28 : return 2;
135 : }
136 :
137 : static GEN FpX_roots_i(GEN f, GEN p);
138 : static GEN Flx_roots_i(GEN f, ulong p);
139 :
140 : static int
141 3212314 : cmpGuGu(GEN a, GEN b) { return (ulong)a < (ulong)b? -1: (a == b? 0: 1); }
142 :
143 : /* Generic driver to computes the roots of f modulo pp, using 'Roots' when
144 : * pp is a small prime.
145 : * if (gpwrap), check types thoroughly and return t_INTMODs, otherwise
146 : * assume that f is an FpX, pp a prime and return t_INTs */
147 : static GEN
148 70939 : rootmod_aux(GEN f, GEN pp)
149 : {
150 : GEN y;
151 70939 : switch(lg(f))
152 : {
153 14 : case 2: pari_err_ROOTS0("rootmod");
154 49 : case 3: return cgetg(1,t_COL);
155 : }
156 70876 : if (typ(f) == t_VECSMALL)
157 : {
158 67766 : ulong p = pp[2];
159 67766 : if (p == 2)
160 8881 : y = Flx_root_mod_2(f);
161 : else
162 : {
163 58885 : if (!odd(p)) pari_err_PRIME("rootmod",utoi(p));
164 58885 : y = Flx_roots_i(f, p);
165 : }
166 67759 : y = Flc_to_ZC(y);
167 : }
168 : else
169 3110 : y = FpX_roots_i(f, pp);
170 70862 : return y;
171 : }
172 : /* assume that f is a ZX and p a prime */
173 : GEN
174 70939 : FpX_roots(GEN f, GEN p)
175 : {
176 70939 : pari_sp av = avma;
177 70939 : GEN y; ZX_rootmod_init(&f, p); y = rootmod_aux(f, p);
178 70911 : return gerepileupto(av, y);
179 : }
180 :
181 : /* assume x reduced mod p > 2, monic. */
182 : static int
183 21 : FpX_quad_factortype(GEN x, GEN p)
184 : {
185 21 : GEN b = gel(x,3), c = gel(x,2);
186 21 : GEN D = subii(sqri(b), shifti(c,2));
187 21 : return kronecker(D,p);
188 : }
189 : /* assume x reduced mod p, monic. Return one root, or NULL if irreducible */
190 : static GEN
191 7532 : FpX_quad_root(GEN x, GEN p, int unknown)
192 : {
193 7532 : GEN s, D, b = gel(x,3), c = gel(x,2);
194 :
195 7532 : if (absequaliu(p, 2)) {
196 0 : if (!signe(b)) return c;
197 0 : return signe(c)? NULL: gen_1;
198 : }
199 7532 : D = subii(sqri(b), shifti(c,2));
200 7532 : D = remii(D,p);
201 7532 : if (unknown && kronecker(D,p) == -1) return NULL;
202 :
203 7019 : s = Fp_sqrt(D,p);
204 : /* p is not prime, go on and give e.g. maxord a chance to recover */
205 7019 : if (!s) return NULL;
206 7011 : return Fp_halve(Fp_sub(s,b, p), p);
207 : }
208 : static GEN
209 3146 : FpX_otherroot(GEN x, GEN r, GEN p)
210 3146 : { return Fp_neg(Fp_add(gel(x,3), r, p), p); }
211 :
212 : /* disc(x^2+bx+c) = b^2 - 4c */
213 : static ulong
214 21200559 : Fl_disc_bc(ulong b, ulong c, ulong p)
215 21200559 : { return Fl_sub(Fl_sqr(b,p), Fl_double(Fl_double(c,p),p), p); }
216 : /* p > 2 */
217 : static ulong
218 20466578 : Flx_quad_root(GEN x, ulong p, int unknown)
219 : {
220 20466578 : ulong s, b = x[3], c = x[2];
221 20466578 : ulong D = Fl_disc_bc(b, c, p);
222 20462537 : if (unknown && krouu(D,p) == -1) return p;
223 13692784 : s = Fl_sqrt(D,p);
224 13686700 : if (s==~0UL) return p;
225 13686687 : return Fl_halve(Fl_sub(s,b, p), p);
226 : }
227 : static ulong
228 12262873 : Flx_otherroot(GEN x, ulong r, ulong p)
229 12262873 : { return Fl_neg(Fl_add(x[3], r, p), p); }
230 :
231 :
232 : /* 'todo' contains the list of factors to be split.
233 : * 'done' the list of finished factors, no longer touched */
234 : struct split_t { GEN todo, done; };
235 : static void
236 4930559 : split_init(struct split_t *S, long max)
237 : {
238 4930559 : S->todo = vectrunc_init(max);
239 4930555 : S->done = vectrunc_init(max);
240 4930409 : }
241 : #if 0
242 : /* move todo[i] to done */
243 : static void
244 : split_convert(struct split_t *S, long i)
245 : {
246 : long n = lg(S->todo)-1;
247 : vectrunc_append(S->done, gel(S->todo,i));
248 : if (n) gel(S->todo,i) = gel(S->todo, n);
249 : setlg(S->todo, n);
250 : }
251 : #endif
252 : /* append t to todo */
253 : static void
254 5177019 : split_add(struct split_t *S, GEN t) { vectrunc_append(S->todo, t); }
255 : /* delete todo[i], add t to done */
256 : static void
257 5177052 : split_moveto_done(struct split_t *S, long i, GEN t)
258 : {
259 5177052 : long n = lg(S->todo)-1;
260 5177052 : vectrunc_append(S->done, t);
261 5177023 : if (n) gel(S->todo,i) = gel(S->todo, n);
262 5177023 : setlg(S->todo, n);
263 :
264 5177561 : }
265 : /* append t to done */
266 : static void
267 391697 : split_add_done(struct split_t *S, GEN t)
268 391697 : { vectrunc_append(S->done, t); }
269 : /* split todo[i] into a and b */
270 : static void
271 337162 : split_todo(struct split_t *S, long i, GEN a, GEN b)
272 : {
273 337162 : gel(S->todo, i) = a;
274 337162 : split_add(S, b);
275 337161 : }
276 : /* split todo[i] into a and b, moved to done */
277 : static void
278 372319 : split_done(struct split_t *S, long i, GEN a, GEN b)
279 : {
280 372319 : split_moveto_done(S, i, a);
281 372321 : split_add_done(S, b);
282 372322 : }
283 :
284 : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
285 : static GEN
286 3110 : FpX_roots_i(GEN f, GEN p)
287 : {
288 : GEN pol, pol0, a, q;
289 : struct split_t S;
290 :
291 3110 : split_init(&S, lg(f)-1);
292 3110 : settyp(S.done, t_COL);
293 3110 : if (ZX_valrem(f, &f)) split_add_done(&S, gen_0);
294 3110 : switch(degpol(f))
295 : {
296 7 : case 0: return ZC_copy(S.done);
297 14 : case 1: split_add_done(&S, subii(p, gel(f,2))); return ZC_copy(S.done);
298 : case 2: {
299 1743 : GEN s, r = FpX_quad_root(f, p, 1);
300 1743 : if (r) {
301 1743 : split_add_done(&S, r);
302 1743 : s = FpX_otherroot(f,r, p);
303 : /* f not known to be square free yet */
304 1743 : if (!equalii(r, s)) split_add_done(&S, s);
305 : }
306 1743 : return sort(S.done);
307 : }
308 : }
309 :
310 1346 : a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
311 1346 : if (lg(a) < 3) pari_err_PRIME("rootmod",p);
312 1346 : a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
313 1346 : a = FpX_gcd(f,a, p);
314 1346 : if (!degpol(a)) return ZC_copy(S.done);
315 1346 : split_add(&S, FpX_normalize(a,p));
316 :
317 1346 : q = shifti(p,-1);
318 1346 : pol0 = icopy(gen_1); /* constant term, will vary in place */
319 1346 : pol = deg1pol_shallow(gen_1, pol0, varn(f));
320 2819 : for (pol0[2] = 1;; pol0[2]++)
321 1473 : {
322 2819 : long j, l = lg(S.todo);
323 2819 : if (l == 1) return sort(S.done);
324 1480 : if (pol0[2] == 100 && !BPSW_psp(p)) pari_err_PRIME("polrootsmod",p);
325 3074 : for (j = 1; j < l; j++)
326 : {
327 1601 : GEN b, r, s, c = gel(S.todo,j);
328 1601 : switch(degpol(c))
329 : { /* convert linear and quadratics to roots, try to split the rest */
330 : case 1:
331 195 : split_moveto_done(&S, j, subii(p, gel(c,2)));
332 195 : j--; l--; break;
333 : case 2:
334 1272 : r = FpX_quad_root(c, p, 0);
335 1272 : if (!r) pari_err_PRIME("polrootsmod",p);
336 1265 : s = FpX_otherroot(c,r, p);
337 1265 : split_done(&S, j, r, s);
338 1265 : j--; l--; break;
339 : default:
340 134 : b = FpXQ_pow(pol,q, c,p);
341 134 : if (degpol(b) <= 0) continue;
342 121 : b = FpX_gcd(c,FpX_Fp_sub_shallow(b,gen_1,p), p);
343 121 : if (!degpol(b)) continue;
344 121 : b = FpX_normalize(b, p);
345 121 : c = FpX_div(c,b, p);
346 121 : split_todo(&S, j, b, c);
347 : }
348 : }
349 : }
350 : }
351 :
352 : /* Assume f is normalized */
353 : static ulong
354 174315 : Flx_cubic_root(GEN ff, ulong p)
355 : {
356 174315 : GEN f = Flx_normalize(ff,p);
357 174316 : ulong pi = get_Fl_red(p);
358 174317 : ulong a = f[4], b=f[3], c=f[2], p3 = p%3==1 ? (2*p+1)/3 :(p+1)/3;
359 174317 : ulong t = Fl_mul_pre(a, p3, p, pi), t2 = Fl_sqr_pre(t, p, pi);
360 174318 : ulong A = Fl_sub(b, Fl_triple(t2, p), p);
361 174316 : ulong B = Fl_addmul_pre(c, t, Fl_sub(Fl_double(t2, p), b, p), p, pi);
362 174318 : ulong A3 = Fl_mul_pre(A, p3, p, pi);
363 174318 : ulong A32 = Fl_sqr_pre(A3, p, pi), A33 = Fl_mul_pre(A3, A32, p, pi);
364 174317 : ulong S = Fl_neg(B,p), P = Fl_neg(A3,p);
365 174317 : ulong D = Fl_add(Fl_sqr_pre(S, p, pi), Fl_double(Fl_double(A33, p), p), p);
366 174315 : ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
367 174317 : if (s!=~0UL)
368 : {
369 109726 : ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
370 109726 : if (p%3==2) /* 1 solutions */
371 19356 : vS1 = Fl_powu_pre(S1, (2*p-1)/3, p, pi);
372 : else
373 : {
374 90370 : vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
375 90370 : if (vS1==~0UL) return p; /*0 solutions*/
376 : /*3 solutions*/
377 : }
378 79212 : vS2 = P? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): 0;
379 79211 : return Fl_sub(Fl_add(vS1,vS2, p), t, p);
380 : }
381 : else
382 : {
383 64591 : pari_sp av = avma;
384 64591 : GEN S1 = mkvecsmall2(Fl_halve(S, p), Fl_halve(1UL, p));
385 64591 : GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
386 : ulong Sa;
387 64592 : if (!vS1) return p; /*0 solutions, p%3==2*/
388 64592 : Sa = vS1[1];
389 64592 : if (p%3==1) /*1 solutions*/
390 : {
391 24086 : ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
392 24086 : if (Fa!=P)
393 16026 : Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
394 : }
395 64592 : avma = av;
396 64592 : return Fl_sub(Fl_double(Sa,p),t,p);
397 : }
398 : }
399 :
400 : /* assume p > 2 prime */
401 : static ulong
402 3133256 : Flx_oneroot_i(GEN f, ulong p, long fl)
403 : {
404 : GEN pol, a;
405 : ulong q;
406 : long da;
407 :
408 3133256 : if (Flx_val(f)) return 0;
409 3132534 : switch(degpol(f))
410 : {
411 12188 : case 1: return Fl_neg(f[2], p);
412 2683232 : case 2: return Flx_quad_root(f, p, 1);
413 158508 : case 3: if (p>3) return Flx_cubic_root(f, p); /*FALL THROUGH*/
414 : }
415 :
416 278605 : if (!fl)
417 : {
418 247223 : a = Flxq_powu(polx_Flx(f[1]), p - 1, f,p);
419 247223 : if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
420 247223 : a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
421 247223 : a = Flx_gcd(f,a, p);
422 31382 : } else a = f;
423 278605 : da = degpol(a);
424 278598 : if (!da) return p;
425 199088 : a = Flx_normalize(a,p);
426 :
427 199099 : q = p >> 1;
428 199099 : pol = polx_Flx(f[1]);
429 299221 : for(pol[2] = 1;; pol[2]++)
430 : {
431 399345 : if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
432 299233 : switch(da)
433 : {
434 124300 : case 1: return Fl_neg(a[2], p);
435 59023 : case 2: return Flx_quad_root(a, p, 0);
436 15814 : case 3: if (p>3) return Flx_cubic_root(a, p); /*FALL THROUGH*/
437 : default: {
438 100096 : GEN b = Flxq_powu(pol,q, a,p);
439 : long db;
440 100120 : if (degpol(b) <= 0) continue;
441 95038 : b = Flx_gcd(a,Flx_Fl_add(b,p-1,p), p);
442 95052 : db = degpol(b); if (!db) continue;
443 95045 : b = Flx_normalize(b, p);
444 95053 : if (db <= (da >> 1)) {
445 57911 : a = b;
446 57911 : da = db;
447 : } else {
448 37142 : a = Flx_div(a,b, p);
449 37139 : da -= db;
450 : }
451 : }
452 : }
453 : }
454 : }
455 :
456 : /* assume p > 2 prime */
457 : static GEN
458 4848 : FpX_oneroot_i(GEN f, GEN p)
459 : {
460 : GEN pol, pol0, a, q;
461 : long da;
462 :
463 4848 : if (ZX_val(f)) return gen_0;
464 4582 : switch(degpol(f))
465 : {
466 675 : case 1: return subii(p, gel(f,2));
467 3837 : case 2: return FpX_quad_root(f, p, 1);
468 : }
469 :
470 70 : a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
471 70 : if (lg(a) < 3) pari_err_PRIME("rootmod",p);
472 70 : a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
473 70 : a = FpX_gcd(f,a, p);
474 70 : da = degpol(a);
475 70 : if (!da) return NULL;
476 70 : a = FpX_normalize(a,p);
477 :
478 70 : q = shifti(p,-1);
479 70 : pol0 = icopy(gen_1); /* constant term, will vary in place */
480 70 : pol = deg1pol_shallow(gen_1, pol0, varn(f));
481 224 : for (pol0[2]=1; ; pol0[2]++)
482 : {
483 378 : if (pol0[2] == 1000 && !BPSW_psp(p)) pari_err_PRIME("FpX_oneroot",p);
484 224 : switch(da)
485 : {
486 42 : case 1: return subii(p, gel(a,2));
487 28 : case 2: return FpX_quad_root(a, p, 0);
488 : default: {
489 154 : GEN b = FpXQ_pow(pol,q, a,p);
490 : long db;
491 154 : if (degpol(b) <= 0) continue;
492 147 : b = FpX_gcd(a,FpX_Fp_sub_shallow(b,gen_1,p), p);
493 147 : db = degpol(b); if (!db) continue;
494 147 : b = FpX_normalize(b, p);
495 147 : if (db <= (da >> 1)) {
496 105 : a = b;
497 105 : da = db;
498 : } else {
499 42 : a = FpX_div(a,b, p);
500 42 : da -= db;
501 : }
502 : }
503 : }
504 : }
505 : }
506 :
507 : ulong
508 3093602 : Flx_oneroot(GEN f, ulong p)
509 : {
510 3093602 : pari_sp av = avma;
511 : ulong r;
512 3093602 : switch(lg(f))
513 : {
514 0 : case 2: return 0;
515 0 : case 3: avma = av; return p;
516 : }
517 3093602 : if (p == 2) return Flx_oneroot_mod_2(f);
518 3093602 : r = Flx_oneroot_i(Flx_normalize(f, p), p, 0);
519 3093602 : avma = av; return r;
520 : }
521 :
522 : ulong
523 32245 : Flx_oneroot_split(GEN f, ulong p)
524 : {
525 32245 : pari_sp av = avma;
526 : ulong r;
527 32245 : switch(lg(f))
528 : {
529 0 : case 2: return 0;
530 0 : case 3: avma = av; return p;
531 : }
532 32245 : if (p == 2) return Flx_oneroot_mod_2(f);
533 32245 : r = Flx_oneroot_i(Flx_normalize(f, p), p, 1);
534 32283 : avma = av; return r;
535 : }
536 :
537 : /* assume that p is prime */
538 : GEN
539 12313 : FpX_oneroot(GEN f, GEN pp) {
540 12313 : pari_sp av = avma;
541 12313 : ZX_rootmod_init(&f, pp);
542 12313 : switch(lg(f))
543 : {
544 0 : case 2: avma = av; return gen_0;
545 0 : case 3: avma = av; return NULL;
546 : }
547 12313 : if (typ(f) == t_VECSMALL)
548 : {
549 7465 : ulong r, p = pp[2];
550 7465 : if (p == 2)
551 56 : r = Flx_oneroot_mod_2(f);
552 : else
553 7409 : r = Flx_oneroot_i(f, p, 0);
554 7465 : avma = av;
555 7465 : return (r == p)? NULL: utoi(r);
556 : }
557 4848 : f = FpX_oneroot_i(f, pp);
558 4848 : if (!f) { avma = av; return NULL; }
559 4848 : return gerepileuptoint(av, f);
560 : }
561 :
562 : /* returns a root of unity in F_p that is suitable for finding a factor */
563 : /* of degree deg_factor of a polynomial of degree deg; the order is */
564 : /* returned in n */
565 : /* A good choice seems to be n close to deg/deg_factor; we choose n */
566 : /* twice as big and decrement until it divides p-1. */
567 : static GEN
568 217 : good_root_of_unity(GEN p, long deg, long deg_factor, long *pt_n)
569 : {
570 217 : pari_sp ltop = avma;
571 : GEN pm, factn, power, base, zeta;
572 : long n;
573 :
574 217 : pm = subis (p, 1ul);
575 217 : for (n = deg / 2 / deg_factor + 1; !dvdiu (pm, n); n--);
576 217 : factn = Z_factor(stoi(n));
577 217 : power = diviuexact (pm, n);
578 217 : base = gen_1;
579 : do {
580 378 : base = addis (base, 1l);
581 378 : zeta = Fp_pow (base, power, p);
582 : }
583 378 : while (!equaliu (Fp_order (zeta, factn, p), n));
584 217 : *pt_n = n;
585 217 : return gerepileuptoint (ltop, zeta);
586 : }
587 :
588 : GEN
589 924 : FpX_oneroot_split(GEN fact, GEN p)
590 : {
591 924 : pari_sp av = avma;
592 : long n, deg_f, i, dmin;
593 : GEN prim, expo, minfactor, xplusa, zeta, xpow;
594 924 : fact = FpX_normalize(fact, p);
595 924 : deg_f = degpol(fact);
596 924 : if (deg_f<=2) return FpX_oneroot(fact, p);
597 217 : minfactor = fact; /* factor of minimal degree found so far */
598 217 : dmin = degpol(minfactor);
599 217 : prim = good_root_of_unity(p, deg_f, 1, &n);
600 217 : expo = diviuexact(subiu(p, 1), n);
601 217 : xplusa = pol_x(varn(fact));
602 217 : zeta = gen_1;
603 1071 : while (dmin != 1)
604 : {
605 : /* split minfactor by computing its gcd with (X+a)^exp-zeta, where */
606 : /* zeta varies over the roots of unity in F_p */
607 637 : fact = minfactor; deg_f = dmin;
608 : /* update X+a, avoid a=0 */
609 637 : gel (xplusa, 2) = addis (gel (xplusa, 2), 1);
610 637 : xpow = FpXQ_pow (xplusa, expo, fact, p);
611 1260 : for (i = 0; i < n; i++)
612 : {
613 980 : GEN tmp = FpX_gcd(FpX_Fp_sub(xpow, zeta, p), fact, p);
614 980 : long dtmp = degpol(tmp);
615 980 : if (dtmp > 0 && dtmp < deg_f)
616 : {
617 427 : fact = FpX_div(fact, tmp, p); deg_f = degpol(fact);
618 427 : if (dtmp < dmin)
619 : {
620 427 : minfactor = FpX_normalize (tmp, p);
621 427 : dmin = dtmp;
622 427 : if (dmin == 1 || dmin <= deg_f / (n / 2) + 1)
623 : /* stop early to avoid too many gcds */
624 : break;
625 : }
626 : }
627 623 : zeta = Fp_mul (zeta, prim, p);
628 : }
629 : }
630 217 : return gerepileuptoint(av, Fp_neg(gel(minfactor,2), p));
631 : }
632 :
633 : /*******************************************************************/
634 : /* */
635 : /* FACTORISATION MODULO p */
636 : /* */
637 : /*******************************************************************/
638 :
639 : /* F / E a vector of vectors of factors / exponents of virtual length l
640 : * (their real lg may be larger). Set their lg to j, concat and return [F,E] */
641 : static GEN
642 482087 : FE_concat(GEN F, GEN E, long l)
643 : {
644 482087 : setlg(E,l); E = shallowconcat1(E);
645 482087 : setlg(F,l); F = shallowconcat1(F); return mkvec2(F,E);
646 : }
647 :
648 : static GEN
649 14 : ddf_to_ddf2_i(GEN V, long fl)
650 : {
651 : GEN F, D;
652 14 : long i, j, l = lg(V);
653 14 : F = cgetg(l, t_VEC);
654 14 : D = cgetg(l, t_VECSMALL);
655 112 : for (i = j = 1; i < l; i++)
656 : {
657 98 : GEN Vi = gel(V,i);
658 98 : if ((fl==2 && F2x_degree(Vi) == 0)
659 84 : ||(fl==0 && degpol(Vi) == 0)) continue;
660 35 : gel(F,j) = Vi;
661 35 : uel(D,j) = i; j++;
662 : }
663 14 : setlg(F,j);
664 14 : setlg(D,j); return mkvec2(F,D);
665 : }
666 :
667 : GEN
668 7 : ddf_to_ddf2(GEN V)
669 7 : { return ddf_to_ddf2_i(V, 0); }
670 :
671 : static GEN
672 7 : F2x_ddf_to_ddf2(GEN V)
673 7 : { return ddf_to_ddf2_i(V, 2); }
674 :
675 : GEN
676 354779 : vddf_to_simplefact(GEN V, long d)
677 : {
678 : GEN E, F;
679 354779 : long i, j, c, l = lg(V);
680 354779 : F = cgetg(d+1, t_VECSMALL);
681 354779 : E = cgetg(d+1, t_VECSMALL);
682 712611 : for (i = c = 1; i < l; i++)
683 : {
684 357832 : GEN Vi = gel(V,i);
685 357832 : long l = lg(Vi);
686 2463663 : for (j = 1; j < l; j++)
687 : {
688 2105831 : long k, n = degpol(gel(Vi,j)) / j;
689 2105831 : for (k = 1; k <= n; k++) { uel(F,c) = j; uel(E,c) = i; c++; }
690 : }
691 : }
692 354779 : setlg(F,c);
693 354779 : setlg(E,c);
694 354779 : return sort_factor(mkvec2(F,E), (void*)&cmpGuGu, cmp_nodata);
695 : }
696 :
697 : /* product of terms of degree 1 in factorization of f */
698 : GEN
699 142823 : FpX_split_part(GEN f, GEN p)
700 : {
701 142823 : long n = degpol(f);
702 142823 : GEN z, X = pol_x(varn(f));
703 142823 : if (n <= 1) return f;
704 142508 : f = FpX_red(f, p);
705 142508 : z = FpX_sub(FpX_Frobenius(f, p), X, p);
706 142508 : return FpX_gcd(z,f,p);
707 : }
708 :
709 : /* Compute the number of roots in Fp without counting multiplicity
710 : * return -1 for 0 polynomial. lc(f) must be prime to p. */
711 : long
712 99148 : FpX_nbroots(GEN f, GEN p)
713 : {
714 99148 : pari_sp av = avma;
715 99148 : GEN z = FpX_split_part(f, p);
716 99148 : avma = av; return degpol(z);
717 : }
718 :
719 : int
720 0 : FpX_is_totally_split(GEN f, GEN p)
721 : {
722 0 : long n=degpol(f);
723 0 : pari_sp av = avma;
724 0 : if (n <= 1) return 1;
725 0 : if (abscmpui(n, p) > 0) return 0;
726 0 : f = FpX_red(f, p);
727 0 : avma = av; return gequalX(FpX_Frobenius(f, p));
728 : }
729 :
730 : long
731 2406352 : Flx_nbroots(GEN f, ulong p)
732 : {
733 2406352 : long n = degpol(f);
734 2406352 : pari_sp av = avma;
735 : GEN z;
736 2406352 : if (n <= 1) return n;
737 2404147 : if (n == 2)
738 : {
739 : ulong D;
740 11844 : if (p==2) return (f[2]==0) + (f[2]!=f[3]);
741 10962 : D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
742 10962 : return 1 + krouu(D,p);
743 : }
744 2392303 : z = Flx_sub(Flx_Frobenius(f, p), polx_Flx(f[1]), p);
745 2392303 : z = Flx_gcd(z, f, p);
746 2392303 : avma = av; return degpol(z);
747 : }
748 :
749 : long
750 4389 : FpX_ddf_degree(GEN T, GEN XP, GEN p)
751 : {
752 4389 : pari_sp av = avma;
753 : GEN X, b, g, xq;
754 : long i, j, n, v, B, l, m;
755 : pari_timer ti;
756 : hashtable h;
757 :
758 4389 : n = get_FpX_degree(T); v = get_FpX_var(T);
759 4389 : X = pol_x(v);
760 4389 : if (ZX_equal(X,XP)) return 1;
761 4389 : B = n/2;
762 4389 : l = usqrt(B);
763 4389 : m = (B+l-1)/l;
764 4389 : T = FpX_get_red(T, p);
765 4389 : hash_init_GEN(&h, l+2, ZX_equal, 1);
766 4389 : hash_insert_long(&h, X, 0);
767 4389 : hash_insert_long(&h, XP, 1);
768 4389 : if (DEBUGLEVEL>=7) timer_start(&ti);
769 4389 : b = XP;
770 4389 : xq = FpXQ_powers(b, brent_kung_optpow(n, l-1, 1), T, p);
771 4389 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq baby");
772 10367 : for (i = 3; i <= l+1; i++)
773 : {
774 6776 : b = FpX_FpXQV_eval(b, xq, T, p);
775 6776 : if (gequalX(b)) { avma = av; return i-1; }
776 5978 : hash_insert_long(&h, b, i-1);
777 : }
778 3591 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: baby");
779 3591 : g = b;
780 3591 : xq = FpXQ_powers(g, brent_kung_optpow(n, m, 1), T, p);
781 3591 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq giant");
782 12278 : for(i = 2; i <= m+1; i++)
783 : {
784 10696 : g = FpX_FpXQV_eval(g, xq, T, p);
785 10696 : if (hash_haskey_long(&h, g, &j)) { avma=av; return l*i-j; }
786 : }
787 1582 : avma = av; return n;
788 : }
789 :
790 : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
791 : static GEN
792 669 : FpX_ddf_Shoup(GEN T, GEN XP, GEN p)
793 : {
794 : GEN b, g, h, F, f, Tr, xq;
795 : long i, j, n, v, B, l, m;
796 : pari_timer ti;
797 :
798 669 : n = get_FpX_degree(T); v = get_FpX_var(T);
799 669 : if (n == 0) return cgetg(1, t_VEC);
800 669 : if (n == 1) return mkvec(get_FpX_mod(T));
801 623 : B = n/2;
802 623 : l = usqrt(B);
803 623 : m = (B+l-1)/l;
804 623 : T = FpX_get_red(T, p);
805 623 : b = cgetg(l+2, t_VEC);
806 623 : gel(b, 1) = pol_x(v);
807 623 : gel(b, 2) = XP;
808 623 : if (DEBUGLEVEL>=7) timer_start(&ti);
809 623 : xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1), T, p);
810 623 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq baby");
811 811 : for (i = 3; i <= l+1; i++)
812 188 : gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
813 623 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: baby");
814 623 : xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1), T, p);
815 623 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq giant");
816 623 : g = cgetg(m+1, t_VEC);
817 623 : gel(g, 1) = gel(xq, 2);
818 623 : for(i = 2; i <= m; i++) gel(g, i) = FpX_FpXQV_eval(gel(g, i-1), xq, T, p);
819 623 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: giant");
820 623 : h = cgetg(m+1, t_VEC);
821 1945 : for (j = 1; j <= m; j++)
822 : {
823 1322 : pari_sp av = avma;
824 1322 : GEN gj = gel(g,j), e = FpX_sub(gj, gel(b,1), p);
825 1322 : for (i = 2; i <= l; i++) e = FpXQ_mul(e, FpX_sub(gj, gel(b,i), p), T, p);
826 1322 : gel(h,j) = gerepileupto(av, e);
827 : }
828 623 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: diff");
829 623 : Tr = get_FpX_mod(T);
830 623 : F = cgetg(m+1, t_VEC);
831 1945 : for (j = 1; j <= m; j++)
832 : {
833 1322 : GEN u = FpX_gcd(Tr, gel(h,j), p);
834 1322 : if (degpol(u))
835 : {
836 226 : u = FpX_normalize(u, p);
837 226 : Tr = FpX_div(Tr, u, p);
838 : }
839 1322 : gel(F,j) = u;
840 : }
841 623 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: F");
842 623 : f = const_vec(n, pol_1(v));
843 1945 : for (j = 1; j <= m; j++)
844 : {
845 1322 : GEN e = gel(F, j);
846 1357 : for (i=l-1; i >= 0; i--)
847 : {
848 1357 : GEN u = FpX_gcd(e, FpX_sub(gel(g, j), gel(b, i+1), p), p);
849 1357 : if (degpol(u))
850 : {
851 233 : u = FpX_normalize(u, p);
852 233 : gel(f, l*j-i) = u;
853 233 : e = FpX_div(e, u, p);
854 : }
855 1357 : if (!degpol(e)) break;
856 : }
857 : }
858 623 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: f");
859 623 : if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
860 623 : return f;
861 : }
862 :
863 : static void
864 0 : FpX_edf_simple(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
865 : {
866 0 : long n = degpol(Tp), r = n/d, ct = 0;
867 : GEN T, f, ff, p2;
868 0 : if (r==1) { gel(V, idx) = Tp; return; }
869 0 : p2 = shifti(p,-1);
870 0 : T = FpX_get_red(Tp, p);
871 0 : XP = FpX_rem(XP, T, p);
872 : while (1)
873 0 : {
874 0 : pari_sp btop = avma;
875 : long i;
876 0 : GEN g = random_FpX(n, varn(Tp), p);
877 0 : GEN t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
878 0 : if (signe(t) == 0) continue;
879 0 : for(i=1; i<=10; i++)
880 : {
881 0 : pari_sp btop2 = avma;
882 0 : GEN R = FpXQ_pow(FpX_Fp_add(t, randomi(p), p), p2, T, p);
883 0 : f = FpX_gcd(FpX_Fp_sub(R, gen_1, p), Tp, p);
884 0 : if (degpol(f) > 0 && degpol(f) < n) break;
885 0 : avma = btop2;
886 : }
887 0 : if (degpol(f) > 0 && degpol(f) < n) break;
888 0 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_simple",p);
889 0 : avma = btop;
890 : }
891 0 : f = FpX_normalize(f, p);
892 0 : ff = FpX_div(Tp, f ,p);
893 0 : FpX_edf_simple(f, XP, d, p, V, idx);
894 0 : FpX_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
895 : }
896 :
897 : static void
898 462 : FpX_edf_rec(GEN T, GEN hp, GEN t, long d, GEN p2, GEN p, GEN V, long idx)
899 : {
900 : pari_sp av;
901 462 : GEN Tp = get_FpX_mod(T);
902 462 : long n = degpol(hp), vT = varn(Tp), ct = 0;
903 : GEN u1, u2, f1, f2, R, h;
904 462 : h = FpX_get_red(hp, p);
905 462 : t = FpX_rem(t, T, p);
906 462 : av = avma;
907 : do
908 : {
909 761 : avma = av;
910 761 : R = FpXQ_pow(deg1pol(gen_1, randomi(p), vT), p2, h, p);
911 761 : u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
912 761 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_rec",p);
913 761 : } while (degpol(u1)==0 || degpol(u1)==n);
914 462 : f1 = FpX_gcd(FpX_FpXQ_eval(u1, t, T, p), Tp, p);
915 462 : f1 = FpX_normalize(f1, p);
916 462 : u2 = FpX_div(hp, u1, p);
917 462 : f2 = FpX_div(Tp, f1, p);
918 462 : if (degpol(u1)==1)
919 336 : gel(V, idx) = f1;
920 : else
921 126 : FpX_edf_rec(FpX_get_red(f1, p), u1, t, d, p2, p, V, idx);
922 462 : idx += degpol(f1)/d;
923 462 : if (degpol(u2)==1)
924 350 : gel(V, idx) = f2;
925 : else
926 112 : FpX_edf_rec(FpX_get_red(f2, p), u2, t, d, p2, p, V, idx);
927 462 : }
928 :
929 : /* assume Tp a squarefree product of r > 1 irred. factors of degree d */
930 : static void
931 224 : FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
932 : {
933 224 : long n = degpol(Tp), r = n/d, vT = varn(Tp), ct = 0;
934 : GEN T, h, t;
935 : pari_timer ti;
936 :
937 224 : T = FpX_get_red(Tp, p);
938 224 : XP = FpX_rem(XP, T, p);
939 224 : if (DEBUGLEVEL>=7) timer_start(&ti);
940 : do
941 : {
942 224 : GEN g = random_FpX(n, vT, p);
943 224 : t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
944 224 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
945 224 : h = FpXQ_minpoly(t, T, p);
946 224 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
947 224 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf",p);
948 224 : } while (degpol(h) != r);
949 224 : FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
950 224 : }
951 :
952 : static GEN
953 655 : FpX_factor_Shoup(GEN T, GEN p)
954 : {
955 655 : long i, n, s = 0;
956 : GEN XP, D, V;
957 655 : long e = expi(p);
958 : pari_timer ti;
959 655 : n = get_FpX_degree(T);
960 655 : T = FpX_get_red(T, p);
961 655 : if (DEBUGLEVEL>=6) timer_start(&ti);
962 655 : XP = FpX_Frobenius(T, p);
963 655 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
964 655 : D = FpX_ddf_Shoup(T, XP, p);
965 655 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
966 655 : s = ddf_to_nbfact(D);
967 655 : V = cgetg(s+1, t_COL);
968 5507 : for (i = 1, s = 1; i <= n; i++)
969 : {
970 4852 : GEN Di = gel(D,i);
971 4852 : long ni = degpol(Di), ri = ni/i;
972 4852 : if (ni == 0) continue;
973 664 : Di = FpX_normalize(Di, p);
974 664 : if (ni == i) { gel(V, s++) = Di; continue; }
975 224 : if (ri <= e*expu(e))
976 224 : FpX_edf(Di, XP, i, p, V, s);
977 : else
978 0 : FpX_edf_simple(Di, XP, i, p, V, s);
979 224 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
980 224 : s += ri;
981 : }
982 655 : return V;
983 : }
984 :
985 : long
986 536358 : ddf_to_nbfact(GEN D)
987 : {
988 536358 : long l = lg(D), i, s = 0;
989 536358 : for(i = 1; i < l; i++) s += degpol(gel(D,i))/i;
990 536358 : return s;
991 : }
992 :
993 : /* Yun algorithm: Assume p > degpol(T) */
994 : static GEN
995 672 : FpX_factor_Yun(GEN T, GEN p)
996 : {
997 672 : long n = degpol(T), i = 1;
998 672 : GEN a, b, c, d = FpX_deriv(T, p);
999 672 : GEN V = cgetg(n+1,t_VEC);
1000 672 : a = FpX_gcd(T, d, p);
1001 672 : if (degpol(a) == 0) return mkvec(T);
1002 488 : b = FpX_div(T, a, p);
1003 : do
1004 : {
1005 2432 : c = FpX_div(d, a, p);
1006 2432 : d = FpX_sub(c, FpX_deriv(b, p), p);
1007 2432 : a = FpX_normalize(FpX_gcd(b, d, p), p);
1008 2432 : gel(V, i++) = a;
1009 2432 : b = FpX_div(b, a, p);
1010 2432 : } while (degpol(b));
1011 488 : setlg(V, i); return V;
1012 : }
1013 : GEN
1014 777 : FpX_factor_squarefree(GEN T, GEN p)
1015 : {
1016 777 : if (lgefint(p)==3)
1017 : {
1018 777 : ulong pp = (ulong)p[2];
1019 777 : GEN u = Flx_factor_squarefree(ZX_to_Flx(T,pp), pp);
1020 777 : return FlxV_to_ZXV(u);
1021 : }
1022 0 : return FpX_factor_Yun(T, p);
1023 : }
1024 :
1025 : long
1026 168 : FpX_ispower(GEN f, ulong k, GEN p, GEN *pt_r)
1027 : {
1028 168 : pari_sp av = avma;
1029 : GEN lc, F;
1030 168 : long i, l, n = degpol(f), v = varn(f);
1031 168 : if (n % k) return 0;
1032 168 : if (lgefint(p)==3)
1033 : {
1034 126 : ulong pp = p[2];
1035 126 : GEN fp = ZX_to_Flx(f, pp);
1036 126 : if (!Flx_ispower(fp, k, pp, pt_r)) { avma = av; return 0; }
1037 105 : if (pt_r) *pt_r = gerepileupto(av, Flx_to_ZX(*pt_r)); else avma = av;
1038 105 : return 1;
1039 : }
1040 42 : lc = Fp_sqrtn(leading_coeff(f), stoi(k), p, NULL);
1041 42 : if (!lc) { av = avma; return 0; }
1042 42 : F = FpX_factor_Yun(f, p); l = lg(F)-1;
1043 1491 : for(i=1; i <= l; i++)
1044 1456 : if (i%k && degpol(gel(F,i))) { avma = av; return 0; }
1045 35 : if (pt_r)
1046 : {
1047 35 : GEN r = scalarpol(lc, v), s = pol_1(v);
1048 1484 : for (i=l; i>=1; i--)
1049 : {
1050 1449 : if (i%k) continue;
1051 294 : s = FpX_mul(s, gel(F,i), p);
1052 294 : r = FpX_mul(r, s, p);
1053 : }
1054 35 : *pt_r = gerepileupto(av, r);
1055 0 : } else av = avma;
1056 35 : return 1;
1057 : }
1058 :
1059 : static GEN
1060 616 : FpX_factor_Cantor(GEN T, GEN p)
1061 : {
1062 616 : GEN E, F, V = FpX_factor_Yun(T, p);
1063 616 : long i, j, l = lg(V);
1064 616 : F = cgetg(l, t_VEC);
1065 616 : E = cgetg(l, t_VEC);
1066 1748 : for (i=1, j=1; i < l; i++)
1067 1132 : if (degpol(gel(V,i)))
1068 : {
1069 655 : GEN Fj = FpX_factor_Shoup(gel(V,i), p);
1070 655 : gel(F, j) = Fj;
1071 655 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1072 655 : j++;
1073 : }
1074 616 : return sort_factor_pol(FE_concat(F,E,j), cmpii);
1075 : }
1076 :
1077 : static GEN
1078 0 : FpX_ddf_i(GEN T, GEN p)
1079 : {
1080 : GEN XP;
1081 0 : T = FpX_get_red(T, p);
1082 0 : XP = FpX_Frobenius(T, p);
1083 0 : return ddf_to_ddf2(FpX_ddf_Shoup(T, XP, p));
1084 : }
1085 :
1086 : GEN
1087 7 : FpX_ddf(GEN f, GEN p)
1088 : {
1089 7 : pari_sp av = avma;
1090 : GEN F;
1091 7 : switch(ZX_factmod_init(&f, p))
1092 : {
1093 7 : case 0: F = F2x_ddf(f);
1094 7 : F2xV_to_ZXV_inplace(gel(F,1)); break;
1095 0 : case 1: F = Flx_ddf(f,p[2]);
1096 0 : FlxV_to_ZXV_inplace(gel(F,1)); break;
1097 0 : default: F = FpX_ddf_i(f,p); break;
1098 : }
1099 7 : return gerepilecopy(av, F);
1100 : }
1101 :
1102 : static GEN Flx_simplefact_Cantor(GEN T, ulong p);
1103 : static GEN
1104 14 : FpX_simplefact_Cantor(GEN T, GEN p)
1105 : {
1106 : GEN V, XP;
1107 : long i, l;
1108 14 : if (lgefint(p) == 3)
1109 : {
1110 0 : ulong pp = p[2];
1111 0 : return Flx_simplefact_Cantor(ZX_to_Flx(T,pp), pp);
1112 : }
1113 14 : T = FpX_get_red(T, p);
1114 14 : XP = FpX_Frobenius(T, p);
1115 14 : V = FpX_factor_Yun(get_FpX_mod(T), p); l = lg(V);
1116 14 : for (i=1; i < l; i++) gel(V,i) = FpX_ddf_Shoup(gel(V,i), XP, p);
1117 14 : return vddf_to_simplefact(V, get_FpX_degree(T));
1118 : }
1119 :
1120 : static int
1121 0 : FpX_isirred_Cantor(GEN Tp, GEN p)
1122 : {
1123 0 : pari_sp av = avma;
1124 : pari_timer ti;
1125 : long n, d;
1126 0 : GEN T = get_FpX_mod(Tp);
1127 0 : GEN dT = FpX_deriv(T, p);
1128 : GEN XP, D;
1129 0 : if (degpol(FpX_gcd(T, dT, p)) != 0) { avma = av; return 0; }
1130 0 : n = get_FpX_degree(T);
1131 0 : T = FpX_get_red(Tp, p);
1132 0 : if (DEBUGLEVEL>=6) timer_start(&ti);
1133 0 : XP = FpX_Frobenius(T, p);
1134 0 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
1135 0 : D = FpX_ddf_Shoup(T, XP, p);
1136 0 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
1137 0 : d = degpol(gel(D, n));
1138 0 : avma = av; return d==n;
1139 : }
1140 :
1141 : static GEN FpX_factor_deg2(GEN f, GEN p, long d, long flag);
1142 :
1143 : /*Assume that p is large and odd*/
1144 : static GEN
1145 1385 : FpX_factor_i(GEN f, GEN pp, long flag)
1146 : {
1147 1385 : long d = degpol(f);
1148 1385 : if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
1149 630 : switch(flag)
1150 : {
1151 616 : default: return FpX_factor_Cantor(f, pp);
1152 14 : case 1: return FpX_simplefact_Cantor(f, pp);
1153 0 : case 2: return FpX_isirred_Cantor(f, pp)? gen_1: NULL;
1154 : }
1155 : }
1156 :
1157 : long
1158 0 : FpX_nbfact_Frobenius(GEN T, GEN XP, GEN p)
1159 : {
1160 0 : pari_sp av = avma;
1161 0 : long s = ddf_to_nbfact(FpX_ddf_Shoup(T, XP, p));
1162 0 : avma = av; return s;
1163 : }
1164 :
1165 : long
1166 0 : FpX_nbfact(GEN T, GEN p)
1167 : {
1168 0 : pari_sp av = avma;
1169 0 : GEN XP = FpX_Frobenius(T, p);
1170 0 : long n = FpX_nbfact_Frobenius(T, XP, p);
1171 0 : avma = av; return n;
1172 : }
1173 :
1174 : /* p > 2 */
1175 : static GEN
1176 7 : FpX_is_irred_2(GEN f, GEN p, long d)
1177 : {
1178 7 : switch(d)
1179 : {
1180 : case -1:
1181 0 : case 0: return NULL;
1182 0 : case 1: return gen_1;
1183 : }
1184 7 : return FpX_quad_factortype(f, p) == -1? gen_1: NULL;
1185 : }
1186 : /* p > 2 */
1187 : static GEN
1188 14 : FpX_degfact_2(GEN f, GEN p, long d)
1189 : {
1190 14 : switch(d)
1191 : {
1192 0 : case -1:retmkvec2(mkvecsmall(-1),mkvecsmall(1));
1193 0 : case 0: return trivial_fact();
1194 0 : case 1: retmkvec2(mkvecsmall(1), mkvecsmall(1));
1195 : }
1196 14 : switch(FpX_quad_factortype(f, p)) {
1197 7 : case 1: retmkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1198 7 : case -1: retmkvec2(mkvecsmall(2), mkvecsmall(1));
1199 0 : default: retmkvec2(mkvecsmall(1), mkvecsmall(2));
1200 : }
1201 : }
1202 :
1203 : GEN
1204 70 : prime_fact(GEN x) { retmkmat2(mkcolcopy(x), mkcol(gen_1)); }
1205 : GEN
1206 67720 : trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
1207 :
1208 : /* not gerepile safe */
1209 : static GEN
1210 734 : FpX_factor_2(GEN f, GEN p, long d)
1211 : {
1212 : GEN r, s, R, S;
1213 : long v;
1214 : int sgn;
1215 734 : switch(d)
1216 : {
1217 7 : case -1: retmkvec2(mkcol(pol_0(varn(f))), mkvecsmall(1));
1218 30 : case 0: retmkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1219 45 : case 1: retmkvec2(mkcol(f), mkvecsmall(1));
1220 : }
1221 652 : r = FpX_quad_root(f, p, 1);
1222 652 : if (!r) return mkvec2(mkcol(f), mkvecsmall(1));
1223 138 : v = varn(f);
1224 138 : s = FpX_otherroot(f, r, p);
1225 138 : if (signe(r)) r = subii(p, r);
1226 138 : if (signe(s)) s = subii(p, s);
1227 138 : sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
1228 138 : R = deg1pol_shallow(gen_1, r, v);
1229 138 : if (!sgn) return mkvec2(mkcol(R), mkvecsmall(2));
1230 47 : S = deg1pol_shallow(gen_1, s, v);
1231 47 : return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
1232 : }
1233 : static GEN
1234 755 : FpX_factor_deg2(GEN f, GEN p, long d, long flag)
1235 : {
1236 755 : switch(flag) {
1237 7 : case 2: return FpX_is_irred_2(f, p, d);
1238 14 : case 1: return FpX_degfact_2(f, p, d);
1239 734 : default: return FpX_factor_2(f, p, d);
1240 : }
1241 : }
1242 :
1243 : static int
1244 58734 : F2x_quad_factortype(GEN x)
1245 58734 : { return x[2] == 7 ? -1: x[2] == 6 ? 1 :0; }
1246 :
1247 : static GEN
1248 7 : F2x_is_irred_2(GEN f, long d)
1249 7 : { return d == 1 || (d==2 && F2x_quad_factortype(f) == -1)? gen_1: NULL; }
1250 :
1251 : static GEN
1252 7140 : F2x_degfact_2(GEN f, long d)
1253 : {
1254 7140 : if (!d) return trivial_fact();
1255 7140 : if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1256 6937 : switch(F2x_quad_factortype(f)) {
1257 2268 : case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1258 2219 : case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1259 2450 : default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1260 : }
1261 : }
1262 :
1263 : static GEN
1264 105616 : F2x_factor_2(GEN f, long d)
1265 : {
1266 105616 : long v = f[1];
1267 105616 : if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1268 97734 : if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1269 50442 : switch(F2x_quad_factortype(f))
1270 : {
1271 12208 : case -1: return mkvec2(mkcol(f), mkvecsmall(1));
1272 23464 : case 0: return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
1273 14770 : default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
1274 : }
1275 : }
1276 : static GEN
1277 112763 : F2x_factor_deg2(GEN f, long d, long flag)
1278 : {
1279 112763 : switch(flag) {
1280 7 : case 2: return F2x_is_irred_2(f, d);
1281 7140 : case 1: return F2x_degfact_2(f, d);
1282 105616 : default: return F2x_factor_2(f, d);
1283 : }
1284 : }
1285 :
1286 : /* xt = NULL or x^(p-1)/2 mod g */
1287 : static void
1288 4809 : split_squares(struct split_t *S, GEN g, ulong p, GEN xt)
1289 : {
1290 4809 : ulong q = p >> 1;
1291 4809 : GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
1292 4809 : long d = degpol(a);
1293 4809 : if (d < 0)
1294 : {
1295 : ulong i;
1296 42 : split_add_done(S, (GEN)1);
1297 42 : for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr(i,p));
1298 : } else {
1299 4767 : if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
1300 4767 : if (d)
1301 : {
1302 4767 : if (xt) xt = Flx_Fl_add(xt, p-1, p); else xt = Flx_Xnm1(g[1], q, p);
1303 4767 : a = Flx_gcd(a, xt, p);
1304 4767 : if (degpol(a)) split_add(S, Flx_normalize(a, p));
1305 : }
1306 : }
1307 4809 : }
1308 : static void
1309 4809 : split_nonsquares(struct split_t *S, GEN g, ulong p, GEN xt)
1310 : {
1311 4809 : ulong q = p >> 1;
1312 4809 : GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
1313 4809 : long d = degpol(a);
1314 4809 : if (d < 0)
1315 : {
1316 28 : ulong i, z = Fl_nonsquare(p);
1317 28 : split_add_done(S, (GEN)z);
1318 28 : for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_mul(z, Fl_sqr(i,p), p));
1319 : } else {
1320 4781 : if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
1321 4781 : if (d)
1322 : {
1323 4781 : if (xt) xt = Flx_Fl_add(xt, 1, p); else xt = Flx_Xn1(g[1], q, p);
1324 4781 : a = Flx_gcd(a, xt, p);
1325 4781 : if (degpol(a)) split_add(S, Flx_normalize(a, p));
1326 : }
1327 : }
1328 4809 : }
1329 : /* p > 2. f monic Flx, f(0) != 0. Add to split_t structs coprime factors
1330 : * of g = \prod_{f(a) = 0} (X - a). Return 0 when f(x) = 0 for all x in Fp* */
1331 : static int
1332 4926904 : split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p)
1333 : {
1334 4926904 : GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
1335 4926690 : long d = degpol(g);
1336 4926603 : if (d < 0) return 0;
1337 4926571 : if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
1338 4926607 : if (!d) return 1;
1339 4910983 : if ((p >> 4) <= (ulong)d)
1340 : { /* small p; split directly using x^((p-1)/2) +/- 1 */
1341 14259 : GEN xt = ((ulong)d < (p>>1))? Flx_rem(monomial_Flx(1, p>>1, g[1]), g, p)
1342 9450 : : NULL;
1343 4809 : split_squares(S, g, p, xt);
1344 4809 : split_nonsquares(S, g, p, xt);
1345 : } else { /* large p; use x^(p-1) - 1 directly */
1346 4906174 : a = Flxq_powu(polx_Flx(f[1]), p-1, g,p);
1347 4903474 : if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
1348 4903474 : a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
1349 4903747 : g = Flx_gcd(g,a, p);
1350 4904510 : if (degpol(g)) split_add(S, Flx_normalize(g,p));
1351 : }
1352 4912230 : return 1;
1353 : }
1354 :
1355 : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
1356 : static GEN
1357 21999947 : Flx_roots_i(GEN f, ulong p)
1358 : {
1359 : GEN pol, g;
1360 21999947 : long v = Flx_valrem(f, &g);
1361 : ulong q;
1362 : struct split_t S;
1363 :
1364 : /* optimization: test for small degree first */
1365 21999750 : switch(degpol(g))
1366 : {
1367 : case 1: {
1368 25515 : ulong r = p - g[2];
1369 25515 : return v? mkvecsmall2(0, r): mkvecsmall(r);
1370 : }
1371 : case 2: {
1372 17040574 : ulong r = Flx_quad_root(g, p, 1), s;
1373 17067628 : if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
1374 11655793 : s = Flx_otherroot(g,r, p);
1375 11703986 : if (r < s)
1376 2921934 : return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
1377 8782052 : else if (r > s)
1378 8781891 : return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
1379 : else
1380 161 : return v? mkvecsmall2(0, s): mkvecsmall(s);
1381 : }
1382 : }
1383 4927525 : q = p >> 1;
1384 4927525 : split_init(&S, lg(f)-1);
1385 4926933 : settyp(S.done, t_VECSMALL);
1386 4926933 : if (v) split_add_done(&S, (GEN)0);
1387 4926933 : if (! split_Flx_cut_out_roots(&S, g, p))
1388 42 : return all_roots_mod_p(p, lg(S.done) == 1);
1389 4927694 : pol = polx_Flx(f[1]);
1390 10186195 : for (pol[2]=1; ; pol[2]++)
1391 5258047 : {
1392 10186195 : long j, l = lg(S.todo);
1393 10186195 : if (l == 1) { vecsmall_sort(S.done); return S.done; }
1394 5257915 : if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
1395 10872105 : for (j = 1; j < l; j++)
1396 : {
1397 5614058 : GEN b, c = gel(S.todo,j);
1398 : ulong r, s;
1399 5614058 : switch(degpol(c))
1400 : {
1401 : case 1:
1402 4804548 : split_moveto_done(&S, j, (GEN)(p - c[2]));
1403 4805068 : j--; l--; break;
1404 : case 2:
1405 371057 : r = Flx_quad_root(c, p, 0);
1406 371062 : if (r == p) pari_err_PRIME("polrootsmod",utoipos(p));
1407 371055 : s = Flx_otherroot(c,r, p);
1408 371055 : split_done(&S, j, (GEN)r, (GEN)s);
1409 371057 : j--; l--; break;
1410 : default:
1411 438427 : b = Flxq_powu(pol,q, c,p); /* pol^(p-1)/2 */
1412 438418 : if (degpol(b) <= 0) continue;
1413 337174 : b = Flx_gcd(c,Flx_Fl_add(b,p-1,p), p);
1414 337177 : if (!degpol(b)) continue;
1415 337043 : b = Flx_normalize(b, p);
1416 337045 : c = Flx_div(c,b, p);
1417 337040 : split_todo(&S, j, b, c);
1418 : }
1419 : }
1420 : }
1421 : }
1422 :
1423 : GEN
1424 21956158 : Flx_roots(GEN f, ulong p)
1425 : {
1426 21956158 : pari_sp av = avma;
1427 21956158 : switch(lg(f))
1428 : {
1429 0 : case 2: pari_err_ROOTS0("Flx_roots");
1430 0 : case 3: avma = av; return cgetg(1, t_VECSMALL);
1431 : }
1432 21957323 : if (p == 2) return Flx_root_mod_2(f);
1433 21957323 : return gerepileuptoleaf(av, Flx_roots_i(Flx_normalize(f, p), p));
1434 : }
1435 :
1436 : /* assume x reduced mod p, monic. */
1437 : static int
1438 734573 : Flx_quad_factortype(GEN x, ulong p)
1439 : {
1440 734573 : ulong b = x[3], c = x[2];
1441 734573 : return krouu(Fl_disc_bc(b, c, p), p);
1442 : }
1443 : static GEN
1444 56 : Flx_is_irred_2(GEN f, ulong p, long d)
1445 : {
1446 56 : if (!d) return NULL;
1447 56 : if (d == 1) return gen_1;
1448 56 : return Flx_quad_factortype(f, p) == -1? gen_1: NULL;
1449 : }
1450 : static GEN
1451 751660 : Flx_degfact_2(GEN f, ulong p, long d)
1452 : {
1453 751660 : if (!d) return trivial_fact();
1454 751660 : if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1455 734517 : switch(Flx_quad_factortype(f, p)) {
1456 349076 : case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1457 376978 : case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1458 8463 : default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1459 : }
1460 : }
1461 : /* p > 2 */
1462 : static GEN
1463 432348 : Flx_factor_2(GEN f, ulong p, long d)
1464 : {
1465 : ulong r, s;
1466 : GEN R,S;
1467 432348 : long v = f[1];
1468 432348 : if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1469 418175 : if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1470 318902 : r = Flx_quad_root(f, p, 1);
1471 318902 : if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
1472 210978 : s = Flx_otherroot(f, r, p);
1473 210978 : r = Fl_neg(r, p);
1474 210978 : s = Fl_neg(s, p);
1475 210978 : if (s < r) lswap(s,r);
1476 210978 : R = mkvecsmall3(v,r,1);
1477 210978 : if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
1478 183208 : S = mkvecsmall3(v,s,1);
1479 183208 : return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
1480 : }
1481 : static GEN
1482 1184064 : Flx_factor_deg2(GEN f, ulong p, long d, long flag)
1483 : {
1484 1184064 : switch(flag) {
1485 56 : case 2: return Flx_is_irred_2(f, p, d);
1486 751660 : case 1: return Flx_degfact_2(f, p, d);
1487 432348 : default: return Flx_factor_2(f, p, d);
1488 : }
1489 : }
1490 :
1491 : void
1492 19957 : F2xV_to_FlxV_inplace(GEN v)
1493 : {
1494 : long i;
1495 19957 : for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_Flx(gel(v,i));
1496 19957 : }
1497 : void
1498 730859 : FlxV_to_ZXV_inplace(GEN v)
1499 : {
1500 : long i;
1501 730859 : for(i=1;i<lg(v);i++) gel(v,i)= Flx_to_ZX(gel(v,i));
1502 730859 : }
1503 : void
1504 144681 : F2xV_to_ZXV_inplace(GEN v)
1505 : {
1506 : long i;
1507 144681 : for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_ZX(gel(v,i));
1508 144681 : }
1509 :
1510 : static GEN
1511 10608 : F2x_Berlekamp_ker(GEN u)
1512 : {
1513 10608 : pari_sp ltop=avma;
1514 10608 : long j,N = F2x_degree(u);
1515 : GEN Q;
1516 : pari_timer T;
1517 10607 : timer_start(&T);
1518 10609 : Q = F2x_matFrobenius(u);
1519 249234 : for (j=1; j<=N; j++)
1520 238625 : F2m_flip(Q,j,j);
1521 10609 : if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
1522 10609 : Q = F2m_ker_sp(Q,0);
1523 10609 : if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
1524 10609 : return gerepileupto(ltop,Q);
1525 : }
1526 : #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
1527 : static long
1528 14567 : F2x_split_Berlekamp(GEN *t)
1529 : {
1530 14567 : GEN u = *t, a, b, vker;
1531 14567 : long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
1532 :
1533 14570 : if (du == 1) return 1;
1534 11212 : if (du == 2)
1535 : {
1536 604 : if (F2x_quad_factortype(u) == 1) /* 0 is a root: shouldn't occur */
1537 : {
1538 0 : t[0] = mkvecsmall2(sv, 2);
1539 0 : t[1] = mkvecsmall2(sv, 3);
1540 0 : return 2;
1541 : }
1542 604 : return 1;
1543 : }
1544 :
1545 10608 : vker = F2x_Berlekamp_ker(u);
1546 10609 : lb = lgcols(vker);
1547 10628 : d = lg(vker)-1;
1548 10628 : ir = 0;
1549 : /* t[i] irreducible for i < ir, still to be treated for i < L */
1550 52986 : for (L=1; L<d; )
1551 : {
1552 : GEN pol;
1553 31751 : if (d == 2)
1554 1744 : pol = F2v_to_F2x(gel(vker,2), sv);
1555 : else
1556 : {
1557 30007 : GEN v = zero_zv(lb);
1558 30014 : v[1] = du;
1559 30014 : v[2] = random_Fl(2); /*Assume vker[1]=1*/
1560 112207 : for (i=2; i<=d; i++)
1561 82192 : if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
1562 30015 : pol = F2v_to_F2x(v, sv);
1563 : }
1564 98654 : for (i=ir; i<L && L<d; i++)
1565 : {
1566 66924 : a = t[i]; la = F2x_degree(a);
1567 66897 : if (la == 1) { set_irred(i); }
1568 66735 : else if (la == 2)
1569 : {
1570 744 : if (F2x_quad_factortype(a) == 1) /* 0 is a root: shouldn't occur */
1571 : {
1572 0 : t[i] = mkvecsmall2(sv, 2);
1573 0 : t[L] = mkvecsmall2(sv, 3); L++;
1574 : }
1575 744 : set_irred(i);
1576 : }
1577 : else
1578 : {
1579 65991 : pari_sp av = avma;
1580 : long lb;
1581 65991 : b = F2x_rem(pol, a);
1582 65977 : if (F2x_degree(b) <= 0) { avma=av; continue; }
1583 21694 : b = F2x_gcd(a,b); lb = F2x_degree(b);
1584 21701 : if (lb && lb < la)
1585 : {
1586 21701 : t[L] = F2x_div(a,b);
1587 21694 : t[i]= b; L++;
1588 : }
1589 0 : else avma = av;
1590 : }
1591 : }
1592 : }
1593 10607 : return d;
1594 : }
1595 : /* assume deg f > 2 */
1596 : static GEN
1597 11699 : F2x_Berlekamp_i(GEN f, long flag)
1598 : {
1599 11699 : long lfact, val, d = F2x_degree(f), j, k, lV;
1600 : GEN y, E, t, V;
1601 :
1602 11700 : val = F2x_valrem(f, &f);
1603 11700 : if (flag == 2 && val) return NULL;
1604 11693 : V = F2x_factor_squarefree(f); lV = lg(V);
1605 11693 : if (flag == 2 && lV > 2) return NULL;
1606 :
1607 : /* to hold factors and exponents */
1608 11623 : t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
1609 11623 : E = cgetg(d+1,t_VECSMALL);
1610 11624 : lfact = 1;
1611 11624 : if (val) {
1612 3482 : if (flag == 1) t[1] = 1; else gel(t,1) = polx_F2x(f[1]);
1613 3482 : E[1] = val; lfact++;
1614 : }
1615 :
1616 49054 : for (k=1; k<lV; k++)
1617 : {
1618 37586 : if (F2x_degree(gel(V, k))==0) continue;
1619 14568 : gel(t,lfact) = gel(V, k);
1620 14568 : d = F2x_split_Berlekamp(&gel(t,lfact));
1621 14568 : if (flag == 2 && d != 1) return NULL;
1622 14414 : if (flag == 1)
1623 1533 : for (j=0; j<d; j++) t[lfact+j] = F2x_degree(gel(t,lfact+j));
1624 14414 : for (j=0; j<d; j++) E[lfact+j] = k;
1625 14414 : lfact += d;
1626 : }
1627 11468 : if (flag == 2) return gen_1; /* irreducible */
1628 11454 : setlg(t, lfact);
1629 11454 : setlg(E, lfact); y = mkvec2(t,E);
1630 : return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
1631 11454 : : sort_factor_pol(y, cmpGuGu);
1632 : }
1633 :
1634 : /* Adapted from Shoup NTL */
1635 : GEN
1636 76228 : F2x_factor_squarefree(GEN f)
1637 : {
1638 : GEN r, t, v, tv;
1639 76228 : long i, q, n = F2x_degree(f);
1640 76229 : GEN u = const_vec(n+1, pol1_F2x(f[1]));
1641 135235 : for(q = 1;;q *= 2)
1642 : {
1643 194241 : r = F2x_gcd(f, F2x_deriv(f));
1644 135228 : if (F2x_degree(r) == 0)
1645 : {
1646 62418 : gel(u, q) = f;
1647 62418 : break;
1648 : }
1649 72814 : t = F2x_div(f, r);
1650 72814 : if (F2x_degree(t) > 0)
1651 : {
1652 : long j;
1653 69267 : for(j = 1;;j++)
1654 : {
1655 107319 : v = F2x_gcd(r, t);
1656 69266 : tv = F2x_div(t, v);
1657 69266 : if (F2x_degree(tv) > 0)
1658 32511 : gel(u, j*q) = tv;
1659 69266 : if (F2x_degree(v) <= 0) break;
1660 38052 : r = F2x_div(r, v);
1661 38052 : t = v;
1662 : }
1663 31214 : if (F2x_degree(r) == 0) break;
1664 : }
1665 59005 : f = F2x_sqrt(r);
1666 : }
1667 560996 : for (i = n; i; i--)
1668 560683 : if (F2x_degree(gel(u,i))) break;
1669 76277 : setlg(u,i+1); return u;
1670 : }
1671 :
1672 : static GEN
1673 79965 : F2x_ddf_simple(GEN T, GEN XP)
1674 : {
1675 79965 : pari_sp av = avma, av2;
1676 : GEN f, z, Tr, X;
1677 79965 : long j, n = F2x_degree(T), v = T[1], B = n/2;
1678 79967 : if (n == 0) return cgetg(1, t_VEC);
1679 79967 : if (n == 1) return mkvec(T);
1680 37482 : z = XP; Tr = T; X = polx_F2x(v);
1681 37484 : f = const_vec(n, pol1_F2x(v));
1682 37484 : av2 = avma;
1683 133515 : for (j = 1; j <= B; j++)
1684 : {
1685 100568 : GEN u = F2x_gcd(Tr, F2x_add(z, X));
1686 100589 : if (F2x_degree(u))
1687 : {
1688 22044 : gel(f, j) = u;
1689 22044 : Tr = F2x_div(Tr, u);
1690 22049 : av2 = avma;
1691 78552 : } else z = gerepileuptoleaf(av2, z);
1692 100615 : if (!F2x_degree(Tr)) break;
1693 96067 : z = F2xq_sqr(z, Tr);
1694 : }
1695 37475 : if (F2x_degree(Tr)) gel(f, F2x_degree(Tr)) = Tr;
1696 37484 : return gerepilecopy(av, f);
1697 : }
1698 :
1699 : GEN
1700 7 : F2x_ddf(GEN T)
1701 : {
1702 : GEN XP;
1703 7 : T = F2x_get_red(T);
1704 7 : XP = F2x_Frobenius(T);
1705 7 : return F2x_ddf_to_ddf2(F2x_ddf_simple(T, XP));
1706 : }
1707 :
1708 : static GEN
1709 6574 : F2xq_frobtrace(GEN a, long d, GEN T)
1710 : {
1711 6574 : pari_sp av = avma;
1712 : long i;
1713 6574 : GEN x = a;
1714 27106 : for(i=1; i<d; i++)
1715 : {
1716 20530 : x = F2x_add(a, F2xq_sqr(x,T));
1717 20532 : if (gc_needed(av, 2))
1718 0 : x = gerepileuptoleaf(av, x);
1719 : }
1720 6576 : return x;
1721 : }
1722 :
1723 : static void
1724 9782 : F2x_edf_simple(GEN Tp, GEN XP, long d, GEN V, long idx)
1725 : {
1726 9782 : long n = F2x_degree(Tp), r = n/d;
1727 : GEN T, f, ff;
1728 9782 : if (r==1) { gel(V, idx) = Tp; return; }
1729 3333 : T = Tp;
1730 3333 : XP = F2x_rem(XP, T);
1731 : while (1)
1732 3243 : {
1733 6576 : pari_sp btop = avma;
1734 : long df;
1735 6576 : GEN g = random_F2x(n, Tp[1]);
1736 6576 : GEN t = F2xq_frobtrace(g, d, T);
1737 6576 : if (lgpol(t) == 0) continue;
1738 4883 : f = F2x_gcd(t, Tp); df = F2x_degree(f);
1739 4883 : if (df > 0 && df < n) break;
1740 1550 : avma = btop;
1741 : }
1742 3333 : ff = F2x_div(Tp, f);
1743 3333 : F2x_edf_simple(f, XP, d, V, idx);
1744 3333 : F2x_edf_simple(ff, XP, d, V, idx+F2x_degree(f)/d);
1745 : }
1746 :
1747 : static GEN
1748 79963 : F2x_factor_Shoup(GEN T)
1749 : {
1750 79963 : long i, n, s = 0;
1751 : GEN XP, D, V;
1752 : pari_timer ti;
1753 79963 : n = F2x_degree(T);
1754 79962 : if (DEBUGLEVEL>=6) timer_start(&ti);
1755 79962 : XP = F2x_Frobenius(T);
1756 79958 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1757 79958 : D = F2x_ddf_simple(T, XP);
1758 79963 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1759 348114 : for (i = 1; i <= n; i++)
1760 268152 : s += F2x_degree(gel(D,i))/i;
1761 79962 : V = cgetg(s+1, t_COL);
1762 348083 : for (i = 1, s = 1; i <= n; i++)
1763 : {
1764 268122 : GEN Di = gel(D,i);
1765 268122 : long ni = F2x_degree(Di), ri = ni/i;
1766 268115 : if (ni == 0) continue;
1767 97473 : if (ni == i) { gel(V, s++) = Di; continue; }
1768 3116 : F2x_edf_simple(Di, XP, i, V, s);
1769 3116 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_edf(%ld)",i);
1770 3116 : s += ri;
1771 : }
1772 79961 : return V;
1773 : }
1774 :
1775 : static GEN
1776 64535 : F2x_factor_Cantor(GEN T)
1777 : {
1778 64535 : GEN E, F, V = F2x_factor_squarefree(T);
1779 64535 : long i, j, l = lg(V);
1780 64535 : E = cgetg(l, t_VEC);
1781 64536 : F = cgetg(l, t_VEC);
1782 256622 : for (i=1, j=1; i < l; i++)
1783 192086 : if (F2x_degree(gel(V,i)))
1784 : {
1785 79963 : GEN Fj = F2x_factor_Shoup(gel(V,i));
1786 79961 : gel(F, j) = Fj;
1787 79961 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1788 79963 : j++;
1789 : }
1790 64536 : return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
1791 : }
1792 :
1793 : #if 0
1794 : static GEN
1795 : F2x_simplefact_Shoup(GEN T)
1796 : {
1797 : long i, n, s = 0, j = 1, k;
1798 : GEN XP, D, V;
1799 : pari_timer ti;
1800 : n = F2x_degree(T);
1801 : if (DEBUGLEVEL>=6) timer_start(&ti);
1802 : XP = F2x_Frobenius(T);
1803 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1804 : D = F2x_ddf_simple(T, XP);
1805 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1806 : for (i = 1; i <= n; i++)
1807 : s += F2x_degree(gel(D,i))/i;
1808 : V = cgetg(s+1, t_VECSMALL);
1809 : for (i = 1; i <= n; i++)
1810 : {
1811 : long ni = F2x_degree(gel(D,i)), ri = ni/i;
1812 : if (ni == 0) continue;
1813 : for (k = 1; k <= ri; k++)
1814 : V[j++] = i;
1815 : }
1816 : return V;
1817 : }
1818 : static GEN
1819 : F2x_simplefact_Cantor(GEN T)
1820 : {
1821 : GEN E, F, V = F2x_factor_squarefree(T);
1822 : long i, j, l = lg(V);
1823 : F = cgetg(l, t_VEC);
1824 : E = cgetg(l, t_VEC);
1825 : for (i=1, j=1; i < l; i++)
1826 : if (F2x_degree(gel(V,i)))
1827 : {
1828 : GEN Fj = F2x_simplefact_Shoup(gel(V,i));
1829 : gel(F, j) = Fj;
1830 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1831 : j++;
1832 : }
1833 : return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
1834 : }
1835 : static int
1836 : F2x_isirred_Cantor(GEN T)
1837 : {
1838 : pari_sp av = avma;
1839 : pari_timer ti;
1840 : long n, d;
1841 : GEN dT = F2x_deriv(T);
1842 : GEN XP, D;
1843 : if (F2x_degree(F2x_gcd(T, dT)) != 0) { avma = av; return 0; }
1844 : n = F2x_degree(T);
1845 : if (DEBUGLEVEL>=6) timer_start(&ti);
1846 : XP = F2x_Frobenius(T);
1847 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1848 : D = F2x_ddf_simple(T, XP);
1849 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1850 : d = F2x_degree(gel(D, n));
1851 : avma = av; return d==n;
1852 : }
1853 : #endif
1854 :
1855 : /* driver for Cantor factorization, assume deg f > 2; not competitive for
1856 : * flag != 0, or as deg f increases */
1857 : static GEN
1858 64536 : F2x_Cantor_i(GEN f, long flag)
1859 : {
1860 : switch(flag)
1861 : {
1862 64536 : default: return F2x_factor_Cantor(f);
1863 : #if 0
1864 : case 1: return F2x_simplefact_Cantor(f);
1865 : case 2: return F2x_isirred_Cantor(f)? gen_1: NULL;
1866 : #endif
1867 : }
1868 : }
1869 : static GEN
1870 188999 : F2x_factor_i(GEN f, long flag)
1871 : {
1872 188999 : long d = F2x_degree(f);
1873 188998 : if (d <= 2) return F2x_factor_deg2(f,d,flag);
1874 74219 : return (flag == 0 && d <= 20)? F2x_Cantor_i(f, flag)
1875 140771 : : F2x_Berlekamp_i(f, flag);
1876 : }
1877 :
1878 : GEN
1879 0 : F2x_degfact(GEN f)
1880 : {
1881 0 : pari_sp av = avma;
1882 0 : GEN z = F2x_factor_i(f, 1);
1883 0 : return gerepilecopy(av, z);
1884 : }
1885 :
1886 : int
1887 238 : F2x_is_irred(GEN f) { return !!F2x_factor_i(f, 2); }
1888 :
1889 : /* Adapted from Shoup NTL */
1890 : GEN
1891 772475 : Flx_factor_squarefree(GEN f, ulong p)
1892 : {
1893 772475 : long i, q, n = degpol(f);
1894 772476 : GEN u = const_vec(n+1, pol1_Flx(f[1]));
1895 826792 : for(q = 1;;q *= p)
1896 54316 : {
1897 826792 : GEN t, v, tv, r = Flx_gcd(f, Flx_deriv(f, p), p);
1898 826793 : if (degpol(r) == 0) { gel(u, q) = f; break; }
1899 99718 : t = Flx_div(f, r, p);
1900 99718 : if (degpol(t) > 0)
1901 : {
1902 : long j;
1903 146861 : for(j = 1;;j++)
1904 : {
1905 247535 : v = Flx_gcd(r, t, p);
1906 146861 : tv = Flx_div(t, v, p);
1907 146860 : if (degpol(tv) > 0)
1908 70504 : gel(u, j*q) = Flx_normalize(tv, p);
1909 146860 : if (degpol(v) <= 0) break;
1910 100674 : r = Flx_div(r, v, p);
1911 100674 : t = v;
1912 : }
1913 46187 : if (degpol(r) == 0) break;
1914 : }
1915 54316 : f = Flx_normalize(Flx_deflate(r, p), p);
1916 : }
1917 4216880 : for (i = n; i; i--)
1918 4216880 : if (degpol(gel(u,i))) break;
1919 772477 : setlg(u,i+1); return u;
1920 : }
1921 :
1922 : long
1923 126 : Flx_ispower(GEN f, ulong k, ulong p, GEN *pt_r)
1924 : {
1925 126 : pari_sp av = avma;
1926 : ulong lc;
1927 : GEN F;
1928 126 : long i, n = degpol(f), v = f[1], l;
1929 126 : if (n % k) return 0;
1930 126 : lc = Fl_sqrtn(Flx_lead(f), k, p, NULL);
1931 126 : if (lc == ULONG_MAX) { av = avma; return 0; }
1932 126 : F = Flx_factor_squarefree(f, p); l = lg(F)-1;
1933 6405 : for (i = 1; i <= l; i++)
1934 6300 : if (i%k && degpol(gel(F,i))) { avma = av; return 0; }
1935 105 : if (pt_r)
1936 : {
1937 105 : GEN r = Fl_to_Flx(lc, v), s = pol1_Flx(v);
1938 6384 : for(i = l; i >= 1; i--)
1939 : {
1940 6279 : if (i%k) continue;
1941 1274 : s = Flx_mul(s, gel(F,i), p);
1942 1274 : r = Flx_mul(r, s, p);
1943 : }
1944 105 : *pt_r = gerepileuptoleaf(av, r);
1945 0 : } else av = avma;
1946 105 : return 1;
1947 : }
1948 :
1949 : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
1950 : static GEN
1951 1085569 : Flx_ddf_Shoup(GEN T, GEN XP, ulong p)
1952 : {
1953 1085569 : pari_sp av = avma;
1954 : GEN b, g, h, F, f, Tr, xq;
1955 : long i, j, n, v, bo, ro;
1956 : long B, l, m;
1957 : pari_timer ti;
1958 1085569 : n = get_Flx_degree(T); v = get_Flx_var(T);
1959 1085569 : if (n == 0) return cgetg(1, t_VEC);
1960 1083554 : if (n == 1) return mkvec(get_Flx_mod(T));
1961 986536 : B = n/2;
1962 986536 : l = usqrt(B);
1963 986536 : m = (B+l-1)/l;
1964 986536 : T = Flx_get_red(T, p);
1965 986535 : b = cgetg(l+2, t_VEC);
1966 986535 : gel(b, 1) = polx_Flx(v);
1967 986535 : gel(b, 2) = XP;
1968 986535 : bo = brent_kung_optpow(n, l-1, 1);
1969 986535 : ro = l<=1 ? 0:(bo-1)/(l-1) + ((n-1)/bo);
1970 986535 : if (DEBUGLEVEL>=7) timer_start(&ti);
1971 986535 : if (expu(p) <= ro)
1972 237627 : for (i = 3; i <= l+1; i++)
1973 133915 : gel(b, i) = Flxq_powu(gel(b, i-1), p, T, p);
1974 : else
1975 : {
1976 882823 : xq = Flxq_powers(gel(b, 2), bo, T, p);
1977 882824 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq baby");
1978 1027391 : for (i = 3; i <= l+1; i++)
1979 144567 : gel(b, i) = Flx_FlxqV_eval(gel(b, i-1), xq, T, p);
1980 : }
1981 986536 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: baby");
1982 986536 : xq = Flxq_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1), T, p);
1983 986536 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq giant");
1984 986536 : g = cgetg(m+1, t_VEC);
1985 986535 : gel(g, 1) = gel(xq, 2);
1986 1847628 : for(i = 2; i <= m; i++)
1987 861093 : gel(g, i) = Flx_FlxqV_eval(gel(g, i-1), xq, T, p);
1988 986535 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: giant");
1989 986535 : h = cgetg(m+1, t_VEC);
1990 2834165 : for (j = 1; j <= m; j++)
1991 : {
1992 1847629 : pari_sp av = avma;
1993 1847629 : GEN gj = gel(g, j);
1994 1847629 : GEN e = Flx_sub(gj, gel(b, 1), p);
1995 2883889 : for (i = 2; i <= l; i++)
1996 1036269 : e = Flxq_mul(e, Flx_sub(gj, gel(b, i), p), T, p);
1997 1847620 : gel(h, j) = gerepileupto(av, e);
1998 : }
1999 986536 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: diff");
2000 986536 : Tr = get_Flx_mod(T);
2001 986536 : F = cgetg(m+1, t_VEC);
2002 2834160 : for (j = 1; j <= m; j++)
2003 : {
2004 1847625 : GEN u = Flx_gcd(Tr, gel(h, j), p);
2005 1847626 : if (degpol(u))
2006 : {
2007 762574 : u = Flx_normalize(u, p);
2008 762574 : Tr = Flx_div(Tr, u, p);
2009 : }
2010 1847623 : gel(F, j) = u;
2011 : }
2012 986535 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: F");
2013 986535 : f = const_vec(n, pol1_Flx(v));
2014 2834159 : for (j = 1; j <= m; j++)
2015 : {
2016 1847623 : GEN e = gel(F, j);
2017 2038234 : for (i=l-1; i >= 0; i--)
2018 : {
2019 2038230 : GEN u = Flx_gcd(e, Flx_sub(gel(g, j), gel(b, i+1), p), p);
2020 2038233 : if (degpol(u))
2021 : {
2022 783756 : gel(f, l*j-i) = u;
2023 783756 : e = Flx_div(e, u, p);
2024 : }
2025 2038227 : if (!degpol(e)) break;
2026 : }
2027 : }
2028 986536 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: f");
2029 986536 : if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
2030 986536 : return gerepilecopy(av, f);
2031 : }
2032 :
2033 : static void
2034 128014 : Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
2035 : {
2036 128014 : long n = degpol(Tp), r = n/d;
2037 : GEN T, f, ff;
2038 : ulong p2;
2039 128014 : if (r==1) { gel(V, idx) = Tp; return; }
2040 55422 : p2 = p>>1;
2041 55422 : T = Flx_get_red(Tp, p);
2042 55422 : XP = Flx_rem(XP, T, p);
2043 : while (1)
2044 6253 : {
2045 61675 : pari_sp btop = avma;
2046 : long i;
2047 61675 : GEN g = random_Flx(n, Tp[1], p);
2048 61675 : GEN t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
2049 61674 : if (lgpol(t) == 0) continue;
2050 134633 : for(i=1; i<=10; i++)
2051 : {
2052 129968 : pari_sp btop2 = avma;
2053 129968 : GEN R = Flxq_powu(Flx_Fl_add(t, random_Fl(p), p), p2, T, p);
2054 129967 : f = Flx_gcd(Flx_Fl_add(R, p-1, p), Tp, p);
2055 129963 : if (degpol(f) > 0 && degpol(f) < n) break;
2056 74547 : avma = btop2;
2057 : }
2058 60087 : if (degpol(f) > 0 && degpol(f) < n) break;
2059 4665 : avma = btop;
2060 : }
2061 55422 : f = Flx_normalize(f, p);
2062 55422 : ff = Flx_div(Tp, f ,p);
2063 55422 : Flx_edf_simple(f, XP, d, p, V, idx);
2064 55422 : Flx_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
2065 : }
2066 : static void
2067 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx);
2068 :
2069 : static void
2070 348902 : Flx_edf_rec(GEN T, GEN XP, GEN hp, GEN t, long d, ulong p, GEN V, long idx)
2071 : {
2072 : pari_sp av;
2073 348902 : GEN Tp = get_Flx_mod(T);
2074 348902 : long n = degpol(hp), vT = Tp[1];
2075 : GEN u1, u2, f1, f2;
2076 348902 : ulong p2 = p>>1;
2077 : GEN R, h;
2078 348902 : h = Flx_get_red(hp, p);
2079 348902 : t = Flx_rem(t, T, p);
2080 348902 : av = avma;
2081 : do
2082 : {
2083 572502 : avma = av;
2084 572502 : R = Flxq_powu(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p);
2085 572502 : u1 = Flx_gcd(Flx_Fl_add(R, p-1, p), hp, p);
2086 572502 : } while (degpol(u1)==0 || degpol(u1)==n);
2087 348902 : f1 = Flx_gcd(Flx_Flxq_eval(u1, t, T, p), Tp, p);
2088 348902 : f1 = Flx_normalize(f1, p);
2089 348902 : u2 = Flx_div(hp, u1, p);
2090 348902 : f2 = Flx_div(Tp, f1, p);
2091 348902 : if (degpol(u1)==1)
2092 : {
2093 252665 : if (degpol(f1)==d)
2094 248470 : gel(V, idx) = f1;
2095 : else
2096 4195 : Flx_edf(f1, XP, d, p, V, idx);
2097 : }
2098 : else
2099 96237 : Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, V, idx);
2100 348902 : idx += degpol(f1)/d;
2101 348902 : if (degpol(u2)==1)
2102 : {
2103 250323 : if (degpol(f2)==d)
2104 246468 : gel(V, idx) = f2;
2105 : else
2106 3855 : Flx_edf(f2, XP, d, p, V, idx);
2107 : }
2108 : else
2109 98579 : Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, V, idx);
2110 348902 : }
2111 :
2112 : static void
2113 154086 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
2114 : {
2115 154086 : long n = degpol(Tp), r = n/d, vT = Tp[1];
2116 : GEN T, h, t;
2117 : pari_timer ti;
2118 154086 : if (r==1) { gel(V, idx) = Tp; return; }
2119 154086 : T = Flx_get_red(Tp, p);
2120 154086 : XP = Flx_rem(XP, T, p);
2121 154086 : if (DEBUGLEVEL>=7) timer_start(&ti);
2122 : do
2123 : {
2124 158993 : GEN g = random_Flx(n, vT, p);
2125 158993 : t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
2126 158993 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
2127 158993 : h = Flxq_minpoly(t, T, p);
2128 158993 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
2129 158993 : } while (degpol(h) <= 1);
2130 154086 : Flx_edf_rec(T, XP, h, t, d, p, V, idx);
2131 : }
2132 :
2133 : static GEN
2134 440537 : Flx_factor_Shoup(GEN T, ulong p)
2135 : {
2136 440537 : long i, n, s = 0;
2137 : GEN XP, D, V;
2138 440537 : long e = expu(p);
2139 : pari_timer ti;
2140 440537 : n = get_Flx_degree(T);
2141 440537 : T = Flx_get_red(T, p);
2142 440537 : if (DEBUGLEVEL>=6) timer_start(&ti);
2143 440537 : XP = Flx_Frobenius(T, p);
2144 440538 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2145 440538 : D = Flx_ddf_Shoup(T, XP, p);
2146 440538 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2147 440538 : s = ddf_to_nbfact(D);
2148 440538 : V = cgetg(s+1, t_COL);
2149 2481985 : for (i = 1, s = 1; i <= n; i++)
2150 : {
2151 2041447 : GEN Di = gel(D,i);
2152 2041447 : long ni = degpol(Di), ri = ni/i;
2153 2041447 : if (ni == 0) continue;
2154 555913 : Di = Flx_normalize(Di, p);
2155 555913 : if (ni == i) { gel(V, s++) = Di; continue; }
2156 163206 : if (ri <= e*expu(e))
2157 146036 : Flx_edf(Di, XP, i, p, V, s);
2158 : else
2159 17170 : Flx_edf_simple(Di, XP, i, p, V, s);
2160 163206 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
2161 163206 : s += ri;
2162 : }
2163 440538 : return V;
2164 : }
2165 :
2166 : static GEN
2167 416934 : Flx_factor_Cantor(GEN T, ulong p)
2168 : {
2169 416934 : GEN E, F, V = Flx_factor_squarefree(get_Flx_mod(T), p);
2170 416935 : long i, j, l = lg(V);
2171 416935 : F = cgetg(l, t_VEC);
2172 416935 : E = cgetg(l, t_VEC);
2173 1051205 : for (i=1, j=1; i < l; i++)
2174 634270 : if (degpol(gel(V,i)))
2175 : {
2176 440537 : GEN Fj = Flx_factor_Shoup(gel(V,i), p);
2177 440538 : gel(F, j) = Fj;
2178 440538 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
2179 440538 : j++;
2180 : }
2181 416935 : return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
2182 : }
2183 :
2184 : GEN
2185 0 : Flx_ddf(GEN T, ulong p)
2186 : {
2187 : GEN XP;
2188 0 : T = Flx_get_red(T, p);
2189 0 : XP = Flx_Frobenius(T, p);
2190 0 : return ddf_to_ddf2(Flx_ddf_Shoup(T, XP, p));
2191 : }
2192 :
2193 : static GEN
2194 354639 : Flx_simplefact_Cantor(GEN T, ulong p)
2195 : {
2196 : GEN XP, V;
2197 : long i, l;
2198 354639 : T = Flx_get_red(T, p);
2199 354639 : XP = Flx_Frobenius(T, p);
2200 354639 : V = Flx_factor_squarefree(get_Flx_mod(T), p); l = lg(V);
2201 354639 : for (i=1; i < l; i++) gel(V,i) = Flx_ddf_Shoup(gel(V,i), XP, p);
2202 354639 : return vddf_to_simplefact(V, get_Flx_degree(T));
2203 : }
2204 :
2205 : static int
2206 1057 : Flx_isirred_Cantor(GEN Tp, ulong p)
2207 : {
2208 1057 : pari_sp av = avma;
2209 : pari_timer ti;
2210 : long n, d;
2211 1057 : GEN T = get_Flx_mod(Tp);
2212 1057 : GEN dT = Flx_deriv(T, p);
2213 : GEN XP, D;
2214 1057 : if (degpol(Flx_gcd(T, dT, p)) != 0) { avma = av; return 0; }
2215 770 : n = get_Flx_degree(T);
2216 770 : T = Flx_get_red(Tp, p);
2217 770 : if (DEBUGLEVEL>=6) timer_start(&ti);
2218 770 : XP = Flx_Frobenius(T, p);
2219 770 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2220 770 : D = Flx_ddf_Shoup(T, XP, p);
2221 770 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2222 770 : d = degpol(gel(D, n));
2223 770 : avma = av; return d==n;
2224 : }
2225 :
2226 : /* f monic */
2227 : static GEN
2228 1985555 : Flx_factor_i(GEN f, ulong pp, long flag)
2229 : {
2230 : long d;
2231 1985555 : if (pp==2) { /*We need to handle 2 specially */
2232 28861 : GEN F = F2x_factor_i(Flx_to_F2x(f),flag);
2233 28861 : if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
2234 28861 : return F;
2235 : }
2236 1956694 : d = degpol(f);
2237 1956694 : if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
2238 772630 : switch(flag)
2239 : {
2240 416934 : default: return Flx_factor_Cantor(f, pp);
2241 354639 : case 1: return Flx_simplefact_Cantor(f, pp);
2242 1057 : case 2: return Flx_isirred_Cantor(f, pp)? gen_1: NULL;
2243 : }
2244 : }
2245 :
2246 : GEN
2247 1115196 : Flx_degfact(GEN f, ulong p)
2248 : {
2249 1115196 : pari_sp av = avma;
2250 1115196 : GEN z = Flx_factor_i(Flx_normalize(f,p),p,1);
2251 1115196 : return gerepilecopy(av, z);
2252 : }
2253 :
2254 : /* T must be squarefree mod p*/
2255 : GEN
2256 191992 : Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
2257 : {
2258 : GEN XP, D;
2259 : pari_timer ti;
2260 191992 : long i, s, n = get_Flx_degree(T);
2261 191992 : GEN V = const_vecsmall(n, 0);
2262 191992 : pari_sp av = avma;
2263 191992 : T = Flx_get_red(T, p);
2264 191992 : if (DEBUGLEVEL>=6) timer_start(&ti);
2265 191992 : XP = Flx_Frobenius(T, p);
2266 191992 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2267 191992 : D = Flx_ddf_Shoup(T, XP, p);
2268 191992 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2269 1228714 : for (i = 1, s = 0; i <= n; i++)
2270 : {
2271 1036722 : V[i] = degpol(gel(D,i))/i;
2272 1036722 : s += V[i];
2273 : }
2274 191992 : *nb = s;
2275 191992 : avma = av; return V;
2276 : }
2277 :
2278 : long
2279 94626 : Flx_nbfact_Frobenius(GEN T, GEN XP, ulong p)
2280 : {
2281 94626 : pari_sp av = avma;
2282 94626 : long s = ddf_to_nbfact(Flx_ddf_Shoup(T, XP, p));
2283 94626 : avma = av; return s;
2284 : }
2285 :
2286 : /* T must be squarefree mod p*/
2287 : long
2288 94626 : Flx_nbfact(GEN T, ulong p)
2289 : {
2290 94626 : pari_sp av = avma;
2291 94626 : GEN XP = Flx_Frobenius(T, p);
2292 94626 : long n = Flx_nbfact_Frobenius(T, XP, p);
2293 94626 : avma = av; return n;
2294 : }
2295 :
2296 : int
2297 1057 : Flx_is_irred(GEN f, ulong p)
2298 : {
2299 1057 : pari_sp av = avma;
2300 1057 : int z = !!Flx_factor_i(Flx_normalize(f,p),p,2);
2301 1057 : avma = av; return z;
2302 : }
2303 :
2304 : /* Use this function when you think f is reducible, and that there are lots of
2305 : * factors. If you believe f has few factors, use FpX_nbfact(f,p)==1 instead */
2306 : int
2307 77 : FpX_is_irred(GEN f, GEN p)
2308 : {
2309 77 : pari_sp av = avma;
2310 : int z;
2311 77 : switch(ZX_factmod_init(&f,p))
2312 : {
2313 14 : case 0: z = !!F2x_factor_i(f,2); break;
2314 56 : case 1: z = !!Flx_factor_i(f,p[2],2); break;
2315 7 : default: z = !!FpX_factor_i(f,p,2); break;
2316 : }
2317 77 : avma = av; return z;
2318 : }
2319 : GEN
2320 42 : FpX_degfact(GEN f, GEN p) {
2321 42 : pari_sp av = avma;
2322 : GEN F;
2323 42 : switch(ZX_factmod_init(&f,p))
2324 : {
2325 7 : case 0: F = F2x_factor_i(f,1); break;
2326 7 : case 1: F = Flx_factor_i(f,p[2],1); break;
2327 28 : default: F = FpX_factor_i(f,p,1); break;
2328 : }
2329 42 : return gerepilecopy(av, F);
2330 : }
2331 :
2332 : #if 0
2333 : /* set x <-- x + c*y mod p */
2334 : /* x is not required to be normalized.*/
2335 : static void
2336 : Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
2337 : {
2338 : long i, lx, ly;
2339 : ulong *x=(ulong *)gx;
2340 : ulong *y=(ulong *)gy;
2341 : if (!c) return;
2342 : lx = lg(gx);
2343 : ly = lg(gy);
2344 : if (lx<ly) pari_err_BUG("lx<ly in Flx_addmul_inplace");
2345 : if (SMALL_ULONG(p))
2346 : for (i=2; i<ly; i++) x[i] = (x[i] + c*y[i]) % p;
2347 : else
2348 : for (i=2; i<ly; i++) x[i] = Fl_add(x[i], Fl_mul(c,y[i],p),p);
2349 : }
2350 : #endif
2351 :
2352 : GEN
2353 874481 : FpX_factor(GEN f, GEN p)
2354 : {
2355 874481 : pari_sp av = avma;
2356 : GEN F;
2357 874481 : switch(ZX_factmod_init(&f, p))
2358 : {
2359 144674 : case 0: F = F2x_factor_i(f,0);
2360 144674 : F2xV_to_ZXV_inplace(gel(F,1)); break;
2361 728457 : case 1: F = Flx_factor_i(f,p[2],0);
2362 728457 : FlxV_to_ZXV_inplace(gel(F,1)); break;
2363 1350 : default: F = FpX_factor_i(f,p,0); break;
2364 : }
2365 874481 : return gerepilecopy(av, F);
2366 : }
2367 :
2368 : GEN
2369 140781 : Flx_factor(GEN f, ulong p)
2370 : {
2371 140781 : pari_sp av = avma;
2372 140781 : return gerepilecopy(av, Flx_factor_i(Flx_normalize(f,p),p,0));
2373 : }
2374 : GEN
2375 15205 : F2x_factor(GEN f)
2376 : {
2377 15205 : pari_sp av = avma;
2378 15205 : return gerepilecopy(av, F2x_factor_i(f,0));
2379 : }
|
