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