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 4479171 : ZX_factmod_init(GEN *F, GEN p)
35 : {
36 4479171 : if (lgefint(p) == 3)
37 : {
38 4476742 : ulong pp = p[2];
39 4476742 : if (pp == 2) { *F = ZX_to_F2x(*F); return 0; }
40 3023157 : *F = ZX_to_Flx(*F, pp);
41 3023216 : if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
42 3023219 : return 1;
43 : }
44 2429 : *F = FpX_red(*F, p);
45 2429 : if (lg(*F) > 3) *F = FpX_normalize(*F, p);
46 2429 : return 2;
47 : }
48 : static GEN
49 670510 : ZX_rootmod_init(GEN F, GEN p)
50 670510 : { 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 2277 : Flx_Xnm1(long sv, long n, ulong p)
71 : {
72 2277 : GEN t = cgetg(n+3, t_VECSMALL);
73 : long i;
74 2277 : t[1] = sv;
75 2277 : t[2] = p - 1;
76 8236 : for (i = 3; i <= n+1; i++) t[i] = 0;
77 2277 : t[i] = 1; return t;
78 : }
79 : /* X^n + 1 */
80 : static GEN
81 2181 : Flx_Xn1(long sv, long n, ulong p)
82 : {
83 2181 : GEN t = cgetg(n+3, t_VECSMALL);
84 : long i;
85 : (void) p;
86 2181 : t[1] = sv;
87 2181 : t[2] = 1;
88 8080 : for (i = 3; i <= n+1; i++) t[i] = 0;
89 2181 : t[i] = 1; return t;
90 : }
91 :
92 : /* assume lg(f) > 3 */
93 : static GEN
94 96014 : Flx_root_mod_2(GEN f)
95 : {
96 96014 : long i, n = lg(f)-1, c = f[2];
97 96014 : int z0 = !c;
98 96014 : c ^= 1; /* c = f[2] + f[n] mod 2, we know f[n] is odd */
99 173078 : for (i=3; i < n; i++) c ^= f[i];
100 : /* c = 0 iff f(1) = 0 (mod 2) */
101 96014 : if (z0) return c? mkvecsmall(0): mkvecsmall2(0, 1);
102 14106 : return c? cgetg(1, t_VECSMALL): mkvecsmall(1);
103 : }
104 : /* assume lg(f) > 3 */
105 : static ulong
106 98 : Flx_oneroot_mod_2(GEN f)
107 : {
108 98 : long i, n, c = f[2];
109 98 : if (!c) return 0;
110 98 : n = lg(f)-1; c = 0; /* = f[2] + f[n] (mod 2); both are odd */
111 196 : for (i=3; i < n; i++) c ^= f[i];
112 98 : return c? 2: 1;
113 : }
114 :
115 : static GEN FpX_roots_i(GEN f, GEN p);
116 :
117 : static int
118 16628207 : 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 436381 : FpX_roots(GEN f, GEN p)
123 : {
124 436381 : pari_sp av = avma;
125 436381 : GEN y; f = ZX_rootmod_init(f, p);
126 436382 : switch(lg(f))
127 : {
128 14 : case 2: pari_err_ROOTS0("FpX_roots");
129 48761 : case 3: return cgetg(1,t_COL);
130 : }
131 387607 : if (typ(f) == t_VECSMALL)
132 : {
133 379360 : ulong pp = p[2];
134 379360 : if (pp == 2)
135 96007 : y = Flx_root_mod_2(f);
136 : else
137 : {
138 283353 : if (!odd(pp)) pari_err_PRIME("FpX_roots", p);
139 283353 : y = Flx_roots_pre(f, pp, SMALL_ULONG(pp)? 0: get_Fl_red(pp));
140 : }
141 379353 : y = Flc_to_ZC(y);
142 : }
143 : else
144 8247 : y = FpX_roots_i(f, p);
145 387593 : return gc_upto(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 14676 : FpX_quad_root(GEN x, GEN p, int unknown)
159 : {
160 14676 : GEN s, D, b = gel(x,3), c = gel(x,2);
161 :
162 14676 : if (absequaliu(p, 2)) {
163 0 : if (!signe(b)) return c;
164 0 : return signe(c)? NULL: gen_1;
165 : }
166 14676 : D = subii(sqri(b), shifti(c,2));
167 14676 : D = remii(D,p);
168 14676 : if (unknown && kronecker(D,p) == -1) return NULL;
169 :
170 14093 : s = Fp_sqrt(D,p);
171 : /* p is not prime, go on and give e.g. maxord a chance to recover */
172 14093 : if (!s) return NULL;
173 14085 : 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 32687176 : Fl_disc_bc(ulong b, ulong c, ulong p)
182 32687176 : { 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 30205119 : Flx_quad_root(GEN x, ulong p, ulong pi, int unknown)
186 : {
187 30205119 : ulong s, b = x[3], c = x[2];
188 30205119 : ulong D = Fl_disc_bc(b, c, p);
189 30178658 : if (unknown && krouu(D,p) == -1) return p;
190 19911316 : s = Fl_sqrt_pre(D, p, pi);
191 20064872 : if (s==~0UL) return p;
192 20064859 : return Fl_halve(Fl_sub(s,b, p), p);
193 : }
194 : static ulong
195 18047563 : Flx_otherroot(GEN x, ulong r, ulong p)
196 18047563 : { 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 4892675 : split_init(struct split_t *S, long max)
203 : {
204 4892675 : S->todo = vectrunc_init(max);
205 4892350 : S->done = vectrunc_init(max);
206 4892126 : }
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 5275489 : 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 5276062 : split_moveto_done(struct split_t *S, long i, GEN t)
224 : {
225 5276062 : long n = lg(S->todo)-1;
226 5276062 : vectrunc_append(S->done, t);
227 5276217 : if (n) gel(S->todo,i) = gel(S->todo, n);
228 5276217 : setlg(S->todo, n);
229 :
230 5276074 : }
231 : /* append t to done */
232 : static void
233 526828 : split_add_done(struct split_t *S, GEN t)
234 526828 : { vectrunc_append(S->done, t); }
235 : /* split todo[i] into a and b */
236 : static void
237 432248 : split_todo(struct split_t *S, long i, GEN a, GEN b)
238 : {
239 432248 : gel(S->todo, i) = a;
240 432248 : split_add(S, b);
241 432242 : }
242 : /* split todo[i] into a and b, moved to done */
243 : static void
244 482546 : split_done(struct split_t *S, long i, GEN a, GEN b)
245 : {
246 482546 : split_moveto_done(S, i, a);
247 482549 : split_add_done(S, b);
248 482550 : }
249 :
250 : /* by splitting, assume p > 2 prime, deg(f) > 0 */
251 : static GEN
252 8247 : FpX_roots_i(GEN f, GEN p)
253 : {
254 : GEN pol, pol0, a, q;
255 : struct split_t S;
256 :
257 8247 : f = FpX_normalize(f, p);
258 8247 : split_init(&S, lg(f)-1);
259 8247 : settyp(S.done, t_COL);
260 8247 : if (ZX_valrem(f, &f)) split_add_done(&S, gen_0);
261 8247 : switch(degpol(f))
262 : {
263 14 : 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 431267 : Flx_cubic_root(GEN ff, ulong p, ulong pi)
322 : {
323 431267 : GEN f = Flx_normalize(ff,p);
324 431264 : 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 431264 : if (pi)
327 : {
328 431267 : t = Fl_mul_pre(a, p3, p, pi);
329 431275 : t2 = Fl_sqr_pre(t, p, pi);
330 431274 : A = Fl_sub(b, Fl_triple(t2, p), p);
331 431274 : B = Fl_sub(c, Fl_mul_pre(t, Fl_add(A, t2, p), p, pi), p);
332 431272 : A3 = Fl_mul_pre(A, p3, p, pi);
333 431272 : 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 431272 : A33 = Fl_powu_pre(A3, 3, p, pi);
345 431265 : D = Fl_add(B2, Fl_double(Fl_double(A33, p), p), p);
346 431259 : S = Fl_neg(B,p);
347 431259 : P = Fl_neg(A3,p);
348 431257 : if (krouu(D,p) >= 0)
349 : {
350 354151 : ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
351 354159 : ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
352 354160 : if (p%3==2) /* 1 solutions */
353 138408 : vS1 = Fl_powu_pre(S1, p - p3, p, pi);
354 : else
355 : {
356 215752 : vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
357 215754 : if (vS1==~0UL) return p; /*0 solutions*/
358 : /*3 solutions*/
359 : }
360 231595 : if (!P) return Fl_sub(vS1, t, p);
361 103762 : vS2 = pi? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): Fl_div(P, vS1, p);
362 103761 : return Fl_sub(Fl_add(vS1,vS2, p), t, p);
363 : }
364 : else
365 : {
366 77107 : pari_sp av = avma;
367 77107 : GEN S1 = mkvecsmall2(Fl_halve(S, p), (p + 1UL) >> 1);
368 77107 : GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
369 : ulong Sa;
370 77110 : if (!vS1) return p; /*0 solutions, p%3==2*/
371 77110 : Sa = vS1[1];
372 77110 : if (p%3==1) /*1 solutions*/
373 : {
374 32959 : ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
375 32959 : if (Fa!=P) Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
376 : }
377 77110 : set_avma(av);
378 77110 : 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 4580626 : 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 4580626 : if (Flx_val(f)) return 0;
441 4579036 : da = degpol(f); f = Flx_normalize(f, p);
442 4578575 : if (da == 1) return Fl_neg(f[2], p);
443 4559224 : PI = pi? pi: get_Fl_red(p); /* PI for Fp, pi for Fp[x] */
444 4559927 : switch(da)
445 : {
446 3789305 : case 2: return Flx_quad_root(f, p, PI, 1);
447 416028 : case 3: if (p>3) return Flx_cubic_root(f, p, PI); /*FALL THROUGH*/
448 : }
449 361272 : if (SMALL_ULONG(p)) pi = 0; /* bilinear ops faster without Fl_*_pre */
450 361272 : if (!fl)
451 : {
452 319626 : a = Flxq_powu_pre(polx_Flx(f[1]), p - 1, f,p,pi);
453 319555 : if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
454 319555 : a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
455 319548 : a = Flx_gcd_pre(f,a, p, pi);
456 41646 : } else a = f;
457 361226 : da = degpol(a);
458 361222 : if (!da) return p;
459 259722 : a = Flx_normalize(a,p);
460 :
461 259752 : q = p >> 1;
462 259752 : pol = polx_Flx(f[1]);
463 393207 : for(pol[2] = 1;; pol[2]++)
464 : {
465 393207 : if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
466 393221 : switch(da)
467 : {
468 162991 : case 1: return Fl_neg(a[2], p);
469 74886 : case 2: return Flx_quad_root(a, p, PI, 0);
470 21917 : case 3: if (p>3) return Flx_cubic_root(a, p, PI); /*FALL THROUGH*/
471 : default: {
472 133427 : GEN b = Flxq_powu_pre(pol,q, a,p,pi);
473 : long db;
474 133432 : if (degpol(b) <= 0) continue;
475 126983 : b = Flx_gcd_pre(a,Flx_Fl_add(b,p-1,p), p, pi);
476 126994 : db = degpol(b); if (!db) continue;
477 126993 : b = Flx_normalize(b, p);
478 126998 : if (db <= (da >> 1)) {
479 77554 : a = b;
480 77554 : da = db;
481 : } else {
482 49444 : a = Flx_div_pre(a,b, p, pi);
483 49449 : da -= db;
484 : }
485 : }
486 : }
487 : }
488 : }
489 : ulong
490 4537614 : Flx_oneroot_pre(GEN f, ulong p, ulong pi)
491 4537614 : { return Flx_oneroot_pre_i(f, p, pi, 0); }
492 : ulong
493 42962 : Flx_oneroot_split_pre(GEN f, ulong p, ulong pi)
494 42962 : { return Flx_oneroot_pre_i(f, p, pi, 1); }
495 :
496 : /* assume p > 3 prime */
497 : static GEN
498 5323 : FpX_oneroot_i(GEN f, GEN p)
499 : {
500 : GEN pol, pol0, a, q;
501 : long da;
502 :
503 5323 : if (ZX_val(f)) return gen_0;
504 4973 : f = FpX_normalize(f, p);
505 4973 : switch(degpol(f))
506 : {
507 737 : case 1: return subii(p, gel(f,2));
508 4047 : 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 2520 : Flx_oneroot(GEN f, ulong p)
551 : {
552 2520 : pari_sp av = avma;
553 2520 : switch(lg(f))
554 : {
555 0 : case 2: return 0;
556 0 : case 3: return p;
557 : }
558 2520 : if (p == 2) return Flx_oneroot_mod_2(f);
559 2520 : 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 234129 : FpX_oneroot(GEN f, GEN p)
578 : {
579 234129 : pari_sp av = avma;
580 234129 : f = ZX_rootmod_init(f, p);
581 234129 : switch(lg(f))
582 : {
583 0 : case 2: set_avma(av); return gen_0;
584 0 : case 3: return gc_NULL(av);
585 : }
586 234129 : if (typ(f) == t_VECSMALL)
587 : {
588 228806 : ulong r, pp = p[2];
589 228806 : if (pp == 2)
590 98 : r = Flx_oneroot_mod_2(f);
591 : else
592 228708 : r = Flx_oneroot_pre(f, pp, SMALL_ULONG(pp)? 0: get_Fl_red(pp));
593 228806 : set_avma(av);
594 228806 : return (r == pp)? NULL: utoi(r);
595 : }
596 5323 : f = FpX_oneroot_i(f, p);
597 5322 : if (!f) return gc_NULL(av);
598 5322 : return gc_INT(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 105 : good_root_of_unity(GEN p, long deg, long deg_factor, long *pt_n)
608 : {
609 105 : pari_sp ltop = avma;
610 : GEN pm, factn, power, base, zeta;
611 : long n;
612 :
613 105 : pm = subis (p, 1ul);
614 224 : for (n = deg / 2 / deg_factor + 1; !dvdiu (pm, n); n--);
615 105 : factn = Z_factor(stoi(n));
616 105 : power = diviuexact (pm, n);
617 105 : base = gen_1;
618 : do {
619 140 : base = addis (base, 1l);
620 140 : zeta = Fp_pow (base, power, p);
621 : }
622 140 : while (!equaliu (Fp_order (zeta, factn, p), n));
623 105 : *pt_n = n;
624 105 : return gc_INT (ltop, zeta);
625 : }
626 :
627 : GEN
628 1043 : FpX_oneroot_split(GEN fact, GEN p)
629 : {
630 1043 : pari_sp av = avma;
631 : long n, deg_f, i, dmin;
632 : GEN prim, expo, minfactor, xplusa, zeta, xpow;
633 1043 : fact = FpX_normalize(fact, p);
634 1043 : deg_f = degpol(fact);
635 1043 : if (deg_f <= 3) return FpX_oneroot(fact, p);
636 77 : minfactor = fact; /* factor of minimal degree found so far */
637 77 : dmin = degpol(minfactor);
638 77 : xplusa = pol_x(varn(fact));
639 182 : 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 105 : fact = minfactor; deg_f = dmin;
644 105 : zeta = gen_1;
645 105 : prim = good_root_of_unity(p, deg_f, 1, &n);
646 105 : expo = diviuexact(subiu(p, 1), n);
647 : /* update X+a, avoid a=0 */
648 105 : gel (xplusa, 2) = addis (gel (xplusa, 2), 1);
649 105 : xpow = FpXQ_pow (xplusa, expo, fact, p);
650 266 : for (i = 0; i < n; i++)
651 : {
652 196 : GEN tmp = FpX_gcd(FpX_Fp_sub(xpow, zeta, p), fact, p);
653 196 : long dtmp = degpol(tmp);
654 196 : if (dtmp > 0 && dtmp < deg_f)
655 : {
656 112 : fact = FpX_div(fact, tmp, p); deg_f = degpol(fact);
657 112 : if (dtmp < dmin)
658 : {
659 105 : minfactor = FpX_normalize (tmp, p);
660 105 : dmin = dtmp;
661 105 : if (dmin == 1 || dmin <= (2 * deg_f) / n - 1)
662 : /* stop early to avoid too many gcds */
663 : break;
664 : }
665 : }
666 161 : zeta = Fp_mul (zeta, prim, p);
667 : }
668 : }
669 77 : return gc_INT(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 1922348 : FE_concat(GEN F, GEN E, long l)
682 : {
683 1922348 : setlg(E,l); E = shallowconcat1(E);
684 1922345 : 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 5351289 : vddf_to_simplefact(GEN V, long d)
716 : {
717 : GEN E, F;
718 5351289 : long i, j, c, l = lg(V);
719 5351289 : F = cgetg(d+1, t_VECSMALL);
720 5350281 : E = cgetg(d+1, t_VECSMALL);
721 10783350 : for (i = c = 1; i < l; i++)
722 : {
723 5433507 : GEN Vi = gel(V,i);
724 5433507 : long l = lg(Vi);
725 27385091 : for (j = 1; j < l; j++)
726 : {
727 21952092 : long k, n = degpol(gel(Vi,j)) / j;
728 33252821 : for (k = 1; k <= n; k++) { uel(F,c) = j; uel(E,c) = i; c++; }
729 : }
730 : }
731 5349843 : setlg(F,c);
732 5349835 : setlg(E,c);
733 5349941 : 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 251953 : FpX_split_part(GEN f, GEN p)
739 : {
740 251953 : long n = degpol(f);
741 251953 : GEN z, X = pol_x(varn(f));
742 251953 : if (n <= 1) return f;
743 249390 : f = FpX_red(f, p);
744 249390 : z = FpX_sub(FpX_Frobenius(f, p), X, p);
745 249388 : 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 154016 : FpX_nbroots(GEN f, GEN p)
752 : {
753 154016 : pari_sp av = avma;
754 154016 : GEN z = FpX_split_part(f, p);
755 154016 : return gc_long(av, degpol(z));
756 : }
757 :
758 : /* 1 < deg(f) <= p */
759 : static int
760 82165 : Flx_is_totally_split_i(GEN f, ulong p)
761 : {
762 82165 : GEN F = Flx_Frobenius(f, p);
763 82165 : return degpol(F)==1 && uel(F,2)==0UL && uel(F,3)==1UL;
764 : }
765 : int
766 82172 : Flx_is_totally_split(GEN f, ulong p)
767 : {
768 82172 : pari_sp av = avma;
769 82172 : ulong n = degpol(f);
770 82172 : if (n <= 1) return 1;
771 82165 : if (n > p) return 0; /* includes n < 0 */
772 82165 : 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 4955421 : Flx_nbroots(GEN f, ulong p)
794 : {
795 4955421 : long n = degpol(f);
796 : ulong pi;
797 4955415 : pari_sp av = avma;
798 : GEN z;
799 4955415 : if (n <= 1) return n;
800 4937488 : if (n == 2)
801 : {
802 : ulong D;
803 926279 : if (p==2) return (f[2]==0) + (f[2]!=f[3]);
804 691822 : D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
805 691823 : return 1 + krouu(D,p);
806 : }
807 4011209 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
808 4011209 : z = Flx_sub(Flx_Frobenius_pre(f, p, pi), polx_Flx(f[1]), p);
809 4011191 : z = Flx_gcd_pre(z, f, p, pi);
810 4011218 : 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 1022 : 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 1022 : n = get_FpX_degree(T); v = get_FpX_var(T);
863 1022 : if (n == 0) return cgetg(1, t_VEC);
864 1022 : if (n == 1) return mkvec(get_FpX_mod(T));
865 845 : B = n/2;
866 845 : l = usqrt(B);
867 845 : m = (B+l-1)/l;
868 845 : T = FpX_get_red(T, p);
869 845 : b = cgetg(l+2, t_VEC);
870 845 : gel(b, 1) = pol_x(v);
871 845 : gel(b, 2) = XP;
872 845 : if (DEBUGLEVEL>=7) timer_start(&ti);
873 845 : xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1), T, p);
874 845 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq baby");
875 1071 : for (i = 3; i <= l+1; i++)
876 226 : gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
877 845 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: baby");
878 845 : xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1), T, p);
879 845 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq giant");
880 845 : g = cgetg(m+1, t_VEC);
881 845 : gel(g, 1) = gel(xq, 2);
882 1836 : for(i = 2; i <= m; i++) gel(g, i) = FpX_FpXQV_eval(gel(g, i-1), xq, T, p);
883 845 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: giant");
884 845 : h = cgetg(m+1, t_VEC);
885 2681 : for (j = 1; j <= m; j++)
886 : {
887 1836 : pari_sp av = avma;
888 1836 : GEN gj = gel(g,j), e = FpX_sub(gj, gel(b,1), p);
889 3063 : for (i = 2; i <= l; i++) e = FpXQ_mul(e, FpX_sub(gj, gel(b,i), p), T, p);
890 1836 : gel(h,j) = gc_upto(av, e);
891 : }
892 845 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: diff");
893 845 : Tr = get_FpX_mod(T);
894 845 : F = cgetg(m+1, t_VEC);
895 2681 : for (j = 1; j <= m; j++)
896 : {
897 1836 : GEN u = FpX_gcd(Tr, gel(h,j), p);
898 1836 : if (degpol(u))
899 : {
900 478 : u = FpX_normalize(u, p);
901 478 : Tr = FpX_div(Tr, u, p);
902 : }
903 1836 : gel(F,j) = u;
904 : }
905 845 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: F");
906 845 : f = const_vec(n, pol_1(v));
907 2681 : for (j = 1; j <= m; j++)
908 : {
909 1836 : GEN e = gel(F, j);
910 1908 : for (i=l-1; i >= 0; i--)
911 : {
912 1908 : GEN u = FpX_gcd(e, FpX_sub(gel(g, j), gel(b, i+1), p), p);
913 1908 : if (degpol(u))
914 : {
915 508 : u = FpX_normalize(u, p);
916 508 : gel(f, l*j-i) = u;
917 508 : e = FpX_div(e, u, p);
918 : }
919 1908 : if (!degpol(e)) break;
920 : }
921 : }
922 845 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: f");
923 845 : if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
924 845 : 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 1051 : 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 1051 : GEN Tp = get_FpX_mod(T);
966 1051 : long n = degpol(hp), vT = varn(Tp), ct = 0;
967 : GEN u1, u2, f1, f2, R, h;
968 1051 : h = FpX_get_red(hp, p);
969 1051 : t = FpX_rem(t, T, p);
970 1051 : av = avma;
971 : do
972 : {
973 1755 : set_avma(av);
974 1755 : R = FpXQ_pow(deg1pol_shallow(gen_1, randomi(p), vT), p2, h, p);
975 1755 : u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
976 1755 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_rec",p);
977 1755 : } while (degpol(u1)==0 || degpol(u1)==n);
978 1051 : f1 = FpX_gcd(FpX_FpXQ_eval(u1, t, T, p), Tp, p);
979 1051 : f1 = FpX_normalize(f1, p);
980 1051 : u2 = FpX_div(hp, u1, p);
981 1051 : f2 = FpX_div(Tp, f1, p);
982 1051 : if (degpol(u1)==1)
983 796 : gel(V, idx) = f1;
984 : else
985 255 : FpX_edf_rec(FpX_get_red(f1, p), u1, t, d, p2, p, V, idx);
986 1051 : idx += degpol(f1)/d;
987 1051 : if (degpol(u2)==1)
988 689 : gel(V, idx) = f2;
989 : else
990 362 : FpX_edf_rec(FpX_get_red(f2, p), u2, t, d, p2, p, V, idx);
991 1051 : }
992 :
993 : /* assume Tp a squarefree product of r > 1 irred. factors of degree d */
994 : static void
995 434 : FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
996 : {
997 434 : long n = degpol(Tp), r = n/d, vT = varn(Tp), ct = 0;
998 : GEN T, h, t;
999 : pari_timer ti;
1000 :
1001 434 : T = FpX_get_red(Tp, p);
1002 434 : XP = FpX_rem(XP, T, p);
1003 434 : if (DEBUGLEVEL>=7) timer_start(&ti);
1004 : do
1005 : {
1006 434 : GEN g = random_FpX(n, vT, p);
1007 434 : t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
1008 434 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
1009 434 : h = FpXQ_minpoly(t, T, p);
1010 434 : if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
1011 434 : if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf",p);
1012 434 : } while (degpol(h) != r);
1013 434 : FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
1014 434 : }
1015 :
1016 : static GEN
1017 1008 : FpX_factor_Shoup(GEN T, GEN p)
1018 : {
1019 1008 : long i, n, s = 0;
1020 : GEN XP, D, V;
1021 1008 : long e = expi(p);
1022 : pari_timer ti;
1023 1008 : n = get_FpX_degree(T);
1024 1008 : T = FpX_get_red(T, p);
1025 1008 : if (DEBUGLEVEL>=6) timer_start(&ti);
1026 1008 : XP = FpX_Frobenius(T, p);
1027 1008 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
1028 1008 : D = FpX_ddf_Shoup(T, XP, p);
1029 1008 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
1030 1008 : s = ddf_to_nbfact(D);
1031 1008 : V = cgetg(s+1, t_COL);
1032 7272 : for (i = 1, s = 1; i <= n; i++)
1033 : {
1034 6264 : GEN Di = gel(D,i);
1035 6264 : long ni = degpol(Di), ri = ni/i;
1036 6264 : if (ni == 0) continue;
1037 1126 : Di = FpX_normalize(Di, p);
1038 1126 : if (ni == i) { gel(V, s++) = Di; continue; }
1039 434 : if (ri <= e*expu(e))
1040 434 : FpX_edf(Di, XP, i, p, V, s);
1041 : else
1042 0 : FpX_edf_simple(Di, XP, i, p, V, s);
1043 434 : if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
1044 434 : s += ri;
1045 : }
1046 1008 : return V;
1047 : }
1048 :
1049 : long
1050 2372334 : ddf_to_nbfact(GEN D)
1051 : {
1052 2372334 : long l = lg(D), i, s = 0;
1053 14970927 : for(i = 1; i < l; i++) s += degpol(gel(D,i))/i;
1054 2372329 : return s;
1055 : }
1056 :
1057 : /* Yun algorithm: Assume p > degpol(T) */
1058 : static GEN
1059 1713 : FpX_factor_Yun(GEN T, GEN p)
1060 : {
1061 1713 : long n = degpol(T), i = 1;
1062 1713 : GEN a, b, c, d = FpX_deriv(T, p);
1063 1713 : GEN V = cgetg(n+1,t_VEC);
1064 1713 : a = FpX_gcd(T, d, p);
1065 1713 : if (degpol(a) == 0) return mkvec(T);
1066 603 : b = FpX_div(T, a, p);
1067 : do
1068 : {
1069 2739 : c = FpX_div(d, a, p);
1070 2739 : d = FpX_sub(c, FpX_deriv(b, p), p);
1071 2739 : a = FpX_normalize(FpX_gcd(b, d, p), p);
1072 2739 : gel(V, i++) = a;
1073 2739 : b = FpX_div(b, a, p);
1074 2739 : } while (degpol(b));
1075 603 : setlg(V, i); return V;
1076 : }
1077 : GEN
1078 349888 : FpX_factor_squarefree(GEN T, GEN p)
1079 : {
1080 349888 : if (lgefint(p)==3)
1081 : {
1082 349140 : ulong pp = (ulong)p[2];
1083 349140 : GEN u = Flx_factor_squarefree(ZX_to_Flx(T,pp), pp);
1084 349140 : return FlxV_to_ZXV(u);
1085 : }
1086 748 : return FpX_factor_Yun(T, p);
1087 : }
1088 :
1089 : GEN
1090 227563 : FpX_roots_mult(GEN T, long n, GEN p)
1091 : {
1092 227563 : pari_sp av = avma;
1093 227563 : GEN V = FpX_factor_squarefree(T, p), W;
1094 227563 : long l = lg(V), i;
1095 227563 : if (l <= n) retgc_const(av, cgetg(1, t_COL));
1096 40864 : W = cgetg(l-n+1,t_VEC);
1097 129824 : for (i = n; i < l; i++)
1098 88960 : gel(W,i-n+1) = FpX_roots(gel(V,i), p);
1099 40864 : return gc_upto(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 = gc_upto(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 = gc_upto(av, r);
1132 0 : } else av = avma;
1133 35 : return 1;
1134 : }
1135 :
1136 : static GEN
1137 909 : FpX_factor_Cantor(GEN T, GEN p)
1138 : {
1139 909 : GEN E, F, V = FpX_factor_Yun(T, p);
1140 909 : long i, j, l = lg(V);
1141 909 : F = cgetg(l, t_VEC);
1142 909 : E = cgetg(l, t_VEC);
1143 2444 : for (i=1, j=1; i < l; i++)
1144 1535 : if (degpol(gel(V,i)))
1145 : {
1146 1008 : GEN Fj = FpX_factor_Shoup(gel(V,i), p);
1147 1008 : gel(F, j) = Fj;
1148 1008 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1149 1008 : j++;
1150 : }
1151 909 : 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 gc_GEN(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 2429 : FpX_factor_i(GEN f, GEN pp, long flag)
1222 : {
1223 2429 : long d = degpol(f);
1224 2429 : if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
1225 923 : switch(flag)
1226 : {
1227 909 : 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 354268 : trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
1283 :
1284 : /* not GC 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 472467 : F2x_quad_factortype(GEN x)
1321 472467 : { 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 19684 : F2x_degfact_2(GEN f, long d)
1329 : {
1330 19684 : if (!d) return trivial_fact();
1331 19684 : if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1332 19488 : switch(F2x_quad_factortype(f)) {
1333 5964 : case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1334 6237 : case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1335 7287 : default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1336 : }
1337 : }
1338 :
1339 : static GEN
1340 1188717 : F2x_factor_2(GEN f, long d)
1341 : {
1342 1188717 : long v = f[1];
1343 1188717 : if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1344 932375 : if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1345 448101 : switch(F2x_quad_factortype(f))
1346 : {
1347 72019 : case -1: return mkvec2(mkcol(f), mkvecsmall(1));
1348 329537 : case 0: return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
1349 46559 : default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
1350 : }
1351 : }
1352 : static GEN
1353 1208416 : F2x_factor_deg2(GEN f, long d, long flag)
1354 : {
1355 1208416 : switch(flag) {
1356 14 : case 2: return F2x_is_irred_2(f, d);
1357 19684 : case 1: return F2x_degfact_2(f, d);
1358 1188718 : default: return F2x_factor_2(f, d);
1359 : }
1360 : }
1361 :
1362 : /* xt = NULL or x^(p-1)/2 mod g */
1363 : static void
1364 31667 : split_squares(struct split_t *S, GEN g, ulong p, ulong pi, GEN xt)
1365 : {
1366 31667 : ulong q = p >> 1;
1367 31667 : GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
1368 31667 : long d = degpol(a);
1369 31667 : if (d < 0)
1370 : {
1371 : ulong i;
1372 428 : split_add_done(S, (GEN)1);
1373 428 : if (!pi)
1374 1618 : 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 31239 : if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
1379 31239 : if (d)
1380 : {
1381 31150 : if (xt) xt = Flx_Fl_add(xt, p-1, p); else xt = Flx_Xnm1(g[1], q, p);
1382 31150 : a = Flx_gcd_pre(a, xt, p, pi);
1383 31150 : if (degpol(a)) split_add(S, Flx_normalize(a, p));
1384 : }
1385 : }
1386 31667 : }
1387 : static void
1388 31667 : split_nonsquares(struct split_t *S, GEN g, ulong p, ulong pi, GEN xt)
1389 : {
1390 31667 : ulong q = p >> 1;
1391 31667 : GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
1392 31667 : long d = degpol(a);
1393 31667 : 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 31418 : if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
1405 31418 : if (d)
1406 : {
1407 31054 : if (xt) xt = Flx_Fl_add(xt, 1, p); else xt = Flx_Xn1(g[1], q, p);
1408 31054 : a = Flx_gcd_pre(a, xt, p, pi);
1409 31054 : if (degpol(a)) split_add(S, Flx_normalize(a, p));
1410 : }
1411 : }
1412 31667 : }
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 4883964 : split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p, ulong pi)
1417 : {
1418 4883964 : GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
1419 4883972 : long d = degpol(g);
1420 4883978 : if (d < 0) return 0;
1421 4883378 : if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
1422 4883377 : if (!d) return 1;
1423 4849107 : if ((p >> 4) <= (ulong)d)
1424 : { /* small p; split directly using x^((p-1)/2) +/- 1 */
1425 28873 : GEN xt = ((ulong)d < (p>>1))? Flx_rem_pre(monomial_Flx(1, p>>1, g[1]), g, p, pi)
1426 31667 : : NULL;
1427 31667 : split_squares(S, g, p, pi, xt);
1428 31667 : split_nonsquares(S, g, p, pi, xt);
1429 : } else { /* large p; use x^(p-1) - 1 directly */
1430 4817440 : a = Flxq_powu_pre(polx_Flx(f[1]), p-1, g, p, pi);
1431 4806599 : if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
1432 4806599 : a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
1433 4807036 : g = Flx_gcd_pre(g,a, p,pi);
1434 4812141 : if (degpol(g)) split_add(S, Flx_normalize(g,p));
1435 : }
1436 4848537 : return 1;
1437 : }
1438 :
1439 : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
1440 : GEN
1441 29225657 : Flx_roots_pre(GEN f, ulong p, ulong pi)
1442 : {
1443 : GEN pol;
1444 29225657 : long v = Flx_valrem(f, &f), n = degpol(f);
1445 : ulong q, PI;
1446 : struct split_t S;
1447 :
1448 29216023 : f = Flx_normalize(f, p);
1449 : /* optimization: test for small degree first */
1450 29237172 : if (n == 1)
1451 : {
1452 115418 : q = p - f[2];
1453 115418 : return v? mkvecsmall2(0, q): mkvecsmall(q);
1454 : }
1455 29121754 : PI = pi? pi: get_Fl_red(p); /* PI for Fp, pi for Fp[x] */
1456 29125393 : if (n == 2)
1457 : {
1458 24240520 : ulong r = Flx_quad_root(f, p, PI, 1), s;
1459 24326685 : if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
1460 16537163 : s = Flx_otherroot(f,r, p);
1461 16593100 : if (r < s)
1462 4169956 : return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
1463 12423144 : else if (r > s)
1464 12464081 : return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
1465 : else
1466 4886 : return v? mkvecsmall2(0, s): mkvecsmall(s);
1467 : }
1468 4884873 : if (SMALL_ULONG(p)) pi = 0; /* bilinear ops faster without Fl_*_pre */
1469 4884873 : q = p >> 1;
1470 4884873 : split_init(&S, lg(f)-1);
1471 4883946 : settyp(S.done, t_VECSMALL);
1472 4883946 : if (v) split_add_done(&S, (GEN)0);
1473 4883946 : if (! split_Flx_cut_out_roots(&S, f, p, pi))
1474 600 : return all_roots_mod_p(p, lg(S.done) == 1);
1475 4882696 : pol = polx_Flx(f[1]);
1476 4883451 : for (pol[2]=1; ; pol[2]++)
1477 5303740 : {
1478 10187191 : long j, l = lg(S.todo);
1479 10187191 : if (l == 1) { vecsmall_sort(S.done); return S.done; }
1480 5303734 : if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
1481 11113537 : for (j = 1; j < l; j++)
1482 : {
1483 5809797 : GEN b, c = gel(S.todo,j);
1484 : ulong r, s;
1485 5809797 : switch(degpol(c))
1486 : {
1487 4789051 : case 1:
1488 4789051 : split_moveto_done(&S, j, (GEN)(p - c[2]));
1489 4789205 : j--; l--; break;
1490 476318 : case 2:
1491 476318 : r = Flx_quad_root(c, p, PI, 0);
1492 476352 : if (r == p) pari_err_PRIME("polrootsmod",utoipos(p));
1493 476345 : s = Flx_otherroot(c,r, p);
1494 476337 : split_done(&S, j, (GEN)r, (GEN)s);
1495 476338 : j--; l--; break;
1496 544314 : default:
1497 544314 : b = Flxq_powu_pre(pol,q, c,p,pi); /* pol^(p-1)/2 */
1498 544483 : if (degpol(b) <= 0) continue;
1499 426491 : b = Flx_gcd_pre(c,Flx_Fl_add(b,p-1,p), p, pi);
1500 426539 : if (!degpol(b)) continue;
1501 425988 : b = Flx_normalize(b, p);
1502 426027 : c = Flx_div_pre(c,b, p,pi);
1503 426007 : 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: retgc_const(av, 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 gc_leaf(av, Flx_roots_pre(f, p, pi));
1522 : }
1523 :
1524 : /* assume x reduced mod p, monic. */
1525 : static int
1526 2489970 : Flx_quad_factortype(GEN x, ulong p)
1527 : {
1528 2489970 : ulong b = x[3], c = x[2];
1529 2489970 : 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 2516017 : Flx_degfact_2(GEN f, ulong p, long d)
1540 : {
1541 2516017 : if (!d) return trivial_fact();
1542 2516017 : if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1543 2489977 : switch(Flx_quad_factortype(f, p)) {
1544 1180453 : case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1545 1279232 : case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1546 30120 : default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1547 : }
1548 : }
1549 : /* p > 2 */
1550 : static GEN
1551 2043190 : Flx_factor_2(GEN f, ulong p, long d)
1552 : {
1553 : ulong r, s;
1554 : GEN R,S;
1555 2043190 : long v = f[1];
1556 2043190 : if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1557 1965735 : if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1558 1644774 : r = Flx_quad_root(f, p, get_Fl_red(p), 1);
1559 1644823 : if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
1560 964723 : s = Flx_otherroot(f, r, p);
1561 964719 : r = Fl_neg(r, p);
1562 964716 : s = Fl_neg(s, p);
1563 964714 : if (s < r) lswap(s,r);
1564 964714 : R = mkvecsmall3(v,r,1);
1565 964718 : if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
1566 851513 : S = mkvecsmall3(v,s,1);
1567 851514 : return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
1568 : }
1569 : static GEN
1570 4559119 : Flx_factor_deg2(GEN f, ulong p, long d, long flag)
1571 : {
1572 4559119 : switch(flag) {
1573 56 : case 2: return Flx_is_irred_2(f, p, d);
1574 2516022 : case 1: return Flx_degfact_2(f, p, d);
1575 2043041 : default: return Flx_factor_2(f, p, d);
1576 : }
1577 : }
1578 :
1579 : static GEN
1580 23447 : F2x_Berlekamp_ker(GEN u)
1581 : {
1582 23447 : pari_sp ltop=avma;
1583 23447 : long j,N = F2x_degree(u);
1584 : GEN Q;
1585 : pari_timer T;
1586 23447 : timer_start(&T);
1587 23447 : Q = F2x_matFrobenius(u);
1588 318661 : for (j=1; j<=N; j++)
1589 295214 : F2m_flip(Q,j,j);
1590 23447 : if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
1591 23447 : Q = F2m_ker_sp(Q,0);
1592 23447 : if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
1593 23447 : return gc_upto(ltop,Q);
1594 : }
1595 : #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
1596 : static long
1597 35966 : F2x_split_Berlekamp(GEN *t)
1598 : {
1599 35966 : GEN u = *t, a, b, vker;
1600 35966 : long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
1601 :
1602 35965 : if (du == 1) return 1;
1603 27499 : if (du == 2)
1604 : {
1605 4053 : 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 4053 : return 1;
1612 : }
1613 :
1614 23446 : vker = F2x_Berlekamp_ker(u);
1615 23447 : lb = lgcols(vker);
1616 23461 : d = lg(vker)-1;
1617 23461 : ir = 0;
1618 : /* t[i] irreducible for i < ir, still to be treated for i < L */
1619 60282 : for (L=1; L<d; )
1620 : {
1621 : GEN pol;
1622 36837 : if (d == 2)
1623 5554 : pol = F2v_to_F2x(gel(vker,2), sv);
1624 : else
1625 : {
1626 31283 : GEN v = zero_zv(lb);
1627 31292 : v[1] = du;
1628 31292 : v[2] = random_Fl(2); /*Assume vker[1]=1*/
1629 117795 : for (i=2; i<=d; i++)
1630 86501 : if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
1631 31294 : pol = F2v_to_F2x(v, sv);
1632 : }
1633 110572 : for (i=ir; i<L && L<d; i++)
1634 : {
1635 73751 : a = t[i]; la = F2x_degree(a);
1636 73741 : if (la == 1) { set_irred(i); }
1637 73511 : else if (la == 2)
1638 : {
1639 805 : 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 805 : set_irred(i);
1645 : }
1646 : else
1647 : {
1648 72706 : pari_sp av = avma;
1649 : long lb;
1650 72706 : b = F2x_rem(pol, a);
1651 72688 : if (F2x_degree(b) <= 0) { set_avma(av); continue; }
1652 26593 : b = F2x_gcd(a,b); lb = F2x_degree(b);
1653 26595 : if (lb && lb < la)
1654 : {
1655 26595 : t[L] = F2x_div(a,b);
1656 26594 : t[i]= b; L++;
1657 : }
1658 0 : else set_avma(av);
1659 : }
1660 : }
1661 : }
1662 23445 : return d;
1663 : }
1664 : /* assume deg f > 2 */
1665 : static GEN
1666 35522 : F2x_Berlekamp_i(GEN f, long flag)
1667 : {
1668 35522 : long lfact, val, d = F2x_degree(f), j, k, lV;
1669 : GEN y, E, t, V;
1670 :
1671 35522 : val = F2x_valrem(f, &f);
1672 35523 : if (flag == 2 && val) return NULL;
1673 35509 : V = F2x_factor_squarefree(f); lV = lg(V);
1674 35508 : if (flag == 2 && lV > 2) return NULL;
1675 :
1676 : /* to hold factors and exponents */
1677 35438 : t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
1678 35439 : E = cgetg(d+1,t_VECSMALL);
1679 35439 : lfact = 1;
1680 35439 : if (val) {
1681 11274 : if (flag == 1) t[1] = 1; else gel(t,1) = polx_F2x(f[1]);
1682 11274 : E[1] = val; lfact++;
1683 : }
1684 :
1685 118306 : for (k=1; k<lV; k++)
1686 : {
1687 83023 : if (F2x_degree(gel(V, k))==0) continue;
1688 35966 : gel(t,lfact) = gel(V, k);
1689 35966 : d = F2x_split_Berlekamp(&gel(t,lfact));
1690 35964 : if (flag == 2 && d != 1) return NULL;
1691 35810 : if (flag == 1)
1692 47895 : for (j=0; j<d; j++) t[lfact+j] = F2x_degree(gel(t,lfact+j));
1693 97803 : for (j=0; j<d; j++) E[lfact+j] = k;
1694 35810 : lfact += d;
1695 : }
1696 35283 : if (flag == 2) return gen_1; /* irreducible */
1697 35269 : setlg(t, lfact);
1698 35268 : setlg(E, lfact); y = mkvec2(t,E);
1699 24399 : return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
1700 59669 : : sort_factor_pol(y, cmpGuGu);
1701 : }
1702 :
1703 : /* Adapted from Shoup NTL */
1704 : GEN
1705 339939 : F2x_factor_squarefree(GEN f)
1706 : {
1707 : GEN r, t, v, tv;
1708 339939 : long i, q, n = F2x_degree(f);
1709 339939 : GEN u = const_vec(n+1, pol1_F2x(f[1]));
1710 339939 : for(q = 1;;q *= 2)
1711 : {
1712 556047 : r = F2x_gcd(f, F2x_deriv(f));
1713 556041 : if (F2x_degree(r) == 0)
1714 : {
1715 245508 : gel(u, q) = f;
1716 245508 : break;
1717 : }
1718 310538 : t = F2x_div(f, r);
1719 310541 : if (F2x_degree(t) > 0)
1720 : {
1721 : long j;
1722 146854 : for(j = 1;;j++)
1723 : {
1724 352518 : v = F2x_gcd(r, t);
1725 352518 : tv = F2x_div(t, v);
1726 352516 : if (F2x_degree(tv) > 0)
1727 148826 : gel(u, j*q) = tv;
1728 352517 : if (F2x_degree(v) <= 0) break;
1729 205663 : r = F2x_div(r, v);
1730 205664 : t = v;
1731 : }
1732 146853 : if (F2x_degree(r) == 0) break;
1733 : }
1734 216105 : f = F2x_sqrt(r);
1735 : }
1736 1209893 : for (i = n; i; i--)
1737 1207217 : if (F2x_degree(gel(u,i))) break;
1738 339919 : setlg(u,i+1); return u;
1739 : }
1740 :
1741 : static GEN
1742 355537 : F2x_ddf_simple(GEN T, GEN XP)
1743 : {
1744 355537 : pari_sp av = avma, av2;
1745 : GEN f, z, Tr, X;
1746 355537 : long j, n = F2x_degree(T), v = T[1], B = n/2;
1747 355537 : if (n == 0) return cgetg(1, t_VEC);
1748 355539 : if (n == 1) return mkvec(T);
1749 139888 : z = XP; Tr = T; X = polx_F2x(v);
1750 139888 : f = const_vec(n, pol1_F2x(v));
1751 139891 : av2 = avma;
1752 360816 : for (j = 1; j <= B; j++)
1753 : {
1754 233707 : GEN u = F2x_gcd(Tr, F2x_add(z, X));
1755 233681 : if (F2x_degree(u))
1756 : {
1757 52355 : gel(f, j) = u;
1758 52355 : Tr = F2x_div(Tr, u);
1759 52349 : av2 = avma;
1760 181386 : } else z = gc_leaf(av2, z);
1761 233758 : if (!F2x_degree(Tr)) break;
1762 220922 : z = F2xq_sqr(z, Tr);
1763 : }
1764 139887 : if (F2x_degree(Tr)) gel(f, F2x_degree(Tr)) = Tr;
1765 139886 : return gc_GEN(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 23295 : F2xq_frobtrace(GEN a, long d, GEN T)
1779 : {
1780 23295 : pari_sp av = avma;
1781 : long i;
1782 23295 : GEN x = a;
1783 61919 : for(i=1; i<d; i++)
1784 : {
1785 38626 : x = F2x_add(a, F2xq_sqr(x,T));
1786 38624 : if (gc_needed(av, 2))
1787 0 : x = gc_leaf(av, x);
1788 : }
1789 23293 : return x;
1790 : }
1791 :
1792 : static void
1793 34506 : F2x_edf_simple(GEN Tp, GEN XP, long d, GEN V, long idx)
1794 : {
1795 34506 : long n = F2x_degree(Tp), r = n/d;
1796 : GEN T, f, ff;
1797 34506 : if (r==1) { gel(V, idx) = Tp; return; }
1798 11597 : T = Tp;
1799 11597 : XP = F2x_rem(XP, T);
1800 : while (1)
1801 11697 : {
1802 23295 : pari_sp btop = avma;
1803 : long df;
1804 23295 : GEN g = random_F2x(n, Tp[1]);
1805 23295 : GEN t = F2xq_frobtrace(g, d, T);
1806 23295 : if (lgpol(t) == 0) continue;
1807 17432 : f = F2x_gcd(t, Tp); df = F2x_degree(f);
1808 17432 : if (df > 0 && df < n) break;
1809 5834 : set_avma(btop);
1810 : }
1811 11598 : ff = F2x_div(Tp, f);
1812 11598 : F2x_edf_simple(f, XP, d, V, idx);
1813 11598 : F2x_edf_simple(ff, XP, d, V, idx+F2x_degree(f)/d);
1814 : }
1815 :
1816 : static GEN
1817 355533 : F2x_factor_Shoup(GEN T)
1818 : {
1819 355533 : long i, n, s = 0;
1820 : GEN XP, D, V;
1821 : pari_timer ti;
1822 355533 : n = F2x_degree(T);
1823 355533 : if (DEBUGLEVEL>=6) timer_start(&ti);
1824 355533 : XP = F2x_Frobenius(T);
1825 355530 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1826 355530 : D = F2x_ddf_simple(T, XP);
1827 355534 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1828 1107516 : for (i = 1; i <= n; i++)
1829 751983 : s += F2x_degree(gel(D,i))/i;
1830 355533 : V = cgetg(s+1, t_COL);
1831 1107454 : for (i = 1, s = 1; i <= n; i++)
1832 : {
1833 751909 : GEN Di = gel(D,i);
1834 751909 : long ni = F2x_degree(Di), ri = ni/i;
1835 751906 : if (ni == 0) continue;
1836 395066 : if (ni == i) { gel(V, s++) = Di; continue; }
1837 11293 : F2x_edf_simple(Di, XP, i, V, s);
1838 11311 : if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_edf(%ld)",i);
1839 11311 : s += ri;
1840 : }
1841 355545 : return V;
1842 : }
1843 :
1844 : static GEN
1845 304432 : F2x_factor_Cantor(GEN T)
1846 : {
1847 304432 : GEN E, F, V = F2x_factor_squarefree(T);
1848 304433 : long i, j, l = lg(V);
1849 304433 : E = cgetg(l, t_VEC);
1850 304433 : F = cgetg(l, t_VEC);
1851 1099283 : for (i=1, j=1; i < l; i++)
1852 794848 : if (F2x_degree(gel(V,i)))
1853 : {
1854 355533 : GEN Fj = F2x_factor_Shoup(gel(V,i));
1855 355532 : gel(F, j) = Fj;
1856 355532 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1857 355531 : j++;
1858 : }
1859 304435 : 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 304432 : F2x_Cantor_i(GEN f, long flag)
1927 : {
1928 : switch(flag)
1929 : {
1930 304432 : 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 1548378 : F2x_factor_i(GEN f, long flag)
1939 : {
1940 1548378 : long d = F2x_degree(f);
1941 1548367 : if (d <= 2) return F2x_factor_deg2(f,d,flag);
1942 315299 : return (flag == 0 && d <= 20)? F2x_Cantor_i(f, flag)
1943 655248 : : 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 gc_GEN(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 7318582 : Flx_factor_squarefree_pre(GEN f, ulong p, ulong pi)
1960 : {
1961 7318582 : long i, q, n = degpol(f);
1962 7318342 : GEN u = const_vec(n+1, pol1_Flx(f[1]));
1963 7319166 : for(q = 1;;q *= p)
1964 161094 : {
1965 7480260 : GEN t, v, tv, r = Flx_gcd_pre(f, Flx_deriv(f, p), p, pi);
1966 7476324 : if (degpol(r) == 0) { gel(u, q) = f; break; }
1967 665370 : t = Flx_div_pre(f, r, p, pi);
1968 665370 : if (degpol(t) > 0)
1969 : {
1970 : long j;
1971 511089 : for(j = 1;;j++)
1972 : {
1973 1266071 : v = Flx_gcd_pre(r, t, p, pi);
1974 1266066 : tv = Flx_div_pre(t, v, p, pi);
1975 1266071 : if (degpol(tv) > 0)
1976 760935 : gel(u, j*q) = Flx_normalize(tv, p);
1977 1266072 : if (degpol(v) <= 0) break;
1978 754975 : r = Flx_div_pre(r, v, p, pi);
1979 754982 : t = v;
1980 : }
1981 511094 : if (degpol(r) == 0) break;
1982 : }
1983 161094 : f = Flx_normalize(Flx_deflate(r, p), p);
1984 : }
1985 30789213 : for (i = n; i; i--)
1986 30787231 : if (degpol(gel(u,i))) break;
1987 7316876 : setlg(u,i+1); return u;
1988 : }
1989 : GEN
1990 349140 : Flx_factor_squarefree(GEN f, ulong p)
1991 349140 : { 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 = gc_leaf(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 9256651 : Flx_ddf_Shoup(GEN T, GEN XP, ulong p, ulong pi)
2024 : {
2025 9256651 : 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 9256651 : n = get_Flx_degree(T); v = get_Flx_var(T);
2031 9257854 : if (n == 0) return cgetg(1, t_VEC);
2032 9213655 : if (n == 1) return mkvec(get_Flx_mod(T));
2033 8825355 : B = n/2;
2034 8825355 : l = usqrt(B);
2035 8825969 : m = (B+l-1)/l;
2036 8825969 : T = Flx_get_red(T, p);
2037 8825539 : b = cgetg(l+2, t_VEC);
2038 8825333 : gel(b, 1) = polx_Flx(v);
2039 8826048 : gel(b, 2) = XP;
2040 8826048 : bo = brent_kung_optpow(n, l-1, 1);
2041 8826995 : ro = l<=1 ? 0:(bo-1)/(l-1) + ((n-1)/bo);
2042 8826995 : if (DEBUGLEVEL>=7) timer_start(&ti);
2043 8826995 : if (expu(p) <= ro)
2044 880555 : for (i = 3; i <= l+1; i++)
2045 488806 : gel(b, i) = Flxq_powu_pre(gel(b, i-1), p, T, p, pi);
2046 : else
2047 : {
2048 8435168 : xq = Flxq_powers_pre(gel(b, 2), bo, T, p, pi);
2049 8432805 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq baby");
2050 9159958 : for (i = 3; i <= l+1; i++)
2051 726916 : gel(b, i) = Flx_FlxqV_eval_pre(gel(b, i-1), xq, T, p, pi);
2052 : }
2053 8824791 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: baby");
2054 8824791 : xq = Flxq_powers_pre(gel(b, l+1), brent_kung_optpow(n, m-1, 1), T, p, pi);
2055 8824603 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq giant");
2056 8824603 : g = cgetg(m+1, t_VEC);
2057 8826085 : gel(g, 1) = gel(xq, 2);
2058 14309204 : for(i = 2; i <= m; i++)
2059 5482810 : gel(g, i) = Flx_FlxqV_eval_pre(gel(g, i-1), xq, T, p, pi);
2060 8826394 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: giant");
2061 8826394 : h = cgetg(m+1, t_VEC);
2062 23133594 : for (j = 1; j <= m; j++)
2063 : {
2064 14308490 : pari_sp av = avma;
2065 14308490 : GEN gj = gel(g, j);
2066 14308490 : GEN e = Flx_sub(gj, gel(b, 1), p);
2067 18232856 : for (i = 2; i <= l; i++)
2068 3927159 : e = Flxq_mul_pre(e, Flx_sub(gj, gel(b, i), p), T, p, pi);
2069 14305697 : gel(h, j) = gc_upto(av, e);
2070 : }
2071 8825104 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: diff");
2072 8825104 : Tr = get_Flx_mod(T);
2073 8825960 : F = cgetg(m+1, t_VEC);
2074 23134754 : for (j = 1; j <= m; j++)
2075 : {
2076 14307980 : GEN u = Flx_gcd_pre(Tr, gel(h, j), p, pi);
2077 14306433 : if (degpol(u))
2078 : {
2079 6619975 : u = Flx_normalize(u, p);
2080 6620834 : Tr = Flx_div_pre(Tr, u, p, pi);
2081 : }
2082 14307269 : gel(F, j) = u;
2083 : }
2084 8826774 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: F");
2085 8826774 : f = const_vec(n, pol1_Flx(v));
2086 23136002 : for (j = 1; j <= m; j++)
2087 : {
2088 14309897 : GEN e = gel(F, j);
2089 15194366 : for (i=l-1; i >= 0; i--)
2090 : {
2091 15194701 : GEN u = Flx_gcd_pre(e, Flx_sub(gel(g, j), gel(b, i+1), p), p, pi);
2092 15191495 : if (degpol(u))
2093 : {
2094 6720301 : gel(f, l*j-i) = u;
2095 6720301 : e = Flx_div_pre(e, u, p, pi);
2096 : }
2097 15190993 : if (!degpol(e)) break;
2098 : }
2099 : }
2100 8826105 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: f");
2101 8826105 : if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
2102 8826143 : return gc_GEN(av, f);
2103 : }
2104 :
2105 : static void
2106 434914 : Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, ulong pi, GEN V, long idx)
2107 : {
2108 434914 : long n = degpol(Tp), r = n/d;
2109 : GEN T, f, ff;
2110 : ulong p2;
2111 434914 : if (r==1) { gel(V, idx) = Tp; return; }
2112 187003 : p2 = p>>1;
2113 187003 : T = Flx_get_red_pre(Tp, p, pi);
2114 187003 : XP = Flx_rem_pre(XP, T, p, pi);
2115 : while (1)
2116 21493 : {
2117 208496 : pari_sp btop = avma;
2118 : long i;
2119 208496 : GEN g = random_Flx(n, Tp[1], p);
2120 208496 : GEN t = gel(Flxq_auttrace_pre(mkvec2(XP, g), d, T, p, pi), 2);
2121 208496 : if (lgpol(t) == 0) continue;
2122 458392 : for(i=1; i<=10; i++)
2123 : {
2124 442341 : pari_sp btop2 = avma;
2125 442341 : GEN R = Flxq_powu_pre(Flx_Fl_add(t, random_Fl(p), p), p2, T, p, pi);
2126 442341 : f = Flx_gcd_pre(Flx_Fl_add(R, p-1, p), Tp, p, pi);
2127 442338 : if (degpol(f) > 0 && degpol(f) < n) break;
2128 255338 : set_avma(btop2);
2129 : }
2130 203053 : if (degpol(f) > 0 && degpol(f) < n) break;
2131 16049 : set_avma(btop);
2132 : }
2133 187002 : f = Flx_normalize(f, p);
2134 187003 : ff = Flx_div_pre(Tp, f, p, pi);
2135 187003 : Flx_edf_simple(f, XP, d, p, pi, V, idx);
2136 187003 : 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 1408357 : 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 1408357 : GEN Tp = get_Flx_mod(T);
2147 1408361 : long n = degpol(hp), vT = Tp[1];
2148 : GEN u1, u2, f1, f2;
2149 1408363 : ulong p2 = p>>1;
2150 : GEN R, h;
2151 1408363 : h = Flx_get_red_pre(hp, p, pi);
2152 1408358 : t = Flx_rem_pre(t, T, p, pi);
2153 1408270 : av = avma;
2154 : do
2155 : {
2156 2369954 : set_avma(av);
2157 2369979 : R = Flxq_powu_pre(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p, pi);
2158 2369797 : u1 = Flx_gcd_pre(Flx_Fl_add(R, p-1, p), hp, p, pi);
2159 2369955 : } while (degpol(u1)==0 || degpol(u1)==n);
2160 1408288 : f1 = Flx_gcd_pre(Flx_Flxq_eval_pre(u1, t, T, p, pi), Tp, p, pi);
2161 1408285 : f1 = Flx_normalize(f1, p);
2162 1408352 : u2 = Flx_div_pre(hp, u1, p, pi);
2163 1408357 : f2 = Flx_div_pre(Tp, f1, p, pi);
2164 1408372 : if (degpol(u1)==1)
2165 : {
2166 1080624 : if (degpol(f1)==d)
2167 1065993 : gel(V, idx) = f1;
2168 : else
2169 14627 : Flx_edf(f1, XP, d, p, pi, V, idx);
2170 : }
2171 : else
2172 327777 : Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, pi, V, idx);
2173 1408405 : idx += degpol(f1)/d;
2174 1408390 : if (degpol(u2)==1)
2175 : {
2176 1074709 : if (degpol(f2)==d)
2177 1060165 : gel(V, idx) = f2;
2178 : else
2179 14541 : Flx_edf(f2, XP, d, p, pi, V, idx);
2180 : }
2181 : else
2182 333696 : Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, pi, V, idx);
2183 1408401 : }
2184 :
2185 : static void
2186 746940 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, ulong pi, GEN V, long idx)
2187 : {
2188 746940 : long n = degpol(Tp), r = n/d, vT = Tp[1];
2189 : GEN T, h, t;
2190 : pari_timer ti;
2191 746940 : if (r==1) { gel(V, idx) = Tp; return; }
2192 746940 : T = Flx_get_red_pre(Tp, p, pi);
2193 746938 : XP = Flx_rem_pre(XP, T, p, pi);
2194 746927 : if (DEBUGLEVEL>=7) timer_start(&ti);
2195 : do
2196 : {
2197 769509 : GEN g = random_Flx(n, vT, p);
2198 769525 : t = gel(Flxq_auttrace_pre(mkvec2(XP, g), d, T, p, pi), 2);
2199 769531 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
2200 769531 : h = Flxq_minpoly_pre(t, T, p, pi);
2201 769520 : if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
2202 769520 : } while (degpol(h) <= 1);
2203 746939 : Flx_edf_rec(T, XP, h, t, d, p, pi, V, idx);
2204 : }
2205 :
2206 : static GEN
2207 1690362 : Flx_factor_Shoup(GEN T, ulong p, ulong pi)
2208 : {
2209 1690362 : long i, n, s = 0, e = expu(p);
2210 : GEN XP, D, V;
2211 : pari_timer ti;
2212 1690360 : n = get_Flx_degree(T);
2213 1690356 : T = Flx_get_red_pre(T, p, pi);
2214 1690350 : if (DEBUGLEVEL>=6) timer_start(&ti);
2215 1690350 : XP = Flx_Frobenius_pre(T, p, pi);
2216 1690314 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2217 1690314 : D = Flx_ddf_Shoup(T, XP, p, pi);
2218 1690384 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2219 1690384 : s = ddf_to_nbfact(D);
2220 1690375 : V = cgetg(s+1, t_COL);
2221 9280321 : for (i = 1, s = 1; i <= n; i++)
2222 : {
2223 7589952 : GEN Di = gel(D,i);
2224 7589952 : long ni = degpol(Di), ri = ni/i;
2225 7589943 : if (ni == 0) continue;
2226 2188501 : Di = Flx_normalize(Di, p);
2227 2188521 : if (ni == i) { gel(V, s++) = Di; continue; }
2228 778674 : if (ri <= e*expu(e))
2229 717773 : Flx_edf(Di, XP, i, p, pi, V, s);
2230 : else
2231 60908 : Flx_edf_simple(Di, XP, i, p, pi, V, s);
2232 778682 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
2233 778682 : s += ri;
2234 : }
2235 1690369 : return V;
2236 : }
2237 :
2238 : static GEN
2239 1616992 : Flx_factor_Cantor(GEN T, ulong p)
2240 : {
2241 1616992 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2242 1616992 : GEN E, F, V = Flx_factor_squarefree_pre(get_Flx_mod(T), p, pi);
2243 1616992 : long i, j, l = lg(V);
2244 1616992 : F = cgetg(l, t_VEC);
2245 1616997 : E = cgetg(l, t_VEC);
2246 3860289 : for (i=1, j=1; i < l; i++)
2247 2243282 : if (degpol(gel(V,i)))
2248 : {
2249 1690362 : GEN Fj = Flx_factor_Shoup(gel(V,i), p, pi);
2250 1690370 : gel(F, j) = Fj;
2251 1690370 : gel(E, j) = const_vecsmall(lg(Fj)-1, i);
2252 1690371 : j++;
2253 : }
2254 1617007 : 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 5349817 : Flx_simplefact_Cantor(GEN T, ulong p)
2271 : {
2272 5349817 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2273 : long i, l;
2274 : GEN V;
2275 5349817 : T = Flx_get_red_pre(T, p, pi);
2276 5349366 : V = Flx_factor_squarefree_pre(get_Flx_mod(T), p, pi); l = lg(V);
2277 10782808 : for (i=1; i < l; i++)
2278 5431231 : gel(V,i) = Flx_ddf_Shoup(gel(V,i), Flx_Frobenius_pre(gel(V,i), p,pi), p,pi);
2279 5351577 : 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 11604926 : Flx_factor_i(GEN f, ulong pp, long flag)
2304 : {
2305 : long d;
2306 11604926 : if (pp==2) { /*We need to handle 2 specially */
2307 77845 : GEN F = F2x_factor_i(Flx_to_F2x(f),flag);
2308 77843 : if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
2309 77843 : return F;
2310 : }
2311 11527081 : d = degpol(f);
2312 11526901 : if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
2313 6967943 : switch(flag)
2314 : {
2315 1616992 : default: return Flx_factor_Cantor(f, pp);
2316 5349873 : 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 7862439 : Flx_degfact(GEN f, ulong p)
2323 : {
2324 7862439 : pari_sp av = avma;
2325 7862439 : GEN z = Flx_factor_i(Flx_normalize(f,p),p,1);
2326 7861389 : return gc_GEN(av, z);
2327 : }
2328 :
2329 : /* T must be squarefree mod p*/
2330 : GEN
2331 1456630 : Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
2332 : {
2333 : GEN XP, D;
2334 : pari_timer ti;
2335 1456630 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2336 1456630 : long i, s, n = get_Flx_degree(T);
2337 1456605 : GEN V = const_vecsmall(n, 0);
2338 1456599 : pari_sp av = avma;
2339 1456599 : T = Flx_get_red_pre(T, p, pi);
2340 1456616 : if (DEBUGLEVEL>=6) timer_start(&ti);
2341 1456616 : XP = Flx_Frobenius_pre(T, p, pi);
2342 1456532 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2343 1456532 : D = Flx_ddf_Shoup(T, XP, p, pi);
2344 1456838 : if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2345 7544288 : for (i = 1, s = 0; i <= n; i++) { V[i] = degpol(gel(D,i))/i; s += V[i]; }
2346 1456791 : *nb = s; set_avma(av); return V;
2347 : }
2348 :
2349 : long
2350 680323 : Flx_nbfact_Frobenius_pre(GEN T, GEN XP, ulong p, ulong pi)
2351 : {
2352 680323 : pari_sp av = avma;
2353 680323 : long s = ddf_to_nbfact(Flx_ddf_Shoup(T, XP, p, pi));
2354 680323 : 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 680325 : Flx_nbfact_pre(GEN T, ulong p, ulong pi)
2363 : {
2364 680325 : pari_sp av = avma;
2365 680325 : GEN XP = Flx_Frobenius_pre(T, p, pi);
2366 680323 : long n = Flx_nbfact_Frobenius_pre(T, XP, p, pi);
2367 680319 : return gc_long(av,n);
2368 : }
2369 : long
2370 680325 : Flx_nbfact(GEN T, ulong p)
2371 680325 : { 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 gc_GEN(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 4431099 : FpX_factor(GEN f, GEN p)
2431 : {
2432 4431099 : pari_sp av = avma;
2433 : GEN F;
2434 4431099 : switch(ZX_factmod_init(&f, p))
2435 : {
2436 1453096 : case 0: F = F2x_factor_i(f,0);
2437 1453093 : F2xV_to_ZXV_inplace(gel(F,1)); break;
2438 2975697 : case 1: F = Flx_factor_i(f,p[2],0);
2439 2975721 : FlxV_to_ZXV_inplace(gel(F,1)); break;
2440 2365 : default: F = FpX_factor_i(f,p,0); break;
2441 : }
2442 4431188 : return gc_GEN(av, F);
2443 : }
2444 :
2445 : GEN
2446 718736 : Flx_factor(GEN f, ulong p)
2447 : {
2448 718736 : pari_sp av = avma;
2449 718736 : return gc_GEN(av, Flx_factor_i(Flx_normalize(f,p),p,0));
2450 : }
2451 : GEN
2452 16697 : F2x_factor(GEN f)
2453 : {
2454 16697 : pari_sp av = avma;
2455 16697 : return gc_GEN(av, F2x_factor_i(f,0));
2456 : }
|
