1 /* $Id: polarit1.c 11851 2009-07-21 21:52:51Z bill $
2
3 Copyright (C) 2000-2004 The PARI group.
4
5 This file is part of the PARI/GP package.
6
7 PARI/GP is free software; you can redistribute it and/or modify it under the
8 terms of the GNU General Public License as published by the Free Software
9 Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10 ANY WARRANTY WHATSOEVER.
11
12 Check the License for details. You should have received a copy of it, along
13 with the package; see the file 'COPYING'. If not, write to the Free Software
14 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15
16 /***********************************************************************/
17 /** **/
18 /** ARITHMETIC OPERATIONS ON POLYNOMIALS **/
19 /** (first part) **/
20 /** **/
21 /***********************************************************************/
22 #include "pari.h"
23 #include "paripriv.h"
24
25 /*******************************************************************/
26 /* */
27 /* DIVISIBILITY */
28 /* Return 1 if y | x, 0 otherwise */
29 /* */
30 /*******************************************************************/
31
32 int
gdvd(GEN x,GEN y)33 gdvd(GEN x, GEN y)
34 {
35 pari_sp av=avma;
36 x=gmod(x,y); avma=av; return gcmp0(x);
37 }
38
39 int
poldvd(GEN x,GEN y,GEN * z)40 poldvd(GEN x, GEN y, GEN *z)
41 {
42 pari_sp av = avma;
43 GEN p1 = poldivrem(x,y,ONLY_DIVIDES);
44 if (p1) { *z = p1; return 1; }
45 avma=av; return 0;
46 }
47
48 /*******************************************************************/
49 /* */
50 /* POLYNOMIAL EUCLIDEAN DIVISION */
51 /* */
52 /*******************************************************************/
53 GEN
poldivrem(GEN x,GEN y,GEN * pr)54 poldivrem(GEN x, GEN y, GEN *pr)
55 {
56 long ty = typ(y), tx, vx = gvar(x), vy = gvar(y);
57 GEN p1;
58
59 if (is_scalar_t(ty) || varncmp(vx, vy) < 0)
60 {
61 if (pr == ONLY_REM) {
62 if (gcmp0(y)) pari_err(gdiver);
63 return gen_0;
64 }
65 if (pr && pr != ONLY_DIVIDES) *pr=gen_0;
66 return gdiv(x,y);
67 }
68 if (ty != t_POL) pari_err(typeer,"euclidean division (poldivrem)");
69 tx = typ(x);
70 if (is_scalar_t(tx) || varncmp(vx, vy) > 0)
71 {
72 if (!signe(y)) pari_err(gdiver);
73 if (!degpol(y)) /* constant */
74 {
75 if (pr == ONLY_REM) return zeropol(vy);
76 p1 = gdiv(x, gel(y,2));
77 if (pr == ONLY_DIVIDES) return p1;
78 if (pr) *pr = zeropol(vy);
79 return p1;
80 }
81 if (pr == ONLY_REM) return gcopy(x);
82 if (pr == ONLY_DIVIDES) return gcmp0(x)? gen_0: NULL;
83 if (pr) *pr = gcopy(x);
84 return gen_0;
85 }
86 if (tx != t_POL) pari_err(typeer,"euclidean division (poldivrem)");
87
88 if (varncmp(vx, vy) < 0)
89 {
90 if (pr && pr != ONLY_DIVIDES)
91 {
92 p1 = zeropol(vx); if (pr == ONLY_REM) return p1;
93 *pr = p1;
94 }
95 return gdiv(x,y);
96 }
97 return RgX_divrem(x, y, pr);
98 }
99
100 /*******************************************************************/
101 /* */
102 /* ROOTS MODULO a prime p (no multiplicities) */
103 /* */
104 /*******************************************************************/
105 static ulong
init_p(GEN pp)106 init_p(GEN pp)
107 {
108 ulong p;
109 if ((ulong)expi(pp) > BITS_IN_LONG - 3) p = 0;
110 else
111 {
112 p = itou(pp);
113 if (p < 2 || signe(pp) < 0) pari_err(talker,"not a prime in factmod");
114 }
115 return p;
116 }
117
118 static long
factmod_init(GEN * F,GEN p)119 factmod_init(GEN *F, GEN p)
120 {
121 long d;
122 if (typ(*F)!=t_POL || typ(p)!=t_INT) pari_err(typeer,"factmod");
123 *F = FpX_normalize(RgX_to_FpX(*F, p), p);
124 d = degpol(*F); if (d < 0) pari_err(zeropoler,"factmod");
125 return d;
126 }
127
128 static GEN
root_mod_2(GEN f)129 root_mod_2(GEN f)
130 {
131 int z1, z0 = !signe(constant_term(f));
132 long i,n;
133 GEN y;
134
135 for (i=2, n=1; i < lg(f); i++)
136 if (signe(f[i])) n++;
137 z1 = n & 1;
138 y = cgetg(z0+z1+1, t_COL); i = 1;
139 if (z0) gel(y,i++) = gen_0;
140 if (z1) gel(y,i) = gen_1;
141 return y;
142 }
143
144 #define i_mod4(x) (signe(x)? mod4((GEN)(x)): 0)
145 static GEN
root_mod_4(GEN f)146 root_mod_4(GEN f)
147 {
148 long i, no, ne;
149 GEN y, t;
150 int z0, z1, z2, z3;
151
152 t = constant_term(f);
153 z0 = !signe(t);
154 z2 = ((i_mod4(t) + 2*i_mod4(f[3])) & 3) == 0;
155
156 for (ne=0,i=2; i<lg(f); i+=2)
157 {
158 t = gel(f,i);
159 if (signe(t)) ne += *int_LSW(t);
160 }
161 for (no=0,i=3; i<lg(f); i+=2)
162 {
163 t = gel(f,i);
164 if (signe(t)) ne += *int_LSW(t);
165 }
166 no &= 3; ne &= 3;
167 z3 = (no == ne);
168 z1 = (no == ((4-ne)&3));
169 y=cgetg(1+z0+z1+z2+z3,t_COL); i = 1;
170 if (z0) gel(y,i++) = gen_0;
171 if (z1) gel(y,i++) = gen_1;
172 if (z2) gel(y,i++) = gen_2;
173 if (z3) gel(y,i) = utoipos(3);
174 return y;
175 }
176 #undef i_mod4
177
178 /* p even, accept p = 4 for p-adic stuff */
179 INLINE GEN
root_mod_even(GEN f,ulong p)180 root_mod_even(GEN f, ulong p)
181 {
182 switch(p)
183 {
184 case 2: return root_mod_2(f);
185 case 4: return root_mod_4(f);
186 }
187 pari_err(talker,"not a prime in rootmod");
188 return NULL; /* not reached */
189 }
190
191 /* by checking f(0..p-1) */
192 static GEN
Flx_roots_naive(GEN f,ulong p)193 Flx_roots_naive(GEN f, ulong p)
194 {
195 long d = degpol(f), n = 0;
196 ulong s = 1UL, r;
197 GEN q, y = cgetg(d + 1, t_VECSMALL);
198 pari_sp av = avma;
199
200 if (!f[2]) y[++n] = 0;
201 do
202 {
203 q = Flx_div_by_X_x(f, s, p, &r); /* TODO: FFT-type multi-evaluation */
204 if (r) avma = av;
205 else
206 {
207 y[++n] = s;
208 f = q; av = avma;
209 }
210 s++;
211 }
212 while (n < d-1 && p > s);
213 if (n == d-1 && p != s) /* -f[2]/f[3] */
214 y[++n] = Fl_mul(p - Fl_inv((ulong)f[3], p), (ulong)f[2], p);
215 setlg(y, n+1); return y;
216 }
217
218 GEN
rootmod2(GEN f,GEN pp)219 rootmod2(GEN f, GEN pp)
220 {
221 pari_sp av = avma;
222 ulong p;
223 GEN y;
224
225 if (!factmod_init(&f, pp)) { avma = av; return cgetg(1,t_COL); }
226 p = init_p(pp); if (!p) pari_err(talker,"prime too big in rootmod2");
227 if (p & 1)
228 y = Flc_to_ZC(Flx_roots_naive(ZX_to_Flx(f,p), p));
229 else
230 y = root_mod_even(f,p);
231 return gerepileupto(av, FpC_to_mod(y, pp));
232 }
233
234 /* assume x reduced mod p, monic. */
235 static int
FpX_quad_factortype(GEN x,GEN p)236 FpX_quad_factortype(GEN x, GEN p)
237 {
238 GEN b = gel(x,3), c = gel(x,2);
239 GEN D;
240
241 if (equaliu(p, 2)) {
242 if (!signe(b)) return 0;
243 return signe(c)? -1: 1;
244 }
245 D = subii(sqri(b), shifti(c,2));
246 return kronecker(D,p);
247 }
248 /* assume x reduced mod p, monic. Return one root, or NULL if irreducible */
249 GEN
FpX_quad_root(GEN x,GEN p,int unknown)250 FpX_quad_root(GEN x, GEN p, int unknown)
251 {
252 GEN s, u, D, b = gel(x,3), c = gel(x,2);
253
254 if (equaliu(p, 2)) {
255 if (!signe(b)) return c;
256 return signe(c)? NULL: gen_1;
257 }
258 D = subii(sqri(b), shifti(c,2));
259 D = remii(D,p);
260 if (unknown && kronecker(D,p) == -1) return NULL;
261
262 s = Fp_sqrt(D,p);
263 /* p is not prime, go on and give e.g. maxord a chance to recover */
264 if (!s) return NULL;
265 u = addis(shifti(p,-1), 1); /* = 1/2 */
266 return modii(mulii(u, subii(s,b)), p);
267 }
268 static GEN
otherroot(GEN x,GEN r,GEN p)269 otherroot(GEN x, GEN r, GEN p)
270 {
271 GEN s = addii(gel(x,3), r);
272 if (!signe(s)) return s;
273 s = subii(p, s); if (signe(s) < 0) s = addii(s,p);
274 return s;
275 }
276
277 /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
278 static GEN
FpX_roots_i(GEN f,GEN p)279 FpX_roots_i(GEN f, GEN p)
280 {
281 long n, j, da, db;
282 GEN y, pol, pol0, a, b, q = shifti(p,-1);
283
284 n = ZX_valuation(f, &f)? 1: 0;
285 y = cgetg(lg(f), t_COL);
286 j = 1;
287 if (n) {
288 gel(y, j++) = gen_0;
289 if (lg(f) <= 3) { setlg(y, 2); return y; }
290 n = 1;
291 }
292 da = degpol(f);
293 if (da == 1) { gel(y,j++) = subii(p, gel(f,2)); setlg(y,j); return y; }
294 if (da == 2) {
295 GEN s, r = FpX_quad_root(f, p, 1);
296 if (r) {
297 gel(y, j++) = r;
298 s = otherroot(f,r, p);
299 if (!equalii(r, s)) gel(y, j++) = s;
300 }
301 setlg(y, j); return sort(y);
302 }
303
304 /* take gcd(x^(p-1) - 1, f) by splitting (x^q-1) * (x^q+1) */
305 b = FpXQ_pow(pol_x[varn(f)],q, f,p);
306 if (lg(b) < 3) pari_err(talker,"not a prime in rootmod");
307 b = ZX_Z_add(b, gen_m1); /* b = x^((p-1)/2) - 1 mod f */
308 a = FpX_gcd(f,b, p);
309 b = ZX_Z_add(b, gen_2); /* b = x^((p-1)/2) + 1 mod f */
310 b = FpX_gcd(f,b, p);
311 da = degpol(a);
312 db = degpol(b); n += da + db; setlg(y, n+1);
313 if (db) gel(y,j) = FpX_normalize(b,p);
314 if (da) gel(y,j+db) = FpX_normalize(a,p);
315 pol = gadd(pol_x[varn(f)], gen_1); pol0 = constant_term(pol);
316 while (j <= n)
317 { /* cf FpX_split_Berlekamp */
318 a = gel(y,j); da = degpol(a);
319 if (da==1)
320 gel(y,j++) = subii(p, gel(a,2));
321 else if (da==2)
322 {
323 GEN r = FpX_quad_root(a, p, 0);
324 gel(y, j++) = r;
325 gel(y, j++) = otherroot(a,r, p);
326 }
327 else for (pol0[2]=1; ; pol0[2]++)
328 {
329 b = ZX_Z_add(FpXQ_pow(pol,q, a,p), gen_m1); /* pol^(p-1)/2 - 1 */
330 b = FpX_gcd(a,b, p); db = degpol(b);
331 if (db && db < da)
332 {
333 b = FpX_normalize(b, p);
334 gel(y,j+db) = FpX_div(a,b, p);
335 gel(y,j) = b; break;
336 }
337 if (pol0[2] == 100 && !BSW_psp(p))
338 pari_err(talker, "not a prime in polrootsmod");
339 }
340 }
341 return sort(y);
342 }
343
344 static GEN
FpX_factmod_init(GEN f,GEN p)345 FpX_factmod_init(GEN f, GEN p) { return FpX_normalize(FpX_red(f,p), p); }
346 GEN
FpX_roots(GEN f,GEN p)347 FpX_roots(GEN f, GEN p) {
348 pari_sp av = avma;
349 long q = modBIL(p);
350 GEN F = FpX_factmod_init(f,p);
351 switch(degpol(F))
352 {
353 case -1: pari_err(zeropoler,"factmod");
354 case 0: avma = av; return cgetg(1, t_COL);
355 }
356 return gerepileupto(av, odd(q)? FpX_roots_i(F, p): root_mod_even(F,q));
357 }
358
359 GEN
rootmod(GEN f,GEN p)360 rootmod(GEN f, GEN p)
361 {
362 ulong q;
363 pari_sp av = avma;
364 GEN y;
365
366 if (!factmod_init(&f, p)) { avma=av; return cgetg(1,t_COL); }
367 q = init_p(p); if (!q) q = modBIL(p);
368 if (q & 1)
369 y = FpX_roots_i(f, p);
370 else
371 y = root_mod_even(f,q);
372 return gerepileupto(av, FpC_to_mod(y, p));
373 }
374
375 GEN
rootmod0(GEN f,GEN p,long flag)376 rootmod0(GEN f, GEN p, long flag)
377 {
378 switch(flag)
379 {
380 case 0: return rootmod(f,p);
381 case 1: return rootmod2(f,p);
382 default: pari_err(flagerr,"polrootsmod");
383 }
384 return NULL; /* not reached */
385 }
386
387 /*******************************************************************/
388 /* */
389 /* FACTORISATION MODULO p */
390 /* */
391 /*******************************************************************/
392 static GEN spec_FpXQ_pow(GEN x, GEN p, GEN S);
393
394 /* Functions giving information on the factorisation. */
395
396 /*TODO: merge with FpXQ_powers */
397 /* u in Z[X], return kernel of (Frob - Id) over Fp[X] / u */
398 GEN
FpX_Berlekamp_ker(GEN u,GEN p)399 FpX_Berlekamp_ker(GEN u, GEN p)
400 {
401 long j,N = degpol(u);
402 GEN v,w,Q,p1;
403 Q = cgetg(N+1,t_MAT); gel(Q,1) = zerocol(N);
404 w = v = FpXQ_pow(pol_x[varn(u)],p,u,p);
405 for (j=2; j<=N; j++)
406 {
407 p1 = RgX_to_RgV(w, N);
408 gel(p1,j) = addis(gel(p1,j), -1);
409 gel(Q,j) = p1;
410 if (j < N)
411 {
412 pari_sp av = avma; /*FpXQ_mul is not stack clean*/
413 w = gerepileupto(av, FpXQ_mul(w, v, u, p));
414 }
415 }
416 return FpM_ker(Q,p);
417 }
418
419 GEN
FqX_Berlekamp_ker(GEN u,GEN T,GEN q,GEN p)420 FqX_Berlekamp_ker(GEN u, GEN T, GEN q, GEN p)
421 {
422 pari_sp ltop=avma;
423 long j,N = degpol(u);
424 GEN v,w,Q,p1;
425 pari_timer Ti;
426 if (DEBUGLEVEL>=4) TIMERstart(&Ti);
427 Q = cgetg(N+1,t_MAT); gel(Q,1) = zerocol(N);
428 w = v = FpXQYQ_pow(pol_x[varn(u)], q, u, T, p);
429 if (DEBUGLEVEL>=4) msgTIMER(&Ti, "FpXQYQ_pow");
430 for (j=2; j<=N; j++)
431 {
432 p1 = RgX_to_RgV(w, N);
433 gel(p1,j) = gaddgs(gel(p1,j), -1);
434 gel(Q,j) = p1;
435 if (j < N)
436 {
437 pari_sp av = avma;
438 w = gerepileupto(av, FpXQX_divrem(FpXQX_mul(w,v, T,p), u,T,p,ONLY_REM));
439 }
440 }
441 if (DEBUGLEVEL>=4) msgTIMER(&Ti, "Berlekamp_matrix");
442 p1 = FqM_ker(Q,T,p);
443 if (DEBUGLEVEL>=4) msgTIMER(&Ti, "Berlekamp_ker");
444 return gerepileupto(ltop,p1);
445 }
446
447 GEN
Flx_Berlekamp_ker(GEN u,ulong p)448 Flx_Berlekamp_ker(GEN u, ulong p)
449 {
450 pari_sp ltop=avma;
451 long j,N = degpol(u);
452 GEN v,w,Q,p1;
453 pari_timer T;
454 TIMERstart(&T);
455 Q = cgetg(N+1,t_VEC); gel(Q,1) = const_vecsmall(N,0);
456 w = v = Flxq_pow(polx_Flx(u[1]), utoipos(p), u, p);
457 for (j=2; j<=N; j++)
458 {
459 p1 = Flx_to_Flv(w, N);
460 p1[j] = Fl_sub((ulong)p1[j], 1, p);
461 gel(Q,j) = p1;
462 if (j < N)
463 {
464 pari_sp av = avma; /*Flxq_mul is not stack clean*/
465 w = gerepileupto(av, Flxq_mul(w, v, u, p));
466 }
467 }
468 if(DEBUGLEVEL>=9) msgTIMER(&T,"Berlekamp matrix");
469 Q=Flm_ker_sp(Q,p,0);
470 if(DEBUGLEVEL>=9) msgTIMER(&T,"kernel");
471 return gerepileupto(ltop,Q);
472 }
473
474 GEN
ZX_deriv(GEN x)475 ZX_deriv(GEN x)
476 {
477 long i,lx = lg(x)-1;
478 GEN y;
479
480 if (lx<3) return zeropol(varn(x));
481 y = cgetg(lx,t_POL);
482 for (i=2; i<lx ; i++) gel(y,i) = mulsi(i-1,gel(x,i+1));
483 y[1] = x[1]; return y;
484 }
485
486 GEN
FpX_deriv(GEN f,GEN p)487 FpX_deriv(GEN f, GEN p) { return FpX_red(ZX_deriv(f), p); }
488 GEN
FqX_deriv(GEN f,GEN T,GEN p)489 FqX_deriv(GEN f, /*unused*/GEN T, GEN p) {
490 (void)T; return FpXX_red(derivpol(f), p);
491 }
492
493 /* f in ZZ[X] and p a prime number. */
494 long
FpX_is_squarefree(GEN f,GEN p)495 FpX_is_squarefree(GEN f, GEN p)
496 {
497 pari_sp av = avma;
498 GEN z;
499 z = FpX_gcd(f,ZX_deriv(f),p);
500 avma = av;
501 return lg(z)==3;
502 }
503
504 /* product of terms of degree 1 in factorization of f */
505 static GEN
FpX_split_part(GEN f,GEN p)506 FpX_split_part(GEN f, GEN p)
507 {
508 long n = degpol(f);
509 GEN z, X = pol_x[varn(f)];
510 if (n <= 1) return f;
511 f = FpX_red(f, p);
512 z = FpXQ_pow(X, p, f, p);
513 z = FpX_sub(z, X, p);
514 return FpX_gcd(z,f,p);
515 }
516
517 static GEN
FqX_split_part(GEN f,GEN T,GEN p)518 FqX_split_part(GEN f, GEN T, GEN p)
519 {
520 long n = degpol(f);
521 GEN z, X = pol_x[varn(f)];
522 if (n <= 1) return f;
523 f = FpXQX_red(f, T, p);
524 z = FpXQYQ_pow(X, powiu(p, degpol(T)), f, T, p);
525 z = gsub(z, X);
526 return FqX_gcd(z, f, T, p);
527 }
528
529 /* Compute the number of roots in Fq without counting multiplicity
530 * return -1 for 0 polynomial. lc(f) must be prime to p. */
531 long
FpX_nbroots(GEN f,GEN p)532 FpX_nbroots(GEN f, GEN p)
533 {
534 pari_sp av = avma;
535 GEN z = FpX_split_part(f, p);
536 avma = av; return degpol(z);
537 }
538
539 long
FqX_nbroots(GEN f,GEN T,GEN p)540 FqX_nbroots(GEN f, GEN T, GEN p)
541 {
542 pari_sp av = avma;
543 GEN z = FqX_split_part(f, T, p);
544 avma = av; return degpol(z);
545 }
546
547 long
FpX_is_totally_split(GEN f,GEN p)548 FpX_is_totally_split(GEN f, GEN p)
549 {
550 long n=degpol(f);
551 pari_sp av = avma;
552 GEN z;
553 if (n <= 1) return 1;
554 if (cmpui(n, p) > 0) return 0;
555 f = FpX_red(f, p);
556 z = FpXQ_pow(pol_x[varn(f)], p, f, p);
557 avma = av;
558 return degpol(z) == 1 && gcmp1(gel(z,3)) && !signe(z[2]); /* x^p = x ? */
559 }
560
561 /* Flv_Flx( Flm_Flc_mul(x, Flx_Flv(y), p) ) */
562 static GEN
Flm_Flx_mul(GEN x,GEN y,ulong p)563 Flm_Flx_mul(GEN x, GEN y, ulong p)
564 {
565 long i,k,l, ly = lg(y)-1;
566 GEN z;
567 long vs=y[1];
568 if (ly==1) return zero_Flx(vs);
569 l = lg(x[1]);
570 y++;
571 z = const_vecsmall(l,0) + 1;
572 if (u_OK_ULONG(p))
573 {
574 for (k=1; k<ly; k++)
575 {
576 GEN c;
577 if (!y[k]) continue;
578 c = gel(x,k);
579 if (y[k] == 1)
580 for (i=1; i<l; i++)
581 {
582 z[i] += c[i];
583 if (z[i] & HIGHBIT) z[i] %= p;
584 }
585 else
586 for (i=1; i<l; i++)
587 {
588 z[i] += c[i] * y[k];
589 if (z[i] & HIGHBIT) z[i] %= p;
590 }
591 }
592 for (i=1; i<l; i++) z[i] %= p;
593 }
594 else
595 {
596 for (k=1; k<ly; k++)
597 {
598 GEN c;
599 if (!y[k]) continue;
600 c = gel(x,k);
601 if (y[k] == 1)
602 for (i=1; i<l; i++)
603 z[i] = Fl_add(z[i], c[i], p);
604 else
605 for (i=1; i<l; i++)
606 z[i] = Fl_add(z[i], Fl_mul(c[i],y[k],p), p);
607 }
608 }
609 while (--l && !z[l]);
610 if (!l) return zero_Flx(vs);
611 *z-- = vs; return z;
612 }
613
614 /* assume deg(u) > 0 */
615 static GEN
Flx_Frobenius(GEN u,ulong p)616 Flx_Frobenius(GEN u, ulong p)
617 {
618 long j,N = degpol(u);
619 GEN v,w,Q;
620 pari_timer T;
621
622 if (DEBUGLEVEL > 7) TIMERstart(&T);
623 Q = cgetg(N+1,t_MAT);
624 gel(Q,1) = const_vecsmall(N, 0);
625 coeff(Q,1,1) = 1;
626 w = v = Flxq_pow(polx_Flx(u[1]), utoipos(p), u, p);
627 for (j=2; j<=N; j++)
628 {
629 gel(Q,j) = Flx_to_Flv(w, N);
630 if (j < N)
631 {
632 pari_sp av = avma;
633 w = gerepileupto(av, Flxq_mul(w, v, u, p));
634 }
635 }
636 if (DEBUGLEVEL > 7) msgTIMER(&T, "frobenius");
637 return Q;
638 }
639
640 /* z must be squarefree mod p*/
641 long
Flx_nbfact(GEN z,ulong p)642 Flx_nbfact(GEN z, ulong p)
643 {
644 long lgg, nfacp = 0, d = 0, e = degpol(z);
645 GEN g, w, MP = Flx_Frobenius(z, p), PolX = polx_Flx(z[1]);
646
647 w = PolX;
648 while (d < (e>>1))
649 { /* here e = degpol(z) */
650 d++;
651 w = Flm_Flx_mul(MP, w, p); /* w^p mod (z,p) */
652 g = Flx_gcd(z, Flx_sub(w, PolX, p), p);
653 lgg = degpol(g);
654 if (!lgg) continue;
655
656 e -= lgg;
657 nfacp += lgg/d;
658 if (DEBUGLEVEL>5)
659 fprintferr(" %3ld fact. of degree %3ld\n", lgg/d, d);
660 if (!e) break;
661 z = Flx_div(z, g, p);
662 w = Flx_rem(w, z, p);
663 }
664 if (e)
665 {
666 if (DEBUGLEVEL>5) fprintferr(" %3ld factor of degree %3ld\n",1,e);
667 nfacp++;
668 }
669 return nfacp;
670 }
671
672 long
Flx_nbroots(GEN f,ulong p)673 Flx_nbroots(GEN f, ulong p)
674 {
675 long n = degpol(f);
676 pari_sp av = avma;
677 GEN z, X;
678 if (n <= 1) return n;
679 X = polx_Flx(f[1]);
680 z = Flxq_pow(X, utoipos(p), f, p);
681 z = Flx_sub(z, X, p);
682 z = Flx_gcd(z, f, p);
683 avma = av; return degpol(z);
684 }
685
686 long
FpX_nbfact(GEN u,GEN p)687 FpX_nbfact(GEN u, GEN p)
688 {
689 pari_sp av = avma;
690 GEN vker = FpX_Berlekamp_ker(u, p);
691 avma = av; return lg(vker)-1;
692 }
693
694 long
FqX_nbfact(GEN u,GEN T,GEN p)695 FqX_nbfact(GEN u, GEN T, GEN p)
696 {
697 pari_sp av = avma;
698 GEN vker;
699 if (!T) return FpX_nbfact(u, p);
700 vker = FqX_Berlekamp_ker(u, T, powiu(p, degpol(T)), p);
701 avma = av; return lg(vker)-1;
702 }
703
704 /************************************************************/
705 GEN
trivfact(void)706 trivfact(void)
707 {
708 GEN y = cgetg(3,t_MAT);
709 gel(y,1) = cgetg(1,t_COL);
710 gel(y,2) = cgetg(1,t_COL); return y;
711 }
712
713 static GEN
try_pow(GEN w0,GEN pol,GEN p,GEN q,long r)714 try_pow(GEN w0, GEN pol, GEN p, GEN q, long r)
715 {
716 GEN w2, w = FpXQ_pow(w0,q, pol,p);
717 long s;
718 if (gcmp1(w)) return w0;
719 for (s=1; s<r; s++,w=w2)
720 {
721 w2 = gsqr(w);
722 w2 = FpX_rem(w2, pol, p);
723 if (gcmp1(w2)) break;
724 }
725 return gcmp_1(w)? NULL: w;
726 }
727
728 /* INPUT:
729 * m integer (converted to polynomial w in Z[X] by stopoly)
730 * p prime; q = (p^d-1) / 2^r
731 * t[0] polynomial of degree k*d product of k irreducible factors of degree d
732 * t[0] is expected to be normalized (leading coeff = 1)
733 * OUTPUT:
734 * t[0],t[1]...t[k-1] the k factors, normalized */
735 static void
split(ulong m,GEN * t,long d,GEN p,GEN q,long r,GEN S)736 split(ulong m, GEN *t, long d, GEN p, GEN q, long r, GEN S)
737 {
738 long l, v, dv;
739 ulong ps;
740 pari_sp av0, av;
741 GEN w,w0;
742
743 dv=degpol(*t); if (dv==d) return;
744 v=varn(*t); av0=avma; ps = (ulong)p[2];
745 for(av=avma;;avma=av)
746 {
747 if (ps==2)
748 {
749 w0 = w = FpXQ_pow(pol_x[v], utoi(m-1), *t, gen_2); m += 2;
750 for (l=1; l<d; l++)
751 w = gadd(w0, spec_FpXQ_pow(w, p, S));
752 }
753 else
754 {
755 w = FpX_rem(stopoly(m,ps,v),*t, p);
756 m++; w = try_pow(w,*t,p,q,r);
757 if (!w) continue;
758 w = ZX_Z_add(w, gen_m1);
759 }
760 w = FpX_gcd(*t,w, p);
761 l = degpol(w); if (l && l!=dv) break;
762 }
763 w = FpX_normalize(w, p);
764 w = gerepileupto(av0, w);
765 l /= d; t[l]=FpX_div(*t,w,p); *t=w;
766 split(m,t+l,d,p,q,r,S);
767 split(m,t, d,p,q,r,S);
768 }
769
770 static void
splitgen(GEN m,GEN * t,long d,GEN p,GEN q,long r)771 splitgen(GEN m, GEN *t, long d, GEN p, GEN q, long r)
772 {
773 long l, v, dv = degpol(*t);
774 pari_sp av;
775 GEN w;
776
777 if (dv==d) return;
778 v = varn(*t);
779 m = setloop(m);
780 av = avma;
781 for(;; avma = av)
782 {
783 m = incloop(m);
784 w = FpX_rem(stopoly_gen(m,p,v),*t, p);
785 w = try_pow(w,*t,p,q,r);
786 if (!w) continue;
787 w = ZX_Z_add(w, gen_m1);
788 w = FpX_gcd(*t,w, p); l=degpol(w);
789 if (l && l!=dv) break;
790 }
791 w = FpX_normalize(w, p);
792 w = gerepileupto(av, w);
793 l /= d; t[l]=FpX_div(*t,w,p); *t=w;
794 splitgen(m,t+l,d,p,q,r);
795 splitgen(m,t, d,p,q,r);
796 }
797
798 /* return S = [ x^p, x^2p, ... x^(n-1)p ] mod (p, T), n = degree(T) > 0 */
799 static GEN
init_spec_FpXQ_pow(GEN p,GEN T)800 init_spec_FpXQ_pow(GEN p, GEN T)
801 {
802 long i, n = degpol(T), v = varn(T);
803 GEN S = cgetg(n, t_VEC), x;
804 if (n == 1) return S;
805 x = FpXQ_pow(pol_x[v], p, T, p);
806 gel(S,1) = x;
807 if ((degpol(x)<<1) < degpol(T)) {
808 for (i=2; i < n; i++)
809 gel(S,i) = FpXQ_mul(gel(S,i-1), x, T,p);
810 } else {
811 for (i=2; i < n; i++)
812 gel(S,i) = (i&1)? FpXQ_mul(gel(S,i-1), x, T,p)
813 : FpXQ_sqr(gel(S,i>>1), T, p);
814 }
815 return S;
816 }
817
818 /* compute x^p, S is as above */
819 static GEN
spec_FpXQ_pow(GEN x,GEN p,GEN S)820 spec_FpXQ_pow(GEN x, GEN p, GEN S)
821 {
822 long i, dx = degpol(x);
823 pari_sp av = avma, lim = stack_lim(av, 1);
824 GEN x0 = x+2, z = gel(x0,0);
825 if (dx < 0) pari_err(talker, "zero polynomial in FpXQ_pow. %Z not prime", p);
826 for (i = 1; i <= dx; i++)
827 {
828 GEN d, c = gel(x0,i); /* assume coeffs in [0, p-1] */
829 if (!signe(c)) continue;
830 d = gel(S,i); if (!gcmp1(c)) d = gmul(c,d);
831 z = gadd(z, d);
832 if (low_stack(lim, stack_lim(av,1)))
833 {
834 if(DEBUGMEM>1) pari_warn(warnmem,"spec_FpXQ_pow");
835 z = gerepileupto(av, z);
836 }
837 }
838 return gerepileupto(av, FpX_red(z, p));
839 }
840
841 static int
cmpGsGs(GEN a,GEN b)842 cmpGsGs(GEN a, GEN b) { return (long)a - (long)b; }
843
844 static GEN
FpX_is_irred_2(GEN f,GEN p,long d)845 FpX_is_irred_2(GEN f, GEN p, long d)
846 {
847 if (!d) return NULL;
848 if (d == 1) return gen_1;
849 return FpX_quad_factortype(f, p) == -1? gen_1: NULL;
850 }
851 static GEN
FpX_degfact_2(GEN f,GEN p,long d)852 FpX_degfact_2(GEN f, GEN p, long d)
853 {
854 if (!d) return trivfact();
855 if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
856 switch(FpX_quad_factortype(f, p)) {
857 case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
858 case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
859 default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
860 }
861 }
862 static GEN
FpX_factor_2(GEN f,GEN p,long d)863 FpX_factor_2(GEN f, GEN p, long d)
864 {
865 GEN r, s, R, S;
866 long v;
867 int sgn;
868 if (d < 0) pari_err(zeropoler,"FpX_factor_2");
869 if (!d) return trivfact();
870 if (d == 1) return mkmat2(mkcol(f), mkvecsmall(1));
871 r = FpX_quad_root(f, p, 1);
872 if (!r) return mkmat2(mkcol(f), mkvecsmall(1));
873 v = varn(f);
874 s = otherroot(f, r, p);
875 if (signe(r)) r = subii(p, r);
876 if (signe(s)) s = subii(p, s);
877 sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
878 R = deg1pol_i(gen_1, r, v);
879 if (!sgn) return mkmat2(mkcol(R), mkvecsmall(2));
880 S = deg1pol_i(gen_1, s, v);
881 return mkmat2(mkcol2(R,S), mkvecsmall2(1,1));
882 }
883
884 /* factor f mod pp.
885 * flag = 1: return the degrees, not the factors
886 * flag = 2: return NULL if f is not irreducible */
887 static GEN
FpX_factcantor_i(GEN f,GEN pp,long flag)888 FpX_factcantor_i(GEN f, GEN pp, long flag)
889 {
890 long j, e, vf, nbfact, d = degpol(f);
891 ulong p, k;
892 GEN E,y,f2,g,g1,u,v,pd,q;
893 GEN *t;
894
895 if (d <= 2) switch(flag) {
896 case 2: return FpX_is_irred_2(f, pp, d);
897 case 1: return FpX_degfact_2(f, pp, d);
898 default: return FpX_factor_2(f, pp, d);
899 }
900 p = init_p(pp);
901
902 /* to hold factors and exponents */
903 t = (GEN*)cgetg(d+1,t_VEC);
904 E = cgetg(d+1, t_VECSMALL);
905 vf=varn(f); e = nbfact = 1;
906 for(;;)
907 {
908 f2 = FpX_gcd(f,ZX_deriv(f), pp);
909 if (flag > 1 && lg(f2) > 3) return NULL;
910 g1 = FpX_div(f,f2,pp);
911 k = 0;
912 while (lg(g1)>3)
913 {
914 long du,dg;
915 GEN S;
916 k++; if (p && !(k%p)) { k++; f2 = FpX_div(f2,g1,pp); }
917 u = g1; g1 = FpX_gcd(f2,g1, pp);
918 if (lg(g1)>3)
919 {
920 u = FpX_div( u,g1,pp);
921 f2= FpX_div(f2,g1,pp);
922 }
923 du = degpol(u);
924 if (du <= 0) continue;
925
926 /* here u is square-free (product of irred. of multiplicity e * k) */
927 S = init_spec_FpXQ_pow(pp, u);
928 pd=gen_1; v=pol_x[vf];
929 for (d=1; d <= du>>1; d++)
930 {
931 if (!flag) pd = mulii(pd,pp);
932 v = spec_FpXQ_pow(v, pp, S);
933 g = FpX_gcd(gadd(v, gneg(pol_x[vf])), u, pp);
934 dg = degpol(g);
935 if (dg <= 0) continue;
936
937 /* g is a product of irred. pols, all of which have degree d */
938 j = nbfact+dg/d;
939 if (flag)
940 {
941 if (flag > 1) return NULL;
942 for ( ; nbfact<j; nbfact++) { t[nbfact]=(GEN)d; E[nbfact]=e*k; }
943 }
944 else
945 {
946 long r;
947 g = FpX_normalize(g, pp);
948 t[nbfact]=g; q = subis(pd,1); /* also ok for p=2: unused */
949 r = vali(q); q = shifti(q,-r);
950 /* First parameter is an integer m, converted to polynomial w_m, whose
951 * coeffs are its digits in base p (initially m = p --> w_m = X). Take
952 * gcd(g, w_m^(p^d-1)/2), m++, until a factor is found. p = 2 is
953 * treated separately */
954 if (p)
955 split(p,t+nbfact,d,pp,q,r,S);
956 else
957 splitgen(pp,t+nbfact,d,pp,q,r);
958 for (; nbfact<j; nbfact++) E[nbfact]=e*k;
959 }
960 du -= dg;
961 u = FpX_div(u,g,pp);
962 v = FpX_rem(v,u,pp);
963 }
964 if (du)
965 {
966 t[nbfact] = flag? (GEN)du: FpX_normalize(u, pp);
967 E[nbfact++]=e*k;
968 }
969 }
970 j = lg(f2); if (j==3) break;
971 e *= p; f = poldeflate_i(f2, p);
972 }
973 if (flag > 1) return gen_1; /* irreducible */
974 setlg(t, nbfact);
975 setlg(E, nbfact); y = mkvec2((GEN)t, E);
976 return flag ? sort_factor_gen(y, cmpGsGs)
977 : sort_factor(y, cmpii);
978 }
979 GEN
FpX_factcantor(GEN f,GEN pp,long flag)980 FpX_factcantor(GEN f, GEN pp, long flag)
981 {
982 pari_sp av = avma;
983 GEN z = FpX_factcantor_i(FpX_factmod_init(f,pp),pp,flag);
984 if (flag == 2) { avma = av; return z; }
985 return gerepileupto(av, z);
986 }
987 GEN
factcantor0(GEN f,GEN pp,long flag)988 factcantor0(GEN f, GEN pp, long flag)
989 {
990 pari_sp av = avma;
991 long j, nbfact;
992 GEN z, y, t, E, u, v;
993 if (! factmod_init(&f, pp)) { avma=av; return trivfact(); }
994 z = FpX_factcantor_i(f,pp,flag); t = gel(z,1); E = gel(z,2);
995 y = cgetg(3, t_MAT); nbfact = lg(t);
996 u = cgetg(nbfact,t_COL); gel(y,1) = u;
997 v = cgetg(nbfact,t_COL); gel(y,2) = v;
998 if (flag)
999 for (j=1; j<nbfact; j++)
1000 {
1001 gel(u,j) = utoi((ulong)t[j]);
1002 gel(v,j) = utoi((ulong)E[j]);
1003 }
1004 else
1005 for (j=1; j<nbfact; j++)
1006 {
1007 gel(u,j) = FpX_to_mod(gel(t,j), pp);
1008 gel(v,j) = utoi((ulong)E[j]);
1009 }
1010 return gerepileupto(av, y);
1011 }
1012
1013 /* Use this function when you think f is reducible, and that there are lots of
1014 * factors. If you believe f has few factors, use FpX_nbfact(f,p)==1 instead */
1015 long
FpX_is_irred(GEN f,GEN p)1016 FpX_is_irred(GEN f, GEN p) {
1017 return !!FpX_factcantor_i(FpX_factmod_init(f,p),p,2);
1018 }
1019 GEN
FpX_degfact(GEN f,GEN p)1020 FpX_degfact(GEN f, GEN p) {
1021 pari_sp av = avma;
1022 GEN z = FpX_factcantor_i(FpX_factmod_init(f,p),p,1);
1023 settyp(z[1], t_VECSMALL);
1024 settyp(z, t_MAT); return gerepilecopy(av, z);
1025 }
1026 GEN
factcantor(GEN f,GEN p)1027 factcantor(GEN f, GEN p) { return factcantor0(f,p,0); }
1028 GEN
simplefactmod(GEN f,GEN p)1029 simplefactmod(GEN f, GEN p) { return factcantor0(f,p,1); }
1030
1031 GEN
col_to_ff(GEN x,long v)1032 col_to_ff(GEN x, long v)
1033 {
1034 long i, k = lg(x);
1035 GEN p;
1036
1037 while (--k && gcmp0(gel(x,k)));
1038 if (k <= 1) return k? gel(x,1): gen_0;
1039 i = k+2; p = cgetg(i,t_POL);
1040 p[1] = evalsigne(1) | evalvarn(v);
1041 x--; for (k=2; k<i; k++) p[k] = x[k];
1042 return p;
1043 }
1044
1045 /* set x <-- x + c*y mod p */
1046 /* x is not required to be normalized.*/
1047 static void
Flx_addmul_inplace(GEN gx,GEN gy,ulong c,ulong p)1048 Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
1049 {
1050 long i, lx, ly;
1051 ulong *x=(ulong *)gx;
1052 ulong *y=(ulong *)gy;
1053 if (!c) return;
1054 lx = lg(gx);
1055 ly = lg(gy);
1056 if (lx<ly) pari_err(bugparier,"lx<ly in Flx_addmul_inplace");
1057 if (u_OK_ULONG(p))
1058 for (i=2; i<ly; i++) x[i] = (x[i] + c*y[i]) % p;
1059 else
1060 for (i=2; i<ly; i++) x[i] = Fl_add(x[i], Fl_mul(c,y[i],p),p);
1061 }
1062
1063 static long
small_rand(ulong p)1064 small_rand(ulong p)
1065 {
1066 return (p==2)? ((pari_rand31() & 0x1000) == 0): pari_rand31() % p;
1067 }
1068
1069 GEN
FpX_rand(long d1,long v,GEN p)1070 FpX_rand(long d1, long v, GEN p)
1071 {
1072 long i, d = d1+2;
1073 GEN y = cgetg(d,t_POL); y[1] = evalsigne(1) | evalvarn(v);
1074 for (i=2; i<d; i++) gel(y,i) = genrand(p);
1075 (void)normalizepol_i(y,d); return y;
1076 }
1077
1078 /* return a random polynomial in F_q[v], degree < d1 */
1079 GEN
FqX_rand(long d1,long v,GEN T,GEN p)1080 FqX_rand(long d1, long v, GEN T, GEN p)
1081 {
1082 long i, d = d1+2, k = degpol(T), w = varn(T);
1083 GEN y = cgetg(d,t_POL); y[1] = evalsigne(1) | evalvarn(v);
1084 for (i=2; i<d; i++) gel(y,i) = FpX_rand(k, w, p);
1085 (void)normalizepol_i(y,d); return y;
1086 }
1087
1088 #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
1089
1090 long
FpX_split_Berlekamp(GEN * t,GEN p)1091 FpX_split_Berlekamp(GEN *t, GEN p)
1092 {
1093 GEN u = *t, a,b,po2,vker;
1094 long d, i, ir, L, la, lb, vu = varn(u);
1095 long l = lg(u);
1096 ulong ps = itou_or_0(p);
1097 if (ps)
1098 {
1099 vker = Flx_Berlekamp_ker(ZX_to_Flx(u,ps),ps);
1100 vker = Flm_to_FlxV(vker, u[1]);
1101 }
1102 else
1103 {
1104 vker = FpX_Berlekamp_ker(u,p);
1105 vker = RgM_to_RgXV(vker,vu);
1106 }
1107 d = lg(vker)-1;
1108 po2 = shifti(p, -1); /* (p-1) / 2 */
1109 ir = 0;
1110 /* t[i] irreducible for i <= ir, still to be treated for i < L */
1111 for (L=1; L<d; )
1112 {
1113 GEN polt;
1114 if (ps)
1115 {
1116 GEN pol = const_vecsmall(l-2,0);
1117 pol[1] = u[1];
1118 pol[2] = small_rand(ps); /*Assume vker[1]=1*/
1119 for (i=2; i<=d; i++)
1120 Flx_addmul_inplace(pol, gel(vker,i), (ulong)small_rand(ps), ps);
1121 (void)Flx_renormalize((GEN)pol,l-1);
1122
1123 polt = Flx_to_ZX(pol);
1124 }
1125 else
1126 {
1127 GEN pol = scalarpol(genrand(p), vu);
1128 for (i=2; i<=d; i++)
1129 pol = ZX_add(pol, ZX_Z_mul(gel(vker,i), randomi(p)));
1130 polt = FpX_red(pol,p);
1131 }
1132 for (i=ir; i<L && L<d; i++)
1133 {
1134 a = t[i]; la = degpol(a);
1135 if (la == 1) { set_irred(i); }
1136 else if (la == 2)
1137 {
1138 GEN r = FpX_quad_root(a,p,1);
1139 if (r)
1140 {
1141 t[i] = deg1pol_i(gen_1, subii(p,r), vu); r = otherroot(a,r,p);
1142 t[L] = deg1pol_i(gen_1, subii(p,r), vu); L++;
1143 }
1144 set_irred(i);
1145 }
1146 else
1147 {
1148 pari_sp av = avma;
1149 b = FpX_rem(polt, a, p);
1150 if (degpol(b) <= 0) { avma=av; continue; }
1151 b = ZX_Z_add(FpXQ_pow(b,po2, a,p), gen_m1);
1152 b = FpX_gcd(a,b, p); lb = degpol(b);
1153 if (lb && lb < la)
1154 {
1155 b = FpX_normalize(b, p);
1156 t[L] = FpX_div(a,b,p);
1157 t[i]= b; L++;
1158 }
1159 else avma = av;
1160 }
1161 }
1162 }
1163 return d;
1164 }
1165
1166 GEN
FqX_gcd(GEN P,GEN Q,GEN T,GEN p)1167 FqX_gcd(GEN P,GEN Q,GEN T,GEN p) {return T? FpXQX_gcd(P,Q,T,p): FpX_gcd(P,Q,p);}
1168 long
FqX_is_squarefree(GEN P,GEN T,GEN p)1169 FqX_is_squarefree(GEN P, GEN T, GEN p)
1170 {
1171 pari_sp av = avma;
1172 GEN z = FqX_gcd(P, FqX_deriv(P, T, p), T, p);
1173 avma = av;
1174 return degpol(z)==0;
1175 }
1176
1177 long
FqX_split_Berlekamp(GEN * t,GEN q,GEN T,GEN p)1178 FqX_split_Berlekamp(GEN *t, GEN q, GEN T, GEN p)
1179 {
1180 GEN u = *t, a,b,qo2,vker,pol;
1181 long N = degpol(u), vu = varn(u), vT = varn(T), dT = degpol(T);
1182 long d, i, ir, L, la, lb;
1183
1184 vker = FqX_Berlekamp_ker(u,T,q,p);
1185 vker = RgM_to_RgXV(vker,vu);
1186 d = lg(vker)-1;
1187 qo2 = shifti(q, -1); /* (q-1) / 2 */
1188 pol = cgetg(N+3,t_POL);
1189 ir = 0;
1190 /* t[i] irreducible for i < ir, still to be treated for i < L */
1191 for (L=1; L<d; )
1192 {
1193 GEN polt;
1194 gel(pol,2) = FpX_rand(dT,vT,p);
1195 setlg(pol, signe(pol[2])? 3: 2);
1196 pol[1] = u[1];
1197 for (i=2; i<=d; i++)
1198 pol = gadd(pol, gmul(gel(vker,i), FpX_rand(dT,vT,p)));
1199 polt = FpXQX_red(pol,T,p);
1200 for (i=ir; i<L && L<d; i++)
1201 {
1202 a = t[i]; la = degpol(a);
1203 if (la == 1) { set_irred(i); }
1204 else
1205 {
1206 pari_sp av = avma;
1207 b = FqX_rem(polt, a, T,p);
1208 if (!degpol(b)) { avma=av; continue; }
1209 b = FpXQYQ_pow(b,qo2, a,T,p);
1210 if (!degpol(b)) { avma=av; continue; }
1211 gel(b,2) = gadd(gel(b,2), gen_1);
1212 b = FqX_gcd(a,b, T,p); lb = degpol(b);
1213 if (lb && lb < la)
1214 {
1215 b = FqX_normalize(b, T,p);
1216 t[L] = FqX_div(a,b,T,p);
1217 t[i]= b; L++;
1218 }
1219 else avma = av;
1220 }
1221 }
1222 }
1223 return d;
1224 }
1225
1226 static GEN
FpX_factor_i(GEN f,GEN pp)1227 FpX_factor_i(GEN f, GEN pp)
1228 {
1229 long e, N, nbfact, val, d = degpol(f);
1230 ulong p, k, j;
1231 GEN E, f2, g1, u, *t;
1232
1233 if (d <= 2) return FpX_factor_2(f, pp, d);
1234 p = init_p(pp);
1235
1236 /* to hold factors and exponents */
1237 t = (GEN*)cgetg(d+1,t_COL); E = cgetg(d+1,t_VECSMALL);
1238 val = ZX_valuation(f, &f);
1239 e = nbfact = 1;
1240 if (val) { t[1] = pol_x[varn(f)]; E[1] = val; nbfact++; }
1241
1242 for(;;)
1243 {
1244 f2 = FpX_gcd(f,ZX_deriv(f), pp);
1245 g1 = lg(f2)==3? f: FpX_div(f,f2,pp); /* is squarefree */
1246 k = 0;
1247 while (lg(g1)>3)
1248 {
1249 k++; if (p && !(k%p)) { k++; f2 = FpX_div(f2,g1,pp); }
1250 u = g1;
1251 if (!degpol(f2)) g1 = pol_1[0]; /* only its degree (= 0) matters */
1252 else
1253 {
1254 g1 = FpX_gcd(f2,g1, pp);
1255 if (degpol(g1))
1256 {
1257 u = FpX_div( u,g1,pp);
1258 f2= FpX_div(f2,g1,pp);
1259 }
1260 }
1261 /* u is square-free (product of irred. of multiplicity e * k) */
1262 N = degpol(u);
1263 if (N > 0)
1264 {
1265 t[nbfact] = FpX_normalize(u,pp);
1266 d = (N==1)? 1: FpX_split_Berlekamp(t+nbfact, pp);
1267 for (j=0; j<(ulong)d; j++) E[nbfact+j] = e*k;
1268 nbfact += d;
1269 }
1270 }
1271 if (!p) break;
1272 j = degpol(f2); if (!j) break;
1273 if (j % p) pari_err(talker, "factmod: %lu is not prime", p);
1274 e *= p; f = poldeflate_i(f2, p);
1275 }
1276 setlg(t, nbfact);
1277 setlg(E, nbfact); return sort_factor(mkvec2((GEN)t,E), cmpii);
1278 }
1279 GEN
FpX_factor(GEN f,GEN p)1280 FpX_factor(GEN f, GEN p)
1281 {
1282 pari_sp av = avma;
1283 return gerepilecopy(av, FpX_factor_i(FpX_factmod_init(f, p), p));
1284 }
1285
1286 GEN
factmod(GEN f,GEN p)1287 factmod(GEN f, GEN p)
1288 {
1289 pari_sp av = avma;
1290 long nbfact;
1291 long j;
1292 GEN y, u, v, z, t, E;
1293
1294 if (!factmod_init(&f, p)) { avma = av; return trivfact(); }
1295 z = FpX_factor_i(f, p); t = gel(z,1); E = gel(z,2);
1296 y = cgetg(3,t_MAT); nbfact = lg(t);
1297 u = cgetg(nbfact,t_COL); gel(y,1) = u;
1298 v = cgetg(nbfact,t_COL); gel(y,2) = v;
1299 for (j=1; j<nbfact; j++)
1300 {
1301 gel(u,j) = FpX_to_mod(gel(t,j), p);
1302 gel(v,j) = utoi((ulong)E[j]);
1303 }
1304 return gerepileupto(av, y);
1305 }
1306 GEN
factormod0(GEN f,GEN p,long flag)1307 factormod0(GEN f, GEN p, long flag)
1308 {
1309 switch(flag)
1310 {
1311 case 0: return factmod(f,p);
1312 case 1: return simplefactmod(f,p);
1313 default: pari_err(flagerr,"factormod");
1314 }
1315 return NULL; /* not reached */
1316 }
1317
1318 /*******************************************************************/
1319 /* */
1320 /* CONVERSIONS RELATED TO p-ADICS */
1321 /* */
1322 /*******************************************************************/
1323 static GEN
Zp_to_Z(GEN x)1324 Zp_to_Z(GEN x) {
1325 switch(typ(x))
1326 {
1327 case t_INT: break;
1328 case t_PADIC: x = gtrunc(x); break;
1329 default: pari_err(typeer,"QpX_to_ZX");
1330 }
1331 return x;
1332 }
1333 static GEN
ZpX_to_ZX(GEN f)1334 ZpX_to_ZX(GEN f) {
1335 long i, l = lg(f);
1336 GEN F = cgetg(l, t_POL); F[1] = f[1];
1337 for (i=2; i<l; i++) gel(f,i) = Zp_to_Z(gel(f,i));
1338 return f;
1339 }
1340 /* make f suitable for [root|factor]padic */
1341 static GEN
QpX_to_ZX(GEN f)1342 QpX_to_ZX(GEN f)
1343 {
1344 GEN c = content(f);
1345 if (gcmp0(c)) /* O(p^n) can occur */
1346 {
1347 if (typ(c) != t_PADIC) pari_err(typeer,"QpX_to_ZX");
1348 f = gdiv(f, gpowgs(gel(c,2), valp(c)));
1349 }
1350 else
1351 f = gdiv(f,c);
1352 return ZpX_to_ZX(f);
1353 }
1354
1355 /* x in Z return x + O(pr), pr = p^r. Return gen_0 instead of zeropadic */
1356 static GEN
Z_to_Zp(GEN x,GEN p,GEN pr,long r)1357 Z_to_Zp(GEN x, GEN p, GEN pr, long r)
1358 {
1359 GEN y;
1360 long v, sx = signe(x);
1361
1362 if (!sx) return gen_0;
1363 v = Z_pvalrem(x,p,&x);
1364 r -= v; if (r <= 0) return gen_0;
1365 y = cgetg(5,t_PADIC);
1366 y[1] = evalprecp(r)|evalvalp(v);
1367 gel(y,2) = p;
1368 gel(y,3) = pr;
1369 gel(y,4) = modii(x,pr); return y;
1370 }
1371 static GEN
ZV_to_ZpV(GEN z,GEN p,long prec)1372 ZV_to_ZpV(GEN z, GEN p, long prec)
1373 {
1374 long i, l = lg(z);
1375 GEN Z = cgetg(l, typ(z)), q = powiu(p, prec);
1376 for (i=1; i<lg(z); i++) gel(Z,i) = Z_to_Zp(gel(z,i),p,q,prec);
1377 return Z;
1378 }
1379 static GEN
ZX_to_ZpX(GEN z,GEN p,GEN q,long prec)1380 ZX_to_ZpX(GEN z, GEN p, GEN q, long prec)
1381 {
1382 long i, l = lg(z);
1383 GEN Z = cgetg(l, t_POL); Z[1] = z[1];
1384 for (i=2; i<lg(z); i++) gel(Z,i) = Z_to_Zp(gel(z,i),p,q,prec);
1385 return Z;
1386 }
1387 /* return (x + O(p^r)) normalized (multiply by a unit such that leading coeff
1388 * is a power of p), x in Z[X] (or Z_p[X]) */
1389 static GEN
ZX_to_ZpX_normalized(GEN x,GEN p,GEN pr,long r)1390 ZX_to_ZpX_normalized(GEN x, GEN p, GEN pr, long r)
1391 {
1392 long i, lx = lg(x);
1393 GEN z, lead = leading_term(x);
1394
1395 if (gcmp1(lead)) return ZX_to_ZpX(x, p, pr, r);
1396 (void)Z_pvalrem(lead, p, &lead); lead = Fp_inv(lead, pr);
1397 z = cgetg(lx,t_POL);
1398 for (i=2; i < lx; i++) gel(z,i) = Z_to_Zp(mulii(lead,gel(x,i)),p,pr,r);
1399 z[1] = x[1]; return z;
1400 }
1401 static GEN
ZXV_to_ZpXQV(GEN z,GEN T,GEN p,long prec)1402 ZXV_to_ZpXQV(GEN z, GEN T, GEN p, long prec)
1403 {
1404 long i, l = lg(z);
1405 GEN Z = cgetg(l, typ(z)), q = powiu(p, prec);
1406 T = gcopy(T);
1407 for (i=1; i<lg(z); i++) gel(Z,i) = mkpolmod(ZX_to_ZpX(gel(z,i),p,q,prec),T);
1408 return Z;
1409 }
1410
1411 static GEN
QpXQ_to_ZXY(GEN f)1412 QpXQ_to_ZXY(GEN f)
1413 {
1414 GEN c = content(f);
1415 long i, l = lg(f);
1416 if (gcmp0(c)) /* O(p^n) can occur */
1417 {
1418 if (typ(c) != t_PADIC) pari_err(typeer,"QpXQ_to_ZXY");
1419 f = gdiv(f, gpowgs(gel(c,2), valp(c)));
1420 }
1421 else
1422 f = gdiv(f,c);
1423 for (i=2; i<l; i++)
1424 {
1425 GEN t = gel(f,i);
1426 switch(typ(t))
1427 {
1428 case t_POL: t = ZpX_to_ZX(t); break;
1429 default: t = Zp_to_Z(t); break;
1430 }
1431 gel(f,i) = t;
1432 }
1433 return f;
1434 }
1435
1436 /*******************************************************************/
1437 /* */
1438 /* p-ADIC ROOTS */
1439 /* */
1440 /*******************************************************************/
1441
1442 /* f primitive ZX, squarefree, leading term prime to p. a in Z such that
1443 * f(a) = 0 mod p. Return p-adic roots of f equal to a mod p, in
1444 * precision >= prec */
1445 static GEN
ZX_Zp_root(GEN f,GEN a,GEN p,long prec)1446 ZX_Zp_root(GEN f, GEN a, GEN p, long prec)
1447 {
1448 GEN z, R, a0 = modii(a, p);
1449 long i, j, k;
1450
1451 if (signe(FpX_eval(FpX_deriv(f, p), a0, p)))
1452 { /* simple zero mod p, go all the way to p^prec */
1453 if (prec > 1) a0 = ZpX_liftroot(f, a0, p, prec);
1454 return mkcol(a0);
1455 }
1456
1457 f = poleval(f, gadd(a, gmul(p,pol_x[varn(f)])));
1458 f = gdivexact(f, powiu(p,ggval(f, p)));
1459 z = cgetg(degpol(f)+1,t_COL);
1460
1461 R = FpX_roots(f, p);
1462 for (j=i=1; i<lg(R); i++)
1463 {
1464 GEN u = ZX_Zp_root(f, gel(R,i), p, prec-1);
1465 for (k=1; k<lg(u); k++) gel(z,j++) = gadd(a, gmul(p, gel(u,k)));
1466 }
1467 setlg(z,j); return z;
1468 }
1469
1470 /* a t_PADIC, return vector of p-adic roots of f equal to a (mod p)
1471 * We assume 1) f(a) = 0 mod p (mod 4 if p=2).
1472 * 2) leading coeff prime to p. */
1473 GEN
Zp_appr(GEN f,GEN a)1474 Zp_appr(GEN f, GEN a)
1475 {
1476 pari_sp av = avma;
1477 long prec;
1478 GEN z, p;
1479 if (typ(f) != t_POL) pari_err(notpoler,"Zp_appr");
1480 if (gcmp0(f)) pari_err(zeropoler,"Zp_appr");
1481 if (typ(a) != t_PADIC) pari_err(typeer,"Zp_appr");
1482 p = gel(a,2); prec = gcmp0(a)? valp(a): precp(a);
1483 f = QpX_to_ZX(f);
1484 z = modulargcd(f, ZX_deriv(f));
1485 if (degpol(z) > 0) f = RgX_div(f,z);
1486 z = ZX_Zp_root(f, gtrunc(a), p, prec);
1487 return gerepilecopy(av, ZV_to_ZpV(z, p, prec));
1488 }
1489 /* vector of p-adic roots of the ZX f, leading term prime to p */
1490 static GEN
ZX_Zp_roots(GEN f,GEN p,long prec)1491 ZX_Zp_roots(GEN f, GEN p, long prec)
1492 {
1493 GEN y, z, rac;
1494 long lx, i, j, k;
1495
1496 z = modulargcd(f, ZX_deriv(f));
1497 if (degpol(z) > 0) f = RgX_div(f,z);
1498 rac = FpX_roots(f, p);
1499 lx = lg(rac); if (lx == 1) return rac;
1500 y = cgetg(degpol(f)+1,t_COL);
1501 for (j=i=1; i<lx; i++)
1502 {
1503 z = ZX_Zp_root(f, gel(rac,i), p, prec);
1504 for (k=1; k<lg(z); k++,j++) gel(y,j) = gel(z,k);
1505 }
1506 setlg(y,j); return ZV_to_ZpV(y, p, prec);
1507 }
1508
1509 /* f a ZX */
1510 static GEN
pnormalize(GEN f,GEN p,long prec,long n,GEN * plead,long * pprec,int * prev)1511 pnormalize(GEN f, GEN p, long prec, long n, GEN *plead, long *pprec, int *prev)
1512 {
1513 *plead = leading_term(f);
1514 *pprec = prec;
1515 *prev = 0;
1516 if (!is_pm1(*plead))
1517 {
1518 long v = ggval(*plead,p), v1 = ggval(constant_term(f),p);
1519 if (v1 < v)
1520 {
1521 *prev = 1; f = polrecip_i(f);
1522 /* beware loss of precision from lc(factor), whose valuation is <= v */
1523 *pprec += v; v = v1;
1524 }
1525 *pprec += v * n;
1526 }
1527 return pol_to_monic(f, plead);
1528 }
1529
1530 /* return p-adic roots of f, precision prec */
1531 GEN
rootpadic(GEN f,GEN p,long prec)1532 rootpadic(GEN f, GEN p, long prec)
1533 {
1534 pari_sp av = avma;
1535 GEN lead,y;
1536 long PREC,i,k;
1537 int reverse;
1538
1539 if (typ(p)!=t_INT) pari_err(typeer,"rootpadic");
1540 if (typ(f)!=t_POL) pari_err(notpoler,"rootpadic");
1541 if (gcmp0(f)) pari_err(zeropoler,"rootpadic");
1542 if (prec <= 0) pari_err(talker,"non-positive precision in rootpadic");
1543 f = QpX_to_ZX(f);
1544 f = pnormalize(f, p, prec, 1, &lead, &PREC, &reverse);
1545 y = ZX_Zp_roots(f,p,PREC);
1546 k = lg(y);
1547 if (lead)
1548 for (i=1; i<k; i++) gel(y,i) = gdiv(gel(y,i), lead);
1549 if (reverse)
1550 for (i=1; i<k; i++) gel(y,i) = ginv(gel(y,i));
1551 return gerepilecopy(av, y);
1552 }
1553 /*************************************************************************/
1554 /* rootpadicfast */
1555 /*************************************************************************/
1556
1557 /* lift accelerator */
1558 long
hensel_lift_accel(long n,long * pmask)1559 hensel_lift_accel(long n, long *pmask)
1560 {
1561 long a, j, mask = 0;
1562 for(j=BITS_IN_LONG-1, a=n ;; j--)
1563 {
1564 mask |= (a&1)<<j;
1565 a = (a>>1) + (a&1);
1566 if (a==1) break;
1567 }
1568 *pmask = mask>>j;
1569 return BITS_IN_LONG-j;
1570 }
1571 /* SPEC:
1572 * p is a t_INT > 1, e >= 0
1573 * f is a ZX with leading term prime to p.
1574 * a is a simple root mod l for all l|p.
1575 * Return roots of f mod p^e, as integers (implicitly mod p^e)
1576 * STANDARD USE: p is a prime power */
1577 GEN
ZpX_liftroot(GEN f,GEN a,GEN p,long e)1578 ZpX_liftroot(GEN f, GEN a, GEN p, long e)
1579 {
1580 pari_sp ltop=avma;
1581 GEN qold, q, qm1;
1582 GEN W, fr, ar, Wr = gen_0;
1583 long i, nb, mask;
1584 qold = q = p; qm1 = gen_1;
1585 nb = hensel_lift_accel(e, &mask);
1586 fr = FpX_red(f,q);
1587 a = modii(a,q);
1588 W = FpX_eval(ZX_deriv(fr), a, q);
1589 W = Fp_inv(W,q);
1590 for(i=0;i<nb;i++)
1591 {
1592 qm1 = (mask&(1<<i))?sqri(qm1):mulii(qm1, q);
1593 q = mulii(qm1, p);
1594 fr = FpX_red(f,q);
1595 ar = a;
1596 if (i)
1597 {
1598 W = modii(mulii(Wr,FpX_eval(ZX_deriv(fr),ar,qold)), qold);
1599 W = modii(mulii(Wr, subsi(2,W)), qold);
1600 }
1601 Wr = W;
1602 a = subii(ar, mulii(Wr, FpX_eval(fr, ar,q)));
1603 a = modii(a,q);
1604 qold = q;
1605 }
1606 return gerepileupto(ltop,a);
1607 }
1608 GEN
ZpXQX_liftroot(GEN f,GEN a,GEN T,GEN p,long e)1609 ZpXQX_liftroot(GEN f, GEN a, GEN T, GEN p, long e)
1610 {
1611 pari_sp ltop=avma;
1612 GEN qold, q, qm1;
1613 GEN W, fr, ar, Wr = gen_0;
1614 long i, nb, mask;
1615 qold = p ; q = p; qm1 = gen_1;
1616 nb=hensel_lift_accel(e, &mask);
1617 fr = FpXQX_red(f, T, q);
1618 a = Fq_red(a, T, q);
1619 W = FqX_eval(derivpol(fr), a, T, q);
1620 W = Fq_inv(W,T,q);
1621 for(i=0;i<nb;i++)
1622 {
1623 qm1 = (mask&(1<<i))?sqri(qm1):mulii(qm1, q);
1624 q = mulii(qm1, p);
1625 fr = FpXQX_red(f,T,q);
1626 ar = a;
1627 if (i)
1628 {
1629 W = Fq_red(gmul(Wr,FqX_eval(derivpol(fr),ar,T,qold)), T, qold);
1630 W = Fq_red(gmul(Wr, gadd(gen_2, gneg(W))), T, qold);
1631 }
1632 Wr = W;
1633 a = gadd(ar, gmul(gneg(Wr), FqX_eval(fr, ar, T, q)));
1634 a = Fq_red(a, T, q);
1635 qold = q;
1636 }
1637 return gerepileupto(ltop,a);
1638 }
1639 /* Apply ZpX_liftroot to all roots in S and trace trick.
1640 * Elements of S must be distinct simple roots mod p for all p|q. */
1641 GEN
ZpX_liftroots(GEN f,GEN S,GEN q,long e)1642 ZpX_liftroots(GEN f, GEN S, GEN q, long e)
1643 {
1644 long i, d, l = lg(S), n = l-1;
1645 GEN y = cgetg(l, typ(S));
1646 if (!n) return y;
1647 for (i=1; i<n; i++)
1648 gel(y,i) = ZpX_liftroot(f, gel(S,i), q, e);
1649 d = degpol(f);
1650 if (n != d) /* not totally split*/
1651 gel(y,n) = ZpX_liftroot(f, gel(S,n), q, e);
1652 else
1653 { /* totally split: use trace trick */
1654 pari_sp av = avma;
1655 GEN z = gel(f, d+1);/* -trace(roots) */
1656 for(i=1; i<n;i++) z = addii(z, gel(y,i));
1657 z = modii(negi(z), powiu(q,e));
1658 gel(y,n) = gerepileuptoint(av,z);
1659 }
1660 return y;
1661 }
1662 /* p is prime
1663 * f in a ZX, with leading term prime to p.
1664 * f must have no multiple roots mod p.
1665 *
1666 * return p-adics roots of f with prec p^e, as integers (implicitly mod p^e) */
1667 GEN
rootpadicfast(GEN f,GEN p,long e)1668 rootpadicfast(GEN f, GEN p, long e)
1669 {
1670 pari_sp av = avma;
1671 GEN y, S = FpX_roots(f, p); /*no multiplicity*/
1672 if (lg(S)==1) { avma = av; return cgetg(1,t_COL); }
1673 S = gclone(S); avma = av;
1674 y = ZpX_liftroots(f,S,p,e);
1675 gunclone(S); return y;
1676 }
1677 /* Same as ZpX_liftroot for the polynomial X^n-T
1678 * TODO: generalize to sparse polynomials. */
1679 GEN
padicsqrtnlift(GEN T,GEN n,GEN a,GEN p,long e)1680 padicsqrtnlift(GEN T, GEN n, GEN a, GEN p, long e)
1681 {
1682 pari_sp ltop=avma;
1683 GEN qold, q, qm1;
1684 GEN W, ar, Wr = gen_0;
1685 long i, nb, mask;
1686 qold = q = p; qm1 = gen_1;
1687 nb = hensel_lift_accel(e, &mask);
1688 W = modii(mulii(n,Fp_pow(a,subis(n,1),q)),q);
1689 W = Fp_inv(W,q);
1690 for(i=0;i<nb;i++)
1691 {
1692 qm1 = (mask&(1<<i))?sqri(qm1):mulii(qm1, q);
1693 q = mulii(qm1, p);
1694 ar = a;
1695 if (i)
1696 {
1697 W = modii(mulii(Wr,mulii(n,Fp_pow(ar,subis(n,1),qold))),qold);
1698 W = modii(mulii(Wr, subsi(2,W)),qold);
1699 }
1700 Wr = W;
1701 a = subii(ar, mulii(Wr, subii(Fp_pow(ar,n,q),T)));
1702 a = modii(a,q);
1703 qold = q;
1704 }
1705 return gerepileupto(ltop,a);
1706 }
1707 /**************************************************************************/
1708
1709 static void
scalar_getprec(GEN x,long * pprec,GEN * pp)1710 scalar_getprec(GEN x, long *pprec, GEN *pp)
1711 {
1712 if (typ(x)==t_PADIC)
1713 {
1714 long e = valp(x); if (signe(x[4])) e += precp(x);
1715 if (e < *pprec) *pprec = e;
1716 if (*pp && !equalii(*pp, gel(x,2))) pari_err(consister,"apprpadic");
1717 *pp = gel(x,2);
1718 }
1719 }
1720
1721 static void
getprec(GEN x,long * pprec,GEN * pp)1722 getprec(GEN x, long *pprec, GEN *pp)
1723 {
1724 long i;
1725 if (typ(x) != t_POL) scalar_getprec(x, pprec, pp);
1726 else
1727 for (i = lg(x)-1; i>1; i--) scalar_getprec(gel(x,i), pprec, pp);
1728 }
1729 static GEN
ZXY_ZpQ_root(GEN f,GEN a,GEN T,GEN p,long prec)1730 ZXY_ZpQ_root(GEN f, GEN a, GEN T, GEN p, long prec)
1731 {
1732 GEN z, R;
1733 long i, j, k, lR;
1734 if (signe(FqX_eval(FqX_deriv(f,T,p), a, T,p)))
1735 { /* simple zero mod (T,p), go all the way to p^prec */
1736 if (prec > 1) a = ZpXQX_liftroot(f, a, T, p, prec);
1737 return mkcol(a);
1738 }
1739 /* TODO: need RgX_RgYQX_compo ? */
1740 f = lift_intern(poleval(f, gadd(mkpolmod(a,T), gmul(p, pol_x[varn(f)]))));
1741 f = gdiv(f, powiu(p, ggval(f,p)));
1742 z = cgetg(degpol(f)+1,t_COL);
1743
1744 #if 1 /* TODO: need a public FqX_roots */
1745 lR = FqX_split_deg1(&R, FqX_red(f, T, p), powiu(p, degpol(T)), T, p) + 1;
1746 R = roots_from_deg1(FqX_split_roots(R, T, p, NULL));
1747 #else
1748 R = FqX_factor(FqX_red(f, T, p), T, p);
1749 R = gel(R,1); lR = lg(R);
1750 for (i=1; i<lR; i++)
1751 {
1752 GEN u = gel(R,i); if (degpol(u) > 1) break;
1753 gel(R,i) = gneg(gel(u,2));
1754 }
1755 lR = i; /* "truncate" R to the deg 1 factors, demoted to roots above */
1756 #endif
1757 for(j=i=1; i<lR; i++)
1758 {
1759 GEN u = ZXY_ZpQ_root(f, gel(R,i), T, p, prec-1);
1760 for (k=1; k<lg(u); k++) gel(z,j++) = gadd(a, gmul(p, gel(u,k)));
1761 }
1762 setlg(z,j); return z;
1763 }
1764
1765 /* a belongs to finite extension of Q_p, return all roots of f equal to a
1766 * mod p. We assume f(a) = 0 (mod p) [mod 4 if p=2] */
1767 GEN
padicappr(GEN f,GEN a)1768 padicappr(GEN f, GEN a)
1769 {
1770 GEN p, z, T;
1771 long prec;
1772 pari_sp av = avma;
1773
1774 switch(typ(a)) {
1775 case t_PADIC: return Zp_appr(f,a);
1776 case t_POLMOD: break;
1777 default: pari_err(typeer,"padicappr");
1778 }
1779 if (typ(f)!=t_POL) pari_err(notpoler,"padicappr");
1780 if (gcmp0(f)) pari_err(zeropoler,"padicappr");
1781 z = ggcd(f, derivpol(f));
1782 if (degpol(z) > 0) f = RgX_div(f,z);
1783 T = gel(a,1); a = gel(a,2);
1784 p = NULL; prec = BIGINT;
1785 getprec(a, &prec, &p);
1786 getprec(T, &prec, &p); if (!p) pari_err(typeer,"padicappr");
1787 f = QpXQ_to_ZXY(lift_intern(f));
1788 a = QpX_to_ZX(a);
1789 T = QpX_to_ZX(T);
1790 z = ZXY_ZpQ_root(f, a, T, p, prec);
1791 return gerepilecopy(av, ZXV_to_ZpXQV(z, T, p, prec));
1792 }
1793
1794 /*******************************************************************/
1795 /* */
1796 /* FACTORIZATION in Zp[X] */
1797 /* */
1798 /*******************************************************************/
1799 int
cmp_padic(GEN x,GEN y)1800 cmp_padic(GEN x, GEN y)
1801 {
1802 long vx, vy;
1803 if (x == gen_0) return -1;
1804 if (y == gen_0) return 1;
1805 vx = valp(x);
1806 vy = valp(y);
1807 if (varncmp(vx, vy) < 0) return 1;
1808 if (varncmp(vx, vy) > 0) return -1;
1809 return cmpii(gel(x,4), gel(y,4));
1810 }
1811
1812 /*******************************/
1813 /* Using Buchmann--Lenstra */
1814 /*******************************/
1815
1816 /* factor T = nf[1] in Zp to precision k */
1817 static GEN
padicff2(GEN nf,GEN p,long k)1818 padicff2(GEN nf,GEN p,long k)
1819 {
1820 GEN z, mat, D, U,fa, pk, dec_p, Ui, mulx;
1821 long i, l;
1822
1823 mulx = eltmul_get_table(nf, gmael(nf,8,2)); /* mul by (x mod T) */
1824 dec_p = primedec(nf,p);
1825 l = lg(dec_p); fa = cgetg(l,t_COL);
1826 D = NULL; /* -Wall */
1827 for (i=1; i<l; i++)
1828 {
1829 GEN P = gel(dec_p,i);
1830 long e = itos(gel(P,3)), ef = e * itos(gel(P,4));
1831 D = smithall(idealpows(nf,P, k*e), &U, NULL);
1832 Ui= ginv(U); setlg(Ui, ef+1); /* cf smithrel */
1833 U = rowslice(U, 1, ef);
1834 mat = gmul(U, gmul(mulx, Ui));
1835 gel(fa,i) = caradj(mat,0,NULL);
1836 }
1837 pk = gcoeff(D,1,1); /* = p^k */
1838 z = cgetg(l,t_COL); pk = icopy(pk);
1839 for (i=1; i<l; i++)
1840 gel(z,i) = ZX_to_ZpX_normalized(gel(fa,i),p,pk,k);
1841 return z;
1842 }
1843
1844 static GEN
padicff(GEN x,GEN p,long pr)1845 padicff(GEN x,GEN p,long pr)
1846 {
1847 pari_sp av = avma;
1848 GEN q, bas, invbas, mul, dK, nf, g, e, dx = absi(ZX_disc(x));
1849 long n = degpol(x), v = Z_pvalrem(dx,p,&q);
1850
1851 nf = cgetg(10,t_VEC); gel(nf,1) = x;
1852 if (is_pm1(q)) {
1853 e = mkcol(utoi(v));
1854 g = mkcol(p);
1855 } else {
1856 e = mkcol2(utoi(v), gen_1);
1857 g = mkcol2(p, q);
1858 }
1859 bas = nfbasis(x, &dK, 0, mkmat2(g,e));
1860 gel(nf,3) = dK;
1861 gel(nf,4) = dvdii( diviiexact(dx, dK), p )? p: gen_1;
1862 invbas = QM_inv(RgXV_to_RgM(bas,n), gen_1);
1863 mul = get_mul_table(x,bas,invbas);
1864 gel(nf,7) = bas;
1865 gel(nf,8) = invbas;
1866 gel(nf,9) = mul;
1867 gel(nf,2) = gel(nf,5) = gel(nf,6) = gen_0;
1868 return gerepileupto(av,padicff2(nf,p,pr));
1869 }
1870
1871 static GEN
padic_trivfact(GEN x,GEN p,long r)1872 padic_trivfact(GEN x, GEN p, long r)
1873 {
1874 return mkmat2(mkcol(ZX_to_ZpX_normalized(x, p, powiu(p,r), r)),
1875 mkcol(gen_1));
1876 }
1877
1878 GEN
factorpadic2(GEN f,GEN p,long prec)1879 factorpadic2(GEN f, GEN p, long prec)
1880 {
1881 pari_sp av = avma;
1882 GEN fa,ex,y;
1883 long n,i,l;
1884
1885 if (typ(f)!=t_POL) pari_err(notpoler,"factorpadic");
1886 if (typ(p)!=t_INT) pari_err(arither1);
1887 if (gcmp0(f)) pari_err(zeropoler,"factorpadic");
1888 if (prec <= 0) pari_err(talker,"non-positive precision in factorpadic");
1889
1890 n = degpol(f);
1891 if (n==0) return trivfact();
1892 f = QpX_to_ZX(f);
1893 if (n==1) return gerepilecopy(av, padic_trivfact(f,p,prec));
1894 if (!gcmp1(leading_term(f)))
1895 pari_err(impl,"factorpadic2 for non-monic polynomial");
1896
1897 fa = ZX_squff(f, &ex);
1898 l = lg(fa); n = 0;
1899 for (i=1; i<l; i++)
1900 {
1901 gel(fa,i) = padicff(gel(fa,i),p,prec);
1902 n += lg(fa[i])-1;
1903 }
1904
1905 y = fact_from_DDF(fa,ex,n);
1906 return gerepileupto(av, sort_factor(y, cmp_padic));
1907 }
1908
1909 /***********************/
1910 /* Using ROUND 4 */
1911 /***********************/
1912 static int
expo_is_squarefree(GEN e)1913 expo_is_squarefree(GEN e)
1914 {
1915 long i, l = lg(e);
1916 for (i=1; i<l; i++)
1917 if (e[i] != 1) return 0;
1918 return 1;
1919 }
1920
1921 /* assume f a ZX with leading_term 1, degree > 0 */
1922 GEN
ZX_monic_factorpadic(GEN f,GEN p,long prec)1923 ZX_monic_factorpadic(GEN f, GEN p, long prec)
1924 {
1925 GEN w, poly, p1, p2, ex, P, E;
1926 long n=degpol(f), i, k, j, pr;
1927
1928 if (n==1) return mkmat2(mkcol(f), mkcol(gen_1));
1929 pr = prec;
1930
1931 poly = ZX_squff(f,&ex);
1932 P = cgetg(n+1,t_COL);
1933 E = cgetg(n+1,t_COL); n = lg(poly);
1934 for (j=i=1; i<n; i++)
1935 {
1936 pari_sp av1 = avma;
1937 GEN fx = gel(poly,i), fa = FpX_factor(fx,p);
1938 w = gel(fa,1);
1939 if (expo_is_squarefree(gel(fa,2)))
1940 { /* no repeated factors: Hensel lift */
1941 p1 = hensel_lift_fact(fx, w, NULL, p, powiu(p,pr), pr);
1942 p2 = utoipos(ex[i]);
1943 for (k=1; k<lg(p1); k++,j++)
1944 {
1945 P[j] = p1[k];
1946 gel(E,j) = p2;
1947 }
1948 continue;
1949 }
1950 /* use Round 4 */
1951 p2 = maxord_i(p, fx, Z_pval(ZX_disc(fx),p), w, pr);
1952 if (p2)
1953 {
1954 p2 = gerepilecopy(av1,p2);
1955 p1 = gel(p2,1);
1956 p2 = gel(p2,2);
1957 for (k=1; k<lg(p1); k++,j++)
1958 {
1959 P[j] = p1[k];
1960 gel(E,j) = mulis(gel(p2,k),ex[i]);
1961 }
1962 }
1963 else
1964 {
1965 avma = av1;
1966 gel(P,j) = fx;
1967 gel(E,j) = utoipos(ex[i]); j++;
1968 }
1969 }
1970 setlg(P,j);
1971 setlg(E,j); return mkmat2(P, E);
1972 }
1973
1974 GEN
factorpadic4(GEN f,GEN p,long prec)1975 factorpadic4(GEN f,GEN p,long prec)
1976 {
1977 pari_sp av = avma;
1978 GEN y, P, ppow, lead, lt;
1979 long i, l, pr, n = degpol(f);
1980 int reverse = 0;
1981
1982 if (typ(f)!=t_POL) pari_err(notpoler,"factorpadic");
1983 if (typ(p)!=t_INT) pari_err(arither1);
1984 if (gcmp0(f)) pari_err(zeropoler,"factorpadic");
1985 if (prec <= 0) pari_err(talker,"non-positive precision in factorpadic");
1986 if (n == 0) return trivfact();
1987
1988 f = QpX_to_ZX(f); (void)Z_pvalrem(leading_term(f), p, <);
1989 f = pnormalize(f, p, prec, n-1, &lead, &pr, &reverse);
1990 y = ZX_monic_factorpadic(f, p, pr);
1991 P = gel(y,1); l = lg(P);
1992 if (lead)
1993 for (i=1; i<l; i++) gel(P,i) = primpart( RgX_unscale(gel(P,i), lead) );
1994 ppow = powiu(p,prec);
1995 for (i=1; i<l; i++)
1996 {
1997 GEN t = gel(P,i);
1998 if (reverse) t = normalizepol(polrecip_i(t));
1999 gel(P,i) = ZX_to_ZpX_normalized(t,p,ppow,prec);
2000 }
2001 if (!gcmp1(lt)) gel(P,1) = gmul(gel(P,1), lt);
2002 return gerepilecopy(av, sort_factor(y, cmp_padic));
2003 }
2004
2005 GEN
factorpadic0(GEN f,GEN p,long r,long flag)2006 factorpadic0(GEN f,GEN p,long r,long flag)
2007 {
2008 switch(flag)
2009 {
2010 case 0: return factorpadic4(f,p,r);
2011 case 1: return factorpadic2(f,p,r);
2012 default: pari_err(flagerr,"factorpadic");
2013 }
2014 return NULL; /* not reached */
2015 }
2016
2017 /*******************************************************************/
2018 /* */
2019 /* FACTORIZATION IN F_q */
2020 /* */
2021 /*******************************************************************/
2022 static GEN spec_FqXQ_pow(GEN x, GEN S, GEN T, GEN p);
2023
2024 static GEN
to_Fq(GEN x,GEN T,GEN p)2025 to_Fq(GEN x, GEN T, GEN p)
2026 {
2027 long i, lx, tx = typ(x);
2028 GEN y;
2029
2030 if (tx == t_INT)
2031 y = mkintmod(x,p);
2032 else
2033 {
2034 if (tx != t_POL) pari_err(typeer,"to_Fq");
2035 lx = lg(x);
2036 y = cgetg(lx,t_POL); y[1] = x[1];
2037 for (i=2; i<lx; i++) gel(y,i) = mkintmod(gel(x,i), p);
2038 }
2039 return mkpolmod(y, T);
2040 }
2041
2042 static GEN
to_Fq_pol(GEN x,GEN T,GEN p)2043 to_Fq_pol(GEN x, GEN T, GEN p)
2044 {
2045 long i, lx, tx = typ(x);
2046 if (tx != t_POL) pari_err(typeer,"to_Fq_pol");
2047 lx = lg(x);
2048 for (i=2; i<lx; i++) gel(x,i) = to_Fq(gel(x,i),T,p);
2049 return x;
2050 }
2051
2052 static GEN
to_Fq_fact(GEN P,GEN E,GEN T,GEN p,pari_sp av)2053 to_Fq_fact(GEN P, GEN E, GEN T, GEN p, pari_sp av)
2054 {
2055 GEN y = cgetg(3,t_MAT), u, v;
2056 long j, l = lg(P), nbf = lg(P);
2057
2058 u = cgetg(nbf,t_COL); gel(y,1) = u;
2059 v = cgetg(nbf,t_COL); gel(y,2) = v;
2060 for (j=1; j<l; j++)
2061 {
2062 gel(u,j) = simplify_i(gel(P,j)); /* may contain pols of degree 0 */
2063 gel(v,j) = utoi((ulong)E[j]);
2064 }
2065 y = gerepilecopy(av, y); u = gel(y,1);
2066 p = icopy(p);
2067 T = FpX_to_mod(T, p);
2068 for (j=1; j<nbf; j++) gel(u,j) = to_Fq_pol(gel(u,j), T,p);
2069 return y;
2070 }
2071
2072 /* split into r factors of degree d */
2073 static void
FqX_split(GEN * t,long d,GEN q,GEN S,GEN T,GEN p)2074 FqX_split(GEN *t, long d, GEN q, GEN S, GEN T, GEN p)
2075 {
2076 long l, v, is2, cnt, dt = degpol(*t), dT = degpol(T);
2077 pari_sp av;
2078 GEN w,w0;
2079
2080 if (dt == d) return;
2081 v = varn(*t);
2082 if (DEBUGLEVEL > 6) (void)timer2();
2083 av = avma; is2 = equaliu(p, 2);
2084 for(cnt = 1;;cnt++, avma = av)
2085 { /* splits *t with probability ~ 1 - 2^(1-r) */
2086 w = w0 = FqX_rand(dt,v, T,p);
2087 if (degpol(w) <= 0) continue;
2088 for (l=1; l<d; l++) /* sum_{0<i<d} w^(q^i), result in (F_q)^r */
2089 w = gadd(w0, spec_FqXQ_pow(w, S, T, p));
2090 w = FpXQX_red(w, T,p);
2091 if (is2)
2092 {
2093 w0 = w;
2094 for (l=1; l<dT; l++) /* sum_{0<i<k} w^(2^i), result in (F_2)^r */
2095 {
2096 w = FqX_rem(FqX_sqr(w,T,p), *t, T,p);
2097 w = FpXX_red(gadd(w0,w), p);
2098 }
2099 }
2100 else
2101 {
2102 w = FpXQYQ_pow(w, shifti(q,-1), *t, T, p);
2103 /* w in {-1,0,1}^r */
2104 if (degpol(w) <= 0) continue;
2105 gel(w,2) = gadd(gel(w,2), gen_1);
2106 }
2107 w = FqX_gcd(*t,w, T,p); l = degpol(w);
2108 if (l && l != dt) break;
2109 }
2110 w = gerepileupto(av,w);
2111 if (DEBUGLEVEL > 6)
2112 fprintferr("[FqX_split] splitting time: %ld (%ld trials)\n",timer2(),cnt);
2113 l /= d; t[l] = FqX_div(*t,w, T,p); *t = w;
2114 FqX_split(t+l,d,q,S,T,p);
2115 FqX_split(t ,d,q,S,T,p);
2116 }
2117
2118 /* to "compare" (real) scalars and t_INTMODs */
2119 static int
cmp_coeff(GEN x,GEN y)2120 cmp_coeff(GEN x, GEN y)
2121 {
2122 if (typ(x) == t_INTMOD) x = gel(x,2);
2123 if (typ(y) == t_INTMOD) y = gel(y,2);
2124 return gcmp(x,y);
2125 }
2126
2127 int
cmp_pol(GEN x,GEN y)2128 cmp_pol(GEN x, GEN y)
2129 {
2130 long fx[3], fy[3];
2131 long i,lx,ly;
2132 int fl;
2133 if (typ(x) == t_POLMOD) x = gel(x,2);
2134 if (typ(y) == t_POLMOD) y = gel(y,2);
2135 if (typ(x) == t_POL) lx = lg(x); else { lx = 3; gel(fx,2) = x; x = fx; }
2136 if (typ(y) == t_POL) ly = lg(y); else { ly = 3; gel(fy,2) = y; y = fy; }
2137 if (lx > ly) return 1;
2138 if (lx < ly) return -1;
2139 for (i=lx-1; i>1; i--)
2140 if ((fl = cmp_coeff(gel(x,i), gel(y,i)))) return fl;
2141 return 0;
2142 }
2143
2144 /*******************************************************************/
2145 /* */
2146 /* FACTOR USING TRAGER'S TRICK */
2147 /* */
2148 /*******************************************************************/
2149 /* Factor polynomial a on the number field defined by polynomial T */
2150 GEN
polfnf(GEN a,GEN T)2151 polfnf(GEN a, GEN T)
2152 {
2153 pari_sp av = avma;
2154 GEN x0, P, E, u, G, fa, n, unt, dent, A;
2155 long lx, i, k, e;
2156 int sqfree, tmonic;
2157
2158 if (typ(a)!=t_POL || typ(T)!=t_POL) pari_err(typeer,"polfnf");
2159 if (gcmp0(a)) return gcopy(a);
2160 A = lift(fix_relative_pol(T,a,0));
2161 P = content(A); if (!gcmp1(P)) A = gdiv(A, P);
2162 T = primpart(T);
2163 tmonic = is_pm1(leading_term(T));
2164
2165 dent = tmonic? indexpartial(T, NULL): ZX_disc(T);
2166 unt = mkpolmod(gen_1,T);
2167 G = nfgcd(A,derivpol(A), T, dent);
2168 sqfree = gcmp1(G);
2169 u = sqfree? A: Q_primpart(RgXQX_div(A, G, T));
2170 k = 0; n = ZY_ZXY_rnfequation(T, u, &k);
2171 if (DEBUGLEVEL>4) fprintferr("polfnf: choosing k = %ld\n",k);
2172 if (!sqfree)
2173 {
2174 G = poleval(G, gadd(pol_x[varn(A)], gmulsg(k, pol_x[varn(T)])));
2175 G = ZY_ZXY_resultant(T, G);
2176 }
2177 /* n guaranteed to be squarefree */
2178 fa = ZX_DDF(n,0); lx = lg(fa);
2179 P = cgetg(lx,t_COL);
2180 E = cgetg(lx,t_COL);
2181 if (lx == 2)
2182 { /* P^k, k irreducible */
2183 gel(P,1) = gmul(unt,u);
2184 gel(E,1) = utoipos(degpol(A) / degpol(u));
2185 return gerepilecopy(av, mkmat2(P,E));
2186 }
2187 x0 = gadd(pol_x[varn(A)], gmulsg(-k, mkpolmod(pol_x[varn(T)], T)));
2188 for (i=lx-1; i>0; i--)
2189 {
2190 GEN f = gel(fa,i), F = lift_intern(poleval(f, x0));
2191 if (!tmonic) F = primpart(F);
2192 F = nfgcd(u, F, T, dent);
2193 if (typ(F) != t_POL || degpol(F) == 0)
2194 pari_err(talker,"reducible modulus in factornf");
2195 e = 1;
2196 if (!sqfree)
2197 {
2198 while (poldvd(G,f, &G)) e++;
2199 if (degpol(G) == 0) sqfree = 1;
2200 }
2201 gel(P,i) = gdiv(gmul(unt,F), leading_term(F));
2202 gel(E,i) = utoipos(e);
2203 }
2204 return gerepilecopy(av, sort_factor(mkmat2(P,E), cmp_pol));
2205 }
2206
2207 static GEN
FqX_frob_deflate(GEN f,GEN T,GEN p)2208 FqX_frob_deflate(GEN f, GEN T, GEN p)
2209 {
2210 GEN F = poldeflate_i(f, itos(p)), frobinv = powiu(p, degpol(T)-1);
2211 long i, l = lg(F);
2212 for (i=2; i<l; i++) gel(F,i) = Fq_pow(gel(F,i), frobinv, T,p);
2213 return F;
2214 }
2215 /* Factor _sqfree_ polynomial a on the finite field Fp[X]/(T) */
2216 static GEN
FqX_split_Trager(GEN A,GEN T,GEN p)2217 FqX_split_Trager(GEN A, GEN T, GEN p)
2218 {
2219 GEN x0, P, u, fa, n;
2220 long lx, i, k;
2221
2222 u = A;
2223 n = NULL;
2224 for (k = 0; cmpui(k, p) < 0; k++)
2225 {
2226 GEN U = poleval(u, gadd(pol_x[varn(A)], gmulsg(k, pol_x[varn(T)])));
2227 n = FpY_FpXY_resultant(T, U, p);
2228 if (FpX_is_squarefree(n, p)) break;
2229 n = NULL;
2230 }
2231 if (!n) return NULL;
2232 if (DEBUGLEVEL>4) fprintferr("FqX_split_Trager: choosing k = %ld\n",k);
2233 /* n guaranteed to be squarefree */
2234 fa = FpX_factor(n, p); fa = gel(fa,1); lx = lg(fa);
2235 P = cgetg(lx,t_COL);
2236 if (lx == 2)
2237 { /* P^k, k irreducible */
2238 gel(P,1) = u;
2239 return P;
2240 }
2241 x0 = gadd(pol_x[varn(A)], gmulsg(-k, mkpolmod(pol_x[varn(T)], T)));
2242 for (i=lx-1; i>1; i--)
2243 {
2244 GEN f = gel(fa,i), F = lift_intern(poleval(f, x0));
2245 F = FqX_gcd(u, F, T, p);
2246 if (typ(F) != t_POL || degpol(F) == 0)
2247 pari_err(talker,"reducible modulus in factornf");
2248 u = FqX_div(u, F, T, p);
2249 gel(P,i) = F;
2250 }
2251 gel(P,1) = u; return P;
2252 }
2253
2254
2255 /* assume n = deg(u) > 1, X over FqX */
2256 /* return S = [ X^q, X^2q, ... X^(n-1)q ] mod u (in Fq[X]) in Kronecker form */
2257 static GEN
init_spec_FqXQ_pow(GEN X,GEN q,GEN u,GEN T,GEN p)2258 init_spec_FqXQ_pow(GEN X, GEN q, GEN u, GEN T, GEN p)
2259 {
2260 long i, n = degpol(u);
2261 GEN x, S = cgetg(n, t_VEC);
2262
2263 if (n == 1) return S;
2264 x = FpXQYQ_pow(X, q, u, T, p);
2265 gel(S,1) = x;
2266 if ((degpol(x)<<1) < degpol(T)) {
2267 for (i=2; i < n; i++)
2268 gel(S,i) = FqX_rem(gmul(gel(S,i-1), x), u, T,p);
2269 } else {
2270 for (i=2; i < n; i++)
2271 gel(S,i) = (i&1)? FqX_rem(gmul(gel(S,i-1), x), u, T,p)
2272 : FqX_rem(gsqr(gel(S,i>>1)), u, T,p);
2273 }
2274 for (i=1; i < n; i++) gel(S,i) = to_Kronecker(gel(S,i), T);
2275 return S;
2276 }
2277
2278 /* compute x^q, S is as above */
2279 static GEN
spec_FqXQ_pow(GEN x,GEN S,GEN T,GEN p)2280 spec_FqXQ_pow(GEN x, GEN S, GEN T, GEN p)
2281 {
2282 pari_sp av = avma, lim = stack_lim(av, 1);
2283 GEN x0 = x+2, z = gel(x0,0);
2284 long i, dx = degpol(x);
2285
2286 for (i = 1; i <= dx; i++)
2287 {
2288 GEN d, c = gel(x0,i);
2289 if (gcmp0(c)) continue;
2290 d = gel(S,i); if (!gcmp1(c)) d = gmul(c,d);
2291 z = gadd(z, d);
2292 if (low_stack(lim, stack_lim(av,1)))
2293 {
2294 if(DEBUGMEM>1) pari_warn(warnmem,"spec_FqXQ_pow");
2295 z = gerepileupto(av, z);
2296 }
2297 }
2298 z = FpXQX_from_Kronecker(z, T, p);
2299 setvarn(z, varn(x)); return gerepileupto(av, z);
2300 }
2301
2302 static long
isabsolutepol(GEN f)2303 isabsolutepol(GEN f)
2304 {
2305 long i, l = lg(f);
2306 for(i=2; i<l; i++)
2307 {
2308 GEN c = gel(f,i);
2309 if (typ(c) == t_POL && degpol(c) > 0) return 0;
2310 }
2311 return 1;
2312 }
2313
2314 typedef struct {
2315 GEN S, L, Xq;
2316 GEN q; /* p^deg(T) */
2317 GEN p, T; /* split mod (p, T(X)) */
2318 } FqX_split_t;
2319
2320 static void
add(GEN z,GEN g,long d)2321 add(GEN z, GEN g, long d) { appendL(z, mkvec2(utoipos(d), g)); }
2322 /* return number of roots */
2323 long
FqX_split_deg1(GEN * pz,GEN u,GEN q,GEN T,GEN p)2324 FqX_split_deg1(GEN *pz, GEN u, GEN q, GEN T, GEN p)
2325 {
2326 long dg, N = degpol(u);
2327 GEN v, S, g, X, z = cget1(N+1, t_VEC);
2328
2329 *pz = z;
2330 if (N == 1) return 1;
2331 v = X = pol_x[varn(u)];
2332 S = init_spec_FqXQ_pow(X, q, u, T, p);
2333 appendL(z, S);
2334 v = spec_FqXQ_pow(v, S, T, p);
2335 g = FqX_gcd(gsub(v,X),u, T,p);
2336 dg = degpol(g);
2337 if (dg > 0) add(z, g, dg);
2338 return dg;
2339 }
2340
2341 /* return number of factors */
2342 long
FqX_split_by_degree(GEN * pz,GEN u,GEN q,GEN T,GEN p)2343 FqX_split_by_degree(GEN *pz, GEN u, GEN q, GEN T, GEN p)
2344 {
2345 long nb = 0, d, dg, N = degpol(u);
2346 GEN v, S, g, X, z = cget1(N+1, t_VEC);
2347
2348 *pz = z;
2349 if (N == 1) return 1;
2350 v = X = pol_x[varn(u)];
2351 S = init_spec_FqXQ_pow(X, q, u, T, p);
2352 appendL(z, S);
2353 for (d=1; d <= N>>1; d++)
2354 {
2355 v = spec_FqXQ_pow(v, S, T, p);
2356 g = FqX_gcd(gsub(v,X),u, T,p);
2357 dg = degpol(g); if (dg <= 0) continue;
2358 /* all factors of g have degree d */
2359 add(z, g, dg / d); nb += dg / d;
2360 N -= dg;
2361 if (N)
2362 {
2363 u = FqX_div(u,g, T,p);
2364 v = FqX_rem(v,u, T,p);
2365 }
2366 }
2367 if (N) { add(z, u, 1); nb++; }
2368 return nb;
2369 }
2370
2371 static GEN
FqX_split_equal(GEN L,GEN S,GEN T,GEN p)2372 FqX_split_equal(GEN L, GEN S, GEN T, GEN p)
2373 {
2374 long n = itos(gel(L,1));
2375 GEN u = gel(L,2), z = cgetg(n + 1, t_VEC);
2376 gel(z,1) = u;
2377 FqX_split((GEN*)(z+1), degpol(u) / n, powiu(p, degpol(T)), S, T, p);
2378 return z;
2379 }
2380 GEN
FqX_split_roots(GEN z,GEN T,GEN p,GEN pol)2381 FqX_split_roots(GEN z, GEN T, GEN p, GEN pol)
2382 {
2383 GEN S = gel(z,1), L = gel(z,2), rep = FqX_split_equal(L, S, T, p);
2384 if (pol) rep = shallowconcat(rep, FqX_div(pol, gel(L,2), T,p));
2385 return rep;
2386 }
2387 GEN
FqX_split_all(GEN z,GEN T,GEN p)2388 FqX_split_all(GEN z, GEN T, GEN p)
2389 {
2390 GEN S = gel(z,1), rep = cgetg(1, t_VEC);
2391 long i, l = lg(z);
2392 for (i = 2; i < l; i++)
2393 rep = shallowconcat(rep, FqX_split_equal(gel(z,i), S, T, p));
2394 return rep;
2395 }
2396
2397 static long
FqX_sqf_split(GEN * t0,GEN q,GEN T,GEN p)2398 FqX_sqf_split(GEN *t0, GEN q, GEN T, GEN p)
2399 {
2400 GEN *t = t0, u = *t, v, S, g, X;
2401 long d, dg, N = degpol(u);
2402
2403 if (N == 1) return 1;
2404 v = X = pol_x[varn(u)];
2405 S = init_spec_FqXQ_pow(X, q, u, T, p);
2406 for (d=1; d <= N>>1; d++)
2407 {
2408 v = spec_FqXQ_pow(v, S, T, p);
2409 g = FqX_gcd(gsub(v,X),u, T,p);
2410 dg = degpol(g); if (dg <= 0) continue;
2411
2412 /* all factors of g have degree d */
2413 *t = g;
2414 FqX_split(t, d, q, S, T, p);
2415 t += dg / d;
2416 N -= dg;
2417 if (N)
2418 {
2419 u = FqX_div(u,g, T,p);
2420 v = FqX_rem(v,u, T,p);
2421 }
2422 }
2423 if (N) *t++ = u;
2424 return t - t0;
2425 }
2426
2427 /* not memory-clean */
2428 /* TODO: provide a public and clean FpX_factorff */
2429 static GEN
FpX_factorff(GEN P,GEN l,GEN Q)2430 FpX_factorff(GEN P,GEN l, GEN Q)
2431 {
2432 GEN V,E, F = FpX_factor(P,l);
2433 long lfact = 1, nmax = lgpol(P), n = lg(gel(F,1));
2434 long i;
2435 V = cgetg(nmax,t_VEC);
2436 E = cgetg(nmax,t_VECSMALL);
2437 for(i=1;i<n;i++)
2438 {
2439 GEN R = FpX_factorff_irred(gmael(F,1,i),Q,l);
2440 long j, r = lg(R);
2441 for (j=1;j<r;j++)
2442 {
2443 V[lfact] = R[j];
2444 E[lfact] = mael(F,2,i); lfact++;
2445 }
2446 }
2447 setlg(V,lfact);
2448 setlg(E,lfact); return sort_factor(mkvec2(V,E), cmp_pol);
2449 }
2450
2451 static GEN
FqX_factor_i(GEN f,GEN T,GEN p)2452 FqX_factor_i(GEN f, GEN T, GEN p)
2453 {
2454 long pg, j, k, d, e, N, nbfact, pk;
2455 GEN E, f2, f3, df1, df2, g1, u, q, *t;
2456
2457 if (!signe(f)) pari_err(zeropoler,"FqX_factor");
2458 d = degpol(f); if (!d) return trivfact();
2459 T = FpX_normalize(T, p);
2460 f = FqX_normalize(f, T, p);
2461 if (isabsolutepol(f)) return FpX_factorff(simplify_i(f), p, T);
2462
2463 pg = itos_or_0(p);
2464 df2 = NULL; /* gcc -Wall */
2465 t = (GEN*)cgetg(d+1,t_VEC);
2466 E = cgetg(d+1, t_VECSMALL);
2467
2468 q = powiu(p, degpol(T));
2469 e = nbfact = 1;
2470 pk = 1;
2471 f3 = NULL;
2472 df1 = FqX_deriv(f, T, p);
2473 for(;;)
2474 {
2475 long nb0;
2476 while (gcmp0(df1))
2477 { /* needs d >= p: pg = 0 can't happen */
2478 pk *= pg; e = pk;
2479 f = FqX_frob_deflate(f, T, p);
2480 df1 = FqX_deriv(f, T, p); f3 = NULL;
2481 }
2482 f2 = f3? f3: FqX_gcd(f,df1, T,p);
2483 if (!degpol(f2)) u = f;
2484 else
2485 {
2486 g1 = FqX_div(f,f2, T,p);
2487 df2 = FqX_deriv(f2, T,p);
2488 if (gcmp0(df2)) { u = g1; f3 = f2; }
2489 else
2490 {
2491 f3 = FqX_gcd(f2,df2, T,p);
2492 u = degpol(f3)? FqX_div(f2, f3, T,p): f2;
2493 u = FqX_div(g1, u, T,p);
2494 }
2495 }
2496 /* u is square-free (product of irreducibles of multiplicity e) */
2497 nb0 = nbfact; N = degpol(u);
2498 t[nbfact] = FqX_normalize(u, T,p);
2499 if (N) {
2500 nb0 = nbfact;
2501 t[nbfact] = FqX_normalize(u, T,p);
2502 if (N == 1) nbfact++;
2503 else
2504 {
2505 #if 0
2506 nbfact += FqX_split_Berlekamp(t+nbfact, q, T, p);
2507 #else
2508 GEN P = FqX_split_Trager(t[nbfact], T, p);
2509 if (P) {
2510 for (j = 1; j < lg(P); j++) t[nbfact++] = gel(P,j);
2511 } else {
2512 if (DEBUGLEVEL) pari_warn(warner, "FqX_split_Trager failed!");
2513 nbfact += FqX_sqf_split(t+nbfact, q, T, p);
2514 }
2515 #endif
2516 }
2517 for (j = nb0; j < nbfact; j++) E[j] = e;
2518 }
2519 if (!degpol(f2)) break;
2520 f = f2; df1 = df2; e += pk;
2521 }
2522
2523 for (j=1; j<nbfact; j++)
2524 {
2525 t[j] = FqX_normalize(gel(t,j), T,p);
2526 for (k=1; k<j; k++)
2527 if (gequal(t[j],t[k]))
2528 {
2529 E[k] += E[j]; nbfact--;
2530 E[j] = E[nbfact];
2531 t[j] = t[nbfact]; break;
2532 }
2533 }
2534 setlg(t, nbfact);
2535 setlg(E, nbfact); return sort_factor(mkvec2((GEN)t, E), cmp_pol);
2536 }
2537 GEN
factorff(GEN f,GEN p,GEN T)2538 factorff(GEN f, GEN p, GEN T)
2539 {
2540 pari_sp av;
2541 long v;
2542 GEN z;
2543
2544 if (typ(T)!=t_POL || typ(f)!=t_POL || typ(p)!=t_INT) pari_err(typeer,"factorff");
2545 v = varn(T);
2546 if (varncmp(v, varn(f)) <= 0)
2547 pari_err(talker,"polynomial variable must have higher priority in factorff");
2548 T = RgX_to_FpX(T, p); av = avma;
2549 z = FqX_factor_i(RgX_to_FqX(f, T,p), T, p);
2550 return to_Fq_fact(gel(z,1),gel(z,2), T,p,av);
2551 }
2552 /* factorization of x modulo (T,p). Assume x already reduced */
2553 GEN
FqX_factor(GEN x,GEN T,GEN p)2554 FqX_factor(GEN x, GEN T, GEN p)
2555 {
2556 pari_sp av = avma;
2557 if (!T) return FpX_factor(x, p);
2558 return gerepilecopy(av, FqX_factor_i(x, T, p));
2559 }
2560 /* See also: Isomorphisms between finite field and relative
2561 * factorization in polarit3.c */
2562
2563 /*******************************************************************/
2564 /* */
2565 /* COMPLEX ROOTS */
2566 /* */
2567 /*******************************************************************/
2568 static GEN laguer(GEN pol,long N,GEN y0,long EPS,long PREC);
2569 GEN zrhqr(GEN a,long PREC);
2570
2571 GEN
rootsold(GEN x,long prec)2572 rootsold(GEN x, long prec)
2573 {
2574 long i, j, f, real, exact, fr, deg, ln;
2575 pari_sp av=avma, av0, av1, av2, av3;
2576 long exc,expmin,m,deg0,k,ti,h,ii,e;
2577 GEN y,xc,xd0,xd,xdabs,p1,p2,p3,p4,p5,p6,p7;
2578 GEN p11,p12,p1r,p1i,pa,pax,pb,pp,pq,ps, pi;
2579
2580 if (typ(x)!=t_POL) pari_err(typeer,"rootsold");
2581 deg0 = degpol(x); expmin = 12 - bit_accuracy(prec);
2582 if (!signe(x)) pari_err(zeropoler,"rootsold");
2583 y = cgetg(deg0+1,t_COL); if (!deg0) return y;
2584 for (i=1; i<=deg0; i++)
2585 {
2586 p1 = cgetc(prec); gel(y,i) = p1;
2587 for (j=3; j<prec; j++) (gel(p1,2))[j] = (gel(p1,1))[j] = 0;
2588 }
2589 real=1; exact=1;
2590 for (i=2; i<=deg0+2; i++)
2591 {
2592 ti = typ(x[i]);
2593 if (ti==t_REAL) exact = 0;
2594 else if (ti==t_QUAD)
2595 {
2596 p2 = gmael3(x,i,1,2);
2597 if (gsigne(p2) > 0) real = 0;
2598 } else if (ti != t_INT && ti != t_FRAC) real = 0;
2599 }
2600 av1 = avma;
2601 k = polvaluation_inexact(x, &pax);
2602 for (i = 1; i <= k; i++) gaffsg(0,gel(y,i));
2603 if (k == deg0) return y;
2604
2605 pi = mppi(DEFAULTPREC);
2606 p2 = mkcomplex(pi, divrs(pi,10)); /* Pi * (1+I/10) */
2607 p11 = cgetg(4,t_POL); p11[1] = x[1];
2608 gel(p11,3) = gen_1;
2609
2610 p12 = cgetg(5,t_POL); p12[1] = x[1];
2611 gel(p12,4) = gen_1;
2612
2613 xd0 = derivpol(pax); pa = pax;
2614 pq = NULL; /* for lint */
2615 if (exact) { pp = ggcd(pax,xd0); h = degpol(pp); if (h) pq = RgX_div(pax,pp); }
2616 else{ pp = gen_1; h = 0; }
2617 m = 0;
2618 while (k != deg0)
2619 {
2620 m++;
2621 if (h)
2622 {
2623 pa = pp; pb = pq; pp = ggcd(pa,derivpol(pa)); h = degpol(pp);
2624 pq = h? RgX_div(pa,pp): pa;
2625 ps = RgX_div(pb,pq);
2626 }
2627 else ps = pa;
2628 deg = degpol(ps); if (!deg) continue;
2629
2630 /* roots of exact order m */
2631 e = gexpo(ps) - gexpo(leading_term(ps));
2632 if (e < 0) e = 0; if (ps!=pax) xd0 = derivpol(ps);
2633 xdabs = cgetg(deg+2,t_POL); xdabs[1] = xd0[1];
2634 for (i=2; i<deg+2; i++)
2635 {
2636 av3 = avma; p3 = gel(xd0,i);
2637 gel(xdabs,i) = gerepileupto(av3, gadd(gabs(real_i(p3),prec),
2638 gabs(imag_i(p3),prec)));
2639 }
2640 av0 = avma; xc = gcopy(ps); xd = gcopy(xd0); av2 = avma;
2641 for (i=1; i<=deg; i++)
2642 {
2643 if (i == deg)
2644 {
2645 p1 = (GEN)y[k+m*i];
2646 gdivz(gneg_i(gel(xc,2)),gel(xc,3), p1);
2647 p1r = gel(p1,1);
2648 p1i = gel(p1,2);
2649 }
2650 else
2651 {
2652 p3 = gshift(p2,e);
2653 p4 = poleval(xc,p3);
2654 p5 = gnorm(p4);
2655 exc = 0;
2656 while (exc >= -20)
2657 {
2658 p7 = gneg_i(gdiv(p4, poleval(xd,p3)));
2659 av3 = avma;
2660 exc = gcmp0(p5)? -32: expo(gnorm(p7))-expo(gnorm(p3));
2661 avma = av3;
2662 for (j=1; j<=10; j++)
2663 {
2664 GEN p8, p9, p10;
2665 p8 = gadd(p3,p7);
2666 p9 = poleval(xc,p8);
2667 p10= gnorm(p9);
2668 if (exc < -20 || cmprr(p10,p5) < 0)
2669 {
2670 GEN *gptr[3];
2671 p3=p8; p4=p9; p5=p10;
2672 gptr[0]=&p3; gptr[1]=&p4; gptr[2]=&p5;
2673 gerepilemanysp(av2,av3,gptr,3);
2674 break;
2675 }
2676 gshiftz(p7,-2,p7); avma = av3;
2677 }
2678 if (j > 10)
2679 {
2680 if (DEBUGLEVEL)
2681 pari_warn(warner,"too many iterations in rootsold(): using roots2()");
2682 avma = av; return roots2(x,prec);
2683 }
2684 }
2685 p1 = (GEN)y[k+m*i];
2686 p1r = gel(p1,1); setlg(p1r, 3);
2687 p1i = gel(p1,2); setlg(p1i, 3); gaffect(p3, p1); avma = av2;
2688 for (ln = 4; ln <= prec; ln = (ln<<1)-2)
2689 {
2690 setlg(p1r,ln); if (gcmp0(p1r)) gel(p1,1) = gen_0;
2691 setlg(p1i,ln); if (gcmp0(p1i)) gel(p1,2) = gen_0;
2692 p6 = gadd(p1, gneg_i(gdiv(poleval(xc,p1), poleval(xd,p1))));
2693 gel(p1,1) = p1r;
2694 gel(p1,2) = p1i; gaffect(p6, p1); avma = av2;
2695 }
2696 }
2697 setlg(p1r,prec);
2698 setlg(p1i,prec); p7 = gcopy(p1);
2699 p1r = gel(p7,1); setlg(p1r,prec+1);
2700 p1i = gel(p7,2); setlg(p1i,prec+1);
2701 for (ii=1; ii<=5; ii++)
2702 {
2703 if (typ(p7) == t_COMPLEX)
2704 {
2705 if (gcmp0(gel(p7,1))) gel(p7,1) = gen_0;
2706 if (gcmp0(gel(p7,2))) gel(p7,2) = gen_0;
2707 }
2708 p7 = gadd(p7, gneg_i(gdiv(poleval(ps,p7), poleval(xd0,p7))));
2709 }
2710 gaffect(p7, p1);
2711 p6 = gdiv(poleval(ps,p7), poleval(xdabs,gabs(p7,prec)));
2712 if (gexpo(p6) >= expmin)
2713 {
2714 if (DEBUGLEVEL) pari_warn(warner,"error in rootsold(): using roots2()");
2715 avma = av; return roots2(x,prec);
2716 }
2717 avma = av2;
2718 if (expo(p1[2]) < expmin && real)
2719 {
2720 gaffect(gen_0, gel(p1,2));
2721 for (j=1; j<m; j++) gaffect(p1, (GEN)y[k+(i-1)*m+j]);
2722 gel(p11,2) = gneg(gel(p1,1));
2723 xc = gerepileupto(av0, RgX_div(xc,p11));
2724 }
2725 else
2726 {
2727 for (j=1; j<m; j++) gaffect(p1, (GEN)y[k+(i-1)*m+j]);
2728 if (real)
2729 {
2730 p1 = gconj(p1);
2731 for (j=1; j<=m; j++) gaffect(p1, (GEN)y[k+i*m+j]);
2732 i++;
2733 gel(p12,2) = gnorm(p1);
2734 gel(p12,3) = gmulsg(-2,gel(p1,1));
2735 xc = gerepileupto(av0, RgX_div(xc,p12));
2736 }
2737 else
2738 {
2739 gel(p11,2) = gneg(p1);
2740 xc = gerepileupto(av0, RgX_div(xc,p11));
2741 }
2742 }
2743 xd = derivpol(xc); av2 = avma;
2744 }
2745 k += deg*m;
2746 }
2747 avma = av1;
2748 for (j=2; j<=deg0; j++)
2749 {
2750 p1 = gel(y,j);
2751 fr = !gcmp0(gel(p1,2));
2752 for (k=j-1; k>=1; k--)
2753 {
2754 p2 = gel(y,k);
2755 f = !gcmp0(gel(p2,2));
2756 if (!f && fr) break;
2757 if (f == fr && gcmp(gel(p2,1),gel(p1,1)) <= 0) break;
2758 y[k+1] = y[k];
2759 }
2760 gel(y,k+1) = p1;
2761 }
2762 return y;
2763 }
2764
2765 GEN
roots2(GEN pol,long PREC)2766 roots2(GEN pol,long PREC)
2767 {
2768 pari_sp av = avma;
2769 long N,flagexactpol,flagrealpol,flagrealrac,ti,i,j;
2770 long nbpol, k, multiqol, deg, nbroot, fr, f, EPS;
2771 pari_sp av1;
2772 GEN p1,p2,rr,qol,qolbis,x,b,c,ad,v, ex, factors;
2773
2774 if (typ(pol)!=t_POL) pari_err(typeer,"roots2");
2775 if (!signe(pol)) pari_err(zeropoler,"roots2");
2776 N=degpol(pol);
2777 if (!N) return cgetg(1,t_COL);
2778 if (N==1)
2779 {
2780 p1 = gmul(real_1(PREC),gel(pol,3));
2781 p2 = gneg_i(gdiv(gel(pol,2),p1));
2782 return gerepilecopy(av,p2);
2783 }
2784 EPS = 12 - bit_accuracy(PREC); /* 2^EPS is "zero" */
2785 flagrealpol = flagexactpol = 1;
2786 for (i=2; i<=N+2; i++)
2787 {
2788 c = gel(pol,i);
2789 switch (typ(c)) {
2790 case t_INT: case t_FRAC: break;
2791
2792 case t_REAL: flagexactpol = 0; break;
2793
2794 case t_COMPLEX: flagexactpol = flagrealpol = 0; break;
2795
2796 case t_QUAD: flagexactpol = 0;
2797 if (gsigne(gmael(c,1,2)) > 0) flagrealpol = 0;
2798 break;
2799 default: pari_err(typeer, "roots2");
2800 }
2801 }
2802 rr=cgetg(N+1,t_COL);
2803 for (i=1; i<=N; i++)
2804 {
2805 p1 = cgetc(PREC); gel(rr,i) = p1;
2806 for (j=3; j<PREC; j++) mael(p1,2,j) = mael(p1,1,j) = 0;
2807 }
2808 if (flagexactpol) { pol = Q_primpart(pol); factors = ZX_squff(pol, &ex); }
2809 else
2810 {
2811 factors = mkcol(pol);
2812 ex = mkvecsmall(1);
2813 }
2814 nbpol = lg(ex)-1;
2815 nbroot= 0;
2816 for (k=1; k<=nbpol; k++)
2817 {
2818 av1=avma; qol = gel(factors,k); qolbis = gcopy(qol);
2819 multiqol = ex[k]; deg = degpol(qol);
2820 for (j=deg; j>=1; j--)
2821 {
2822 x = gen_0; flagrealrac = 0;
2823 if (j==1) x = gneg_i(gdiv(gel(qolbis,2),gel(qolbis,3)));
2824 else
2825 {
2826 x = laguer(qolbis,j,x,EPS,PREC);
2827 if (x == NULL) goto RLAB;
2828 }
2829 if (flagexactpol)
2830 {
2831 x = gprec_w(x, PREC+2);
2832 x = laguer(qol,deg,x, EPS-BITS_IN_LONG, PREC+1);
2833 }
2834 else
2835 x = laguer(qol,deg,x,EPS,PREC);
2836 if (x == NULL) goto RLAB;
2837
2838 if (typ(x)==t_COMPLEX && gexpo(imag_i(x)) <= gexpo(real_i(x)) + EPS+1)
2839 { gel(x,2) = gen_0; flagrealrac = 1; }
2840 else if (j==1 && flagrealpol)
2841 { gel(x,2) = gen_0; flagrealrac = 1; }
2842 else if (typ(x)!=t_COMPLEX) flagrealrac = 1;
2843
2844 for (i=1; i<=multiqol; i++) gaffect(x, gel(rr,nbroot+i));
2845 nbroot += multiqol;
2846 if (!flagrealpol || flagrealrac)
2847 {
2848 ad = new_chunk(j+1);
2849 for (i=0; i<=j; i++) ad[i] = qolbis[i+2];
2850 b = gel(ad,j);
2851 for (i=j-1; i>=0; i--)
2852 {
2853 c = gel(ad,i); gel(ad,i) = b;
2854 b = gadd(gmul(gel(rr,nbroot),b),c);
2855 }
2856 v = cgetg(j+1,t_VEC); for (i=1; i<=j; i++) v[i] = ad[j-i];
2857 qolbis = gtopoly(v,varn(qolbis));
2858 if (flagrealpol)
2859 for (i=2; i<=j+1; i++)
2860 if (typ(qolbis[i])==t_COMPLEX) gmael(qolbis,i,2)= gen_0;
2861 }
2862 else
2863 {
2864 ad = new_chunk(j-1);
2865 ad[j-2] = qolbis[j+2];
2866 p1 = gmulsg(2,real_i(gel(rr,nbroot)));
2867 p2 = gnorm(gel(rr,nbroot));
2868 gel(ad,j-3) = gadd(gel(qolbis,j+1), gmul(p1,gel(ad,j-2)));
2869 for (i=j-2; i>=2; i--)
2870 gel(ad,i-2) = gadd(gel(qolbis,i+2),
2871 gsub(gmul(p1,gel(ad,i-1)),
2872 gmul(p2,gel(ad,i))));
2873 v = cgetg(j,t_VEC); for (i=1; i<=j-1; i++) v[i] = ad[j-1-i];
2874 qolbis = gtopoly(v,varn(qolbis));
2875 for (i=2; i<=j; i++)
2876 if (typ(qolbis[i])==t_COMPLEX) gmael(qolbis,i,2)= gen_0;
2877 for (i=1; i<=multiqol; i++)
2878 gaffect(gconj(gel(rr,nbroot)), gel(rr,nbroot+i));
2879 nbroot+=multiqol; j--;
2880 }
2881 }
2882 avma=av1;
2883 }
2884 for (j=2; j<=N; j++)
2885 {
2886 x=gel(rr,j); if (gcmp0(gel(x,2))) fr=0; else fr=1;
2887 for (i=j-1; i>=1; i--)
2888 {
2889 if (gcmp0(gmael(rr,i,2))) f=0; else f=1;
2890 if (f<fr) break;
2891 if (f==fr && gcmp(real_i(gel(rr,i)),real_i(x)) <= 0) break;
2892 rr[i+1]=rr[i];
2893 }
2894 gel(rr,i+1) = x;
2895 }
2896 return gerepilecopy(av,rr);
2897
2898 RLAB:
2899 avma = av;
2900 for(i=2;i<=N+2;i++)
2901 {
2902 ti = typ(pol[i]);
2903 if (!is_intreal_t(ti)) pari_err(talker,"too many iterations in roots");
2904 }
2905 if (DEBUGLEVEL)
2906 {
2907 fprintferr("too many iterations in roots2() ( laguer() ):\n");
2908 fprintferr(" real coefficients polynomial, using zrhqr()\n");
2909 }
2910 return zrhqr(pol,PREC);
2911 }
2912
2913 #define MR 8
2914 #define MT 10
2915
2916 static GEN
laguer(GEN pol,long N,GEN y0,long EPS,long PREC)2917 laguer(GEN pol,long N,GEN y0,long EPS,long PREC)
2918 {
2919 long MAXIT, iter, j;
2920 pari_sp av = avma, av1;
2921 GEN rac,erre,I,x,abx,abp,abm,dx,x1,b,d,f,g,h,sq,gp,gm,g2,*ffrac;
2922
2923 MAXIT = MR*MT; rac = cgetc(PREC);
2924 av1 = avma;
2925 I = mkcomplex(gen_1,gen_1);
2926 ffrac = (GEN*)new_chunk(MR+1);
2927 ffrac[0] = dbltor(0.0);
2928 ffrac[1] = dbltor(0.5);
2929 ffrac[2] = dbltor(0.25);
2930 ffrac[3] = dbltor(0.75);
2931 ffrac[4] = dbltor(0.13);
2932 ffrac[5] = dbltor(0.38);
2933 ffrac[6] = dbltor(0.62);
2934 ffrac[7] = dbltor(0.88);
2935 ffrac[8] = dbltor(1.0);
2936 x=y0;
2937 for (iter=1; iter<=MAXIT; iter++)
2938 {
2939 b = gel(pol,N+2); d = f = gen_0;
2940 erre = QuickNormL1(b,PREC);
2941 abx = QuickNormL1(x,PREC);
2942 for (j=N-1; j>=0; j--)
2943 {
2944 f = gadd(gmul(x,f), d);
2945 d = gadd(gmul(x,d), b);
2946 b = gadd(gmul(x,b), gel(pol,j+2));
2947 erre = gadd(QuickNormL1(b,PREC), gmul(abx,erre));
2948 }
2949 erre = gmul2n(erre, EPS);
2950 if (gcmp(QuickNormL1(b,PREC),erre)<=0)
2951 {
2952 gaffect(x,rac); avma = av1; return rac;
2953 }
2954 g = gdiv(d,b);
2955 g2 = gsqr(g); h = gsub(g2, gmul2n(gdiv(f,b),1));
2956 sq = gsqrt(gmulsg(N-1,gsub(gmulsg(N,h),g2)),PREC);
2957 gp = gadd(g,sq); abp = gnorm(gp);
2958 gm = gsub(g,sq); abm = gnorm(gm);
2959 if (gcmp(abp,abm) < 0) gp = gm;
2960 if (gsigne(gmax(abp,abm)) > 0)
2961 dx = gdivsg(N,gp);
2962 else
2963 dx = gmul(gadd(gen_1,abx), gexp(gmulgs(I,iter),PREC));
2964 x1 = gsub(x,dx);
2965 if (gexpo(QuickNormL1(gsub(x,x1),PREC)) < EPS)
2966 {
2967 gaffect(x,rac); avma = av1; return rac;
2968 }
2969 if (iter%MT) x = gcopy(x1); else x = gsub(x, gmul(ffrac[iter/MT],dx));
2970 }
2971 avma = av; return NULL;
2972 }
2973 #undef MR
2974 #undef MT
2975 /***********************************************************************/
2976 /** **/
2977 /** ROOTS of a polynomial with REAL coeffs **/
2978 /** **/
2979 /***********************************************************************/
2980 #define RADIX 1L
2981 #define COF 0.95
2982
2983 /* x t_MAT in M_n(R) : compute a symmetric matrix with the same eigenvalues */
2984 static GEN
balanc(GEN x)2985 balanc(GEN x)
2986 {
2987 pari_sp av = avma;
2988 long last, i, j, sqrdx = (RADIX<<1), n = lg(x);
2989 GEN r, c, cofgen, a = shallowcopy(x);
2990
2991 last = 0; cofgen = dbltor(COF);
2992 while (!last)
2993 {
2994 last = 1;
2995 for (i=1; i<n; i++)
2996 {
2997 r = c = gen_0;
2998 for (j=1; j<n; j++)
2999 if (j!=i)
3000 {
3001 c = gadd(c, gabs(gcoeff(a,j,i),0));
3002 r = gadd(r, gabs(gcoeff(a,i,j),0));
3003 }
3004 if (!gcmp0(r) && !gcmp0(c))
3005 {
3006 GEN g, s = gmul(cofgen, gadd(c,r));
3007 long ex = 0;
3008 g = gmul2n(r,-RADIX); while (gcmp(c,g) < 0) {ex++; c=gmul2n(c, sqrdx);}
3009 g = gmul2n(r, RADIX); while (gcmp(c,g) > 0) {ex--; c=gmul2n(c,-sqrdx);}
3010 if (gcmp(gadd(c,r), gmul2n(s,ex)) < 0)
3011 {
3012 last = 0;
3013 for (j=1; j<n; j++) gcoeff(a,i,j) = gmul2n(gcoeff(a,i,j),-ex);
3014 for (j=1; j<n; j++) gcoeff(a,j,i) = gmul2n(gcoeff(a,j,i), ex);
3015 }
3016 }
3017 }
3018 }
3019 return gerepilecopy(av, a);
3020 }
3021
3022 #define SIGN(a,b) ((b)>=0.0 ? fabs(a) : -fabs(a))
3023 /* find the eigenvalues of the symmetric matrix mat */
3024 static GEN
hqr(GEN mat)3025 hqr(GEN mat)
3026 {
3027 long n = lg(mat)-1, N, m, l, k, j, i, mmin, flj, flk;
3028 double **a, p, q, r, s, t, u, v, w, x, y, z, anorm, *wr, *wi;
3029 const double eps = 0.000001;
3030 GEN eig;
3031
3032 init_dalloc();
3033 a = (double**)stackmalloc(sizeof(double*)*(n+1));
3034 for (i=1; i<=n; i++) a[i] = (double*)stackmalloc(sizeof(double)*(n+1));
3035 for (j=1; j<=n; j++)
3036 for (i=1; i<=n; i++) a[i][j] = gtodouble(gcoeff(mat,i,j));
3037 wr = (double*)stackmalloc(sizeof(double)*(n+1));
3038 wi = (double*)stackmalloc(sizeof(double)*(n+1));
3039
3040 anorm = fabs(a[1][1]);
3041 for (i=2; i<=n; i++) for (j=(i-1); j<=n; j++) anorm += fabs(a[i][j]);
3042 N = n; t = 0.;
3043 p = q = r = 0.; /* -Wall */
3044 if (DEBUGLEVEL>3) { fprintferr("* Finding eigenvalues\n"); flusherr(); }
3045 while (N>=1)
3046 {
3047 long its = 0;
3048 for(;;)
3049 {
3050 for (l=N; l>=2; l--)
3051 {
3052 s = fabs(a[l-1][l-1])+fabs(a[l][l]); if (s==0.) s = anorm;
3053 if (fabs(a[l][l-1])+s == s) break;
3054 }
3055 x = a[N][N];
3056 if (l == N){ wr[N] = x+t; wi[N] = 0.; N--; break; } /* OK */
3057
3058 y = a[N-1][N-1];
3059 w = a[N][N-1]*a[N-1][N];
3060 if (l == N-1)
3061 {
3062 p = 0.5*(y-x); q = p*p+w; z = sqrt(fabs(q)); x += t;
3063 if (q >= 0. || fabs(q) <= eps)
3064 {
3065 z = p + SIGN(z,p);
3066 wr[N-1] = wr[N] = x+z;
3067 if (fabs(z)>eps) wr[N] = x-w/z;
3068 wi[N-1] = wi[N] = 0.;
3069 }
3070 else { wr[N-1] = wr[N]= x+p; wi[N-1] = -z; wi[N] = z; }
3071 N -= 2; break; /* OK */
3072 }
3073
3074 if (its==30) pari_err(talker,"too many iterations in hqr");
3075 if (its==10 || its==20)
3076 {
3077 t += x; for (i=1; i<=N; i++) a[i][i] -= x;
3078 s = fabs(a[N][N-1]) + fabs(a[N-1][N-2]);
3079 y = x = 0.75*s;
3080 w = -0.4375*s*s;
3081 }
3082 its++;
3083 for (m=N-2; m>=l; m--)
3084 {
3085 z = a[m][m]; r = x-z; s = y-z;
3086 p = (r*s-w)/a[m+1][m]+a[m][m+1];
3087 q = a[m+1][m+1]-z-r-s;
3088 r = a[m+2][m+1];
3089 s = fabs(p)+fabs(q)+fabs(r); p/=s; q/=s; r/=s;
3090 if (m==l) break;
3091 u = fabs(a[m][m-1])*(fabs(q)+fabs(r));
3092 v = fabs(p) * (fabs(a[m-1][m-1])+fabs(z)+fabs(a[m+1][m+1]));
3093 if (u+v==v) break;
3094 }
3095 for (i=m+2; i<=N; i++){ a[i][i-2]=0.; if (i!=m+2) a[i][i-3]=0.; }
3096 for (k=m; k<=N-1; k++)
3097 {
3098 if (k!=m)
3099 {
3100 p = a[k][k-1]; q = a[k+1][k-1];
3101 r = (k != N-1)? a[k+2][k-1]: 0.;
3102 x = fabs(p)+fabs(q)+fabs(r);
3103 if (x != 0.) { p/=x; q/=x; r/=x; }
3104 }
3105 s = SIGN(sqrt(p*p+q*q+r*r),p);
3106 if (s == 0.) continue;
3107
3108 if (k==m)
3109 { if (l!=m) a[k][k-1] = -a[k][k-1]; }
3110 else
3111 a[k][k-1] = -s*x;
3112 p+=s; x=p/s; y=q/s; z=r/s; q/=p; r/=p;
3113 for (j=k; j<=N; j++)
3114 {
3115 p = a[k][j]+q*a[k+1][j];
3116 if (k != N-1) { p+=r*a[k+2][j]; a[k+2][j]-=p*z; }
3117 a[k+1][j] -= p*y; a[k][j] -= p*x;
3118 }
3119 mmin = (N < k+3)? N: k+3;
3120 for (i=l; i<=mmin; i++)
3121 {
3122 p = x*a[i][k]+y*a[i][k+1];
3123 if (k != N-1) { p+=z*a[i][k+2]; a[i][k+2]-=p*r; }
3124 a[i][k+1] -= p*q; a[i][k] -= p;
3125 }
3126 }
3127 }
3128 }
3129 for (j=2; j<=n; j++) /* ordering the roots */
3130 {
3131 x = wr[j];
3132 y = wi[j]; flj = (y != 0.);
3133 for (k=j-1; k>=1; k--)
3134 {
3135 flk = (wi[k] != 0.);
3136 if (!flk && flj) break;
3137 if (flk == flj && wr[k] <= x) break;
3138 wr[k+1] = wr[k];
3139 wi[k+1] = wi[k];
3140 }
3141 wr[k+1] = x;
3142 wi[k+1] = y;
3143 }
3144 if (DEBUGLEVEL>3) { fprintferr("* Eigenvalues computed\n"); flusherr(); }
3145 eig = cgetg(n+1,t_COL);
3146 for (i=1; i<=n; i++)
3147 gel(eig,i) = (wi[i] == 0.)? dbltor(wr[i])
3148 : mkcomplex(dbltor(wr[i]), dbltor(wi[i]));
3149 return eig;
3150 }
3151
3152 /* a t_POL in R[X], squarefree: give the roots of the polynomial a (real roots
3153 * first) in increasing order of their real parts. */
3154 GEN
zrhqr(GEN a,long prec)3155 zrhqr(GEN a, long prec)
3156 {
3157 pari_sp av = avma;
3158 long i, prec2, n = degpol(a), ex = -bit_accuracy(prec);
3159 GEN aa, b, rt, rr, x, dx, y, newval, oldval;
3160
3161 rt = hqr(balanc(assmat(a)));
3162 prec2 = 2*prec; /* polishing the roots */
3163 aa = gprec_w(a, prec2);
3164 b = derivpol(aa); rr = cgetg(n+1,t_COL);
3165 for (i=1; i<=n; i++)
3166 {
3167 x = gprec_w(gel(rt,i), prec2);
3168 for (oldval=NULL;; oldval=newval, x=y)
3169 { /* Newton iteration */
3170 dx = poleval(b,x);
3171 if (gexpo(dx) < ex)
3172 pari_err(talker,"polynomial has probably multiple roots in zrhqr");
3173 y = gsub(x, gdiv(poleval(aa,x),dx));
3174 newval = gabs(poleval(aa,y),prec2);
3175 if (gexpo(newval) < ex || (oldval && gcmp(newval,oldval) > 0)) break;
3176 }
3177 if (DEBUGLEVEL>3) fprintferr("%ld ",i);
3178 gel(rr,i) = gtofp(y, prec);
3179 }
3180 if (DEBUGLEVEL>3) { fprintferr("\npolished roots = %Z",rr); flusherr(); }
3181 return gerepilecopy(av, rr);
3182 }
3183