Annotation of OpenXM_contrib/pari-2.2/src/modules/nffactor.c, Revision 1.1.1.1
1.1 noro 1: /* $Id: nffactor.c,v 1.29 2001/10/01 12:11:33 karim Exp $
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: /* POLYNOMIAL FACTORIZATION IN A NUMBER FIELD */
19: /* */
20: /*******************************************************************/
21: #include "pari.h"
22: #include "parinf.h"
23:
24: extern GEN hensel_lift_fact(GEN pol, GEN fact, GEN T, GEN p, GEN pev, long e);
25: extern GEN nf_get_T2(GEN base, GEN polr);
26: extern GEN nfreducemodpr_i(GEN x, GEN prh);
27: extern GEN sort_factor(GEN y, int (*cmp)(GEN,GEN));
28: extern GEN pidealprimeinv(GEN nf, GEN x);
29:
30: static GEN nffactormod2(GEN nf,GEN pol,GEN pr);
31: static GEN nfmod_split2(GEN nf, GEN prhall, GEN polb, GEN v, GEN q);
32: static GEN nf_pol_mul(GEN nf,GEN pol1,GEN pol2);
33: static GEN nf_pol_divres(GEN nf,GEN pol1,GEN pol2, GEN *pr);
34: static GEN nf_pol_subres(GEN nf,GEN pol1,GEN pol2);
35: static GEN nfmod_pol_reduce(GEN nf,GEN prhall,GEN pol);
36: static GEN nfmod_pol_divres(GEN nf,GEN prhall,GEN pol1,GEN pol2, GEN *pr);
37: static GEN nfmod_pol_gcd(GEN nf,GEN prhall,GEN pol1,GEN pol2);
38: static GEN nf_bestlift(GEN id,GEN idinv,GEN den,GEN elt);
39: static GEN nf_pol_lift(GEN id,GEN idinv,GEN den,GEN pol);
40: static GEN nfsqff(GEN nf,GEN pol,long fl);
41: static int nf_combine_factors(GEN nf,long fxn,GEN psf,long dlim,long hint);
42:
43: typedef struct Nfcmbf /* for use in nfsqff */
44: {
45: GEN pol, h, hinv, fact, res, lt, den;
46: long nfact, nfactmod;
47: } Nfcmbf;
48: static Nfcmbf nfcmbf;
49:
50: static GEN
51: unifpol0(GEN nf,GEN pol,long flag)
52: {
53: static long n = 0;
54: static GEN vun = NULL;
55: GEN f = (GEN) nf[1];
56: long i = degpol(f), av;
57:
58: if (i != n)
59: {
60: n=i; if (vun) free(vun);
61: vun = (GEN) gpmalloc((n+1)*sizeof(long));
62: vun[0] = evaltyp(t_COL) | evallg(n+1); vun[1]=un;
63: for ( ; i>=2; i--) vun[i]=zero;
64: }
65:
66: av = avma;
67: switch(typ(pol))
68: {
69: case t_INT: case t_FRAC: case t_RFRAC:
70: pol = gmul(pol,vun); break;
71:
72: case t_POL:
73: pol = gmodulcp(pol,f); /* fall through */
74: case t_POLMOD:
75: pol = algtobasis(nf,pol);
76: }
77: if (flag) pol = basistoalg(nf, lift(pol));
78: return gerepileupto(av,pol);
79: }
80:
81: /* Let pol be a polynomial with coefficients in Z or nf (vectors or polymods)
82: * return the same polynomial with coefficients expressed:
83: * if flag=0: as vectors (on the integral basis).
84: * if flag=1: as polymods.
85: */
86: GEN
87: unifpol(GEN nf,GEN pol,long flag)
88: {
89: if (typ(pol)==t_POL && varn(pol) < varn(nf[1]))
90: {
91: long i, d = lgef(pol);
92: GEN p1 = pol;
93:
94: pol=cgetg(d,t_POL); pol[1]=p1[1];
95: for (i=2; i<d; i++)
96: pol[i] = (long) unifpol0(nf,(GEN) p1[i],flag);
97:
98: return pol;
99: }
100: return unifpol0(nf,(GEN) pol, flag);
101: }
102:
103: #if 0
104: /* return a monic polynomial of degree d with random coefficients in Z_nf */
105: static GEN
106: random_pol(GEN nf,long d)
107: {
108: long i,j, n = degpol(nf[1]);
109: GEN pl,p;
110:
111: pl=cgetg(d+3,t_POL);
112: for (i=2; i<d+2; i++)
113: {
114: p=cgetg(n+1,t_COL); pl[i]=(long)p;
115: for (j=1; j<=n; j++)
116: p[j] = lstoi(mymyrand()%101 - 50);
117: }
118: p=cgetg(n+1,t_COL); pl[i]=(long)p;
119: p[1]=un; for (i=2; i<=n; i++) p[i]=zero;
120:
121: pl[1] = evalsigne(1) | evallgef(d+3) | evalvarn(0);
122: return pl;
123: }
124: #endif
125:
126: /* multiplication of x by y */
127: static GEN
128: nf_pol_mul(GEN nf,GEN x,GEN y)
129: {
130: long tetpil,av=avma;
131: GEN res = gmul(unifpol(nf,x,1), unifpol(nf,y,1));
132:
133: tetpil = avma;
134: return gerepile(av,tetpil,unifpol(nf,res,0));
135: }
136:
137: /* compute x^2 in nf */
138: static GEN
139: nf_pol_sqr(GEN nf,GEN x)
140: {
141: long tetpil,av=avma;
142: GEN res = gsqr(unifpol(nf,x,1));
143:
144: tetpil = avma;
145: return gerepile(av,tetpil,unifpol(nf,res,0));
146: }
147:
148: /* reduce the coefficients of pol modulo prhall */
149: static GEN
150: nfmod_pol_reduce(GEN nf,GEN prhall,GEN pol)
151: {
152: long av=avma,tetpil,i;
153: GEN p1;
154:
155: if (typ(pol)!=t_POL) return nfreducemodpr(nf,pol,prhall);
156: pol=unifpol(nf,pol,0);
157:
158: tetpil=avma; i=lgef(pol);
159: p1=cgetg(i,t_POL); p1[1]=pol[1];
160: for (i--; i>=2; i--)
161: p1[i] = (long) nfreducemodpr(nf,(GEN)pol[i],prhall);
162: return gerepile(av,tetpil, normalizepol(p1));
163: }
164:
165: /* x^2 modulo prhall */
166: static GEN
167: nfmod_pol_sqr(GEN nf,GEN prhall,GEN x)
168: {
169: long av=avma,tetpil;
170: GEN px;
171:
172: px = nfmod_pol_reduce(nf,prhall,x);
173: px = unifpol(nf,lift(px),1);
174: px = unifpol(nf,nf_pol_sqr(nf,px),0);
175: tetpil=avma;
176: return gerepile(av,tetpil,nfmod_pol_reduce(nf,prhall,px));
177: }
178:
179: /* multiplication of x by y modulo prhall */
180: static GEN
181: nfmod_pol_mul(GEN nf,GEN prhall,GEN x,GEN y)
182: {
183: long av=avma,tetpil;
184: GEN px,py;
185:
186: px = nfmod_pol_reduce(nf,prhall,x); px = unifpol(nf,lift(px),1);
187: py = nfmod_pol_reduce(nf,prhall,y); py = unifpol(nf,lift(py),1);
188: px = unifpol(nf,nf_pol_mul(nf,px,py),0);
189: tetpil=avma;
190: return gerepile(av,tetpil,nfmod_pol_reduce(nf,prhall,px));
191: }
192:
193: /* Euclidan division of x by y */
194: static GEN
195: nf_pol_divres(GEN nf,GEN x,GEN y,GEN *pr)
196: {
197: long av = avma,tetpil;
198: GEN nq = poldivres(unifpol(nf,x,1),unifpol(nf,y,1),pr);
199: GEN *gptr[2];
200:
201: tetpil=avma; nq=unifpol(nf,nq,0);
202: if (pr) *pr = unifpol(nf,*pr,0);
203: gptr[0]=&nq; gptr[1]=pr;
204: gerepilemanysp(av,tetpil,gptr,pr ? 2:1);
205: return nq;
206: }
207:
208: /* Euclidan division of x by y modulo prhall */
209: static GEN
210: nfmod_pol_divres(GEN nf,GEN prhall,GEN x,GEN y, GEN *pr)
211: {
212: long av=avma,dx,dy,dz,i,j,k,l,n,tetpil;
213: GEN z,p1,p3,px,py;
214:
215: py = nfmod_pol_reduce(nf,prhall,y);
216: if (gcmp0(py))
217: err(talker, "division by zero in nfmod_pol_divres");
218:
219: tetpil=avma;
220: px=nfmod_pol_reduce(nf,prhall,x);
221: dx=degpol(px); dy=degpol(py); dz=dx-dy;
222: if (dx<dy)
223: {
224: GEN vzero;
225:
226: if (pr) *pr = gerepile(av,tetpil,px);
227: else avma = av;
228:
229: n=degpol(nf[1]);
230: vzero = cgetg(n+1,t_COL);
231: n=degpol(nf[1]);
232: for (i=1; i<=n; i++) vzero[i]=zero;
233:
234: z=cgetg(3,t_POL); z[2]=(long)vzero;
235: z[1]=evallgef(2) | evalvarn(varn(px));
236: return z;
237: }
238: p1 = NULL; /* gcc -Wall */
239:
240: z=cgetg(dz+3,t_POL); z[1]=evalsigne(1) | evallgef(3+dz);
241: setvarn(z,varn(px));
242: z[dz+2] = (long) element_divmodpr(nf,(GEN)px[dx+2],(GEN)py[dy+2],prhall);
243: for (i=dx-1; i>=dy; --i)
244: {
245: l=avma; p1=(GEN)px[i+2];
246: for (j=i-dy+1; j<=i && j<=dz; j++)
247: p1 = gsub(p1, element_mul(nf,(GEN)z[j+2],(GEN)py[i-j+2]));
248: p1 = nfreducemodpr(nf,p1,prhall);
249: tetpil=avma; p3=element_divmodpr(nf,p1,(GEN)py[dy+2],prhall);
250: z[i-dy+2]=lpile(l,tetpil,p3);
251: z[i-dy+2]=(long)nfreducemodpr(nf,(GEN)z[i-dy+2],prhall);
252: }
253: l=avma;
254: for (i=dy-1; i>=0; --i)
255: {
256: l=avma; p1=((GEN)px[i+2]);
257: for (j=0; j<=i && j<=dz; j++)
258: p1 = gsub(p1, element_mul(nf,(GEN)z[j+2],(GEN)py[i-j+2]));
259: p1 = gerepileupto(l, nfreducemodpr(nf,p1,prhall));
260: if (!gcmp0(p1)) break;
261: }
262:
263: if (!pr) { avma = l; return z; }
264:
265: if (i<0)
266: {
267: avma=l;
268: p3 = cgetg(3,t_POL); p3[2]=zero;
269: p3[1] = evallgef(2) | evalvarn(varn(px));
270: *pr=p3; return z;
271: }
272:
273: p3=cgetg(i+3,t_POL);
274: p3[1]=evalsigne(1) | evallgef(i+3) | evalvarn(varn(px));
275: p3[i+2]=(long)nfreducemodpr(nf,p1,prhall);
276: for (k=i-1; k>=0; --k)
277: {
278: l=avma; p1=((GEN)px[k+2]);
279: for (j=0; j<=k && j<=dz; j++)
280: p1 = gsub(p1, element_mul(nf,(GEN)z[j+2],(GEN)py[k-j+2]));
281: p3[k+2]=lpileupto(l,nfreducemodpr(nf,p1,prhall));
282: }
283: *pr=p3; return z;
284: }
285:
286: /* GCD of x and y */
287: static GEN
288: nf_pol_subres(GEN nf,GEN x,GEN y)
289: {
290: long av=avma,tetpil;
291: GEN s = srgcd(unifpol(nf,x,1), unifpol(nf,y,1));
292:
293: tetpil=avma; return gerepile(av,tetpil,unifpol(nf,s,1));
294: }
295:
296: /* GCD of x and y modulo prhall */
297: static GEN
298: nfmod_pol_gcd(GEN nf,GEN prhall,GEN x,GEN y)
299: {
300: long av=avma;
301: GEN p1,p2;
302:
303: if (lgef(x)<lgef(y)) { p1=y; y=x; x=p1; }
304: p1=nfmod_pol_reduce(nf,prhall,x);
305: p2=nfmod_pol_reduce(nf,prhall,y);
306: while (!isexactzero(p2))
307: {
308: GEN p3;
309:
310: nfmod_pol_divres(nf,prhall,p1,p2,&p3);
311: p1=p2; p2=p3;
312: }
313: return gerepileupto(av,p1);
314: }
315:
316: /* return pol^e modulo prhall and the polynomial pmod */
317: static GEN
318: nfmod_pol_pow(GEN nf,GEN prhall,GEN pmod,GEN pol,GEN e)
319: {
320: long i, av = avma, n = degpol(nf[1]);
321: GEN p1,p2,vun;
322:
323: vun=cgetg(n+1,t_COL); vun[1]=un; for (i=2; i<=n; i++) vun[i]=zero;
324: p1=gcopy(polun[varn(pol)]); p1[2]=(long)vun;
325: if (gcmp0(e)) return p1;
326:
327: p2=nfmod_pol_reduce(nf,prhall,pol);
328: for(;;)
329: {
330: if (!vali(e))
331: {
332: p1=nfmod_pol_mul(nf,prhall,p1,p2);
333: nfmod_pol_divres(nf,prhall,p1,pmod,&p1);
334: }
335: if (gcmp1(e)) break;
336:
337: e=shifti(e,-1);
338: p2=nfmod_pol_sqr(nf,prhall,p2);
339: nfmod_pol_divres(nf,prhall,p2,pmod,&p2);
340: }
341: return gerepileupto(av,p1);
342: }
343:
344: static long
345: isdivbyprime(GEN nf, GEN x, GEN pr)
346: {
347: GEN elt, p = (GEN)pr[1], tau = (GEN)pr[5];
348:
349: elt = element_mul(nf, x, tau);
350: if (divise(content(elt), p)) return 1;
351:
352: return 0;
353: }
354:
355: /* return the factor of nf.pol modulo p corresponding to pr */
356: static GEN
357: localpol(GEN nf, GEN pr)
358: {
359: long i, l;
360: GEN fct, pol = (GEN)nf[1], p = (GEN)pr[1];
361:
362: fct = lift((GEN)factmod(pol, p)[1]);
363: l = lg(fct) - 1;
364: for (i = 1; i <= l; i++)
365: if (isdivbyprime(nf, (GEN)fct[i], pr)) return (GEN)fct[i];
366:
367: err(talker, "cannot find a suitable factor in localpol");
368: return NULL; /* not reached */
369: }
370:
371: /* factorization of x modulo pr */
372: static GEN
373: nffactormod0(GEN nf, GEN x, GEN pr)
374: {
375: long av = avma, j, l, vx = varn(x), vn;
376: GEN rep, pol, xrd, prh, p1;
377:
378: nf=checknf(nf);
379: vn = varn((GEN)nf[1]);
380: if (typ(x)!=t_POL) err(typeer,"nffactormod");
381: if (vx >= vn)
382: err(talker,"polynomial variable must have highest priority in nffactormod");
383:
384: prh = nfmodprinit(nf, pr);
385:
386: if (divise((GEN)nf[4], (GEN)pr[1]))
387: return gerepileupto(av, nffactormod2(nf,x,pr));
388:
389: xrd = nfmod_pol_reduce(nf, prh, x);
390: if (gcmp1((GEN)pr[4]))
391: {
392: pol = gun; /* dummy */
393: for (j = 2; j < lg(xrd); j++)
394: xrd[j] = mael(xrd, j, 1);
395: rep = factmod(xrd, (GEN)pr[1]);
396: rep = lift(rep);
397: }
398: else
399: {
400: pol = localpol(nf, pr);
401: xrd = lift(unifpol(nf, xrd, 1));
402: p1 = gun;
403: for (j = 2; j < lg(xrd); j++)
404: {
405: xrd[j] = lmod((GEN)xrd[j], pol);
406: p1 = mpppcm(p1, denom(content((GEN)xrd[j])));
407: }
408: rep = factmod9(gmul(xrd, p1), (GEN)pr[1], pol);
409: rep = lift(lift(rep));
410: }
411:
412: l = lg((GEN)rep[1]);
413: for (j = 1; j < l; j++)
414: mael(rep, 1, j) = (long)unifpol(nf, gmael(rep, 1, j), 1);
415:
416: return gerepilecopy(av, rep);
417: }
418:
419: GEN
420: nffactormod(GEN nf, GEN x, GEN pr)
421: {
422: return nffactormod0(nf, x, pr);
423: }
424:
425: extern GEN trivfact(void);
426:
427: /* factorization of x modulo pr */
428: GEN
429: nffactormod2(GEN nf,GEN pol,GEN pr)
430: {
431: long av = avma, tetpil,lb,nbfact,psim,N,n,i,j,k,d,e,vf,r,kk;
432: GEN y,ex,*t,f1,f2,f3,df1,g1,polb,pold,polu,vker;
433: GEN Q,f,x,u,v,v2,v3,vz,q,vun,vzero,prhall;
434:
435: nf=checknf(nf);
436: if (typ(pol)!=t_POL) err(typeer,"nffactormod");
437: if (varn(pol) >= varn(nf[1]))
438: err(talker,"polynomial variable must have highest priority in nffactormod");
439:
440: prhall=nfmodprinit(nf,pr); n=degpol(nf[1]);
441: vun = gscalcol_i(gun, n);
442: vzero = gscalcol_i(gzero, n);
443: q = vker = NULL; /* gcc -Wall */
444:
445: f=unifpol(nf,pol,0); f=nfmod_pol_reduce(nf,prhall,f);
446: if (!signe(f)) err(zeropoler,"nffactormod");
447: d=degpol(f); vf=varn(f);
448: if (d == 0) return trivfact();
449: t = (GEN*)new_chunk(d+1); ex = cgetg(d+1, t_VECSMALL);
450: x=gcopy(polx[vf]); x[3]=(long)vun; x[2]=(long)vzero;
451: if (d<=1) { nbfact=2; t[1] = f; ex[1]=1; }
452: else
453: {
454: q = (GEN)pr[1]; psim = VERYBIGINT;
455: if (!is_bigint(q)) psim = itos(q);
456: /* psim has an effect only when p is small. If too big, set it to a huge
457: * number (i.e ignore it) to avoid an error in itos on next line.
458: */
459: q=gpuigs(q, itos((GEN)pr[4]));
460: f1=f; e=1; nbfact=1;
461: while (lgef(f1)>3)
462: {
463: df1=deriv(f1,vf); f2=nfmod_pol_gcd(nf,prhall,f1,df1);
464: g1=nfmod_pol_divres(nf,prhall,f1,f2,NULL); k=0;
465: while (lgef(g1)>3)
466: {
467: k++;
468: if (k%psim == 0)
469: {
470: k++; f2=nfmod_pol_divres(nf,prhall,f2,g1,NULL);
471: }
472: f3=nfmod_pol_gcd(nf,prhall,f2,g1);
473: u = nfmod_pol_divres(nf,prhall,g1,f3,NULL);
474: f2= nfmod_pol_divres(nf,prhall,f2,f3,NULL);
475: g1=f3;
476: if (lgef(u)>3)
477: {
478: N=degpol(u); Q=cgetg(N+1,t_MAT);
479: v3=cgetg(N+1,t_COL); Q[1]=(long)v3;
480: v3[1]=(long)vun; for (i=2; i<=N; i++) v3[i]=(long)vzero;
481:
482: v2 = v = nfmod_pol_pow(nf,prhall,u,x,q);
483: for (j=2; j<=N; j++)
484: {
485: v3=cgetg(N+1,t_COL); Q[j]=(long)v3;
486: for (i=1; i<=lgef(v2)-2; i++) v3[i]=v2[i+1];
487: for (; i<=N; i++) v3[i]=(long)vzero;
488: if (j<N)
489: {
490: v2=nfmod_pol_mul(nf,prhall,v2,v);
491: nfmod_pol_divres(nf,prhall,v2,u,&v2);
492: }
493: }
494: for (i=1; i<=N; i++)
495: coeff(Q,i,i)=lsub((GEN)coeff(Q,i,i),vun);
496: v2=nfkermodpr(nf,Q,prhall); r=lg(v2)-1; t[nbfact]=gcopy(u); kk=1;
497: if (r>1)
498: {
499: vker=cgetg(r+1,t_COL);
500: for (i=1; i<=r; i++)
501: {
502: v3=cgetg(N+2,t_POL);
503: v3[1]=evalsigne(1)+evallgef(2+N); setvarn(v3,vf);
504: vker[i]=(long)v3; for (j=1; j<=N; j++) v3[j+1]=coeff(v2,j,i);
505: normalizepol(v3);
506: }
507: }
508: while (kk<r)
509: {
510: v=gcopy(polun[vf]); v[2]=(long)vzero;
511: for (i=1; i<=r; i++)
512: {
513: vz=cgetg(n+1,t_COL);
514: for (j=1; j<=n; j++)
515: vz[j] = lmodsi(mymyrand()>>8, q);
516: vz=nfreducemodpr(nf,vz,prhall);
517: v=gadd(v,nfmod_pol_mul(nf,prhall,vz,(GEN)vker[i]));
518: }
519: for (i=1; i<=kk && kk<r; i++)
520: {
521: polb=t[nbfact+i-1]; lb=lgef(polb);
522: if (lb>4)
523: {
524: if(psim==2)
525: polu=nfmod_split2(nf,prhall,polb,v,q);
526: else
527: {
528: polu=nfmod_pol_pow(nf,prhall,polb,v,shifti(q,-1));
529: polu=gsub(polu,vun);
530: }
531: pold=nfmod_pol_gcd(nf,prhall,polb,polu);
532: if (lgef(pold)>3 && lgef(pold)<lb)
533: {
534: t[nbfact+i-1]=pold; kk++;
535: t[nbfact+kk-1]=nfmod_pol_divres(nf,prhall,polb,pold,NULL);
536: }
537: }
538: }
539: }
540: for (i=nbfact; i<nbfact+r; i++) ex[i]=e*k;
541: nbfact+=r;
542: }
543: }
544: e*=psim; j=(degpol(f2))/psim+3; f1=cgetg(j,t_POL);
545: f1[1] = evalsigne(1) | evallgef(j) | evalvarn(vf);
546: for (i=2; i<j; i++)
547: f1[i]=(long)element_powmodpr(nf,(GEN)f2[psim*(i-2)+2],
548: gdiv(q,(GEN)pr[1]),prhall);
549: }
550: }
551: if (nbfact < 2)
552: err(talker, "%Z not a prime (nffactormod)", q);
553: for (j=1; j<nbfact; j++)
554: {
555: v = element_divmodpr(nf,vun,gmael(t,j,lgef(t[j])-1),prhall);
556: t[j] = unifpol(nf,nfmod_pol_mul(nf,prhall,v,(GEN)t[j]),1);
557: }
558:
559: tetpil=avma; y=cgetg(3,t_MAT);
560: u=cgetg(nbfact,t_COL); y[1]=(long)u;
561: v=cgetg(nbfact,t_COL); y[2]=(long)v;
562: for (j=1,k=0; j<nbfact; j++)
563: if (ex[j])
564: {
565: k++;
566: u[k]=lcopy((GEN)t[j]);
567: v[k]=lstoi(ex[j]);
568: }
569: return gerepile(av,tetpil,y);
570: }
571:
572: /* return pol + pol^2 + ... + pol^(q/2) modulo prhall and
573: the polynomial pmod */
574: static GEN
575: nfmod_split2(GEN nf,GEN prhall,GEN pmod,GEN pol,GEN exp)
576: {
577: long av = avma;
578: GEN p1,p2,q;
579:
580: if (cmpis(exp,2)<=0) return pol;
581: p2=p1=pol; q=shifti(exp,-1);
582: while (!gcmp1(q))
583: {
584: p2=nfmod_pol_sqr(nf,prhall,p2);
585: nfmod_pol_divres(nf,prhall,p2,pmod,&p2);
586: q=shifti(q,-1); p1=gadd(p1,p2);
587: }
588: return gerepileupto(av,p1);
589: }
590:
591: /* If p doesn't divide either a or b and has a divisor of degree 1, return it.
592: * Return NULL otherwise.
593: */
594: static GEN
595: p_ok(GEN nf, GEN p, GEN a)
596: {
597: long av,m,i;
598: GEN dec;
599:
600: if (divise(a,p)) return NULL;
601: av = avma; dec = primedec(nf,p); m=lg(dec);
602: for (i=1; i<m; i++)
603: {
604: GEN pr = (GEN)dec[i];
605: if (is_pm1(pr[4]))
606: return pr;
607: }
608: avma = av; return NULL;
609: }
610:
611: /* for each new prime ct--, if ct = 0, return NULL */
612: static GEN
613: choose_prime(GEN nf, GEN dk, GEN lim, long ct)
614: {
615: GEN p, pr;
616:
617: p = nextprime(lim);
618: for (;;)
619: {
620: if ((pr = p_ok(nf,p,dk))) break;
621: ct--;
622: if (!ct) return NULL;
623: p = nextprime(addis(p,2));
624: }
625:
626: return pr;
627: }
628:
629: /* test if the discriminant of polbase modulo some few primes
630: is non-zero. Return 1 if it is so (=> polbase is square-free)
631: and 0 otherwise (=> polbase may or may not be square-free) */
632: static int
633: is_sqf(GEN nf, GEN polbase)
634: {
635: GEN lt, pr, prh, p2, p;
636: long i, d = lgef(polbase), ct = 5;
637:
638: lt = (GEN)leading_term(polbase)[1];
639: p = stoi(101);
640:
641: while (ct > 0)
642: {
643: /* small primes tend to divide discriminants more often
644: than large ones so we look at primes >= 101 */
645: pr = choose_prime(nf,lt,p,30);
646: if (!pr) break;
647:
648: p=(GEN)pr[1];
649: prh=prime_to_ideal(nf,pr);
650:
651: p2=gcopy(polbase);
652: lt=mpinvmod(lt,p);
653:
654: for (i=2; i<d; i++)
655: p2[i] = nfreducemodpr_i(gmul(lt,(GEN)p2[i]), prh)[1];
656: p2 = normalizepol(p2);
657:
658: /* discriminant is non-zero => polynomial is square-free */
659: if (!gcmp0(p2) && !divise(discsr(p2),p)) { return 1; }
660:
661: ct--;
662: p=addis(p,1);
663: }
664:
665: return 0;
666: }
667:
668: /* rescale p in K[X] (coeffs in algtobasis form) --> primitive in O_K[X] */
669: GEN
670: nf_pol_to_int(GEN p, GEN *den)
671: {
672: long i, l = lgef(p);
673: GEN d = gun;
674: for (i=2; i<l; i++)
675: d = mpppcm(d,denom((GEN)p[i]));
676: if (! gcmp1(d)) p = gmul(p, d);
677: *den = d; return p;
678: }
679:
680: /* return the roots of pol in nf */
681: GEN
682: nfroots(GEN nf,GEN pol)
683: {
684: long av=avma,tetpil,d=lgef(pol),fl;
685: GEN p1,p2,polbase,polmod,den;
686:
687: p2=NULL; /* gcc -Wall */
688: nf=checknf(nf);
689: if (typ(pol)!=t_POL) err(talker,"not a polynomial in nfroots");
690: if (varn(pol) >= varn(nf[1]))
691: err(talker,"polynomial variable must have highest priority in nfroots");
692:
693: polbase=unifpol(nf,pol,0);
694:
695: if (d==3)
696: {
697: tetpil=avma; p1=cgetg(1,t_VEC);
698: return gerepile(av,tetpil,p1);
699: }
700:
701: if (d==4)
702: {
703: tetpil=avma; p1=cgetg(2,t_VEC);
704: p1[1] = (long)basistoalg(nf,gneg_i(
705: element_div(nf,(GEN)polbase[2],(GEN)polbase[3])));
706: return gerepile(av,tetpil,p1);
707: }
708:
709: p1=element_inv(nf,leading_term(polbase));
710: polbase=nf_pol_mul(nf,p1,polbase);
711:
712: polbase = nf_pol_to_int(polbase, &den);
713: polmod=unifpol(nf,polbase,1);
714:
715: if (DEBUGLEVEL>=4)
716: fprintferr("test if the polynomial is square-free\n");
717:
718: fl = is_sqf(nf, polbase);
719:
720: /* the polynomial may not be square-free ... */
721: if (!fl)
722: {
723: p1=derivpol(polmod);
724: p2=nf_pol_subres(nf,polmod,p1);
725: if (degpol(p2) == 0) fl = 1;
726: }
727:
728: if (!fl)
729: {
730: p1=element_inv(nf,leading_term(p2));
731: p2=nf_pol_mul(nf,p1,p2);
732: polmod=nf_pol_divres(nf,polmod,p2,NULL);
733:
734: p1=element_inv(nf,leading_term(polmod));
735: polmod=nf_pol_mul(nf,p1,polmod);
736:
737: polmod = nf_pol_to_int(polmod, &den);
738: polmod=unifpol(nf,polmod,1);
739: }
740:
741: p1 = nfsqff(nf,polmod,1);
742: tetpil=avma; return gerepile(av, tetpil, gen_sort(p1, 0, cmp_pol));
743: }
744:
745: /* return a minimal lift of elt modulo id */
746: static GEN
747: nf_bestlift(GEN id,GEN idinv,GEN den,GEN elt)
748: {
749: GEN u = gmul(idinv,elt);
750: long i, l = lg(u);
751: for (i=1; i<l; i++) u[i] = (long)gdivround((GEN)u[i], den);
752: return gsub(elt, gmul(id, u));
753: }
754:
755: /* return the lift of pol with coefficients of T2-norm <= C (if possible) */
756: static GEN
757: nf_pol_lift(GEN id,GEN idinv,GEN den,GEN pol)
758: {
759: long i, d = lgef(pol);
760: GEN x = cgetg(d,t_POL);
761: x[1] = pol[1];
762: for (i=2; i<d; i++)
763: x[i] = (long) nf_bestlift(id,idinv,den,(GEN)pol[i]);
764: return x;
765: }
766:
767: #if 0
768: /* return pol(elt) */
769: static GEN
770: nf_pol_eval(GEN nf,GEN pol,GEN elt)
771: {
772: long av=avma,tetpil,i;
773: GEN p1;
774:
775: i=lgef(pol)-1; if (i==2) return gcopy((GEN)pol[2]);
776:
777: p1=element_mul(nf,(GEN)pol[i],elt);
778: for (i--; i>=3; i--)
779: p1=element_mul(nf,elt,gadd((GEN)pol[i],p1));
780: tetpil=avma; return gerepile(av,tetpil,gadd(p1,(GEN)pol[2]));
781: }
782: #endif
783:
784: /* return the factorization of x in nf */
785: GEN
786: nffactor(GEN nf,GEN pol)
787: {
788: GEN y,p1,p2,den,p3,p4,quot,rep=cgetg(3,t_MAT);
789: long av = avma,tetpil,i,j,d,fl;
790:
791: if (DEBUGLEVEL >= 4) timer2();
792:
793: p3=NULL; /* gcc -Wall */
794: nf=checknf(nf);
795: if (typ(pol)!=t_POL) err(typeer,"nffactor");
796: if (varn(pol) >= varn(nf[1]))
797: err(talker,"polynomial variable must have highest priority in nffactor");
798:
799: d=lgef(pol);
800: if (d==3)
801: {
802: rep[1]=lgetg(1,t_COL);
803: rep[2]=lgetg(1,t_COL);
804: return rep;
805: }
806: if (d==4)
807: {
808: p1=cgetg(2,t_COL); rep[1]=(long)p1; p1[1]=lcopy(pol);
809: p1=cgetg(2,t_COL); rep[2]=(long)p1; p1[1]=un;
810: return rep;
811: }
812:
813: p1=element_inv(nf,leading_term(pol));
814: pol=nf_pol_mul(nf,p1,pol);
815:
816: p1=unifpol(nf,pol,0);
817: p1 = nf_pol_to_int(p1, &den);
818:
819: if (DEBUGLEVEL>=4)
820: fprintferr("test if the polynomial is square-free\n");
821:
822: fl = is_sqf(nf, p1);
823:
824: /* polynomial may not be square-free ... */
825: if (!fl)
826: {
827: p2=derivpol(p1);
828: p3=nf_pol_subres(nf,p1,p2);
829: if (degpol(p3) == 0) fl = 1;
830: }
831:
832: if (!fl)
833: {
834: p4=element_inv(nf,leading_term(p3));
835: p3=nf_pol_mul(nf,p4,p3);
836:
837: p2=nf_pol_divres(nf,p1,p3,NULL);
838: p4=element_inv(nf,leading_term(p2));
839: p2=nf_pol_mul(nf,p4,p2);
840:
841: p2 = nf_pol_to_int(p2, &den);
842:
843: p2=unifpol(nf,p2,1);
844: tetpil = avma; y = nfsqff(nf,p2,0);
845: i = nfcmbf.nfact;
846:
847: quot=nf_pol_divres(nf,p1,p2,NULL);
848: p3=(GEN)gpmalloc((i+1) * sizeof(long));
849: for (j=i; j>=1; j--)
850: {
851: GEN fact=(GEN)y[j], quo = quot, rem;
852: long e = 0;
853:
854: do
855: {
856: quo = nf_pol_divres(nf,quo,fact,&rem);
857: e++;
858: }
859: while (gcmp0(rem));
860: p3[j]=lstoi(e);
861: }
862: avma = (long)y; y = gerepile(av, tetpil, y);
863: p2=cgetg(i+1, t_COL); for (; i>=1; i--) p2[i]=lcopy((GEN)p3[i]);
864: free(p3);
865: }
866: else
867: {
868: tetpil=avma; y = gerepile(av,tetpil,nfsqff(nf,p1,0));
869: i = nfcmbf.nfact;
870: p2=cgetg(i+1, t_COL); for (; i>=1; i--) p2[i]=un;
871: }
872: if (DEBUGLEVEL>=4)
873: fprintferr("number of factor(s) found: %ld\n", nfcmbf.nfact);
874: rep[1]=(long)y;
875: rep[2]=(long)p2; return sort_factor(rep, cmp_pol);
876: }
877:
878: /* test if the matrix M is suitable */
879: static long
880: test_mat(GEN M, GEN p, GEN C2, long k)
881: {
882: long av = avma, i, N = lg(M);
883: GEN min, prod, L2, R;
884:
885: min = prod = gcoeff(M,1,1);
886: for (i = 2; i < N; i++)
887: {
888: L2 = gcoeff(M,i,i); prod = mpmul(prod,L2);
889: if (mpcmp(L2,min) < 0) min = L2;
890: }
891: R = mpmul(min, gpowgs(p, k<<1));
892: i = mpcmp(mpmul(C2,prod), R); avma = av;
893: return (i < 0);
894: }
895:
896: /* return the matrix corresponding to pr^e with R(pr^e) > C */
897: static GEN
898: T2_matrix_pow(GEN nf, GEN pre, GEN p, GEN C, long *kmax, long prec)
899: {
900: long av=avma,av1,lim, k = *kmax, N = degpol((GEN)nf[1]);
901: int tot_real = !signe(gmael(nf,2,2));
902: GEN p1,p2,p3,u,C2,T2, x = (GEN)nf[1];
903:
904: C2 = gdiv(gmul2n(C,2), absi((GEN)nf[3]));
905: p1 = gmul(pre,lllintpartial(pre)); av1 = avma;
906: T2 = tot_real? gmael(nf,5,4)
907: : nf_get_T2((GEN) nf[7], roots(x,prec));
908: p3 = qf_base_change(T2,p1,1);
909:
910: if (N <= 6 && test_mat(p3,p,C2,k))
911: {
912: avma = av1; return gerepileupto(av,p1);
913: }
914:
915: av1=avma; lim = stack_lim(av1,2);
916: for (;;)
917: {
918: if (DEBUGLEVEL>2) fprintferr("exponent: %ld\n",k);
919:
920: for(;;)
921: {
922: u = tot_real? lllgramall(p3,2,lll_IM) : lllgramintern(p3,2,1,prec);
923: if (u) break;
924:
925: prec=(prec<<1)-2;
926: if (DEBUGLEVEL > 1) err(warnprec,"nffactor[1]",prec);
927: T2 = nf_get_T2((GEN) nf[7], roots(x,prec));
928: p3 = qf_base_change(T2,p1,1);
929: }
930: if (DEBUGLEVEL>2) msgtimer("lllgram + base change");
931: p3 = qf_base_change(p3,u,1);
932: if (test_mat(p3,p,C2,k))
933: {
934: *kmax = k;
935: return gerepileupto(av,gmul(p1,u));
936: }
937:
938: /* we also need to increase the precision */
939: p2=shifti(gceil(mulsr(k, glog(p,DEFAULTPREC))),-1);
940: prec += (long)(itos(p2)*pariK1);
941: if (DEBUGLEVEL > 1) err(warnprec,"nffactor[2]",prec);
942: k = k<<1; p1 = idealmullll(nf,p1,p1);
943: if (low_stack(lim, stack_lim(av1,2)))
944: {
945: if (DEBUGMEM>1) err(warnmem,"T2_matrix_pow");
946: p1 = gerepileupto(av1,p1);
947: }
948: if (!tot_real) T2 = nf_get_T2((GEN) nf[7], roots(x,prec));
949: p3 = qf_base_change(T2,p1,1);
950: }
951: }
952:
953: /* return the factorization of the square-free polynomial x.
954: The coeff of x are in Z_nf and its leading term is a rational
955: integer. If fl = 1,return only the roots of x in nf */
956: static GEN
957: nfsqff(GEN nf,GEN pol, long fl)
958: {
959: long d=lgef(pol),i,k,m,n,av=avma,tetpil,newprec,prec,nbf=BIGINT,anbf,ct=5;
960: GEN p1,pr,p2,rep,k2,C,h,dk,dki,p,prh,p3,T2,polbase,fact,pk,ap,apr;
961: GEN polmod,polred,hinv,lt,minp,den,maxp=shifti(gun,32),run,aprh;
962:
963: if (DEBUGLEVEL>=4) msgtimer("square-free");
964:
965: dk=absi((GEN)nf[3]);
966: dki=mulii(dk,(GEN)nf[4]);
967: n=degpol(nf[1]);
968:
969: polbase = unifpol(nf,pol,0);
970: polmod = unifpol(nf,pol,1);
971: dki=mulii(dki,gnorm((GEN)polmod[d-1]));
972:
973: prec = DEFAULTPREC;
974: for (;;)
975: {
976: if (prec <= gprecision(nf))
977: T2 = gprec_w(gmael(nf,5,3), prec);
978: else
979: T2 = nf_get_T2((GEN)nf[7], roots((GEN)nf[1], prec));
980:
981: run=realun(prec);
982: p1=realzero(prec);
983: for (i=2; i<d; i++)
984: {
985: p2 = gmul(run, (GEN)polbase[i]);
986: p2 = qfeval(T2, p2);
987: p1 = addrr(p1, gdiv(p2, binome(stoi(d), i-2)));
988: if (signe(p1) < 0) break;
989: }
990:
991: if (signe(p1) > 0) break;
992:
993: prec = (prec<<1)-2;
994: if (DEBUGLEVEL > 1) err(warnprec, "nfsqff", prec);
995: }
996:
997: lt = (GEN)leading_term(polbase)[1];
998: p1 = gmul(p1, mulis(sqri(lt), n));
999: C = gpow(stoi(3), gadd(gmulsg(3, ghalf), stoi(d)), prec);
1000: C = gdiv(gmul(C, p1), gmulsg(d, mppi(prec)));
1001:
1002: if (DEBUGLEVEL>=4)
1003: fprintferr("the bound on the T2-norm of the coeff. is: %Z\n", C);
1004:
1005: /* the theoretical bound for the exponent should be:
1006: k2=gadd(glog(gdivgs(C,n),DEFAULTPREC), mulsr(n*(n-1), dbltor(0.347))); */
1007: k2=gadd(glog(gdivgs(C,n),DEFAULTPREC), mulsr(n*(n-1), dbltor(0.15)));
1008: k2=gmul2n(gmulgs(k2,n),-1);
1009:
1010: minp=gmin(gexp(gmul2n(k2,-6),BIGDEFAULTPREC), maxp);
1011: minp=gceil(minp);
1012:
1013: if (DEBUGLEVEL>=4)
1014: {
1015: fprintferr("lower bound for the prime numbers: %Z\n", minp);
1016: msgtimer("bounds computation");
1017: }
1018:
1019: p = rep = polred = NULL; /* gcc -Wall */
1020: pr=NULL;
1021: for (;;)
1022: {
1023: apr = choose_prime(nf,dki,minp, pr?30:0);
1024: if (!apr) break;
1025:
1026: ap=(GEN)apr[1];
1027: aprh=prime_to_ideal(nf,apr);
1028:
1029: polred=gcopy(polbase);
1030: lt=(GEN)leading_term(polbase)[1];
1031: lt=mpinvmod(lt,ap);
1032:
1033: for (i=2; i<d; i++)
1034: polred[i] = nfreducemodpr_i(gmul(lt,(GEN)polbase[i]), aprh)[1];
1035:
1036: if (!divise(discsr(polred),ap))
1037: {
1038: rep=(GEN)simplefactmod(polred,ap)[1];
1039: anbf=lg(rep)-1;
1040: ct--;
1041: if (anbf < nbf)
1042: {
1043: nbf=anbf;
1044: pr=gcopy(apr);
1045: p=gcopy(ap);
1046: if (DEBUGLEVEL>=4)
1047: {
1048: fprintferr("prime ideal considered: %Z\n", pr);
1049: fprintferr("number of irreducible factors: %ld\n", nbf);
1050: }
1051: if (nbf == 1) break;
1052: }
1053: }
1054: if (pr && !ct) break;
1055:
1056: minp = addis(ap,1);
1057: }
1058:
1059: k = itos(gceil(gdiv(k2,glog(p,BIGDEFAULTPREC))));
1060:
1061: if (DEBUGLEVEL>=4)
1062: {
1063: fprintferr("prime ideal chosen: %Z\n", pr);
1064: msgtimer("choice of the prime ideal");
1065: }
1066:
1067: if (lg(rep)==2)
1068: {
1069: if (fl) { avma=av; return cgetg(1,t_VEC); }
1070: rep=cgetg(2,t_VEC); rep[1]=lcopy(polmod);
1071: nfcmbf.nfact=1; return gerepileupto(av, rep);
1072: }
1073:
1074: p2=cgetr(DEFAULTPREC);
1075: affir(mulii(absi(dk),gpowgs(p,k)),p2);
1076: p2=shifti(gceil(mplog(p2)),-1);
1077:
1078: newprec = MEDDEFAULTPREC + (long)(itos(p2)*pariK1);
1079: if (DEBUGLEVEL>=4)
1080: fprintferr("new precision: %ld\n",newprec);
1081:
1082: prh = idealpows(nf,pr,k); m = k;
1083: h = T2_matrix_pow(nf,prh,p,C, &k, newprec);
1084: if (m != k) prh=idealpows(nf,pr,k);
1085:
1086: if (DEBUGLEVEL>=4)
1087: {
1088: fprintferr("a suitable exponent is: %ld\n",(long)k);
1089: msgtimer("computation of H");
1090: }
1091:
1092: pk = gcoeff(prh,1,1);
1093: lt=(GEN)leading_term(polbase)[1];
1094: lt=mpinvmod(lt,pk);
1095:
1096: polred[1] = polbase[1];
1097: for (i=2; i<d; i++)
1098: polred[i] = nfreducemodpr_i(gmul(lt,(GEN)polbase[i]), prh)[1];
1099:
1100: fact = lift_intern((GEN)factmod(polred,p)[1]);
1101: rep = hensel_lift_fact(polred,fact,NULL,p,pk,k);
1102:
1103: if (DEBUGLEVEL >= 4) msgtimer("computation of the p-adic factorization");
1104:
1105: den = det(h); /* |den| = NP^k */
1106: hinv= adj(h);
1107: lt=(GEN)leading_term(polbase)[1];
1108:
1109: if (fl)
1110: {
1111: long x_a[] = { evaltyp(t_POL)|m_evallg(4), 0,0,0 };
1112: GEN mlt = gneg_i(lt), rem;
1113: x_a[1] = polbase[1]; setlgef(x_a, 4);
1114: x_a[3] = un;
1115: p1=cgetg(lg(rep)+1,t_VEC);
1116: polbase = unifpol(nf,polbase,1);
1117: for (m=1,i=1; i<lg(rep); i++)
1118: {
1119: p2=(GEN)rep[i]; if (lgef(p2)!=4) break;
1120:
1121: p3 = algtobasis(nf, gmul(mlt,(GEN)p2[2]));
1122: p3 = nf_bestlift(h,hinv,den,p3);
1123: p3 = gdiv(p3,lt);
1124: x_a[2] = lneg(p3); /* check P(p3) == 0 */
1125: p2 = poldivres(polbase, unifpol(nf,x_a,1), &rem);
1126: if (!signe(rem)) { polbase = p2; p1[m++] = (long)p3; }
1127: }
1128: tetpil=avma; rep=cgetg(m,t_VEC);
1129: for (i=1; i<m; i++) rep[i]=(long)basistoalg(nf,(GEN)p1[i]);
1130: return gerepile(av,tetpil,rep);
1131: }
1132:
1133: for (i=1; i<lg(rep); i++)
1134: rep[i] = (long)unifpol(nf,(GEN)rep[i],0);
1135:
1136: nfcmbf.pol = polmod;
1137: nfcmbf.lt = lt;
1138: nfcmbf.h = h;
1139: nfcmbf.hinv = hinv;
1140: nfcmbf.den = den;
1141: nfcmbf.fact = rep;
1142: nfcmbf.res = cgetg(lg(rep)+1,t_VEC);
1143: nfcmbf.nfact = 0;
1144: nfcmbf.nfactmod = lg(rep)-1;
1145: nf_combine_factors(nf,1,NULL,d-3,1);
1146:
1147: if (DEBUGLEVEL >= 4) msgtimer("computation of the factors");
1148:
1149: i = nfcmbf.nfact;
1150: if (lgef(nfcmbf.pol)>3)
1151: {
1152: nfcmbf.res[++i] = (long) nf_pol_divres(nf,nfcmbf.pol,nfcmbf.lt,NULL);
1153: nfcmbf.nfact = i;
1154: }
1155:
1156: tetpil=avma; rep=cgetg(i+1,t_VEC);
1157: for ( ; i>=1; i--)
1158: rep[i]=(long)unifpol(nf,(GEN)nfcmbf.res[i],1);
1159: return gerepile(av,tetpil,rep);
1160: }
1161:
1162: static int
1163: nf_combine_factors(GEN nf,long fxn,GEN psf,long dlim,long hint)
1164: {
1165: int val = 0; /* assume failure */
1166: GEN newf, newpsf = NULL;
1167: long av,newd,ltop,i;
1168:
1169: /* Assertion: fxn <= nfcmbf.nfactmod && dlim > 0 */
1170:
1171: /* first, try deeper factors without considering the current one */
1172: if (fxn != nfcmbf.nfactmod)
1173: {
1174: val = nf_combine_factors(nf,fxn+1,psf,dlim,hint);
1175: if (val && psf) return 1;
1176: }
1177:
1178: /* second, try including the current modular factor in the product */
1179: newf = (GEN)nfcmbf.fact[fxn];
1180: if (!newf) return val; /* modular factor already used */
1181: newd = degpol(newf);
1182: if (newd > dlim) return val; /* degree of new factor is too large */
1183:
1184: av = avma;
1185: if (newd % hint == 0)
1186: {
1187: GEN p, quot,rem;
1188:
1189: newpsf = nf_pol_mul(nf, psf? psf: nfcmbf.lt, newf);
1190: newpsf = nf_pol_lift(nfcmbf.h,nfcmbf.hinv,nfcmbf.den,newpsf);
1191: /* try out the new combination */
1192: ltop=avma;
1193: quot=nf_pol_divres(nf,nfcmbf.pol,newpsf,&rem);
1194: if (gcmp0(rem)) /* found a factor */
1195: {
1196: p = nf_pol_mul(nf,element_inv(nf,leading_term(newpsf)),newpsf);
1197: nfcmbf.res[++nfcmbf.nfact] = (long) p; /* store factor */
1198: nfcmbf.fact[fxn]=0; /* remove used modular factor */
1199:
1200: /* fix up target */
1201: p=gun; quot=unifpol(nf,quot,0);
1202: for (i=2; i<lgef(quot); i++)
1203: p = mpppcm(p, denom((GEN)quot[i]));
1204:
1205: nfcmbf.pol = nf_pol_mul(nf,p,quot);
1206: nfcmbf.lt = leading_term(nfcmbf.pol);
1207: return 1;
1208: }
1209: avma=ltop;
1210: }
1211:
1212: /* If room in degree limit + more modular factors to try, add more factors to
1213: * newpsf */
1214: if (newd < dlim && fxn < nfcmbf.nfactmod
1215: && nf_combine_factors(nf,fxn+1,newpsf,dlim-newd,hint))
1216: {
1217: nfcmbf.fact[fxn]=0; /* remove used modular factor */
1218: return 1;
1219: }
1220: avma = av; return val;
1221: }
1222:
1223: /* return the characteristic polynomial of alpha over nf, where alpha
1224: is an element of the algebra nf[X]/(T) given as a polynomial in X */
1225: GEN
1226: rnfcharpoly(GEN nf,GEN T,GEN alpha,int v)
1227: {
1228: long av = avma, vnf, vT, lT;
1229: GEN p1;
1230:
1231: nf=checknf(nf); vnf = varn(nf[1]);
1232: if (v<0) v = 0;
1233: T = fix_relative_pol(nf,T,1);
1234: if (typ(alpha) == t_POLMOD) alpha = lift_to_pol(alpha);
1235: lT = lgef(T);
1236: if (typ(alpha) != t_POL || varn(alpha) == vnf)
1237: return gerepileupto(av, gpowgs(gsub(polx[v], alpha), lT - 3));
1238: vT = varn(T);
1239: if (varn(alpha) != vT || v >= vnf)
1240: err(talker,"incorrect variables in rnfcharpoly");
1241: if (lgef(alpha) >= lT) alpha = gmod(alpha,T);
1242: if (lT <= 4)
1243: return gerepileupto(av, gsub(polx[v], alpha));
1244: p1 = caract2(unifpol(nf,T,1), unifpol(nf,alpha,1), v);
1245: return gerepileupto(av, unifpol(nf,p1,1));
1246: }
1247:
1248: #if 0
1249: /* return the minimal polynomial of alpha over nf, where alpha is
1250: an element of the algebra nf[X]/(T) given as a polynomial in X */
1251: GEN
1252: rnfminpoly(GEN nf,GEN T,GEN alpha,int n)
1253: {
1254: long av=avma,tetpil;
1255: GEN p1,p2;
1256:
1257: nf=checknf(nf); p1=rnfcharpoly(nf,T,alpha,n);
1258: tetpil=avma; p2=nf_pol_subres(nf,p1,deriv(p1,varn(T)));
1259: if (lgef(p2)==3) { avma=tetpil; return p1; }
1260:
1261: p1 = nf_pol_divres(nf,p1,p2,NULL);
1262: p2 = element_inv(nf,leading_term(p1));
1263: tetpil=avma; return gerepile(av,tetpil,unifpol(nf,nf_pol_mul(nf,p2,p1),1));
1264: }
1265: #endif
1266:
1267: /* relative Dedekind criterion over nf, applied to the order defined by a
1268: * root of irreducible polynomial T, modulo the prime ideal pr. Returns
1269: * [flag,basis,val], where basis is a pseudo-basis of the enlarged order,
1270: * flag is 1 iff this order is pr-maximal, and val is the valuation in pr of
1271: * the order discriminant
1272: */
1273: GEN
1274: rnfdedekind(GEN nf,GEN T,GEN pr)
1275: {
1276: long av=avma,vt,r,d,da,n,m,i,j;
1277: GEN p1,p2,p,tau,g,vecun,veczero,matid;
1278: GEN prhall,res,h,k,base,Ca;
1279:
1280: nf=checknf(nf); Ca=unifpol(nf,T,0);
1281: res=cgetg(4,t_VEC);
1282: if (typ(pr)==t_VEC && lg(pr)==3)
1283: { prhall = (GEN)pr[2]; pr = (GEN)pr[1]; }
1284: else
1285: prhall=nfmodprinit(nf,pr);
1286: p=(GEN)pr[1]; tau=gdiv((GEN)pr[5],p);
1287: n=degpol(nf[1]); m=degpol(T);
1288:
1289: vecun = gscalcol_i(gun,n);
1290: veczero = zerocol(n);
1291:
1292: p1=(GEN)nffactormod(nf,Ca,pr)[1];
1293: r=lg(p1); if (r < 2) err(constpoler,"rnfdedekind");
1294: g=lift((GEN)p1[1]);
1295: for (i=2; i<r; i++)
1296: g = nf_pol_mul(nf,g, lift((GEN)p1[i]));
1297: h=nfmod_pol_divres(nf,prhall,Ca,g,NULL);
1298: k=nf_pol_mul(nf,tau,gsub(Ca, nf_pol_mul(nf,lift(g),lift(h))));
1299: p2=nfmod_pol_gcd(nf,prhall,g,h);
1300: k= nfmod_pol_gcd(nf,prhall,p2,k);
1301:
1302: d=degpol(k); /* <= m */
1303: vt = idealval(nf,discsr(T),pr) - 2*d;
1304: res[3]=lstoi(vt);
1305: if (!d || vt<=1) res[1]=un; else res[1]=zero;
1306:
1307: base=cgetg(3,t_VEC);
1308: p1=cgetg(m+d+1,t_MAT); base[1]=(long)p1;
1309: p2=cgetg(m+d+1,t_VEC); base[2]=(long)p2;
1310: /* if d > 0, base[2] temporarily multiplied by p, for the final nfhermitemod
1311: * [ which requires integral ideals ] */
1312: matid = gscalmat(d? p: gun, n);
1313: for (j=1; j<=m; j++)
1314: {
1315: p2[j]=(long)matid;
1316: p1[j]=lgetg(m+1,t_COL);
1317: for (i=1; i<=m; i++)
1318: coeff(p1,i,j) = (i==j)?(long)vecun:(long)veczero;
1319: }
1320: if (d)
1321: {
1322: GEN pal = lift(nfmod_pol_divres(nf,prhall,Ca,k,NULL));
1323: GEN prinvp = pidealprimeinv(nf,pr); /* again multiplied by p */
1324: GEN nfx = unifpol(nf,polx[varn(T)],0);
1325:
1326: for ( ; j<=m+d; j++)
1327: {
1328: p1[j]=lgetg(m+1,t_COL);
1329: da=degpol(pal);
1330: for (i=1; i<=da+1; i++) coeff(p1,i,j)=pal[i+1];
1331: for ( ; i<=m; i++) coeff(p1,i,j)=(long)veczero;
1332: p2[j]=(long)prinvp;
1333: nf_pol_divres(nf,nf_pol_mul(nf,pal,nfx),T,&pal);
1334: }
1335: /* the modulus is integral */
1336: base = nfhermitemod(nf,base, gmul(gpowgs(p, m-d),
1337: idealpows(nf, prinvp, d)));
1338: base[2] = ldiv((GEN)base[2], p); /* cancel the factor p */
1339: }
1340: res[2]=(long)base; return gerepilecopy(av, res);
1341: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to nffactor.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 / modules