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; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : #define DEBUGLEVEL DEBUGLEVEL_factormod
19 :
20 : /***********************************************************************/
21 : /** **/
22 : /** Factorisation over finite field **/
23 : /** **/
24 : /***********************************************************************/
25 :
26 : /*******************************************************************/
27 : /* */
28 : /* ROOTS MODULO a prime p (no multiplicities) */
29 : /* */
30 : /*******************************************************************/
31 : /* Replace F by a monic normalized FpX having the same factors;
32 : * assume p prime and *F a ZX */
33 : static int
34 3459405 : ZX_factmod_init(GEN *F, GEN p)
35 : {
36 3459405 : if (lgefint(p) == 3)
37 : {
38 3457715 : ulong pp = p[2];
39 3457715 : if (pp == 2) { *F = ZX_to_F2x(*F); return 0; }
40 2079705 : *F = ZX_to_Flx(*F, pp);
41 2079795 : if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
42 2079786 : return 1;
43 : }
44 1690 : *F = FpX_red(*F, p);
45 1695 : if (lg(*F) > 3) *F = FpX_normalize(*F, p);
46 1695 : return 2;
47 : }
48 : static void
49 367625 : ZX_rootmod_init(GEN *F, GEN p)
50 : {
51 367625 : if (lgefint(p) == 3)
52 : {
53 359229 : ulong pp = p[2];
54 359229 : *F = ZX_to_Flx(*F, pp);
55 359225 : if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
56 : }
57 : else
58 : {
59 8396 : *F = FpX_red(*F, p);
60 8395 : if (lg(*F) > 3) *F = FpX_normalize(*F, p);
61 : }
62 367621 : }
63 :
64 : /* return 1,...,p-1 [not_0 = 1] or 0,...,p [not_0 = 0] */
65 : static GEN
66 434 : all_roots_mod_p(ulong p, int not_0)
67 : {
68 : GEN r;
69 : ulong i;
70 434 : if (not_0) {
71 329 : r = cgetg(p, t_VECSMALL);
72 987 : for (i = 1; i < p; i++) r[i] = i;
73 : } else {
74 105 : r = cgetg(p+1, t_VECSMALL);
75 448 : for (i = 0; i < p; i++) r[i+1] = i;
76 : }
77 434 : return r;
78 : }
79 :
80 : /* X^n - 1 */
81 : static GEN
82 3188 : Flx_Xnm1(long sv, long n, ulong p)
83 : {
84 3188 : GEN t = cgetg(n+3, t_VECSMALL);
85 : long i;
86 3188 : t[1] = sv;
87 3188 : t[2] = p - 1;
88 9824 : for (i = 3; i <= n+1; i++) t[i] = 0;
89 3188 : t[i] = 1; return t;
90 : }
91 : /* X^n + 1 */
92 : static GEN
93 3135 : Flx_Xn1(long sv, long n, ulong p)
94 : {
95 3135 : GEN t = cgetg(n+3, t_VECSMALL);
96 : long i;
97 : (void) p;
98 3135 : t[1] = sv;
99 3135 : t[2] = 1;
100 9718 : for (i = 3; i <= n+1; i++) t[i] = 0;
101 3135 : t[i] = 1; return t;
102 : }
103 :
104 : static GEN
105 89430 : Flx_root_mod_2(GEN f)
106 : {
107 89430 : int z1, z0 = !(f[2] & 1);
108 : long i,n;
109 : GEN y;
110 :
111 343641 : for (i=2, n=1; i < lg(f); i++) n += f[i];
112 89430 : z1 = n & 1;
113 89430 : y = cgetg(z0+z1+1, t_VECSMALL); i = 1;
114 89428 : if (z0) y[i++] = 0;
115 89428 : if (z1) y[i ] = 1;
116 89428 : return y;
117 : }
118 : static ulong
119 91 : Flx_oneroot_mod_2(GEN f)
120 : {
121 : long i,n;
122 91 : if (!(f[2] & 1)) return 0;
123 364 : for (i=2, n=1; i < lg(f); i++) n += f[i];
124 91 : if (n & 1) return 1;
125 28 : return 2;
126 : }
127 :
128 : static GEN FpX_roots_i(GEN f, GEN p);
129 : static GEN Flx_roots_i(GEN f, ulong p);
130 :
131 : static int
132 16049042 : cmpGuGu(GEN a, GEN b) { return (ulong)a < (ulong)b? -1: (a == b? 0: 1); }
133 :
134 : /* Generic driver to computes the roots of f modulo pp, using 'Roots' when
135 : * pp is a small prime.
136 : * if (gpwrap), check types thoroughly and return t_INTMODs, otherwise
137 : * assume that f is an FpX, pp a prime and return t_INTs */
138 : static GEN
139 354476 : rootmod_aux(GEN f, GEN pp)
140 : {
141 : GEN y;
142 354476 : switch(lg(f))
143 : {
144 14 : case 2: pari_err_ROOTS0("rootmod");
145 46221 : case 3: return cgetg(1,t_COL);
146 : }
147 308241 : if (typ(f) == t_VECSMALL)
148 : {
149 305105 : ulong p = pp[2];
150 305105 : if (p == 2)
151 89423 : y = Flx_root_mod_2(f);
152 : else
153 : {
154 215682 : if (!odd(p)) pari_err_PRIME("rootmod",utoi(p));
155 215682 : y = Flx_roots_i(f, p);
156 : }
157 305099 : y = Flc_to_ZC(y);
158 : }
159 : else
160 3136 : y = FpX_roots_i(f, pp);
161 308229 : return y;
162 : }
163 : /* assume that f is a ZX and p a prime */
164 : GEN
165 354479 : FpX_roots(GEN f, GEN p)
166 : {
167 354479 : pari_sp av = avma;
168 354479 : GEN y; ZX_rootmod_init(&f, p); y = rootmod_aux(f, p);
169 354450 : return gerepileupto(av, y);
170 : }
171 :
172 : /* assume x reduced mod p > 2, monic. */
173 : static int
174 21 : FpX_quad_factortype(GEN x, GEN p)
175 : {
176 21 : GEN b = gel(x,3), c = gel(x,2);
177 21 : GEN D = subii(sqri(b), shifti(c,2));
178 21 : return kronecker(D,p);
179 : }
180 : /* assume x reduced mod p, monic. Return one root, or NULL if irreducible */
181 : static GEN
182 7706 : FpX_quad_root(GEN x, GEN p, int unknown)
183 : {
184 7706 : GEN s, D, b = gel(x,3), c = gel(x,2);
185 :
186 7706 : if (absequaliu(p, 2)) {
187 0 : if (!signe(b)) return c;
188 0 : return signe(c)? NULL: gen_1;
189 : }
190 7706 : D = subii(sqri(b), shifti(c,2));
191 7706 : D = remii(D,p);
192 7706 : if (unknown && kronecker(D,p) == -1) return NULL;
193 :
194 7172 : s = Fp_sqrt(D,p);
195 : /* p is not prime, go on and give e.g. maxord a chance to recover */
196 7172 : if (!s) return NULL;
197 7164 : return Fp_halve(Fp_sub(s,b, p), p);
198 : }
199 : static GEN
200 3222 : FpX_otherroot(GEN x, GEN r, GEN p)
201 3222 : { return Fp_neg(Fp_add(gel(x,3), r, p), p); }
202 :
203 : /* disc(x^2+bx+c) = b^2 - 4c */
204 : static ulong
205 25286197 : Fl_disc_bc(ulong b, ulong c, ulong p)
206 25286197 : { return Fl_sub(Fl_sqr(b,p), Fl_double(Fl_double(c,p),p), p); }
207 : /* p > 2 */
208 : static ulong
209 22821997 : Flx_quad_root(GEN x, ulong p, int unknown)
210 : {
211 22821997 : ulong s, b = x[3], c = x[2];
212 22821997 : ulong D = Fl_disc_bc(b, c, p);
213 22797285 : if (unknown && krouu(D,p) == -1) return p;
214 15257236 : s = Fl_sqrt(D,p);
215 15314611 : if (s==~0UL) return p;
216 15314598 : return Fl_halve(Fl_sub(s,b, p), p);
217 : }
218 : static ulong
219 13832601 : Flx_otherroot(GEN x, ulong r, ulong p)
220 13832601 : { return Fl_neg(Fl_add(x[3], r, p), p); }
221 :
222 : /* 'todo' contains the list of factors to be split.
223 : * 'done' the list of finished factors, no longer touched */
224 : struct split_t { GEN todo, done; };
225 : static void
226 5022311 : split_init(struct split_t *S, long max)
227 : {
228 5022311 : S->todo = vectrunc_init(max);
229 5021857 : S->done = vectrunc_init(max);
230 5021446 : }
231 : #if 0
232 : /* move todo[i] to done */
233 : static void
234 : split_convert(struct split_t *S, long i)
235 : {
236 : long n = lg(S->todo)-1;
237 : vectrunc_append(S->done, gel(S->todo,i));
238 : if (n) gel(S->todo,i) = gel(S->todo, n);
239 : setlg(S->todo, n);
240 : }
241 : #endif
242 : /* append t to todo */
243 : static void
244 5336575 : split_add(struct split_t *S, GEN t) { vectrunc_append(S->todo, t); }
245 : /* delete todo[i], add t to done */
246 : static void
247 5336826 : split_moveto_done(struct split_t *S, long i, GEN t)
248 : {
249 5336826 : long n = lg(S->todo)-1;
250 5336826 : vectrunc_append(S->done, t);
251 5337304 : if (n) gel(S->todo,i) = gel(S->todo, n);
252 5337304 : setlg(S->todo, n);
253 :
254 5337193 : }
255 : /* append t to done */
256 : static void
257 484551 : split_add_done(struct split_t *S, GEN t)
258 484551 : { vectrunc_append(S->done, t); }
259 : /* split todo[i] into a and b */
260 : static void
261 404803 : split_todo(struct split_t *S, long i, GEN a, GEN b)
262 : {
263 404803 : gel(S->todo, i) = a;
264 404803 : split_add(S, b);
265 404804 : }
266 : /* split todo[i] into a and b, moved to done */
267 : static void
268 449281 : split_done(struct split_t *S, long i, GEN a, GEN b)
269 : {
270 449281 : split_moveto_done(S, i, a);
271 449287 : split_add_done(S, b);
272 449287 : }
273 :
274 : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
275 : static GEN
276 3136 : FpX_roots_i(GEN f, GEN p)
277 : {
278 : GEN pol, pol0, a, q;
279 : struct split_t S;
280 :
281 3136 : split_init(&S, lg(f)-1);
282 3136 : settyp(S.done, t_COL);
283 3136 : if (ZX_valrem(f, &f)) split_add_done(&S, gen_0);
284 3136 : switch(degpol(f))
285 : {
286 7 : case 0: return ZC_copy(S.done);
287 14 : case 1: split_add_done(&S, subii(p, gel(f,2))); return ZC_copy(S.done);
288 1744 : case 2: {
289 1744 : GEN s, r = FpX_quad_root(f, p, 1);
290 1744 : if (r) {
291 1744 : split_add_done(&S, r);
292 1744 : s = FpX_otherroot(f,r, p);
293 : /* f not known to be square free yet */
294 1744 : if (!equalii(r, s)) split_add_done(&S, s);
295 : }
296 1744 : return sort(S.done);
297 : }
298 : }
299 :
300 1371 : a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
301 1371 : if (lg(a) < 3) pari_err_PRIME("rootmod",p);
302 1371 : a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
303 1371 : a = FpX_gcd(f,a, p);
304 1371 : if (!degpol(a)) return ZC_copy(S.done);
305 1371 : split_add(&S, FpX_normalize(a,p));
306 :
307 1371 : q = shifti(p,-1);
308 1371 : pol0 = icopy(gen_1); /* constant term, will vary in place */
309 1371 : pol = deg1pol_shallow(gen_1, pol0, varn(f));
310 1371 : for (pol0[2] = 1;; pol0[2]++)
311 1501 : {
312 2872 : long j, l = lg(S.todo);
313 2872 : if (l == 1) return sort(S.done);
314 1508 : if (pol0[2] == 100 && !BPSW_psp(p)) pari_err_PRIME("polrootsmod",p);
315 3132 : for (j = 1; j < l; j++)
316 : {
317 1631 : GEN b, r, s, c = gel(S.todo,j);
318 1631 : switch(degpol(c))
319 : { /* convert linear and quadratics to roots, try to split the rest */
320 196 : case 1:
321 196 : split_moveto_done(&S, j, subii(p, gel(c,2)));
322 196 : j--; l--; break;
323 1298 : case 2:
324 1298 : r = FpX_quad_root(c, p, 0);
325 1298 : if (!r) pari_err_PRIME("polrootsmod",p);
326 1291 : s = FpX_otherroot(c,r, p);
327 1291 : split_done(&S, j, r, s);
328 1291 : j--; l--; break;
329 137 : default:
330 137 : b = FpXQ_pow(pol,q, c,p);
331 137 : if (degpol(b) <= 0) continue;
332 123 : b = FpX_gcd(c,FpX_Fp_sub_shallow(b,gen_1,p), p);
333 123 : if (!degpol(b)) continue;
334 123 : b = FpX_normalize(b, p);
335 123 : c = FpX_div(c,b, p);
336 123 : split_todo(&S, j, b, c);
337 : }
338 : }
339 : }
340 : }
341 :
342 : /* Assume f is normalized */
343 : static ulong
344 180228 : Flx_cubic_root(GEN ff, ulong p)
345 : {
346 180228 : GEN f = Flx_normalize(ff,p);
347 180228 : ulong pi = get_Fl_red(p);
348 180228 : ulong a = f[4], b=f[3], c=f[2], p3 = p%3==1 ? (2*p+1)/3 :(p+1)/3;
349 180228 : ulong t = Fl_mul_pre(a, p3, p, pi), t2 = Fl_sqr_pre(t, p, pi);
350 180229 : ulong A = Fl_sub(b, Fl_triple(t2, p), p);
351 180228 : ulong B = Fl_addmul_pre(c, t, Fl_sub(Fl_double(t2, p), b, p), p, pi);
352 180229 : ulong A3 = Fl_mul_pre(A, p3, p, pi);
353 180229 : ulong A32 = Fl_sqr_pre(A3, p, pi), A33 = Fl_mul_pre(A3, A32, p, pi);
354 180229 : ulong S = Fl_neg(B,p), P = Fl_neg(A3,p);
355 180229 : ulong D = Fl_add(Fl_sqr_pre(S, p, pi), Fl_double(Fl_double(A33, p), p), p);
356 180229 : if (krouu(D,p) >= 0)
357 : {
358 115363 : ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
359 115363 : ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
360 115363 : if (p%3==2) /* 1 solutions */
361 20567 : vS1 = Fl_powu_pre(S1, (2*p-1)/3, p, pi);
362 : else
363 : {
364 94796 : vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
365 94794 : if (vS1==~0UL) return p; /*0 solutions*/
366 : /*3 solutions*/
367 : }
368 83875 : vS2 = P? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): 0;
369 83877 : return Fl_sub(Fl_add(vS1,vS2, p), t, p);
370 : }
371 : else
372 : {
373 64866 : pari_sp av = avma;
374 64866 : GEN S1 = mkvecsmall2(Fl_halve(S, p), Fl_halve(1UL, p));
375 64866 : GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
376 : ulong Sa;
377 64866 : if (!vS1) return p; /*0 solutions, p%3==2*/
378 64866 : Sa = vS1[1];
379 64866 : if (p%3==1) /*1 solutions*/
380 : {
381 24235 : ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
382 24235 : if (Fa!=P)
383 15952 : Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
384 : }
385 64866 : set_avma(av);
386 64866 : return Fl_sub(Fl_double(Sa,p),t,p);
387 : }
388 : }
389 :
390 : /* Assume f is normalized */
391 : static GEN
392 126 : FpX_cubic_root(GEN ff, GEN p)
393 : {
394 126 : GEN f = FpX_normalize(ff,p);
395 126 : GEN a = gel(f,4), b = gel(f,3), c = gel(f,2);
396 126 : ulong pm3 = umodiu(p,3);
397 28 : GEN p3 = pm3==1 ? diviuexact(addiu(shifti(p,1),1),3)
398 126 : : diviuexact(addiu(p,1),3);
399 126 : GEN t = Fp_mul(a, p3, p), t2 = Fp_sqr(t, p);
400 126 : GEN A = Fp_sub(b, Fp_mulu(t2, 3, p), p);
401 126 : GEN B = Fp_addmul(c, t, Fp_sub(shifti(t2, 1), b, p), p);
402 126 : GEN A3 = Fp_mul(A, p3, p);
403 126 : GEN A32 = Fp_sqr(A3, p), A33 = Fp_mul(A3, A32, p);
404 126 : GEN S = Fp_neg(B,p), P = Fp_neg(A3,p);
405 126 : GEN D = Fp_add(Fp_sqr(S, p), shifti(A33, 2), p);
406 126 : if (kronecker(D,p) >= 0)
407 : {
408 28 : GEN s = Fp_sqrt(D, p), vS1, vS2;
409 28 : GEN S1 = S==s ? S: Fp_halve(Fp_sub(S, s, p), p);
410 28 : if (pm3 == 2) /* 1 solutions */
411 0 : vS1 = Fp_pow(S1, diviuexact(addiu(shifti(p, 1), 1), 3), p);
412 : else
413 : {
414 28 : vS1 = Fp_sqrtn(S1, utoi(3), p, NULL);
415 28 : if (!vS1) return p; /*0 solutions*/
416 : /*3 solutions*/
417 : }
418 28 : vS2 = P? Fp_mul(P, Fp_inv(vS1, p), p): 0;
419 28 : return Fp_sub(Fp_add(vS1,vS2, p), t, p);
420 : }
421 : else
422 : {
423 98 : pari_sp av = avma;
424 98 : GEN T = deg2pol_shallow(gen_1, gen_0, negi(D), 0);
425 98 : GEN S1 = deg1pol_shallow(Fp_halve(gen_1, p), Fp_halve(S, p), 0);
426 98 : GEN vS1 = FpXQ_sqrtn(S1, utoi(3), T, p, NULL);
427 : GEN Sa;
428 98 : if (!vS1) return p; /*0 solutions, p%3==2*/
429 98 : Sa = gel(vS1,2);
430 98 : if (pm3 == 1) /*1 solutions*/
431 : {
432 0 : GEN Fa = FpXQ_norm(vS1, T, p);
433 0 : if (!equalii(Fa,P))
434 0 : Sa = Fp_mul(Sa, Fp_div(Fa, P, p),p);
435 : }
436 98 : set_avma(av);
437 98 : return Fp_sub(shifti(Sa,1),t,p);
438 : }
439 : }
440 :
441 : /* assume p > 2 prime */
442 : static ulong
443 3248591 : Flx_oneroot_i(GEN f, ulong p, long fl)
444 : {
445 : GEN pol, a;
446 : ulong q;
447 : long da;
448 :
449 3248591 : if (Flx_val(f)) return 0;
450 3247807 : switch(degpol(f))
451 : {
452 13014 : case 1: return Fl_neg(f[2], p);
453 2795644 : case 2: return Flx_quad_root(f, p, 1);
454 164308 : case 3: if (p>3) return Flx_cubic_root(f, p); /*FALL THROUGH*/
455 : }
456 :
457 274841 : if (!fl)
458 : {
459 242297 : a = Flxq_powu(polx_Flx(f[1]), p - 1, f,p);
460 242297 : if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
461 242297 : a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
462 242297 : a = Flx_gcd(f,a, p);
463 32544 : } else a = f;
464 274841 : da = degpol(a);
465 274838 : if (!da) return p;
466 196505 : a = Flx_normalize(a,p);
467 :
468 196505 : q = p >> 1;
469 196505 : pol = polx_Flx(f[1]);
470 196504 : for(pol[2] = 1;; pol[2]++)
471 : {
472 299873 : if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
473 299874 : switch(da)
474 : {
475 123298 : case 1: return Fl_neg(a[2], p);
476 57293 : case 2: return Flx_quad_root(a, p, 0);
477 15920 : case 3: if (p>3) return Flx_cubic_root(a, p); /*FALL THROUGH*/
478 : default: {
479 103363 : GEN b = Flxq_powu(pol,q, a,p);
480 : long db;
481 103347 : if (degpol(b) <= 0) continue;
482 98391 : b = Flx_gcd(a,Flx_Fl_add(b,p-1,p), p);
483 98395 : db = degpol(b); if (!db) continue;
484 98394 : b = Flx_normalize(b, p);
485 98402 : if (db <= (da >> 1)) {
486 60388 : a = b;
487 60388 : da = db;
488 : } else {
489 38014 : a = Flx_div(a,b, p);
490 38021 : da -= db;
491 : }
492 : }
493 : }
494 : }
495 : }
496 :
497 : /* assume p > 3 prime */
498 : static GEN
499 5245 : FpX_oneroot_i(GEN f, GEN p)
500 : {
501 : GEN pol, pol0, a, q;
502 : long da;
503 :
504 5245 : if (ZX_val(f)) return gen_0;
505 4923 : switch(degpol(f))
506 : {
507 814 : case 1: return subii(p, gel(f,2));
508 3914 : case 2: return FpX_quad_root(f, p, 1);
509 126 : case 3: return FpX_cubic_root(f, p);
510 : }
511 :
512 70 : a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
513 70 : if (lg(a) < 3) pari_err_PRIME("rootmod",p);
514 70 : a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
515 70 : a = FpX_gcd(f,a, p);
516 70 : da = degpol(a);
517 70 : if (!da) return NULL;
518 70 : a = FpX_normalize(a,p);
519 :
520 70 : q = shifti(p,-1);
521 70 : pol0 = icopy(gen_1); /* constant term, will vary in place */
522 70 : pol = deg1pol_shallow(gen_1, pol0, varn(f));
523 70 : for (pol0[2]=1; ; pol0[2]++)
524 : {
525 224 : if (pol0[2] == 1000 && !BPSW_psp(p)) pari_err_PRIME("FpX_oneroot",p);
526 224 : switch(da)
527 : {
528 42 : case 1: return subii(p, gel(a,2));
529 28 : case 2: return FpX_quad_root(a, p, 0);
530 154 : default: {
531 154 : GEN b = FpXQ_pow(pol,q, a,p);
532 : long db;
533 154 : if (degpol(b) <= 0) continue;
534 147 : b = FpX_gcd(a,FpX_Fp_sub_shallow(b,gen_1,p), p);
535 147 : db = degpol(b); if (!db) continue;
536 147 : b = FpX_normalize(b, p);
537 147 : if (db <= (da >> 1)) {
538 105 : a = b;
539 105 : da = db;
540 : } else {
541 42 : a = FpX_div(a,b, p);
542 42 : da -= db;
543 : }
544 : }
545 : }
546 : }
547 : }
548 :
549 : ulong
550 3206776 : Flx_oneroot(GEN f, ulong p)
551 : {
552 3206776 : pari_sp av = avma;
553 : ulong r;
554 3206776 : switch(lg(f))
555 : {
556 0 : case 2: return 0;
557 0 : case 3: return p;
558 : }
559 3206776 : if (p == 2) return Flx_oneroot_mod_2(f);
560 3206776 : r = Flx_oneroot_i(Flx_normalize(f, p), p, 0);
561 3206776 : return gc_ulong(av,r);
562 : }
563 :
564 : ulong
565 34006 : Flx_oneroot_split(GEN f, ulong p)
566 : {
567 34006 : pari_sp av = avma;
568 : ulong r;
569 34006 : switch(lg(f))
570 : {
571 0 : case 2: return 0;
572 0 : case 3: return p;
573 : }
574 34006 : if (p == 2) return Flx_oneroot_mod_2(f);
575 34006 : r = Flx_oneroot_i(Flx_normalize(f, p), p, 1);
576 34010 : return gc_ulong(av,r);
577 : }
578 :
579 : /* assume that p is prime */
580 : GEN
581 13146 : FpX_oneroot(GEN f, GEN pp) {
582 13146 : pari_sp av = avma;
583 13146 : ZX_rootmod_init(&f, pp);
584 13145 : switch(lg(f))
585 : {
586 0 : case 2: set_avma(av); return gen_0;
587 0 : case 3: return gc_NULL(av);
588 : }
589 13145 : if (typ(f) == t_VECSMALL)
590 : {
591 7900 : ulong r, p = pp[2];
592 7900 : if (p == 2)
593 91 : r = Flx_oneroot_mod_2(f);
594 : else
595 7809 : r = Flx_oneroot_i(f, p, 0);
596 7900 : set_avma(av);
597 7900 : return (r == p)? NULL: utoi(r);
598 : }
599 5245 : f = FpX_oneroot_i(f, pp);
600 5246 : if (!f) return gc_NULL(av);
601 5246 : return gerepileuptoint(av, f);
602 : }
603 :
604 : /* returns a root of unity in F_p that is suitable for finding a factor */
605 : /* of degree deg_factor of a polynomial of degree deg; the order is */
606 : /* returned in n */
607 : /* A good choice seems to be n close to deg/deg_factor; we choose n */
608 : /* twice as big and decrement until it divides p-1. */
609 : static GEN
610 154 : good_root_of_unity(GEN p, long deg, long deg_factor, long *pt_n)
611 : {
612 154 : pari_sp ltop = avma;
613 : GEN pm, factn, power, base, zeta;
614 : long n;
615 :
616 154 : pm = subis (p, 1ul);
617 336 : for (n = deg / 2 / deg_factor + 1; !dvdiu (pm, n); n--);
618 154 : factn = Z_factor(stoi(n));
619 154 : power = diviuexact (pm, n);
620 154 : base = gen_1;
621 : do {
622 259 : base = addis (base, 1l);
623 259 : zeta = Fp_pow (base, power, p);
624 : }
625 259 : while (!equaliu (Fp_order (zeta, factn, p), n));
626 154 : *pt_n = n;
627 154 : return gerepileuptoint (ltop, zeta);
628 : }
629 :
630 : GEN
631 1092 : FpX_oneroot_split(GEN fact, GEN p)
632 : {
633 1092 : pari_sp av = avma;
634 : long n, deg_f, i, dmin;
635 : GEN prim, expo, minfactor, xplusa, zeta, xpow;
636 1092 : fact = FpX_normalize(fact, p);
637 1092 : deg_f = degpol(fact);
638 1092 : if (deg_f <= 3) return FpX_oneroot(fact, p);
639 133 : minfactor = fact; /* factor of minimal degree found so far */
640 133 : dmin = degpol(minfactor);
641 133 : xplusa = pol_x(varn(fact));
642 287 : while (dmin > 3)
643 : {
644 : /* split minfactor by computing its gcd with (X+a)^exp-zeta, where */
645 : /* zeta varies over the roots of unity in F_p */
646 154 : fact = minfactor; deg_f = dmin;
647 154 : zeta = gen_1;
648 154 : prim = good_root_of_unity(p, deg_f, 1, &n);
649 154 : expo = diviuexact(subiu(p, 1), n);
650 : /* update X+a, avoid a=0 */
651 154 : gel (xplusa, 2) = addis (gel (xplusa, 2), 1);
652 154 : xpow = FpXQ_pow (xplusa, expo, fact, p);
653 329 : for (i = 0; i < n; i++)
654 : {
655 252 : GEN tmp = FpX_gcd(FpX_Fp_sub(xpow, zeta, p), fact, p);
656 252 : long dtmp = degpol(tmp);
657 252 : if (dtmp > 0 && dtmp < deg_f)
658 : {
659 161 : fact = FpX_div(fact, tmp, p); deg_f = degpol(fact);
660 161 : if (dtmp < dmin)
661 : {
662 154 : minfactor = FpX_normalize (tmp, p);
663 154 : dmin = dtmp;
664 154 : if (dmin == 1 || dmin <= (2 * deg_f) / n - 1)
665 : /* stop early to avoid too many gcds */
666 : break;
667 : }
668 : }
669 175 : zeta = Fp_mul (zeta, prim, p);
670 : }
671 : }
672 133 : return gerepileuptoint(av, FpX_oneroot(minfactor, p));
673 : }
674 :
675 : /*******************************************************************/
676 : /* */
677 : /* FACTORISATION MODULO p */
678 : /* */
679 : /*******************************************************************/
680 :
681 : /* F / E a vector of vectors of factors / exponents of virtual length l
682 : * (their real lg may be larger). Set their lg to j, concat and return [F,E] */
683 : static GEN
684 1587632 : FE_concat(GEN F, GEN E, long l)
685 : {
686 1587632 : setlg(E,l); E = shallowconcat1(E);
687 1587628 : setlg(F,l); F = shallowconcat1(F); return mkvec2(F,E);
688 : }
689 :
690 : static GEN
691 14 : ddf_to_ddf2_i(GEN V, long fl)
692 : {
693 : GEN F, D;
694 14 : long i, j, l = lg(V);
695 14 : F = cgetg(l, t_VEC);
696 14 : D = cgetg(l, t_VECSMALL);
697 112 : for (i = j = 1; i < l; i++)
698 : {
699 98 : GEN Vi = gel(V,i);
700 98 : if ((fl==2 && F2x_degree(Vi) == 0)
701 98 : ||(fl==0 && degpol(Vi) == 0)) continue;
702 35 : gel(F,j) = Vi;
703 35 : uel(D,j) = i; j++;
704 : }
705 14 : setlg(F,j);
706 14 : setlg(D,j); return mkvec2(F,D);
707 : }
708 :
709 : GEN
710 7 : ddf_to_ddf2(GEN V)
711 7 : { return ddf_to_ddf2_i(V, 0); }
712 :
713 : static GEN
714 7 : F2x_ddf_to_ddf2(GEN V)
715 7 : { return ddf_to_ddf2_i(V, 2); }
716 :
717 : GEN
718 5368745 : vddf_to_simplefact(GEN V, long d)
719 : {
720 : GEN E, F;
721 5368745 : long i, j, c, l = lg(V);
722 5368745 : F = cgetg(d+1, t_VECSMALL);
723 5368092 : E = cgetg(d+1, t_VECSMALL);
724 10821622 : for (i = c = 1; i < l; i++)
725 : {
726 5453037 : GEN Vi = gel(V,i);
727 5453037 : long l = lg(Vi);
728 27713413 : for (j = 1; j < l; j++)
729 : {
730 22259413 : long k, n = degpol(gel(Vi,j)) / j;
731 33619515 : for (k = 1; k <= n; k++) { uel(F,c) = j; uel(E,c) = i; c++; }
732 : }
733 : }
734 5368585 : setlg(F,c);
735 5368656 : setlg(E,c);
736 5368399 : return sort_factor(mkvec2(F,E), (void*)&cmpGuGu, cmp_nodata);
737 : }
738 :
739 : /* product of terms of degree 1 in factorization of f */
740 : GEN
741 182695 : FpX_split_part(GEN f, GEN p)
742 : {
743 182695 : long n = degpol(f);
744 182695 : GEN z, X = pol_x(varn(f));
745 182695 : if (n <= 1) return f;
746 180762 : f = FpX_red(f, p);
747 180761 : z = FpX_sub(FpX_Frobenius(f, p), X, p);
748 180762 : return FpX_gcd(z,f,p);
749 : }
750 :
751 : /* Compute the number of roots in Fp without counting multiplicity
752 : * return -1 for 0 polynomial. lc(f) must be prime to p. */
753 : long
754 99275 : FpX_nbroots(GEN f, GEN p)
755 : {
756 99275 : pari_sp av = avma;
757 99275 : GEN z = FpX_split_part(f, p);
758 99275 : return gc_long(av, degpol(z));
759 : }
760 :
761 : /* 1 < deg(f) <= p */
762 : static int
763 80569 : Flx_is_totally_split_i(GEN f, ulong p)
764 : {
765 80569 : GEN F = Flx_Frobenius(f, p);
766 80574 : return degpol(F)==1 && uel(F,2)==0UL && uel(F,3)==1UL;
767 : }
768 : int
769 80576 : Flx_is_totally_split(GEN f, ulong p)
770 : {
771 80576 : pari_sp av = avma;
772 80576 : ulong n = degpol(f);
773 80576 : if (n <= 1) return 1;
774 80569 : if (n > p) return 0; /* includes n < 0 */
775 80569 : return gc_bool(av, Flx_is_totally_split_i(f,p));
776 : }
777 : int
778 0 : FpX_is_totally_split(GEN f, GEN p)
779 : {
780 0 : pari_sp av = avma;
781 0 : ulong n = degpol(f);
782 : int u;
783 0 : if (n <= 1) return 1;
784 0 : if (abscmpui(n, p) > 0) return 0; /* includes n < 0 */
785 0 : if (lgefint(p) != 3)
786 0 : u = gequalX(FpX_Frobenius(FpX_red(f,p), p));
787 : else
788 : {
789 0 : ulong pp = (ulong)p[2];
790 0 : u = Flx_is_totally_split_i(ZX_to_Flx(f,pp), pp);
791 : }
792 0 : return gc_bool(av, u);
793 : }
794 :
795 : long
796 2702744 : Flx_nbroots(GEN f, ulong p)
797 : {
798 2702744 : long n = degpol(f);
799 2702745 : pari_sp av = avma;
800 : GEN z;
801 2702745 : if (n <= 1) return n;
802 2689263 : if (n == 2)
803 : {
804 : ulong D;
805 62461 : if (p==2) return (f[2]==0) + (f[2]!=f[3]);
806 61467 : D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
807 61467 : return 1 + krouu(D,p);
808 : }
809 2626802 : z = Flx_sub(Flx_Frobenius(f, p), polx_Flx(f[1]), p);
810 2626799 : z = Flx_gcd(z, f, p);
811 2626806 : return gc_long(av, degpol(z));
812 : }
813 :
814 : long
815 4403 : FpX_ddf_degree(GEN T, GEN XP, GEN p)
816 : {
817 4403 : pari_sp av = avma;
818 : GEN X, b, g, xq;
819 : long i, j, n, v, B, l, m;
820 : pari_timer ti;
821 : hashtable h;
822 :
823 4403 : n = get_FpX_degree(T); v = get_FpX_var(T);
824 4403 : X = pol_x(v);
825 4403 : if (ZX_equal(X,XP)) return 1;
826 4403 : B = n/2;
827 4403 : l = usqrt(B);
828 4403 : m = (B+l-1)/l;
829 4403 : T = FpX_get_red(T, p);
830 4403 : hash_init_GEN(&h, l+2, ZX_equal, 1);
831 4403 : hash_insert_long(&h, X, 0);
832 4403 : hash_insert_long(&h, XP, 1);
833 4403 : if (DEBUGLEVEL>=7) timer_start(&ti);
834 4403 : b = XP;
835 4403 : xq = FpXQ_powers(b, brent_kung_optpow(n, l-1, 1), T, p);
836 4403 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq baby");
837 10402 : for (i = 3; i <= l+1; i++)
838 : {
839 6790 : b = FpX_FpXQV_eval(b, xq, T, p);
840 6790 : if (gequalX(b)) return gc_long(av,i-1);
841 5999 : hash_insert_long(&h, b, i-1);
842 : }
843 3612 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: baby");
844 3612 : g = b;
845 3612 : xq = FpXQ_powers(g, brent_kung_optpow(n, m, 1), T, p);
846 3612 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq giant");
847 12341 : for(i = 2; i <= m+1; i++)
848 : {
849 10745 : g = FpX_FpXQV_eval(g, xq, T, p);
850 10745 : if (hash_haskey_long(&h, g, &j)) return gc_long(av, l*i-j);
851 : }
852 1596 : return gc_long(av,n);
853 : }
854 :
855 : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
856 : static GEN
857 911 : FpX_ddf_Shoup(GEN T, GEN XP, GEN p)
858 : {
859 : GEN b, g, h, F, f, Tr, xq;
860 : long i, j, n, v, B, l, m;
861 : pari_timer ti;
862 :
863 911 : n = get_FpX_degree(T); v = get_FpX_var(T);
864 911 : if (n == 0) return cgetg(1, t_VEC);
865 911 : if (n == 1) return mkvec(get_FpX_mod(T));
866 805 : B = n/2;
867 805 : l = usqrt(B);
868 805 : m = (B+l-1)/l;
869 805 : T = FpX_get_red(T, p);
870 805 : b = cgetg(l+2, t_VEC);
871 805 : gel(b, 1) = pol_x(v);
872 805 : gel(b, 2) = XP;
873 805 : if (DEBUGLEVEL>=7) timer_start(&ti);
874 805 : xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1), T, p);
875 805 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq baby");
876 1065 : for (i = 3; i <= l+1; i++)
877 260 : gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
878 805 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: baby");
879 805 : xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1), T, p);
880 805 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq giant");
881 805 : g = cgetg(m+1, t_VEC);
882 805 : gel(g, 1) = gel(xq, 2);
883 1804 : for(i = 2; i <= m; i++) gel(g, i) = FpX_FpXQV_eval(gel(g, i-1), xq, T, p);
884 805 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: giant");
885 805 : h = cgetg(m+1, t_VEC);
886 2609 : for (j = 1; j <= m; j++)
887 : {
888 1804 : pari_sp av = avma;
889 1804 : GEN gj = gel(g,j), e = FpX_sub(gj, gel(b,1), p);
890 3175 : for (i = 2; i <= l; i++) e = FpXQ_mul(e, FpX_sub(gj, gel(b,i), p), T, p);
891 1804 : gel(h,j) = gerepileupto(av, e);
892 : }
893 805 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: diff");
894 805 : Tr = get_FpX_mod(T);
895 805 : F = cgetg(m+1, t_VEC);
896 2609 : for (j = 1; j <= m; j++)
897 : {
898 1804 : GEN u = FpX_gcd(Tr, gel(h,j), p);
899 1804 : if (degpol(u))
900 : {
901 409 : u = FpX_normalize(u, p);
902 409 : Tr = FpX_div(Tr, u, p);
903 : }
904 1804 : gel(F,j) = u;
905 : }
906 805 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: F");
907 805 : f = const_vec(n, pol_1(v));
908 2609 : for (j = 1; j <= m; j++)
909 : {
910 1804 : GEN e = gel(F, j);
911 1861 : for (i=l-1; i >= 0; i--)
912 : {
913 1861 : GEN u = FpX_gcd(e, FpX_sub(gel(g, j), gel(b, i+1), p), p);
914 1861 : if (degpol(u))
915 : {
916 417 : u = FpX_normalize(u, p);
917 417 : gel(f, l*j-i) = u;
918 417 : e = FpX_div(e, u, p);
919 : }
920 1861 : if (!degpol(e)) break;
921 : }
922 : }
923 805 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: f");
924 805 : if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
925 805 : return f;
926 : }
927 :
928 : static void
929 0 : FpX_edf_simple(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
930 : {
931 0 : long n = degpol(Tp), r = n/d, ct = 0;
932 : GEN T, f, ff, p2;
933 0 : if (r==1) { gel(V, idx) = Tp; return; }
934 0 : p2 = shifti(p,-1);
935 0 : T = FpX_get_red(Tp, p);
936 0 : XP = FpX_rem(XP, T, p);
937 : while (1)
938 0 : {
939 0 : pari_sp btop = avma;
940 : long i;
941 0 : GEN g = random_FpX(n, varn(Tp), p);
942 0 : GEN t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
943 0 : if (signe(t) == 0) continue;
944 0 : for(i=1; i<=10; i++)
945 : {
946 0 : pari_sp btop2 = avma;
947 0 : GEN R = FpXQ_pow(FpX_Fp_add(t, randomi(p), p), p2, T, p);
948 0 : f = FpX_gcd(FpX_Fp_sub(R, gen_1, p), Tp, p);
949 0 : if (degpol(f) > 0 && degpol(f) < n) break;
950 0 : set_avma(btop2);
951 : }
952 0 : if (degpol(f) > 0 && degpol(f) < n) break;
953 0 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_simple",p);
954 0 : set_avma(btop);
955 : }
956 0 : f = FpX_normalize(f, p);
957 0 : ff = FpX_div(Tp, f ,p);
958 0 : FpX_edf_simple(f, XP, d, p, V, idx);
959 0 : FpX_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
960 : }
961 :
962 : static void
963 1007 : FpX_edf_rec(GEN T, GEN hp, GEN t, long d, GEN p2, GEN p, GEN V, long idx)
964 : {
965 : pari_sp av;
966 1007 : GEN Tp = get_FpX_mod(T);
967 1007 : long n = degpol(hp), vT = varn(Tp), ct = 0;
968 : GEN u1, u2, f1, f2, R, h;
969 1007 : h = FpX_get_red(hp, p);
970 1007 : t = FpX_rem(t, T, p);
971 1007 : av = avma;
972 : do
973 : {
974 1682 : set_avma(av);
975 1682 : R = FpXQ_pow(deg1pol(gen_1, randomi(p), vT), p2, h, p);
976 1682 : u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
977 1682 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_rec",p);
978 1682 : } while (degpol(u1)==0 || degpol(u1)==n);
979 1007 : f1 = FpX_gcd(FpX_FpXQ_eval(u1, t, T, p), Tp, p);
980 1007 : f1 = FpX_normalize(f1, p);
981 1007 : u2 = FpX_div(hp, u1, p);
982 1007 : f2 = FpX_div(Tp, f1, p);
983 1007 : if (degpol(u1)==1)
984 669 : gel(V, idx) = f1;
985 : else
986 338 : FpX_edf_rec(FpX_get_red(f1, p), u1, t, d, p2, p, V, idx);
987 1007 : idx += degpol(f1)/d;
988 1007 : if (degpol(u2)==1)
989 730 : gel(V, idx) = f2;
990 : else
991 277 : FpX_edf_rec(FpX_get_red(f2, p), u2, t, d, p2, p, V, idx);
992 1007 : }
993 :
994 : /* assume Tp a squarefree product of r > 1 irred. factors of degree d */
995 : static void
996 392 : FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
997 : {
998 392 : long n = degpol(Tp), r = n/d, vT = varn(Tp), ct = 0;
999 : GEN T, h, t;
1000 : pari_timer ti;
1001 :
1002 392 : T = FpX_get_red(Tp, p);
1003 392 : XP = FpX_rem(XP, T, p);
1004 392 : if (DEBUGLEVEL>=7) timer_start(&ti);
1005 : do
1006 : {
1007 392 : GEN g = random_FpX(n, vT, p);
1008 392 : t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
1009 392 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
1010 392 : h = FpXQ_minpoly(t, T, p);
1011 392 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
1012 392 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf",p);
1013 392 : } while (degpol(h) != r);
1014 392 : FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
1015 392 : }
1016 :
1017 : static GEN
1018 897 : FpX_factor_Shoup(GEN T, GEN p)
1019 : {
1020 897 : long i, n, s = 0;
1021 : GEN XP, D, V;
1022 897 : long e = expi(p);
1023 : pari_timer ti;
1024 897 : n = get_FpX_degree(T);
1025 897 : T = FpX_get_red(T, p);
1026 897 : if (DEBUGLEVEL>=6) timer_start(&ti);
1027 897 : XP = FpX_Frobenius(T, p);
1028 897 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
1029 897 : D = FpX_ddf_Shoup(T, XP, p);
1030 897 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
1031 897 : s = ddf_to_nbfact(D);
1032 897 : V = cgetg(s+1, t_COL);
1033 7279 : for (i = 1, s = 1; i <= n; i++)
1034 : {
1035 6382 : GEN Di = gel(D,i);
1036 6382 : long ni = degpol(Di), ri = ni/i;
1037 6382 : if (ni == 0) continue;
1038 930 : Di = FpX_normalize(Di, p);
1039 930 : if (ni == i) { gel(V, s++) = Di; continue; }
1040 392 : if (ri <= e*expu(e))
1041 392 : FpX_edf(Di, XP, i, p, V, s);
1042 : else
1043 0 : FpX_edf_simple(Di, XP, i, p, V, s);
1044 392 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
1045 392 : s += ri;
1046 : }
1047 897 : return V;
1048 : }
1049 :
1050 : long
1051 2047063 : ddf_to_nbfact(GEN D)
1052 : {
1053 2047063 : long l = lg(D), i, s = 0;
1054 13733030 : for(i = 1; i < l; i++) s += degpol(gel(D,i))/i;
1055 2047064 : return s;
1056 : }
1057 :
1058 : /* Yun algorithm: Assume p > degpol(T) */
1059 : static GEN
1060 1595 : FpX_factor_Yun(GEN T, GEN p)
1061 : {
1062 1595 : long n = degpol(T), i = 1;
1063 1595 : GEN a, b, c, d = FpX_deriv(T, p);
1064 1595 : GEN V = cgetg(n+1,t_VEC);
1065 1595 : a = FpX_gcd(T, d, p);
1066 1595 : if (degpol(a) == 0) return mkvec(T);
1067 535 : b = FpX_div(T, a, p);
1068 : do
1069 : {
1070 2526 : c = FpX_div(d, a, p);
1071 2526 : d = FpX_sub(c, FpX_deriv(b, p), p);
1072 2526 : a = FpX_normalize(FpX_gcd(b, d, p), p);
1073 2526 : gel(V, i++) = a;
1074 2526 : b = FpX_div(b, a, p);
1075 2526 : } while (degpol(b));
1076 535 : setlg(V, i); return V;
1077 : }
1078 : GEN
1079 333676 : FpX_factor_squarefree(GEN T, GEN p)
1080 : {
1081 333676 : if (lgefint(p)==3)
1082 : {
1083 332960 : ulong pp = (ulong)p[2];
1084 332960 : GEN u = Flx_factor_squarefree(ZX_to_Flx(T,pp), pp);
1085 332960 : return FlxV_to_ZXV(u);
1086 : }
1087 716 : return FpX_factor_Yun(T, p);
1088 : }
1089 :
1090 : GEN
1091 218470 : FpX_roots_mult(GEN T, long n, GEN p)
1092 : {
1093 218470 : GEN V = FpX_factor_squarefree(T,p), W;
1094 218470 : long l = lg(V), i;
1095 218470 : if (l<=n) return cgetg(1,t_COL);
1096 35455 : W = cgetg(l-n+1,t_VEC);
1097 116473 : for (i = n; i < l; i++)
1098 81018 : gel(W,i-n+1) = FpX_roots(gel(V,i), p);
1099 35455 : return shallowconcat1(W);
1100 : }
1101 :
1102 : long
1103 168 : FpX_ispower(GEN f, ulong k, GEN p, GEN *pt_r)
1104 : {
1105 168 : pari_sp av = avma;
1106 : GEN lc, F;
1107 168 : long i, l, n = degpol(f), v = varn(f);
1108 168 : if (n % k) return 0;
1109 168 : if (lgefint(p)==3)
1110 : {
1111 126 : ulong pp = p[2];
1112 126 : GEN fp = ZX_to_Flx(f, pp);
1113 126 : if (!Flx_ispower(fp, k, pp, pt_r)) return gc_long(av,0);
1114 105 : if (pt_r) *pt_r = gerepileupto(av, Flx_to_ZX(*pt_r)); else set_avma(av);
1115 105 : return 1;
1116 : }
1117 42 : lc = Fp_sqrtn(leading_coeff(f), stoi(k), p, NULL);
1118 42 : if (!lc) { av = avma; return 0; }
1119 42 : F = FpX_factor_Yun(f, p); l = lg(F)-1;
1120 1491 : for(i=1; i <= l; i++)
1121 1456 : if (i%k && degpol(gel(F,i))) return gc_long(av,0);
1122 35 : if (pt_r)
1123 : {
1124 35 : GEN r = scalarpol(lc, v), s = pol_1(v);
1125 1484 : for (i=l; i>=1; i--)
1126 : {
1127 1449 : if (i%k) continue;
1128 294 : s = FpX_mul(s, gel(F,i), p);
1129 294 : r = FpX_mul(r, s, p);
1130 : }
1131 35 : *pt_r = gerepileupto(av, r);
1132 0 : } else av = avma;
1133 35 : return 1;
1134 : }
1135 :
1136 : static GEN
1137 823 : FpX_factor_Cantor(GEN T, GEN p)
1138 : {
1139 823 : GEN E, F, V = FpX_factor_Yun(T, p);
1140 823 : long i, j, l = lg(V);
1141 823 : F = cgetg(l, t_VEC);
1142 823 : E = cgetg(l, t_VEC);
1143 2197 : for (i=1, j=1; i < l; i++)
1144 1374 : if (degpol(gel(V,i)))
1145 : {
1146 897 : GEN Fj = FpX_factor_Shoup(gel(V,i), p);
1147 897 : gel(F, j) = Fj;
1148 897 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1149 897 : j++;
1150 : }
1151 823 : return sort_factor_pol(FE_concat(F,E,j), cmpii);
1152 : }
1153 :
1154 : static GEN
1155 0 : FpX_ddf_i(GEN T, GEN p)
1156 : {
1157 : GEN XP;
1158 0 : T = FpX_get_red(T, p);
1159 0 : XP = FpX_Frobenius(T, p);
1160 0 : return ddf_to_ddf2(FpX_ddf_Shoup(T, XP, p));
1161 : }
1162 :
1163 : GEN
1164 7 : FpX_ddf(GEN f, GEN p)
1165 : {
1166 7 : pari_sp av = avma;
1167 : GEN F;
1168 7 : switch(ZX_factmod_init(&f, p))
1169 : {
1170 7 : case 0: F = F2x_ddf(f);
1171 7 : F2xV_to_ZXV_inplace(gel(F,1)); break;
1172 0 : case 1: F = Flx_ddf(f,p[2]);
1173 0 : FlxV_to_ZXV_inplace(gel(F,1)); break;
1174 0 : default: F = FpX_ddf_i(f,p); break;
1175 : }
1176 7 : return gerepilecopy(av, F);
1177 : }
1178 :
1179 : static GEN Flx_simplefact_Cantor(GEN T, ulong p);
1180 : static GEN
1181 14 : FpX_simplefact_Cantor(GEN T, GEN p)
1182 : {
1183 : GEN V;
1184 : long i, l;
1185 14 : if (lgefint(p) == 3)
1186 : {
1187 0 : ulong pp = p[2];
1188 0 : return Flx_simplefact_Cantor(ZX_to_Flx(T,pp), pp);
1189 : }
1190 14 : T = FpX_get_red(T, p);
1191 14 : V = FpX_factor_Yun(get_FpX_mod(T), p); l = lg(V);
1192 28 : for (i=1; i < l; i++)
1193 14 : gel(V,i) = FpX_ddf_Shoup(gel(V,i), FpX_Frobenius(gel(V,i), p), p);
1194 14 : return vddf_to_simplefact(V, get_FpX_degree(T));
1195 : }
1196 :
1197 : static int
1198 0 : FpX_isirred_Cantor(GEN Tp, GEN p)
1199 : {
1200 0 : pari_sp av = avma;
1201 : pari_timer ti;
1202 : long n;
1203 0 : GEN T = get_FpX_mod(Tp);
1204 0 : GEN dT = FpX_deriv(T, p);
1205 : GEN XP, D;
1206 0 : if (degpol(FpX_gcd(T, dT, p)) != 0) return gc_bool(av,0);
1207 0 : n = get_FpX_degree(T);
1208 0 : T = FpX_get_red(Tp, p);
1209 0 : if (DEBUGLEVEL>=6) timer_start(&ti);
1210 0 : XP = FpX_Frobenius(T, p);
1211 0 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
1212 0 : D = FpX_ddf_Shoup(T, XP, p);
1213 0 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
1214 0 : return gc_bool(av, degpol(gel(D,n)) == n);
1215 : }
1216 :
1217 : static GEN FpX_factor_deg2(GEN f, GEN p, long d, long flag);
1218 :
1219 : /*Assume that p is large and odd*/
1220 : static GEN
1221 1695 : FpX_factor_i(GEN f, GEN pp, long flag)
1222 : {
1223 1695 : long d = degpol(f);
1224 1695 : if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
1225 837 : switch(flag)
1226 : {
1227 823 : default: return FpX_factor_Cantor(f, pp);
1228 14 : case 1: return FpX_simplefact_Cantor(f, pp);
1229 0 : case 2: return FpX_isirred_Cantor(f, pp)? gen_1: NULL;
1230 : }
1231 : }
1232 :
1233 : long
1234 0 : FpX_nbfact_Frobenius(GEN T, GEN XP, GEN p)
1235 : {
1236 0 : pari_sp av = avma;
1237 0 : long s = ddf_to_nbfact(FpX_ddf_Shoup(T, XP, p));
1238 0 : return gc_long(av,s);
1239 : }
1240 :
1241 : long
1242 0 : FpX_nbfact(GEN T, GEN p)
1243 : {
1244 0 : pari_sp av = avma;
1245 0 : GEN XP = FpX_Frobenius(T, p);
1246 0 : long n = FpX_nbfact_Frobenius(T, XP, p);
1247 0 : return gc_long(av,n);
1248 : }
1249 :
1250 : /* p > 2 */
1251 : static GEN
1252 7 : FpX_is_irred_2(GEN f, GEN p, long d)
1253 : {
1254 7 : switch(d)
1255 : {
1256 0 : case -1:
1257 0 : case 0: return NULL;
1258 0 : case 1: return gen_1;
1259 : }
1260 7 : return FpX_quad_factortype(f, p) == -1? gen_1: NULL;
1261 : }
1262 : /* p > 2 */
1263 : static GEN
1264 14 : FpX_degfact_2(GEN f, GEN p, long d)
1265 : {
1266 14 : switch(d)
1267 : {
1268 0 : case -1:retmkvec2(mkvecsmall(-1),mkvecsmall(1));
1269 0 : case 0: return trivial_fact();
1270 0 : case 1: retmkvec2(mkvecsmall(1), mkvecsmall(1));
1271 : }
1272 14 : switch(FpX_quad_factortype(f, p)) {
1273 7 : case 1: retmkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1274 7 : case -1: retmkvec2(mkvecsmall(2), mkvecsmall(1));
1275 0 : default: retmkvec2(mkvecsmall(1), mkvecsmall(2));
1276 : }
1277 : }
1278 :
1279 : GEN
1280 70 : prime_fact(GEN x) { retmkmat2(mkcolcopy(x), mkcol(gen_1)); }
1281 : GEN
1282 164419 : trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
1283 :
1284 : /* not gerepile safe */
1285 : static GEN
1286 837 : FpX_factor_2(GEN f, GEN p, long d)
1287 : {
1288 : GEN r, s, R, S;
1289 : long v;
1290 : int sgn;
1291 837 : switch(d)
1292 : {
1293 7 : case -1: retmkvec2(mkcol(pol_0(varn(f))), mkvecsmall(1));
1294 30 : case 0: retmkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1295 78 : case 1: retmkvec2(mkcol(f), mkvecsmall(1));
1296 : }
1297 722 : r = FpX_quad_root(f, p, 1);
1298 722 : if (!r) return mkvec2(mkcol(f), mkvecsmall(1));
1299 187 : v = varn(f);
1300 187 : s = FpX_otherroot(f, r, p);
1301 187 : if (signe(r)) r = subii(p, r);
1302 187 : if (signe(s)) s = subii(p, s);
1303 187 : sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
1304 187 : R = deg1pol_shallow(gen_1, r, v);
1305 187 : if (!sgn) return mkvec2(mkcol(R), mkvecsmall(2));
1306 93 : S = deg1pol_shallow(gen_1, s, v);
1307 93 : return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
1308 : }
1309 : static GEN
1310 858 : FpX_factor_deg2(GEN f, GEN p, long d, long flag)
1311 : {
1312 858 : switch(flag) {
1313 7 : case 2: return FpX_is_irred_2(f, p, d);
1314 14 : case 1: return FpX_degfact_2(f, p, d);
1315 837 : default: return FpX_factor_2(f, p, d);
1316 : }
1317 : }
1318 :
1319 : static int
1320 450366 : F2x_quad_factortype(GEN x)
1321 450366 : { return x[2] == 7 ? -1: x[2] == 6 ? 1 :0; }
1322 :
1323 : static GEN
1324 14 : F2x_is_irred_2(GEN f, long d)
1325 14 : { return d == 1 || (d==2 && F2x_quad_factortype(f) == -1)? gen_1: NULL; }
1326 :
1327 : static GEN
1328 19593 : F2x_degfact_2(GEN f, long d)
1329 : {
1330 19593 : if (!d) return trivial_fact();
1331 19593 : if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1332 19404 : switch(F2x_quad_factortype(f)) {
1333 5964 : case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1334 6202 : case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1335 7238 : default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1336 : }
1337 : }
1338 :
1339 : static GEN
1340 1165906 : F2x_factor_2(GEN f, long d)
1341 : {
1342 1165906 : long v = f[1];
1343 1165906 : if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1344 909186 : if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1345 426161 : switch(F2x_quad_factortype(f))
1346 : {
1347 72507 : case -1: return mkvec2(mkcol(f), mkvecsmall(1));
1348 326235 : case 0: return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
1349 27438 : default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
1350 : }
1351 : }
1352 : static GEN
1353 1185515 : F2x_factor_deg2(GEN f, long d, long flag)
1354 : {
1355 1185515 : switch(flag) {
1356 14 : case 2: return F2x_is_irred_2(f, d);
1357 19593 : case 1: return F2x_degfact_2(f, d);
1358 1165908 : default: return F2x_factor_2(f, d);
1359 : }
1360 : }
1361 :
1362 : /* xt = NULL or x^(p-1)/2 mod g */
1363 : static void
1364 27725 : split_squares(struct split_t *S, GEN g, ulong p, GEN xt)
1365 : {
1366 27725 : ulong q = p >> 1;
1367 27725 : GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
1368 27725 : long d = degpol(a);
1369 27725 : if (d < 0)
1370 : {
1371 : ulong i;
1372 611 : split_add_done(S, (GEN)1);
1373 2454 : for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr(i,p));
1374 : } else {
1375 27114 : if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
1376 27114 : if (d)
1377 : {
1378 27079 : if (xt) xt = Flx_Fl_add(xt, p-1, p); else xt = Flx_Xnm1(g[1], q, p);
1379 27079 : a = Flx_gcd(a, xt, p);
1380 27079 : if (degpol(a)) split_add(S, Flx_normalize(a, p));
1381 : }
1382 : }
1383 27725 : }
1384 : static void
1385 27725 : split_nonsquares(struct split_t *S, GEN g, ulong p, GEN xt)
1386 : {
1387 27725 : ulong q = p >> 1;
1388 27725 : GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
1389 27725 : long d = degpol(a);
1390 27725 : if (d < 0)
1391 : {
1392 167 : ulong i, z = nonsquare_Fl(p);
1393 167 : split_add_done(S, (GEN)z);
1394 306 : for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_mul(z, Fl_sqr(i,p), p));
1395 : } else {
1396 27558 : if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
1397 27558 : if (d)
1398 : {
1399 27026 : if (xt) xt = Flx_Fl_add(xt, 1, p); else xt = Flx_Xn1(g[1], q, p);
1400 27026 : a = Flx_gcd(a, xt, p);
1401 27025 : if (degpol(a)) split_add(S, Flx_normalize(a, p));
1402 : }
1403 : }
1404 27725 : }
1405 : /* p > 2. f monic Flx, f(0) != 0. Add to split_t structs coprime factors
1406 : * of g = \prod_{f(a) = 0} (X - a). Return 0 when f(x) = 0 for all x in Fp* */
1407 : static int
1408 5018248 : split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p)
1409 : {
1410 5018248 : GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
1411 5018060 : long d = degpol(g);
1412 5017955 : if (d < 0) return 0;
1413 5017521 : if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
1414 5017533 : if (!d) return 1;
1415 4989018 : if ((p >> 4) <= (ulong)d)
1416 : { /* small p; split directly using x^((p-1)/2) +/- 1 */
1417 23891 : GEN xt = ((ulong)d < (p>>1))? Flx_rem(monomial_Flx(1, p>>1, g[1]), g, p)
1418 27725 : : NULL;
1419 27725 : split_squares(S, g, p, xt);
1420 27725 : split_nonsquares(S, g, p, xt);
1421 : } else { /* large p; use x^(p-1) - 1 directly */
1422 4961293 : a = Flxq_powu(polx_Flx(f[1]), p-1, g,p);
1423 4947902 : if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
1424 4947902 : a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
1425 4949997 : g = Flx_gcd(g,a, p);
1426 4954920 : if (degpol(g)) split_add(S, Flx_normalize(g,p));
1427 : }
1428 4987814 : return 1;
1429 : }
1430 :
1431 : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
1432 : static GEN
1433 23733161 : Flx_roots_i(GEN f, ulong p)
1434 : {
1435 : GEN pol, g;
1436 23733161 : long v = Flx_valrem(f, &g);
1437 : ulong q;
1438 : struct split_t S;
1439 :
1440 : /* optimization: test for small degree first */
1441 23744146 : switch(degpol(g))
1442 : {
1443 99469 : case 1: {
1444 99469 : ulong r = p - g[2];
1445 99469 : return v? mkvecsmall2(0, r): mkvecsmall(r);
1446 : }
1447 18614932 : case 2: {
1448 18614932 : ulong r = Flx_quad_root(g, p, 1), s;
1449 18604954 : if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
1450 12685477 : s = Flx_otherroot(g,r, p);
1451 12771476 : if (r < s)
1452 3196349 : return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
1453 9575127 : else if (r > s)
1454 9598952 : return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
1455 : else
1456 4554 : return v? mkvecsmall2(0, s): mkvecsmall(s);
1457 : }
1458 : }
1459 5019275 : q = p >> 1;
1460 5019275 : split_init(&S, lg(f)-1);
1461 5018253 : settyp(S.done, t_VECSMALL);
1462 5018253 : if (v) split_add_done(&S, (GEN)0);
1463 5018253 : if (! split_Flx_cut_out_roots(&S, g, p))
1464 434 : return all_roots_mod_p(p, lg(S.done) == 1);
1465 5016206 : pol = polx_Flx(f[1]);
1466 5017061 : for (pol[2]=1; ; pol[2]++)
1467 5371602 : {
1468 10388663 : long j, l = lg(S.todo);
1469 10388663 : if (l == 1) { vecsmall_sort(S.done); return S.done; }
1470 5371324 : if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
1471 11224779 : for (j = 1; j < l; j++)
1472 : {
1473 5853177 : GEN b, c = gel(S.todo,j);
1474 : ulong r, s;
1475 5853177 : switch(degpol(c))
1476 : {
1477 4887282 : case 1:
1478 4887282 : split_moveto_done(&S, j, (GEN)(p - c[2]));
1479 4887712 : j--; l--; break;
1480 447988 : case 2:
1481 447988 : r = Flx_quad_root(c, p, 0);
1482 447995 : if (r == p) pari_err_PRIME("polrootsmod",utoipos(p));
1483 447988 : s = Flx_otherroot(c,r, p);
1484 447989 : split_done(&S, j, (GEN)r, (GEN)s);
1485 447996 : j--; l--; break;
1486 517799 : default:
1487 517799 : b = Flxq_powu(pol,q, c,p); /* pol^(p-1)/2 */
1488 517895 : if (degpol(b) <= 0) continue;
1489 405132 : b = Flx_gcd(c,Flx_Fl_add(b,p-1,p), p);
1490 405148 : if (!degpol(b)) continue;
1491 404683 : b = Flx_normalize(b, p);
1492 404692 : c = Flx_div(c,b, p);
1493 404681 : split_todo(&S, j, b, c);
1494 : }
1495 : }
1496 : }
1497 : }
1498 :
1499 : GEN
1500 23526617 : Flx_roots(GEN f, ulong p)
1501 : {
1502 23526617 : pari_sp av = avma;
1503 23526617 : switch(lg(f))
1504 : {
1505 0 : case 2: pari_err_ROOTS0("Flx_roots");
1506 0 : case 3: set_avma(av); return cgetg(1, t_VECSMALL);
1507 : }
1508 23527619 : if (p == 2) return Flx_root_mod_2(f);
1509 23527612 : return gerepileuptoleaf(av, Flx_roots_i(Flx_normalize(f, p), p));
1510 : }
1511 :
1512 : /* assume x reduced mod p, monic. */
1513 : static int
1514 2467329 : Flx_quad_factortype(GEN x, ulong p)
1515 : {
1516 2467329 : ulong b = x[3], c = x[2];
1517 2467329 : return krouu(Fl_disc_bc(b, c, p), p);
1518 : }
1519 : static GEN
1520 56 : Flx_is_irred_2(GEN f, ulong p, long d)
1521 : {
1522 56 : if (!d) return NULL;
1523 56 : if (d == 1) return gen_1;
1524 56 : return Flx_quad_factortype(f, p) == -1? gen_1: NULL;
1525 : }
1526 : static GEN
1527 2491533 : Flx_degfact_2(GEN f, ulong p, long d)
1528 : {
1529 2491533 : if (!d) return trivial_fact();
1530 2491533 : if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1531 2467439 : switch(Flx_quad_factortype(f, p)) {
1532 1171057 : case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1533 1269130 : case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1534 27829 : default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1535 : }
1536 : }
1537 : /* p > 2 */
1538 : static GEN
1539 1221206 : Flx_factor_2(GEN f, ulong p, long d)
1540 : {
1541 : ulong r, s;
1542 : GEN R,S;
1543 1221206 : long v = f[1];
1544 1221206 : if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1545 1144764 : if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1546 916507 : r = Flx_quad_root(f, p, 1);
1547 916537 : if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
1548 615314 : s = Flx_otherroot(f, r, p);
1549 615328 : r = Fl_neg(r, p);
1550 615326 : s = Fl_neg(s, p);
1551 615324 : if (s < r) lswap(s,r);
1552 615324 : R = mkvecsmall3(v,r,1);
1553 615323 : if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
1554 506424 : S = mkvecsmall3(v,s,1);
1555 506428 : return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
1556 : }
1557 : static GEN
1558 3712084 : Flx_factor_deg2(GEN f, ulong p, long d, long flag)
1559 : {
1560 3712084 : switch(flag) {
1561 56 : case 2: return Flx_is_irred_2(f, p, d);
1562 2491619 : case 1: return Flx_degfact_2(f, p, d);
1563 1220409 : default: return Flx_factor_2(f, p, d);
1564 : }
1565 : }
1566 :
1567 : static GEN
1568 23359 : F2x_Berlekamp_ker(GEN u)
1569 : {
1570 23359 : pari_sp ltop=avma;
1571 23359 : long j,N = F2x_degree(u);
1572 : GEN Q;
1573 : pari_timer T;
1574 23356 : timer_start(&T);
1575 23359 : Q = F2x_matFrobenius(u);
1576 322302 : for (j=1; j<=N; j++)
1577 298943 : F2m_flip(Q,j,j);
1578 23359 : if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
1579 23359 : Q = F2m_ker_sp(Q,0);
1580 23359 : if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
1581 23359 : return gerepileupto(ltop,Q);
1582 : }
1583 : #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
1584 : static long
1585 35947 : F2x_split_Berlekamp(GEN *t)
1586 : {
1587 35947 : GEN u = *t, a, b, vker;
1588 35947 : long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
1589 :
1590 35947 : if (du == 1) return 1;
1591 27425 : if (du == 2)
1592 : {
1593 4066 : if (F2x_quad_factortype(u) == 1) /* 0 is a root: shouldn't occur */
1594 : {
1595 0 : t[0] = mkvecsmall2(sv, 2);
1596 0 : t[1] = mkvecsmall2(sv, 3);
1597 0 : return 2;
1598 : }
1599 4066 : return 1;
1600 : }
1601 :
1602 23359 : vker = F2x_Berlekamp_ker(u);
1603 23359 : lb = lgcols(vker);
1604 23377 : d = lg(vker)-1;
1605 23377 : ir = 0;
1606 : /* t[i] irreducible for i < ir, still to be treated for i < L */
1607 60811 : for (L=1; L<d; )
1608 : {
1609 : GEN pol;
1610 37454 : if (d == 2)
1611 5616 : pol = F2v_to_F2x(gel(vker,2), sv);
1612 : else
1613 : {
1614 31838 : GEN v = zero_zv(lb);
1615 31838 : v[1] = du;
1616 31838 : v[2] = random_Fl(2); /*Assume vker[1]=1*/
1617 119792 : for (i=2; i<=d; i++)
1618 87953 : if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
1619 31839 : pol = F2v_to_F2x(v, sv);
1620 : }
1621 112614 : for (i=ir; i<L && L<d; i++)
1622 : {
1623 75180 : a = t[i]; la = F2x_degree(a);
1624 75174 : if (la == 1) { set_irred(i); }
1625 74996 : else if (la == 2)
1626 : {
1627 709 : if (F2x_quad_factortype(a) == 1) /* 0 is a root: shouldn't occur */
1628 : {
1629 0 : t[i] = mkvecsmall2(sv, 2);
1630 0 : t[L] = mkvecsmall2(sv, 3); L++;
1631 : }
1632 709 : set_irred(i);
1633 : }
1634 : else
1635 : {
1636 74287 : pari_sp av = avma;
1637 : long lb;
1638 74287 : b = F2x_rem(pol, a);
1639 74278 : if (F2x_degree(b) <= 0) { set_avma(av); continue; }
1640 27053 : b = F2x_gcd(a,b); lb = F2x_degree(b);
1641 27060 : if (lb && lb < la)
1642 : {
1643 27060 : t[L] = F2x_div(a,b);
1644 27058 : t[i]= b; L++;
1645 : }
1646 0 : else set_avma(av);
1647 : }
1648 : }
1649 : }
1650 23357 : return d;
1651 : }
1652 : /* assume deg f > 2 */
1653 : static GEN
1654 35306 : F2x_Berlekamp_i(GEN f, long flag)
1655 : {
1656 35306 : long lfact, val, d = F2x_degree(f), j, k, lV;
1657 : GEN y, E, t, V;
1658 :
1659 35306 : val = F2x_valrem(f, &f);
1660 35306 : if (flag == 2 && val) return NULL;
1661 35292 : V = F2x_factor_squarefree(f); lV = lg(V);
1662 35291 : if (flag == 2 && lV > 2) return NULL;
1663 :
1664 : /* to hold factors and exponents */
1665 35221 : t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
1666 35223 : E = cgetg(d+1,t_VECSMALL);
1667 35223 : lfact = 1;
1668 35223 : if (val) {
1669 11272 : if (flag == 1) t[1] = 1; else gel(t,1) = polx_F2x(f[1]);
1670 11272 : E[1] = val; lfact++;
1671 : }
1672 :
1673 116286 : for (k=1; k<lV; k++)
1674 : {
1675 81219 : if (F2x_degree(gel(V, k))==0) continue;
1676 35947 : gel(t,lfact) = gel(V, k);
1677 35947 : d = F2x_split_Berlekamp(&gel(t,lfact));
1678 35945 : if (flag == 2 && d != 1) return NULL;
1679 35791 : if (flag == 1)
1680 48066 : for (j=0; j<d; j++) t[lfact+j] = F2x_degree(gel(t,lfact+j));
1681 98231 : for (j=0; j<d; j++) E[lfact+j] = k;
1682 35791 : lfact += d;
1683 : }
1684 35067 : if (flag == 2) return gen_1; /* irreducible */
1685 35053 : setlg(t, lfact);
1686 35054 : setlg(E, lfact); y = mkvec2(t,E);
1687 24352 : return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
1688 59406 : : sort_factor_pol(y, cmpGuGu);
1689 : }
1690 :
1691 : /* Adapted from Shoup NTL */
1692 : GEN
1693 281678 : F2x_factor_squarefree(GEN f)
1694 : {
1695 : GEN r, t, v, tv;
1696 281678 : long i, q, n = F2x_degree(f);
1697 281677 : GEN u = const_vec(n+1, pol1_F2x(f[1]));
1698 281682 : for(q = 1;;q *= 2)
1699 : {
1700 486209 : r = F2x_gcd(f, F2x_deriv(f));
1701 486209 : if (F2x_degree(r) == 0)
1702 : {
1703 232676 : gel(u, q) = f;
1704 232676 : break;
1705 : }
1706 253534 : t = F2x_div(f, r);
1707 253531 : if (F2x_degree(t) > 0)
1708 : {
1709 : long j;
1710 95499 : for(j = 1;;j++)
1711 : {
1712 209289 : v = F2x_gcd(r, t);
1713 209288 : tv = F2x_div(t, v);
1714 209288 : if (F2x_degree(tv) > 0)
1715 97353 : gel(u, j*q) = tv;
1716 209288 : if (F2x_degree(v) <= 0) break;
1717 113789 : r = F2x_div(r, v);
1718 113790 : t = v;
1719 : }
1720 95499 : if (F2x_degree(r) == 0) break;
1721 : }
1722 204529 : f = F2x_sqrt(r);
1723 : }
1724 1136043 : for (i = n; i; i--)
1725 1133429 : if (F2x_degree(gel(u,i))) break;
1726 281668 : setlg(u,i+1); return u;
1727 : }
1728 :
1729 : static GEN
1730 291326 : F2x_ddf_simple(GEN T, GEN XP)
1731 : {
1732 291326 : pari_sp av = avma, av2;
1733 : GEN f, z, Tr, X;
1734 291326 : long j, n = F2x_degree(T), v = T[1], B = n/2;
1735 291326 : if (n == 0) return cgetg(1, t_VEC);
1736 291326 : if (n == 1) return mkvec(T);
1737 135248 : z = XP; Tr = T; X = polx_F2x(v);
1738 135248 : f = const_vec(n, pol1_F2x(v));
1739 135249 : av2 = avma;
1740 355299 : for (j = 1; j <= B; j++)
1741 : {
1742 232370 : GEN u = F2x_gcd(Tr, F2x_add(z, X));
1743 232387 : if (F2x_degree(u))
1744 : {
1745 49140 : gel(f, j) = u;
1746 49140 : Tr = F2x_div(Tr, u);
1747 49135 : av2 = avma;
1748 183261 : } else z = gerepileuptoleaf(av2, z);
1749 232435 : if (!F2x_degree(Tr)) break;
1750 220066 : z = F2xq_sqr(z, Tr);
1751 : }
1752 135240 : if (F2x_degree(Tr)) gel(f, F2x_degree(Tr)) = Tr;
1753 135244 : return gerepilecopy(av, f);
1754 : }
1755 :
1756 : GEN
1757 7 : F2x_ddf(GEN T)
1758 : {
1759 : GEN XP;
1760 7 : T = F2x_get_red(T);
1761 7 : XP = F2x_Frobenius(T);
1762 7 : return F2x_ddf_to_ddf2(F2x_ddf_simple(T, XP));
1763 : }
1764 :
1765 : static GEN
1766 22165 : F2xq_frobtrace(GEN a, long d, GEN T)
1767 : {
1768 22165 : pari_sp av = avma;
1769 : long i;
1770 22165 : GEN x = a;
1771 63818 : for(i=1; i<d; i++)
1772 : {
1773 41653 : x = F2x_add(a, F2xq_sqr(x,T));
1774 41653 : if (gc_needed(av, 2))
1775 0 : x = gerepileuptoleaf(av, x);
1776 : }
1777 22165 : return x;
1778 : }
1779 :
1780 : static void
1781 33112 : F2x_edf_simple(GEN Tp, GEN XP, long d, GEN V, long idx)
1782 : {
1783 33112 : long n = F2x_degree(Tp), r = n/d;
1784 : GEN T, f, ff;
1785 33112 : if (r==1) { gel(V, idx) = Tp; return; }
1786 11133 : T = Tp;
1787 11133 : XP = F2x_rem(XP, T);
1788 : while (1)
1789 11032 : {
1790 22165 : pari_sp btop = avma;
1791 : long df;
1792 22165 : GEN g = random_F2x(n, Tp[1]);
1793 22165 : GEN t = F2xq_frobtrace(g, d, T);
1794 22165 : if (lgpol(t) == 0) continue;
1795 16570 : f = F2x_gcd(t, Tp); df = F2x_degree(f);
1796 16570 : if (df > 0 && df < n) break;
1797 5437 : set_avma(btop);
1798 : }
1799 11133 : ff = F2x_div(Tp, f);
1800 11133 : F2x_edf_simple(f, XP, d, V, idx);
1801 11133 : F2x_edf_simple(ff, XP, d, V, idx+F2x_degree(f)/d);
1802 : }
1803 :
1804 : static GEN
1805 291317 : F2x_factor_Shoup(GEN T)
1806 : {
1807 291317 : long i, n, s = 0;
1808 : GEN XP, D, V;
1809 : pari_timer ti;
1810 291317 : n = F2x_degree(T);
1811 291319 : if (DEBUGLEVEL>=6) timer_start(&ti);
1812 291319 : XP = F2x_Frobenius(T);
1813 291319 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1814 291319 : D = F2x_ddf_simple(T, XP);
1815 291321 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1816 977056 : for (i = 1; i <= n; i++)
1817 685751 : s += F2x_degree(gel(D,i))/i;
1818 291305 : V = cgetg(s+1, t_COL);
1819 977107 : for (i = 1, s = 1; i <= n; i++)
1820 : {
1821 685742 : GEN Di = gel(D,i);
1822 685742 : long ni = F2x_degree(Di), ri = ni/i;
1823 685729 : if (ni == 0) continue;
1824 328077 : if (ni == i) { gel(V, s++) = Di; continue; }
1825 10789 : F2x_edf_simple(Di, XP, i, V, s);
1826 10846 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_edf(%ld)",i);
1827 10846 : s += ri;
1828 : }
1829 291365 : return V;
1830 : }
1831 :
1832 : static GEN
1833 246387 : F2x_factor_Cantor(GEN T)
1834 : {
1835 246387 : GEN E, F, V = F2x_factor_squarefree(T);
1836 246390 : long i, j, l = lg(V);
1837 246390 : E = cgetg(l, t_VEC);
1838 246388 : F = cgetg(l, t_VEC);
1839 881271 : for (i=1, j=1; i < l; i++)
1840 634880 : if (F2x_degree(gel(V,i)))
1841 : {
1842 291318 : GEN Fj = F2x_factor_Shoup(gel(V,i));
1843 291321 : gel(F, j) = Fj;
1844 291321 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1845 291320 : j++;
1846 : }
1847 246391 : return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
1848 : }
1849 :
1850 : #if 0
1851 : static GEN
1852 : F2x_simplefact_Shoup(GEN T)
1853 : {
1854 : long i, n, s = 0, j = 1, k;
1855 : GEN XP, D, V;
1856 : pari_timer ti;
1857 : n = F2x_degree(T);
1858 : if (DEBUGLEVEL>=6) timer_start(&ti);
1859 : XP = F2x_Frobenius(T);
1860 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1861 : D = F2x_ddf_simple(T, XP);
1862 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1863 : for (i = 1; i <= n; i++)
1864 : s += F2x_degree(gel(D,i))/i;
1865 : V = cgetg(s+1, t_VECSMALL);
1866 : for (i = 1; i <= n; i++)
1867 : {
1868 : long ni = F2x_degree(gel(D,i)), ri = ni/i;
1869 : if (ni == 0) continue;
1870 : for (k = 1; k <= ri; k++)
1871 : V[j++] = i;
1872 : }
1873 : return V;
1874 : }
1875 : static GEN
1876 : F2x_simplefact_Cantor(GEN T)
1877 : {
1878 : GEN E, F, V = F2x_factor_squarefree(T);
1879 : long i, j, l = lg(V);
1880 : F = cgetg(l, t_VEC);
1881 : E = cgetg(l, t_VEC);
1882 : for (i=1, j=1; i < l; i++)
1883 : if (F2x_degree(gel(V,i)))
1884 : {
1885 : GEN Fj = F2x_simplefact_Shoup(gel(V,i));
1886 : gel(F, j) = Fj;
1887 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1888 : j++;
1889 : }
1890 : return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
1891 : }
1892 : static int
1893 : F2x_isirred_Cantor(GEN T)
1894 : {
1895 : pari_sp av = avma;
1896 : pari_timer ti;
1897 : long n;
1898 : GEN dT = F2x_deriv(T);
1899 : GEN XP, D;
1900 : if (F2x_degree(F2x_gcd(T, dT)) != 0) return gc_bool(av,0);
1901 : n = F2x_degree(T);
1902 : if (DEBUGLEVEL>=6) timer_start(&ti);
1903 : XP = F2x_Frobenius(T);
1904 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1905 : D = F2x_ddf_simple(T, XP);
1906 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1907 : return gc_bool(av, F2x_degree(gel(D,n)) == n);
1908 : }
1909 : #endif
1910 :
1911 : /* driver for Cantor factorization, assume deg f > 2; not competitive for
1912 : * flag != 0, or as deg f increases */
1913 : static GEN
1914 246387 : F2x_Cantor_i(GEN f, long flag)
1915 : {
1916 : switch(flag)
1917 : {
1918 246387 : default: return F2x_factor_Cantor(f);
1919 : #if 0
1920 : case 1: return F2x_simplefact_Cantor(f);
1921 : case 2: return F2x_isirred_Cantor(f)? gen_1: NULL;
1922 : #endif
1923 : }
1924 : }
1925 : static GEN
1926 1467174 : F2x_factor_i(GEN f, long flag)
1927 : {
1928 1467174 : long d = F2x_degree(f);
1929 1467205 : if (d <= 2) return F2x_factor_deg2(f,d,flag);
1930 257088 : return (flag == 0 && d <= 20)? F2x_Cantor_i(f, flag)
1931 538772 : : F2x_Berlekamp_i(f, flag);
1932 : }
1933 :
1934 : GEN
1935 0 : F2x_degfact(GEN f)
1936 : {
1937 0 : pari_sp av = avma;
1938 0 : GEN z = F2x_factor_i(f, 1);
1939 0 : return gerepilecopy(av, z);
1940 : }
1941 :
1942 : int
1943 238 : F2x_is_irred(GEN f) { return !!F2x_factor_i(f, 2); }
1944 :
1945 : /* Adapted from Shoup NTL */
1946 : GEN
1947 7040909 : Flx_factor_squarefree(GEN f, ulong p)
1948 : {
1949 7040909 : long i, q, n = degpol(f);
1950 7040330 : GEN u = const_vec(n+1, pol1_Flx(f[1]));
1951 7043777 : for(q = 1;;q *= p)
1952 86646 : {
1953 7130423 : GEN t, v, tv, r = Flx_gcd(f, Flx_deriv(f, p), p);
1954 7124275 : if (degpol(r) == 0) { gel(u, q) = f; break; }
1955 565534 : t = Flx_div(f, r, p);
1956 565531 : if (degpol(t) > 0)
1957 : {
1958 : long j;
1959 485820 : for(j = 1;;j++)
1960 : {
1961 1195589 : v = Flx_gcd(r, t, p);
1962 1195572 : tv = Flx_div(t, v, p);
1963 1195566 : if (degpol(tv) > 0)
1964 728406 : gel(u, j*q) = Flx_normalize(tv, p);
1965 1195560 : if (degpol(v) <= 0) break;
1966 709744 : r = Flx_div(r, v, p);
1967 709769 : t = v;
1968 : }
1969 485821 : if (degpol(r) == 0) break;
1970 : }
1971 86644 : f = Flx_normalize(Flx_deflate(r, p), p);
1972 : }
1973 30217568 : for (i = n; i; i--)
1974 30217862 : if (degpol(gel(u,i))) break;
1975 7038604 : setlg(u,i+1); return u;
1976 : }
1977 :
1978 : long
1979 3451 : Flx_ispower(GEN f, ulong k, ulong p, GEN *pt_r)
1980 : {
1981 3451 : pari_sp av = avma;
1982 : ulong lc;
1983 : GEN F;
1984 3451 : long i, n = degpol(f), v = f[1], l;
1985 3451 : if (n % k) return 0;
1986 3451 : lc = Fl_sqrtn(Flx_lead(f), k, p, NULL);
1987 3451 : if (lc == ULONG_MAX) { av = avma; return 0; }
1988 3451 : F = Flx_factor_squarefree(f, p); l = lg(F)-1;
1989 38325 : for (i = 1; i <= l; i++)
1990 34895 : if (i%k && degpol(gel(F,i))) return gc_long(av,0);
1991 3430 : if (pt_r)
1992 : {
1993 3430 : GEN r = Fl_to_Flx(lc, v), s = pol1_Flx(v);
1994 38304 : for(i = l; i >= 1; i--)
1995 : {
1996 34874 : if (i%k) continue;
1997 12453 : s = Flx_mul(s, gel(F,i), p);
1998 12453 : r = Flx_mul(r, s, p);
1999 : }
2000 3430 : *pt_r = gerepileuptoleaf(av, r);
2001 0 : } else set_avma(av);
2002 3430 : return 1;
2003 : }
2004 :
2005 : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
2006 : static GEN
2007 8937695 : Flx_ddf_Shoup(GEN T, GEN XP, ulong p)
2008 : {
2009 8937695 : pari_sp av = avma;
2010 : GEN b, g, h, F, f, Tr, xq;
2011 : long i, j, n, v, bo, ro;
2012 : long B, l, m;
2013 : pari_timer ti;
2014 8937695 : n = get_Flx_degree(T); v = get_Flx_var(T);
2015 8939491 : if (n == 0) return cgetg(1, t_VEC);
2016 8894917 : if (n == 1) return mkvec(get_Flx_mod(T));
2017 8598992 : B = n/2;
2018 8598992 : l = usqrt(B);
2019 8598831 : m = (B+l-1)/l;
2020 8598831 : T = Flx_get_red(T, p);
2021 8596863 : b = cgetg(l+2, t_VEC);
2022 8596567 : gel(b, 1) = polx_Flx(v);
2023 8597352 : gel(b, 2) = XP;
2024 8597352 : bo = brent_kung_optpow(n, l-1, 1);
2025 8598788 : ro = l<=1 ? 0:(bo-1)/(l-1) + ((n-1)/bo);
2026 8598788 : if (DEBUGLEVEL>=7) timer_start(&ti);
2027 8598788 : if (expu(p) <= ro)
2028 877553 : for (i = 3; i <= l+1; i++)
2029 486713 : gel(b, i) = Flxq_powu(gel(b, i-1), p, T, p);
2030 : else
2031 : {
2032 8208437 : xq = Flxq_powers(gel(b, 2), bo, T, p);
2033 8207461 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq baby");
2034 8936954 : for (i = 3; i <= l+1; i++)
2035 729380 : gel(b, i) = Flx_FlxqV_eval(gel(b, i-1), xq, T, p);
2036 : }
2037 8598414 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: baby");
2038 8598414 : xq = Flxq_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1), T, p);
2039 8596851 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq giant");
2040 8596851 : g = cgetg(m+1, t_VEC);
2041 8597676 : gel(g, 1) = gel(xq, 2);
2042 13958064 : for(i = 2; i <= m; i++)
2043 5360380 : gel(g, i) = Flx_FlxqV_eval(gel(g, i-1), xq, T, p);
2044 8597684 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: giant");
2045 8597684 : h = cgetg(m+1, t_VEC);
2046 22557049 : for (j = 1; j <= m; j++)
2047 : {
2048 13959519 : pari_sp av = avma;
2049 13959519 : GEN gj = gel(g, j);
2050 13959519 : GEN e = Flx_sub(gj, gel(b, 1), p);
2051 17903308 : for (i = 2; i <= l; i++)
2052 3953940 : e = Flxq_mul(e, Flx_sub(gj, gel(b, i), p), T, p);
2053 13949368 : gel(h, j) = gerepileupto(av, e);
2054 : }
2055 8597530 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: diff");
2056 8597530 : Tr = get_Flx_mod(T);
2057 8598622 : F = cgetg(m+1, t_VEC);
2058 22552908 : for (j = 1; j <= m; j++)
2059 : {
2060 13954287 : GEN u = Flx_gcd(Tr, gel(h, j), p);
2061 13956354 : if (degpol(u))
2062 : {
2063 6359042 : u = Flx_normalize(u, p);
2064 6358878 : Tr = Flx_div(Tr, u, p);
2065 : }
2066 13953111 : gel(F, j) = u;
2067 : }
2068 8598621 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: F");
2069 8598621 : f = const_vec(n, pol1_Flx(v));
2070 22559273 : for (j = 1; j <= m; j++)
2071 : {
2072 13960855 : GEN e = gel(F, j);
2073 14845075 : for (i=l-1; i >= 0; i--)
2074 : {
2075 14844367 : GEN u = Flx_gcd(e, Flx_sub(gel(g, j), gel(b, i+1), p), p);
2076 14841411 : if (degpol(u))
2077 : {
2078 6455262 : gel(f, l*j-i) = u;
2079 6455262 : e = Flx_div(e, u, p);
2080 : }
2081 14838716 : if (!degpol(e)) break;
2082 : }
2083 : }
2084 8598418 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: f");
2085 8598418 : if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
2086 8598798 : return gerepilecopy(av, f);
2087 : }
2088 :
2089 : static void
2090 441517 : Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
2091 : {
2092 441517 : long n = degpol(Tp), r = n/d;
2093 : GEN T, f, ff;
2094 : ulong p2;
2095 441517 : if (r==1) { gel(V, idx) = Tp; return; }
2096 190217 : p2 = p>>1;
2097 190217 : T = Flx_get_red(Tp, p);
2098 190219 : XP = Flx_rem(XP, T, p);
2099 : while (1)
2100 20691 : {
2101 210911 : pari_sp btop = avma;
2102 : long i;
2103 210911 : GEN g = random_Flx(n, Tp[1], p);
2104 210911 : GEN t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
2105 210911 : if (lgpol(t) == 0) continue;
2106 458467 : for(i=1; i<=10; i++)
2107 : {
2108 442869 : pari_sp btop2 = avma;
2109 442869 : GEN R = Flxq_powu(Flx_Fl_add(t, random_Fl(p), p), p2, T, p);
2110 442854 : f = Flx_gcd(Flx_Fl_add(R, p-1, p), Tp, p);
2111 442869 : if (degpol(f) > 0 && degpol(f) < n) break;
2112 252649 : set_avma(btop2);
2113 : }
2114 205817 : if (degpol(f) > 0 && degpol(f) < n) break;
2115 15597 : set_avma(btop);
2116 : }
2117 190219 : f = Flx_normalize(f, p);
2118 190219 : ff = Flx_div(Tp, f ,p);
2119 190219 : Flx_edf_simple(f, XP, d, p, V, idx);
2120 190220 : Flx_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
2121 : }
2122 : static void
2123 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx);
2124 :
2125 : static void
2126 1229908 : Flx_edf_rec(GEN T, GEN XP, GEN hp, GEN t, long d, ulong p, GEN V, long idx)
2127 : {
2128 : pari_sp av;
2129 1229908 : GEN Tp = get_Flx_mod(T);
2130 1229909 : long n = degpol(hp), vT = Tp[1];
2131 : GEN u1, u2, f1, f2;
2132 1229915 : ulong p2 = p>>1;
2133 : GEN R, h;
2134 1229915 : h = Flx_get_red(hp, p);
2135 1229896 : t = Flx_rem(t, T, p);
2136 1229865 : av = avma;
2137 : do
2138 : {
2139 2052614 : set_avma(av);
2140 2052652 : R = Flxq_powu(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p);
2141 2052468 : u1 = Flx_gcd(Flx_Fl_add(R, p-1, p), hp, p);
2142 2052615 : } while (degpol(u1)==0 || degpol(u1)==n);
2143 1229884 : f1 = Flx_gcd(Flx_Flxq_eval(u1, t, T, p), Tp, p);
2144 1229882 : f1 = Flx_normalize(f1, p);
2145 1229894 : u2 = Flx_div(hp, u1, p);
2146 1229889 : f2 = Flx_div(Tp, f1, p);
2147 1229915 : if (degpol(u1)==1)
2148 : {
2149 926432 : if (degpol(f1)==d)
2150 912086 : gel(V, idx) = f1;
2151 : else
2152 14333 : Flx_edf(f1, XP, d, p, V, idx);
2153 : }
2154 : else
2155 303528 : Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, V, idx);
2156 1229949 : idx += degpol(f1)/d;
2157 1229923 : if (degpol(u2)==1)
2158 : {
2159 921165 : if (degpol(f2)==d)
2160 907528 : gel(V, idx) = f2;
2161 : else
2162 13636 : Flx_edf(f2, XP, d, p, V, idx);
2163 : }
2164 : else
2165 308794 : Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, V, idx);
2166 1229969 : }
2167 :
2168 : static void
2169 617660 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
2170 : {
2171 617660 : long n = degpol(Tp), r = n/d, vT = Tp[1];
2172 : GEN T, h, t;
2173 : pari_timer ti;
2174 617660 : if (r==1) { gel(V, idx) = Tp; return; }
2175 617660 : T = Flx_get_red(Tp, p);
2176 617658 : XP = Flx_rem(XP, T, p);
2177 617659 : if (DEBUGLEVEL>=7) timer_start(&ti);
2178 : do
2179 : {
2180 635838 : GEN g = random_Flx(n, vT, p);
2181 635835 : t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
2182 635854 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
2183 635854 : h = Flxq_minpoly(t, T, p);
2184 635827 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
2185 635827 : } while (degpol(h) <= 1);
2186 617649 : Flx_edf_rec(T, XP, h, t, d, p, V, idx);
2187 : }
2188 :
2189 : static GEN
2190 1410116 : Flx_factor_Shoup(GEN T, ulong p)
2191 : {
2192 1410116 : long i, n, s = 0;
2193 : GEN XP, D, V;
2194 1410116 : long e = expu(p);
2195 : pari_timer ti;
2196 1410114 : n = get_Flx_degree(T);
2197 1410114 : T = Flx_get_red(T, p);
2198 1410105 : if (DEBUGLEVEL>=6) timer_start(&ti);
2199 1410105 : XP = Flx_Frobenius(T, p);
2200 1410095 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2201 1410095 : D = Flx_ddf_Shoup(T, XP, p);
2202 1410145 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2203 1410145 : s = ddf_to_nbfact(D);
2204 1410130 : V = cgetg(s+1, t_COL);
2205 8191211 : for (i = 1, s = 1; i <= n; i++)
2206 : {
2207 6781088 : GEN Di = gel(D,i);
2208 6781088 : long ni = degpol(Di), ri = ni/i;
2209 6781084 : if (ni == 0) continue;
2210 1863281 : Di = Flx_normalize(Di, p);
2211 1863320 : if (ni == i) { gel(V, s++) = Di; continue; }
2212 650754 : if (ri <= e*expu(e))
2213 589684 : Flx_edf(Di, XP, i, p, V, s);
2214 : else
2215 61078 : Flx_edf_simple(Di, XP, i, p, V, s);
2216 650760 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
2217 650760 : s += ri;
2218 : }
2219 1410123 : return V;
2220 : }
2221 :
2222 : static GEN
2223 1340418 : Flx_factor_Cantor(GEN T, ulong p)
2224 : {
2225 1340418 : GEN E, F, V = Flx_factor_squarefree(get_Flx_mod(T), p);
2226 1340418 : long i, j, l = lg(V);
2227 1340418 : F = cgetg(l, t_VEC);
2228 1340416 : E = cgetg(l, t_VEC);
2229 3132598 : for (i=1, j=1; i < l; i++)
2230 1792177 : if (degpol(gel(V,i)))
2231 : {
2232 1410116 : GEN Fj = Flx_factor_Shoup(gel(V,i), p);
2233 1410121 : gel(F, j) = Fj;
2234 1410121 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
2235 1410117 : j++;
2236 : }
2237 1340421 : return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
2238 : }
2239 :
2240 : GEN
2241 0 : Flx_ddf(GEN T, ulong p)
2242 : {
2243 : GEN XP;
2244 0 : T = Flx_get_red(T, p);
2245 0 : XP = Flx_Frobenius(T, p);
2246 0 : return ddf_to_ddf2(Flx_ddf_Shoup(T, XP, p));
2247 : }
2248 :
2249 : static GEN
2250 5367384 : Flx_simplefact_Cantor(GEN T, ulong p)
2251 : {
2252 : GEN V;
2253 : long i, l;
2254 5367384 : T = Flx_get_red(T, p);
2255 5365297 : V = Flx_factor_squarefree(get_Flx_mod(T), p); l = lg(V);
2256 10818934 : for (i=1; i < l; i++)
2257 5448370 : gel(V,i) = Flx_ddf_Shoup(gel(V,i), Flx_Frobenius(gel(V,i), p), p);
2258 5370564 : return vddf_to_simplefact(V, get_Flx_degree(T));
2259 : }
2260 :
2261 : static int
2262 1078 : Flx_isirred_Cantor(GEN Tp, ulong p)
2263 : {
2264 1078 : pari_sp av = avma;
2265 : pari_timer ti;
2266 : long n;
2267 1078 : GEN T = get_Flx_mod(Tp), dT = Flx_deriv(T, p), XP, D;
2268 1078 : if (degpol(Flx_gcd(T, dT, p)) != 0) return gc_bool(av,0);
2269 791 : n = get_Flx_degree(T);
2270 791 : T = Flx_get_red(Tp, p);
2271 791 : if (DEBUGLEVEL>=6) timer_start(&ti);
2272 791 : XP = Flx_Frobenius(T, p);
2273 791 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2274 791 : D = Flx_ddf_Shoup(T, XP, p);
2275 791 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2276 791 : return gc_bool(av, degpol(gel(D,n)) == n);
2277 : }
2278 :
2279 : /* f monic */
2280 : static GEN
2281 10493415 : Flx_factor_i(GEN f, ulong pp, long flag)
2282 : {
2283 : long d;
2284 10493415 : if (pp==2) { /*We need to handle 2 specially */
2285 73672 : GEN F = F2x_factor_i(Flx_to_F2x(f),flag);
2286 73671 : if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
2287 73671 : return F;
2288 : }
2289 10419743 : d = degpol(f);
2290 10420343 : if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
2291 6708352 : switch(flag)
2292 : {
2293 1340418 : default: return Flx_factor_Cantor(f, pp);
2294 5366856 : case 1: return Flx_simplefact_Cantor(f, pp);
2295 1078 : case 2: return Flx_isirred_Cantor(f, pp)? gen_1: NULL;
2296 : }
2297 : }
2298 :
2299 : GEN
2300 7838588 : Flx_degfact(GEN f, ulong p)
2301 : {
2302 7838588 : pari_sp av = avma;
2303 7838588 : GEN z = Flx_factor_i(Flx_normalize(f,p),p,1);
2304 7839084 : return gerepilecopy(av, z);
2305 : }
2306 :
2307 : /* T must be squarefree mod p*/
2308 : GEN
2309 1444644 : Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
2310 : {
2311 : GEN XP, D;
2312 : pari_timer ti;
2313 1444644 : long i, s, n = get_Flx_degree(T);
2314 1444592 : GEN V = const_vecsmall(n, 0);
2315 1444577 : pari_sp av = avma;
2316 1444577 : T = Flx_get_red(T, p);
2317 1444545 : if (DEBUGLEVEL>=6) timer_start(&ti);
2318 1444545 : XP = Flx_Frobenius(T, p);
2319 1444569 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2320 1444569 : D = Flx_ddf_Shoup(T, XP, p);
2321 1444922 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2322 7424834 : for (i = 1, s = 0; i <= n; i++)
2323 : {
2324 5979967 : V[i] = degpol(gel(D,i))/i;
2325 5979913 : s += V[i];
2326 : }
2327 1444867 : *nb = s;
2328 1444867 : set_avma(av); return V;
2329 : }
2330 :
2331 : long
2332 635405 : Flx_nbfact_Frobenius(GEN T, GEN XP, ulong p)
2333 : {
2334 635405 : pari_sp av = avma;
2335 635405 : long s = ddf_to_nbfact(Flx_ddf_Shoup(T, XP, p));
2336 635412 : return gc_long(av,s);
2337 : }
2338 :
2339 : /* T must be squarefree mod p*/
2340 : long
2341 635402 : Flx_nbfact(GEN T, ulong p)
2342 : {
2343 635402 : pari_sp av = avma;
2344 635402 : GEN XP = Flx_Frobenius(T, p);
2345 635404 : long n = Flx_nbfact_Frobenius(T, XP, p);
2346 635412 : return gc_long(av,n);
2347 : }
2348 :
2349 : int
2350 1057 : Flx_is_irred(GEN f, ulong p)
2351 : {
2352 1057 : pari_sp av = avma;
2353 1057 : f = Flx_normalize(f,p);
2354 1057 : return gc_bool(av, !!Flx_factor_i(f,p,2));
2355 : }
2356 :
2357 : /* Use this function when you think f is reducible, and that there are lots of
2358 : * factors. If you believe f has few factors, use FpX_nbfact(f,p)==1 instead */
2359 : int
2360 112 : FpX_is_irred(GEN f, GEN p)
2361 : {
2362 112 : pari_sp av = avma;
2363 : int z;
2364 112 : switch(ZX_factmod_init(&f,p))
2365 : {
2366 28 : case 0: z = !!F2x_factor_i(f,2); break;
2367 77 : case 1: z = !!Flx_factor_i(f,p[2],2); break;
2368 7 : default: z = !!FpX_factor_i(f,p,2); break;
2369 : }
2370 112 : return gc_bool(av,z);
2371 : }
2372 : GEN
2373 63287 : FpX_degfact(GEN f, GEN p) {
2374 63287 : pari_sp av = avma;
2375 : GEN F;
2376 63287 : switch(ZX_factmod_init(&f,p))
2377 : {
2378 567 : case 0: F = F2x_factor_i(f,1); break;
2379 62692 : case 1: F = Flx_factor_i(f,p[2],1); break;
2380 28 : default: F = FpX_factor_i(f,p,1); break;
2381 : }
2382 63287 : return gerepilecopy(av, F);
2383 : }
2384 :
2385 : #if 0
2386 : /* set x <-- x + c*y mod p */
2387 : /* x is not required to be normalized.*/
2388 : static void
2389 : Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
2390 : {
2391 : long i, lx, ly;
2392 : ulong *x=(ulong *)gx;
2393 : ulong *y=(ulong *)gy;
2394 : if (!c) return;
2395 : lx = lg(gx);
2396 : ly = lg(gy);
2397 : if (lx<ly) pari_err_BUG("lx<ly in Flx_addmul_inplace");
2398 : if (SMALL_ULONG(p))
2399 : for (i=2; i<ly; i++) x[i] = (x[i] + c*y[i]) % p;
2400 : else
2401 : for (i=2; i<ly; i++) x[i] = Fl_add(x[i], Fl_mul(c,y[i],p),p);
2402 : }
2403 : #endif
2404 :
2405 : GEN
2406 3395988 : FpX_factor(GEN f, GEN p)
2407 : {
2408 3395988 : pari_sp av = avma;
2409 : GEN F;
2410 3395988 : switch(ZX_factmod_init(&f, p))
2411 : {
2412 1377389 : case 0: F = F2x_factor_i(f,0);
2413 1377428 : F2xV_to_ZXV_inplace(gel(F,1)); break;
2414 2017015 : case 1: F = Flx_factor_i(f,p[2],0);
2415 2017021 : FlxV_to_ZXV_inplace(gel(F,1)); break;
2416 1616 : default: F = FpX_factor_i(f,p,0); break;
2417 : }
2418 3396135 : return gerepilecopy(av, F);
2419 : }
2420 :
2421 : GEN
2422 574920 : Flx_factor(GEN f, ulong p)
2423 : {
2424 574920 : pari_sp av = avma;
2425 574920 : return gerepilecopy(av, Flx_factor_i(Flx_normalize(f,p),p,0));
2426 : }
2427 : GEN
2428 15285 : F2x_factor(GEN f)
2429 : {
2430 15285 : pari_sp av = avma;
2431 15285 : return gerepilecopy(av, F2x_factor_i(f,0));
2432 : }
|
