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 4143513 : ZX_factmod_init(GEN *F, GEN p)
35 : {
36 4143513 : if (lgefint(p) == 3)
37 : {
38 4140997 : ulong pp = p[2];
39 4140997 : if (pp == 2) { *F = ZX_to_F2x(*F); return 0; }
40 2772456 : *F = ZX_to_Flx(*F, pp);
41 2772518 : if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
42 2772512 : return 1;
43 : }
44 2516 : *F = FpX_red(*F, p);
45 2518 : if (lg(*F) > 3) *F = FpX_normalize(*F, p);
46 2518 : return 2;
47 : }
48 : static GEN
49 427285 : ZX_rootmod_init(GEN F, GEN p)
50 427285 : { return lgefint(p) == 3? ZX_to_Flx(F, p[2]): FpX_red(F, p); }
51 :
52 : /* return 1,...,p-1 [not_0 = 1] or 0,...,p [not_0 = 0] */
53 : static GEN
54 600 : all_roots_mod_p(ulong p, int not_0)
55 : {
56 : GEN r;
57 : ulong i;
58 600 : if (not_0) {
59 412 : r = cgetg(p, t_VECSMALL);
60 1276 : for (i = 1; i < p; i++) r[i] = i;
61 : } else {
62 188 : r = cgetg(p+1, t_VECSMALL);
63 780 : for (i = 0; i < p; i++) r[i+1] = i;
64 : }
65 600 : return r;
66 : }
67 :
68 : /* X^n - 1 */
69 : static GEN
70 2139 : Flx_Xnm1(long sv, long n, ulong p)
71 : {
72 2139 : GEN t = cgetg(n+3, t_VECSMALL);
73 : long i;
74 2139 : t[1] = sv;
75 2139 : t[2] = p - 1;
76 6786 : for (i = 3; i <= n+1; i++) t[i] = 0;
77 2139 : t[i] = 1; return t;
78 : }
79 : /* X^n + 1 */
80 : static GEN
81 2027 : Flx_Xn1(long sv, long n, ulong p)
82 : {
83 2027 : GEN t = cgetg(n+3, t_VECSMALL);
84 : long i;
85 : (void) p;
86 2027 : t[1] = sv;
87 2027 : t[2] = 1;
88 6582 : for (i = 3; i <= n+1; i++) t[i] = 0;
89 2027 : t[i] = 1; return t;
90 : }
91 :
92 : /* assume lg(f) > 3 */
93 : static GEN
94 91366 : Flx_root_mod_2(GEN f)
95 : {
96 91366 : long i, n = lg(f)-1, c = f[2];
97 91366 : int z0 = !c;
98 91366 : c ^= 1; /* c = f[2] + f[n] mod 2, we know f[n] is odd */
99 164650 : for (i=3; i < n; i++) c ^= f[i];
100 : /* c = 0 iff f(1) = 0 (mod 2) */
101 91366 : if (z0) return c? mkvecsmall(0): mkvecsmall2(0, 1);
102 13539 : return c? cgetg(1, t_VECSMALL): mkvecsmall(1);
103 : }
104 : /* assume lg(f) > 3 */
105 : static ulong
106 91 : Flx_oneroot_mod_2(GEN f)
107 : {
108 91 : long i, n, c = f[2];
109 91 : if (!c) return 0;
110 91 : n = lg(f)-1; c = 0; /* = f[2] + f[n] (mod 2); both are odd */
111 182 : for (i=3; i < n; i++) c ^= f[i];
112 91 : return c? 2: 1;
113 : }
114 :
115 : static GEN FpX_roots_i(GEN f, GEN p);
116 :
117 : static int
118 15976731 : cmpGuGu(GEN a, GEN b) { return (ulong)a < (ulong)b? -1: (a == b? 0: 1); }
119 :
120 : /* assume that f is a ZX and p a prime */
121 : GEN
122 414167 : FpX_roots(GEN f, GEN p)
123 : {
124 414167 : pari_sp av = avma;
125 414167 : GEN y; f = ZX_rootmod_init(f, p);
126 414167 : switch(lg(f))
127 : {
128 14 : case 2: pari_err_ROOTS0("FpX_roots");
129 46388 : case 3: return cgetg(1,t_COL);
130 : }
131 367765 : if (typ(f) == t_VECSMALL)
132 : {
133 359525 : ulong pp = p[2];
134 359525 : if (pp == 2)
135 91360 : y = Flx_root_mod_2(f);
136 : else
137 : {
138 268165 : if (!odd(pp)) pari_err_PRIME("FpX_roots", p);
139 268167 : y = Flx_roots_pre(f, pp, SMALL_ULONG(pp)? 0: get_Fl_red(pp));
140 : }
141 359520 : y = Flc_to_ZC(y);
142 : }
143 : else
144 8240 : y = FpX_roots_i(f, p);
145 367752 : return gerepileupto(av, y);
146 : }
147 :
148 : /* assume x reduced mod p > 2, monic. */
149 : static int
150 21 : FpX_quad_factortype(GEN x, GEN p)
151 : {
152 21 : GEN b = gel(x,3), c = gel(x,2);
153 21 : GEN D = subii(sqri(b), shifti(c,2));
154 21 : return kronecker(D,p);
155 : }
156 : /* assume x reduced mod p, monic. Return one root, or NULL if irreducible */
157 : static GEN
158 14463 : FpX_quad_root(GEN x, GEN p, int unknown)
159 : {
160 14463 : GEN s, D, b = gel(x,3), c = gel(x,2);
161 :
162 14463 : if (absequaliu(p, 2)) {
163 0 : if (!signe(b)) return c;
164 0 : return signe(c)? NULL: gen_1;
165 : }
166 14463 : D = subii(sqri(b), shifti(c,2));
167 14463 : D = remii(D,p);
168 14463 : if (unknown && kronecker(D,p) == -1) return NULL;
169 :
170 13873 : s = Fp_sqrt(D,p);
171 : /* p is not prime, go on and give e.g. maxord a chance to recover */
172 13873 : if (!s) return NULL;
173 13865 : return Fp_halve(Fp_sub(s,b, p), p);
174 : }
175 : static GEN
176 10000 : FpX_otherroot(GEN x, GEN r, GEN p)
177 10000 : { return Fp_neg(Fp_add(gel(x,3), r, p), p); }
178 :
179 : /* disc(x^2+bx+c) = b^2 - 4c */
180 : static ulong
181 27361858 : Fl_disc_bc(ulong b, ulong c, ulong p)
182 27361858 : { return Fl_sub(Fl_sqr(b,p), Fl_double(Fl_double(c,p),p), p); }
183 : /* p > 2; allow pi = 0 */
184 : static ulong
185 24880426 : Flx_quad_root(GEN x, ulong p, ulong pi, int unknown)
186 : {
187 24880426 : ulong s, b = x[3], c = x[2];
188 24880426 : ulong D = Fl_disc_bc(b, c, p);
189 24863325 : if (unknown && krouu(D,p) == -1) return p;
190 16439854 : s = Fl_sqrt_pre(D, p, pi);
191 16484917 : if (s==~0UL) return p;
192 16484904 : return Fl_halve(Fl_sub(s,b, p), p);
193 : }
194 : static ulong
195 14688093 : Flx_otherroot(GEN x, ulong r, ulong p)
196 14688093 : { return Fl_neg(Fl_add(x[3], r, p), p); }
197 :
198 : /* 'todo' contains the list of factors to be split.
199 : * 'done' the list of finished factors, no longer touched */
200 : struct split_t { GEN todo, done; };
201 : static void
202 5069413 : split_init(struct split_t *S, long max)
203 : {
204 5069413 : S->todo = vectrunc_init(max);
205 5069179 : S->done = vectrunc_init(max);
206 5068916 : }
207 : #if 0
208 : /* move todo[i] to done */
209 : static void
210 : split_convert(struct split_t *S, long i)
211 : {
212 : long n = lg(S->todo)-1;
213 : vectrunc_append(S->done, gel(S->todo,i));
214 : if (n) gel(S->todo,i) = gel(S->todo, n);
215 : setlg(S->todo, n);
216 : }
217 : #endif
218 : /* append t to todo */
219 : static void
220 5431221 : split_add(struct split_t *S, GEN t) { vectrunc_append(S->todo, t); }
221 : /* delete todo[i], add t to done */
222 : static void
223 5432685 : split_moveto_done(struct split_t *S, long i, GEN t)
224 : {
225 5432685 : long n = lg(S->todo)-1;
226 5432685 : vectrunc_append(S->done, t);
227 5432915 : if (n) gel(S->todo,i) = gel(S->todo, n);
228 5432915 : setlg(S->todo, n);
229 :
230 5432818 : }
231 : /* append t to done */
232 : static void
233 500030 : split_add_done(struct split_t *S, GEN t)
234 500030 : { vectrunc_append(S->done, t); }
235 : /* split todo[i] into a and b */
236 : static void
237 411955 : split_todo(struct split_t *S, long i, GEN a, GEN b)
238 : {
239 411955 : gel(S->todo, i) = a;
240 411955 : split_add(S, b);
241 411953 : }
242 : /* split todo[i] into a and b, moved to done */
243 : static void
244 460728 : split_done(struct split_t *S, long i, GEN a, GEN b)
245 : {
246 460728 : split_moveto_done(S, i, a);
247 460731 : split_add_done(S, b);
248 460728 : }
249 :
250 : /* by splitting, assume p > 2 prime, deg(f) > 0 */
251 : static GEN
252 8240 : FpX_roots_i(GEN f, GEN p)
253 : {
254 : GEN pol, pol0, a, q;
255 : struct split_t S;
256 :
257 8240 : f = FpX_normalize(f, p);
258 8240 : split_init(&S, lg(f)-1);
259 8240 : settyp(S.done, t_COL);
260 8240 : if (ZX_valrem(f, &f)) split_add_done(&S, gen_0);
261 8240 : switch(degpol(f))
262 : {
263 7 : case 0: return ZC_copy(S.done);
264 14 : case 1: split_add_done(&S, subii(p, gel(f,2))); return ZC_copy(S.done);
265 3524 : case 2: {
266 3524 : GEN s, r = FpX_quad_root(f, p, 1);
267 3524 : if (r) {
268 3524 : split_add_done(&S, r);
269 3524 : s = FpX_otherroot(f,r, p);
270 : /* f not known to be square free yet */
271 3524 : if (!equalii(r, s)) split_add_done(&S, s);
272 : }
273 3524 : return sort(S.done);
274 : }
275 : }
276 :
277 4695 : a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
278 4695 : if (lg(a) < 3) pari_err_PRIME("rootmod",p);
279 4695 : a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
280 4695 : a = FpX_gcd(f,a, p);
281 4695 : if (!degpol(a)) return ZC_copy(S.done);
282 4315 : split_add(&S, FpX_normalize(a,p));
283 :
284 4315 : q = shifti(p,-1);
285 4315 : pol0 = icopy(gen_1); /* constant term, will vary in place */
286 4315 : pol = deg1pol_shallow(gen_1, pol0, varn(f));
287 4315 : for (pol0[2] = 1;; pol0[2]++)
288 10749 : {
289 15064 : long j, l = lg(S.todo);
290 15064 : if (l == 1) return sort(S.done);
291 10756 : if (pol0[2] == 100 && !BPSW_psp(p)) pari_err_PRIME("polrootsmod",p);
292 28738 : for (j = 1; j < l; j++)
293 : {
294 17989 : GEN b, r, s, c = gel(S.todo,j);
295 17989 : switch(degpol(c))
296 : { /* convert linear and quadratics to roots, try to split the rest */
297 4330 : case 1:
298 4330 : split_moveto_done(&S, j, subii(p, gel(c,2)));
299 4330 : j--; l--; break;
300 6219 : case 2:
301 6219 : r = FpX_quad_root(c, p, 0);
302 6219 : if (!r) pari_err_PRIME("polrootsmod",p);
303 6212 : s = FpX_otherroot(c,r, p);
304 6212 : split_done(&S, j, r, s);
305 6212 : j--; l--; break;
306 7440 : default:
307 7440 : b = FpXQ_pow(pol,q, c,p);
308 7440 : if (degpol(b) <= 0) continue;
309 6234 : b = FpX_gcd(c,FpX_Fp_sub_shallow(b,gen_1,p), p);
310 6234 : if (!degpol(b)) continue;
311 6234 : b = FpX_normalize(b, p);
312 6234 : c = FpX_div(c,b, p);
313 6234 : split_todo(&S, j, b, c);
314 : }
315 : }
316 : }
317 : }
318 :
319 : /* Assume f is normalized; allow pi = 0 */
320 : static ulong
321 200102 : Flx_cubic_root(GEN ff, ulong p, ulong pi)
322 : {
323 200102 : GEN f = Flx_normalize(ff,p);
324 200101 : ulong a = f[4], b=f[3], c=f[2], p3 = p%3==1 ? (2*p+1)/3 :(p+1)/3;
325 : ulong t, t2, A, B2, B, A3, A33, S, P, D;
326 200101 : if (pi)
327 : {
328 200103 : t = Fl_mul_pre(a, p3, p, pi);
329 200108 : t2 = Fl_sqr_pre(t, p, pi);
330 200104 : A = Fl_sub(b, Fl_triple(t2, p), p);
331 200101 : B = Fl_sub(c, Fl_mul_pre(t, Fl_add(A, t2, p), p, pi), p);
332 200108 : A3 = Fl_mul_pre(A, p3, p, pi);
333 200110 : B2 = Fl_sqr_pre(B, p, pi);
334 : }
335 : else
336 : {
337 0 : t = Fl_mul(a, p3, p);
338 0 : t2 = Fl_sqr(t, p);
339 0 : A = Fl_sub(b, Fl_triple(t2, p), p);
340 0 : B = Fl_sub(c, Fl_mul(t, Fl_add(A, t2, p), p), p);
341 0 : A3 = Fl_mul(A, p3, p);
342 0 : B2 = Fl_sqr(B, p);
343 : }
344 200107 : A33 = Fl_powu_pre(A3, 3, p, pi);
345 200100 : D = Fl_add(B2, Fl_double(Fl_double(A33, p), p), p);
346 200100 : S = Fl_neg(B,p);
347 200104 : P = Fl_neg(A3,p);
348 200105 : if (krouu(D,p) >= 0)
349 : {
350 125932 : ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
351 125931 : ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
352 125933 : if (p%3==2) /* 1 solutions */
353 23504 : vS1 = Fl_powu_pre(S1, p - p3, p, pi);
354 : else
355 : {
356 102429 : vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
357 102431 : if (vS1==~0UL) return p; /*0 solutions*/
358 : /*3 solutions*/
359 : }
360 93523 : if (!P) return Fl_sub(vS1, t, p);
361 93160 : vS2 = pi? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): Fl_div(P, vS1, p);
362 93160 : return Fl_sub(Fl_add(vS1,vS2, p), t, p);
363 : }
364 : else
365 : {
366 74174 : pari_sp av = avma;
367 74174 : GEN S1 = mkvecsmall2(Fl_halve(S, p), (p + 1UL) >> 1);
368 74174 : GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
369 : ulong Sa;
370 74173 : if (!vS1) return p; /*0 solutions, p%3==2*/
371 74173 : Sa = vS1[1];
372 74173 : if (p%3==1) /*1 solutions*/
373 : {
374 28705 : ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
375 28705 : if (Fa!=P) Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
376 : }
377 74173 : set_avma(av);
378 74173 : return Fl_sub(Fl_double(Sa,p),t,p);
379 : }
380 : }
381 :
382 : /* Assume f is normalized */
383 : static GEN
384 125 : FpX_cubic_root(GEN ff, GEN p)
385 : {
386 125 : GEN f = FpX_normalize(ff,p);
387 125 : GEN a = gel(f,4), b = gel(f,3), c = gel(f,2);
388 125 : ulong pm3 = umodiu(p,3);
389 28 : GEN p3 = pm3==1 ? diviuexact(addiu(shifti(p,1),1),3)
390 125 : : diviuexact(addiu(p,1),3);
391 125 : GEN t = Fp_mul(a, p3, p), t2 = Fp_sqr(t, p);
392 125 : GEN A = Fp_sub(b, Fp_mulu(t2, 3, p), p);
393 125 : GEN B = Fp_addmul(c, t, Fp_sub(shifti(t2, 1), b, p), p);
394 125 : GEN A3 = Fp_mul(A, p3, p), A33 = Fp_powu(A3, 3, p);
395 125 : GEN S = Fp_neg(B,p), P = Fp_neg(A3,p);
396 125 : GEN D = Fp_add(Fp_sqr(S, p), shifti(A33, 2), p);
397 125 : if (kronecker(D,p) >= 0)
398 : {
399 28 : GEN s = Fp_sqrt(D, p), vS1, vS2;
400 28 : GEN S1 = S==s ? S: Fp_halve(Fp_sub(S, s, p), p);
401 28 : if (pm3 == 2) /* 1 solutions */
402 0 : vS1 = Fp_pow(S1, diviuexact(addiu(shifti(p, 1), 1), 3), p);
403 : else
404 : {
405 28 : vS1 = Fp_sqrtn(S1, utoi(3), p, NULL);
406 28 : if (!vS1) return p; /*0 solutions*/
407 : /*3 solutions*/
408 : }
409 28 : vS2 = P? Fp_mul(P, Fp_inv(vS1, p), p): 0;
410 28 : return Fp_sub(Fp_add(vS1,vS2, p), t, p);
411 : }
412 : else
413 : {
414 97 : pari_sp av = avma;
415 97 : GEN T = deg2pol_shallow(gen_1, gen_0, negi(D), 0);
416 97 : GEN S1 = deg1pol_shallow(Fp_halve(gen_1, p), Fp_halve(S, p), 0);
417 97 : GEN vS1 = FpXQ_sqrtn(S1, utoi(3), T, p, NULL);
418 : GEN Sa;
419 97 : if (!vS1) return p; /*0 solutions, p%3==2*/
420 97 : Sa = gel(vS1,2);
421 97 : if (pm3 == 1) /*1 solutions*/
422 : {
423 0 : GEN Fa = FpXQ_norm(vS1, T, p);
424 0 : if (!equalii(Fa,P))
425 0 : Sa = Fp_mul(Sa, Fp_div(Fa, P, p),p);
426 : }
427 97 : set_avma(av);
428 97 : return Fp_sub(shifti(Sa,1),t,p);
429 : }
430 : }
431 :
432 : /* assume p > 2 prime; if fl is set, assume that f splits mod p */
433 : static ulong
434 3925834 : Flx_oneroot_pre_i(GEN f, ulong p, ulong pi, long fl)
435 : {
436 : GEN pol, a;
437 : ulong q, PI;
438 : long da;
439 :
440 3925834 : if (Flx_val(f)) return 0;
441 3924555 : da = degpol(f); f = Flx_normalize(f, p);
442 3925297 : if (da == 1) return Fl_neg(f[2], p);
443 3911705 : PI = pi? pi: get_Fl_red(p); /* PI for Fp, pi for Fp[x] */
444 3912480 : switch(da)
445 : {
446 3393098 : case 2: return Flx_quad_root(f, p, PI, 1);
447 179699 : case 3: if (p>3) return Flx_cubic_root(f, p, PI); /*FALL THROUGH*/
448 : }
449 339683 : if (SMALL_ULONG(p)) pi = 0; /* bilinear ops faster without Fl_*_pre */
450 339683 : if (!fl)
451 : {
452 303487 : a = Flxq_powu_pre(polx_Flx(f[1]), p - 1, f,p,pi);
453 303451 : if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
454 303451 : a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
455 303422 : a = Flx_gcd_pre(f,a, p, pi);
456 36196 : } else a = f;
457 339652 : da = degpol(a);
458 339654 : if (!da) return p;
459 239612 : a = Flx_normalize(a,p);
460 :
461 239638 : q = p >> 1;
462 239638 : pol = polx_Flx(f[1]);
463 353780 : for(pol[2] = 1;; pol[2]++)
464 : {
465 353780 : if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
466 353789 : switch(da)
467 : {
468 148618 : case 1: return Fl_neg(a[2], p);
469 70634 : case 2: return Flx_quad_root(a, p, PI, 0);
470 20403 : case 3: if (p>3) return Flx_cubic_root(a, p, PI); /*FALL THROUGH*/
471 : default: {
472 114134 : GEN b = Flxq_powu_pre(pol,q, a,p,pi);
473 : long db;
474 114119 : if (degpol(b) <= 0) continue;
475 108491 : b = Flx_gcd_pre(a,Flx_Fl_add(b,p-1,p), p, pi);
476 108497 : db = degpol(b); if (!db) continue;
477 108497 : b = Flx_normalize(b, p);
478 108503 : if (db <= (da >> 1)) {
479 66939 : a = b;
480 66939 : da = db;
481 : } else {
482 41564 : a = Flx_div_pre(a,b, p, pi);
483 41568 : da -= db;
484 : }
485 : }
486 : }
487 : }
488 : }
489 : ulong
490 3888420 : Flx_oneroot_pre(GEN f, ulong p, ulong pi)
491 3888420 : { return Flx_oneroot_pre_i(f, p, pi, 0); }
492 : ulong
493 37659 : Flx_oneroot_split_pre(GEN f, ulong p, ulong pi)
494 37659 : { return Flx_oneroot_pre_i(f, p, pi, 1); }
495 :
496 : /* assume p > 3 prime */
497 : static GEN
498 5169 : FpX_oneroot_i(GEN f, GEN p)
499 : {
500 : GEN pol, pol0, a, q;
501 : long da;
502 :
503 5169 : if (ZX_val(f)) return gen_0;
504 4847 : f = FpX_normalize(f, p);
505 4847 : switch(degpol(f))
506 : {
507 815 : case 1: return subii(p, gel(f,2));
508 3837 : case 2: return FpX_quad_root(f, p, 1);
509 125 : 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 224 : 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 2478 : Flx_oneroot(GEN f, ulong p)
551 : {
552 2478 : pari_sp av = avma;
553 2478 : switch(lg(f))
554 : {
555 0 : case 2: return 0;
556 0 : case 3: return p;
557 : }
558 2478 : if (p == 2) return Flx_oneroot_mod_2(f);
559 2478 : return gc_ulong(av, Flx_oneroot_pre(f, p, SMALL_ULONG(p)? 0: get_Fl_red(p)));
560 : }
561 :
562 : ulong
563 14 : Flx_oneroot_split(GEN f, ulong p)
564 : {
565 14 : pari_sp av = avma;
566 14 : switch(lg(f))
567 : {
568 0 : case 2: return 0;
569 0 : case 3: return p;
570 : }
571 14 : if (p == 2) return Flx_oneroot_mod_2(f);
572 14 : return gc_ulong(av, Flx_oneroot_split_pre(f, p, 0));
573 : }
574 :
575 : /* assume that p is prime */
576 : GEN
577 13118 : FpX_oneroot(GEN f, GEN p)
578 : {
579 13118 : pari_sp av = avma;
580 13118 : f = ZX_rootmod_init(f, p);
581 13118 : switch(lg(f))
582 : {
583 0 : case 2: set_avma(av); return gen_0;
584 0 : case 3: return gc_NULL(av);
585 : }
586 13118 : if (typ(f) == t_VECSMALL)
587 : {
588 7949 : ulong r, pp = p[2];
589 7949 : if (pp == 2)
590 91 : r = Flx_oneroot_mod_2(f);
591 : else
592 7858 : r = Flx_oneroot_pre(f, pp, SMALL_ULONG(pp)? 0: get_Fl_red(pp));
593 7949 : set_avma(av);
594 7949 : return (r == pp)? NULL: utoi(r);
595 : }
596 5169 : f = FpX_oneroot_i(f, p);
597 5169 : if (!f) return gc_NULL(av);
598 5169 : return gerepileuptoint(av, f);
599 : }
600 :
601 : /* returns a root of unity in F_p that is suitable for finding a factor */
602 : /* of degree deg_factor of a polynomial of degree deg; the order is */
603 : /* returned in n */
604 : /* A good choice seems to be n close to deg/deg_factor; we choose n */
605 : /* twice as big and decrement until it divides p-1. */
606 : static GEN
607 154 : good_root_of_unity(GEN p, long deg, long deg_factor, long *pt_n)
608 : {
609 154 : pari_sp ltop = avma;
610 : GEN pm, factn, power, base, zeta;
611 : long n;
612 :
613 154 : pm = subis (p, 1ul);
614 336 : for (n = deg / 2 / deg_factor + 1; !dvdiu (pm, n); n--);
615 154 : factn = Z_factor(stoi(n));
616 154 : power = diviuexact (pm, n);
617 154 : base = gen_1;
618 : do {
619 259 : base = addis (base, 1l);
620 259 : zeta = Fp_pow (base, power, p);
621 : }
622 259 : while (!equaliu (Fp_order (zeta, factn, p), n));
623 154 : *pt_n = n;
624 154 : return gerepileuptoint (ltop, zeta);
625 : }
626 :
627 : GEN
628 1092 : FpX_oneroot_split(GEN fact, GEN p)
629 : {
630 1092 : pari_sp av = avma;
631 : long n, deg_f, i, dmin;
632 : GEN prim, expo, minfactor, xplusa, zeta, xpow;
633 1092 : fact = FpX_normalize(fact, p);
634 1092 : deg_f = degpol(fact);
635 1092 : if (deg_f <= 3) return FpX_oneroot(fact, p);
636 133 : minfactor = fact; /* factor of minimal degree found so far */
637 133 : dmin = degpol(minfactor);
638 133 : xplusa = pol_x(varn(fact));
639 287 : while (dmin > 3)
640 : {
641 : /* split minfactor by computing its gcd with (X+a)^exp-zeta, where */
642 : /* zeta varies over the roots of unity in F_p */
643 154 : fact = minfactor; deg_f = dmin;
644 154 : zeta = gen_1;
645 154 : prim = good_root_of_unity(p, deg_f, 1, &n);
646 154 : expo = diviuexact(subiu(p, 1), n);
647 : /* update X+a, avoid a=0 */
648 154 : gel (xplusa, 2) = addis (gel (xplusa, 2), 1);
649 154 : xpow = FpXQ_pow (xplusa, expo, fact, p);
650 325 : for (i = 0; i < n; i++)
651 : {
652 249 : GEN tmp = FpX_gcd(FpX_Fp_sub(xpow, zeta, p), fact, p);
653 249 : long dtmp = degpol(tmp);
654 249 : if (dtmp > 0 && dtmp < deg_f)
655 : {
656 160 : fact = FpX_div(fact, tmp, p); deg_f = degpol(fact);
657 160 : if (dtmp < dmin)
658 : {
659 154 : minfactor = FpX_normalize (tmp, p);
660 154 : dmin = dtmp;
661 154 : if (dmin == 1 || dmin <= (2 * deg_f) / n - 1)
662 : /* stop early to avoid too many gcds */
663 : break;
664 : }
665 : }
666 171 : zeta = Fp_mul (zeta, prim, p);
667 : }
668 : }
669 133 : return gerepileuptoint(av, FpX_oneroot(minfactor, p));
670 : }
671 :
672 : /*******************************************************************/
673 : /* */
674 : /* FACTORISATION MODULO p */
675 : /* */
676 : /*******************************************************************/
677 :
678 : /* F / E a vector of vectors of factors / exponents of virtual length l
679 : * (their real lg may be larger). Set their lg to j, concat and return [F,E] */
680 : static GEN
681 1697264 : FE_concat(GEN F, GEN E, long l)
682 : {
683 1697264 : setlg(E,l); E = shallowconcat1(E);
684 1697265 : setlg(F,l); F = shallowconcat1(F); return mkvec2(F,E);
685 : }
686 :
687 : static GEN
688 14 : ddf_to_ddf2_i(GEN V, long fl)
689 : {
690 : GEN F, D;
691 14 : long i, j, l = lg(V);
692 14 : F = cgetg(l, t_VEC);
693 14 : D = cgetg(l, t_VECSMALL);
694 112 : for (i = j = 1; i < l; i++)
695 : {
696 98 : GEN Vi = gel(V,i);
697 98 : if ((fl==2 && F2x_degree(Vi) == 0)
698 98 : ||(fl==0 && degpol(Vi) == 0)) continue;
699 35 : gel(F,j) = Vi;
700 35 : uel(D,j) = i; j++;
701 : }
702 14 : setlg(F,j);
703 14 : setlg(D,j); return mkvec2(F,D);
704 : }
705 :
706 : GEN
707 7 : ddf_to_ddf2(GEN V)
708 7 : { return ddf_to_ddf2_i(V, 0); }
709 :
710 : static GEN
711 7 : F2x_ddf_to_ddf2(GEN V)
712 7 : { return ddf_to_ddf2_i(V, 2); }
713 :
714 : GEN
715 5293068 : vddf_to_simplefact(GEN V, long d)
716 : {
717 : GEN E, F;
718 5293068 : long i, j, c, l = lg(V);
719 5293068 : F = cgetg(d+1, t_VECSMALL);
720 5292955 : E = cgetg(d+1, t_VECSMALL);
721 10669483 : for (i = c = 1; i < l; i++)
722 : {
723 5376193 : GEN Vi = gel(V,i);
724 5376193 : long l = lg(Vi);
725 27075731 : for (j = 1; j < l; j++)
726 : {
727 21699167 : long k, n = degpol(gel(Vi,j)) / j;
728 32866783 : for (k = 1; k <= n; k++) { uel(F,c) = j; uel(E,c) = i; c++; }
729 : }
730 : }
731 5293290 : setlg(F,c);
732 5293307 : setlg(E,c);
733 5293148 : return sort_factor(mkvec2(F,E), (void*)&cmpGuGu, cmp_nodata);
734 : }
735 :
736 : /* product of terms of degree 1 in factorization of f */
737 : GEN
738 189919 : FpX_split_part(GEN f, GEN p)
739 : {
740 189919 : long n = degpol(f);
741 189919 : GEN z, X = pol_x(varn(f));
742 189919 : if (n <= 1) return f;
743 187979 : f = FpX_red(f, p);
744 187978 : z = FpX_sub(FpX_Frobenius(f, p), X, p);
745 187980 : return FpX_gcd(z,f,p);
746 : }
747 :
748 : /* Compute the number of roots in Fp without counting multiplicity
749 : * return -1 for 0 polynomial. lc(f) must be prime to p. */
750 : long
751 99541 : FpX_nbroots(GEN f, GEN p)
752 : {
753 99541 : pari_sp av = avma;
754 99541 : GEN z = FpX_split_part(f, p);
755 99541 : return gc_long(av, degpol(z));
756 : }
757 :
758 : /* 1 < deg(f) <= p */
759 : static int
760 81346 : Flx_is_totally_split_i(GEN f, ulong p)
761 : {
762 81346 : GEN F = Flx_Frobenius(f, p);
763 81347 : return degpol(F)==1 && uel(F,2)==0UL && uel(F,3)==1UL;
764 : }
765 : int
766 81353 : Flx_is_totally_split(GEN f, ulong p)
767 : {
768 81353 : pari_sp av = avma;
769 81353 : ulong n = degpol(f);
770 81353 : if (n <= 1) return 1;
771 81346 : if (n > p) return 0; /* includes n < 0 */
772 81346 : return gc_bool(av, Flx_is_totally_split_i(f,p));
773 : }
774 : int
775 0 : FpX_is_totally_split(GEN f, GEN p)
776 : {
777 0 : pari_sp av = avma;
778 0 : ulong n = degpol(f);
779 : int u;
780 0 : if (n <= 1) return 1;
781 0 : if (abscmpui(n, p) > 0) return 0; /* includes n < 0 */
782 0 : if (lgefint(p) != 3)
783 0 : u = gequalX(FpX_Frobenius(FpX_red(f,p), p));
784 : else
785 : {
786 0 : ulong pp = (ulong)p[2];
787 0 : u = Flx_is_totally_split_i(ZX_to_Flx(f,pp), pp);
788 : }
789 0 : return gc_bool(av, u);
790 : }
791 :
792 : long
793 3011192 : Flx_nbroots(GEN f, ulong p)
794 : {
795 3011192 : long n = degpol(f);
796 : ulong pi;
797 3011191 : pari_sp av = avma;
798 : GEN z;
799 3011191 : if (n <= 1) return n;
800 2997660 : if (n == 2)
801 : {
802 : ulong D;
803 65490 : if (p==2) return (f[2]==0) + (f[2]!=f[3]);
804 64475 : D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
805 64475 : return 1 + krouu(D,p);
806 : }
807 2932170 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
808 2932170 : z = Flx_sub(Flx_Frobenius_pre(f, p, pi), polx_Flx(f[1]), p);
809 2932195 : z = Flx_gcd_pre(z, f, p, pi);
810 2932196 : return gc_long(av, degpol(z));
811 : }
812 :
813 : long
814 4256 : FpX_ddf_degree(GEN T, GEN XP, GEN p)
815 : {
816 4256 : pari_sp av = avma;
817 : GEN X, b, g, xq;
818 : long i, j, n, v, B, l, m;
819 : pari_timer ti;
820 : hashtable h;
821 :
822 4256 : n = get_FpX_degree(T); v = get_FpX_var(T);
823 4256 : X = pol_x(v);
824 4256 : if (ZX_equal(X,XP)) return 1;
825 4256 : B = n/2;
826 4256 : l = usqrt(B);
827 4256 : m = (B+l-1)/l;
828 4256 : T = FpX_get_red(T, p);
829 4256 : hash_init_GEN(&h, l+2, ZX_equal, 1);
830 4256 : hash_insert_long(&h, X, 0);
831 4256 : hash_insert_long(&h, XP, 1);
832 4256 : if (DEBUGLEVEL>=7) timer_start(&ti);
833 4256 : b = XP;
834 4256 : xq = FpXQ_powers(b, brent_kung_optpow(n, l-1, 1), T, p);
835 4256 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq baby");
836 10178 : for (i = 3; i <= l+1; i++)
837 : {
838 6601 : b = FpX_FpXQV_eval(b, xq, T, p);
839 6601 : if (gequalX(b)) return gc_long(av,i-1);
840 5922 : hash_insert_long(&h, b, i-1);
841 : }
842 3577 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: baby");
843 3577 : g = b;
844 3577 : xq = FpXQ_powers(g, brent_kung_optpow(n, m, 1), T, p);
845 3577 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq giant");
846 12208 : for(i = 2; i <= m+1; i++)
847 : {
848 10619 : g = FpX_FpXQV_eval(g, xq, T, p);
849 10619 : if (hash_haskey_long(&h, g, &j)) return gc_long(av, l*i-j);
850 : }
851 1589 : return gc_long(av,n);
852 : }
853 :
854 : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
855 : static GEN
856 1030 : FpX_ddf_Shoup(GEN T, GEN XP, GEN p)
857 : {
858 : GEN b, g, h, F, f, Tr, xq;
859 : long i, j, n, v, B, l, m;
860 : pari_timer ti;
861 :
862 1030 : n = get_FpX_degree(T); v = get_FpX_var(T);
863 1030 : if (n == 0) return cgetg(1, t_VEC);
864 1030 : if (n == 1) return mkvec(get_FpX_mod(T));
865 852 : B = n/2;
866 852 : l = usqrt(B);
867 852 : m = (B+l-1)/l;
868 852 : T = FpX_get_red(T, p);
869 852 : b = cgetg(l+2, t_VEC);
870 852 : gel(b, 1) = pol_x(v);
871 852 : gel(b, 2) = XP;
872 852 : if (DEBUGLEVEL>=7) timer_start(&ti);
873 852 : xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1), T, p);
874 852 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq baby");
875 1081 : for (i = 3; i <= l+1; i++)
876 229 : gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
877 852 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: baby");
878 852 : xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1), T, p);
879 852 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq giant");
880 852 : g = cgetg(m+1, t_VEC);
881 852 : gel(g, 1) = gel(xq, 2);
882 1843 : for(i = 2; i <= m; i++) gel(g, i) = FpX_FpXQV_eval(gel(g, i-1), xq, T, p);
883 852 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: giant");
884 852 : h = cgetg(m+1, t_VEC);
885 2695 : for (j = 1; j <= m; j++)
886 : {
887 1843 : pari_sp av = avma;
888 1843 : GEN gj = gel(g,j), e = FpX_sub(gj, gel(b,1), p);
889 3076 : for (i = 2; i <= l; i++) e = FpXQ_mul(e, FpX_sub(gj, gel(b,i), p), T, p);
890 1843 : gel(h,j) = gerepileupto(av, e);
891 : }
892 852 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: diff");
893 852 : Tr = get_FpX_mod(T);
894 852 : F = cgetg(m+1, t_VEC);
895 2695 : for (j = 1; j <= m; j++)
896 : {
897 1843 : GEN u = FpX_gcd(Tr, gel(h,j), p);
898 1843 : if (degpol(u))
899 : {
900 485 : u = FpX_normalize(u, p);
901 485 : Tr = FpX_div(Tr, u, p);
902 : }
903 1843 : gel(F,j) = u;
904 : }
905 852 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: F");
906 852 : f = const_vec(n, pol_1(v));
907 2695 : for (j = 1; j <= m; j++)
908 : {
909 1843 : GEN e = gel(F, j);
910 1915 : for (i=l-1; i >= 0; i--)
911 : {
912 1915 : GEN u = FpX_gcd(e, FpX_sub(gel(g, j), gel(b, i+1), p), p);
913 1915 : if (degpol(u))
914 : {
915 515 : u = FpX_normalize(u, p);
916 515 : gel(f, l*j-i) = u;
917 515 : e = FpX_div(e, u, p);
918 : }
919 1915 : if (!degpol(e)) break;
920 : }
921 : }
922 852 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: f");
923 852 : if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
924 852 : return f;
925 : }
926 :
927 : static void
928 0 : FpX_edf_simple(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
929 : {
930 0 : long n = degpol(Tp), r = n/d, ct = 0;
931 : GEN T, f, ff, p2;
932 0 : if (r==1) { gel(V, idx) = Tp; return; }
933 0 : p2 = shifti(p,-1);
934 0 : T = FpX_get_red(Tp, p);
935 0 : XP = FpX_rem(XP, T, p);
936 : while (1)
937 0 : {
938 0 : pari_sp btop = avma;
939 : long i;
940 0 : GEN g = random_FpX(n, varn(Tp), p);
941 0 : GEN t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
942 0 : if (signe(t) == 0) continue;
943 0 : for(i=1; i<=10; i++)
944 : {
945 0 : pari_sp btop2 = avma;
946 0 : GEN R = FpXQ_pow(FpX_Fp_add(t, randomi(p), p), p2, T, p);
947 0 : f = FpX_gcd(FpX_Fp_sub(R, gen_1, p), Tp, p);
948 0 : if (degpol(f) > 0 && degpol(f) < n) break;
949 0 : set_avma(btop2);
950 : }
951 0 : if (degpol(f) > 0 && degpol(f) < n) break;
952 0 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_simple",p);
953 0 : set_avma(btop);
954 : }
955 0 : f = FpX_normalize(f, p);
956 0 : ff = FpX_div(Tp, f ,p);
957 0 : FpX_edf_simple(f, XP, d, p, V, idx);
958 0 : FpX_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
959 : }
960 :
961 : static void
962 1075 : FpX_edf_rec(GEN T, GEN hp, GEN t, long d, GEN p2, GEN p, GEN V, long idx)
963 : {
964 : pari_sp av;
965 1075 : GEN Tp = get_FpX_mod(T);
966 1075 : long n = degpol(hp), vT = varn(Tp), ct = 0;
967 : GEN u1, u2, f1, f2, R, h;
968 1075 : h = FpX_get_red(hp, p);
969 1075 : t = FpX_rem(t, T, p);
970 1075 : av = avma;
971 : do
972 : {
973 1695 : set_avma(av);
974 1695 : R = FpXQ_pow(deg1pol(gen_1, randomi(p), vT), p2, h, p);
975 1695 : u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
976 1695 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_rec",p);
977 1695 : } while (degpol(u1)==0 || degpol(u1)==n);
978 1075 : f1 = FpX_gcd(FpX_FpXQ_eval(u1, t, T, p), Tp, p);
979 1075 : f1 = FpX_normalize(f1, p);
980 1075 : u2 = FpX_div(hp, u1, p);
981 1075 : f2 = FpX_div(Tp, f1, p);
982 1075 : if (degpol(u1)==1)
983 774 : gel(V, idx) = f1;
984 : else
985 301 : FpX_edf_rec(FpX_get_red(f1, p), u1, t, d, p2, p, V, idx);
986 1075 : idx += degpol(f1)/d;
987 1075 : if (degpol(u2)==1)
988 739 : gel(V, idx) = f2;
989 : else
990 336 : FpX_edf_rec(FpX_get_red(f2, p), u2, t, d, p2, p, V, idx);
991 1075 : }
992 :
993 : /* assume Tp a squarefree product of r > 1 irred. factors of degree d */
994 : static void
995 438 : FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
996 : {
997 438 : long n = degpol(Tp), r = n/d, vT = varn(Tp), ct = 0;
998 : GEN T, h, t;
999 : pari_timer ti;
1000 :
1001 438 : T = FpX_get_red(Tp, p);
1002 438 : XP = FpX_rem(XP, T, p);
1003 438 : if (DEBUGLEVEL>=7) timer_start(&ti);
1004 : do
1005 : {
1006 438 : GEN g = random_FpX(n, vT, p);
1007 438 : t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
1008 438 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
1009 438 : h = FpXQ_minpoly(t, T, p);
1010 438 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
1011 438 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf",p);
1012 438 : } while (degpol(h) != r);
1013 438 : FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
1014 438 : }
1015 :
1016 : static GEN
1017 1016 : FpX_factor_Shoup(GEN T, GEN p)
1018 : {
1019 1016 : long i, n, s = 0;
1020 : GEN XP, D, V;
1021 1016 : long e = expi(p);
1022 : pari_timer ti;
1023 1016 : n = get_FpX_degree(T);
1024 1016 : T = FpX_get_red(T, p);
1025 1016 : if (DEBUGLEVEL>=6) timer_start(&ti);
1026 1016 : XP = FpX_Frobenius(T, p);
1027 1016 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
1028 1016 : D = FpX_ddf_Shoup(T, XP, p);
1029 1016 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
1030 1016 : s = ddf_to_nbfact(D);
1031 1016 : V = cgetg(s+1, t_COL);
1032 7316 : for (i = 1, s = 1; i <= n; i++)
1033 : {
1034 6300 : GEN Di = gel(D,i);
1035 6300 : long ni = degpol(Di), ri = ni/i;
1036 6300 : if (ni == 0) continue;
1037 1137 : Di = FpX_normalize(Di, p);
1038 1137 : if (ni == i) { gel(V, s++) = Di; continue; }
1039 438 : if (ri <= e*expu(e))
1040 438 : FpX_edf(Di, XP, i, p, V, s);
1041 : else
1042 0 : FpX_edf_simple(Di, XP, i, p, V, s);
1043 438 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
1044 438 : s += ri;
1045 : }
1046 1016 : return V;
1047 : }
1048 :
1049 : long
1050 2173363 : ddf_to_nbfact(GEN D)
1051 : {
1052 2173363 : long l = lg(D), i, s = 0;
1053 14233518 : for(i = 1; i < l; i++) s += degpol(gel(D,i))/i;
1054 2173365 : return s;
1055 : }
1056 :
1057 : /* Yun algorithm: Assume p > degpol(T) */
1058 : static GEN
1059 1678 : FpX_factor_Yun(GEN T, GEN p)
1060 : {
1061 1678 : long n = degpol(T), i = 1;
1062 1678 : GEN a, b, c, d = FpX_deriv(T, p);
1063 1678 : GEN V = cgetg(n+1,t_VEC);
1064 1678 : a = FpX_gcd(T, d, p);
1065 1678 : if (degpol(a) == 0) return mkvec(T);
1066 589 : b = FpX_div(T, a, p);
1067 : do
1068 : {
1069 2648 : c = FpX_div(d, a, p);
1070 2648 : d = FpX_sub(c, FpX_deriv(b, p), p);
1071 2648 : a = FpX_normalize(FpX_gcd(b, d, p), p);
1072 2648 : gel(V, i++) = a;
1073 2648 : b = FpX_div(b, a, p);
1074 2648 : } while (degpol(b));
1075 589 : setlg(V, i); return V;
1076 : }
1077 : GEN
1078 335713 : FpX_factor_squarefree(GEN T, GEN p)
1079 : {
1080 335713 : if (lgefint(p)==3)
1081 : {
1082 335001 : ulong pp = (ulong)p[2];
1083 335001 : GEN u = Flx_factor_squarefree(ZX_to_Flx(T,pp), pp);
1084 335001 : return FlxV_to_ZXV(u);
1085 : }
1086 712 : return FpX_factor_Yun(T, p);
1087 : }
1088 :
1089 : GEN
1090 219548 : FpX_roots_mult(GEN T, long n, GEN p)
1091 : {
1092 219548 : GEN V = FpX_factor_squarefree(T,p), W;
1093 219548 : long l = lg(V), i;
1094 219548 : if (l<=n) return cgetg(1,t_COL);
1095 35663 : W = cgetg(l-n+1,t_VEC);
1096 117049 : for (i = n; i < l; i++)
1097 81386 : gel(W,i-n+1) = FpX_roots(gel(V,i), p);
1098 35663 : return shallowconcat1(W);
1099 : }
1100 :
1101 : long
1102 168 : FpX_ispower(GEN f, ulong k, GEN p, GEN *pt_r)
1103 : {
1104 168 : pari_sp av = avma;
1105 : GEN lc, F;
1106 168 : long i, l, n = degpol(f), v = varn(f);
1107 168 : if (n % k) return 0;
1108 168 : if (lgefint(p)==3)
1109 : {
1110 126 : ulong pp = p[2];
1111 126 : GEN fp = ZX_to_Flx(f, pp);
1112 126 : if (!Flx_ispower(fp, k, pp, pt_r)) return gc_long(av,0);
1113 105 : if (pt_r) *pt_r = gerepileupto(av, Flx_to_ZX(*pt_r)); else set_avma(av);
1114 105 : return 1;
1115 : }
1116 42 : lc = Fp_sqrtn(leading_coeff(f), stoi(k), p, NULL);
1117 42 : if (!lc) { av = avma; return 0; }
1118 42 : F = FpX_factor_Yun(f, p); l = lg(F)-1;
1119 1491 : for(i=1; i <= l; i++)
1120 1456 : if (i%k && degpol(gel(F,i))) return gc_long(av,0);
1121 35 : if (pt_r)
1122 : {
1123 35 : GEN r = scalarpol(lc, v), s = pol_1(v);
1124 1484 : for (i=l; i>=1; i--)
1125 : {
1126 1449 : if (i%k) continue;
1127 294 : s = FpX_mul(s, gel(F,i), p);
1128 294 : r = FpX_mul(r, s, p);
1129 : }
1130 35 : *pt_r = gerepileupto(av, r);
1131 0 : } else av = avma;
1132 35 : return 1;
1133 : }
1134 :
1135 : static GEN
1136 910 : FpX_factor_Cantor(GEN T, GEN p)
1137 : {
1138 910 : GEN E, F, V = FpX_factor_Yun(T, p);
1139 910 : long i, j, l = lg(V);
1140 910 : F = cgetg(l, t_VEC);
1141 910 : E = cgetg(l, t_VEC);
1142 2439 : for (i=1, j=1; i < l; i++)
1143 1529 : if (degpol(gel(V,i)))
1144 : {
1145 1016 : GEN Fj = FpX_factor_Shoup(gel(V,i), p);
1146 1016 : gel(F, j) = Fj;
1147 1016 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1148 1016 : j++;
1149 : }
1150 910 : return sort_factor_pol(FE_concat(F,E,j), cmpii);
1151 : }
1152 :
1153 : static GEN
1154 0 : FpX_ddf_i(GEN T, GEN p)
1155 : {
1156 : GEN XP;
1157 0 : T = FpX_get_red(T, p);
1158 0 : XP = FpX_Frobenius(T, p);
1159 0 : return ddf_to_ddf2(FpX_ddf_Shoup(T, XP, p));
1160 : }
1161 :
1162 : GEN
1163 7 : FpX_ddf(GEN f, GEN p)
1164 : {
1165 7 : pari_sp av = avma;
1166 : GEN F;
1167 7 : switch(ZX_factmod_init(&f, p))
1168 : {
1169 7 : case 0: F = F2x_ddf(f);
1170 7 : F2xV_to_ZXV_inplace(gel(F,1)); break;
1171 0 : case 1: F = Flx_ddf(f,p[2]);
1172 0 : FlxV_to_ZXV_inplace(gel(F,1)); break;
1173 0 : default: F = FpX_ddf_i(f,p); break;
1174 : }
1175 7 : return gerepilecopy(av, F);
1176 : }
1177 :
1178 : static GEN Flx_simplefact_Cantor(GEN T, ulong p);
1179 : static GEN
1180 14 : FpX_simplefact_Cantor(GEN T, GEN p)
1181 : {
1182 : GEN V;
1183 : long i, l;
1184 14 : if (lgefint(p) == 3)
1185 : {
1186 0 : ulong pp = p[2];
1187 0 : return Flx_simplefact_Cantor(ZX_to_Flx(T,pp), pp);
1188 : }
1189 14 : T = FpX_get_red(T, p);
1190 14 : V = FpX_factor_Yun(get_FpX_mod(T), p); l = lg(V);
1191 28 : for (i=1; i < l; i++)
1192 14 : gel(V,i) = FpX_ddf_Shoup(gel(V,i), FpX_Frobenius(gel(V,i), p), p);
1193 14 : return vddf_to_simplefact(V, get_FpX_degree(T));
1194 : }
1195 :
1196 : static int
1197 0 : FpX_isirred_Cantor(GEN Tp, GEN p)
1198 : {
1199 0 : pari_sp av = avma;
1200 : pari_timer ti;
1201 : long n;
1202 0 : GEN T = get_FpX_mod(Tp);
1203 0 : GEN dT = FpX_deriv(T, p);
1204 : GEN XP, D;
1205 0 : if (degpol(FpX_gcd(T, dT, p)) != 0) return gc_bool(av,0);
1206 0 : n = get_FpX_degree(T);
1207 0 : T = FpX_get_red(Tp, p);
1208 0 : if (DEBUGLEVEL>=6) timer_start(&ti);
1209 0 : XP = FpX_Frobenius(T, p);
1210 0 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
1211 0 : D = FpX_ddf_Shoup(T, XP, p);
1212 0 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
1213 0 : return gc_bool(av, degpol(gel(D,n)) == n);
1214 : }
1215 :
1216 : static GEN FpX_factor_deg2(GEN f, GEN p, long d, long flag);
1217 :
1218 : /*Assume that p is large and odd*/
1219 : static GEN
1220 2518 : FpX_factor_i(GEN f, GEN pp, long flag)
1221 : {
1222 2518 : long d = degpol(f);
1223 2518 : if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
1224 924 : switch(flag)
1225 : {
1226 910 : default: return FpX_factor_Cantor(f, pp);
1227 14 : case 1: return FpX_simplefact_Cantor(f, pp);
1228 0 : case 2: return FpX_isirred_Cantor(f, pp)? gen_1: NULL;
1229 : }
1230 : }
1231 :
1232 : long
1233 0 : FpX_nbfact_Frobenius(GEN T, GEN XP, GEN p)
1234 : {
1235 0 : pari_sp av = avma;
1236 0 : long s = ddf_to_nbfact(FpX_ddf_Shoup(T, XP, p));
1237 0 : return gc_long(av,s);
1238 : }
1239 :
1240 : long
1241 0 : FpX_nbfact(GEN T, GEN p)
1242 : {
1243 0 : pari_sp av = avma;
1244 0 : GEN XP = FpX_Frobenius(T, p);
1245 0 : long n = FpX_nbfact_Frobenius(T, XP, p);
1246 0 : return gc_long(av,n);
1247 : }
1248 :
1249 : /* p > 2 */
1250 : static GEN
1251 7 : FpX_is_irred_2(GEN f, GEN p, long d)
1252 : {
1253 7 : switch(d)
1254 : {
1255 0 : case -1:
1256 0 : case 0: return NULL;
1257 0 : case 1: return gen_1;
1258 : }
1259 7 : return FpX_quad_factortype(f, p) == -1? gen_1: NULL;
1260 : }
1261 : /* p > 2 */
1262 : static GEN
1263 14 : FpX_degfact_2(GEN f, GEN p, long d)
1264 : {
1265 14 : switch(d)
1266 : {
1267 0 : case -1:retmkvec2(mkvecsmall(-1),mkvecsmall(1));
1268 0 : case 0: return trivial_fact();
1269 0 : case 1: retmkvec2(mkvecsmall(1), mkvecsmall(1));
1270 : }
1271 14 : switch(FpX_quad_factortype(f, p)) {
1272 7 : case 1: retmkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1273 7 : case -1: retmkvec2(mkvecsmall(2), mkvecsmall(1));
1274 0 : default: retmkvec2(mkvecsmall(1), mkvecsmall(2));
1275 : }
1276 : }
1277 :
1278 : GEN
1279 70 : prime_fact(GEN x) { retmkmat2(mkcolcopy(x), mkcol(gen_1)); }
1280 : GEN
1281 164231 : trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
1282 :
1283 : /* not gerepile safe */
1284 : static GEN
1285 1573 : FpX_factor_2(GEN f, GEN p, long d)
1286 : {
1287 : GEN r, s, R, S;
1288 : long v;
1289 : int sgn;
1290 1573 : switch(d)
1291 : {
1292 7 : case -1: retmkvec2(mkcol(pol_0(varn(f))), mkvecsmall(1));
1293 37 : case 0: retmkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1294 674 : case 1: retmkvec2(mkcol(f), mkvecsmall(1));
1295 : }
1296 855 : r = FpX_quad_root(f, p, 1);
1297 855 : if (!r) return mkvec2(mkcol(f), mkvecsmall(1));
1298 264 : v = varn(f);
1299 264 : s = FpX_otherroot(f, r, p);
1300 264 : if (signe(r)) r = subii(p, r);
1301 264 : if (signe(s)) s = subii(p, s);
1302 264 : sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
1303 264 : R = deg1pol_shallow(gen_1, r, v);
1304 264 : if (!sgn) return mkvec2(mkcol(R), mkvecsmall(2));
1305 162 : S = deg1pol_shallow(gen_1, s, v);
1306 162 : return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
1307 : }
1308 : static GEN
1309 1594 : FpX_factor_deg2(GEN f, GEN p, long d, long flag)
1310 : {
1311 1594 : switch(flag) {
1312 7 : case 2: return FpX_is_irred_2(f, p, d);
1313 14 : case 1: return FpX_degfact_2(f, p, d);
1314 1573 : default: return FpX_factor_2(f, p, d);
1315 : }
1316 : }
1317 :
1318 : static int
1319 447624 : F2x_quad_factortype(GEN x)
1320 447624 : { return x[2] == 7 ? -1: x[2] == 6 ? 1 :0; }
1321 :
1322 : static GEN
1323 14 : F2x_is_irred_2(GEN f, long d)
1324 14 : { return d == 1 || (d==2 && F2x_quad_factortype(f) == -1)? gen_1: NULL; }
1325 :
1326 : static GEN
1327 19620 : F2x_degfact_2(GEN f, long d)
1328 : {
1329 19620 : if (!d) return trivial_fact();
1330 19620 : if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1331 19431 : switch(F2x_quad_factortype(f)) {
1332 5964 : case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1333 6223 : case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1334 7244 : default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1335 : }
1336 : }
1337 :
1338 : static GEN
1339 1158748 : F2x_factor_2(GEN f, long d)
1340 : {
1341 1158748 : long v = f[1];
1342 1158748 : if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1343 904602 : if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1344 423428 : switch(F2x_quad_factortype(f))
1345 : {
1346 70912 : case -1: return mkvec2(mkcol(f), mkvecsmall(1));
1347 325117 : case 0: return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
1348 27423 : default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
1349 : }
1350 : }
1351 : static GEN
1352 1178379 : F2x_factor_deg2(GEN f, long d, long flag)
1353 : {
1354 1178379 : switch(flag) {
1355 14 : case 2: return F2x_is_irred_2(f, d);
1356 19620 : case 1: return F2x_degfact_2(f, d);
1357 1158745 : default: return F2x_factor_2(f, d);
1358 : }
1359 : }
1360 :
1361 : /* xt = NULL or x^(p-1)/2 mod g */
1362 : static void
1363 29715 : split_squares(struct split_t *S, GEN g, ulong p, ulong pi, GEN xt)
1364 : {
1365 29715 : ulong q = p >> 1;
1366 29715 : GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
1367 29715 : long d = degpol(a);
1368 29715 : if (d < 0)
1369 : {
1370 : ulong i;
1371 405 : split_add_done(S, (GEN)1);
1372 405 : if (!pi)
1373 1549 : for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr(i,p));
1374 : else
1375 0 : for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr_pre(i,p,pi));
1376 : } else {
1377 29310 : if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
1378 29310 : if (d)
1379 : {
1380 29261 : if (xt) xt = Flx_Fl_add(xt, p-1, p); else xt = Flx_Xnm1(g[1], q, p);
1381 29261 : a = Flx_gcd_pre(a, xt, p, pi);
1382 29261 : if (degpol(a)) split_add(S, Flx_normalize(a, p));
1383 : }
1384 : }
1385 29715 : }
1386 : static void
1387 29715 : split_nonsquares(struct split_t *S, GEN g, ulong p, ulong pi, GEN xt)
1388 : {
1389 29715 : ulong q = p >> 1;
1390 29715 : GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
1391 29715 : long d = degpol(a);
1392 29715 : if (d < 0)
1393 : {
1394 209 : ulong i, z = nonsquare_Fl(p);
1395 209 : split_add_done(S, (GEN)z);
1396 209 : if (!pi)
1397 402 : for (i = 2; i <= q; i++)
1398 193 : split_add_done(S, (GEN)Fl_mul(z, Fl_sqr(i,p), p));
1399 : else
1400 0 : for (i = 2; i <= q; i++)
1401 0 : split_add_done(S, (GEN)Fl_mul_pre(z, Fl_sqr_pre(i,p,pi), p,pi));
1402 : } else {
1403 29506 : if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
1404 29506 : if (d)
1405 : {
1406 29149 : if (xt) xt = Flx_Fl_add(xt, 1, p); else xt = Flx_Xn1(g[1], q, p);
1407 29149 : a = Flx_gcd_pre(a, xt, p, pi);
1408 29149 : if (degpol(a)) split_add(S, Flx_normalize(a, p));
1409 : }
1410 : }
1411 29715 : }
1412 : /* p > 2. f monic Flx, f(0) != 0. Add to split_t structs coprime factors
1413 : * of g = \prod_{f(a) = 0} (X - a). Return 0 when f(x) = 0 for all x in Fp* */
1414 : static int
1415 5060714 : split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p, ulong pi)
1416 : {
1417 5060714 : GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
1418 5060573 : long d = degpol(g);
1419 5060544 : if (d < 0) return 0;
1420 5059944 : if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
1421 5059946 : if (!d) return 1;
1422 5030387 : if ((p >> 4) <= (ulong)d)
1423 : { /* small p; split directly using x^((p-1)/2) +/- 1 */
1424 27122 : GEN xt = ((ulong)d < (p>>1))? Flx_rem_pre(monomial_Flx(1, p>>1, g[1]), g, p, pi)
1425 29715 : : NULL;
1426 29715 : split_squares(S, g, p, pi, xt);
1427 29715 : split_nonsquares(S, g, p, pi, xt);
1428 : } else { /* large p; use x^(p-1) - 1 directly */
1429 5000672 : a = Flxq_powu_pre(polx_Flx(f[1]), p-1, g, p, pi);
1430 4989105 : if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
1431 4989105 : a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
1432 4990249 : g = Flx_gcd_pre(g,a, p,pi);
1433 4994290 : if (degpol(g)) split_add(S, Flx_normalize(g,p));
1434 : }
1435 5028854 : return 1;
1436 : }
1437 :
1438 : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
1439 : GEN
1440 24631714 : Flx_roots_pre(GEN f, ulong p, ulong pi)
1441 : {
1442 : GEN pol;
1443 24631714 : long v = Flx_valrem(f, &f), n = degpol(f);
1444 : ulong q, PI;
1445 : struct split_t S;
1446 :
1447 24612005 : f = Flx_normalize(f, p);
1448 : /* optimization: test for small degree first */
1449 24615835 : if (n == 1)
1450 : {
1451 109298 : q = p - f[2];
1452 109298 : return v? mkvecsmall2(0, q): mkvecsmall(q);
1453 : }
1454 24506537 : PI = pi? pi: get_Fl_red(p); /* PI for Fp, pi for Fp[x] */
1455 24508724 : if (n == 2)
1456 : {
1457 19447234 : ulong r = Flx_quad_root(f, p, PI, 1), s;
1458 19419893 : if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
1459 13228605 : s = Flx_otherroot(f,r, p);
1460 13323366 : if (r < s)
1461 3344045 : return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
1462 9979321 : else if (r > s)
1463 10006363 : return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
1464 : else
1465 4832 : return v? mkvecsmall2(0, s): mkvecsmall(s);
1466 : }
1467 5061490 : if (SMALL_ULONG(p)) pi = 0; /* bilinear ops faster without Fl_*_pre */
1468 5061490 : q = p >> 1;
1469 5061490 : split_init(&S, lg(f)-1);
1470 5060712 : settyp(S.done, t_VECSMALL);
1471 5060712 : if (v) split_add_done(&S, (GEN)0);
1472 5060712 : if (! split_Flx_cut_out_roots(&S, f, p, pi))
1473 600 : return all_roots_mod_p(p, lg(S.done) == 1);
1474 5058299 : pol = polx_Flx(f[1]);
1475 5059770 : for (pol[2]=1; ; pol[2]++)
1476 5455249 : {
1477 10515019 : long j, l = lg(S.todo);
1478 10515019 : if (l == 1) { vecsmall_sort(S.done); return S.done; }
1479 5455374 : if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
1480 11395676 : for (j = 1; j < l; j++)
1481 : {
1482 5940427 : GEN b, c = gel(S.todo,j);
1483 : ulong r, s;
1484 5940427 : switch(degpol(c))
1485 : {
1486 4967576 : case 1:
1487 4967576 : split_moveto_done(&S, j, (GEN)(p - c[2]));
1488 4967793 : j--; l--; break;
1489 454487 : case 2:
1490 454487 : r = Flx_quad_root(c, p, PI, 0);
1491 454521 : if (r == p) pari_err_PRIME("polrootsmod",utoipos(p));
1492 454514 : s = Flx_otherroot(c,r, p);
1493 454515 : split_done(&S, j, (GEN)r, (GEN)s);
1494 454516 : j--; l--; break;
1495 518134 : default:
1496 518134 : b = Flxq_powu_pre(pol,q, c,p,pi); /* pol^(p-1)/2 */
1497 518281 : if (degpol(b) <= 0) continue;
1498 406158 : b = Flx_gcd_pre(c,Flx_Fl_add(b,p-1,p), p, pi);
1499 406214 : if (!degpol(b)) continue;
1500 405721 : b = Flx_normalize(b, p);
1501 405748 : c = Flx_div_pre(c,b, p,pi);
1502 405721 : split_todo(&S, j, b, c);
1503 : }
1504 : }
1505 : }
1506 : }
1507 :
1508 : GEN
1509 1372 : Flx_roots(GEN f, ulong p)
1510 : {
1511 1372 : pari_sp av = avma;
1512 : ulong pi;
1513 1372 : switch(lg(f))
1514 : {
1515 0 : case 2: pari_err_ROOTS0("Flx_roots");
1516 0 : case 3: set_avma(av); return cgetg(1, t_VECSMALL);
1517 : }
1518 1372 : if (p == 2) return Flx_root_mod_2(f);
1519 1365 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
1520 1365 : return gerepileuptoleaf(av, Flx_roots_pre(f, p, pi));
1521 : }
1522 :
1523 : /* assume x reduced mod p, monic. */
1524 : static int
1525 2477403 : Flx_quad_factortype(GEN x, ulong p)
1526 : {
1527 2477403 : ulong b = x[3], c = x[2];
1528 2477403 : return krouu(Fl_disc_bc(b, c, p), p);
1529 : }
1530 : static GEN
1531 56 : Flx_is_irred_2(GEN f, ulong p, long d)
1532 : {
1533 56 : if (!d) return NULL;
1534 56 : if (d == 1) return gen_1;
1535 56 : return Flx_quad_factortype(f, p) == -1? gen_1: NULL;
1536 : }
1537 : static GEN
1538 2501528 : Flx_degfact_2(GEN f, ulong p, long d)
1539 : {
1540 2501528 : if (!d) return trivial_fact();
1541 2501528 : if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1542 2477434 : switch(Flx_quad_factortype(f, p)) {
1543 1174673 : case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1544 1273357 : case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1545 29955 : default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1546 : }
1547 : }
1548 : /* p > 2 */
1549 : static GEN
1550 1920949 : Flx_factor_2(GEN f, ulong p, long d)
1551 : {
1552 : ulong r, s;
1553 : GEN R,S;
1554 1920949 : long v = f[1];
1555 1920949 : if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1556 1844018 : if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1557 1526359 : r = Flx_quad_root(f, p, get_Fl_red(p), 1);
1558 1526413 : if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
1559 915388 : s = Flx_otherroot(f, r, p);
1560 915383 : r = Fl_neg(r, p);
1561 915376 : s = Fl_neg(s, p);
1562 915374 : if (s < r) lswap(s,r);
1563 915374 : R = mkvecsmall3(v,r,1);
1564 915379 : if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
1565 807044 : S = mkvecsmall3(v,s,1);
1566 807052 : return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
1567 : }
1568 : static GEN
1569 4422379 : Flx_factor_deg2(GEN f, ulong p, long d, long flag)
1570 : {
1571 4422379 : switch(flag) {
1572 56 : case 2: return Flx_is_irred_2(f, p, d);
1573 2501573 : case 1: return Flx_degfact_2(f, p, d);
1574 1920750 : default: return Flx_factor_2(f, p, d);
1575 : }
1576 : }
1577 :
1578 : static GEN
1579 23226 : F2x_Berlekamp_ker(GEN u)
1580 : {
1581 23226 : pari_sp ltop=avma;
1582 23226 : long j,N = F2x_degree(u);
1583 : GEN Q;
1584 : pari_timer T;
1585 23226 : timer_start(&T);
1586 23226 : Q = F2x_matFrobenius(u);
1587 321013 : for (j=1; j<=N; j++)
1588 297787 : F2m_flip(Q,j,j);
1589 23226 : if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
1590 23226 : Q = F2m_ker_sp(Q,0);
1591 23226 : if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
1592 23226 : return gerepileupto(ltop,Q);
1593 : }
1594 : #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
1595 : static long
1596 35718 : F2x_split_Berlekamp(GEN *t)
1597 : {
1598 35718 : GEN u = *t, a, b, vker;
1599 35718 : long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
1600 :
1601 35717 : if (du == 1) return 1;
1602 27245 : if (du == 2)
1603 : {
1604 4019 : if (F2x_quad_factortype(u) == 1) /* 0 is a root: shouldn't occur */
1605 : {
1606 0 : t[0] = mkvecsmall2(sv, 2);
1607 0 : t[1] = mkvecsmall2(sv, 3);
1608 0 : return 2;
1609 : }
1610 4019 : return 1;
1611 : }
1612 :
1613 23226 : vker = F2x_Berlekamp_ker(u);
1614 23226 : lb = lgcols(vker);
1615 23239 : d = lg(vker)-1;
1616 23239 : ir = 0;
1617 : /* t[i] irreducible for i < ir, still to be treated for i < L */
1618 60555 : for (L=1; L<d; )
1619 : {
1620 : GEN pol;
1621 37330 : if (d == 2)
1622 5593 : pol = F2v_to_F2x(gel(vker,2), sv);
1623 : else
1624 : {
1625 31737 : GEN v = zero_zv(lb);
1626 31740 : v[1] = du;
1627 31740 : v[2] = random_Fl(2); /*Assume vker[1]=1*/
1628 119521 : for (i=2; i<=d; i++)
1629 87780 : if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
1630 31741 : pol = F2v_to_F2x(v, sv);
1631 : }
1632 112253 : for (i=ir; i<L && L<d; i++)
1633 : {
1634 74937 : a = t[i]; la = F2x_degree(a);
1635 74931 : if (la == 1) { set_irred(i); }
1636 74737 : else if (la == 2)
1637 : {
1638 716 : if (F2x_quad_factortype(a) == 1) /* 0 is a root: shouldn't occur */
1639 : {
1640 0 : t[i] = mkvecsmall2(sv, 2);
1641 0 : t[L] = mkvecsmall2(sv, 3); L++;
1642 : }
1643 716 : set_irred(i);
1644 : }
1645 : else
1646 : {
1647 74021 : pari_sp av = avma;
1648 : long lb;
1649 74021 : b = F2x_rem(pol, a);
1650 74014 : if (F2x_degree(b) <= 0) { set_avma(av); continue; }
1651 27034 : b = F2x_gcd(a,b); lb = F2x_degree(b);
1652 27037 : if (lb && lb < la)
1653 : {
1654 27037 : t[L] = F2x_div(a,b);
1655 27036 : t[i]= b; L++;
1656 : }
1657 0 : else set_avma(av);
1658 : }
1659 : }
1660 : }
1661 23225 : return d;
1662 : }
1663 : /* assume deg f > 2 */
1664 : static GEN
1665 35162 : F2x_Berlekamp_i(GEN f, long flag)
1666 : {
1667 35162 : long lfact, val, d = F2x_degree(f), j, k, lV;
1668 : GEN y, E, t, V;
1669 :
1670 35160 : val = F2x_valrem(f, &f);
1671 35161 : if (flag == 2 && val) return NULL;
1672 35147 : V = F2x_factor_squarefree(f); lV = lg(V);
1673 35148 : if (flag == 2 && lV > 2) return NULL;
1674 :
1675 : /* to hold factors and exponents */
1676 35078 : t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
1677 35078 : E = cgetg(d+1,t_VECSMALL);
1678 35078 : lfact = 1;
1679 35078 : if (val) {
1680 11286 : if (flag == 1) t[1] = 1; else gel(t,1) = polx_F2x(f[1]);
1681 11286 : E[1] = val; lfact++;
1682 : }
1683 :
1684 114811 : for (k=1; k<lV; k++)
1685 : {
1686 79889 : if (F2x_degree(gel(V, k))==0) continue;
1687 35717 : gel(t,lfact) = gel(V, k);
1688 35717 : d = F2x_split_Berlekamp(&gel(t,lfact));
1689 35715 : if (flag == 2 && d != 1) return NULL;
1690 35561 : if (flag == 1)
1691 47641 : for (j=0; j<d; j++) t[lfact+j] = F2x_degree(gel(t,lfact+j));
1692 97746 : for (j=0; j<d; j++) E[lfact+j] = k;
1693 35561 : lfact += d;
1694 : }
1695 34922 : if (flag == 2) return gen_1; /* irreducible */
1696 34908 : setlg(t, lfact);
1697 34908 : setlg(E, lfact); y = mkvec2(t,E);
1698 24234 : return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
1699 59144 : : sort_factor_pol(y, cmpGuGu);
1700 : }
1701 :
1702 : /* Adapted from Shoup NTL */
1703 : GEN
1704 279385 : F2x_factor_squarefree(GEN f)
1705 : {
1706 : GEN r, t, v, tv;
1707 279385 : long i, q, n = F2x_degree(f);
1708 279386 : GEN u = const_vec(n+1, pol1_F2x(f[1]));
1709 279392 : for(q = 1;;q *= 2)
1710 : {
1711 481401 : r = F2x_gcd(f, F2x_deriv(f));
1712 481391 : if (F2x_degree(r) == 0)
1713 : {
1714 230940 : gel(u, q) = f;
1715 230940 : break;
1716 : }
1717 250461 : t = F2x_div(f, r);
1718 250456 : if (F2x_degree(t) > 0)
1719 : {
1720 : long j;
1721 93649 : for(j = 1;;j++)
1722 : {
1723 206053 : v = F2x_gcd(r, t);
1724 206051 : tv = F2x_div(t, v);
1725 206055 : if (F2x_degree(tv) > 0)
1726 95639 : gel(u, j*q) = tv;
1727 206053 : if (F2x_degree(v) <= 0) break;
1728 112403 : r = F2x_div(r, v);
1729 112404 : t = v;
1730 : }
1731 93646 : if (F2x_degree(r) == 0) break;
1732 : }
1733 202009 : f = F2x_sqrt(r);
1734 : }
1735 1116672 : for (i = n; i; i--)
1736 1114042 : if (F2x_degree(gel(u,i))) break;
1737 279368 : setlg(u,i+1); return u;
1738 : }
1739 :
1740 : static GEN
1741 288072 : F2x_ddf_simple(GEN T, GEN XP)
1742 : {
1743 288072 : pari_sp av = avma, av2;
1744 : GEN f, z, Tr, X;
1745 288072 : long j, n = F2x_degree(T), v = T[1], B = n/2;
1746 288074 : if (n == 0) return cgetg(1, t_VEC);
1747 288077 : if (n == 1) return mkvec(T);
1748 133284 : z = XP; Tr = T; X = polx_F2x(v);
1749 133286 : f = const_vec(n, pol1_F2x(v));
1750 133286 : av2 = avma;
1751 348805 : for (j = 1; j <= B; j++)
1752 : {
1753 226808 : GEN u = F2x_gcd(Tr, F2x_add(z, X));
1754 226778 : if (F2x_degree(u))
1755 : {
1756 47902 : gel(f, j) = u;
1757 47902 : Tr = F2x_div(Tr, u);
1758 47904 : av2 = avma;
1759 178911 : } else z = gerepileuptoleaf(av2, z);
1760 226860 : if (!F2x_degree(Tr)) break;
1761 215551 : z = F2xq_sqr(z, Tr);
1762 : }
1763 133274 : if (F2x_degree(Tr)) gel(f, F2x_degree(Tr)) = Tr;
1764 133284 : return gerepilecopy(av, f);
1765 : }
1766 :
1767 : GEN
1768 7 : F2x_ddf(GEN T)
1769 : {
1770 : GEN XP;
1771 7 : T = F2x_get_red(T);
1772 7 : XP = F2x_Frobenius(T);
1773 7 : return F2x_ddf_to_ddf2(F2x_ddf_simple(T, XP));
1774 : }
1775 :
1776 : static GEN
1777 20163 : F2xq_frobtrace(GEN a, long d, GEN T)
1778 : {
1779 20163 : pari_sp av = avma;
1780 : long i;
1781 20163 : GEN x = a;
1782 59599 : for(i=1; i<d; i++)
1783 : {
1784 39436 : x = F2x_add(a, F2xq_sqr(x,T));
1785 39436 : if (gc_needed(av, 2))
1786 0 : x = gerepileuptoleaf(av, x);
1787 : }
1788 20163 : return x;
1789 : }
1790 :
1791 : static void
1792 29988 : F2x_edf_simple(GEN Tp, GEN XP, long d, GEN V, long idx)
1793 : {
1794 29988 : long n = F2x_degree(Tp), r = n/d;
1795 : GEN T, f, ff;
1796 29988 : if (r==1) { gel(V, idx) = Tp; return; }
1797 10092 : T = Tp;
1798 10092 : XP = F2x_rem(XP, T);
1799 : while (1)
1800 10071 : {
1801 20163 : pari_sp btop = avma;
1802 : long df;
1803 20163 : GEN g = random_F2x(n, Tp[1]);
1804 20163 : GEN t = F2xq_frobtrace(g, d, T);
1805 20163 : if (lgpol(t) == 0) continue;
1806 15281 : f = F2x_gcd(t, Tp); df = F2x_degree(f);
1807 15281 : if (df > 0 && df < n) break;
1808 5189 : set_avma(btop);
1809 : }
1810 10092 : ff = F2x_div(Tp, f);
1811 10092 : F2x_edf_simple(f, XP, d, V, idx);
1812 10092 : F2x_edf_simple(ff, XP, d, V, idx+F2x_degree(f)/d);
1813 : }
1814 :
1815 : static GEN
1816 288077 : F2x_factor_Shoup(GEN T)
1817 : {
1818 288077 : long i, n, s = 0;
1819 : GEN XP, D, V;
1820 : pari_timer ti;
1821 288077 : n = F2x_degree(T);
1822 288071 : if (DEBUGLEVEL>=6) timer_start(&ti);
1823 288071 : XP = F2x_Frobenius(T);
1824 288066 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1825 288066 : D = F2x_ddf_simple(T, XP);
1826 288075 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1827 961422 : for (i = 1; i <= n; i++)
1828 673362 : s += F2x_degree(gel(D,i))/i;
1829 288060 : V = cgetg(s+1, t_COL);
1830 961512 : for (i = 1, s = 1; i <= n; i++)
1831 : {
1832 673396 : GEN Di = gel(D,i);
1833 673396 : long ni = F2x_degree(Di), ri = ni/i;
1834 673388 : if (ni == 0) continue;
1835 324644 : if (ni == i) { gel(V, s++) = Di; continue; }
1836 9755 : F2x_edf_simple(Di, XP, i, V, s);
1837 9804 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_edf(%ld)",i);
1838 9804 : s += ri;
1839 : }
1840 288116 : return V;
1841 : }
1842 :
1843 : static GEN
1844 244240 : F2x_factor_Cantor(GEN T)
1845 : {
1846 244240 : GEN E, F, V = F2x_factor_squarefree(T);
1847 244246 : long i, j, l = lg(V);
1848 244246 : E = cgetg(l, t_VEC);
1849 244247 : F = cgetg(l, t_VEC);
1850 867696 : for (i=1, j=1; i < l; i++)
1851 623451 : if (F2x_degree(gel(V,i)))
1852 : {
1853 288077 : GEN Fj = F2x_factor_Shoup(gel(V,i));
1854 288074 : gel(F, j) = Fj;
1855 288074 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1856 288074 : j++;
1857 : }
1858 244245 : return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
1859 : }
1860 :
1861 : #if 0
1862 : static GEN
1863 : F2x_simplefact_Shoup(GEN T)
1864 : {
1865 : long i, n, s = 0, j = 1, k;
1866 : GEN XP, D, V;
1867 : pari_timer ti;
1868 : n = F2x_degree(T);
1869 : if (DEBUGLEVEL>=6) timer_start(&ti);
1870 : XP = F2x_Frobenius(T);
1871 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1872 : D = F2x_ddf_simple(T, XP);
1873 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1874 : for (i = 1; i <= n; i++)
1875 : s += F2x_degree(gel(D,i))/i;
1876 : V = cgetg(s+1, t_VECSMALL);
1877 : for (i = 1; i <= n; i++)
1878 : {
1879 : long ni = F2x_degree(gel(D,i)), ri = ni/i;
1880 : if (ni == 0) continue;
1881 : for (k = 1; k <= ri; k++)
1882 : V[j++] = i;
1883 : }
1884 : return V;
1885 : }
1886 : static GEN
1887 : F2x_simplefact_Cantor(GEN T)
1888 : {
1889 : GEN E, F, V = F2x_factor_squarefree(T);
1890 : long i, j, l = lg(V);
1891 : F = cgetg(l, t_VEC);
1892 : E = cgetg(l, t_VEC);
1893 : for (i=1, j=1; i < l; i++)
1894 : if (F2x_degree(gel(V,i)))
1895 : {
1896 : GEN Fj = F2x_simplefact_Shoup(gel(V,i));
1897 : gel(F, j) = Fj;
1898 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1899 : j++;
1900 : }
1901 : return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
1902 : }
1903 : static int
1904 : F2x_isirred_Cantor(GEN T)
1905 : {
1906 : pari_sp av = avma;
1907 : pari_timer ti;
1908 : long n;
1909 : GEN dT = F2x_deriv(T);
1910 : GEN XP, D;
1911 : if (F2x_degree(F2x_gcd(T, dT)) != 0) return gc_bool(av,0);
1912 : n = F2x_degree(T);
1913 : if (DEBUGLEVEL>=6) timer_start(&ti);
1914 : XP = F2x_Frobenius(T);
1915 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1916 : D = F2x_ddf_simple(T, XP);
1917 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1918 : return gc_bool(av, F2x_degree(gel(D,n)) == n);
1919 : }
1920 : #endif
1921 :
1922 : /* driver for Cantor factorization, assume deg f > 2; not competitive for
1923 : * flag != 0, or as deg f increases */
1924 : static GEN
1925 244239 : F2x_Cantor_i(GEN f, long flag)
1926 : {
1927 : switch(flag)
1928 : {
1929 244239 : default: return F2x_factor_Cantor(f);
1930 : #if 0
1931 : case 1: return F2x_simplefact_Cantor(f);
1932 : case 2: return F2x_isirred_Cantor(f)? gen_1: NULL;
1933 : #endif
1934 : }
1935 : }
1936 : static GEN
1937 1457744 : F2x_factor_i(GEN f, long flag)
1938 : {
1939 1457744 : long d = F2x_degree(f);
1940 1457761 : if (d <= 2) return F2x_factor_deg2(f,d,flag);
1941 254914 : return (flag == 0 && d <= 20)? F2x_Cantor_i(f, flag)
1942 534308 : : F2x_Berlekamp_i(f, flag);
1943 : }
1944 :
1945 : GEN
1946 0 : F2x_degfact(GEN f)
1947 : {
1948 0 : pari_sp av = avma;
1949 0 : GEN z = F2x_factor_i(f, 1);
1950 0 : return gerepilecopy(av, z);
1951 : }
1952 :
1953 : int
1954 238 : F2x_is_irred(GEN f) { return !!F2x_factor_i(f, 2); }
1955 :
1956 : /* Adapted from Shoup NTL */
1957 : GEN
1958 7079610 : Flx_factor_squarefree_pre(GEN f, ulong p, ulong pi)
1959 : {
1960 7079610 : long i, q, n = degpol(f);
1961 7079332 : GEN u = const_vec(n+1, pol1_Flx(f[1]));
1962 7083398 : for(q = 1;;q *= p)
1963 86369 : {
1964 7169767 : GEN t, v, tv, r = Flx_gcd_pre(f, Flx_deriv(f, p), p, pi);
1965 7163435 : if (degpol(r) == 0) { gel(u, q) = f; break; }
1966 565196 : t = Flx_div_pre(f, r, p, pi);
1967 565198 : if (degpol(t) > 0)
1968 : {
1969 : long j;
1970 485561 : for(j = 1;;j++)
1971 : {
1972 1193017 : v = Flx_gcd_pre(r, t, p, pi);
1973 1192998 : tv = Flx_div_pre(t, v, p, pi);
1974 1193002 : if (degpol(tv) > 0)
1975 725861 : gel(u, j*q) = Flx_normalize(tv, p);
1976 1193005 : if (degpol(v) <= 0) break;
1977 707443 : r = Flx_div_pre(r, v, p, pi);
1978 707456 : t = v;
1979 : }
1980 485561 : if (degpol(r) == 0) break;
1981 : }
1982 86369 : f = Flx_normalize(Flx_deflate(r, p), p);
1983 : }
1984 30024852 : for (i = n; i; i--)
1985 30026143 : if (degpol(gel(u,i))) break;
1986 7077370 : setlg(u,i+1); return u;
1987 : }
1988 : GEN
1989 335001 : Flx_factor_squarefree(GEN f, ulong p)
1990 335001 : { return Flx_factor_squarefree_pre(f, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1991 :
1992 : long
1993 3451 : Flx_ispower(GEN f, ulong k, ulong p, GEN *pt_r)
1994 : {
1995 3451 : pari_sp av = avma;
1996 : ulong lc, pi;
1997 : GEN F;
1998 3451 : long i, n = degpol(f), v = f[1], l;
1999 3451 : if (n % k) return 0;
2000 3451 : lc = Fl_sqrtn(Flx_lead(f), k, p, NULL);
2001 3451 : if (lc == ULONG_MAX) { av = avma; return 0; }
2002 3451 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2003 3451 : F = Flx_factor_squarefree_pre(f, p, pi); l = lg(F)-1;
2004 38325 : for (i = 1; i <= l; i++)
2005 34895 : if (i%k && degpol(gel(F,i))) return gc_long(av,0);
2006 3430 : if (pt_r)
2007 : {
2008 3430 : GEN r = Fl_to_Flx(lc, v), s = pol1_Flx(v);
2009 38304 : for(i = l; i >= 1; i--)
2010 : {
2011 34874 : if (i%k) continue;
2012 12453 : s = Flx_mul_pre(s, gel(F,i), p, pi);
2013 12453 : r = Flx_mul_pre(r, s, p, pi);
2014 : }
2015 3430 : *pt_r = gerepileuptoleaf(av, r);
2016 0 : } else set_avma(av);
2017 3430 : return 1;
2018 : }
2019 :
2020 : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
2021 : static GEN
2022 8993232 : Flx_ddf_Shoup(GEN T, GEN XP, ulong p, ulong pi)
2023 : {
2024 8993232 : pari_sp av = avma;
2025 : GEN b, g, h, F, f, Tr, xq;
2026 : long i, j, n, v, bo, ro;
2027 : long B, l, m;
2028 : pari_timer ti;
2029 8993232 : n = get_Flx_degree(T); v = get_Flx_var(T);
2030 8994977 : if (n == 0) return cgetg(1, t_VEC);
2031 8951408 : if (n == 1) return mkvec(get_Flx_mod(T));
2032 8652118 : B = n/2;
2033 8652118 : l = usqrt(B);
2034 8651870 : m = (B+l-1)/l;
2035 8651870 : T = Flx_get_red(T, p);
2036 8650788 : b = cgetg(l+2, t_VEC);
2037 8651650 : gel(b, 1) = polx_Flx(v);
2038 8652661 : gel(b, 2) = XP;
2039 8652661 : bo = brent_kung_optpow(n, l-1, 1);
2040 8653352 : ro = l<=1 ? 0:(bo-1)/(l-1) + ((n-1)/bo);
2041 8653352 : if (DEBUGLEVEL>=7) timer_start(&ti);
2042 8653352 : if (expu(p) <= ro)
2043 871556 : for (i = 3; i <= l+1; i++)
2044 483341 : gel(b, i) = Flxq_powu_pre(gel(b, i-1), p, T, p, pi);
2045 : else
2046 : {
2047 8265157 : xq = Flxq_powers_pre(gel(b, 2), bo, T, p, pi);
2048 8263371 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq baby");
2049 8975311 : for (i = 3; i <= l+1; i++)
2050 711837 : gel(b, i) = Flx_FlxqV_eval_pre(gel(b, i-1), xq, T, p, pi);
2051 : }
2052 8651689 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: baby");
2053 8651689 : xq = Flxq_powers_pre(gel(b, l+1), brent_kung_optpow(n, m-1, 1), T, p, pi);
2054 8650088 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq giant");
2055 8650088 : g = cgetg(m+1, t_VEC);
2056 8652389 : gel(g, 1) = gel(xq, 2);
2057 14017194 : for(i = 2; i <= m; i++)
2058 5365632 : gel(g, i) = Flx_FlxqV_eval_pre(gel(g, i-1), xq, T, p, pi);
2059 8651562 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: giant");
2060 8651562 : h = cgetg(m+1, t_VEC);
2061 22669917 : for (j = 1; j <= m; j++)
2062 : {
2063 14018619 : pari_sp av = avma;
2064 14018619 : GEN gj = gel(g, j);
2065 14018619 : GEN e = Flx_sub(gj, gel(b, 1), p);
2066 17858873 : for (i = 2; i <= l; i++)
2067 3850657 : e = Flxq_mul_pre(e, Flx_sub(gj, gel(b, i), p), T, p, pi);
2068 14008216 : gel(h, j) = gerepileupto(av, e);
2069 : }
2070 8651298 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: diff");
2071 8651298 : Tr = get_Flx_mod(T);
2072 8652157 : F = cgetg(m+1, t_VEC);
2073 22665307 : for (j = 1; j <= m; j++)
2074 : {
2075 14014832 : GEN u = Flx_gcd_pre(Tr, gel(h, j), p, pi);
2076 14011764 : if (degpol(u))
2077 : {
2078 6470387 : u = Flx_normalize(u, p);
2079 6470491 : Tr = Flx_div_pre(Tr, u, p, pi);
2080 : }
2081 14010991 : gel(F, j) = u;
2082 : }
2083 8650475 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: F");
2084 8650475 : f = const_vec(n, pol1_Flx(v));
2085 22667590 : for (j = 1; j <= m; j++)
2086 : {
2087 14018301 : GEN e = gel(F, j);
2088 14888835 : for (i=l-1; i >= 0; i--)
2089 : {
2090 14888479 : GEN u = Flx_gcd_pre(e, Flx_sub(gel(g, j), gel(b, i+1), p), p, pi);
2091 14882926 : if (degpol(u))
2092 : {
2093 6569258 : gel(f, l*j-i) = u;
2094 6569258 : e = Flx_div_pre(e, u, p, pi);
2095 : }
2096 14881424 : if (!degpol(e)) break;
2097 : }
2098 : }
2099 8649289 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: f");
2100 8649289 : if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
2101 8649551 : return gerepilecopy(av, f);
2102 : }
2103 :
2104 : static void
2105 404666 : Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, ulong pi, GEN V, long idx)
2106 : {
2107 404666 : long n = degpol(Tp), r = n/d;
2108 : GEN T, f, ff;
2109 : ulong p2;
2110 404668 : if (r==1) { gel(V, idx) = Tp; return; }
2111 172969 : p2 = p>>1;
2112 172969 : T = Flx_get_red_pre(Tp, p, pi);
2113 172971 : XP = Flx_rem_pre(XP, T, p, pi);
2114 : while (1)
2115 20109 : {
2116 193080 : pari_sp btop = avma;
2117 : long i;
2118 193080 : GEN g = random_Flx(n, Tp[1], p);
2119 193080 : GEN t = gel(Flxq_auttrace_pre(mkvec2(XP, g), d, T, p, pi), 2);
2120 193080 : if (lgpol(t) == 0) continue;
2121 426922 : for(i=1; i<=10; i++)
2122 : {
2123 411770 : pari_sp btop2 = avma;
2124 411770 : GEN R = Flxq_powu_pre(Flx_Fl_add(t, random_Fl(p), p), p2, T, p, pi);
2125 411764 : f = Flx_gcd_pre(Flx_Fl_add(R, p-1, p), Tp, p, pi);
2126 411767 : if (degpol(f) > 0 && degpol(f) < n) break;
2127 238797 : set_avma(btop2);
2128 : }
2129 188123 : if (degpol(f) > 0 && degpol(f) < n) break;
2130 15150 : set_avma(btop);
2131 : }
2132 172971 : f = Flx_normalize(f, p);
2133 172970 : ff = Flx_div_pre(Tp, f, p, pi);
2134 172969 : Flx_edf_simple(f, XP, d, p, pi, V, idx);
2135 172971 : Flx_edf_simple(ff, XP, d, p, pi, V, idx+degpol(f)/d);
2136 : }
2137 : static void
2138 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, ulong pi, GEN V, long idx);
2139 :
2140 : static void
2141 1352105 : Flx_edf_rec(GEN T, GEN XP, GEN hp, GEN t, long d, ulong p, ulong pi,
2142 : GEN V, long idx)
2143 : {
2144 : pari_sp av;
2145 1352105 : GEN Tp = get_Flx_mod(T);
2146 1352106 : long n = degpol(hp), vT = Tp[1];
2147 : GEN u1, u2, f1, f2;
2148 1352105 : ulong p2 = p>>1;
2149 : GEN R, h;
2150 1352105 : h = Flx_get_red_pre(hp, p, pi);
2151 1352099 : t = Flx_rem_pre(t, T, p, pi);
2152 1352032 : av = avma;
2153 : do
2154 : {
2155 2277073 : set_avma(av);
2156 2277091 : R = Flxq_powu_pre(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p, pi);
2157 2276926 : u1 = Flx_gcd_pre(Flx_Fl_add(R, p-1, p), hp, p, pi);
2158 2277070 : } while (degpol(u1)==0 || degpol(u1)==n);
2159 1352053 : f1 = Flx_gcd_pre(Flx_Flxq_eval_pre(u1, t, T, p, pi), Tp, p, pi);
2160 1352045 : f1 = Flx_normalize(f1, p);
2161 1352083 : u2 = Flx_div_pre(hp, u1, p, pi);
2162 1352105 : f2 = Flx_div_pre(Tp, f1, p, pi);
2163 1352093 : if (degpol(u1)==1)
2164 : {
2165 1035681 : if (degpol(f1)==d)
2166 1021041 : gel(V, idx) = f1;
2167 : else
2168 14634 : Flx_edf(f1, XP, d, p, pi, V, idx);
2169 : }
2170 : else
2171 316454 : Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, pi, V, idx);
2172 1352133 : idx += degpol(f1)/d;
2173 1352107 : if (degpol(u2)==1)
2174 : {
2175 1030876 : if (degpol(f2)==d)
2176 1017131 : gel(V, idx) = f2;
2177 : else
2178 13746 : Flx_edf(f2, XP, d, p, pi, V, idx);
2179 : }
2180 : else
2181 321270 : Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, pi, V, idx);
2182 1352147 : }
2183 :
2184 : static void
2185 714442 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, ulong pi, GEN V, long idx)
2186 : {
2187 714442 : long n = degpol(Tp), r = n/d, vT = Tp[1];
2188 : GEN T, h, t;
2189 : pari_timer ti;
2190 714442 : if (r==1) { gel(V, idx) = Tp; return; }
2191 714442 : T = Flx_get_red_pre(Tp, p, pi);
2192 714440 : XP = Flx_rem_pre(XP, T, p, pi);
2193 714432 : if (DEBUGLEVEL>=7) timer_start(&ti);
2194 : do
2195 : {
2196 735963 : GEN g = random_Flx(n, vT, p);
2197 735973 : t = gel(Flxq_auttrace_pre(mkvec2(XP, g), d, T, p, pi), 2);
2198 735979 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
2199 735979 : h = Flxq_minpoly_pre(t, T, p, pi);
2200 735964 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
2201 735964 : } while (degpol(h) <= 1);
2202 714434 : Flx_edf_rec(T, XP, h, t, d, p, pi, V, idx);
2203 : }
2204 :
2205 : static GEN
2206 1519239 : Flx_factor_Shoup(GEN T, ulong p, ulong pi)
2207 : {
2208 1519239 : long i, n, s = 0, e = expu(p);
2209 : GEN XP, D, V;
2210 : pari_timer ti;
2211 1519235 : n = get_Flx_degree(T);
2212 1519231 : T = Flx_get_red_pre(T, p, pi);
2213 1519220 : if (DEBUGLEVEL>=6) timer_start(&ti);
2214 1519220 : XP = Flx_Frobenius_pre(T, p, pi);
2215 1519202 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2216 1519202 : D = Flx_ddf_Shoup(T, XP, p, pi);
2217 1519255 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2218 1519255 : s = ddf_to_nbfact(D);
2219 1519245 : V = cgetg(s+1, t_COL);
2220 8673154 : for (i = 1, s = 1; i <= n; i++)
2221 : {
2222 7153913 : GEN Di = gel(D,i);
2223 7153913 : long ni = degpol(Di), ri = ni/i;
2224 7153902 : if (ni == 0) continue;
2225 1979233 : Di = Flx_normalize(Di, p);
2226 1979273 : if (ni == i) { gel(V, s++) = Di; continue; }
2227 744781 : if (ri <= e*expu(e))
2228 686060 : Flx_edf(Di, XP, i, p, pi, V, s);
2229 : else
2230 58727 : Flx_edf_simple(Di, XP, i, p, pi, V, s);
2231 744781 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
2232 744782 : s += ri;
2233 : }
2234 1519241 : return V;
2235 : }
2236 :
2237 : static GEN
2238 1452094 : Flx_factor_Cantor(GEN T, ulong p)
2239 : {
2240 1452094 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2241 1452094 : GEN E, F, V = Flx_factor_squarefree_pre(get_Flx_mod(T), p, pi);
2242 1452098 : long i, j, l = lg(V);
2243 1452098 : F = cgetg(l, t_VEC);
2244 1452105 : E = cgetg(l, t_VEC);
2245 3353310 : for (i=1, j=1; i < l; i++)
2246 1901197 : if (degpol(gel(V,i)))
2247 : {
2248 1519241 : GEN Fj = Flx_factor_Shoup(gel(V,i), p, pi);
2249 1519240 : gel(F, j) = Fj;
2250 1519240 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
2251 1519242 : j++;
2252 : }
2253 1452113 : return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
2254 : }
2255 :
2256 : GEN
2257 0 : Flx_ddf_pre(GEN T, ulong p, ulong pi)
2258 : {
2259 : GEN XP;
2260 0 : T = Flx_get_red_pre(T, p, pi);
2261 0 : XP = Flx_Frobenius_pre(T, p, pi);
2262 0 : return ddf_to_ddf2(Flx_ddf_Shoup(T, XP, p, pi));
2263 : }
2264 : GEN
2265 0 : Flx_ddf(GEN T, ulong p)
2266 0 : { return Flx_ddf_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2267 :
2268 : static GEN
2269 5291110 : Flx_simplefact_Cantor(GEN T, ulong p)
2270 : {
2271 5291110 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2272 : long i, l;
2273 : GEN V;
2274 5291110 : T = Flx_get_red_pre(T, p, pi);
2275 5289714 : V = Flx_factor_squarefree_pre(get_Flx_mod(T), p, pi); l = lg(V);
2276 10666753 : for (i=1; i < l; i++)
2277 5372580 : gel(V,i) = Flx_ddf_Shoup(gel(V,i), Flx_Frobenius_pre(gel(V,i), p,pi), p,pi);
2278 5294173 : return vddf_to_simplefact(V, get_Flx_degree(T));
2279 : }
2280 :
2281 : static int
2282 1078 : Flx_isirred_Cantor(GEN Tp, ulong p)
2283 : {
2284 1078 : pari_sp av = avma;
2285 : pari_timer ti;
2286 1078 : GEN T = get_Flx_mod(Tp), dT = Flx_deriv(T, p), XP, D;
2287 1078 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2288 : long n;
2289 1078 : if (degpol(Flx_gcd_pre(T, dT, p, pi)) != 0) return gc_bool(av,0);
2290 791 : n = get_Flx_degree(T);
2291 791 : T = Flx_get_red_pre(Tp, p, pi);
2292 791 : if (DEBUGLEVEL>=6) timer_start(&ti);
2293 791 : XP = Flx_Frobenius_pre(T, p, pi);
2294 791 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2295 791 : D = Flx_ddf_Shoup(T, XP, p, pi);
2296 791 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2297 791 : return gc_bool(av, degpol(gel(D,n)) == n);
2298 : }
2299 :
2300 : /* f monic */
2301 : static GEN
2302 11238907 : Flx_factor_i(GEN f, ulong pp, long flag)
2303 : {
2304 : long d;
2305 11238907 : if (pp==2) { /*We need to handle 2 specially */
2306 73697 : GEN F = F2x_factor_i(Flx_to_F2x(f),flag);
2307 73698 : if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
2308 73698 : return F;
2309 : }
2310 11165210 : d = degpol(f);
2311 11165707 : if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
2312 6743723 : switch(flag)
2313 : {
2314 1452092 : default: return Flx_factor_Cantor(f, pp);
2315 5290553 : case 1: return Flx_simplefact_Cantor(f, pp);
2316 1078 : case 2: return Flx_isirred_Cantor(f, pp)? gen_1: NULL;
2317 : }
2318 : }
2319 :
2320 : GEN
2321 7787993 : Flx_degfact(GEN f, ulong p)
2322 : {
2323 7787993 : pari_sp av = avma;
2324 7787993 : GEN z = Flx_factor_i(Flx_normalize(f,p),p,1);
2325 7790803 : return gerepilecopy(av, z);
2326 : }
2327 :
2328 : /* T must be squarefree mod p*/
2329 : GEN
2330 1449389 : Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
2331 : {
2332 : GEN XP, D;
2333 : pari_timer ti;
2334 1449389 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2335 1449389 : long i, s, n = get_Flx_degree(T);
2336 1449346 : GEN V = const_vecsmall(n, 0);
2337 1449479 : pari_sp av = avma;
2338 1449479 : T = Flx_get_red_pre(T, p, pi);
2339 1449444 : if (DEBUGLEVEL>=6) timer_start(&ti);
2340 1449444 : XP = Flx_Frobenius_pre(T, p, pi);
2341 1449413 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2342 1449413 : D = Flx_ddf_Shoup(T, XP, p, pi);
2343 1449658 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2344 7452958 : for (i = 1, s = 0; i <= n; i++) { V[i] = degpol(gel(D,i))/i; s += V[i]; }
2345 1449614 : *nb = s; set_avma(av); return V;
2346 : }
2347 :
2348 : long
2349 652405 : Flx_nbfact_Frobenius_pre(GEN T, GEN XP, ulong p, ulong pi)
2350 : {
2351 652405 : pari_sp av = avma;
2352 652405 : long s = ddf_to_nbfact(Flx_ddf_Shoup(T, XP, p, pi));
2353 652415 : return gc_long(av,s);
2354 : }
2355 : long
2356 0 : Flx_nbfact_Frobenius(GEN T, GEN XP, ulong p)
2357 0 : { return Flx_nbfact_Frobenius_pre(T, XP, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2358 :
2359 : /* T must be squarefree mod p*/
2360 : long
2361 652404 : Flx_nbfact_pre(GEN T, ulong p, ulong pi)
2362 : {
2363 652404 : pari_sp av = avma;
2364 652404 : GEN XP = Flx_Frobenius_pre(T, p, pi);
2365 652404 : long n = Flx_nbfact_Frobenius_pre(T, XP, p, pi);
2366 652413 : return gc_long(av,n);
2367 : }
2368 : long
2369 652405 : Flx_nbfact(GEN T, ulong p)
2370 652405 : { return Flx_nbfact_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2371 :
2372 : int
2373 1057 : Flx_is_irred(GEN f, ulong p)
2374 : {
2375 1057 : pari_sp av = avma;
2376 1057 : f = Flx_normalize(f,p);
2377 1057 : return gc_bool(av, !!Flx_factor_i(f,p,2));
2378 : }
2379 :
2380 : /* Use this function when you think f is reducible, and that there are lots of
2381 : * factors. If you believe f has few factors, use FpX_nbfact(f,p)==1 instead */
2382 : int
2383 112 : FpX_is_irred(GEN f, GEN p)
2384 : {
2385 112 : pari_sp av = avma;
2386 : int z;
2387 112 : switch(ZX_factmod_init(&f,p))
2388 : {
2389 28 : case 0: z = !!F2x_factor_i(f,2); break;
2390 77 : case 1: z = !!Flx_factor_i(f,p[2],2); break;
2391 7 : default: z = !!FpX_factor_i(f,p,2); break;
2392 : }
2393 112 : return gc_bool(av,z);
2394 : }
2395 : GEN
2396 47502 : FpX_degfact(GEN f, GEN p) {
2397 47502 : pari_sp av = avma;
2398 : GEN F;
2399 47502 : switch(ZX_factmod_init(&f,p))
2400 : {
2401 462 : case 0: F = F2x_factor_i(f,1); break;
2402 47012 : case 1: F = Flx_factor_i(f,p[2],1); break;
2403 28 : default: F = FpX_factor_i(f,p,1); break;
2404 : }
2405 47502 : return gerepilecopy(av, F);
2406 : }
2407 :
2408 : #if 0
2409 : /* set x <-- x + c*y mod p */
2410 : /* x is not required to be normalized.*/
2411 : static void
2412 : Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
2413 : {
2414 : long i, lx, ly;
2415 : ulong *x=(ulong *)gx;
2416 : ulong *y=(ulong *)gy;
2417 : if (!c) return;
2418 : lx = lg(gx);
2419 : ly = lg(gy);
2420 : if (lx<ly) pari_err_BUG("lx<ly in Flx_addmul_inplace");
2421 : if (SMALL_ULONG(p))
2422 : for (i=2; i<ly; i++) x[i] = (x[i] + c*y[i]) % p;
2423 : else
2424 : for (i=2; i<ly; i++) x[i] = Fl_add(x[i], Fl_mul(c,y[i],p),p);
2425 : }
2426 : #endif
2427 :
2428 : GEN
2429 4095881 : FpX_factor(GEN f, GEN p)
2430 : {
2431 4095881 : pari_sp av = avma;
2432 : GEN F;
2433 4095881 : switch(ZX_factmod_init(&f, p))
2434 : {
2435 1368020 : case 0: F = F2x_factor_i(f,0);
2436 1368080 : F2xV_to_ZXV_inplace(gel(F,1)); break;
2437 2725420 : case 1: F = Flx_factor_i(f,p[2],0);
2438 2725476 : FlxV_to_ZXV_inplace(gel(F,1)); break;
2439 2419 : default: F = FpX_factor_i(f,p,0); break;
2440 : }
2441 4096042 : return gerepilecopy(av, F);
2442 : }
2443 :
2444 : GEN
2445 677931 : Flx_factor(GEN f, ulong p)
2446 : {
2447 677931 : pari_sp av = avma;
2448 677931 : return gerepilecopy(av, Flx_factor_i(Flx_normalize(f,p),p,0));
2449 : }
2450 : GEN
2451 15309 : F2x_factor(GEN f)
2452 : {
2453 15309 : pari_sp av = avma;
2454 15309 : return gerepilecopy(av, F2x_factor_i(f,0));
2455 : }
|
