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 4343434 : ZX_factmod_init(GEN *F, GEN p)
35 : {
36 4343434 : if (lgefint(p) == 3)
37 : {
38 4340927 : ulong pp = p[2];
39 4340927 : if (pp == 2) { *F = ZX_to_F2x(*F); return 0; }
40 2910635 : *F = ZX_to_Flx(*F, pp);
41 2910689 : if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
42 2910687 : return 1;
43 : }
44 2507 : *F = FpX_red(*F, p);
45 2508 : if (lg(*F) > 3) *F = FpX_normalize(*F, p);
46 2508 : return 2;
47 : }
48 : static GEN
49 659518 : ZX_rootmod_init(GEN F, GEN p)
50 659518 : { 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 2186 : Flx_Xnm1(long sv, long n, ulong p)
71 : {
72 2186 : GEN t = cgetg(n+3, t_VECSMALL);
73 : long i;
74 2186 : t[1] = sv;
75 2186 : t[2] = p - 1;
76 6927 : for (i = 3; i <= n+1; i++) t[i] = 0;
77 2186 : t[i] = 1; return t;
78 : }
79 : /* X^n + 1 */
80 : static GEN
81 2076 : Flx_Xn1(long sv, long n, ulong p)
82 : {
83 2076 : GEN t = cgetg(n+3, t_VECSMALL);
84 : long i;
85 : (void) p;
86 2076 : t[1] = sv;
87 2076 : t[2] = 1;
88 6729 : for (i = 3; i <= n+1; i++) t[i] = 0;
89 2076 : t[i] = 1; return t;
90 : }
91 :
92 : /* assume lg(f) > 3 */
93 : static GEN
94 92715 : Flx_root_mod_2(GEN f)
95 : {
96 92715 : long i, n = lg(f)-1, c = f[2];
97 92715 : int z0 = !c;
98 92715 : c ^= 1; /* c = f[2] + f[n] mod 2, we know f[n] is odd */
99 166795 : for (i=3; i < n; i++) c ^= f[i];
100 : /* c = 0 iff f(1) = 0 (mod 2) */
101 92715 : 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 16299319 : 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 425634 : FpX_roots(GEN f, GEN p)
123 : {
124 425634 : pari_sp av = avma;
125 425634 : GEN y; f = ZX_rootmod_init(f, p);
126 425634 : switch(lg(f))
127 : {
128 14 : case 2: pari_err_ROOTS0("FpX_roots");
129 47438 : case 3: return cgetg(1,t_COL);
130 : }
131 378182 : if (typ(f) == t_VECSMALL)
132 : {
133 369942 : ulong pp = p[2];
134 369942 : if (pp == 2)
135 92708 : y = Flx_root_mod_2(f);
136 : else
137 : {
138 277234 : if (!odd(pp)) pari_err_PRIME("FpX_roots", p);
139 277235 : y = Flx_roots_pre(f, pp, SMALL_ULONG(pp)? 0: get_Fl_red(pp));
140 : }
141 369936 : y = Flc_to_ZC(y);
142 : }
143 : else
144 8240 : y = FpX_roots_i(f, p);
145 378168 : 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 14457 : FpX_quad_root(GEN x, GEN p, int unknown)
159 : {
160 14457 : GEN s, D, b = gel(x,3), c = gel(x,2);
161 :
162 14457 : if (absequaliu(p, 2)) {
163 0 : if (!signe(b)) return c;
164 0 : return signe(c)? NULL: gen_1;
165 : }
166 14457 : D = subii(sqri(b), shifti(c,2));
167 14457 : D = remii(D,p);
168 14457 : if (unknown && kronecker(D,p) == -1) return NULL;
169 :
170 13867 : s = Fp_sqrt(D,p);
171 : /* p is not prime, go on and give e.g. maxord a chance to recover */
172 13867 : if (!s) return NULL;
173 13859 : return Fp_halve(Fp_sub(s,b, p), p);
174 : }
175 : static GEN
176 9994 : FpX_otherroot(GEN x, GEN r, GEN p)
177 9994 : { 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 27443341 : Fl_disc_bc(ulong b, ulong c, ulong p)
182 27443341 : { 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 24957013 : Flx_quad_root(GEN x, ulong p, ulong pi, int unknown)
186 : {
187 24957013 : ulong s, b = x[3], c = x[2];
188 24957013 : ulong D = Fl_disc_bc(b, c, p);
189 24929133 : if (unknown && krouu(D,p) == -1) return p;
190 16473032 : s = Fl_sqrt_pre(D, p, pi);
191 16568896 : if (s==~0UL) return p;
192 16568883 : return Fl_halve(Fl_sub(s,b, p), p);
193 : }
194 : static ulong
195 14755316 : Flx_otherroot(GEN x, ulong r, ulong p)
196 14755316 : { 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 5076198 : split_init(struct split_t *S, long max)
203 : {
204 5076198 : S->todo = vectrunc_init(max);
205 5075745 : S->done = vectrunc_init(max);
206 5075243 : }
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 5442861 : 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 5443937 : split_moveto_done(struct split_t *S, long i, GEN t)
224 : {
225 5443937 : long n = lg(S->todo)-1;
226 5443937 : vectrunc_append(S->done, t);
227 5444208 : if (n) gel(S->todo,i) = gel(S->todo, n);
228 5444208 : setlg(S->todo, n);
229 :
230 5444115 : }
231 : /* append t to done */
232 : static void
233 506301 : split_add_done(struct split_t *S, GEN t)
234 506301 : { vectrunc_append(S->done, t); }
235 : /* split todo[i] into a and b */
236 : static void
237 416671 : split_todo(struct split_t *S, long i, GEN a, GEN b)
238 : {
239 416671 : gel(S->todo, i) = a;
240 416671 : split_add(S, b);
241 416669 : }
242 : /* split todo[i] into a and b, moved to done */
243 : static void
244 466435 : split_done(struct split_t *S, long i, GEN a, GEN b)
245 : {
246 466435 : split_moveto_done(S, i, a);
247 466434 : split_add_done(S, b);
248 466433 : }
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 416347 : Flx_cubic_root(GEN ff, ulong p, ulong pi)
322 : {
323 416347 : GEN f = Flx_normalize(ff,p);
324 416342 : 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 416342 : if (pi)
327 : {
328 416344 : t = Fl_mul_pre(a, p3, p, pi);
329 416353 : t2 = Fl_sqr_pre(t, p, pi);
330 416346 : A = Fl_sub(b, Fl_triple(t2, p), p);
331 416344 : B = Fl_sub(c, Fl_mul_pre(t, Fl_add(A, t2, p), p, pi), p);
332 416350 : A3 = Fl_mul_pre(A, p3, p, pi);
333 416351 : 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 416351 : A33 = Fl_powu_pre(A3, 3, p, pi);
345 416348 : D = Fl_add(B2, Fl_double(Fl_double(A33, p), p), p);
346 416347 : S = Fl_neg(B,p);
347 416346 : P = Fl_neg(A3,p);
348 416346 : if (krouu(D,p) >= 0)
349 : {
350 341663 : ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
351 341670 : ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
352 341669 : if (p%3==2) /* 1 solutions */
353 138580 : vS1 = Fl_powu_pre(S1, p - p3, p, pi);
354 : else
355 : {
356 203089 : vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
357 203085 : if (vS1==~0UL) return p; /*0 solutions*/
358 : /*3 solutions*/
359 : }
360 222609 : if (!P) return Fl_sub(vS1, t, p);
361 94768 : vS2 = pi? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): Fl_div(P, vS1, p);
362 94769 : return Fl_sub(Fl_add(vS1,vS2, p), t, p);
363 : }
364 : else
365 : {
366 74685 : pari_sp av = avma;
367 74685 : GEN S1 = mkvecsmall2(Fl_halve(S, p), (p + 1UL) >> 1);
368 74686 : GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
369 : ulong Sa;
370 74687 : if (!vS1) return p; /*0 solutions, p%3==2*/
371 74687 : Sa = vS1[1];
372 74687 : if (p%3==1) /*1 solutions*/
373 : {
374 29233 : ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
375 29234 : if (Fa!=P) Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
376 : }
377 74688 : set_avma(av);
378 74688 : 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 4152864 : 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 4152864 : if (Flx_val(f)) return 0;
441 4151313 : da = degpol(f); f = Flx_normalize(f, p);
442 4150939 : if (da == 1) return Fl_neg(f[2], p);
443 4137345 : PI = pi? pi: get_Fl_red(p); /* PI for Fp, pi for Fp[x] */
444 4138206 : switch(da)
445 : {
446 3397166 : case 2: return Flx_quad_root(f, p, PI, 1);
447 402347 : case 3: if (p>3) return Flx_cubic_root(f, p, PI); /*FALL THROUGH*/
448 : }
449 345371 : if (SMALL_ULONG(p)) pi = 0; /* bilinear ops faster without Fl_*_pre */
450 345371 : if (!fl)
451 : {
452 308063 : a = Flxq_powu_pre(polx_Flx(f[1]), p - 1, f,p,pi);
453 307988 : if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
454 307988 : a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
455 307973 : a = Flx_gcd_pre(f,a, p, pi);
456 37308 : } else a = f;
457 345346 : da = degpol(a);
458 345349 : if (!da) return p;
459 246541 : a = Flx_normalize(a,p);
460 :
461 246572 : q = p >> 1;
462 246572 : pol = polx_Flx(f[1]);
463 363055 : for(pol[2] = 1;; pol[2]++)
464 : {
465 363055 : if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
466 363065 : switch(da)
467 : {
468 155184 : case 1: return Fl_neg(a[2], p);
469 70727 : case 2: return Flx_quad_root(a, p, PI, 0);
470 20680 : case 3: if (p>3) return Flx_cubic_root(a, p, PI); /*FALL THROUGH*/
471 : default: {
472 116474 : GEN b = Flxq_powu_pre(pol,q, a,p,pi);
473 : long db;
474 116473 : if (degpol(b) <= 0) continue;
475 110705 : b = Flx_gcd_pre(a,Flx_Fl_add(b,p-1,p), p, pi);
476 110705 : db = degpol(b); if (!db) continue;
477 110706 : b = Flx_normalize(b, p);
478 110711 : if (db <= (da >> 1)) {
479 68177 : a = b;
480 68177 : da = db;
481 : } else {
482 42534 : a = Flx_div_pre(a,b, p, pi);
483 42536 : da -= db;
484 : }
485 : }
486 : }
487 : }
488 : }
489 : ulong
490 4114298 : Flx_oneroot_pre(GEN f, ulong p, ulong pi)
491 4114298 : { return Flx_oneroot_pre_i(f, p, pi, 0); }
492 : ulong
493 38777 : Flx_oneroot_split_pre(GEN f, ulong p, ulong pi)
494 38777 : { 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 821 : case 1: return subii(p, gel(f,2));
508 3837 : 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 2485 : Flx_oneroot(GEN f, ulong p)
551 : {
552 2485 : pari_sp av = avma;
553 2485 : switch(lg(f))
554 : {
555 0 : case 2: return 0;
556 0 : case 3: return p;
557 : }
558 2485 : if (p == 2) return Flx_oneroot_mod_2(f);
559 2485 : 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 233884 : FpX_oneroot(GEN f, GEN p)
578 : {
579 233884 : pari_sp av = avma;
580 233884 : f = ZX_rootmod_init(f, p);
581 233884 : switch(lg(f))
582 : {
583 0 : case 2: set_avma(av); return gen_0;
584 0 : case 3: return gc_NULL(av);
585 : }
586 233884 : if (typ(f) == t_VECSMALL)
587 : {
588 228715 : ulong r, pp = p[2];
589 228715 : if (pp == 2)
590 91 : r = Flx_oneroot_mod_2(f);
591 : else
592 228624 : r = Flx_oneroot_pre(f, pp, SMALL_ULONG(pp)? 0: get_Fl_red(pp));
593 228715 : set_avma(av);
594 228715 : 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 301 : for (i = 0; i < n; i++)
651 : {
652 231 : GEN tmp = FpX_gcd(FpX_Fp_sub(xpow, zeta, p), fact, p);
653 231 : long dtmp = degpol(tmp);
654 231 : if (dtmp > 0 && dtmp < deg_f)
655 : {
656 154 : fact = FpX_div(fact, tmp, p); deg_f = degpol(fact);
657 154 : 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 147 : 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 1882510 : FE_concat(GEN F, GEN E, long l)
682 : {
683 1882510 : setlg(E,l); E = shallowconcat1(E);
684 1882529 : 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 5318000 : vddf_to_simplefact(GEN V, long d)
716 : {
717 : GEN E, F;
718 5318000 : long i, j, c, l = lg(V);
719 5318000 : F = cgetg(d+1, t_VECSMALL);
720 5317009 : E = cgetg(d+1, t_VECSMALL);
721 10716266 : for (i = c = 1; i < l; i++)
722 : {
723 5399813 : GEN Vi = gel(V,i);
724 5399813 : long l = lg(Vi);
725 27224199 : for (j = 1; j < l; j++)
726 : {
727 21824461 : long k, n = degpol(gel(Vi,j)) / j;
728 33048145 : for (k = 1; k <= n; k++) { uel(F,c) = j; uel(E,c) = i; c++; }
729 : }
730 : }
731 5316453 : setlg(F,c);
732 5316383 : setlg(E,c);
733 5316366 : 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 234191 : FpX_split_part(GEN f, GEN p)
739 : {
740 234191 : long n = degpol(f);
741 234191 : GEN z, X = pol_x(varn(f));
742 234191 : if (n <= 1) return f;
743 232223 : f = FpX_red(f, p);
744 232220 : z = FpX_sub(FpX_Frobenius(f, p), X, p);
745 232224 : 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 138146 : FpX_nbroots(GEN f, GEN p)
752 : {
753 138146 : pari_sp av = avma;
754 138146 : GEN z = FpX_split_part(f, p);
755 138146 : return gc_long(av, degpol(z));
756 : }
757 :
758 : /* 1 < deg(f) <= p */
759 : static int
760 81347 : Flx_is_totally_split_i(GEN f, ulong p)
761 : {
762 81347 : GEN F = Flx_Frobenius(f, p);
763 81347 : return degpol(F)==1 && uel(F,2)==0UL && uel(F,3)==1UL;
764 : }
765 : int
766 81354 : Flx_is_totally_split(GEN f, ulong p)
767 : {
768 81354 : pari_sp av = avma;
769 81354 : ulong n = degpol(f);
770 81354 : if (n <= 1) return 1;
771 81347 : if (n > p) return 0; /* includes n < 0 */
772 81347 : 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 4388719 : Flx_nbroots(GEN f, ulong p)
794 : {
795 4388719 : long n = degpol(f);
796 : ulong pi;
797 4388718 : pari_sp av = avma;
798 : GEN z;
799 4388718 : if (n <= 1) return n;
800 4374991 : if (n == 2)
801 : {
802 : ulong D;
803 578544 : if (p==2) return (f[2]==0) + (f[2]!=f[3]);
804 433787 : D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
805 433791 : return 1 + krouu(D,p);
806 : }
807 3796447 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
808 3796447 : z = Flx_sub(Flx_Frobenius_pre(f, p, pi), polx_Flx(f[1]), p);
809 3796442 : z = Flx_gcd_pre(z, f, p, pi);
810 3796443 : 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 1021 : 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 1021 : n = get_FpX_degree(T); v = get_FpX_var(T);
863 1021 : if (n == 0) return cgetg(1, t_VEC);
864 1021 : if (n == 1) return mkvec(get_FpX_mod(T));
865 848 : B = n/2;
866 848 : l = usqrt(B);
867 848 : m = (B+l-1)/l;
868 848 : T = FpX_get_red(T, p);
869 848 : b = cgetg(l+2, t_VEC);
870 848 : gel(b, 1) = pol_x(v);
871 848 : gel(b, 2) = XP;
872 848 : if (DEBUGLEVEL>=7) timer_start(&ti);
873 848 : xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1), T, p);
874 848 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq baby");
875 1077 : for (i = 3; i <= l+1; i++)
876 229 : gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
877 848 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: baby");
878 848 : xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1), T, p);
879 848 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq giant");
880 848 : g = cgetg(m+1, t_VEC);
881 848 : 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 848 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: giant");
884 848 : h = cgetg(m+1, t_VEC);
885 2691 : 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 848 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: diff");
893 848 : Tr = get_FpX_mod(T);
894 848 : F = cgetg(m+1, t_VEC);
895 2691 : 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 481 : u = FpX_normalize(u, p);
901 481 : Tr = FpX_div(Tr, u, p);
902 : }
903 1843 : gel(F,j) = u;
904 : }
905 848 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: F");
906 848 : f = const_vec(n, pol_1(v));
907 2691 : 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 511 : u = FpX_normalize(u, p);
916 511 : gel(f, l*j-i) = u;
917 511 : e = FpX_div(e, u, p);
918 : }
919 1915 : if (!degpol(e)) break;
920 : }
921 : }
922 848 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: f");
923 848 : if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
924 848 : 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 1660 : set_avma(av);
974 1660 : R = FpXQ_pow(deg1pol(gen_1, randomi(p), vT), p2, h, p);
975 1660 : u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
976 1660 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_rec",p);
977 1660 : } 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 784 : gel(V, idx) = f1;
984 : else
985 291 : 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 728 : gel(V, idx) = f2;
989 : else
990 347 : 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 437 : FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
996 : {
997 437 : long n = degpol(Tp), r = n/d, vT = varn(Tp), ct = 0;
998 : GEN T, h, t;
999 : pari_timer ti;
1000 :
1001 437 : T = FpX_get_red(Tp, p);
1002 437 : XP = FpX_rem(XP, T, p);
1003 437 : if (DEBUGLEVEL>=7) timer_start(&ti);
1004 : do
1005 : {
1006 437 : GEN g = random_FpX(n, vT, p);
1007 437 : t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
1008 437 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
1009 437 : h = FpXQ_minpoly(t, T, p);
1010 437 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
1011 437 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf",p);
1012 437 : } while (degpol(h) != r);
1013 437 : FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
1014 437 : }
1015 :
1016 : static GEN
1017 1007 : FpX_factor_Shoup(GEN T, GEN p)
1018 : {
1019 1007 : long i, n, s = 0;
1020 : GEN XP, D, V;
1021 1007 : long e = expi(p);
1022 : pari_timer ti;
1023 1007 : n = get_FpX_degree(T);
1024 1007 : T = FpX_get_red(T, p);
1025 1007 : if (DEBUGLEVEL>=6) timer_start(&ti);
1026 1007 : XP = FpX_Frobenius(T, p);
1027 1007 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
1028 1007 : D = FpX_ddf_Shoup(T, XP, p);
1029 1007 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
1030 1007 : s = ddf_to_nbfact(D);
1031 1007 : V = cgetg(s+1, t_COL);
1032 7292 : for (i = 1, s = 1; i <= n; i++)
1033 : {
1034 6285 : GEN Di = gel(D,i);
1035 6285 : long ni = degpol(Di), ri = ni/i;
1036 6285 : if (ni == 0) continue;
1037 1125 : Di = FpX_normalize(Di, p);
1038 1125 : if (ni == i) { gel(V, s++) = Di; continue; }
1039 437 : if (ri <= e*expu(e))
1040 437 : FpX_edf(Di, XP, i, p, V, s);
1041 : else
1042 0 : FpX_edf_simple(Di, XP, i, p, V, s);
1043 437 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
1044 437 : s += ri;
1045 : }
1046 1007 : return V;
1047 : }
1048 :
1049 : long
1050 2315096 : ddf_to_nbfact(GEN D)
1051 : {
1052 2315096 : long l = lg(D), i, s = 0;
1053 14646347 : for(i = 1; i < l; i++) s += degpol(gel(D,i))/i;
1054 2315091 : return s;
1055 : }
1056 :
1057 : /* Yun algorithm: Assume p > degpol(T) */
1058 : static GEN
1059 1696 : FpX_factor_Yun(GEN T, GEN p)
1060 : {
1061 1696 : long n = degpol(T), i = 1;
1062 1696 : GEN a, b, c, d = FpX_deriv(T, p);
1063 1696 : GEN V = cgetg(n+1,t_VEC);
1064 1696 : a = FpX_gcd(T, d, p);
1065 1696 : if (degpol(a) == 0) return mkvec(T);
1066 584 : b = FpX_div(T, a, p);
1067 : do
1068 : {
1069 2638 : c = FpX_div(d, a, p);
1070 2638 : d = FpX_sub(c, FpX_deriv(b, p), p);
1071 2638 : a = FpX_normalize(FpX_gcd(b, d, p), p);
1072 2638 : gel(V, i++) = a;
1073 2638 : b = FpX_div(b, a, p);
1074 2638 : } while (degpol(b));
1075 584 : setlg(V, i); return V;
1076 : }
1077 : GEN
1078 341852 : FpX_factor_squarefree(GEN T, GEN p)
1079 : {
1080 341852 : if (lgefint(p)==3)
1081 : {
1082 341118 : ulong pp = (ulong)p[2];
1083 341118 : GEN u = Flx_factor_squarefree(ZX_to_Flx(T,pp), pp);
1084 341118 : return FlxV_to_ZXV(u);
1085 : }
1086 734 : return FpX_factor_Yun(T, p);
1087 : }
1088 :
1089 : GEN
1090 221697 : FpX_roots_mult(GEN T, long n, GEN p)
1091 : {
1092 221697 : pari_sp av = avma;
1093 221697 : GEN V = FpX_factor_squarefree(T, p), W;
1094 221697 : long l = lg(V), i;
1095 221697 : if (l <= n) { set_avma(av); return cgetg(1,t_COL); }
1096 36776 : W = cgetg(l-n+1,t_VEC);
1097 120325 : for (i = n; i < l; i++)
1098 83549 : gel(W,i-n+1) = FpX_roots(gel(V,i), p);
1099 36776 : 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 906 : FpX_factor_Cantor(GEN T, GEN p)
1138 : {
1139 906 : GEN E, F, V = FpX_factor_Yun(T, p);
1140 906 : long i, j, l = lg(V);
1141 906 : F = cgetg(l, t_VEC);
1142 906 : E = cgetg(l, t_VEC);
1143 2426 : for (i=1, j=1; i < l; i++)
1144 1520 : if (degpol(gel(V,i)))
1145 : {
1146 1007 : GEN Fj = FpX_factor_Shoup(gel(V,i), p);
1147 1007 : gel(F, j) = Fj;
1148 1007 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1149 1007 : j++;
1150 : }
1151 906 : 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 2508 : FpX_factor_i(GEN f, GEN pp, long flag)
1222 : {
1223 2508 : long d = degpol(f);
1224 2508 : if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
1225 920 : switch(flag)
1226 : {
1227 906 : 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 349935 : trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
1283 :
1284 : /* not gerepile safe */
1285 : static GEN
1286 1567 : FpX_factor_2(GEN f, GEN p, long d)
1287 : {
1288 : GEN r, s, R, S;
1289 : long v;
1290 : int sgn;
1291 1567 : 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 674 : case 1: retmkvec2(mkcol(f), mkvecsmall(1));
1296 : }
1297 849 : r = FpX_quad_root(f, p, 1);
1298 849 : if (!r) return mkvec2(mkcol(f), mkvecsmall(1));
1299 258 : v = varn(f);
1300 258 : s = FpX_otherroot(f, r, p);
1301 258 : if (signe(r)) r = subii(p, r);
1302 258 : if (signe(s)) s = subii(p, s);
1303 258 : sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
1304 258 : R = deg1pol_shallow(gen_1, r, v);
1305 258 : if (!sgn) return mkvec2(mkcol(R), mkvecsmall(2));
1306 156 : S = deg1pol_shallow(gen_1, s, v);
1307 156 : return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
1308 : }
1309 : static GEN
1310 1588 : FpX_factor_deg2(GEN f, GEN p, long d, long flag)
1311 : {
1312 1588 : 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 1567 : default: return FpX_factor_2(f, p, d);
1316 : }
1317 : }
1318 :
1319 : static int
1320 451466 : F2x_quad_factortype(GEN x)
1321 451466 : { 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 19677 : F2x_degfact_2(GEN f, long d)
1329 : {
1330 19677 : if (!d) return trivial_fact();
1331 19677 : if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1332 19481 : switch(F2x_quad_factortype(f)) {
1333 5964 : case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1334 6244 : case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1335 7273 : default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1336 : }
1337 : }
1338 :
1339 : static GEN
1340 1164909 : F2x_factor_2(GEN f, long d)
1341 : {
1342 1164909 : long v = f[1];
1343 1164909 : if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1344 909934 : if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1345 427176 : switch(F2x_quad_factortype(f))
1346 : {
1347 71702 : case -1: return mkvec2(mkcol(f), mkvecsmall(1));
1348 327982 : case 0: return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
1349 27515 : default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
1350 : }
1351 : }
1352 : static GEN
1353 1184600 : F2x_factor_deg2(GEN f, long d, long flag)
1354 : {
1355 1184600 : switch(flag) {
1356 14 : case 2: return F2x_is_irred_2(f, d);
1357 19677 : case 1: return F2x_degfact_2(f, d);
1358 1164909 : default: return F2x_factor_2(f, d);
1359 : }
1360 : }
1361 :
1362 : /* xt = NULL or x^(p-1)/2 mod g */
1363 : static void
1364 30282 : split_squares(struct split_t *S, GEN g, ulong p, ulong pi, GEN xt)
1365 : {
1366 30282 : ulong q = p >> 1;
1367 30282 : GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
1368 30282 : long d = degpol(a);
1369 30282 : if (d < 0)
1370 : {
1371 : ulong i;
1372 407 : split_add_done(S, (GEN)1);
1373 407 : if (!pi)
1374 1555 : 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 29875 : if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
1379 29875 : if (d)
1380 : {
1381 29786 : if (xt) xt = Flx_Fl_add(xt, p-1, p); else xt = Flx_Xnm1(g[1], q, p);
1382 29786 : a = Flx_gcd_pre(a, xt, p, pi);
1383 29786 : if (degpol(a)) split_add(S, Flx_normalize(a, p));
1384 : }
1385 : }
1386 30282 : }
1387 : static void
1388 30282 : split_nonsquares(struct split_t *S, GEN g, ulong p, ulong pi, GEN xt)
1389 : {
1390 30282 : ulong q = p >> 1;
1391 30282 : GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
1392 30282 : long d = degpol(a);
1393 30282 : if (d < 0)
1394 : {
1395 249 : ulong i, z = nonsquare_Fl(p);
1396 249 : split_add_done(S, (GEN)z);
1397 249 : if (!pi)
1398 522 : for (i = 2; i <= q; i++)
1399 273 : 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 30033 : if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
1405 30033 : if (d)
1406 : {
1407 29676 : if (xt) xt = Flx_Fl_add(xt, 1, p); else xt = Flx_Xn1(g[1], q, p);
1408 29676 : a = Flx_gcd_pre(a, xt, p, pi);
1409 29676 : if (degpol(a)) split_add(S, Flx_normalize(a, p));
1410 : }
1411 : }
1412 30282 : }
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 5067125 : split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p, ulong pi)
1417 : {
1418 5067125 : GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
1419 5067099 : long d = degpol(g);
1420 5067012 : if (d < 0) return 0;
1421 5066412 : if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
1422 5066412 : if (!d) return 1;
1423 5036440 : if ((p >> 4) <= (ulong)d)
1424 : { /* small p; split directly using x^((p-1)/2) +/- 1 */
1425 27600 : GEN xt = ((ulong)d < (p>>1))? Flx_rem_pre(monomial_Flx(1, p>>1, g[1]), g, p, pi)
1426 30282 : : NULL;
1427 30282 : split_squares(S, g, p, pi, xt);
1428 30282 : split_nonsquares(S, g, p, pi, xt);
1429 : } else { /* large p; use x^(p-1) - 1 directly */
1430 5006158 : a = Flxq_powu_pre(polx_Flx(f[1]), p-1, g, p, pi);
1431 4994178 : if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
1432 4994178 : a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
1433 4995580 : g = Flx_gcd_pre(g,a, p,pi);
1434 4999904 : if (degpol(g)) split_add(S, Flx_normalize(g,p));
1435 : }
1436 5035150 : return 1;
1437 : }
1438 :
1439 : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
1440 : GEN
1441 24679759 : Flx_roots_pre(GEN f, ulong p, ulong pi)
1442 : {
1443 : GEN pol;
1444 24679759 : long v = Flx_valrem(f, &f), n = degpol(f);
1445 : ulong q, PI;
1446 : struct split_t S;
1447 :
1448 24676256 : f = Flx_normalize(f, p);
1449 : /* optimization: test for small degree first */
1450 24690984 : if (n == 1)
1451 : {
1452 114919 : q = p - f[2];
1453 114919 : return v? mkvecsmall2(0, q): mkvecsmall(q);
1454 : }
1455 24576065 : PI = pi? pi: get_Fl_red(p); /* PI for Fp, pi for Fp[x] */
1456 24577467 : if (n == 2)
1457 : {
1458 19509090 : ulong r = Flx_quad_root(f, p, PI, 1), s;
1459 19527640 : if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
1460 13316203 : s = Flx_otherroot(f,r, p);
1461 13364262 : if (r < s)
1462 3357632 : return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
1463 10006630 : else if (r > s)
1464 10039200 : return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
1465 : else
1466 4832 : return v? mkvecsmall2(0, s): mkvecsmall(s);
1467 : }
1468 5068377 : if (SMALL_ULONG(p)) pi = 0; /* bilinear ops faster without Fl_*_pre */
1469 5068377 : q = p >> 1;
1470 5068377 : split_init(&S, lg(f)-1);
1471 5067120 : settyp(S.done, t_VECSMALL);
1472 5067120 : if (v) split_add_done(&S, (GEN)0);
1473 5067120 : if (! split_Flx_cut_out_roots(&S, f, p, pi))
1474 600 : return all_roots_mod_p(p, lg(S.done) == 1);
1475 5064981 : pol = polx_Flx(f[1]);
1476 5066466 : for (pol[2]=1; ; pol[2]++)
1477 5467475 : {
1478 10533941 : long j, l = lg(S.todo);
1479 10533941 : if (l == 1) { vecsmall_sort(S.done); return S.done; }
1480 5467655 : if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
1481 11425244 : for (j = 1; j < l; j++)
1482 : {
1483 5957769 : GEN b, c = gel(S.todo,j);
1484 : ulong r, s;
1485 5957769 : switch(degpol(c))
1486 : {
1487 4973088 : case 1:
1488 4973088 : split_moveto_done(&S, j, (GEN)(p - c[2]));
1489 4973347 : j--; l--; break;
1490 460203 : case 2:
1491 460203 : r = Flx_quad_root(c, p, PI, 0);
1492 460236 : if (r == p) pari_err_PRIME("polrootsmod",utoipos(p));
1493 460229 : s = Flx_otherroot(c,r, p);
1494 460224 : split_done(&S, j, (GEN)r, (GEN)s);
1495 460221 : j--; l--; break;
1496 524062 : default:
1497 524062 : b = Flxq_powu_pre(pol,q, c,p,pi); /* pol^(p-1)/2 */
1498 524247 : if (degpol(b) <= 0) continue;
1499 410868 : b = Flx_gcd_pre(c,Flx_Fl_add(b,p-1,p), p, pi);
1500 410932 : if (!degpol(b)) continue;
1501 410415 : b = Flx_normalize(b, p);
1502 410443 : c = Flx_div_pre(c,b, p,pi);
1503 410434 : 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 2490066 : Flx_quad_factortype(GEN x, ulong p)
1527 : {
1528 2490066 : ulong b = x[3], c = x[2];
1529 2490066 : 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 2516113 : Flx_degfact_2(GEN f, ulong p, long d)
1540 : {
1541 2516113 : if (!d) return trivial_fact();
1542 2516113 : if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1543 2490073 : switch(Flx_quad_factortype(f, p)) {
1544 1180453 : case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1545 1279147 : case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1546 30358 : default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1547 : }
1548 : }
1549 : /* p > 2 */
1550 : static GEN
1551 1932488 : Flx_factor_2(GEN f, ulong p, long d)
1552 : {
1553 : ulong r, s;
1554 : GEN R,S;
1555 1932488 : long v = f[1];
1556 1932488 : if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1557 1855886 : if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1558 1536432 : r = Flx_quad_root(f, p, get_Fl_red(p), 1);
1559 1536483 : if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
1560 923205 : s = Flx_otherroot(f, r, p);
1561 923202 : r = Fl_neg(r, p);
1562 923201 : s = Fl_neg(s, p);
1563 923198 : if (s < r) lswap(s,r);
1564 923198 : R = mkvecsmall3(v,r,1);
1565 923198 : if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
1566 812381 : S = mkvecsmall3(v,s,1);
1567 812388 : return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
1568 : }
1569 : static GEN
1570 4448538 : Flx_factor_deg2(GEN f, ulong p, long d, long flag)
1571 : {
1572 4448538 : switch(flag) {
1573 56 : case 2: return Flx_is_irred_2(f, p, d);
1574 2516132 : case 1: return Flx_degfact_2(f, p, d);
1575 1932350 : default: return Flx_factor_2(f, p, d);
1576 : }
1577 : }
1578 :
1579 : static GEN
1580 23299 : F2x_Berlekamp_ker(GEN u)
1581 : {
1582 23299 : pari_sp ltop=avma;
1583 23299 : long j,N = F2x_degree(u);
1584 : GEN Q;
1585 : pari_timer T;
1586 23298 : timer_start(&T);
1587 23299 : Q = F2x_matFrobenius(u);
1588 317949 : for (j=1; j<=N; j++)
1589 294650 : F2m_flip(Q,j,j);
1590 23299 : if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
1591 23299 : Q = F2m_ker_sp(Q,0);
1592 23299 : if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
1593 23299 : return gerepileupto(ltop,Q);
1594 : }
1595 : #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
1596 : static long
1597 35790 : F2x_split_Berlekamp(GEN *t)
1598 : {
1599 35790 : GEN u = *t, a, b, vker;
1600 35790 : long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
1601 :
1602 35789 : if (du == 1) return 1;
1603 27345 : if (du == 2)
1604 : {
1605 4046 : 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 4046 : return 1;
1612 : }
1613 :
1614 23299 : vker = F2x_Berlekamp_ker(u);
1615 23299 : lb = lgcols(vker);
1616 23312 : d = lg(vker)-1;
1617 23312 : ir = 0;
1618 : /* t[i] irreducible for i < ir, still to be treated for i < L */
1619 60317 : for (L=1; L<d; )
1620 : {
1621 : GEN pol;
1622 37019 : if (d == 2)
1623 5548 : pol = F2v_to_F2x(gel(vker,2), sv);
1624 : else
1625 : {
1626 31471 : GEN v = zero_zv(lb);
1627 31481 : v[1] = du;
1628 31481 : v[2] = random_Fl(2); /*Assume vker[1]=1*/
1629 118525 : for (i=2; i<=d; i++)
1630 87044 : if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
1631 31481 : pol = F2v_to_F2x(v, sv);
1632 : }
1633 111498 : for (i=ir; i<L && L<d; i++)
1634 : {
1635 74493 : a = t[i]; la = F2x_degree(a);
1636 74480 : if (la == 1) { set_irred(i); }
1637 74268 : else if (la == 2)
1638 : {
1639 733 : 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 733 : set_irred(i);
1645 : }
1646 : else
1647 : {
1648 73535 : pari_sp av = avma;
1649 : long lb;
1650 73535 : b = F2x_rem(pol, a);
1651 73520 : if (F2x_degree(b) <= 0) { set_avma(av); continue; }
1652 26554 : b = F2x_gcd(a,b); lb = F2x_degree(b);
1653 26554 : if (lb && lb < la)
1654 : {
1655 26554 : t[L] = F2x_div(a,b);
1656 26555 : t[i]= b; L++;
1657 : }
1658 0 : else set_avma(av);
1659 : }
1660 : }
1661 : }
1662 23298 : return d;
1663 : }
1664 : /* assume deg f > 2 */
1665 : static GEN
1666 35366 : F2x_Berlekamp_i(GEN f, long flag)
1667 : {
1668 35366 : long lfact, val, d = F2x_degree(f), j, k, lV;
1669 : GEN y, E, t, V;
1670 :
1671 35366 : val = F2x_valrem(f, &f);
1672 35365 : if (flag == 2 && val) return NULL;
1673 35351 : V = F2x_factor_squarefree(f); lV = lg(V);
1674 35352 : if (flag == 2 && lV > 2) return NULL;
1675 :
1676 : /* to hold factors and exponents */
1677 35282 : t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
1678 35282 : E = cgetg(d+1,t_VECSMALL);
1679 35282 : lfact = 1;
1680 35282 : if (val) {
1681 11268 : if (flag == 1) t[1] = 1; else gel(t,1) = polx_F2x(f[1]);
1682 11268 : E[1] = val; lfact++;
1683 : }
1684 :
1685 117158 : for (k=1; k<lV; k++)
1686 : {
1687 82034 : if (F2x_degree(gel(V, k))==0) continue;
1688 35790 : gel(t,lfact) = gel(V, k);
1689 35790 : d = F2x_split_Berlekamp(&gel(t,lfact));
1690 35786 : if (flag == 2 && d != 1) return NULL;
1691 35632 : if (flag == 1)
1692 47783 : for (j=0; j<d; j++) t[lfact+j] = F2x_degree(gel(t,lfact+j));
1693 97404 : for (j=0; j<d; j++) E[lfact+j] = k;
1694 35632 : lfact += d;
1695 : }
1696 35124 : if (flag == 2) return gen_1; /* irreducible */
1697 35110 : setlg(t, lfact);
1698 35110 : setlg(E, lfact); y = mkvec2(t,E);
1699 24343 : return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
1700 59455 : : sort_factor_pol(y, cmpGuGu);
1701 : }
1702 :
1703 : /* Adapted from Shoup NTL */
1704 : GEN
1705 335024 : F2x_factor_squarefree(GEN f)
1706 : {
1707 : GEN r, t, v, tv;
1708 335024 : long i, q, n = F2x_degree(f);
1709 335023 : GEN u = const_vec(n+1, pol1_F2x(f[1]));
1710 335028 : for(q = 1;;q *= 2)
1711 : {
1712 543324 : r = F2x_gcd(f, F2x_deriv(f));
1713 543306 : if (F2x_degree(r) == 0)
1714 : {
1715 240933 : gel(u, q) = f;
1716 240933 : break;
1717 : }
1718 302381 : t = F2x_div(f, r);
1719 302382 : if (F2x_degree(t) > 0)
1720 : {
1721 : long j;
1722 146467 : for(j = 1;;j++)
1723 : {
1724 351437 : v = F2x_gcd(r, t);
1725 351439 : tv = F2x_div(t, v);
1726 351436 : if (F2x_degree(tv) > 0)
1727 148439 : gel(u, j*q) = tv;
1728 351436 : if (F2x_degree(v) <= 0) break;
1729 204968 : r = F2x_div(r, v);
1730 204970 : t = v;
1731 : }
1732 146466 : if (F2x_degree(r) == 0) break;
1733 : }
1734 208292 : f = F2x_sqrt(r);
1735 : }
1736 1193385 : for (i = n; i; i--)
1737 1190721 : if (F2x_degree(gel(u,i))) break;
1738 334990 : setlg(u,i+1); return u;
1739 : }
1740 :
1741 : static GEN
1742 350753 : F2x_ddf_simple(GEN T, GEN XP)
1743 : {
1744 350753 : pari_sp av = avma, av2;
1745 : GEN f, z, Tr, X;
1746 350753 : long j, n = F2x_degree(T), v = T[1], B = n/2;
1747 350753 : if (n == 0) return cgetg(1, t_VEC);
1748 350753 : if (n == 1) return mkvec(T);
1749 136912 : z = XP; Tr = T; X = polx_F2x(v);
1750 136917 : f = const_vec(n, pol1_F2x(v));
1751 136918 : av2 = avma;
1752 356080 : for (j = 1; j <= B; j++)
1753 : {
1754 230395 : GEN u = F2x_gcd(Tr, F2x_add(z, X));
1755 230398 : if (F2x_degree(u))
1756 : {
1757 50706 : gel(f, j) = u;
1758 50706 : Tr = F2x_div(Tr, u);
1759 50707 : av2 = avma;
1760 179741 : } else z = gerepileuptoleaf(av2, z);
1761 230452 : if (!F2x_degree(Tr)) break;
1762 219192 : z = F2xq_sqr(z, Tr);
1763 : }
1764 136903 : if (F2x_degree(Tr)) gel(f, F2x_degree(Tr)) = Tr;
1765 136914 : 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 20061 : F2xq_frobtrace(GEN a, long d, GEN T)
1779 : {
1780 20061 : pari_sp av = avma;
1781 : long i;
1782 20061 : GEN x = a;
1783 59056 : for(i=1; i<d; i++)
1784 : {
1785 38998 : x = F2x_add(a, F2xq_sqr(x,T));
1786 38995 : if (gc_needed(av, 2))
1787 0 : x = gerepileuptoleaf(av, x);
1788 : }
1789 20058 : return x;
1790 : }
1791 :
1792 : static void
1793 29845 : F2x_edf_simple(GEN Tp, GEN XP, long d, GEN V, long idx)
1794 : {
1795 29845 : long n = F2x_degree(Tp), r = n/d;
1796 : GEN T, f, ff;
1797 29845 : if (r==1) { gel(V, idx) = Tp; return; }
1798 10044 : T = Tp;
1799 10044 : XP = F2x_rem(XP, T);
1800 : while (1)
1801 10016 : {
1802 20060 : pari_sp btop = avma;
1803 : long df;
1804 20060 : GEN g = random_F2x(n, Tp[1]);
1805 20061 : GEN t = F2xq_frobtrace(g, d, T);
1806 20060 : if (lgpol(t) == 0) continue;
1807 15126 : f = F2x_gcd(t, Tp); df = F2x_degree(f);
1808 15125 : if (df > 0 && df < n) break;
1809 5082 : set_avma(btop);
1810 : }
1811 10043 : ff = F2x_div(Tp, f);
1812 10044 : F2x_edf_simple(f, XP, d, V, idx);
1813 10044 : F2x_edf_simple(ff, XP, d, V, idx+F2x_degree(f)/d);
1814 : }
1815 :
1816 : static GEN
1817 350750 : F2x_factor_Shoup(GEN T)
1818 : {
1819 350750 : long i, n, s = 0;
1820 : GEN XP, D, V;
1821 : pari_timer ti;
1822 350750 : n = F2x_degree(T);
1823 350747 : if (DEBUGLEVEL>=6) timer_start(&ti);
1824 350747 : XP = F2x_Frobenius(T);
1825 350746 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1826 350746 : D = F2x_ddf_simple(T, XP);
1827 350749 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1828 1093955 : for (i = 1; i <= n; i++)
1829 743206 : s += F2x_degree(gel(D,i))/i;
1830 350749 : V = cgetg(s+1, t_COL);
1831 1093950 : for (i = 1, s = 1; i <= n; i++)
1832 : {
1833 743183 : GEN Di = gel(D,i);
1834 743183 : long ni = F2x_degree(Di), ri = ni/i;
1835 743187 : if (ni == 0) continue;
1836 390218 : if (ni == i) { gel(V, s++) = Di; continue; }
1837 9744 : F2x_edf_simple(Di, XP, i, V, s);
1838 9757 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_edf(%ld)",i);
1839 9757 : s += ri;
1840 : }
1841 350767 : return V;
1842 : }
1843 :
1844 : static GEN
1845 299673 : F2x_factor_Cantor(GEN T)
1846 : {
1847 299673 : GEN E, F, V = F2x_factor_squarefree(T);
1848 299671 : long i, j, l = lg(V);
1849 299671 : E = cgetg(l, t_VEC);
1850 299675 : F = cgetg(l, t_VEC);
1851 1076016 : for (i=1, j=1; i < l; i++)
1852 776341 : if (F2x_degree(gel(V,i)))
1853 : {
1854 350750 : GEN Fj = F2x_factor_Shoup(gel(V,i));
1855 350749 : gel(F, j) = Fj;
1856 350749 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1857 350748 : j++;
1858 : }
1859 299675 : 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 299673 : F2x_Cantor_i(GEN f, long flag)
1927 : {
1928 : switch(flag)
1929 : {
1930 299673 : 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 1519643 : F2x_factor_i(GEN f, long flag)
1939 : {
1940 1519643 : long d = F2x_degree(f);
1941 1519633 : if (d <= 2) return F2x_factor_deg2(f,d,flag);
1942 310442 : return (flag == 0 && d <= 20)? F2x_Cantor_i(f, flag)
1943 645474 : : 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 7241677 : Flx_factor_squarefree_pre(GEN f, ulong p, ulong pi)
1960 : {
1961 7241677 : long i, q, n = degpol(f);
1962 7241486 : GEN u = const_vec(n+1, pol1_Flx(f[1]));
1963 7242849 : for(q = 1;;q *= p)
1964 160861 : {
1965 7403710 : GEN t, v, tv, r = Flx_gcd_pre(f, Flx_deriv(f, p), p, pi);
1966 7399756 : if (degpol(r) == 0) { gel(u, q) = f; break; }
1967 653375 : t = Flx_div_pre(f, r, p, pi);
1968 653379 : if (degpol(t) > 0)
1969 : {
1970 : long j;
1971 499242 : for(j = 1;;j++)
1972 : {
1973 1234516 : v = Flx_gcd_pre(r, t, p, pi);
1974 1234505 : tv = Flx_div_pre(t, v, p, pi);
1975 1234511 : if (degpol(tv) > 0)
1976 744650 : gel(u, j*q) = Flx_normalize(tv, p);
1977 1234513 : if (degpol(v) <= 0) break;
1978 735267 : r = Flx_div_pre(r, v, p, pi);
1979 735274 : t = v;
1980 : }
1981 499245 : if (degpol(r) == 0) break;
1982 : }
1983 160861 : f = Flx_normalize(Flx_deflate(r, p), p);
1984 : }
1985 30447099 : for (i = n; i; i--)
1986 30442092 : if (degpol(gel(u,i))) break;
1987 7243578 : setlg(u,i+1); return u;
1988 : }
1989 : GEN
1990 341118 : Flx_factor_squarefree(GEN f, ulong p)
1991 341118 : { return Flx_factor_squarefree_pre(f, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1992 :
1993 : long
1994 3451 : Flx_ispower(GEN f, ulong k, ulong p, GEN *pt_r)
1995 : {
1996 3451 : pari_sp av = avma;
1997 : ulong lc, pi;
1998 : GEN F;
1999 3451 : long i, n = degpol(f), v = f[1], l;
2000 3451 : if (n % k) return 0;
2001 3451 : lc = Fl_sqrtn(Flx_lead(f), k, p, NULL);
2002 3451 : if (lc == ULONG_MAX) { av = avma; return 0; }
2003 3451 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2004 3451 : F = Flx_factor_squarefree_pre(f, p, pi); l = lg(F)-1;
2005 38325 : for (i = 1; i <= l; i++)
2006 34895 : if (i%k && degpol(gel(F,i))) return gc_long(av,0);
2007 3430 : if (pt_r)
2008 : {
2009 3430 : GEN r = Fl_to_Flx(lc, v), s = pol1_Flx(v);
2010 38304 : for(i = l; i >= 1; i--)
2011 : {
2012 34874 : if (i%k) continue;
2013 12453 : s = Flx_mul_pre(s, gel(F,i), p, pi);
2014 12453 : r = Flx_mul_pre(r, s, p, pi);
2015 : }
2016 3430 : *pt_r = gerepileuptoleaf(av, r);
2017 0 : } else set_avma(av);
2018 3430 : return 1;
2019 : }
2020 :
2021 : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
2022 : static GEN
2023 9161279 : Flx_ddf_Shoup(GEN T, GEN XP, ulong p, ulong pi)
2024 : {
2025 9161279 : 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 9161279 : n = get_Flx_degree(T); v = get_Flx_var(T);
2031 9162662 : if (n == 0) return cgetg(1, t_VEC);
2032 9119002 : if (n == 1) return mkvec(get_Flx_mod(T));
2033 8734802 : B = n/2;
2034 8734802 : l = usqrt(B);
2035 8735280 : m = (B+l-1)/l;
2036 8735280 : T = Flx_get_red(T, p);
2037 8734935 : b = cgetg(l+2, t_VEC);
2038 8734857 : gel(b, 1) = polx_Flx(v);
2039 8735702 : gel(b, 2) = XP;
2040 8735702 : bo = brent_kung_optpow(n, l-1, 1);
2041 8736437 : ro = l<=1 ? 0:(bo-1)/(l-1) + ((n-1)/bo);
2042 8736437 : if (DEBUGLEVEL>=7) timer_start(&ti);
2043 8736437 : if (expu(p) <= ro)
2044 875032 : for (i = 3; i <= l+1; i++)
2045 485085 : gel(b, i) = Flxq_powu_pre(gel(b, i-1), p, T, p, pi);
2046 : else
2047 : {
2048 8346395 : xq = Flxq_powers_pre(gel(b, 2), bo, T, p, pi);
2049 8344344 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq baby");
2050 9065315 : for (i = 3; i <= l+1; i++)
2051 720733 : gel(b, i) = Flx_FlxqV_eval_pre(gel(b, i-1), xq, T, p, pi);
2052 : }
2053 8734529 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: baby");
2054 8734529 : xq = Flxq_powers_pre(gel(b, l+1), brent_kung_optpow(n, m-1, 1), T, p, pi);
2055 8734229 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq giant");
2056 8734229 : g = cgetg(m+1, t_VEC);
2057 8735322 : gel(g, 1) = gel(xq, 2);
2058 14137897 : for(i = 2; i <= m; i++)
2059 5402366 : gel(g, i) = Flx_FlxqV_eval_pre(gel(g, i-1), xq, T, p, pi);
2060 8735531 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: giant");
2061 8735531 : h = cgetg(m+1, t_VEC);
2062 22872860 : for (j = 1; j <= m; j++)
2063 : {
2064 14138202 : pari_sp av = avma;
2065 14138202 : GEN gj = gel(g, j);
2066 14138202 : GEN e = Flx_sub(gj, gel(b, 1), p);
2067 18007863 : for (i = 2; i <= l; i++)
2068 3872937 : e = Flxq_mul_pre(e, Flx_sub(gj, gel(b, i), p), T, p, pi);
2069 14134926 : gel(h, j) = gerepileupto(av, e);
2070 : }
2071 8734658 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: diff");
2072 8734658 : Tr = get_Flx_mod(T);
2073 8735606 : F = cgetg(m+1, t_VEC);
2074 22872976 : for (j = 1; j <= m; j++)
2075 : {
2076 14137256 : GEN u = Flx_gcd_pre(Tr, gel(h, j), p, pi);
2077 14135622 : if (degpol(u))
2078 : {
2079 6537913 : u = Flx_normalize(u, p);
2080 6538545 : Tr = Flx_div_pre(Tr, u, p, pi);
2081 : }
2082 14136087 : gel(F, j) = u;
2083 : }
2084 8735720 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: F");
2085 8735720 : f = const_vec(n, pol1_Flx(v));
2086 22874520 : for (j = 1; j <= m; j++)
2087 : {
2088 14139085 : GEN e = gel(F, j);
2089 15017099 : for (i=l-1; i >= 0; i--)
2090 : {
2091 15016952 : GEN u = Flx_gcd_pre(e, Flx_sub(gel(g, j), gel(b, i+1), p), p, pi);
2092 15013993 : if (degpol(u))
2093 : {
2094 6638083 : gel(f, l*j-i) = u;
2095 6638083 : e = Flx_div_pre(e, u, p, pi);
2096 : }
2097 15013518 : if (!degpol(e)) break;
2098 : }
2099 : }
2100 8735435 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: f");
2101 8735435 : if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
2102 8735540 : return gerepilecopy(av, f);
2103 : }
2104 :
2105 : static void
2106 408843 : Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, ulong pi, GEN V, long idx)
2107 : {
2108 408843 : long n = degpol(Tp), r = n/d;
2109 : GEN T, f, ff;
2110 : ulong p2;
2111 408843 : if (r==1) { gel(V, idx) = Tp; return; }
2112 174687 : p2 = p>>1;
2113 174687 : T = Flx_get_red_pre(Tp, p, pi);
2114 174688 : XP = Flx_rem_pre(XP, T, p, pi);
2115 : while (1)
2116 20665 : {
2117 195349 : pari_sp btop = avma;
2118 : long i;
2119 195349 : GEN g = random_Flx(n, Tp[1], p);
2120 195352 : GEN t = gel(Flxq_auttrace_pre(mkvec2(XP, g), d, T, p, pi), 2);
2121 195353 : if (lgpol(t) == 0) continue;
2122 434089 : for(i=1; i<=10; i++)
2123 : {
2124 418501 : pari_sp btop2 = avma;
2125 418501 : GEN R = Flxq_powu_pre(Flx_Fl_add(t, random_Fl(p), p), p2, T, p, pi);
2126 418500 : f = Flx_gcd_pre(Flx_Fl_add(R, p-1, p), Tp, p, pi);
2127 418497 : if (degpol(f) > 0 && degpol(f) < n) break;
2128 243812 : set_avma(btop2);
2129 : }
2130 190274 : if (degpol(f) > 0 && degpol(f) < n) break;
2131 15588 : set_avma(btop);
2132 : }
2133 174686 : f = Flx_normalize(f, p);
2134 174688 : ff = Flx_div_pre(Tp, f, p, pi);
2135 174688 : Flx_edf_simple(f, XP, d, p, pi, V, idx);
2136 174688 : 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 1357048 : 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 1357048 : GEN Tp = get_Flx_mod(T);
2147 1357048 : long n = degpol(hp), vT = Tp[1];
2148 : GEN u1, u2, f1, f2;
2149 1357049 : ulong p2 = p>>1;
2150 : GEN R, h;
2151 1357049 : h = Flx_get_red_pre(hp, p, pi);
2152 1357047 : t = Flx_rem_pre(t, T, p, pi);
2153 1356955 : av = avma;
2154 : do
2155 : {
2156 2284765 : set_avma(av);
2157 2284764 : R = Flxq_powu_pre(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p, pi);
2158 2284615 : u1 = Flx_gcd_pre(Flx_Fl_add(R, p-1, p), hp, p, pi);
2159 2284766 : } while (degpol(u1)==0 || degpol(u1)==n);
2160 1356968 : f1 = Flx_gcd_pre(Flx_Flxq_eval_pre(u1, t, T, p, pi), Tp, p, pi);
2161 1356985 : f1 = Flx_normalize(f1, p);
2162 1357030 : u2 = Flx_div_pre(hp, u1, p, pi);
2163 1357054 : f2 = Flx_div_pre(Tp, f1, p, pi);
2164 1357072 : if (degpol(u1)==1)
2165 : {
2166 1040117 : if (degpol(f1)==d)
2167 1025738 : gel(V, idx) = f1;
2168 : else
2169 14372 : Flx_edf(f1, XP, d, p, pi, V, idx);
2170 : }
2171 : else
2172 316988 : Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, pi, V, idx);
2173 1357096 : idx += degpol(f1)/d;
2174 1357080 : if (degpol(u2)==1)
2175 : {
2176 1034819 : if (degpol(f2)==d)
2177 1020834 : gel(V, idx) = f2;
2178 : else
2179 13984 : Flx_edf(f2, XP, d, p, pi, V, idx);
2180 : }
2181 : else
2182 322277 : Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, pi, V, idx);
2183 1357083 : }
2184 :
2185 : static void
2186 717852 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, ulong pi, GEN V, long idx)
2187 : {
2188 717852 : long n = degpol(Tp), r = n/d, vT = Tp[1];
2189 : GEN T, h, t;
2190 : pari_timer ti;
2191 717850 : if (r==1) { gel(V, idx) = Tp; return; }
2192 717850 : T = Flx_get_red_pre(Tp, p, pi);
2193 717849 : XP = Flx_rem_pre(XP, T, p, pi);
2194 717843 : if (DEBUGLEVEL>=7) timer_start(&ti);
2195 : do
2196 : {
2197 739391 : GEN g = random_Flx(n, vT, p);
2198 739403 : t = gel(Flxq_auttrace_pre(mkvec2(XP, g), d, T, p, pi), 2);
2199 739406 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
2200 739406 : h = Flxq_minpoly_pre(t, T, p, pi);
2201 739386 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
2202 739386 : } while (degpol(h) <= 1);
2203 717838 : Flx_edf_rec(T, XP, h, t, d, p, pi, V, idx);
2204 : }
2205 :
2206 : static GEN
2207 1653901 : Flx_factor_Shoup(GEN T, ulong p, ulong pi)
2208 : {
2209 1653901 : long i, n, s = 0, e = expu(p);
2210 : GEN XP, D, V;
2211 : pari_timer ti;
2212 1653901 : n = get_Flx_degree(T);
2213 1653898 : T = Flx_get_red_pre(T, p, pi);
2214 1653895 : if (DEBUGLEVEL>=6) timer_start(&ti);
2215 1653895 : XP = Flx_Frobenius_pre(T, p, pi);
2216 1653846 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2217 1653846 : D = Flx_ddf_Shoup(T, XP, p, pi);
2218 1653924 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2219 1653924 : s = ddf_to_nbfact(D);
2220 1653907 : V = cgetg(s+1, t_COL);
2221 9050674 : for (i = 1, s = 1; i <= n; i++)
2222 : {
2223 7396765 : GEN Di = gel(D,i);
2224 7396765 : long ni = degpol(Di), ri = ni/i;
2225 7396758 : if (ni == 0) continue;
2226 2149683 : Di = Flx_normalize(Di, p);
2227 2149718 : if (ni == i) { gel(V, s++) = Di; continue; }
2228 748950 : if (ri <= e*expu(e))
2229 689495 : Flx_edf(Di, XP, i, p, pi, V, s);
2230 : else
2231 59465 : Flx_edf_simple(Di, XP, i, p, pi, V, s);
2232 748959 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
2233 748959 : s += ri;
2234 : }
2235 1653909 : return V;
2236 : }
2237 :
2238 : static GEN
2239 1581948 : Flx_factor_Cantor(GEN T, ulong p)
2240 : {
2241 1581948 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2242 1581948 : GEN E, F, V = Flx_factor_squarefree_pre(get_Flx_mod(T), p, pi);
2243 1581940 : long i, j, l = lg(V);
2244 1581940 : F = cgetg(l, t_VEC);
2245 1581945 : E = cgetg(l, t_VEC);
2246 3784061 : for (i=1, j=1; i < l; i++)
2247 2202109 : if (degpol(gel(V,i)))
2248 : {
2249 1653900 : GEN Fj = Flx_factor_Shoup(gel(V,i), p, pi);
2250 1653908 : gel(F, j) = Fj;
2251 1653908 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
2252 1653901 : j++;
2253 : }
2254 1581952 : 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 5315894 : Flx_simplefact_Cantor(GEN T, ulong p)
2271 : {
2272 5315894 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2273 : long i, l;
2274 : GEN V;
2275 5315894 : T = Flx_get_red_pre(T, p, pi);
2276 5315520 : V = Flx_factor_squarefree_pre(get_Flx_mod(T), p, pi); l = lg(V);
2277 10715635 : for (i=1; i < l; i++)
2278 5397239 : gel(V,i) = Flx_ddf_Shoup(gel(V,i), Flx_Frobenius_pre(gel(V,i), p,pi), p,pi);
2279 5318396 : 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 11422069 : Flx_factor_i(GEN f, ulong pp, long flag)
2304 : {
2305 : long d;
2306 11422069 : if (pp==2) { /*We need to handle 2 specially */
2307 74053 : GEN F = F2x_factor_i(Flx_to_F2x(f),flag);
2308 74053 : if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
2309 74053 : return F;
2310 : }
2311 11348016 : d = degpol(f);
2312 11347328 : if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
2313 6898976 : switch(flag)
2314 : {
2315 1581948 : default: return Flx_factor_Cantor(f, pp);
2316 5315950 : 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 7829189 : Flx_degfact(GEN f, ulong p)
2323 : {
2324 7829189 : pari_sp av = avma;
2325 7829189 : GEN z = Flx_factor_i(Flx_normalize(f,p),p,1);
2326 7828099 : return gerepilecopy(av, z);
2327 : }
2328 :
2329 : /* T must be squarefree mod p*/
2330 : GEN
2331 1452489 : Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
2332 : {
2333 : GEN XP, D;
2334 : pari_timer ti;
2335 1452489 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2336 1452489 : long i, s, n = get_Flx_degree(T);
2337 1452470 : GEN V = const_vecsmall(n, 0);
2338 1452490 : pari_sp av = avma;
2339 1452490 : T = Flx_get_red_pre(T, p, pi);
2340 1452520 : if (DEBUGLEVEL>=6) timer_start(&ti);
2341 1452520 : XP = Flx_Frobenius_pre(T, p, pi);
2342 1452437 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2343 1452437 : D = Flx_ddf_Shoup(T, XP, p, pi);
2344 1452648 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2345 7477243 : for (i = 1, s = 0; i <= n; i++) { V[i] = degpol(gel(D,i))/i; s += V[i]; }
2346 1452616 : *nb = s; set_avma(av); return V;
2347 : }
2348 :
2349 : long
2350 659476 : Flx_nbfact_Frobenius_pre(GEN T, GEN XP, ulong p, ulong pi)
2351 : {
2352 659476 : pari_sp av = avma;
2353 659476 : long s = ddf_to_nbfact(Flx_ddf_Shoup(T, XP, p, pi));
2354 659480 : 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 659471 : Flx_nbfact_pre(GEN T, ulong p, ulong pi)
2363 : {
2364 659471 : pari_sp av = avma;
2365 659471 : GEN XP = Flx_Frobenius_pre(T, p, pi);
2366 659476 : long n = Flx_nbfact_Frobenius_pre(T, XP, p, pi);
2367 659480 : return gc_long(av,n);
2368 : }
2369 : long
2370 659471 : Flx_nbfact(GEN T, ulong p)
2371 659471 : { 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 4295369 : FpX_factor(GEN f, GEN p)
2431 : {
2432 4295369 : pari_sp av = avma;
2433 : GEN F;
2434 4295369 : switch(ZX_factmod_init(&f, p))
2435 : {
2436 1429803 : case 0: F = F2x_factor_i(f,0);
2437 1429802 : F2xV_to_ZXV_inplace(gel(F,1)); break;
2438 2863165 : case 1: F = Flx_factor_i(f,p[2],0);
2439 2863181 : FlxV_to_ZXV_inplace(gel(F,1)); break;
2440 2431 : default: F = FpX_factor_i(f,p,0); break;
2441 : }
2442 4295441 : return gerepilecopy(av, F);
2443 : }
2444 :
2445 : GEN
2446 681800 : Flx_factor(GEN f, ulong p)
2447 : {
2448 681800 : pari_sp av = avma;
2449 681800 : return gerepilecopy(av, Flx_factor_i(Flx_normalize(f,p),p,0));
2450 : }
2451 : GEN
2452 15045 : F2x_factor(GEN f)
2453 : {
2454 15045 : pari_sp av = avma;
2455 15045 : return gerepilecopy(av, F2x_factor_i(f,0));
2456 : }
|
