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 4351068 : ZX_factmod_init(GEN *F, GEN p)
35 : {
36 4351068 : if (lgefint(p) == 3)
37 : {
38 4348649 : ulong pp = p[2];
39 4348649 : if (pp == 2) { *F = ZX_to_F2x(*F); return 0; }
40 2913069 : *F = ZX_to_Flx(*F, pp);
41 2913139 : if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
42 2913138 : return 1;
43 : }
44 2419 : *F = FpX_red(*F, p);
45 2425 : if (lg(*F) > 3) *F = FpX_normalize(*F, p);
46 2425 : return 2;
47 : }
48 : static GEN
49 665904 : ZX_rootmod_init(GEN F, GEN p)
50 665904 : { 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 2174 : Flx_Xnm1(long sv, long n, ulong p)
71 : {
72 2174 : GEN t = cgetg(n+3, t_VECSMALL);
73 : long i;
74 2174 : t[1] = sv;
75 2174 : t[2] = p - 1;
76 6896 : for (i = 3; i <= n+1; i++) t[i] = 0;
77 2174 : t[i] = 1; return t;
78 : }
79 : /* X^n + 1 */
80 : static GEN
81 2069 : Flx_Xn1(long sv, long n, ulong p)
82 : {
83 2069 : GEN t = cgetg(n+3, t_VECSMALL);
84 : long i;
85 : (void) p;
86 2069 : t[1] = sv;
87 2069 : t[2] = 1;
88 6708 : for (i = 3; i <= n+1; i++) t[i] = 0;
89 2069 : t[i] = 1; return t;
90 : }
91 :
92 : /* assume lg(f) > 3 */
93 : static GEN
94 93051 : Flx_root_mod_2(GEN f)
95 : {
96 93051 : long i, n = lg(f)-1, c = f[2];
97 93051 : int z0 = !c;
98 93051 : c ^= 1; /* c = f[2] + f[n] mod 2, we know f[n] is odd */
99 167467 : for (i=3; i < n; i++) c ^= f[i];
100 : /* c = 0 iff f(1) = 0 (mod 2) */
101 93051 : if (z0) return c? mkvecsmall(0): mkvecsmall2(0, 1);
102 13931 : 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 16287741 : 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 431992 : FpX_roots(GEN f, GEN p)
123 : {
124 431992 : pari_sp av = avma;
125 431992 : GEN y; f = ZX_rootmod_init(f, p);
126 431992 : switch(lg(f))
127 : {
128 14 : case 2: pari_err_ROOTS0("FpX_roots");
129 48747 : case 3: return cgetg(1,t_COL);
130 : }
131 383231 : if (typ(f) == t_VECSMALL)
132 : {
133 374991 : ulong pp = p[2];
134 374991 : if (pp == 2)
135 93044 : y = Flx_root_mod_2(f);
136 : else
137 : {
138 281947 : if (!odd(pp)) pari_err_PRIME("FpX_roots", p);
139 281947 : y = Flx_roots_pre(f, pp, SMALL_ULONG(pp)? 0: get_Fl_red(pp));
140 : }
141 374984 : y = Flc_to_ZC(y);
142 : }
143 : else
144 8240 : y = FpX_roots_i(f, p);
145 383218 : 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 14529 : FpX_quad_root(GEN x, GEN p, int unknown)
159 : {
160 14529 : GEN s, D, b = gel(x,3), c = gel(x,2);
161 :
162 14529 : if (absequaliu(p, 2)) {
163 0 : if (!signe(b)) return c;
164 0 : return signe(c)? NULL: gen_1;
165 : }
166 14529 : D = subii(sqri(b), shifti(c,2));
167 14529 : D = remii(D,p);
168 14529 : if (unknown && kronecker(D,p) == -1) return NULL;
169 :
170 13946 : s = Fp_sqrt(D,p);
171 : /* p is not prime, go on and give e.g. maxord a chance to recover */
172 13946 : if (!s) return NULL;
173 13938 : return Fp_halve(Fp_sub(s,b, p), p);
174 : }
175 : static GEN
176 10010 : FpX_otherroot(GEN x, GEN r, GEN p)
177 10010 : { 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 27434627 : Fl_disc_bc(ulong b, ulong c, ulong p)
182 27434627 : { 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 24952015 : Flx_quad_root(GEN x, ulong p, ulong pi, int unknown)
186 : {
187 24952015 : ulong s, b = x[3], c = x[2];
188 24952015 : ulong D = Fl_disc_bc(b, c, p);
189 24916040 : if (unknown && krouu(D,p) == -1) return p;
190 16482433 : s = Fl_sqrt_pre(D, p, pi);
191 16541613 : if (s==~0UL) return p;
192 16541600 : return Fl_halve(Fl_sub(s,b, p), p);
193 : }
194 : static ulong
195 14743967 : Flx_otherroot(GEN x, ulong r, ulong p)
196 14743967 : { 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 5058977 : split_init(struct split_t *S, long max)
203 : {
204 5058977 : S->todo = vectrunc_init(max);
205 5058486 : S->done = vectrunc_init(max);
206 5057997 : }
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 5421889 : 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 5422978 : split_moveto_done(struct split_t *S, long i, GEN t)
224 : {
225 5422978 : long n = lg(S->todo)-1;
226 5422978 : vectrunc_append(S->done, t);
227 5422919 : if (n) gel(S->todo,i) = gel(S->todo, n);
228 5422919 : setlg(S->todo, n);
229 :
230 5422843 : }
231 : /* append t to done */
232 : static void
233 511745 : split_add_done(struct split_t *S, GEN t)
234 511745 : { vectrunc_append(S->done, t); }
235 : /* split todo[i] into a and b */
236 : static void
237 418576 : split_todo(struct split_t *S, long i, GEN a, GEN b)
238 : {
239 418576 : gel(S->todo, i) = a;
240 418576 : split_add(S, b);
241 418569 : }
242 : /* split todo[i] into a and b, moved to done */
243 : static void
244 467730 : split_done(struct split_t *S, long i, GEN a, GEN b)
245 : {
246 467730 : split_moveto_done(S, i, a);
247 467737 : split_add_done(S, b);
248 467734 : }
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 417717 : Flx_cubic_root(GEN ff, ulong p, ulong pi)
322 : {
323 417717 : GEN f = Flx_normalize(ff,p);
324 417715 : 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 417715 : if (pi)
327 : {
328 417715 : t = Fl_mul_pre(a, p3, p, pi);
329 417720 : t2 = Fl_sqr_pre(t, p, pi);
330 417717 : A = Fl_sub(b, Fl_triple(t2, p), p);
331 417718 : B = Fl_sub(c, Fl_mul_pre(t, Fl_add(A, t2, p), p, pi), p);
332 417718 : A3 = Fl_mul_pre(A, p3, p, pi);
333 417719 : 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 417718 : A33 = Fl_powu_pre(A3, 3, p, pi);
345 417706 : D = Fl_add(B2, Fl_double(Fl_double(A33, p), p), p);
346 417703 : S = Fl_neg(B,p);
347 417702 : P = Fl_neg(A3,p);
348 417701 : if (krouu(D,p) >= 0)
349 : {
350 343739 : ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
351 343745 : ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
352 343747 : if (p%3==2) /* 1 solutions */
353 138827 : vS1 = Fl_powu_pre(S1, p - p3, p, pi);
354 : else
355 : {
356 204920 : vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
357 204918 : if (vS1==~0UL) return p; /*0 solutions*/
358 : /*3 solutions*/
359 : }
360 223950 : if (!P) return Fl_sub(vS1, t, p);
361 96156 : vS2 = pi? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): Fl_div(P, vS1, p);
362 96159 : return Fl_sub(Fl_add(vS1,vS2, p), t, p);
363 : }
364 : else
365 : {
366 73973 : pari_sp av = avma;
367 73973 : GEN S1 = mkvecsmall2(Fl_halve(S, p), (p + 1UL) >> 1);
368 73975 : GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
369 : ulong Sa;
370 73976 : if (!vS1) return p; /*0 solutions, p%3==2*/
371 73976 : Sa = vS1[1];
372 73976 : if (p%3==1) /*1 solutions*/
373 : {
374 29298 : ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
375 29297 : if (Fa!=P) Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
376 : }
377 73975 : set_avma(av);
378 73975 : return Fl_sub(Fl_double(Sa,p),t,p);
379 : }
380 : }
381 :
382 : /* Assume f is normalized */
383 : static GEN
384 119 : FpX_cubic_root(GEN ff, GEN p)
385 : {
386 119 : GEN f = FpX_normalize(ff,p);
387 119 : GEN a = gel(f,4), b = gel(f,3), c = gel(f,2);
388 119 : ulong pm3 = umodiu(p,3);
389 28 : GEN p3 = pm3==1 ? diviuexact(addiu(shifti(p,1),1),3)
390 119 : : diviuexact(addiu(p,1),3);
391 119 : GEN t = Fp_mul(a, p3, p), t2 = Fp_sqr(t, p);
392 119 : GEN A = Fp_sub(b, Fp_mulu(t2, 3, p), p);
393 119 : GEN B = Fp_addmul(c, t, Fp_sub(shifti(t2, 1), b, p), p);
394 119 : GEN A3 = Fp_mul(A, p3, p), A33 = Fp_powu(A3, 3, p);
395 119 : GEN S = Fp_neg(B,p), P = Fp_neg(A3,p);
396 119 : GEN D = Fp_add(Fp_sqr(S, p), shifti(A33, 2), p);
397 119 : 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 91 : pari_sp av = avma;
415 91 : GEN T = deg2pol_shallow(gen_1, gen_0, negi(D), 0);
416 91 : GEN S1 = deg1pol_shallow(Fp_halve(gen_1, p), Fp_halve(S, p), 0);
417 91 : GEN vS1 = FpXQ_sqrtn(S1, utoi(3), T, p, NULL);
418 : GEN Sa;
419 91 : if (!vS1) return p; /*0 solutions, p%3==2*/
420 91 : Sa = gel(vS1,2);
421 91 : 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 91 : set_avma(av);
428 91 : 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 4144737 : 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 4144737 : if (Flx_val(f)) return 0;
441 4143093 : da = degpol(f); f = Flx_normalize(f, p);
442 4142622 : if (da == 1) return Fl_neg(f[2], p);
443 4126488 : PI = pi? pi: get_Fl_red(p); /* PI for Fp, pi for Fp[x] */
444 4127472 : switch(da)
445 : {
446 3392233 : case 2: return Flx_quad_root(f, p, PI, 1);
447 404011 : case 3: if (p>3) return Flx_cubic_root(f, p, PI); /*FALL THROUGH*/
448 : }
449 337906 : if (SMALL_ULONG(p)) pi = 0; /* bilinear ops faster without Fl_*_pre */
450 337906 : if (!fl)
451 : {
452 300630 : a = Flxq_powu_pre(polx_Flx(f[1]), p - 1, f,p,pi);
453 300515 : if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
454 300515 : a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
455 300512 : a = Flx_gcd_pre(f,a, p, pi);
456 37276 : } else a = f;
457 337873 : da = degpol(a);
458 337867 : if (!da) return p;
459 241973 : a = Flx_normalize(a,p);
460 :
461 242010 : q = p >> 1;
462 242010 : pol = polx_Flx(f[1]);
463 358940 : for(pol[2] = 1;; pol[2]++)
464 : {
465 358940 : if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
466 358947 : switch(da)
467 : {
468 152375 : case 1: return Fl_neg(a[2], p);
469 69287 : case 2: return Flx_quad_root(a, p, PI, 0);
470 20384 : case 3: if (p>3) return Flx_cubic_root(a, p, PI); /*FALL THROUGH*/
471 : default: {
472 116901 : GEN b = Flxq_powu_pre(pol,q, a,p,pi);
473 : long db;
474 116878 : if (degpol(b) <= 0) continue;
475 111039 : b = Flx_gcd_pre(a,Flx_Fl_add(b,p-1,p), p, pi);
476 111068 : db = degpol(b); if (!db) continue;
477 111070 : b = Flx_normalize(b, p);
478 111068 : if (db <= (da >> 1)) {
479 68368 : a = b;
480 68368 : da = db;
481 : } else {
482 42700 : a = Flx_div_pre(a,b, p, pi);
483 42712 : da -= db;
484 : }
485 : }
486 : }
487 : }
488 : }
489 : ulong
490 4106044 : Flx_oneroot_pre(GEN f, ulong p, ulong pi)
491 4106044 : { return Flx_oneroot_pre_i(f, p, pi, 0); }
492 : ulong
493 38689 : Flx_oneroot_split_pre(GEN f, ulong p, ulong pi)
494 38689 : { return Flx_oneroot_pre_i(f, p, pi, 1); }
495 :
496 : /* assume p > 3 prime */
497 : static GEN
498 5183 : FpX_oneroot_i(GEN f, GEN p)
499 : {
500 : GEN pol, pol0, a, q;
501 : long da;
502 :
503 5183 : if (ZX_val(f)) return gen_0;
504 4812 : f = FpX_normalize(f, p);
505 4812 : switch(degpol(f))
506 : {
507 723 : case 1: return subii(p, gel(f,2));
508 3900 : case 2: return FpX_quad_root(f, p, 1);
509 119 : 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 2492 : Flx_oneroot(GEN f, ulong p)
551 : {
552 2492 : pari_sp av = avma;
553 2492 : switch(lg(f))
554 : {
555 0 : case 2: return 0;
556 0 : case 3: return p;
557 : }
558 2492 : if (p == 2) return Flx_oneroot_mod_2(f);
559 2492 : 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 233912 : FpX_oneroot(GEN f, GEN p)
578 : {
579 233912 : pari_sp av = avma;
580 233912 : f = ZX_rootmod_init(f, p);
581 233912 : switch(lg(f))
582 : {
583 0 : case 2: set_avma(av); return gen_0;
584 0 : case 3: return gc_NULL(av);
585 : }
586 233912 : if (typ(f) == t_VECSMALL)
587 : {
588 228729 : ulong r, pp = p[2];
589 228729 : if (pp == 2)
590 91 : r = Flx_oneroot_mod_2(f);
591 : else
592 228638 : r = Flx_oneroot_pre(f, pp, SMALL_ULONG(pp)? 0: get_Fl_red(pp));
593 228729 : set_avma(av);
594 228729 : return (r == pp)? NULL: utoi(r);
595 : }
596 5183 : f = FpX_oneroot_i(f, p);
597 5183 : if (!f) return gc_NULL(av);
598 5183 : 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 112 : good_root_of_unity(GEN p, long deg, long deg_factor, long *pt_n)
608 : {
609 112 : pari_sp ltop = avma;
610 : GEN pm, factn, power, base, zeta;
611 : long n;
612 :
613 112 : pm = subis (p, 1ul);
614 252 : for (n = deg / 2 / deg_factor + 1; !dvdiu (pm, n); n--);
615 112 : factn = Z_factor(stoi(n));
616 112 : power = diviuexact (pm, n);
617 112 : base = gen_1;
618 : do {
619 161 : base = addis (base, 1l);
620 161 : zeta = Fp_pow (base, power, p);
621 : }
622 161 : while (!equaliu (Fp_order (zeta, factn, p), n));
623 112 : *pt_n = n;
624 112 : return gerepileuptoint (ltop, zeta);
625 : }
626 :
627 : GEN
628 1036 : FpX_oneroot_split(GEN fact, GEN p)
629 : {
630 1036 : pari_sp av = avma;
631 : long n, deg_f, i, dmin;
632 : GEN prim, expo, minfactor, xplusa, zeta, xpow;
633 1036 : fact = FpX_normalize(fact, p);
634 1036 : deg_f = degpol(fact);
635 1036 : if (deg_f <= 3) return FpX_oneroot(fact, p);
636 98 : minfactor = fact; /* factor of minimal degree found so far */
637 98 : dmin = degpol(minfactor);
638 98 : xplusa = pol_x(varn(fact));
639 210 : 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 112 : fact = minfactor; deg_f = dmin;
644 112 : zeta = gen_1;
645 112 : prim = good_root_of_unity(p, deg_f, 1, &n);
646 112 : expo = diviuexact(subiu(p, 1), n);
647 : /* update X+a, avoid a=0 */
648 112 : gel (xplusa, 2) = addis (gel (xplusa, 2), 1);
649 112 : xpow = FpXQ_pow (xplusa, expo, fact, p);
650 245 : for (i = 0; i < n; i++)
651 : {
652 182 : GEN tmp = FpX_gcd(FpX_Fp_sub(xpow, zeta, p), fact, p);
653 182 : long dtmp = degpol(tmp);
654 182 : if (dtmp > 0 && dtmp < deg_f)
655 : {
656 119 : fact = FpX_div(fact, tmp, p); deg_f = degpol(fact);
657 119 : if (dtmp < dmin)
658 : {
659 112 : minfactor = FpX_normalize (tmp, p);
660 112 : dmin = dtmp;
661 112 : if (dmin == 1 || dmin <= (2 * deg_f) / n - 1)
662 : /* stop early to avoid too many gcds */
663 : break;
664 : }
665 : }
666 133 : zeta = Fp_mul (zeta, prim, p);
667 : }
668 : }
669 98 : 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 1884687 : FE_concat(GEN F, GEN E, long l)
682 : {
683 1884687 : setlg(E,l); E = shallowconcat1(E);
684 1884685 : 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 5314845 : vddf_to_simplefact(GEN V, long d)
716 : {
717 : GEN E, F;
718 5314845 : long i, j, c, l = lg(V);
719 5314845 : F = cgetg(d+1, t_VECSMALL);
720 5313915 : E = cgetg(d+1, t_VECSMALL);
721 10710210 : for (i = c = 1; i < l; i++)
722 : {
723 5396752 : GEN Vi = gel(V,i);
724 5396752 : long l = lg(Vi);
725 27209173 : for (j = 1; j < l; j++)
726 : {
727 21812558 : long k, n = degpol(gel(Vi,j)) / j;
728 33028803 : for (k = 1; k <= n; k++) { uel(F,c) = j; uel(E,c) = i; c++; }
729 : }
730 : }
731 5313458 : setlg(F,c);
732 5313417 : setlg(E,c);
733 5313451 : 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 234324 : FpX_split_part(GEN f, GEN p)
739 : {
740 234324 : long n = degpol(f);
741 234324 : GEN z, X = pol_x(varn(f));
742 234324 : if (n <= 1) return f;
743 232356 : f = FpX_red(f, p);
744 232352 : z = FpX_sub(FpX_Frobenius(f, p), X, p);
745 232355 : 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 138259 : FpX_nbroots(GEN f, GEN p)
752 : {
753 138259 : pari_sp av = avma;
754 138259 : GEN z = FpX_split_part(f, p);
755 138259 : 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 4391284 : Flx_nbroots(GEN f, ulong p)
794 : {
795 4391284 : long n = degpol(f);
796 : ulong pi;
797 4391283 : pari_sp av = avma;
798 : GEN z;
799 4391283 : if (n <= 1) return n;
800 4377556 : if (n == 2)
801 : {
802 : ulong D;
803 580879 : if (p==2) return (f[2]==0) + (f[2]!=f[3]);
804 432243 : D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
805 432242 : return 1 + krouu(D,p);
806 : }
807 3796677 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
808 3796677 : z = Flx_sub(Flx_Frobenius_pre(f, p, pi), polx_Flx(f[1]), p);
809 3796662 : z = Flx_gcd_pre(z, f, p, pi);
810 3796666 : 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 1018 : 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 1018 : n = get_FpX_degree(T); v = get_FpX_var(T);
863 1018 : if (n == 0) return cgetg(1, t_VEC);
864 1018 : if (n == 1) return mkvec(get_FpX_mod(T));
865 842 : B = n/2;
866 842 : l = usqrt(B);
867 842 : m = (B+l-1)/l;
868 842 : T = FpX_get_red(T, p);
869 842 : b = cgetg(l+2, t_VEC);
870 842 : gel(b, 1) = pol_x(v);
871 842 : gel(b, 2) = XP;
872 842 : if (DEBUGLEVEL>=7) timer_start(&ti);
873 842 : xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1), T, p);
874 842 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq baby");
875 1068 : for (i = 3; i <= l+1; i++)
876 226 : gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
877 842 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: baby");
878 842 : xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1), T, p);
879 842 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq giant");
880 842 : g = cgetg(m+1, t_VEC);
881 842 : gel(g, 1) = gel(xq, 2);
882 1827 : for(i = 2; i <= m; i++) gel(g, i) = FpX_FpXQV_eval(gel(g, i-1), xq, T, p);
883 842 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: giant");
884 842 : h = cgetg(m+1, t_VEC);
885 2669 : for (j = 1; j <= m; j++)
886 : {
887 1827 : pari_sp av = avma;
888 1827 : GEN gj = gel(g,j), e = FpX_sub(gj, gel(b,1), p);
889 3054 : for (i = 2; i <= l; i++) e = FpXQ_mul(e, FpX_sub(gj, gel(b,i), p), T, p);
890 1827 : gel(h,j) = gerepileupto(av, e);
891 : }
892 842 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: diff");
893 842 : Tr = get_FpX_mod(T);
894 842 : F = cgetg(m+1, t_VEC);
895 2669 : for (j = 1; j <= m; j++)
896 : {
897 1827 : GEN u = FpX_gcd(Tr, gel(h,j), p);
898 1827 : if (degpol(u))
899 : {
900 475 : u = FpX_normalize(u, p);
901 475 : Tr = FpX_div(Tr, u, p);
902 : }
903 1827 : gel(F,j) = u;
904 : }
905 842 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: F");
906 842 : f = const_vec(n, pol_1(v));
907 2669 : for (j = 1; j <= m; j++)
908 : {
909 1827 : GEN e = gel(F, j);
910 1899 : for (i=l-1; i >= 0; i--)
911 : {
912 1899 : GEN u = FpX_gcd(e, FpX_sub(gel(g, j), gel(b, i+1), p), p);
913 1899 : if (degpol(u))
914 : {
915 505 : u = FpX_normalize(u, p);
916 505 : gel(f, l*j-i) = u;
917 505 : e = FpX_div(e, u, p);
918 : }
919 1899 : if (!degpol(e)) break;
920 : }
921 : }
922 842 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: f");
923 842 : if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
924 842 : 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 1036 : 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 1036 : GEN Tp = get_FpX_mod(T);
966 1036 : long n = degpol(hp), vT = varn(Tp), ct = 0;
967 : GEN u1, u2, f1, f2, R, h;
968 1036 : h = FpX_get_red(hp, p);
969 1036 : t = FpX_rem(t, T, p);
970 1036 : av = avma;
971 : do
972 : {
973 1729 : set_avma(av);
974 1729 : R = FpXQ_pow(deg1pol(gen_1, randomi(p), vT), p2, h, p);
975 1729 : u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
976 1729 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_rec",p);
977 1729 : } while (degpol(u1)==0 || degpol(u1)==n);
978 1036 : f1 = FpX_gcd(FpX_FpXQ_eval(u1, t, T, p), Tp, p);
979 1036 : f1 = FpX_normalize(f1, p);
980 1036 : u2 = FpX_div(hp, u1, p);
981 1036 : f2 = FpX_div(Tp, f1, p);
982 1036 : if (degpol(u1)==1)
983 786 : gel(V, idx) = f1;
984 : else
985 250 : FpX_edf_rec(FpX_get_red(f1, p), u1, t, d, p2, p, V, idx);
986 1036 : idx += degpol(f1)/d;
987 1036 : if (degpol(u2)==1)
988 681 : gel(V, idx) = f2;
989 : else
990 355 : FpX_edf_rec(FpX_get_red(f2, p), u2, t, d, p2, p, V, idx);
991 1036 : }
992 :
993 : /* assume Tp a squarefree product of r > 1 irred. factors of degree d */
994 : static void
995 431 : FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
996 : {
997 431 : long n = degpol(Tp), r = n/d, vT = varn(Tp), ct = 0;
998 : GEN T, h, t;
999 : pari_timer ti;
1000 :
1001 431 : T = FpX_get_red(Tp, p);
1002 431 : XP = FpX_rem(XP, T, p);
1003 431 : if (DEBUGLEVEL>=7) timer_start(&ti);
1004 : do
1005 : {
1006 431 : GEN g = random_FpX(n, vT, p);
1007 431 : t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
1008 431 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
1009 431 : h = FpXQ_minpoly(t, T, p);
1010 431 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
1011 431 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf",p);
1012 431 : } while (degpol(h) != r);
1013 431 : FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
1014 431 : }
1015 :
1016 : static GEN
1017 1004 : FpX_factor_Shoup(GEN T, GEN p)
1018 : {
1019 1004 : long i, n, s = 0;
1020 : GEN XP, D, V;
1021 1004 : long e = expi(p);
1022 : pari_timer ti;
1023 1004 : n = get_FpX_degree(T);
1024 1004 : T = FpX_get_red(T, p);
1025 1004 : if (DEBUGLEVEL>=6) timer_start(&ti);
1026 1004 : XP = FpX_Frobenius(T, p);
1027 1004 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
1028 1004 : D = FpX_ddf_Shoup(T, XP, p);
1029 1004 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
1030 1004 : s = ddf_to_nbfact(D);
1031 1004 : V = cgetg(s+1, t_COL);
1032 7249 : for (i = 1, s = 1; i <= n; i++)
1033 : {
1034 6245 : GEN Di = gel(D,i);
1035 6245 : long ni = degpol(Di), ri = ni/i;
1036 6245 : if (ni == 0) continue;
1037 1122 : Di = FpX_normalize(Di, p);
1038 1122 : if (ni == i) { gel(V, s++) = Di; continue; }
1039 431 : if (ri <= e*expu(e))
1040 431 : FpX_edf(Di, XP, i, p, V, s);
1041 : else
1042 0 : FpX_edf_simple(Di, XP, i, p, V, s);
1043 431 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
1044 431 : s += ri;
1045 : }
1046 1004 : return V;
1047 : }
1048 :
1049 : long
1050 2315722 : ddf_to_nbfact(GEN D)
1051 : {
1052 2315722 : long l = lg(D), i, s = 0;
1053 14647045 : for(i = 1; i < l; i++) s += degpol(gel(D,i))/i;
1054 2315716 : return s;
1055 : }
1056 :
1057 : /* Yun algorithm: Assume p > degpol(T) */
1058 : static GEN
1059 1695 : FpX_factor_Yun(GEN T, GEN p)
1060 : {
1061 1695 : long n = degpol(T), i = 1;
1062 1695 : GEN a, b, c, d = FpX_deriv(T, p);
1063 1695 : GEN V = cgetg(n+1,t_VEC);
1064 1695 : a = FpX_gcd(T, d, p);
1065 1695 : if (degpol(a) == 0) return mkvec(T);
1066 588 : b = FpX_div(T, a, p);
1067 : do
1068 : {
1069 2646 : c = FpX_div(d, a, p);
1070 2646 : d = FpX_sub(c, FpX_deriv(b, p), p);
1071 2646 : a = FpX_normalize(FpX_gcd(b, d, p), p);
1072 2646 : gel(V, i++) = a;
1073 2646 : b = FpX_div(b, a, p);
1074 2646 : } while (degpol(b));
1075 588 : setlg(V, i); return V;
1076 : }
1077 : GEN
1078 347291 : FpX_factor_squarefree(GEN T, GEN p)
1079 : {
1080 347291 : if (lgefint(p)==3)
1081 : {
1082 346557 : ulong pp = (ulong)p[2];
1083 346557 : GEN u = Flx_factor_squarefree(ZX_to_Flx(T,pp), pp);
1084 346557 : return FlxV_to_ZXV(u);
1085 : }
1086 734 : return FpX_factor_Yun(T, p);
1087 : }
1088 :
1089 : GEN
1090 226464 : FpX_roots_mult(GEN T, long n, GEN p)
1091 : {
1092 226464 : pari_sp av = avma;
1093 226464 : GEN V = FpX_factor_squarefree(T, p), W;
1094 226464 : long l = lg(V), i;
1095 226464 : if (l <= n) { set_avma(av); return cgetg(1,t_COL); }
1096 40836 : W = cgetg(l-n+1,t_VEC);
1097 129754 : for (i = n; i < l; i++)
1098 88918 : gel(W,i-n+1) = FpX_roots(gel(V,i), p);
1099 40836 : return gerepileupto(av, sort(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 905 : FpX_factor_Cantor(GEN T, GEN p)
1138 : {
1139 905 : GEN E, F, V = FpX_factor_Yun(T, p);
1140 905 : long i, j, l = lg(V);
1141 905 : F = cgetg(l, t_VEC);
1142 905 : E = cgetg(l, t_VEC);
1143 2428 : for (i=1, j=1; i < l; i++)
1144 1523 : if (degpol(gel(V,i)))
1145 : {
1146 1004 : GEN Fj = FpX_factor_Shoup(gel(V,i), p);
1147 1004 : gel(F, j) = Fj;
1148 1004 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1149 1004 : j++;
1150 : }
1151 905 : 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 2425 : FpX_factor_i(GEN f, GEN pp, long flag)
1222 : {
1223 2425 : long d = degpol(f);
1224 2425 : if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
1225 919 : switch(flag)
1226 : {
1227 905 : 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 350032 : trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
1283 :
1284 : /* not gerepile safe */
1285 : static GEN
1286 1485 : FpX_factor_2(GEN f, GEN p, long d)
1287 : {
1288 : GEN r, s, R, S;
1289 : long v;
1290 : int sgn;
1291 1485 : switch(d)
1292 : {
1293 7 : case -1: retmkvec2(mkcol(pol_0(varn(f))), mkvecsmall(1));
1294 37 : case 0: retmkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1295 583 : case 1: retmkvec2(mkcol(f), mkvecsmall(1));
1296 : }
1297 858 : r = FpX_quad_root(f, p, 1);
1298 858 : if (!r) return mkvec2(mkcol(f), mkvecsmall(1));
1299 274 : v = varn(f);
1300 274 : s = FpX_otherroot(f, r, p);
1301 274 : if (signe(r)) r = subii(p, r);
1302 274 : if (signe(s)) s = subii(p, s);
1303 274 : sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
1304 274 : R = deg1pol_shallow(gen_1, r, v);
1305 274 : if (!sgn) return mkvec2(mkcol(R), mkvecsmall(2));
1306 162 : S = deg1pol_shallow(gen_1, s, v);
1307 162 : return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
1308 : }
1309 : static GEN
1310 1506 : FpX_factor_deg2(GEN f, GEN p, long d, long flag)
1311 : {
1312 1506 : 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 1485 : default: return FpX_factor_2(f, p, d);
1316 : }
1317 : }
1318 :
1319 : static int
1320 452911 : F2x_quad_factortype(GEN x)
1321 452911 : { 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 19676 : F2x_degfact_2(GEN f, long d)
1329 : {
1330 19676 : if (!d) return trivial_fact();
1331 19676 : if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1332 19480 : switch(F2x_quad_factortype(f)) {
1333 5957 : case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1334 6244 : case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1335 7279 : default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1336 : }
1337 : }
1338 :
1339 : static GEN
1340 1168109 : F2x_factor_2(GEN f, long d)
1341 : {
1342 1168109 : long v = f[1];
1343 1168109 : if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1344 912226 : if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1345 428579 : switch(F2x_quad_factortype(f))
1346 : {
1347 71643 : case -1: return mkvec2(mkcol(f), mkvecsmall(1));
1348 329249 : case 0: return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
1349 27704 : default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
1350 : }
1351 : }
1352 : static GEN
1353 1187801 : F2x_factor_deg2(GEN f, long d, long flag)
1354 : {
1355 1187801 : switch(flag) {
1356 14 : case 2: return F2x_is_irred_2(f, d);
1357 19676 : case 1: return F2x_degfact_2(f, d);
1358 1168111 : default: return F2x_factor_2(f, d);
1359 : }
1360 : }
1361 :
1362 : /* xt = NULL or x^(p-1)/2 mod g */
1363 : static void
1364 30100 : split_squares(struct split_t *S, GEN g, ulong p, ulong pi, GEN xt)
1365 : {
1366 30100 : ulong q = p >> 1;
1367 30100 : GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
1368 30100 : long d = degpol(a);
1369 30100 : if (d < 0)
1370 : {
1371 : ulong i;
1372 414 : split_add_done(S, (GEN)1);
1373 414 : if (!pi)
1374 1576 : for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr(i,p));
1375 : else
1376 0 : for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr_pre(i,p,pi));
1377 : } else {
1378 29686 : if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
1379 29686 : if (d)
1380 : {
1381 29597 : if (xt) xt = Flx_Fl_add(xt, p-1, p); else xt = Flx_Xnm1(g[1], q, p);
1382 29597 : a = Flx_gcd_pre(a, xt, p, pi);
1383 29597 : if (degpol(a)) split_add(S, Flx_normalize(a, p));
1384 : }
1385 : }
1386 30100 : }
1387 : static void
1388 30100 : split_nonsquares(struct split_t *S, GEN g, ulong p, ulong pi, GEN xt)
1389 : {
1390 30100 : ulong q = p >> 1;
1391 30100 : GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
1392 30100 : long d = degpol(a);
1393 30100 : if (d < 0)
1394 : {
1395 244 : ulong i, z = nonsquare_Fl(p);
1396 244 : split_add_done(S, (GEN)z);
1397 244 : if (!pi)
1398 512 : for (i = 2; i <= q; i++)
1399 268 : split_add_done(S, (GEN)Fl_mul(z, Fl_sqr(i,p), p));
1400 : else
1401 0 : for (i = 2; i <= q; i++)
1402 0 : split_add_done(S, (GEN)Fl_mul_pre(z, Fl_sqr_pre(i,p,pi), p,pi));
1403 : } else {
1404 29856 : if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
1405 29856 : if (d)
1406 : {
1407 29492 : if (xt) xt = Flx_Fl_add(xt, 1, p); else xt = Flx_Xn1(g[1], q, p);
1408 29492 : a = Flx_gcd_pre(a, xt, p, pi);
1409 29492 : if (degpol(a)) split_add(S, Flx_normalize(a, p));
1410 : }
1411 : }
1412 30101 : }
1413 : /* p > 2. f monic Flx, f(0) != 0. Add to split_t structs coprime factors
1414 : * of g = \prod_{f(a) = 0} (X - a). Return 0 when f(x) = 0 for all x in Fp* */
1415 : static int
1416 5049779 : split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p, ulong pi)
1417 : {
1418 5049779 : GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
1419 5049725 : long d = degpol(g);
1420 5049694 : if (d < 0) return 0;
1421 5049094 : if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
1422 5049097 : if (!d) return 1;
1423 5015009 : if ((p >> 4) <= (ulong)d)
1424 : { /* small p; split directly using x^((p-1)/2) +/- 1 */
1425 27423 : GEN xt = ((ulong)d < (p>>1))? Flx_rem_pre(monomial_Flx(1, p>>1, g[1]), g, p, pi)
1426 30100 : : NULL;
1427 30100 : split_squares(S, g, p, pi, xt);
1428 30100 : split_nonsquares(S, g, p, pi, xt);
1429 : } else { /* large p; use x^(p-1) - 1 directly */
1430 4984909 : a = Flxq_powu_pre(polx_Flx(f[1]), p-1, g, p, pi);
1431 4972203 : if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
1432 4972203 : a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
1433 4973551 : g = Flx_gcd_pre(g,a, p,pi);
1434 4977664 : if (degpol(g)) split_add(S, Flx_normalize(g,p));
1435 : }
1436 5013810 : return 1;
1437 : }
1438 :
1439 : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
1440 : GEN
1441 24663918 : Flx_roots_pre(GEN f, ulong p, ulong pi)
1442 : {
1443 : GEN pol;
1444 24663918 : long v = Flx_valrem(f, &f), n = degpol(f);
1445 : ulong q, PI;
1446 : struct split_t S;
1447 :
1448 24660823 : f = Flx_normalize(f, p);
1449 : /* optimization: test for small degree first */
1450 24671965 : if (n == 1)
1451 : {
1452 115190 : q = p - f[2];
1453 115190 : return v? mkvecsmall2(0, q): mkvecsmall(q);
1454 : }
1455 24556775 : PI = pi? pi: get_Fl_red(p); /* PI for Fp, pi for Fp[x] */
1456 24560800 : if (n == 2)
1457 : {
1458 19509505 : ulong r = Flx_quad_root(f, p, PI, 1), s;
1459 19522075 : if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
1460 13308884 : s = Flx_otherroot(f,r, p);
1461 13352171 : if (r < s)
1462 3355046 : return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
1463 9997125 : else if (r > s)
1464 10031070 : return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
1465 : else
1466 4886 : return v? mkvecsmall2(0, s): mkvecsmall(s);
1467 : }
1468 5051295 : if (SMALL_ULONG(p)) pi = 0; /* bilinear ops faster without Fl_*_pre */
1469 5051295 : q = p >> 1;
1470 5051295 : split_init(&S, lg(f)-1);
1471 5049780 : settyp(S.done, t_VECSMALL);
1472 5049780 : if (v) split_add_done(&S, (GEN)0);
1473 5049780 : if (! split_Flx_cut_out_roots(&S, f, p, pi))
1474 600 : return all_roots_mod_p(p, lg(S.done) == 1);
1475 5047746 : pol = polx_Flx(f[1]);
1476 5049061 : for (pol[2]=1; ; pol[2]++)
1477 5446549 : {
1478 10495610 : long j, l = lg(S.todo);
1479 10495610 : if (l == 1) { vecsmall_sort(S.done); return S.done; }
1480 5446815 : if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
1481 11384814 : for (j = 1; j < l; j++)
1482 : {
1483 5938265 : GEN b, c = gel(S.todo,j);
1484 : ulong r, s;
1485 5938265 : switch(degpol(c))
1486 : {
1487 4950645 : case 1:
1488 4950645 : split_moveto_done(&S, j, (GEN)(p - c[2]));
1489 4950796 : j--; l--; break;
1490 461478 : case 2:
1491 461478 : r = Flx_quad_root(c, p, PI, 0);
1492 461534 : if (r == p) pari_err_PRIME("polrootsmod",utoipos(p));
1493 461527 : s = Flx_otherroot(c,r, p);
1494 461518 : split_done(&S, j, (GEN)r, (GEN)s);
1495 461521 : j--; l--; break;
1496 525784 : default:
1497 525784 : b = Flxq_powu_pre(pol,q, c,p,pi); /* pol^(p-1)/2 */
1498 526006 : if (degpol(b) <= 0) continue;
1499 412796 : b = Flx_gcd_pre(c,Flx_Fl_add(b,p-1,p), p, pi);
1500 412882 : if (!degpol(b)) continue;
1501 412314 : b = Flx_normalize(b, p);
1502 412360 : c = Flx_div_pre(c,b, p,pi);
1503 412333 : split_todo(&S, j, b, c);
1504 : }
1505 : }
1506 : }
1507 : }
1508 :
1509 : GEN
1510 1372 : Flx_roots(GEN f, ulong p)
1511 : {
1512 1372 : pari_sp av = avma;
1513 : ulong pi;
1514 1372 : switch(lg(f))
1515 : {
1516 0 : case 2: pari_err_ROOTS0("Flx_roots");
1517 0 : case 3: set_avma(av); return cgetg(1, t_VECSMALL);
1518 : }
1519 1372 : if (p == 2) return Flx_root_mod_2(f);
1520 1365 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
1521 1365 : return gerepileuptoleaf(av, Flx_roots_pre(f, p, pi));
1522 : }
1523 :
1524 : /* assume x reduced mod p, monic. */
1525 : static int
1526 2487543 : Flx_quad_factortype(GEN x, ulong p)
1527 : {
1528 2487543 : ulong b = x[3], c = x[2];
1529 2487543 : return krouu(Fl_disc_bc(b, c, p), p);
1530 : }
1531 : static GEN
1532 56 : Flx_is_irred_2(GEN f, ulong p, long d)
1533 : {
1534 56 : if (!d) return NULL;
1535 56 : if (d == 1) return gen_1;
1536 56 : return Flx_quad_factortype(f, p) == -1? gen_1: NULL;
1537 : }
1538 : static GEN
1539 2513596 : Flx_degfact_2(GEN f, ulong p, long d)
1540 : {
1541 2513596 : if (!d) return trivial_fact();
1542 2513596 : if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1543 2487556 : switch(Flx_quad_factortype(f, p)) {
1544 1179484 : case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1545 1278105 : case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1546 29915 : default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1547 : }
1548 : }
1549 : /* p > 2 */
1550 : static GEN
1551 1933283 : Flx_factor_2(GEN f, ulong p, long d)
1552 : {
1553 : ulong r, s;
1554 : GEN R,S;
1555 1933283 : long v = f[1];
1556 1933283 : if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1557 1856074 : if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1558 1535417 : r = Flx_quad_root(f, p, get_Fl_red(p), 1);
1559 1535491 : if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
1560 922420 : s = Flx_otherroot(f, r, p);
1561 922416 : r = Fl_neg(r, p);
1562 922413 : s = Fl_neg(s, p);
1563 922410 : if (s < r) lswap(s,r);
1564 922410 : R = mkvecsmall3(v,r,1);
1565 922411 : if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
1566 811735 : S = mkvecsmall3(v,s,1);
1567 811734 : return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
1568 : }
1569 : static GEN
1570 4446844 : Flx_factor_deg2(GEN f, ulong p, long d, long flag)
1571 : {
1572 4446844 : switch(flag) {
1573 56 : case 2: return Flx_is_irred_2(f, p, d);
1574 2513596 : case 1: return Flx_degfact_2(f, p, d);
1575 1933192 : default: return Flx_factor_2(f, p, d);
1576 : }
1577 : }
1578 :
1579 : static GEN
1580 23319 : F2x_Berlekamp_ker(GEN u)
1581 : {
1582 23319 : pari_sp ltop=avma;
1583 23319 : long j,N = F2x_degree(u);
1584 : GEN Q;
1585 : pari_timer T;
1586 23318 : timer_start(&T);
1587 23319 : Q = F2x_matFrobenius(u);
1588 318461 : for (j=1; j<=N; j++)
1589 295142 : F2m_flip(Q,j,j);
1590 23319 : if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
1591 23319 : Q = F2m_ker_sp(Q,0);
1592 23319 : if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
1593 23319 : return gerepileupto(ltop,Q);
1594 : }
1595 : #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
1596 : static long
1597 35818 : F2x_split_Berlekamp(GEN *t)
1598 : {
1599 35818 : GEN u = *t, a, b, vker;
1600 35818 : long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
1601 :
1602 35818 : if (du == 1) return 1;
1603 27354 : if (du == 2)
1604 : {
1605 4035 : if (F2x_quad_factortype(u) == 1) /* 0 is a root: shouldn't occur */
1606 : {
1607 0 : t[0] = mkvecsmall2(sv, 2);
1608 0 : t[1] = mkvecsmall2(sv, 3);
1609 0 : return 2;
1610 : }
1611 4035 : return 1;
1612 : }
1613 :
1614 23319 : vker = F2x_Berlekamp_ker(u);
1615 23318 : lb = lgcols(vker);
1616 23332 : d = lg(vker)-1;
1617 23332 : ir = 0;
1618 : /* t[i] irreducible for i < ir, still to be treated for i < L */
1619 60169 : for (L=1; L<d; )
1620 : {
1621 : GEN pol;
1622 36851 : if (d == 2)
1623 5549 : pol = F2v_to_F2x(gel(vker,2), sv);
1624 : else
1625 : {
1626 31302 : GEN v = zero_zv(lb);
1627 31313 : v[1] = du;
1628 31313 : v[2] = random_Fl(2); /*Assume vker[1]=1*/
1629 117926 : for (i=2; i<=d; i++)
1630 86610 : if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
1631 31316 : pol = F2v_to_F2x(v, sv);
1632 : }
1633 110577 : for (i=ir; i<L && L<d; i++)
1634 : {
1635 73740 : a = t[i]; la = F2x_degree(a);
1636 73728 : if (la == 1) { set_irred(i); }
1637 73500 : else if (la == 2)
1638 : {
1639 793 : if (F2x_quad_factortype(a) == 1) /* 0 is a root: shouldn't occur */
1640 : {
1641 0 : t[i] = mkvecsmall2(sv, 2);
1642 0 : t[L] = mkvecsmall2(sv, 3); L++;
1643 : }
1644 793 : set_irred(i);
1645 : }
1646 : else
1647 : {
1648 72707 : pari_sp av = avma;
1649 : long lb;
1650 72707 : b = F2x_rem(pol, a);
1651 72689 : if (F2x_degree(b) <= 0) { set_avma(av); continue; }
1652 26630 : b = F2x_gcd(a,b); lb = F2x_degree(b);
1653 26632 : if (lb && lb < la)
1654 : {
1655 26632 : t[L] = F2x_div(a,b);
1656 26635 : t[i]= b; L++;
1657 : }
1658 0 : else set_avma(av);
1659 : }
1660 : }
1661 : }
1662 23318 : return d;
1663 : }
1664 : /* assume deg f > 2 */
1665 : static GEN
1666 35374 : F2x_Berlekamp_i(GEN f, long flag)
1667 : {
1668 35374 : long lfact, val, d = F2x_degree(f), j, k, lV;
1669 : GEN y, E, t, V;
1670 :
1671 35375 : val = F2x_valrem(f, &f);
1672 35376 : if (flag == 2 && val) return NULL;
1673 35362 : V = F2x_factor_squarefree(f); lV = lg(V);
1674 35361 : if (flag == 2 && lV > 2) return NULL;
1675 :
1676 : /* to hold factors and exponents */
1677 35291 : t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
1678 35291 : E = cgetg(d+1,t_VECSMALL);
1679 35292 : lfact = 1;
1680 35292 : if (val) {
1681 11267 : if (flag == 1) t[1] = 1; else gel(t,1) = polx_F2x(f[1]);
1682 11267 : E[1] = val; lfact++;
1683 : }
1684 :
1685 117182 : for (k=1; k<lV; k++)
1686 : {
1687 82045 : if (F2x_degree(gel(V, k))==0) continue;
1688 35818 : gel(t,lfact) = gel(V, k);
1689 35818 : d = F2x_split_Berlekamp(&gel(t,lfact));
1690 35817 : if (flag == 2 && d != 1) return NULL;
1691 35663 : if (flag == 1)
1692 47763 : for (j=0; j<d; j++) t[lfact+j] = F2x_degree(gel(t,lfact+j));
1693 97548 : for (j=0; j<d; j++) E[lfact+j] = k;
1694 35663 : lfact += d;
1695 : }
1696 35137 : if (flag == 2) return gen_1; /* irreducible */
1697 35123 : setlg(t, lfact);
1698 35123 : setlg(E, lfact); y = mkvec2(t,E);
1699 24333 : return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
1700 59457 : : sort_factor_pol(y, cmpGuGu);
1701 : }
1702 :
1703 : /* Adapted from Shoup NTL */
1704 : GEN
1705 336990 : F2x_factor_squarefree(GEN f)
1706 : {
1707 : GEN r, t, v, tv;
1708 336990 : long i, q, n = F2x_degree(f);
1709 336990 : GEN u = const_vec(n+1, pol1_F2x(f[1]));
1710 336991 : for(q = 1;;q *= 2)
1711 : {
1712 548789 : r = F2x_gcd(f, F2x_deriv(f));
1713 548783 : if (F2x_degree(r) == 0)
1714 : {
1715 242539 : gel(u, q) = f;
1716 242539 : break;
1717 : }
1718 306247 : t = F2x_div(f, r);
1719 306247 : if (F2x_degree(t) > 0)
1720 : {
1721 : long j;
1722 146854 : for(j = 1;;j++)
1723 : {
1724 352547 : v = F2x_gcd(r, t);
1725 352548 : tv = F2x_div(t, v);
1726 352545 : if (F2x_degree(tv) > 0)
1727 148818 : gel(u, j*q) = tv;
1728 352546 : if (F2x_degree(v) <= 0) break;
1729 205691 : r = F2x_div(r, v);
1730 205693 : t = v;
1731 : }
1732 146853 : if (F2x_degree(r) == 0) break;
1733 : }
1734 211796 : f = F2x_sqrt(r);
1735 : }
1736 1199540 : for (i = n; i; i--)
1737 1196898 : if (F2x_degree(gel(u,i))) break;
1738 336926 : setlg(u,i+1); return u;
1739 : }
1740 :
1741 : static GEN
1742 352711 : F2x_ddf_simple(GEN T, GEN XP)
1743 : {
1744 352711 : pari_sp av = avma, av2;
1745 : GEN f, z, Tr, X;
1746 352711 : long j, n = F2x_degree(T), v = T[1], B = n/2;
1747 352709 : if (n == 0) return cgetg(1, t_VEC);
1748 352710 : if (n == 1) return mkvec(T);
1749 137872 : z = XP; Tr = T; X = polx_F2x(v);
1750 137872 : f = const_vec(n, pol1_F2x(v));
1751 137877 : av2 = avma;
1752 358057 : for (j = 1; j <= B; j++)
1753 : {
1754 231380 : GEN u = F2x_gcd(Tr, F2x_add(z, X));
1755 231410 : if (F2x_degree(u))
1756 : {
1757 50729 : gel(f, j) = u;
1758 50729 : Tr = F2x_div(Tr, u);
1759 50725 : av2 = avma;
1760 180718 : } else z = gerepileuptoleaf(av2, z);
1761 231445 : if (!F2x_degree(Tr)) break;
1762 220226 : z = F2xq_sqr(z, Tr);
1763 : }
1764 137862 : if (F2x_degree(Tr)) gel(f, F2x_degree(Tr)) = Tr;
1765 137876 : return gerepilecopy(av, f);
1766 : }
1767 :
1768 : GEN
1769 7 : F2x_ddf(GEN T)
1770 : {
1771 : GEN XP;
1772 7 : T = F2x_get_red(T);
1773 7 : XP = F2x_Frobenius(T);
1774 7 : return F2x_ddf_to_ddf2(F2x_ddf_simple(T, XP));
1775 : }
1776 :
1777 : static GEN
1778 20368 : F2xq_frobtrace(GEN a, long d, GEN T)
1779 : {
1780 20368 : pari_sp av = avma;
1781 : long i;
1782 20368 : GEN x = a;
1783 59871 : for(i=1; i<d; i++)
1784 : {
1785 39504 : x = F2x_add(a, F2xq_sqr(x,T));
1786 39503 : if (gc_needed(av, 2))
1787 0 : x = gerepileuptoleaf(av, x);
1788 : }
1789 20367 : return x;
1790 : }
1791 :
1792 : static void
1793 29749 : F2x_edf_simple(GEN Tp, GEN XP, long d, GEN V, long idx)
1794 : {
1795 29749 : long n = F2x_degree(Tp), r = n/d;
1796 : GEN T, f, ff;
1797 29749 : if (r==1) { gel(V, idx) = Tp; return; }
1798 10012 : T = Tp;
1799 10012 : XP = F2x_rem(XP, T);
1800 : while (1)
1801 10356 : {
1802 20368 : pari_sp btop = avma;
1803 : long df;
1804 20368 : GEN g = random_F2x(n, Tp[1]);
1805 20368 : GEN t = F2xq_frobtrace(g, d, T);
1806 20367 : if (lgpol(t) == 0) continue;
1807 15262 : f = F2x_gcd(t, Tp); df = F2x_degree(f);
1808 15263 : if (df > 0 && df < n) break;
1809 5251 : set_avma(btop);
1810 : }
1811 10012 : ff = F2x_div(Tp, f);
1812 10012 : F2x_edf_simple(f, XP, d, V, idx);
1813 10012 : F2x_edf_simple(ff, XP, d, V, idx+F2x_degree(f)/d);
1814 : }
1815 :
1816 : static GEN
1817 352707 : F2x_factor_Shoup(GEN T)
1818 : {
1819 352707 : long i, n, s = 0;
1820 : GEN XP, D, V;
1821 : pari_timer ti;
1822 352707 : n = F2x_degree(T);
1823 352703 : if (DEBUGLEVEL>=6) timer_start(&ti);
1824 352703 : XP = F2x_Frobenius(T);
1825 352704 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1826 352704 : D = F2x_ddf_simple(T, XP);
1827 352710 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1828 1099099 : for (i = 1; i <= n; i++)
1829 746396 : s += F2x_degree(gel(D,i))/i;
1830 352703 : V = cgetg(s+1, t_COL);
1831 1099064 : for (i = 1, s = 1; i <= n; i++)
1832 : {
1833 746345 : GEN Di = gel(D,i);
1834 746345 : long ni = F2x_degree(Di), ri = ni/i;
1835 746342 : if (ni == 0) continue;
1836 392223 : if (ni == i) { gel(V, s++) = Di; continue; }
1837 9712 : F2x_edf_simple(Di, XP, i, V, s);
1838 9725 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_edf(%ld)",i);
1839 9725 : s += ri;
1840 : }
1841 352719 : return V;
1842 : }
1843 :
1844 : static GEN
1845 301628 : F2x_factor_Cantor(GEN T)
1846 : {
1847 301628 : GEN E, F, V = F2x_factor_squarefree(T);
1848 301622 : long i, j, l = lg(V);
1849 301622 : E = cgetg(l, t_VEC);
1850 301625 : F = cgetg(l, t_VEC);
1851 1087375 : for (i=1, j=1; i < l; i++)
1852 785745 : if (F2x_degree(gel(V,i)))
1853 : {
1854 352707 : GEN Fj = F2x_factor_Shoup(gel(V,i));
1855 352707 : gel(F, j) = Fj;
1856 352707 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1857 352706 : j++;
1858 : }
1859 301630 : return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
1860 : }
1861 :
1862 : #if 0
1863 : static GEN
1864 : F2x_simplefact_Shoup(GEN T)
1865 : {
1866 : long i, n, s = 0, j = 1, k;
1867 : GEN XP, D, V;
1868 : pari_timer ti;
1869 : n = F2x_degree(T);
1870 : if (DEBUGLEVEL>=6) timer_start(&ti);
1871 : XP = F2x_Frobenius(T);
1872 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1873 : D = F2x_ddf_simple(T, XP);
1874 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1875 : for (i = 1; i <= n; i++)
1876 : s += F2x_degree(gel(D,i))/i;
1877 : V = cgetg(s+1, t_VECSMALL);
1878 : for (i = 1; i <= n; i++)
1879 : {
1880 : long ni = F2x_degree(gel(D,i)), ri = ni/i;
1881 : if (ni == 0) continue;
1882 : for (k = 1; k <= ri; k++)
1883 : V[j++] = i;
1884 : }
1885 : return V;
1886 : }
1887 : static GEN
1888 : F2x_simplefact_Cantor(GEN T)
1889 : {
1890 : GEN E, F, V = F2x_factor_squarefree(T);
1891 : long i, j, l = lg(V);
1892 : F = cgetg(l, t_VEC);
1893 : E = cgetg(l, t_VEC);
1894 : for (i=1, j=1; i < l; i++)
1895 : if (F2x_degree(gel(V,i)))
1896 : {
1897 : GEN Fj = F2x_simplefact_Shoup(gel(V,i));
1898 : gel(F, j) = Fj;
1899 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1900 : j++;
1901 : }
1902 : return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
1903 : }
1904 : static int
1905 : F2x_isirred_Cantor(GEN T)
1906 : {
1907 : pari_sp av = avma;
1908 : pari_timer ti;
1909 : long n;
1910 : GEN dT = F2x_deriv(T);
1911 : GEN XP, D;
1912 : if (F2x_degree(F2x_gcd(T, dT)) != 0) return gc_bool(av,0);
1913 : n = F2x_degree(T);
1914 : if (DEBUGLEVEL>=6) timer_start(&ti);
1915 : XP = F2x_Frobenius(T);
1916 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1917 : D = F2x_ddf_simple(T, XP);
1918 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1919 : return gc_bool(av, F2x_degree(gel(D,n)) == n);
1920 : }
1921 : #endif
1922 :
1923 : /* driver for Cantor factorization, assume deg f > 2; not competitive for
1924 : * flag != 0, or as deg f increases */
1925 : static GEN
1926 301628 : F2x_Cantor_i(GEN f, long flag)
1927 : {
1928 : switch(flag)
1929 : {
1930 301628 : default: return F2x_factor_Cantor(f);
1931 : #if 0
1932 : case 1: return F2x_simplefact_Cantor(f);
1933 : case 2: return F2x_isirred_Cantor(f)? gen_1: NULL;
1934 : #endif
1935 : }
1936 : }
1937 : static GEN
1938 1524805 : F2x_factor_i(GEN f, long flag)
1939 : {
1940 1524805 : long d = F2x_degree(f);
1941 1524801 : if (d <= 2) return F2x_factor_deg2(f,d,flag);
1942 312419 : return (flag == 0 && d <= 20)? F2x_Cantor_i(f, flag)
1943 649421 : : F2x_Berlekamp_i(f, flag);
1944 : }
1945 :
1946 : GEN
1947 0 : F2x_degfact(GEN f)
1948 : {
1949 0 : pari_sp av = avma;
1950 0 : GEN z = F2x_factor_i(f, 1);
1951 0 : return gerepilecopy(av, z);
1952 : }
1953 :
1954 : int
1955 238 : F2x_is_irred(GEN f) { return !!F2x_factor_i(f, 2); }
1956 :
1957 : /* Adapted from Shoup NTL */
1958 : GEN
1959 7244851 : Flx_factor_squarefree_pre(GEN f, ulong p, ulong pi)
1960 : {
1961 7244851 : long i, q, n = degpol(f);
1962 7244614 : GEN u = const_vec(n+1, pol1_Flx(f[1]));
1963 7245356 : for(q = 1;;q *= p)
1964 160937 : {
1965 7406293 : GEN t, v, tv, r = Flx_gcd_pre(f, Flx_deriv(f, p), p, pi);
1966 7402293 : if (degpol(r) == 0) { gel(u, q) = f; break; }
1967 659029 : t = Flx_div_pre(f, r, p, pi);
1968 659028 : if (degpol(t) > 0)
1969 : {
1970 : long j;
1971 504812 : for(j = 1;;j++)
1972 : {
1973 1252903 : v = Flx_gcd_pre(r, t, p, pi);
1974 1252896 : tv = Flx_div_pre(t, v, p, pi);
1975 1252899 : if (degpol(tv) > 0)
1976 750873 : gel(u, j*q) = Flx_normalize(tv, p);
1977 1252895 : if (degpol(v) <= 0) break;
1978 748082 : r = Flx_div_pre(r, v, p, pi);
1979 748091 : t = v;
1980 : }
1981 504812 : if (degpol(r) == 0) break;
1982 : }
1983 160937 : f = Flx_normalize(Flx_deflate(r, p), p);
1984 : }
1985 30438224 : for (i = n; i; i--)
1986 30436067 : if (degpol(gel(u,i))) break;
1987 7242824 : setlg(u,i+1); return u;
1988 : }
1989 : GEN
1990 346557 : Flx_factor_squarefree(GEN f, ulong p)
1991 346557 : { return Flx_factor_squarefree_pre(f, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1992 :
1993 : long
1994 3444 : Flx_ispower(GEN f, ulong k, ulong p, GEN *pt_r)
1995 : {
1996 3444 : pari_sp av = avma;
1997 : ulong lc, pi;
1998 : GEN F;
1999 3444 : long i, n = degpol(f), v = f[1], l;
2000 3444 : if (n % k) return 0;
2001 3444 : lc = Fl_sqrtn(Flx_lead(f), k, p, NULL);
2002 3444 : if (lc == ULONG_MAX) { av = avma; return 0; }
2003 3444 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2004 3444 : F = Flx_factor_squarefree_pre(f, p, pi); l = lg(F)-1;
2005 38290 : for (i = 1; i <= l; i++)
2006 34867 : if (i%k && degpol(gel(F,i))) return gc_long(av,0);
2007 3423 : if (pt_r)
2008 : {
2009 3423 : GEN r = Fl_to_Flx(lc, v), s = pol1_Flx(v);
2010 38269 : for(i = l; i >= 1; i--)
2011 : {
2012 34846 : if (i%k) continue;
2013 12439 : s = Flx_mul_pre(s, gel(F,i), p, pi);
2014 12439 : r = Flx_mul_pre(r, s, p, pi);
2015 : }
2016 3423 : *pt_r = gerepileuptoleaf(av, r);
2017 0 : } else set_avma(av);
2018 3423 : return 1;
2019 : }
2020 :
2021 : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
2022 : static GEN
2023 9159590 : Flx_ddf_Shoup(GEN T, GEN XP, ulong p, ulong pi)
2024 : {
2025 9159590 : pari_sp av = avma;
2026 : GEN b, g, h, F, f, Tr, xq;
2027 : long i, j, n, v, bo, ro;
2028 : long B, l, m;
2029 : pari_timer ti;
2030 9159590 : n = get_Flx_degree(T); v = get_Flx_var(T);
2031 9160967 : if (n == 0) return cgetg(1, t_VEC);
2032 9117358 : if (n == 1) return mkvec(get_Flx_mod(T));
2033 8732920 : B = n/2;
2034 8732920 : l = usqrt(B);
2035 8733306 : m = (B+l-1)/l;
2036 8733306 : T = Flx_get_red(T, p);
2037 8732965 : b = cgetg(l+2, t_VEC);
2038 8732960 : gel(b, 1) = polx_Flx(v);
2039 8733571 : gel(b, 2) = XP;
2040 8733571 : bo = brent_kung_optpow(n, l-1, 1);
2041 8734150 : ro = l<=1 ? 0:(bo-1)/(l-1) + ((n-1)/bo);
2042 8734150 : if (DEBUGLEVEL>=7) timer_start(&ti);
2043 8734150 : if (expu(p) <= ro)
2044 875097 : for (i = 3; i <= l+1; i++)
2045 485069 : gel(b, i) = Flxq_powu_pre(gel(b, i-1), p, T, p, pi);
2046 : else
2047 : {
2048 8344004 : xq = Flxq_powers_pre(gel(b, 2), bo, T, p, pi);
2049 8342636 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq baby");
2050 9063427 : for (i = 3; i <= l+1; i++)
2051 720580 : gel(b, i) = Flx_FlxqV_eval_pre(gel(b, i-1), xq, T, p, pi);
2052 : }
2053 8732875 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: baby");
2054 8732875 : xq = Flxq_powers_pre(gel(b, l+1), brent_kung_optpow(n, m-1, 1), T, p, pi);
2055 8732755 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq giant");
2056 8732755 : g = cgetg(m+1, t_VEC);
2057 8733858 : gel(g, 1) = gel(xq, 2);
2058 14134345 : for(i = 2; i <= m; i++)
2059 5400285 : gel(g, i) = Flx_FlxqV_eval_pre(gel(g, i-1), xq, T, p, pi);
2060 8734060 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: giant");
2061 8734060 : h = cgetg(m+1, t_VEC);
2062 22867188 : for (j = 1; j <= m; j++)
2063 : {
2064 14134039 : pari_sp av = avma;
2065 14134039 : GEN gj = gel(g, j);
2066 14134039 : GEN e = Flx_sub(gj, gel(b, 1), p);
2067 18003413 : for (i = 2; i <= l; i++)
2068 3872217 : e = Flxq_mul_pre(e, Flx_sub(gj, gel(b, i), p), T, p, pi);
2069 14131196 : gel(h, j) = gerepileupto(av, e);
2070 : }
2071 8733149 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: diff");
2072 8733149 : Tr = get_Flx_mod(T);
2073 8733890 : F = cgetg(m+1, t_VEC);
2074 22866704 : for (j = 1; j <= m; j++)
2075 : {
2076 14133032 : GEN u = Flx_gcd_pre(Tr, gel(h, j), p, pi);
2077 14131534 : if (degpol(u))
2078 : {
2079 6535604 : u = Flx_normalize(u, p);
2080 6536322 : Tr = Flx_div_pre(Tr, u, p, pi);
2081 : }
2082 14131695 : gel(F, j) = u;
2083 : }
2084 8733672 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: F");
2085 8733672 : f = const_vec(n, pol1_Flx(v));
2086 22868064 : for (j = 1; j <= m; j++)
2087 : {
2088 14134571 : GEN e = gel(F, j);
2089 15012551 : for (i=l-1; i >= 0; i--)
2090 : {
2091 15012668 : GEN u = Flx_gcd_pre(e, Flx_sub(gel(g, j), gel(b, i+1), p), p, pi);
2092 15010187 : if (degpol(u))
2093 : {
2094 6636014 : gel(f, l*j-i) = u;
2095 6636014 : e = Flx_div_pre(e, u, p, pi);
2096 : }
2097 15009654 : if (!degpol(e)) break;
2098 : }
2099 : }
2100 8733493 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: f");
2101 8733493 : if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
2102 8733522 : return gerepilecopy(av, f);
2103 : }
2104 :
2105 : static void
2106 408990 : Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, ulong pi, GEN V, long idx)
2107 : {
2108 408990 : long n = degpol(Tp), r = n/d;
2109 : GEN T, f, ff;
2110 : ulong p2;
2111 408990 : if (r==1) { gel(V, idx) = Tp; return; }
2112 174676 : p2 = p>>1;
2113 174676 : T = Flx_get_red_pre(Tp, p, pi);
2114 174675 : XP = Flx_rem_pre(XP, T, p, pi);
2115 : while (1)
2116 20823 : {
2117 195498 : pari_sp btop = avma;
2118 : long i;
2119 195498 : GEN g = random_Flx(n, Tp[1], p);
2120 195499 : GEN t = gel(Flxq_auttrace_pre(mkvec2(XP, g), d, T, p, pi), 2);
2121 195499 : if (lgpol(t) == 0) continue;
2122 435252 : for(i=1; i<=10; i++)
2123 : {
2124 419689 : pari_sp btop2 = avma;
2125 419689 : GEN R = Flxq_powu_pre(Flx_Fl_add(t, random_Fl(p), p), p2, T, p, pi);
2126 419687 : f = Flx_gcd_pre(Flx_Fl_add(R, p-1, p), Tp, p, pi);
2127 419689 : if (degpol(f) > 0 && degpol(f) < n) break;
2128 245012 : set_avma(btop2);
2129 : }
2130 190239 : if (degpol(f) > 0 && degpol(f) < n) break;
2131 15563 : set_avma(btop);
2132 : }
2133 174676 : f = Flx_normalize(f, p);
2134 174676 : ff = Flx_div_pre(Tp, f, p, pi);
2135 174676 : Flx_edf_simple(f, XP, d, p, pi, V, idx);
2136 174676 : Flx_edf_simple(ff, XP, d, p, pi, V, idx+degpol(f)/d);
2137 : }
2138 : static void
2139 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, ulong pi, GEN V, long idx);
2140 :
2141 : static void
2142 1355686 : Flx_edf_rec(GEN T, GEN XP, GEN hp, GEN t, long d, ulong p, ulong pi,
2143 : GEN V, long idx)
2144 : {
2145 : pari_sp av;
2146 1355686 : GEN Tp = get_Flx_mod(T);
2147 1355691 : long n = degpol(hp), vT = Tp[1];
2148 : GEN u1, u2, f1, f2;
2149 1355692 : ulong p2 = p>>1;
2150 : GEN R, h;
2151 1355692 : h = Flx_get_red_pre(hp, p, pi);
2152 1355687 : t = Flx_rem_pre(t, T, p, pi);
2153 1355613 : av = avma;
2154 : do
2155 : {
2156 2280421 : set_avma(av);
2157 2280428 : R = Flxq_powu_pre(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p, pi);
2158 2280089 : u1 = Flx_gcd_pre(Flx_Fl_add(R, p-1, p), hp, p, pi);
2159 2280426 : } while (degpol(u1)==0 || degpol(u1)==n);
2160 1355627 : f1 = Flx_gcd_pre(Flx_Flxq_eval_pre(u1, t, T, p, pi), Tp, p, pi);
2161 1355566 : f1 = Flx_normalize(f1, p);
2162 1355633 : u2 = Flx_div_pre(hp, u1, p, pi);
2163 1355672 : f2 = Flx_div_pre(Tp, f1, p, pi);
2164 1355666 : if (degpol(u1)==1)
2165 : {
2166 1040348 : if (degpol(f1)==d)
2167 1026071 : gel(V, idx) = f1;
2168 : else
2169 14269 : Flx_edf(f1, XP, d, p, pi, V, idx);
2170 : }
2171 : else
2172 315370 : Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, pi, V, idx);
2173 1355704 : idx += degpol(f1)/d;
2174 1355674 : if (degpol(u2)==1)
2175 : {
2176 1033082 : if (degpol(f2)==d)
2177 1018947 : gel(V, idx) = f2;
2178 : else
2179 14135 : Flx_edf(f2, XP, d, p, pi, V, idx);
2180 : }
2181 : else
2182 322615 : Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, pi, V, idx);
2183 1355686 : }
2184 :
2185 : static void
2186 717771 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, ulong pi, GEN V, long idx)
2187 : {
2188 717771 : long n = degpol(Tp), r = n/d, vT = Tp[1];
2189 : GEN T, h, t;
2190 : pari_timer ti;
2191 717770 : if (r==1) { gel(V, idx) = Tp; return; }
2192 717770 : T = Flx_get_red_pre(Tp, p, pi);
2193 717770 : XP = Flx_rem_pre(XP, T, p, pi);
2194 717761 : if (DEBUGLEVEL>=7) timer_start(&ti);
2195 : do
2196 : {
2197 739851 : GEN g = random_Flx(n, vT, p);
2198 739854 : t = gel(Flxq_auttrace_pre(mkvec2(XP, g), d, T, p, pi), 2);
2199 739863 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
2200 739863 : h = Flxq_minpoly_pre(t, T, p, pi);
2201 739848 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
2202 739848 : } while (degpol(h) <= 1);
2203 717758 : Flx_edf_rec(T, XP, h, t, d, p, pi, V, idx);
2204 : }
2205 :
2206 : static GEN
2207 1653770 : Flx_factor_Shoup(GEN T, ulong p, ulong pi)
2208 : {
2209 1653770 : long i, n, s = 0, e = expu(p);
2210 : GEN XP, D, V;
2211 : pari_timer ti;
2212 1653770 : n = get_Flx_degree(T);
2213 1653767 : T = Flx_get_red_pre(T, p, pi);
2214 1653759 : if (DEBUGLEVEL>=6) timer_start(&ti);
2215 1653759 : XP = Flx_Frobenius_pre(T, p, pi);
2216 1653761 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2217 1653761 : D = Flx_ddf_Shoup(T, XP, p, pi);
2218 1653817 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2219 1653817 : s = ddf_to_nbfact(D);
2220 1653805 : V = cgetg(s+1, t_COL);
2221 9049271 : for (i = 1, s = 1; i <= n; i++)
2222 : {
2223 7395465 : GEN Di = gel(D,i);
2224 7395465 : long ni = degpol(Di), ri = ni/i;
2225 7395472 : if (ni == 0) continue;
2226 2149953 : Di = Flx_normalize(Di, p);
2227 2149979 : if (ni == i) { gel(V, s++) = Di; continue; }
2228 748998 : if (ri <= e*expu(e))
2229 689366 : Flx_edf(Di, XP, i, p, pi, V, s);
2230 : else
2231 59637 : Flx_edf_simple(Di, XP, i, p, pi, V, s);
2232 748992 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
2233 748992 : s += ri;
2234 : }
2235 1653806 : return V;
2236 : }
2237 :
2238 : static GEN
2239 1582129 : Flx_factor_Cantor(GEN T, ulong p)
2240 : {
2241 1582129 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2242 1582129 : GEN E, F, V = Flx_factor_squarefree_pre(get_Flx_mod(T), p, pi);
2243 1582111 : long i, j, l = lg(V);
2244 1582111 : F = cgetg(l, t_VEC);
2245 1582115 : E = cgetg(l, t_VEC);
2246 3786520 : for (i=1, j=1; i < l; i++)
2247 2204369 : if (degpol(gel(V,i)))
2248 : {
2249 1653770 : GEN Fj = Flx_factor_Shoup(gel(V,i), p, pi);
2250 1653806 : gel(F, j) = Fj;
2251 1653806 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
2252 1653805 : j++;
2253 : }
2254 1582151 : return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
2255 : }
2256 :
2257 : GEN
2258 0 : Flx_ddf_pre(GEN T, ulong p, ulong pi)
2259 : {
2260 : GEN XP;
2261 0 : T = Flx_get_red_pre(T, p, pi);
2262 0 : XP = Flx_Frobenius_pre(T, p, pi);
2263 0 : return ddf_to_ddf2(Flx_ddf_Shoup(T, XP, p, pi));
2264 : }
2265 : GEN
2266 0 : Flx_ddf(GEN T, ulong p)
2267 0 : { return Flx_ddf_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2268 :
2269 : static GEN
2270 5313261 : Flx_simplefact_Cantor(GEN T, ulong p)
2271 : {
2272 5313261 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2273 : long i, l;
2274 : GEN V;
2275 5313261 : T = Flx_get_red_pre(T, p, pi);
2276 5313032 : V = Flx_factor_squarefree_pre(get_Flx_mod(T), p, pi); l = lg(V);
2277 10709151 : for (i=1; i < l; i++)
2278 5393935 : gel(V,i) = Flx_ddf_Shoup(gel(V,i), Flx_Frobenius_pre(gel(V,i), p,pi), p,pi);
2279 5315216 : return vddf_to_simplefact(V, get_Flx_degree(T));
2280 : }
2281 :
2282 : static int
2283 1078 : Flx_isirred_Cantor(GEN Tp, ulong p)
2284 : {
2285 1078 : pari_sp av = avma;
2286 : pari_timer ti;
2287 1078 : GEN T = get_Flx_mod(Tp), dT = Flx_deriv(T, p), XP, D;
2288 1078 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2289 : long n;
2290 1078 : if (degpol(Flx_gcd_pre(T, dT, p, pi)) != 0) return gc_bool(av,0);
2291 791 : n = get_Flx_degree(T);
2292 791 : T = Flx_get_red_pre(Tp, p, pi);
2293 791 : if (DEBUGLEVEL>=6) timer_start(&ti);
2294 791 : XP = Flx_Frobenius_pre(T, p, pi);
2295 791 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2296 791 : D = Flx_ddf_Shoup(T, XP, p, pi);
2297 791 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2298 791 : return gc_bool(av, degpol(gel(D,n)) == n);
2299 : }
2300 :
2301 : /* f monic */
2302 : static GEN
2303 11417541 : Flx_factor_i(GEN f, ulong pp, long flag)
2304 : {
2305 : long d;
2306 11417541 : if (pp==2) { /*We need to handle 2 specially */
2307 73899 : GEN F = F2x_factor_i(Flx_to_F2x(f),flag);
2308 73899 : if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
2309 73899 : return F;
2310 : }
2311 11343642 : d = degpol(f);
2312 11343195 : if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
2313 6896506 : switch(flag)
2314 : {
2315 1582130 : default: return Flx_factor_Cantor(f, pp);
2316 5313298 : case 1: return Flx_simplefact_Cantor(f, pp);
2317 1078 : case 2: return Flx_isirred_Cantor(f, pp)? gen_1: NULL;
2318 : }
2319 : }
2320 :
2321 : GEN
2322 7823750 : Flx_degfact(GEN f, ulong p)
2323 : {
2324 7823750 : pari_sp av = avma;
2325 7823750 : GEN z = Flx_factor_i(Flx_normalize(f,p),p,1);
2326 7822488 : return gerepilecopy(av, z);
2327 : }
2328 :
2329 : /* T must be squarefree mod p*/
2330 : GEN
2331 1452788 : Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
2332 : {
2333 : GEN XP, D;
2334 : pari_timer ti;
2335 1452788 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2336 1452788 : long i, s, n = get_Flx_degree(T);
2337 1452766 : GEN V = const_vecsmall(n, 0);
2338 1452766 : pari_sp av = avma;
2339 1452766 : T = Flx_get_red_pre(T, p, pi);
2340 1452789 : if (DEBUGLEVEL>=6) timer_start(&ti);
2341 1452789 : XP = Flx_Frobenius_pre(T, p, pi);
2342 1452699 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2343 1452699 : D = Flx_ddf_Shoup(T, XP, p, pi);
2344 1452953 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2345 7478760 : for (i = 1, s = 0; i <= n; i++) { V[i] = degpol(gel(D,i))/i; s += V[i]; }
2346 1452921 : *nb = s; set_avma(av); return V;
2347 : }
2348 :
2349 : long
2350 660268 : Flx_nbfact_Frobenius_pre(GEN T, GEN XP, ulong p, ulong pi)
2351 : {
2352 660268 : pari_sp av = avma;
2353 660268 : long s = ddf_to_nbfact(Flx_ddf_Shoup(T, XP, p, pi));
2354 660284 : return gc_long(av,s);
2355 : }
2356 : long
2357 0 : Flx_nbfact_Frobenius(GEN T, GEN XP, ulong p)
2358 0 : { return Flx_nbfact_Frobenius_pre(T, XP, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2359 :
2360 : /* T must be squarefree mod p*/
2361 : long
2362 660263 : Flx_nbfact_pre(GEN T, ulong p, ulong pi)
2363 : {
2364 660263 : pari_sp av = avma;
2365 660263 : GEN XP = Flx_Frobenius_pre(T, p, pi);
2366 660268 : long n = Flx_nbfact_Frobenius_pre(T, XP, p, pi);
2367 660284 : return gc_long(av,n);
2368 : }
2369 : long
2370 660262 : Flx_nbfact(GEN T, ulong p)
2371 660262 : { return Flx_nbfact_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2372 :
2373 : int
2374 1057 : Flx_is_irred(GEN f, ulong p)
2375 : {
2376 1057 : pari_sp av = avma;
2377 1057 : f = Flx_normalize(f,p);
2378 1057 : return gc_bool(av, !!Flx_factor_i(f,p,2));
2379 : }
2380 :
2381 : /* Use this function when you think f is reducible, and that there are lots of
2382 : * factors. If you believe f has few factors, use FpX_nbfact(f,p)==1 instead */
2383 : int
2384 112 : FpX_is_irred(GEN f, GEN p)
2385 : {
2386 112 : pari_sp av = avma;
2387 : int z;
2388 112 : switch(ZX_factmod_init(&f,p))
2389 : {
2390 28 : case 0: z = !!F2x_factor_i(f,2); break;
2391 77 : case 1: z = !!Flx_factor_i(f,p[2],2); break;
2392 7 : default: z = !!FpX_factor_i(f,p,2); break;
2393 : }
2394 112 : return gc_bool(av,z);
2395 : }
2396 : GEN
2397 47947 : FpX_degfact(GEN f, GEN p) {
2398 47947 : pari_sp av = avma;
2399 : GEN F;
2400 47947 : switch(ZX_factmod_init(&f,p))
2401 : {
2402 476 : case 0: F = F2x_factor_i(f,1); break;
2403 47443 : case 1: F = Flx_factor_i(f,p[2],1); break;
2404 28 : default: F = FpX_factor_i(f,p,1); break;
2405 : }
2406 47947 : return gerepilecopy(av, F);
2407 : }
2408 :
2409 : #if 0
2410 : /* set x <-- x + c*y mod p */
2411 : /* x is not required to be normalized.*/
2412 : static void
2413 : Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
2414 : {
2415 : long i, lx, ly;
2416 : ulong *x=(ulong *)gx;
2417 : ulong *y=(ulong *)gy;
2418 : if (!c) return;
2419 : lx = lg(gx);
2420 : ly = lg(gy);
2421 : if (lx<ly) pari_err_BUG("lx<ly in Flx_addmul_inplace");
2422 : if (SMALL_ULONG(p))
2423 : for (i=2; i<ly; i++) x[i] = (x[i] + c*y[i]) % p;
2424 : else
2425 : for (i=2; i<ly; i++) x[i] = Fl_add(x[i], Fl_mul(c,y[i],p),p);
2426 : }
2427 : #endif
2428 :
2429 : GEN
2430 4302999 : FpX_factor(GEN f, GEN p)
2431 : {
2432 4302999 : pari_sp av = avma;
2433 : GEN F;
2434 4302999 : switch(ZX_factmod_init(&f, p))
2435 : {
2436 1435103 : case 0: F = F2x_factor_i(f,0);
2437 1435107 : F2xV_to_ZXV_inplace(gel(F,1)); break;
2438 2865615 : case 1: F = Flx_factor_i(f,p[2],0);
2439 2865634 : FlxV_to_ZXV_inplace(gel(F,1)); break;
2440 2348 : default: F = FpX_factor_i(f,p,0); break;
2441 : }
2442 4303119 : return gerepilecopy(av, F);
2443 : }
2444 :
2445 : GEN
2446 680207 : Flx_factor(GEN f, ulong p)
2447 : {
2448 680207 : pari_sp av = avma;
2449 680207 : return gerepilecopy(av, Flx_factor_i(Flx_normalize(f,p),p,0));
2450 : }
2451 : GEN
2452 15062 : F2x_factor(GEN f)
2453 : {
2454 15062 : pari_sp av = avma;
2455 15062 : return gerepilecopy(av, F2x_factor_i(f,0));
2456 : }
|
