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