Annotation of OpenXM_contrib/pari-2.2/src/basemath/base2.c, Revision 1.1.1.1
1.1 noro 1: /* $Id: base2.c,v 1.87 2001/10/01 12:11:28 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: /* MAXIMAL ORDERS */
19: /* */
20: /*******************************************************************/
21: #include "pari.h"
22:
23: extern GEN caractducos(GEN p, GEN x, int v);
24: extern GEN element_muli(GEN nf, GEN x, GEN y);
25: extern GEN element_mulid(GEN nf, GEN x, long i);
26: extern GEN eleval(GEN f,GEN h,GEN a);
27: extern GEN ideal_better_basis(GEN nf, GEN x, GEN M);
28: extern long int_elt_val(GEN nf, GEN x, GEN p, GEN bp, GEN *t, long v);
29: extern GEN mat_to_vecpol(GEN x, long v);
30: extern GEN nfidealdet1(GEN nf, GEN a, GEN b);
31: extern GEN nfsuppl(GEN nf, GEN x, long n, GEN prhall);
32: extern GEN pol_to_monic(GEN pol, GEN *lead);
33: extern GEN pol_to_vec(GEN x, long N);
34: extern GEN quicktrace(GEN x, GEN sym);
35: extern GEN respm(GEN f1,GEN f2,GEN pm);
36:
37: static void
38: allbase_check_args(GEN f, long code, GEN *y, GEN *ptw1, GEN *ptw2)
39: {
40: GEN w;
41: if (typ(f)!=t_POL) err(notpoler,"allbase");
42: if (degpol(f) <= 0) err(constpoler,"allbase");
43: if (DEBUGLEVEL) timer2();
44: switch(code)
45: {
46: case 0: case 1:
47: *y = ZX_disc(f);
48: if (!signe(*y)) err(talker,"reducible polynomial in allbase");
49: w = auxdecomp(absi(*y),1-code);
50: break;
51: default:
52: w = (GEN)code;
53: *y = factorback(w, NULL);
54: }
55: if (DEBUGLEVEL) msgtimer("disc. factorisation");
56: *ptw1 = (GEN)w[1];
57: *ptw2 = (GEN)w[2];
58: }
59:
60: /*******************************************************************/
61: /* */
62: /* ROUND 2 */
63: /* */
64: /*******************************************************************/
65: /* Normalized quotient and remainder ( -1/2 |y| < r = x-q*y <= 1/2 |y| ) */
66: static GEN
67: rquot(GEN x, GEN y)
68: {
69: long av=avma,av1;
70: GEN u,v,w,p;
71:
72: u=absi(y); v=shifti(x,1); w=shifti(y,1);
73: if (cmpii(u,v)>0) p=subii(v,u);
74: else p=addsi(-1,addii(u,v));
75: av1=avma; return gerepile(av,av1,divii(p,w));
76: }
77:
78: /* space needed lx + 2*ly */
79: static GEN
80: rrmdr(GEN x, GEN y)
81: {
82: long av=avma,tetpil,k;
83: GEN r,ys2;
84:
85: if (!signe(x)) return gzero;
86: r = resii(x,y); tetpil = avma;
87: ys2 = shifti(y,-1);
88: k = absi_cmp(r, ys2);
89: if (k>0 || (k==0 && signe(r)>0))
90: {
91: avma = tetpil;
92: if (signe(y) == signe(r)) r = subii(r,y); else r = addii(r,y);
93: return gerepile(av,tetpil,r);
94: }
95: avma = tetpil; return r;
96: }
97:
98: /* companion matrix of unitary polynomial x */
99: static GEN
100: companion(GEN x) /* cf assmat */
101: {
102: long i,j,l;
103: GEN y;
104:
105: l=degpol(x)+1; y=cgetg(l,t_MAT);
106: for (j=1; j<l; j++)
107: {
108: y[j] = lgetg(l,t_COL);
109: for (i=1; i<l-1; i++)
110: coeff(y,i,j)=(i+1==j)? un: zero;
111: coeff(y,i,j) = lneg((GEN)x[j+1]);
112: }
113: return y;
114: }
115:
116: /* assume x, y are square integer matrices of same dim. Multiply them */
117: static GEN
118: mulmati(GEN x, GEN y)
119: {
120: long n = lg(x),i,j,k,av;
121: GEN z = cgetg(n,t_MAT),p1,p2;
122:
123: for (j=1; j<n; j++)
124: {
125: z[j] = lgetg(n,t_COL);
126: for (i=1; i<n; i++)
127: {
128: p1=gzero; av=avma;
129: for (k=1; k<n; k++)
130: {
131: p2=mulii(gcoeff(x,i,k),gcoeff(y,k,j));
132: if (p2 != gzero) p1=addii(p1,p2);
133: }
134: coeff(z,i,j)=lpileupto(av,p1);
135: }
136: }
137: return z;
138: }
139:
140: static GEN
141: powmati(GEN x, long m)
142: {
143: long av=avma,j;
144: GEN y = x;
145:
146: j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
147: for (; j; m<<=1,j--)
148: {
149: y=mulmati(y,y);
150: if (m<0) y=mulmati(y,x);
151: }
152: return gerepileupto(av,y);
153: }
154:
155: static GEN
156: rtran(GEN v, GEN w, GEN q)
157: {
158: long av,tetpil;
159: GEN p1;
160:
161: if (signe(q))
162: {
163: av=avma; p1=gneg(gmul(q,w)); tetpil=avma;
164: return gerepile(av,tetpil,gadd(v,p1));
165: }
166: return v;
167: }
168:
169: /* return (v - qw) mod m (only compute entries k0,..,n)
170: * v and w are expected to have entries smaller than m */
171: static GEN
172: mtran(GEN v, GEN w, GEN q, GEN m, long k0)
173: {
174: long k,l;
175: GEN p1;
176:
177: if (signe(q))
178: {
179: l = lgefint(m) << 2;
180: for (k=lg(v)-1; k>= k0; k--)
181: {
182: long av = avma; (void)new_chunk(l);
183: p1 = subii((GEN)v[k], mulii(q,(GEN)w[k]));
184: avma = av; v[k]=(long)rrmdr(p1, m);
185: }
186: }
187: return v;
188: }
189:
190: /* entries of v and w are C small integers */
191: static GEN
192: mtran_long(GEN v, GEN w, long q, long m, long k0)
193: {
194: long k, p1;
195:
196: if (q)
197: {
198: for (k=lg(v)-1; k>= k0; k--)
199: {
200: p1 = v[k] - q * w[k];
201: v[k] = p1 % m;
202: }
203: }
204: return v;
205: }
206:
207: /* coeffs of a are C-long integers */
208: static void
209: rowred_long(GEN a, long rmod)
210: {
211: long q,j,k,pro, c = lg(a), r = lg(a[1]);
212:
213: for (j=1; j<r; j++)
214: {
215: for (k=j+1; k<c; k++)
216: while (coeff(a,j,k))
217: {
218: q = coeff(a,j,j) / coeff(a,j,k);
219: pro=(long)mtran_long((GEN)a[j],(GEN)a[k],q,rmod, j);
220: a[j]=a[k]; a[k]=pro;
221: }
222: if (coeff(a,j,j) < 0)
223: for (k=j; k<r; k++) coeff(a,k,j)=-coeff(a,k,j);
224: for (k=1; k<j; k++)
225: {
226: q = coeff(a,j,k) / coeff(a,j,j);
227: a[k]=(long)mtran_long((GEN)a[k],(GEN)a[j],q,rmod, k);
228: }
229: }
230: /* don't update the 0s in the last columns */
231: for (j=1; j<r; j++)
232: for (k=1; k<r; k++) coeff(a,j,k) = lstoi(coeff(a,j,k));
233: }
234:
235: static void
236: rowred(GEN a, GEN rmod)
237: {
238: long j,k,pro, c = lg(a), r = lg(a[1]);
239: long av=avma, lim=stack_lim(av,1);
240: GEN q;
241:
242: for (j=1; j<r; j++)
243: {
244: for (k=j+1; k<c; k++)
245: while (signe(gcoeff(a,j,k)))
246: {
247: q=rquot(gcoeff(a,j,j),gcoeff(a,j,k));
248: pro=(long)mtran((GEN)a[j],(GEN)a[k],q,rmod, j);
249: a[j]=a[k]; a[k]=pro;
250: }
251: if (signe(gcoeff(a,j,j)) < 0)
252: for (k=j; k<r; k++) coeff(a,k,j)=lnegi(gcoeff(a,k,j));
253: for (k=1; k<j; k++)
254: {
255: q=rquot(gcoeff(a,j,k),gcoeff(a,j,j));
256: a[k]=(long)mtran((GEN)a[k],(GEN)a[j],q,rmod, k);
257: }
258: if (low_stack(lim, stack_lim(av,1)))
259: {
260: long j1,k1;
261: GEN p1 = a;
262: if(DEBUGMEM>1) err(warnmem,"rowred j=%ld", j);
263: p1 = gerepilecopy(av,a);
264: for (j1=1; j1<r; j1++)
265: for (k1=1; k1<c; k1++) coeff(a,j1,k1) = coeff(p1,j1,k1);
266: }
267: }
268: }
269:
270: /* Calcule d/x ou d est entier et x matrice triangulaire inferieure
271: * entiere dont les coeff diagonaux divisent d (resultat entier).
272: */
273: static GEN
274: matinv(GEN x, GEN d, long n)
275: {
276: long i,j,k,av,av1;
277: GEN y,h;
278:
279: y=idmat(n);
280: for (i=1; i<=n; i++)
281: coeff(y,i,i)=ldivii(d,gcoeff(x,i,i));
282: av=avma;
283: for (i=2; i<=n; i++)
284: for (j=i-1; j; j--)
285: {
286: for (h=gzero,k=j+1; k<=i; k++)
287: {
288: GEN p1 = mulii(gcoeff(y,i,k),gcoeff(x,k,j));
289: if (p1 != gzero) h=addii(h,p1);
290: }
291: setsigne(h,-signe(h)); av1=avma;
292: coeff(y,i,j) = lpile(av,av1,divii(h,gcoeff(x,j,j)));
293: av = avma;
294: }
295: return y;
296: }
297:
298: static GEN
299: ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
300: {
301: long sp,hard_case_exponent,i,n=lg(cf)-1,av=avma, av2,limit;
302: GEN T,T2,Tn,m,v,delta, *w;
303: const GEN pp = sqri(p);
304: const long pps = (2*expi(pp)+2<BITS_IN_LONG)? pp[2]: 0;
305:
306: if (cmpis(p,n) > 0)
307: {
308: hard_case_exponent = 0;
309: sp = 0; /* gcc -Wall */
310: }
311: else
312: {
313: long k;
314: k = sp = itos(p);
315: i=1; while (k < n) { k *= sp; i++; }
316: hard_case_exponent = i;
317: }
318: T=cgetg(n+1,t_MAT); for (i=1; i<=n; i++) T[i]=lgetg(n+1,t_COL);
319: T2=cgetg(2*n+1,t_MAT); for (i=1; i<=2*n; i++) T2[i]=lgetg(n+1,t_COL);
320: Tn=cgetg(n*n+1,t_MAT); for (i=1; i<=n*n; i++) Tn[i]=lgetg(n+1,t_COL);
321: v = new_chunk(n+1);
322: w = (GEN*)new_chunk(n+1);
323:
324: av2 = avma; limit = stack_lim(av2,1);
325: delta=gun; m=idmat(n);
326:
327: for(;;)
328: {
329: long j,k,h, av0 = avma;
330: GEN t,b,jp,hh,index,p1, dd = sqri(delta), ppdd = mulii(dd,pp);
331:
332: if (DEBUGLEVEL > 3)
333: fprintferr("ROUND2: epsilon = %ld\tavma = %ld\n",epsilon,avma);
334:
335: b=matinv(m,delta,n);
336: for (i=1; i<=n; i++)
337: {
338: for (j=1; j<=n; j++)
339: for (k=1; k<=n; k++)
340: {
341: p1 = j==k? gcoeff(m,i,1): gzero;
342: for (h=2; h<=n; h++)
343: {
344: GEN p2 = mulii(gcoeff(m,i,h),gcoeff(cf[h],j,k));
345: if (p2!=gzero) p1 = addii(p1,p2);
346: }
347: coeff(T,j,k) = (long)rrmdr(p1, ppdd);
348: }
349: p1 = mulmati(m, mulmati(T,b));
350: for (j=1; j<=n; j++)
351: for (k=1; k<=n; k++)
352: coeff(p1,j,k)=(long)rrmdr(divii(gcoeff(p1,j,k),dd),pp);
353: w[i] = p1;
354: }
355:
356: if (hard_case_exponent)
357: {
358: for (j=1; j<=n; j++)
359: {
360: for (i=1; i<=n; i++) coeff(T,i,j) = coeff(w[j],1,i);
361: /* ici la boucle en k calcule la puissance p mod p de w[j] */
362: for (k=1; k<sp; k++)
363: {
364: for (i=1; i<=n; i++)
365: {
366: p1 = gzero;
367: for (h=1; h<=n; h++)
368: {
369: GEN p2=mulii(gcoeff(T,h,j),gcoeff(w[j],h,i));
370: if (p2!=gzero) p1 = addii(p1,p2);
371: }
372: v[i] = lmodii(p1, p);
373: }
374: for (i=1; i<=n; i++) coeff(T,i,j)=v[i];
375: }
376: }
377: t = powmati(T, hard_case_exponent);
378: }
379: else
380: {
381: for (i=1; i<=n; i++)
382: for (j=1; j<=n; j++)
383: {
384: long av1 = avma;
385: p1 = gzero;
386: for (k=1; k<=n; k++)
387: for (h=1; h<=n; h++)
388: {
389: const GEN r=modii(gcoeff(w[i],k,h),p);
390: const GEN s=modii(gcoeff(w[j],h,k),p);
391: const GEN p2 = mulii(r,s);
392: if (p2!=gzero) p1 = addii(p1,p2);
393: }
394: coeff(T,i,j) = lpileupto(av1,p1);
395: }
396: t = T;
397: }
398:
399: if (pps)
400: {
401: long ps = p[2];
402: for (i=1; i<=n; i++)
403: for (j=1; j<=n; j++)
404: {
405: coeff(T2,j,i)=(i==j)? ps: 0;
406: coeff(T2,j,n+i)=smodis(gcoeff(t,i,j),ps);
407: }
408: rowred_long(T2,pps);
409: }
410: else
411: {
412: for (i=1; i<=n; i++)
413: for (j=1; j<=n; j++)
414: {
415: coeff(T2,j,i)=(i==j)? (long)p: zero;
416: coeff(T2,j,n+i)=lmodii(gcoeff(t,i,j),p);
417: }
418: rowred(T2,pp);
419: }
420: jp=matinv(T2,p,n);
421: if (pps)
422: {
423: for (k=1; k<=n; k++)
424: {
425: long av1=avma;
426: t = mulmati(mulmati(jp,w[k]), T2);
427: for (h=i=1; i<=n; i++)
428: for (j=1; j<=n; j++)
429: { coeff(Tn,k,h) = itos(divii(gcoeff(t,i,j), p)) % pps; h++; }
430: avma=av1;
431: }
432: avma = av0;
433: rowred_long(Tn,pps);
434: }
435: else
436: {
437: for (k=1; k<=n; k++)
438: {
439: t = mulmati(mulmati(jp,w[k]), T2);
440: for (h=i=1; i<=n; i++)
441: for (j=1; j<=n; j++)
442: { coeff(Tn,k,h) = ldivii(gcoeff(t,i,j), p); h++; }
443: }
444: rowred(Tn,pp);
445: }
446: for (index=gun,i=1; i<=n; i++)
447: index = mulii(index,gcoeff(Tn,i,i));
448: if (gcmp1(index)) break;
449:
450: m = mulmati(matinv(Tn,index,n), m);
451: hh = delta = mulii(index,delta);
452: for (i=1; i<=n; i++)
453: for (j=1; j<=n; j++)
454: hh = mppgcd(gcoeff(m,i,j),hh);
455: if (!is_pm1(hh))
456: {
457: m = gdiv(m,hh);
458: delta = divii(delta,hh);
459: }
460: epsilon -= 2 * ggval(index,p);
461: if (epsilon < 2) break;
462: if (low_stack(limit,stack_lim(av2,1)))
463: {
464: GEN *gptr[3]; gptr[0]=&m; gptr[1]=δ
465: if(DEBUGMEM>1) err(warnmem,"ordmax");
466: gerepilemany(av2, gptr,2);
467: }
468: }
469: {
470: GEN *gptr[2]; gptr[0]=&m; gptr[1]=δ
471: gerepilemany(av,gptr,2);
472: }
473: *ptdelta=delta; return m;
474: }
475:
476: /* Input:
477: * x normalized integral polynomial of degree n, defining K=Q(theta).
478: *
479: * code 0, 1 or (long)p if we want base, smallbase ou factoredbase (resp.).
480: * y is GEN *, which will receive the discriminant of K.
481: *
482: * Output
483: * 1) A t_COL whose n components are rationnal polynomials (with degree
484: * 0,1...n-1) : integral basis for K (putting x=theta).
485: * Rem: common denominator is in da.
486: *
487: * 2) discriminant of K (in *y).
488: */
489: GEN
490: allbase(GEN f, long code, GEN *y)
491: {
492: GEN w1,w2,a,pro,at,bt,b,da,db,q, *cf,*gptr[2];
493: long av=avma,tetpil,n,h,j,i,k,r,s,t,v,mf;
494:
495: allbase_check_args(f,code,y, &w1,&w2);
496: v = varn(f); n = degpol(f); h = lg(w1)-1;
497: cf = (GEN*)cgetg(n+1,t_VEC);
498: cf[2]=companion(f);
499: for (i=3; i<=n; i++) cf[i]=mulmati(cf[2],cf[i-1]);
500:
501: a=idmat(n); da=gun;
502: for (i=1; i<=h; i++)
503: {
504: long av1 = avma;
505: mf=itos((GEN)w2[i]); if (mf==1) continue;
506: if (DEBUGLEVEL) fprintferr("Treating p^k = %Z^%ld\n",w1[i],mf);
507:
508: b=ordmax(cf,(GEN)w1[i],mf,&db);
509: a=gmul(db,a); b=gmul(da,b);
510: da=mulii(db,da);
511: at=gtrans(a); bt=gtrans(b);
512: for (r=n; r; r--)
513: for (s=r; s; s--)
514: while (signe(gcoeff(bt,s,r)))
515: {
516: q=rquot(gcoeff(at,s,s),gcoeff(bt,s,r));
517: pro=rtran((GEN)at[s],(GEN)bt[r],q);
518: for (t=s-1; t; t--)
519: {
520: q=rquot(gcoeff(at,t,s),gcoeff(at,t,t));
521: pro=rtran(pro,(GEN)at[t],q);
522: }
523: at[s]=bt[r]; bt[r]=(long)pro;
524: }
525: for (j=n; j; j--)
526: {
527: for (k=1; k<j; k++)
528: {
529: while (signe(gcoeff(at,j,k)))
530: {
531: q=rquot(gcoeff(at,j,j),gcoeff(at,j,k));
532: pro=rtran((GEN)at[j],(GEN)at[k],q);
533: at[j]=at[k]; at[k]=(long)pro;
534: }
535: }
536: if (signe(gcoeff(at,j,j))<0)
537: for (k=1; k<=j; k++) coeff(at,k,j)=lnegi(gcoeff(at,k,j));
538: for (k=j+1; k<=n; k++)
539: {
540: q=rquot(gcoeff(at,j,k),gcoeff(at,j,j));
541: at[k]=(long)rtran((GEN)at[k],(GEN)at[j],q);
542: }
543: }
544: for (j=2; j<=n; j++)
545: if (egalii(gcoeff(at,j,j), gcoeff(at,j-1,j-1)))
546: {
547: coeff(at,1,j)=zero;
548: for (k=2; k<=j; k++) coeff(at,k,j)=coeff(at,k-1,j-1);
549: }
550: tetpil=avma; a=gtrans(at);
551: {
552: GEN *gptr[2];
553: da = icopy(da); gptr[0]=&a; gptr[1]=&da;
554: gerepilemanysp(av1,tetpil,gptr,2);
555: }
556: }
557: for (j=1; j<=n; j++)
558: *y = divii(mulii(*y,sqri(gcoeff(a,j,j))), sqri(da));
559: tetpil=avma; *y=icopy(*y);
560: at=cgetg(n+1,t_VEC); v=varn(f);
561: for (k=1; k<=n; k++)
562: {
563: q=cgetg(k+2,t_POL); at[k]=(long)q;
564: q[1] = evalsigne(1) | evallgef(2+k) | evalvarn(v);
565: for (j=1; j<=k; j++) q[j+1]=ldiv(gcoeff(a,k,j),da);
566: }
567: gptr[0]=&at; gptr[1]=y;
568: gerepilemanysp(av,tetpil,gptr,2);
569: return at;
570: }
571:
572: GEN
573: base2(GEN x, GEN *y)
574: {
575: return allbase(x,0,y);
576: }
577:
578: GEN
579: discf2(GEN x)
580: {
581: GEN y;
582: long av=avma,tetpil;
583:
584: allbase(x,0,&y); tetpil=avma;
585: return gerepile(av,tetpil,icopy(y));
586: }
587:
588: /*******************************************************************/
589: /* */
590: /* ROUND 4 */
591: /* */
592: /*******************************************************************/
593:
594: GEN nilord(GEN p,GEN fx,long mf,GEN gx,long flag);
595: GEN Decomp(GEN p,GEN f,long mf,GEN theta,GEN chi,GEN nu,long r);
596: static GEN dbasis(GEN p, GEN f, long mf, GEN alpha, GEN U);
597: static GEN maxord(GEN p,GEN f,long mf);
598: static GEN nbasis(GEN ibas,GEN pd);
599: static GEN testb2(GEN p,GEN fa,long Fa,GEN theta,GEN pmf,long Ft,GEN ns);
600: static GEN testc2(GEN p,GEN fa,GEN pmr,GEN pmf,GEN alph2,
601: long Ea,GEN thet2,long Et,GEN ns);
602:
603: static int
604: fnz(GEN x,long j)
605: {
606: long i;
607: for (i=1; i<j; i++)
608: if (signe(x[i])) return 0;
609: return 1;
610: }
611:
612: /* retourne la base, dans y le discf et dans ptw la factorisation (peut
613: etre partielle) de discf */
614: GEN
615: allbase4(GEN f,long code, GEN *y, GEN *ptw)
616: {
617: GEN w,w1,w2,a,da,b,db,bas,q,p1,*gptr[3];
618: long v,n,mf,h,lfa,i,j,k,l,tetpil,av = avma;
619:
620: allbase_check_args(f,code,y, &w1,&w2);
621: v = varn(f); n = degpol(f); h = lg(w1)-1;
622: a = NULL; /* gcc -Wall */
623: da= NULL;
624: for (i=1; i<=h; i++)
625: {
626: mf=itos((GEN)w2[i]); if (mf == 1) continue;
627: if (DEBUGLEVEL) fprintferr("Treating p^k = %Z^%ld\n",w1[i],mf);
628:
629: b = maxord((GEN)w1[i],f,mf); db = gun;
630: for (j=1; j<=n; j++)
631: {
632: p1 = denom(gcoeff(b,j,j));
633: if (cmpii(p1,db) > 0) db = p1;
634: }
635: if (db != gun)
636: { /* db = denom(diag(b)), (da,db) = 1 */
637: b = gmul(b,db);
638: if (!da) { da=db; a=b; }
639: else
640: {
641: j=1; while (j<=n && fnz((GEN)a[j],j) && fnz((GEN)b[j],j)) j++;
642: b = gmul(da,b); a = gmul(db,a);
643: k=j-1; p1=cgetg(2*n-k+1,t_MAT);
644: for (j=1; j<=k; j++)
645: {
646: p1[j] = a[j];
647: coeff(p1,j,j) = lmppgcd(gcoeff(a,j,j),gcoeff(b,j,j));
648: }
649: for ( ; j<=n; j++) p1[j] = a[j];
650: for ( ; j<=2*n-k; j++) p1[j] = b[j+k-n];
651: da = mulii(da,db); a = hnfmodid(p1, da);
652: }
653: }
654: if (DEBUGLEVEL>5)
655: fprintferr("Result for prime %Z is:\n%Z\n",w1[i],b);
656: }
657: if (da)
658: {
659: for (j=1; j<=n; j++)
660: *y = mulii(divii(*y,sqri(da)),sqri(gcoeff(a,j,j)));
661: for (j=n-1; j; j--)
662: if (cmpis(gcoeff(a,j,j),2) > 0)
663: {
664: p1=shifti(gcoeff(a,j,j),-1);
665: for (k=j+1; k<=n; k++)
666: if (cmpii(gcoeff(a,j,k),p1) > 0)
667: for (l=1; l<=j; l++)
668: coeff(a,l,k)=lsubii(gcoeff(a,l,k),gcoeff(a,l,j));
669: }
670: }
671: lfa = 0;
672: if (ptw)
673: {
674: for (j=1; j<=h; j++)
675: {
676: k=ggval(*y,(GEN)w1[j]);
677: if (k) { lfa++; w1[lfa]=w1[j]; w2[lfa]=k; }
678: }
679: }
680: tetpil=avma; *y=icopy(*y);
681: bas=cgetg(n+1,t_VEC); v=varn(f);
682: for (k=1; k<=n; k++)
683: {
684: q=cgetg(k+2,t_POL); bas[k]=(long)q;
685: q[1] = evalsigne(1) | evallgef(k+2) | evalvarn(v);
686: if (da)
687: for (j=1; j<=k; j++) q[j+1]=ldiv(gcoeff(a,j,k),da);
688: else
689: {
690: for (j=2; j<=k; j++) q[j]=zero;
691: q[j]=un;
692: }
693: }
694: if (ptw)
695: {
696: *ptw=w=cgetg(3,t_MAT);
697: w[1]=lgetg(lfa+1,t_COL);
698: w[2]=lgetg(lfa+1,t_COL);
699: for (j=1; j<=lfa; j++)
700: {
701: coeff(w,j,1)=(long)icopy((GEN)w1[j]);
702: coeff(w,j,2)=lstoi(w2[j]);
703: }
704: gptr[2]=ptw;
705: }
706: gptr[0]=&bas; gptr[1]=y;
707: gerepilemanysp(av,tetpil,gptr, ptw?3:2);
708: return bas;
709: }
710:
711: extern GEN merge_factor_i(GEN f, GEN g);
712:
713: static GEN
714: update_fact(GEN x, GEN f)
715: {
716: GEN e,q,d = ZX_disc(x), g = cgetg(3, t_MAT), p = (GEN)f[1];
717: long iq,i,k,l;
718: if (typ(f)!=t_MAT || lg(f)!=3)
719: err(talker,"not a factorisation in nfbasis");
720: l = lg(p);
721: q = cgetg(l,t_COL); g[1]=(long)q;
722: e = cgetg(l,t_COL); g[2]=(long)e; iq = 1;
723: for (i=1; i<l; i++)
724: {
725: k = pvaluation(d, (GEN)p[i], &d);
726: if (k) { q[iq] = p[i]; e[iq] = lstoi(k); iq++; }
727: }
728: setlg(q,iq); setlg(e,iq);
729: return merge_factor_i(decomp(d), g);
730: }
731:
732: /* if y is non-NULL, it receives the discriminant
733: * return basis if (ret_basis != 0), discriminant otherwise
734: */
735: static GEN
736: nfbasis00(GEN x0, long flag, GEN p, long ret_basis, GEN *y)
737: {
738: GEN x, disc, basis, lead;
739: GEN *gptr[2];
740: long k, tetpil, av = avma, l = lgef(x0), smll;
741:
742: if (typ(x0)!=t_POL) err(typeer,"nfbasis00");
743: if (l<=3) err(zeropoler,"nfbasis00");
744: for (k=2; k<l; k++)
745: if (typ(x0[k])!=t_INT) err(talker,"polynomial not in Z[X] in nfbasis");
746:
747: x = pol_to_monic(x0,&lead);
748:
749: if (!p || gcmp0(p))
750: smll = (flag & 1); /* small basis */
751: else
752: {
753: if (lead) p = update_fact(x,p);
754: smll = (long) p; /* factored basis */
755: }
756:
757: if (flag & 2)
758: basis = allbase(x,smll,&disc); /* round 2 */
759: else
760: basis = allbase4(x,smll,&disc,NULL); /* round 4 */
761:
762: if (!ret_basis) return gerepilecopy(av,disc);
763:
764: tetpil=avma;
765: if (!lead) basis = gcopy(basis);
766: else
767: {
768: long v = varn(x);
769: GEN pol = gmul(polx[v],lead);
770:
771: tetpil = avma; basis = gsubst(basis,v,pol);
772: }
773: if (!y)
774: return gerepile(av,tetpil,basis);
775:
776: *y = gcopy(disc);
777: gptr[0]=&basis; gptr[1]=y;
778: gerepilemanysp(av,tetpil,gptr,2);
779: return basis;
780: }
781:
782: GEN
783: nfbasis(GEN x, GEN *y, long flag, GEN p)
784: {
785: return nfbasis00(x,flag,p,1,y);
786: }
787:
788: GEN
789: nfbasis0(GEN x, long flag, GEN p)
790: {
791: return nfbasis00(x,flag,p,1,NULL);
792: }
793:
794: GEN
795: nfdiscf0(GEN x, long flag, GEN p)
796: {
797: return nfbasis00(x,flag,p,0,&p);
798: }
799:
800: GEN
801: base(GEN x, GEN *y)
802: {
803: return allbase4(x,0,y,NULL);
804: }
805:
806: GEN
807: smallbase(GEN x, GEN *y)
808: {
809: return allbase4(x,1,y,NULL);
810: }
811:
812: GEN
813: factoredbase(GEN x, GEN p, GEN *y)
814: {
815: return allbase4(x,(long)p,y,NULL);
816: }
817:
818: GEN
819: discf(GEN x)
820: {
821: GEN y;
822: long av=avma,tetpil;
823:
824: allbase4(x,0,&y,NULL); tetpil=avma;
825: return gerepile(av,tetpil,icopy(y));
826: }
827:
828: GEN
829: smalldiscf(GEN x)
830: {
831: GEN y;
832: long av=avma,tetpil;
833:
834: allbase4(x,1,&y,NULL); tetpil=avma;
835: return gerepile(av,tetpil,icopy(y));
836: }
837:
838: GEN
839: factoreddiscf(GEN x, GEN p)
840: {
841: GEN y;
842: long av=avma,tetpil;
843:
844: allbase4(x,(long)p,&y,NULL); tetpil=avma;
845: return gerepile(av,tetpil,icopy(y));
846: }
847:
848: /* return U if Z[alpha] is not maximal or 2*dU < m-1; else return NULL */
849: static GEN
850: dedek(GEN f, long mf, GEN p,GEN g)
851: {
852: GEN k,h;
853: long dk;
854:
855: if (DEBUGLEVEL>=3)
856: {
857: fprintferr(" entering dedek ");
858: if (DEBUGLEVEL>5)
859: fprintferr("with parameters p=%Z,\n f=%Z",p,f);
860: fprintferr("\n");
861: }
862: h = FpX_div(f,g,p);
863: k = gdivexact(gadd(f, gneg_i(gmul(g,h))), p);
864: k = FpX_gcd(k, FpX_gcd(g,h, p), p);
865:
866: dk = degpol(k);
867: if (DEBUGLEVEL>=3) fprintferr(" gcd has degree %ld\n", dk);
868: if (2*dk >= mf-1) return FpX_div(f,k,p);
869: return dk? (GEN)NULL: f;
870: }
871:
872: /* p-maximal order of Af; mf = v_p(Disc(f)) */
873: static GEN
874: maxord(GEN p,GEN f,long mf)
875: {
876: long j,r, av = avma, flw = (cmpsi(degpol(f),p) < 0);
877: GEN w,g,h,res;
878:
879: if (flw)
880: {
881: h = NULL; r = 0; /* gcc -Wall */
882: g = FpX_div(f, FpX_gcd(f,derivpol(f), p), p);
883: }
884: else
885: {
886: w=(GEN)factmod(f,p)[1]; r=lg(w)-1;
887: g = h = lift_intern((GEN)w[r]); /* largest factor */
888: for (j=1; j<r; j++) g = FpX_red(gmul(g, lift_intern((GEN)w[j])), p);
889: }
890: res = dedek(f,mf,p,g);
891: if (res)
892: res = dbasis(p,f,mf,polx[varn(f)],res);
893: else
894: {
895: if (flw) { w=(GEN)factmod(f,p)[1]; r=lg(w)-1; h=lift_intern((GEN)w[r]); }
896: res = (r==1)? nilord(p,f,mf,h,0): Decomp(p,f,mf,polx[varn(f)],f,h,0);
897: }
898: return gerepileupto(av,res);
899: }
900:
901: /* do a centermod on integer or rational number */
902: static GEN
903: polmodiaux(GEN x, GEN y, GEN ys2)
904: {
905: if (typ(x)!=t_INT)
906: x = mulii((GEN)x[1], mpinvmod((GEN)x[2],y));
907: x = modii(x,y);
908: if (cmpii(x,ys2) > 0) x = subii(x,y);
909: return x;
910: }
911:
912: /* x polynomial with integer or rational coeff. Reduce them mod y IN PLACE */
913: GEN
914: polmodi(GEN x, GEN y)
915: {
916: long lx=lgef(x), i;
917: GEN ys2 = shifti(y,-1);
918: for (i=2; i<lx; i++) x[i]=(long)polmodiaux((GEN)x[i],y,ys2);
919: return normalizepol_i(x, lx);
920: }
921:
922: /* same but not in place */
923: GEN
924: polmodi_keep(GEN x, GEN y)
925: {
926: long lx=lgef(x), i;
927: GEN ys2 = shifti(y,-1);
928: GEN z = cgetg(lx,t_POL);
929: for (i=2; i<lx; i++) z[i]=(long)polmodiaux((GEN)x[i],y,ys2);
930: z[1]=x[1]; return normalizepol_i(z, lx);
931: }
932:
933: static GEN
934: dbasis(GEN p, GEN f, long mf, GEN alpha, GEN U)
935: {
936: long n=degpol(f),dU,c;
937: GEN b,ha,pd,pdp;
938:
939: if (n == 1) return gscalmat(gun, 1);
940: if (DEBUGLEVEL>=3)
941: {
942: fprintferr(" entering Dedekind Basis ");
943: if (DEBUGLEVEL>5)
944: {
945: fprintferr("with parameters p=%Z\n",p);
946: fprintferr(" f = %Z,\n alpha = %Z",f,alpha);
947: }
948: fprintferr("\n");
949: }
950: ha = pd = gpuigs(p,mf/2); pdp = mulii(pd,p);
951: dU = typ(U)==t_POL? degpol(U): 0;
952: b = cgetg(n,t_MAT); /* Z[a] + U/p Z[a] is maximal */
953: /* skip first column = gscalcol(pd,n) */
954: for (c=1; c<n; c++)
955: {
956: if (c == dU)
957: {
958: ha = gdiv(gmul(pd,eleval(f,U,alpha)),p);
959: ha = polmodi(ha,pdp);
960: }
961: else
962: {
963: GEN p2, mod;
964: ha = gmul(ha,alpha);
965: p2 = content(ha); /* to cancel denominator */
966: if (gcmp1(p2)) { p2 = NULL; mod = pdp; }
967: else
968: {
969: ha = gdiv(ha,p2);
970: if (typ(p2)==t_INT)
971: mod = divii(pdp, mppgcd(pdp,p2));
972: else
973: mod = mulii(pdp, (GEN)p2[2]); /* p2 = a / p^e */
974: }
975: ha = FpX_res(ha, f, mod);
976: if (p2) ha = gmul(ha,p2);
977: }
978: b[c] = (long)pol_to_vec(ha,n);
979: }
980: b = hnfmodid(b,pd);
981: if (DEBUGLEVEL>5) fprintferr(" new order: %Z\n",b);
982: return gdiv(b,pd);
983: }
984:
985: static GEN
986: get_partial_order_as_pols(GEN p, GEN f)
987: {
988: long i,j, n = degpol(f), vf = varn(f);
989: GEN b,ib,h,col;
990:
991: b = maxord(p,f, ggval(ZX_disc(f),p));
992: ib = cgetg(n+1,t_VEC);
993: for (i=1; i<=n; i++)
994: {
995: h=cgetg(i+2,t_POL); ib[i]=(long)h; col=(GEN)b[i];
996: h[1]=evalsigne(1)|evallgef(i+2)|evalvarn(vf);
997: for (j=1;j<=i;j++) h[j+1]=col[j];
998: }
999: return ib;
1000: }
1001:
1002: /* if flag != 0, factorization to precision r (maximal order otherwise) */
1003: GEN
1004: Decomp(GEN p,GEN f,long mf,GEN theta,GEN chi,GEN nu,long flag)
1005: {
1006: GEN res,pr,pk,ph,pdr,unmodp,b1,b2,b3,a1,e,f1,f2;
1007:
1008: if (DEBUGLEVEL>2)
1009: {
1010: fprintferr(" entering Decomp ");
1011: if (DEBUGLEVEL>5)
1012: {
1013: fprintferr("with parameters: p=%Z, expo=%ld\n",p,mf);
1014: if (flag) fprintferr("precision = %ld\n",flag);
1015: fprintferr(" f=%Z",f);
1016: }
1017: fprintferr("\n");
1018: }
1019: unmodp = gmodulsg(1,p);
1020: b1=lift_intern(gmul(chi,unmodp));
1021: a1=gun; b2=gun;
1022: b3=lift_intern(gmul(nu,unmodp));
1023: while (degpol(b3) > 0)
1024: {
1025: GEN p1;
1026: b1 = FpX_div(b1,b3, p);
1027: b2 = FpX_red(gmul(b2,b3), p);
1028: b3 = FpX_extgcd(b2,b1, p, &a1,&p1); /* p1 = junk */
1029: p1 = leading_term(b3);
1030: if (!gcmp1(p1))
1031: { /* FpX_extgcd does not return normalized gcd */
1032: p1 = mpinvmod(p1,p);
1033: b3 = gmul(b3,p1);
1034: a1 = gmul(a1,p1);
1035: }
1036: }
1037: pdr = respm(f,derivpol(f),gpuigs(p,mf+1));
1038: e = eleval(f,FpX_red(gmul(a1,b2), p),theta);
1039: e = gdiv(polmodi(gmul(pdr,e), mulii(pdr,p)),pdr);
1040:
1041: pr = flag? gpowgs(p,flag): mulii(p,sqri(pdr));
1042: pk=p; ph=mulii(pdr,pr);
1043: /* E(t) - e(t) belongs to p^k Op, which is contained in p^(k-df)*Zp[xi] */
1044: while (cmpii(pk,ph) < 0)
1045: {
1046: e = gmul(gsqr(e), gsubsg(3,gmul2n(e,1)));
1047: e = gres(e,f); pk = sqri(pk);
1048: e = gdiv(polmodi(gmul(pdr,e), mulii(pdr,pk)), pdr);
1049: }
1050: f1 = gcdpm(f,gmul(pdr,gsubsg(1,e)), ph);
1051: f1 = FpX_res(f1,f, pr);
1052: f2 = FpX_res(FpX_div(f,f1, pr), f, pr);
1053:
1054: if (DEBUGLEVEL>2)
1055: {
1056: fprintferr(" leaving Decomp");
1057: if (DEBUGLEVEL>5)
1058: fprintferr(" with parameters: f1 = %Z\nf2 = %Z\ne = %Z\n", f1,f2,e);
1059: fprintferr("\n");
1060: }
1061:
1062: if (flag)
1063: {
1064: b1=factorpadic4(f1,p,flag);
1065: b2=factorpadic4(f2,p,flag); res=cgetg(3,t_MAT);
1066: res[1]=lconcat((GEN)b1[1],(GEN)b2[1]);
1067: res[2]=lconcat((GEN)b1[2],(GEN)b2[2]); return res;
1068: }
1069: else
1070: {
1071: GEN ib1,ib2;
1072: long n1,n2,i;
1073: ib1 = get_partial_order_as_pols(p,f1); n1=lg(ib1)-1;
1074: ib2 = get_partial_order_as_pols(p,f2); n2=lg(ib2)-1;
1075: res=cgetg(n1+n2+1,t_VEC);
1076: for (i=1; i<=n1; i++)
1077: res[i]=(long)polmodi(gmod(gmul(gmul(pdr,(GEN)ib1[i]),e),f), pdr);
1078: e=gsubsg(1,e); ib2 -= n1;
1079: for ( ; i<=n1+n2; i++)
1080: res[i]=(long)polmodi(gmod(gmul(gmul(pdr,(GEN)ib2[i]),e),f), pdr);
1081: return nbasis(res,pdr);
1082: }
1083: }
1084:
1085: /* minimum extension valuation: res[0]/res[1] (both are longs) */
1086: static long *
1087: vstar(GEN p,GEN h)
1088: {
1089: static long res[2];
1090: long m,first,j,k,v,w;
1091:
1092: m=degpol(h); first=1; k=1; v=0;
1093: for (j=1; j<=m; j++)
1094: if (! gcmp0((GEN)h[m-j+2]))
1095: {
1096: w = ggval((GEN)h[m-j+2],p);
1097: if (first || w*k < v*j) { v=w; k=j; }
1098: first=0;
1099: }
1100: m = cgcd(v,k);
1101: res[0]=v/m; res[1]=k/m; return res;
1102: }
1103:
1104: /* reduce the element elt modulo rd, taking first care of the denominators */
1105: static GEN
1106: redelt(GEN elt, GEN rd, GEN pd)
1107: {
1108: GEN den, nelt, nrd, relt;
1109:
1110: den = ggcd(denom(content(elt)), pd);
1111: nelt = gmul(den, elt);
1112: nrd = gmul(den, rd);
1113:
1114: if (typ(elt) == t_POL)
1115: relt = polmodi(nelt, nrd);
1116: else
1117: relt = centermod(nelt, nrd);
1118:
1119: return gdiv(relt, den);
1120: }
1121:
1122: /* compute the Newton sums of g(x) mod pp from its coefficients */
1123: GEN
1124: polsymmodpp(GEN g, GEN pp)
1125: {
1126: long av1, av2, d = degpol(g), i, k;
1127: GEN s , y;
1128:
1129: y = cgetg(d + 1, t_COL);
1130: y[1] = lstoi(d);
1131: for (k = 1; k < d; k++)
1132: {
1133: av1 = avma;
1134: s = centermod(gmulsg(k, polcoeff0(g,d-k,-1)), pp);
1135: for (i = 1; i < k; i++)
1136: s = gadd(s, gmul((GEN)y[k-i+1], polcoeff0(g,d-i,-1)));
1137: av2 = avma;
1138: y[k + 1] = lpile(av1, av2, centermod(gneg(s), pp));
1139: }
1140:
1141: return y;
1142: }
1143:
1144: /* no GC */
1145: static GEN
1146: manage_cache(GEN chi, GEN pp, GEN ns)
1147: {
1148: long j, n = degpol(chi);
1149: GEN ns2, npp = (GEN)ns[n+1];
1150:
1151: if (gcmp(pp, npp) > 0)
1152: {
1153: if (DEBUGLEVEL > 4)
1154: fprintferr("newtonsums: result too large to fit in cache\n");
1155: return polsymmodpp(chi, pp);
1156: }
1157:
1158: if (!signe((GEN)ns[1]))
1159: {
1160: ns2 = polsymmodpp(chi, pp);
1161: for (j = 1; j <= n; j++)
1162: affii((GEN)ns2[j], (GEN)ns[j]);
1163: }
1164:
1165: return ns;
1166: }
1167:
1168: /* compute the Newton sums modulo pp of the characteristic
1169: polynomial of a(x) mod g(x) */
1170: static GEN
1171: newtonsums(GEN a, GEN chi, GEN pp, GEN ns)
1172: {
1173: GEN va, pa, s, ns2;
1174: long j, k, n = degpol(chi), av2, lim;
1175:
1176: ns2 = manage_cache(chi, pp, ns);
1177:
1178: av2 = avma;
1179: lim = stack_lim(av2, 1);
1180:
1181: pa = gun;
1182: va = zerovec(n);
1183:
1184: for (j = 1; j <= n; j++)
1185: {
1186: pa = gmul(pa, a);
1187: if (pp) pa = polmodi(pa, pp);
1188: pa = gmod(pa, chi);
1189: if (pp) pa = polmodi(pa, pp);
1190:
1191: s = gzero;
1192:
1193: for (k = 0; k <= n-1; k++)
1194: s = addii(s, mulii(polcoeff0(pa, k, -1), (GEN)ns2[k+1]));
1195:
1196: if (pp) va[j] = (long)centermod(s, pp);
1197:
1198: if (low_stack(lim, stack_lim(av2, 1)))
1199: {
1200: GEN *gptr[2];
1201: gptr[0]=&pa; gptr[1]=&va;
1202: if(DEBUGMEM>1) err(warnmem, "newtonsums");
1203: gerepilemany(av2, gptr, 2);
1204: }
1205: }
1206:
1207: return va;
1208: }
1209:
1210: /* compute the characteristic polynomial of a mod g
1211: to a precision of pp using Newton sums */
1212: static GEN
1213: newtoncharpoly(GEN a, GEN chi, GEN pp, GEN ns)
1214: {
1215: GEN v, c, s, t;
1216: long n = degpol(chi), j, k, vn = varn(chi), av = avma, av2, lim;
1217:
1218: v = newtonsums(a, chi, pp, ns);
1219: av2 = avma;
1220: lim = stack_lim(av2, 1);
1221: c = cgetg(n + 2, t_VEC);
1222: c[1] = un;
1223: if (n%2) c[1] = lneg((GEN)c[1]);
1224: for (k = 2; k <= n+1; k++) c[k] = zero;
1225:
1226: for (k = 2; k <= n+1; k++)
1227: {
1228: s = gzero;
1229: for (j = 1; j < k; j++)
1230: {
1231: t = gmul((GEN)v[j], (GEN)c[k-j]);
1232: if (!(j%2)) t = gneg(t);
1233: s = gadd(s, t);
1234: }
1235: c[k] = ldiv(s, stoi(k - 1));
1236:
1237: if (low_stack(lim, stack_lim(av2, 1)))
1238: {
1239: if(DEBUGMEM>1) err(warnmem, "newtoncharpoly");
1240: c = gerepilecopy(av2, c);
1241: }
1242: }
1243:
1244: k = (n%2)? 1: 2;
1245: for ( ; k <= n+1; k += 2)
1246: c[k] = lneg((GEN)c[k]);
1247:
1248: return gerepileupto(av, gtopoly(c, vn));
1249: }
1250:
1251: static GEN
1252: mycaract(GEN f, GEN beta, GEN p, GEN pp, GEN ns)
1253: {
1254: GEN p1, p2, chi, chi2, npp;
1255: long j, a, v = varn(f), n = degpol(f);
1256:
1257: if (gcmp0(beta)) return zeropol(v);
1258:
1259: p1 = content(beta);
1260: if (gcmp1(p1)) p1 = NULL;
1261: else beta = gdiv(beta, p1);
1262:
1263: if (!pp)
1264: chi = caractducos(f, beta, v);
1265: else
1266: {
1267: a = 0;
1268: for (j = 1; j <= n; j++) /* compute the extra precision needed */
1269: a += ggval(stoi(j), p);
1270: npp = mulii(pp, gpowgs(p, a));
1271: if (p1) npp = gmul(npp, gpowgs(denom(p1), n));
1272:
1273: chi = newtoncharpoly(beta, f, npp, ns);
1274: }
1275:
1276: if (p1)
1277: {
1278: chi2 = cgetg(n+3, t_POL);
1279: chi2[1] = chi[1];
1280: p2 = gun;
1281: for (j = n+2; j >= 2; j--)
1282: {
1283: chi2[j] = lmul((GEN)chi[j], p2);
1284: p2 = gmul(p2, p1);
1285: }
1286: }
1287: else
1288: chi2 = chi;
1289:
1290: if (!pp) return chi2;
1291:
1292: /* this can happen only if gamma is incorrect (not an integer) */
1293: if (divise(denom(content(chi2)), p)) return NULL;
1294:
1295: return redelt(chi2, pp, pp);
1296: }
1297:
1298: /* Factor characteristic polynomial of beta mod p */
1299: static GEN
1300: factcp(GEN p, GEN f, GEN beta, GEN pp, GEN ns)
1301: {
1302: long av,l;
1303: GEN chi,nu, b = cgetg(4,t_VEC);
1304:
1305: chi=mycaract(f,beta,p,pp,ns);
1306: av=avma; nu=(GEN)factmod(chi,p)[1]; l=lg(nu)-1;
1307: nu=lift_intern((GEN)nu[1]);
1308: b[1]=(long)chi;
1309: b[2]=lpilecopy(av,nu);
1310: b[3]=lstoi(l); return b;
1311: }
1312:
1313: /* return the prime element in Zp[phi] */
1314: static GEN
1315: getprime(GEN p, GEN chi, GEN phi, GEN chip, GEN nup, long *Lp, long *Ep)
1316: {
1317: long v = varn(chi), L, E, r, s;
1318: GEN chin, pip, pp, vn;
1319:
1320: if (gegal(nup, polx[v]))
1321: chin = chip;
1322: else
1323: chin = mycaract(chip, nup, p, NULL, NULL);
1324:
1325: vn = vstar(p, chin);
1326: L = vn[0];
1327: E = vn[1];
1328:
1329: cbezout(L, -E, &r, &s);
1330:
1331: if (r <= 0)
1332: {
1333: long q = (-r) / E;
1334: q++;
1335: r += q*E;
1336: s += q*L;
1337: }
1338:
1339: pip = eleval(chi, nup, phi);
1340: pip = lift_intern(gpuigs(gmodulcp(pip, chi), r));
1341: pp = gpuigs(p, s);
1342:
1343: *Lp = L;
1344: *Ep = E;
1345: return gdiv(pip, pp);
1346: }
1347:
1348: static GEN
1349: update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr, GEN pmf, long mf,
1350: GEN ns)
1351: {
1352: long l, v = varn(fx);
1353: GEN nalph = NULL, nchi, w, nnu, pdr, npmr, rep;
1354:
1355: affii(gzero, (GEN)ns[1]); /* kill cache */
1356:
1357: if (!chi)
1358: nchi = mycaract(fx, alph, p, pmf, ns);
1359: else
1360: {
1361: nchi = chi;
1362: nalph = alph;
1363: }
1364:
1365: pdr = respm(nchi, derivpol(nchi), pmr);
1366: for (;;)
1367: {
1368: if (signe(pdr)) break;
1369: if (!nalph) nalph = gadd(alph, gmul(p, polx[v]));
1370: /* nchi was too reduced at this point */
1371: nchi = mycaract(fx, nalph, NULL, NULL, ns);
1372: pdr = respm(nchi, derivpol(nchi), pmf);
1373: if (signe(pdr)) break;
1374: if (DEBUGLEVEL >= 6)
1375: fprintferr(" non separable polynomial in update_alpha!\n");
1376: /* at this point, we assume that chi is not square-free */
1377: nalph = gadd(nalph, gmul(p, polx[v]));
1378: w = factcp(p, fx, nalph, NULL, ns);
1379: nchi = (GEN)w[1];
1380: nnu = (GEN)w[2];
1381: l = itos((GEN)w[3]);
1382: if (l > 1) return Decomp(p, fx, mf, nalph, nchi, nnu, 0);
1383: pdr = respm(nchi, derivpol(nchi), pmr);
1384: }
1385:
1386: if (is_pm1(pdr))
1387: npmr = gun;
1388: else
1389: {
1390: npmr = mulii(sqri(pdr), p);
1391: nchi = polmodi(nchi, npmr);
1392: if (!nalph) nalph = redelt(alph, npmr, pmf);
1393: else nalph = redelt(nalph, npmr, pmf);
1394: }
1395:
1396: affii(gzero, (GEN)ns[1]); /* kill cache again (contains data for fx) */
1397:
1398: rep = cgetg(5, t_VEC);
1399: rep[1] = (long)nalph;
1400: rep[2] = (long)nchi;
1401: rep[3] = (long)npmr;
1402: rep[4] = lmulii(p, pdr);
1403:
1404: return rep;
1405: }
1406:
1407: extern GEN Fp_factor_irred(GEN P,GEN l, GEN Q);
1408:
1409: /* flag != 0 iff we're looking for the p-adic factorization,
1410: in which case it is the p-adic precision we want */
1411: GEN
1412: nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
1413: {
1414: long Fa, La, Ea, oE, Fg, eq, er, v = varn(fx), i, nv, Le, Ee, N, l, vn;
1415: GEN p1, alph, chi, nu, w, phi, pmf, pdr, pmr, kapp, pie, chib, ns;
1416: GEN gamm, chig, nug, delt, beta, eta, chie, nue, pia, vb, opa;
1417:
1418: if (DEBUGLEVEL >= 3)
1419: {
1420: if (flag)
1421: fprintferr(" entering Nilord2 (factorization)");
1422: else
1423: fprintferr(" entering Nilord2 (basis/discriminant)");
1424: if (DEBUGLEVEL >= 5)
1425: {
1426: fprintferr(" with parameters: p = %Z, expo = %ld\n", p, mf);
1427: fprintferr(" fx = %Z, gx = %Z", fx, gx);
1428: }
1429: fprintferr("\n");
1430: }
1431:
1432: pmf = gpowgs(p, mf + 1);
1433: pdr = respm(fx, derivpol(fx), pmf);
1434: pmr = mulii(sqri(pdr), p);
1435: pdr = mulii(p, pdr);
1436: chi = polmodi_keep(fx, pmr);
1437:
1438: alph = polx[v];
1439: nu = gx;
1440: N = degpol(fx);
1441: oE = 0;
1442: opa = NULL;
1443:
1444: /* used to cache the newton sums of chi */
1445: ns = cgetg(N + 2, t_COL);
1446: p1 = powgi(p, gceil(gdivsg(N, mulii(p, subis(p, 1))))); /* p^(N/(p(p-1))) */
1447: p1 = mulii(p1, mulii(pmf, gpowgs(pmr, N)));
1448: l = lg(p1); /* enough in general... */
1449: for (i = 1; i <= N + 1; i++) ns[i] = lgeti(l);
1450: ns[N+1] = (long)p1;
1451: affii(gzero, (GEN)ns[1]); /* zero means: need to be computed */
1452:
1453: for(;;)
1454: {
1455: /* kappa need to be recomputed */
1456: kapp = NULL;
1457: Fa = degpol(nu);
1458: /* the prime element in Zp[alpha] */
1459: pia = getprime(p, chi, polx[v], chi, nu, &La, &Ea);
1460: pia = redelt(pia, pmr, pmf);
1461:
1462: if (Ea < oE)
1463: {
1464: alph = gadd(alph, opa);
1465: w = update_alpha(p, fx, alph, NULL, pmr, pmf, mf, ns);
1466: alph = (GEN)w[1];
1467: chi = (GEN)w[2];
1468: pmr = (GEN)w[3];
1469: pdr = (GEN)w[4];
1470: kapp = NULL;
1471: pia = getprime(p, chi, polx[v], chi, nu, &La, &Ea);
1472: pia = redelt(pia, pmr, pmf);
1473: }
1474:
1475: oE = Ea; opa = eleval(fx, pia, alph);
1476:
1477: if (DEBUGLEVEL >= 5)
1478: fprintferr(" Fa = %ld and Ea = %ld \n", Fa, Ea);
1479:
1480: /* we change alpha such that nu = pia */
1481: if (La > 1)
1482: {
1483: alph = gadd(alph, eleval(fx, pia, alph));
1484:
1485: w = update_alpha(p, fx, alph, NULL, pmr, pmf, mf, ns);
1486: alph = (GEN)w[1];
1487: chi = (GEN)w[2];
1488: pmr = (GEN)w[3];
1489: pdr = (GEN)w[4];
1490: }
1491:
1492: /* if Ea*Fa == N then O = Zp[alpha] */
1493: if (Ea*Fa == N)
1494: {
1495: if (flag) return NULL;
1496:
1497: alph = redelt(alph, sqri(p), pmf);
1498: return dbasis(p, fx, mf, alph, p);
1499: }
1500:
1501: /* during the process beta tends to a factor of chi */
1502: beta = lift_intern(gpowgs(gmodulcp(nu, chi), Ea));
1503:
1504: for (;;)
1505: {
1506: if (DEBUGLEVEL >= 5)
1507: fprintferr(" beta = %Z\n", beta);
1508:
1509: /* if pmf divides norm(beta) then it's useless */
1510: p1 = gmod(gnorm(gmodulcp(beta, chi)), pmf);
1511: if (signe(p1))
1512: {
1513: chib = NULL;
1514: vn = ggval(p1, p);
1515: eq = (long)(vn / N);
1516: er = (long)(vn*Ea/N - eq*Ea);
1517: }
1518: else
1519: {
1520: chib = mycaract(chi, beta, NULL, NULL, ns);
1521: vb = vstar(p, chib);
1522: eq = (long)(vb[0] / vb[1]);
1523: er = (long)(vb[0]*Ea / vb[1] - eq*Ea);
1524: }
1525:
1526: /* eq and er are such that gamma = beta.p^-eq.nu^-er is a unit */
1527: if (eq) gamm = gdiv(beta, gpowgs(p, eq));
1528: else gamm = beta;
1529:
1530: if (er)
1531: {
1532: /* kappa = nu^-1 in Zp[alpha] */
1533: if (!kapp)
1534: {
1535: kapp = ginvmod(nu, chi);
1536: kapp = redelt(kapp, pmr, pmf);
1537: kapp = gmodulcp(kapp, chi);
1538: }
1539: gamm = lift(gmul(gamm, gpowgs(kapp, er)));
1540: gamm = redelt(gamm, p, pmf);
1541: }
1542:
1543: if (DEBUGLEVEL >= 6)
1544: fprintferr(" gamma = %Z\n", gamm);
1545:
1546: if (er || !chib)
1547: chig = mycaract(chi, gamm, p, pmf, ns);
1548: else
1549: {
1550: chig = poleval(chib, gmul(polx[v], gpowgs(p, eq)));
1551: chig = gdiv(chig, gpowgs(p, N*eq));
1552: chig = polmodi(chig, pmf);
1553: }
1554:
1555: if (!chig || !gcmp1(denom(content(chig))))
1556: {
1557: /* the valuation of beta was wrong... This also means
1558: that chi_gamma has more than one factor modulo p */
1559: if (!chig) chig = mycaract(chi, gamm, p, NULL, NULL);
1560:
1561: vb = vstar(p, chig);
1562: eq = (long)(-vb[0] / vb[1]);
1563: er = (long)(-vb[0]*Ea / vb[1] - eq*Ea);
1564: if (eq) gamm = gmul(gamm, gpowgs(p, eq));
1565: if (er)
1566: {
1567: gamm = gmul(gamm, gpowgs(nu, er));
1568: gamm = gmod(gamm, chi);
1569: gamm = redelt(gamm, p, pmr);
1570: }
1571: if (eq || er) chig = mycaract(chi, gamm, p, pmf, ns);
1572: }
1573:
1574: nug = (GEN)factmod(chig, p)[1];
1575: l = lg(nug) - 1;
1576: nug = lift((GEN)nug[l]);
1577:
1578: if (l > 1)
1579: {
1580: /* there are at least 2 factors mod. p => chi can be split */
1581: phi = eleval(fx, gamm, alph);
1582: phi = redelt(phi, p, pmf);
1583: if (flag) mf += 3;
1584: return Decomp(p, fx, mf, phi, chig, nug, flag);
1585: }
1586:
1587: Fg = degpol(nug);
1588: if (Fa%Fg)
1589: {
1590: if (DEBUGLEVEL >= 5)
1591: fprintferr(" Increasing Fa\n");
1592: /* we compute a new element such F = lcm(Fa, Fg) */
1593: w = testb2(p, chi, Fa, gamm, pmf, Fg, ns);
1594: if (gcmp1((GEN)w[1]))
1595: {
1596: /* there are at least 2 factors mod. p => chi can be split */
1597: phi = eleval(fx, (GEN)w[2], alph);
1598: phi = redelt(phi, p, pmf);
1599: if (flag) mf += 3;
1600: return Decomp(p, fx, mf, phi, (GEN)w[3], (GEN)w[4], flag);
1601: }
1602: break;
1603: }
1604:
1605: /* we look for a root delta of nug in Fp[alpha] such that
1606: vp(gamma - delta) > 0. This root can then be used to
1607: improved the approximation given by beta */
1608: nv = fetch_var();
1609: w = Fp_factor_irred(nug, p, gsubst(nu, varn(nu), polx[nv]));
1610: if (degpol(w[1]) != 1)
1611: err(talker,"bug in nilord (no root). Is p a prime ?");
1612:
1613: for (i = 1;; i++)
1614: {
1615: if (i >= lg(w))
1616: err(talker, "bug in nilord (no suitable root), is p = %Z a prime?",
1617: (long)p);
1618: delt = gneg_i(gsubst(gcoeff(w, 2, i), nv, polx[v]));
1619: eta = gsub(gamm, delt);
1620: if (typ(delt) == t_INT)
1621: {
1622: chie = poleval(chig, gadd(polx[v], delt));
1623: nue = (GEN)factmod(chie, p)[1];
1624: l = lg(nue) - 1;
1625: nue = lift((GEN)nue[l]);
1626: }
1627: else
1628: {
1629: p1 = factcp(p, chi, eta, pmf, ns);
1630: chie = (GEN)p1[1];
1631: nue = (GEN)p1[2];
1632: l = itos((GEN)p1[3]);
1633: }
1634: if (l > 1)
1635: {
1636: /* there are at least 2 factors mod. p => chi can be split */
1637: delete_var();
1638: phi = eleval(fx, eta, alph);
1639: phi = redelt(phi, p, pmf);
1640: if (flag) mf += 3;
1641: return Decomp(p, fx, mf, phi, chie, nue, flag);
1642: }
1643:
1644: /* if vp(eta) = vp(gamma - delta) > 0 */
1645: if (gegal(nue, polx[v])) break;
1646: }
1647: delete_var();
1648:
1649: if (!signe(modii((GEN)chie[2], pmr)))
1650: chie = mycaract(chi, eta, p, pmf, ns);
1651:
1652: pie = getprime(p, chi, eta, chie, nue, &Le, &Ee);
1653: if (Ea%Ee)
1654: {
1655: if (DEBUGLEVEL >= 5)
1656: fprintferr(" Increasing Ea\n");
1657: pie = redelt(pie, p, pmf);
1658: /* we compute a new element such E = lcm(Ea, Ee) */
1659: w = testc2(p, chi, pmr, pmf, nu, Ea, pie, Ee, ns);
1660: if (gcmp1((GEN)w[1]))
1661: {
1662: /* there are at least 2 factors mod. p => chi can be split */
1663: phi = eleval(fx, (GEN)w[2], alph);
1664: phi = redelt(phi, p, pmf);
1665: if (flag) mf += 3;
1666: return Decomp(p, fx, mf, phi, (GEN)w[3], (GEN)w[4], flag);
1667: }
1668: break;
1669: }
1670:
1671: if (eq) delt = gmul(delt, gpowgs(p, eq));
1672: if (er) delt = gmul(delt, gpowgs(nu, er));
1673: beta = gsub(beta, delt);
1674: }
1675:
1676: /* we replace alpha by a new alpha with a larger F or E */
1677: alph = eleval(fx, (GEN)w[2], alph);
1678: chi = (GEN)w[3];
1679: nu = (GEN)w[4];
1680:
1681: w = update_alpha(p, fx, alph, chi, pmr, pmf, mf, ns);
1682: alph = (GEN)w[1];
1683: chi = (GEN)w[2];
1684: pmr = (GEN)w[3];
1685: pdr = (GEN)w[4];
1686:
1687: /* that can happen if p does not divide the field discriminant! */
1688: if (is_pm1(pmr))
1689: {
1690: if (flag)
1691: {
1692: p1 = lift((GEN)factmod(chi, p)[1]);
1693: l = lg(p1) - 1;
1694: if (l == 1) return NULL;
1695: phi = redelt(alph, p, pmf);
1696: return Decomp(p, fx, ggval(pmf, p), phi, chi, (GEN)p1[l], flag);
1697: }
1698: else
1699: return dbasis(p, fx, mf, alph, chi);
1700: }
1701: }
1702: }
1703:
1704: /* Returns [1,phi,chi,nu] if phi non-primary
1705: * [2,phi,chi,nu] with F_phi = lcm (F_alpha, F_theta)
1706: * and E_phi = E_alpha
1707: */
1708: static GEN
1709: testb2(GEN p, GEN fa, long Fa, GEN theta, GEN pmf, long Ft, GEN ns)
1710: {
1711: long m, Dat, t, v = varn(fa);
1712: GEN b, w, phi, h;
1713:
1714: Dat = clcm(Fa, Ft) + 3;
1715: b = cgetg(5, t_VEC);
1716: m = p[2];
1717: if (degpol(p) > 0 || m < 0) m = 0;
1718:
1719: for (t = 1;; t++)
1720: {
1721: h = m? stopoly(t, m, v): scalarpol(stoi(t), v);
1722: phi = gadd(theta, gmod(h, fa));
1723: w = factcp(p, fa, phi, pmf, ns);
1724: h = (GEN)w[3];
1725: if (h[2] > 1) { b[1] = un; break; }
1726: if (lgef(w[2]) == Dat) { b[1] = deux; break; }
1727: }
1728:
1729: b[2] = (long)phi;
1730: b[3] = w[1];
1731: b[4] = w[2];
1732: return b;
1733: }
1734:
1735: /* Returns [1, phi, chi, nu] if phi non-primary
1736: * [2, phi, chi, nu] if E_phi = lcm (E_alpha, E_theta)
1737: */
1738: static GEN
1739: testc2(GEN p, GEN fa, GEN pmr, GEN pmf, GEN alph2, long Ea, GEN thet2,
1740: long Et, GEN ns)
1741: {
1742: GEN b, c1, c2, c3, psi, phi, w, h;
1743: long r, s, t, v = varn(fa);
1744:
1745: b=cgetg(5, t_VEC);
1746:
1747: cbezout(Ea, Et, &r, &s); t = 0;
1748: while (r < 0) { r = r + Et; t++; }
1749: while (s < 0) { s = s + Ea; t++; }
1750:
1751: c1 = lift_intern(gpuigs(gmodulcp(alph2, fa), s));
1752: c2 = lift_intern(gpuigs(gmodulcp(thet2, fa), r));
1753: c3 = gdiv(gmod(gmul(c1, c2), fa), gpuigs(p, t));
1754:
1755: psi = redelt(c3, pmr, pmr);
1756: phi = gadd(polx[v], psi);
1757:
1758: w = factcp(p, fa, phi, pmf, ns);
1759: h = (GEN)w[3];
1760: b[1] = (h[2] > 1)? un: deux;
1761: b[2] = (long)phi;
1762: b[3] = w[1];
1763: b[4] = w[2];
1764: return b;
1765: }
1766:
1767: /* evaluate h(a) mod f */
1768: GEN
1769: eleval(GEN f,GEN h,GEN a)
1770: {
1771: long n,av,tetpil;
1772: GEN y;
1773:
1774: if (typ(h) != t_POL) return gcopy(h);
1775: av = tetpil = avma;
1776: n=lgef(h)-1; y=(GEN)h[n];
1777: for (n--; n>=2; n--)
1778: {
1779: y = gadd(gmul(y,a),(GEN)h[n]);
1780: tetpil=avma; y = gmod(y,f);
1781: }
1782: return gerepile(av,tetpil,y);
1783: }
1784:
1785: GEN addshiftw(GEN x, GEN y, long d);
1786:
1787: static GEN
1788: shiftpol(GEN x, long v)
1789: {
1790: x = addshiftw(x, zeropol(v), 1);
1791: setvarn(x,v); return x;
1792: }
1793:
1794: /* Sylvester's matrix, mod p^m (assumes f1 monic) */
1795: static GEN
1796: sylpm(GEN f1,GEN f2,GEN pm)
1797: {
1798: long n,j,v=varn(f1);
1799: GEN a,h;
1800:
1801: n=degpol(f1); a=cgetg(n+1,t_MAT);
1802: h = FpX_res(f2,f1,pm);
1803: for (j=1;; j++)
1804: {
1805: a[j] = (long)pol_to_vec(h,n);
1806: if (j == n) break;
1807: h = FpX_res(shiftpol(h,v),f1,pm);
1808: }
1809: return hnfmodid(a,pm);
1810: }
1811:
1812: /* polynomial gcd mod p^m (assumes f1 monic) */
1813: GEN
1814: gcdpm(GEN f1,GEN f2,GEN pm)
1815: {
1816: long n,c,v=varn(f1),av=avma,tetpil;
1817: GEN a,col;
1818:
1819: n=degpol(f1); a=sylpm(f1,f2,pm);
1820: for (c=1; c<=n; c++)
1821: if (signe(resii(gcoeff(a,c,c),pm))) break;
1822: if (c > n) { avma=av; return zeropol(v); }
1823:
1824: col = gdiv((GEN)a[c], gcoeff(a,c,c)); tetpil=avma;
1825: return gerepile(av,tetpil, gtopolyrev(col,v));
1826: }
1827:
1828: /* reduced resultant mod p^m (assumes x monic) */
1829: GEN
1830: respm(GEN x,GEN y,GEN pm)
1831: {
1832: long av = avma;
1833: GEN p1 = sylpm(x,y,pm);
1834:
1835: p1 = gcoeff(p1,1,1);
1836: if (egalii(p1,pm)) { avma = av; return gzero; }
1837: return gerepileuptoint(av, icopy(p1));
1838: }
1839:
1840: /* Normalized integral basis */
1841: static GEN
1842: nbasis(GEN ibas,GEN pd)
1843: {
1844: long k, n = lg(ibas)-1;
1845: GEN a = cgetg(n+1,t_MAT);
1846: for (k=1; k<=n; k++) a[k] = (long)pol_to_vec((GEN)ibas[k],n);
1847: return gdiv(hnfmodid(a,pd), pd);
1848: }
1849:
1850: /*******************************************************************/
1851: /* */
1852: /* BUCHMANN-LENSTRA ALGORITHM */
1853: /* */
1854: /*******************************************************************/
1855: static GEN lens(GEN nf,GEN p,GEN a);
1856: GEN element_powid_mod_p(GEN nf, long I, GEN n, GEN p);
1857:
1858: /* return a Z basis of Z_K's p-radical, modfrob = x--> x^p-x */
1859: static GEN
1860: pradical(GEN nf, GEN p, GEN *modfrob)
1861: {
1862: long i,N = degpol(nf[1]);
1863: GEN p1,m,frob,rad;
1864:
1865: frob = cgetg(N+1,t_MAT);
1866: for (i=1; i<=N; i++)
1867: frob[i] = (long) element_powid_mod_p(nf,i,p, p);
1868:
1869: /* p1 = smallest power of p st p^k >= N */
1870: p1=p; while (cmpis(p1,N)<0) p1=mulii(p1,p);
1871: if (p1==p) m = frob;
1872: else
1873: {
1874: m=cgetg(N+1,t_MAT); p1 = divii(p1,p);
1875: for (i=1; i<=N; i++)
1876: m[i]=(long)element_pow_mod_p(nf,(GEN)frob[i],p1, p);
1877: }
1878: rad = FpM_ker(m, p);
1879: for (i=1; i<=N; i++)
1880: coeff(frob,i,i) = lsubis(gcoeff(frob,i,i), 1);
1881: *modfrob = frob; return rad;
1882: }
1883:
1884: static GEN
1885: project(GEN algebre, GEN x, long k, long kbar)
1886: {
1887: x = inverseimage(algebre,x);
1888: x += k; x[0] = evaltyp(t_COL) | evallg(kbar+1);
1889: return x;
1890: }
1891:
1892: /* Calcule le polynome minimal de alpha dans algebre (coeffs dans Z) */
1893: static GEN
1894: pol_min(GEN alpha,GEN nf,GEN algebre,long kbar,GEN p)
1895: {
1896: long av=avma,tetpil,i,N,k;
1897: GEN p1,puiss;
1898:
1899: N = lg(nf[1])-3; puiss=cgetg(N+2,t_MAT);
1900: k = N-kbar; p1=alpha;
1901: puiss[1] = (long)gscalcol_i(gun,kbar);
1902: for (i=2; i<=N+1; i++)
1903: {
1904: if (i>2) p1 = element_mul(nf,p1,alpha);
1905: puiss[i] = (long) project(algebre,p1,k,kbar);
1906: }
1907: puiss = lift_intern(puiss);
1908: p1 = (GEN)FpM_ker(puiss, p)[1]; tetpil=avma;
1909: return gerepile(av,tetpil,gtopolyrev(p1,0));
1910: }
1911:
1912: /* Evalue le polynome pol en alpha,element de nf */
1913: static GEN
1914: eval_pol(GEN nf,GEN pol,GEN alpha,GEN algebre,GEN algebre1)
1915: {
1916: long av=avma,tetpil,i,kbar,k, lx = lgef(pol)-1, N = degpol(nf[1]);
1917: GEN res;
1918:
1919: kbar = lg(algebre1)-1; k = N-kbar;
1920: res = gscalcol_i((GEN)pol[lx], N);
1921: for (i=2; i<lx; i++)
1922: {
1923: res = element_mul(nf,alpha,res);
1924: res[1] = ladd((GEN)res[1],(GEN)pol[i]);
1925: }
1926: res = project(algebre,res,k,kbar); tetpil=avma;
1927: return gerepile(av,tetpil,gmul(algebre1,res));
1928: }
1929:
1930: static GEN
1931: kerlens2(GEN x, GEN p)
1932: {
1933: long i,j,k,t,nbc,nbl,av,av1;
1934: GEN a,c,l,d,y,q;
1935:
1936: av=avma; a=gmul(x,gmodulsg(1,p));
1937: nbl=nbc=lg(x)-1;
1938: c=new_chunk(nbl+1); for (i=1; i<=nbl; i++) c[i]=0;
1939: l=new_chunk(nbc+1);
1940: d=new_chunk(nbc+1);
1941: k = t = 1;
1942: while (t<=nbl && k<=nbc)
1943: {
1944: for (j=1; j<k; j++)
1945: for (i=1; i<=nbl; i++)
1946: if (i!=l[j])
1947: coeff(a,i,k) = lsub(gmul((GEN)d[j],gcoeff(a,i,k)),
1948: gmul(gcoeff(a,l[j],k),gcoeff(a,i,j)));
1949: t=1; while (t<=nbl && (c[t] || gcmp0(gcoeff(a,t,k)))) t++;
1950: if (t<=nbl) { d[k]=coeff(a,t,k); c[t]=k; l[k]=t; k++; }
1951: }
1952: if (k>nbc) err(bugparier,"kerlens2");
1953: y=cgetg(nbc+1,t_COL);
1954: y[1]=(k>1)?coeff(a,l[1],k):un;
1955: for (q=gun,j=2; j<k; j++)
1956: {
1957: q=gmul(q,(GEN)d[j-1]);
1958: y[j]=lmul(gcoeff(a,l[j],k),q);
1959: }
1960: if (k>1) y[k]=lneg(gmul(q,(GEN)d[k-1]));
1961: for (j=k+1; j<=nbc; j++) y[j]=zero;
1962: av1=avma; return gerepile(av,av1,lift(y));
1963: }
1964:
1965: static GEN
1966: kerlens(GEN x, GEN pgen)
1967: {
1968: long av = avma, i,j,k,t,nbc,nbl,p,q,*c,*l,*d,**a;
1969: GEN y;
1970:
1971: if (cmpis(pgen, MAXHALFULONG>>1) > 0)
1972: return kerlens2(x,pgen);
1973: /* ici p <= (MAXHALFULONG>>1) ==> long du C */
1974: p=itos(pgen); nbl=nbc=lg(x)-1;
1975: a=(long**)new_chunk(nbc+1);
1976: for (j=1; j<=nbc; j++)
1977: {
1978: c=a[j]=new_chunk(nbl+1);
1979: for (i=1; i<=nbl; i++) c[i]=smodis(gcoeff(x,i,j),p);
1980: }
1981: c=new_chunk(nbl+1); for (i=1; i<=nbl; i++) c[i]=0;
1982: l=new_chunk(nbc+1);
1983: d=new_chunk(nbc+1);
1984: k = t = 1;
1985: while (t<=nbl && k<=nbc)
1986: {
1987: for (j=1; j<k; j++)
1988: for (i=1; i<=nbl; i++)
1989: if (i!=l[j])
1990: a[k][i] = (d[j]*a[k][i] - a[j][i]*a[k][l[j]]) % p;
1991: t=1; while (t<=nbl && (c[t] || !a[k][t])) t++;
1992: if (t<=nbl) { d[k]=a[k][t]; c[t]=k; l[k++]=t; }
1993: }
1994: if (k>nbc) err(bugparier,"kerlens");
1995: avma=av; y=cgetg(nbc+1,t_COL);
1996: t=(k>1) ? a[k][l[1]]:1;
1997: y[1]=(t>0)? lstoi(t):lstoi(t+p);
1998: for (q=1,j=2; j<k; j++)
1999: {
2000: q = (q*d[j-1]) % p;
2001: t = (a[k][l[j]]*q) % p;
2002: y[j] = (t>0) ? lstoi(t) : lstoi(t+p);
2003: }
2004: if (k>1)
2005: {
2006: t = (q*d[k-1]) % p;
2007: y[k] = (t>0) ? lstoi(p-t) : lstoi(-t);
2008: }
2009: for (j=k+1; j<=nbc; j++) y[j]=zero;
2010: return y;
2011: }
2012:
2013: /* Calcule la constante de lenstra de l'ideal p.Z_K+a.Z_K ou a est un
2014: vecteur sur la base d'entiers */
2015: static GEN
2016: lens(GEN nf, GEN p, GEN a)
2017: {
2018: long av=avma,tetpil,N=degpol(nf[1]),j;
2019: GEN mat=cgetg(N+1,t_MAT);
2020: for (j=1; j<=N; j++) mat[j]=(long)element_mulid(nf,a,j);
2021: tetpil=avma; return gerepile(av,tetpil,kerlens(mat,p));
2022: }
2023:
2024: extern GEN det_mod_P_n(GEN a, GEN N, GEN P);
2025: GEN sylvestermatrix_i(GEN x, GEN y);
2026:
2027: /* check if p^va doesnt divide norm x (or norm(x+p)) */
2028: #if 0
2029: /* compute norm mod p^whatneeded using Sylvester's matrix */
2030: /* looks slower than the new subresultant. Have to re-check this */
2031: static GEN
2032: prime_check_elt(GEN a, GEN pol, GEN p, GEN pf)
2033: {
2034: GEN M,mod,x, c = denom(content(a));
2035: long v = pvaluation(c, p, &x); /* x is junk */
2036:
2037: mod = mulii(pf, gpowgs(p, degpol(pol)*v + 1));
2038:
2039: x = FpX_red(gmul(a,c), mod);
2040: M = sylvestermatrix_i(pol,x);
2041: if (det_mod_P_n(M,mod,p) == gzero)
2042: {
2043: x[2] = ladd((GEN)x[2], mulii(p,c));
2044: M = sylvestermatrix_i(pol,x);
2045: if (det_mod_P_n(M,mod,p) == gzero) return NULL;
2046: a[2] = ladd((GEN)a[2], p);
2047: }
2048: return a;
2049: }
2050: #else
2051: /* use subres to compute norm */
2052: static GEN
2053: prime_check_elt(GEN a, GEN pol, GEN p, GEN pf)
2054: {
2055: GEN norme=subres(pol,a);
2056: if (resii(divii(norme,pf),p) != gzero) return a;
2057: /* Note: a+p can't succeed if e > 1, can we know this at this point ? */
2058: a=gadd(a,p); norme=subres(pol,a);
2059: if (resii(divii(norme,pf),p) != gzero) return a;
2060: return NULL;
2061: }
2062: #endif
2063:
2064: #if 0
2065: static GEN
2066: prime_two_elt_loop(GEN beta, GEN pol, GEN p, GEN pf)
2067: {
2068: long av, m = lg(beta)-1;
2069: int i,j,K, *x = (int*)new_chunk(m+1);
2070: GEN a;
2071:
2072: K = 1; av = avma;
2073: for(;;)
2074: { /* x runs through strictly increasing sequences of length K,
2075: * 1 <= x[i] <= m */
2076: nextK:
2077: if (DEBUGLEVEL) fprintferr("K = %d\n", K);
2078: for (i=1; i<=K; i++) x[i] = i;
2079: for(;;)
2080: {
2081: if (DEBUGLEVEL > 1)
2082: {
2083: for (i=1; i<=K; i++) fprintferr("%d ",x[i]);
2084: fprintferr("\n"); flusherr();
2085: }
2086: a = (GEN)beta[x[1]];
2087: for (i=2; i<=K; i++) a = gadd(a, (GEN)beta[x[i]]);
2088: if ((a = prime_check_elt(a,pol,p,pf))) return a;
2089: avma = av;
2090:
2091: /* start: i = K+1; */
2092: do
2093: {
2094: if (--i == 0)
2095: {
2096: if (++K > m) return NULL; /* fail */
2097: goto nextK;
2098: }
2099: x[i]++;
2100: } while (x[i] > m - K + i);
2101: for (j=i; j<K; j++) x[j+1] = x[j]+1;
2102: }
2103: }
2104: }
2105: #endif
2106:
2107: static GEN
2108: random_prime_two_elt_loop(GEN beta, GEN pol, GEN p, GEN pf)
2109: {
2110: long av = avma, z,i, m = lg(beta)-1;
2111: long keep = getrand();
2112: int c = 0;
2113: GEN a;
2114:
2115: for(i=1; i<=m; i++)
2116: if ((a = prime_check_elt((GEN)beta[i],pol,p,pf))) return a;
2117: (void)setrand(1);
2118: if (DEBUGLEVEL) fprintferr("prime_two_elt_loop, hard case: ");
2119: for(;;avma=av)
2120: {
2121: if (DEBUGLEVEL) fprintferr("%d ", ++c);
2122: a = gzero;
2123: for (i=1; i<=m; i++)
2124: {
2125: z = mymyrand() >> (BITS_IN_RANDOM-5); /* in [0,15] */
2126: if (z >= 9) z -= 7;
2127: a = gadd(a,gmulsg(z,(GEN)beta[i]));
2128: }
2129: if ((a = prime_check_elt(a,pol,p,pf)))
2130: {
2131: if (DEBUGLEVEL) fprintferr("\n");
2132: (void)setrand(keep); return a;
2133: }
2134: }
2135: }
2136:
2137: /* Input: an ideal mod p (!= Z_K)
2138: * Output: a 2-elt representation [p, x] */
2139: static GEN
2140: prime_two_elt(GEN nf, GEN p, GEN ideal)
2141: {
2142: GEN beta,a,pf, pol = (GEN)nf[1];
2143: long f, N=degpol(pol), m=lg(ideal)-1;
2144: ulong av;
2145:
2146: if (!m) return gscalcol_i(p,N);
2147:
2148: /* we want v_p(Norm(beta)) = p^f, f = N-m */
2149: av = avma; f = N-m; pf = gpuigs(p,f);
2150: ideal = centerlift(ideal);
2151: ideal = concatsp(gscalcol(p,N), ideal);
2152: ideal = ideal_better_basis(nf, ideal, p);
2153: beta = gmul((GEN)nf[7], ideal);
2154:
2155: #if 0
2156: a = prime_two_elt_loop(beta,pol,p,pf);
2157: if (!a) err(bugparier, "prime_two_elt (failed)");
2158: #else
2159: a = random_prime_two_elt_loop(beta,pol,p,pf);
2160: #endif
2161:
2162: a = centermod(algtobasis_intern(nf,a), p);
2163: if (resii(divii(subres(gmul((GEN)nf[7],a),pol),pf),p) == gzero)
2164: a[1] = laddii((GEN)a[1],p);
2165: return gerepilecopy(av,a);
2166: }
2167:
2168: static GEN
2169: apply_kummer(GEN nf,GEN pol,GEN e,GEN p,long N,GEN *beta)
2170: {
2171: GEN T,p1, res = cgetg(6,t_VEC);
2172: long f = degpol(pol);
2173:
2174: res[1]=(long)p;
2175: res[3]=(long)e;
2176: res[4]=lstoi(f);
2177: if (f == N) /* inert */
2178: {
2179: res[2]=(long)gscalcol_i(p,N);
2180: res[5]=(long)gscalcol_i(gun,N);
2181: }
2182: else
2183: {
2184: T = (GEN) nf[1];
2185: if (ggval(subres(pol,T),p) > f)
2186: pol[2] = laddii((GEN)pol[2],p);
2187: res[2] = (long) algtobasis_intern(nf,pol);
2188:
2189: p1 = FpX_div(T,pol,p);
2190: res[5] = (long) centermod(algtobasis_intern(nf,p1), p);
2191:
2192: if (beta)
2193: *beta = *beta? FpX_div(*beta,pol,p): p1;
2194: }
2195: return res;
2196: }
2197:
2198: /* prime ideal decomposition of p sorted by increasing residual degree */
2199: GEN
2200: primedec(GEN nf, GEN p)
2201: {
2202: long av=avma,tetpil,i,j,k,kbar,np,c,indice,N,lp;
2203: GEN ex,f,list,ip,elth,h;
2204: GEN modfrob,algebre,algebre1,b,mat1,T;
2205: GEN alpha,beta,p1,p2,unmodp,zmodp,idmodp;
2206:
2207: if (DEBUGLEVEL>=3) timer2();
2208: nf=checknf(nf); T=(GEN)nf[1]; N=degpol(T);
2209: f=factmod(T,p); ex=(GEN)f[2];
2210: f=centerlift((GEN)f[1]); np=lg(f);
2211: if (DEBUGLEVEL>=6) msgtimer("factmod");
2212:
2213: if (signe(modii((GEN)nf[4],p))) /* p doesn't divide index */
2214: {
2215: list=cgetg(np,t_VEC);
2216: for (i=1; i<np; i++)
2217: list[i]=(long)apply_kummer(nf,(GEN)f[i],(GEN)ex[i],p,N, NULL);
2218: if (DEBUGLEVEL>=6) msgtimer("simple primedec");
2219: p1=stoi(4); tetpil=avma;
2220: return gerepile(av,tetpil,vecsort(list,p1));
2221: }
2222:
2223: p1 = (GEN)f[1];
2224: for (i=2; i<np; i++)
2225: p1 = FpX_red(gmul(p1, (GEN)f[i]), p);
2226: p1 = FpX_red(gdiv(gadd(gmul(p1, FpX_div(T,p1,p)), gneg(T)), p), p);
2227: list = cgetg(N+1,t_VEC);
2228: indice=1; beta=NULL;
2229: for (i=1; i<np; i++) /* e = 1 or f[i] does not divide p1 (mod p) */
2230: if (is_pm1(ex[i]) || signe(FpX_res(p1,(GEN)f[i],p)))
2231: list[indice++] = (long)apply_kummer(nf,(GEN)f[i],(GEN)ex[i],p,N,&beta);
2232: if (DEBUGLEVEL>=3) msgtimer("unramified factors");
2233:
2234: /* modfrob = modified Frobenius: x -> x^p - x mod p */
2235: ip = pradical(nf,p,&modfrob);
2236: if (DEBUGLEVEL>=3) msgtimer("pradical");
2237:
2238: if (beta)
2239: {
2240: beta = algtobasis_intern(nf,beta);
2241: lp=lg(ip)-1; p1=cgetg(2*lp+N+1,t_MAT);
2242: for (i=1; i<=N; i++) p1[i]=(long)element_mulid(nf,beta,i);
2243: for ( ; i<=N+lp; i++)
2244: {
2245: p2 = (GEN) ip[i-N];
2246: p1[i+lp] = (long) p2;
2247: p1[i] = ldiv(element_mul(nf,lift(p2),beta),p);
2248: }
2249: ip = FpM_image(p1, p);
2250: if (lg(ip)>N) err(bugparier,"primedec (bad pradical)");
2251: }
2252: unmodp=gmodulsg(1,p); zmodp=gmodulsg(0,p);
2253: idmodp = idmat_intern(N,unmodp,zmodp);
2254: ip = gmul(ip, unmodp);
2255:
2256: h=cgetg(N+1,t_VEC); h[1]=(long)ip;
2257: for (c=1; c; c--)
2258: {
2259: elth=(GEN)h[c]; k=lg(elth)-1; kbar=N-k;
2260: p1 = concatsp(elth,(GEN)idmodp[1]);
2261: algebre = suppl_intern(p1,idmodp);
2262: algebre1 = cgetg(kbar+1,t_MAT);
2263: for (i=1; i<=kbar; i++) algebre1[i]=algebre[i+k];
2264: b = gmul(modfrob,algebre1);
2265: for (i=1;i<=kbar;i++)
2266: b[i] = (long) project(algebre,(GEN) b[i],k,kbar);
2267: mat1 = FpM_ker(lift_intern(b), p);
2268: if (lg(mat1)>2)
2269: {
2270: GEN mat2 = cgetg(k+N+1,t_MAT);
2271: for (i=1; i<=k; i++) mat2[i]=elth[i];
2272: alpha=gmul(algebre1,(GEN)mat1[2]);
2273: p1 = pol_min(alpha,nf,algebre,kbar,p);
2274: p1 = (GEN)factmod(p1,p)[1];
2275: for (i=1; i<lg(p1); i++)
2276: {
2277: beta = eval_pol(nf,(GEN)p1[i],alpha,algebre,algebre1);
2278: beta = lift_intern(beta);
2279: for (j=1; j<=N; j++)
2280: mat2[k+j] = (long)FpV(element_mulid(nf,beta,j), p);
2281: h[c] = (long) image(mat2); c++;
2282: }
2283: }
2284: else
2285: {
2286: long av1; p1 = cgetg(6,t_VEC);
2287: list[indice++] = (long)p1;
2288: p1[1]=(long)p; p1[4]=lstoi(kbar);
2289: p1[2]=(long)prime_two_elt(nf,p,elth);
2290: p1[5]=(long)lens(nf,p,(GEN)p1[2]);
2291: av1=avma;
2292: i = int_elt_val(nf,(GEN)p1[5],p,(GEN)p1[5],NULL,N-1);
2293: avma=av1;
2294: p1[3]=lstoi(i+1);
2295: }
2296: if (DEBUGLEVEL>=3) msgtimer("h[%ld]",c);
2297: }
2298: setlg(list, indice); tetpil=avma;
2299: return gerepile(av,tetpil,gen_sort(list,0,cmp_prime_over_p));
2300: }
2301:
2302: /* REDUCTION Modulo a prime ideal */
2303:
2304: /* x integral, reduce mod prh in HNF */
2305: GEN
2306: nfreducemodpr_i(GEN x, GEN prh)
2307: {
2308: GEN p = gcoeff(prh,1,1);
2309: long i,j;
2310:
2311: x = dummycopy(x);
2312: for (i=lg(x)-1; i>=2; i--)
2313: {
2314: GEN t = (GEN)prh[i], p1 = resii((GEN)x[i], p);
2315: x[i] = (long)p1;
2316: if (signe(p1) && is_pm1(t[i]))
2317: {
2318: for (j=1; j<i; j++)
2319: x[j] = lsubii((GEN)x[j], mulii(p1, (GEN)t[j]));
2320: x[i] = zero;
2321: }
2322: }
2323: x[1] = lresii((GEN)x[1], p); return x;
2324: }
2325:
2326: /* for internal use */
2327: GEN
2328: nfreducemodpr(GEN nf, GEN x, GEN prhall)
2329: {
2330: long i,v;
2331: GEN p,prh,den;
2332:
2333: for (i=lg(x)-1; i>0; i--)
2334: if (typ(x[i]) == t_INTMOD) { x=lift_intern(x); break; }
2335: prh=(GEN)prhall[1]; p=gcoeff(prh,1,1);
2336: den=denom(x);
2337: if (!gcmp1(den))
2338: {
2339: v=ggval(den,p);
2340: if (v) x=element_mul(nf,x,element_pow(nf,(GEN)prhall[2],stoi(v)));
2341: x = gmod(x,p);
2342: }
2343: return FpV(nfreducemodpr_i(x, prh), p);
2344: }
2345:
2346: /* public function */
2347: GEN
2348: nfreducemodpr2(GEN nf, GEN x, GEN prhall)
2349: {
2350: long av = avma; checkprhall(prhall);
2351: if (typ(x) != t_COL) x = algtobasis(nf,x);
2352: return gerepileupto(av, nfreducemodpr(nf,x,prhall));
2353: }
2354:
2355: /* relative ROUND 2
2356: *
2357: * input: nf = base field K
2358: * x monic polynomial, coefficients in Z_K, degree n defining a relative
2359: * extension L=K(theta).
2360: * One MUST have varn(x) < varn(nf[1]).
2361: * output: a pseudo-basis [A,I] of Z_L, where A is in M_n(K) in HNF form and
2362: * I a vector of n ideals.
2363: */
2364:
2365: /* given MODULES x and y by their pseudo-bases in HNF, gives a
2366: * pseudo-basis of the module generated by x and y. For internal use.
2367: */
2368: static GEN
2369: rnfjoinmodules(GEN nf, GEN x, GEN y)
2370: {
2371: long i,lx,ly;
2372: GEN p1,p2,z,Hx,Hy,Ix,Iy;
2373:
2374: if (x == NULL) return y;
2375: Hx=(GEN)x[1]; lx=lg(Hx); Ix=(GEN)x[2];
2376: Hy=(GEN)y[1]; ly=lg(Hy); Iy=(GEN)y[2];
2377: i = lx+ly-1;
2378: z = (GEN)gpmalloc(sizeof(long*)*(3+2*i));
2379: *z = evaltyp(t_VEC)|evallg(3);
2380: p1 = z+3; z[1]=(long)p1; *p1 = evaltyp(t_MAT)|evallg(i);
2381: p2 = p1+i; z[2]=(long)p2; *p2 = evaltyp(t_VEC)|evallg(i);
2382:
2383: for (i=1; i<lx; i++) { p1[i]=Hx[i]; p2[i]=Ix[i]; }
2384: for ( ; i<lx+ly-1; i++) { p1[i]=Hy[i-lx+1]; p2[i]=Iy[i-lx+1]; }
2385: x = nfhermite(nf,z); free(z); return x;
2386: }
2387:
2388: /* a usage interne, pas de gestion de pile : x et y sont des vecteurs dont
2389: * les coefficients sont les composantes sur nf[7]; avec reduction mod pr sauf
2390: * si prhall=NULL
2391: */
2392: static GEN
2393: rnfelement_mulidmod(GEN nf, GEN multab, GEN unnf, GEN x, long h, GEN prhall)
2394: {
2395: long j,k,N;
2396: GEN p1,c,v,s,znf;
2397:
2398: if (h==1) return gcopy(x);
2399: N = lg(x)-1; multab += (h-1)*N;
2400: x = lift(x); v = cgetg(N+1,t_COL);
2401: znf = gscalcol_i(gzero,lg(unnf)-1);
2402: for (k=1; k<=N; k++)
2403: {
2404: s = gzero;
2405: for (j=1; j<=N; j++)
2406: if (!gcmp0(p1 = (GEN)x[j]) && !gcmp0(c = gcoeff(multab,k,j)))
2407: {
2408: if (!gegal(c,unnf)) p1 = element_mul(nf,p1,c);
2409: s = gadd(s,p1);
2410: }
2411: if (s == gzero) s = znf;
2412: else
2413: if (prhall) s = nfreducemodpr(nf,s,prhall);
2414: v[k] = (long)s;
2415: }
2416: return v;
2417: }
2418:
2419: /* a usage interne, pas de gestion de pile : x est un vecteur dont
2420: * les coefficients sont les composantes sur nf[7]
2421: */
2422: static GEN
2423: rnfelement_sqrmod(GEN nf, GEN multab, GEN unnf, GEN x, GEN prhall)
2424: {
2425: long i,j,k,n;
2426: GEN p1,c,z,s;
2427:
2428: n=lg(x)-1; x=lift(x); z=cgetg(n+1,t_COL);
2429: for (k=1; k<=n; k++)
2430: {
2431: if (k == 1)
2432: s = element_sqr(nf,(GEN)x[1]);
2433: else
2434: s = gmul2n(element_mul(nf,(GEN)x[1],(GEN)x[k]), 1);
2435: for (i=2; i<=n; i++)
2436: {
2437: c = gcoeff(multab,k,(i-1)*n+i);
2438: if (!gcmp0(c))
2439: {
2440: p1=element_sqr(nf,(GEN)x[i]);
2441: if (!gegal(c,unnf)) p1 = element_mul(nf,p1,c);
2442: s = gadd(s,p1);
2443: }
2444: for (j=i+1; j<=n; j++)
2445: {
2446: c = gcoeff(multab,k,(i-1)*n+j);
2447: if (!gcmp0(c))
2448: {
2449: p1=gmul2n(element_mul(nf,(GEN)x[i],(GEN)x[j]),1);
2450: if (!gegal(c,unnf)) p1 = element_mul(nf,p1,c);
2451: s = gadd(s,p1);
2452: }
2453: }
2454: }
2455: if (prhall) s = nfreducemodpr(nf,s,prhall);
2456: z[k]=(long)s;
2457: }
2458: return z;
2459: }
2460:
2461: /* Compute x^n mod pr in the extension, assume n >= 0 [cf puissii()]*/
2462: static GEN
2463: rnfelementid_powmod(GEN nf, GEN multab, GEN matId, long h, GEN n, GEN prhall)
2464: {
2465: ulong av = avma;
2466: long i,k,m;
2467: GEN y,p1, unrnf=(GEN)matId[1], unnf=(GEN)unrnf[1];
2468:
2469: if (!signe(n)) return unrnf;
2470: y = (GEN)matId[h]; p1 = n+2; m = *p1;
2471: k = 1+bfffo(m); m<<=k; k = BITS_IN_LONG-k;
2472: for (i=lgefint(n)-2;;)
2473: {
2474: for (; k; m<<=1,k--)
2475: {
2476: y = rnfelement_sqrmod(nf,multab,unnf,y,prhall);
2477: if (m < 0) y = rnfelement_mulidmod(nf,multab,unnf,y,h,prhall);
2478: }
2479: if (--i == 0) break;
2480: m = *++p1; k = BITS_IN_LONG;
2481: }
2482: return gerepilecopy(av, y);
2483: }
2484:
2485: GEN
2486: mymod(GEN x, GEN p)
2487: {
2488: long i,lx, tx = typ(x);
2489: GEN y,p1;
2490:
2491: if (tx == t_INT) return resii(x,p);
2492: if (tx == t_FRAC)
2493: {
2494: p1 = resii((GEN)x[2], p);
2495: if (p1 != gzero) { cgiv(p1); return gmod(x,p); }
2496: return x;
2497: }
2498: if (!is_matvec_t(tx))
2499: err(bugparier, "mymod (missing type)");
2500: lx = lg(x); y = cgetg(lx,tx);
2501: for (i=1; i<lx; i++) y[i] = (long)mymod((GEN)x[i],p);
2502: return y;
2503: }
2504:
2505: static GEN
2506: rnfordmax(GEN nf, GEN pol, GEN pr, GEN unnf, GEN id, GEN matId)
2507: {
2508: long av=avma,tetpil,av1,lim,i,j,k,n,v1,v2,vpol,m,cmpt,sep;
2509: GEN p,q,q1,prhall,A,Aa,Aaa,A1,I,R,p1,p2,p3,multab,multabmod,Aainv;
2510: GEN pip,baseIp,baseOp,alpha,matprod,alphainv,matC,matG,vecpro,matH;
2511: GEN neworder,H,Hid,alphalistinv,alphalist,prhinv;
2512:
2513: if (DEBUGLEVEL>1) fprintferr(" treating %Z\n",pr);
2514: prhall=nfmodprinit(nf,pr);
2515: q=cgetg(3,t_VEC); q[1]=(long)pr;q[2]=(long)prhall;
2516: p1=rnfdedekind(nf,pol,q);
2517: if (gcmp1((GEN)p1[1])) return gerepilecopy(av,(GEN)p1[2]);
2518:
2519: sep=itos((GEN)p1[3]);
2520: A=gmael(p1,2,1);
2521: I=gmael(p1,2,2);
2522:
2523: n=degpol(pol); vpol=varn(pol);
2524: p=(GEN)pr[1]; q=powgi(p,(GEN)pr[4]); pip=(GEN)pr[2];
2525: q1=q; while (cmpis(q1,n)<0) q1=mulii(q1,q);
2526:
2527: multab=cgetg(n*n+1,t_MAT);
2528: for (j=1; j<=n*n; j++) multab[j]=lgetg(n+1,t_COL);
2529: prhinv = idealinv(nf,(GEN)prhall[1]);
2530: alphalistinv=cgetg(n+1,t_VEC);
2531: alphalist=cgetg(n+1,t_VEC);
2532: A1=cgetg(n+1,t_MAT);
2533: av1=avma; lim=stack_lim(av1,1);
2534: for(cmpt=1; ; cmpt++)
2535: {
2536: if (DEBUGLEVEL>1)
2537: {
2538: fprintferr(" %ld%s pass\n",cmpt,eng_ord(cmpt));
2539: flusherr();
2540: }
2541: for (i=1; i<=n; i++)
2542: {
2543: if (gegal((GEN)I[i],id)) alphalist[i]=alphalistinv[i]=(long)unnf;
2544: else
2545: {
2546: p1=ideal_two_elt(nf,(GEN)I[i]);
2547: v1=gcmp0((GEN)p1[1])? EXP220
2548: : element_val(nf,(GEN)p1[1],pr);
2549: v2=element_val(nf,(GEN)p1[2],pr);
2550: if (v1>v2) p2=(GEN)p1[2]; else p2=(GEN)p1[1];
2551: alphalist[i]=(long)p2;
2552: alphalistinv[i]=(long)element_inv(nf,p2);
2553: }
2554: }
2555: for (j=1; j<=n; j++)
2556: {
2557: p1=cgetg(n+1,t_COL); A1[j]=(long)p1;
2558: for (i=1; i<=n; i++)
2559: p1[i] = (long)element_mul(nf,gcoeff(A,i,j),(GEN)alphalist[j]);
2560: }
2561: Aa=basistoalg(nf,A1); Aainv=lift_intern(ginv(Aa));
2562: Aaa=mat_to_vecpol(Aa,vpol);
2563: for (i=1; i<=n; i++) for (j=i; j<=n; j++)
2564: {
2565: long tp;
2566: p1 = lift_intern(gres(gmul((GEN)Aaa[i],(GEN)Aaa[j]), pol));
2567: tp = typ(p1);
2568: if (is_scalar_t(tp) || (tp==t_POL && varn(p1)>vpol))
2569: p2 = gmul(p1, (GEN)Aainv[1]);
2570: else
2571: p2 = gmul(Aainv, pol_to_vec(p1, n));
2572:
2573: p3 = algtobasis(nf,p2);
2574: for (k=1; k<=n; k++)
2575: {
2576: coeff(multab,k,(i-1)*n+j) = p3[k];
2577: coeff(multab,k,(j-1)*n+i) = p3[k];
2578: }
2579: }
2580: R=cgetg(n+1,t_MAT); multabmod = mymod(multab,p);
2581: R[1] = matId[1];
2582: for (j=2; j<=n; j++)
2583: R[j] = (long) rnfelementid_powmod(nf,multabmod,matId, j,q1,prhall);
2584: baseIp = nfkermodpr(nf,R,prhall);
2585: baseOp = lift_intern(nfsuppl(nf,baseIp,n,prhall));
2586: alpha=cgetg(n+1,t_MAT);
2587: for (j=1; j<lg(baseIp); j++) alpha[j]=baseOp[j];
2588: for ( ; j<=n; j++)
2589: {
2590: p1=cgetg(n+1,t_COL); alpha[j]=(long)p1;
2591: for (i=1; i<=n; i++)
2592: p1[i]=(long)element_mul(nf,pip,gcoeff(baseOp,i,j));
2593: }
2594: matprod=cgetg(n+1,t_MAT);
2595: for (j=1; j<=n; j++)
2596: {
2597: p1=cgetg(n+1,t_COL); matprod[j]=(long)p1;
2598: for (i=1; i<=n; i++)
2599: {
2600: p2 = rnfelement_mulidmod(nf,multab,unnf, (GEN)alpha[i],j, NULL);
2601: for (k=1; k<=n; k++)
2602: p2[k] = lmul((GEN)nf[7], (GEN)p2[k]);
2603: p1[i] = (long)p2;
2604: }
2605: }
2606: alphainv = lift_intern(ginv(basistoalg(nf,alpha)));
2607: matC = cgetg(n+1,t_MAT);
2608: for (j=1; j<=n; j++)
2609: {
2610: p1=cgetg(n*n+1,t_COL); matC[j]=(long)p1;
2611: for (i=1; i<=n; i++)
2612: {
2613: p2 = gmul(alphainv, gcoeff(matprod,i,j));
2614: for (k=1; k<=n; k++)
2615: p1[(i-1)*n+k]=(long)nfreducemodpr(nf,algtobasis(nf,(GEN)p2[k]),prhall);
2616: }
2617: }
2618: matG=nfkermodpr(nf,matC,prhall); m=lg(matG)-1;
2619: vecpro=cgetg(3,t_VEC);
2620: p1=cgetg(n+m+1,t_MAT); vecpro[1]=(long)p1;
2621: p2=cgetg(n+m+1,t_VEC); vecpro[2]=(long)p2;
2622: for (j=1; j<=m; j++)
2623: {
2624: p1[j] = llift((GEN)matG[j]);
2625: p2[j] = (long)prhinv;
2626: }
2627: p1 += m;
2628: p2 += m;
2629: for (j=1; j<=n; j++)
2630: {
2631: p1[j] = matId[j];
2632: p2[j] = (long)idealmul(nf,(GEN)I[j],(GEN)alphalistinv[j]);
2633: }
2634: matH=nfhermite(nf,vecpro);
2635: p1=algtobasis(nf,gmul(Aa,basistoalg(nf,(GEN)matH[1])));
2636: p2=(GEN)matH[2];
2637:
2638: tetpil=avma; neworder=cgetg(3,t_VEC);
2639: H=cgetg(n+1,t_MAT); Hid=cgetg(n+1,t_VEC);
2640: for (j=1; j<=n; j++)
2641: {
2642: p3=cgetg(n+1,t_COL); H[j]=(long)p3;
2643: for (i=1; i<=n; i++)
2644: p3[i]=(long)element_mul(nf,gcoeff(p1,i,j),(GEN)alphalistinv[j]);
2645: Hid[j]=(long)idealmul(nf,(GEN)p2[j],(GEN)alphalist[j]);
2646: }
2647: if (DEBUGLEVEL>3)
2648: { fprintferr(" new order:\n"); outerr(H); outerr(Hid); }
2649: if (sep == 2 || gegal(I,Hid))
2650: {
2651: neworder[1]=(long)H; neworder[2]=(long)Hid;
2652: return gerepile(av,tetpil,neworder);
2653: }
2654:
2655: A=H; I=Hid;
2656: if (low_stack(lim, stack_lim(av1,1)) || (cmpt & 3) == 0)
2657: {
2658: GEN *gptr[2]; gptr[0]=&A; gptr[1]=&I;
2659: if(DEBUGMEM>1) err(warnmem,"rnfordmax");
2660: gerepilemany(av1,gptr,2);
2661: }
2662: }
2663: }
2664:
2665: static void
2666: check_pol(GEN x, long v)
2667: {
2668: long i,tx, lx = lg(x);
2669: if (varn(x) != v)
2670: err(talker,"incorrect variable in rnf function");
2671: for (i=2; i<lx; i++)
2672: {
2673: tx = typ(x[i]);
2674: if (!is_scalar_t(tx) || tx == t_POLMOD)
2675: err(talker,"incorrect polcoeff in rnf function");
2676: }
2677: }
2678:
2679: GEN
2680: fix_relative_pol(GEN nf, GEN x, int chk_lead)
2681: {
2682: GEN xnf = (typ(nf) == t_POL)? nf: (GEN)nf[1];
2683: long i, vnf = varn(xnf), lx = lg(x);
2684: if (typ(x) != t_POL || varn(x) >= vnf)
2685: err(talker,"incorrect polynomial in rnf function");
2686: x = dummycopy(x);
2687: for (i=2; i<lx; i++)
2688: switch(typ(x[i]))
2689: {
2690: case t_POL:
2691: check_pol((GEN)x[i], vnf);
2692: x[i] = lmodulcp((GEN)x[i], xnf); break;
2693: case t_POLMOD:
2694: if (!gegal(gmael(x,i,1), xnf)) err(consister,"rnf function");
2695: break;
2696: }
2697:
2698: if (chk_lead && !gcmp1(leading_term(x)))
2699: err(impl,"non-monic relative polynomials");
2700: return x;
2701: }
2702:
2703: static GEN
2704: rnfround2all(GEN nf, GEN pol, long all)
2705: {
2706: long av=avma,tetpil,i,j,n,N,nbidp,vpol,*ep;
2707: GEN p1,p2,p3,p4,polnf,list,unnf,id,matId,I,W,pseudo,y,discpol,d,D,sym;
2708:
2709: nf=checknf(nf); polnf=(GEN)nf[1]; vpol=varn(pol);
2710: pol = fix_relative_pol(nf,pol,1);
2711: N=degpol(polnf); n=degpol(pol); discpol=discsr(pol);
2712: list=idealfactor(nf,discpol); ep=(long*)list[2]; list=(GEN)list[1];
2713: nbidp=lg(list)-1; for(i=1;i<=nbidp;i++) ep[i]=itos((GEN)ep[i]);
2714: if (DEBUGLEVEL>1)
2715: {
2716: fprintferr("Ideals to consider:\n");
2717: for (i=1; i<=nbidp; i++)
2718: if (ep[i]>1) fprintferr("%Z^%ld\n",list[i],ep[i]);
2719: flusherr();
2720: }
2721: id=idmat(N); unnf=gscalcol_i(gun,N);
2722: matId=idmat_intern(n,unnf, gscalcol_i(gzero,N));
2723: pseudo = NULL;
2724: for (i=1; i<=nbidp; i++)
2725: if (ep[i] > 1)
2726: {
2727: y=rnfordmax(nf,pol,(GEN)list[i],unnf,id,matId);
2728: pseudo = rnfjoinmodules(nf,pseudo,y);
2729: }
2730: if (!pseudo)
2731: {
2732: I=cgetg(n+1,t_VEC); for (i=1; i<=n; i++) I[i]=(long)id;
2733: pseudo=cgetg(3,t_VEC); pseudo[1]=(long)matId; pseudo[2]=(long)I;
2734: }
2735: W=gmodulcp(mat_to_vecpol(basistoalg(nf,(GEN)pseudo[1]),vpol),pol);
2736: p2=cgetg(n+1,t_MAT); for (j=1; j<=n; j++) p2[j]=lgetg(n+1,t_COL);
2737: sym=polsym(pol,n-1);
2738: for (j=1; j<=n; j++)
2739: for (i=j; i<=n; i++)
2740: {
2741: p1 = lift_intern(gmul((GEN)W[i],(GEN)W[j]));
2742: coeff(p2,j,i)=coeff(p2,i,j)=(long)quicktrace(p1,sym);
2743: }
2744: d = algtobasis_intern(nf,det(p2));
2745:
2746: I=(GEN)pseudo[2];
2747: i=1; while (i<=n && gegal((GEN)I[i],id)) i++;
2748: if (i>n) D=id;
2749: else
2750: {
2751: D = (GEN)I[i];
2752: for (i++; i<=n; i++)
2753: if (!gegal((GEN)I[i],id)) D = idealmul(nf,D,(GEN)I[i]);
2754: D = idealpow(nf,D,gdeux);
2755: }
2756: p4=gun; p3=auxdecomp(content(d),0);
2757: for (i=1; i<lg(p3[1]); i++)
2758: p4 = gmul(p4, gpuigs(gcoeff(p3,i,1), itos(gcoeff(p3,i,2))>>1));
2759: p4 = gsqr(p4); tetpil=avma;
2760: i = all? 2: 0;
2761: p1=cgetg(3 + i,t_VEC);
2762: if (i) { p1[1]=lcopy((GEN)pseudo[1]); p1[2]=lcopy(I); }
2763: p1[1+i] = (long)idealmul(nf,D,d);
2764: p1[2+i] = ldiv(d,p4);
2765: return gerepile(av,tetpil,p1);
2766: }
2767:
2768: GEN
2769: rnfpseudobasis(GEN nf, GEN pol)
2770: {
2771: return rnfround2all(nf,pol,1);
2772: }
2773:
2774: GEN
2775: rnfdiscf(GEN nf, GEN pol)
2776: {
2777: return rnfround2all(nf,pol,0);
2778: }
2779:
2780: /* given bnf as output by buchinit and a pseudo-basis of an order
2781: * in HNF [A,I] (or [A,I,D,d] it does not matter), tries to simplify the
2782: * HNF as much as possible. The resulting matrix will be upper triangular
2783: * but the diagonal coefficients will not be equal to 1. The ideals
2784: * are guaranteed to be integral and primitive.
2785: */
2786: GEN
2787: rnfsimplifybasis(GEN bnf, GEN order)
2788: {
2789: long av=avma,tetpil,j,N,n;
2790: GEN p1,id,Az,Iz,nf,A,I;
2791:
2792: bnf = checkbnf(bnf);
2793: if (typ(order)!=t_VEC || lg(order)<3)
2794: err(talker,"not a pseudo-basis in nfsimplifybasis");
2795: A=(GEN)order[1]; I=(GEN)order[2]; n=lg(A)-1; nf=(GEN)bnf[7];
2796: N=degpol(nf[1]); id=idmat(N); Iz=cgetg(n+1,t_VEC); Az=cgetg(n+1,t_MAT);
2797: for (j=1; j<=n; j++)
2798: {
2799: if (gegal((GEN)I[j],id)) { Iz[j]=(long)id; Az[j]=A[j]; }
2800: else
2801: {
2802: p1=content((GEN)I[j]);
2803: if (!gcmp1(p1))
2804: {
2805: Iz[j]=(long)gdiv((GEN)I[j],p1); Az[j]=lmul((GEN)A[j],p1);
2806: }
2807: else Az[j]=A[j];
2808: if (!gegal((GEN)Iz[j],id))
2809: {
2810: p1=isprincipalgen(bnf,(GEN)Iz[j]);
2811: if (gcmp0((GEN)p1[1]))
2812: {
2813: p1=(GEN)p1[2]; Iz[j]=(long)id;
2814: Az[j]=(long)element_mulvec(nf,p1,(GEN)Az[j]);
2815: }
2816: }
2817: }
2818: }
2819: tetpil=avma; p1=cgetg(lg(order),t_VEC); p1[1]=lcopy(Az); p1[2]=lcopy(Iz);
2820: for (j=3; j<lg(order); j++) p1[j]=lcopy((GEN)order[j]);
2821: return gerepile(av,tetpil,p1);
2822: }
2823:
2824: GEN
2825: rnfdet2(GEN nf, GEN A, GEN I)
2826: {
2827: long av,tetpil,i;
2828: GEN p1;
2829:
2830: nf=checknf(nf); av = tetpil = avma;
2831: p1=idealhermite(nf,det(matbasistoalg(nf,A)));
2832: for(i=1;i<lg(I);i++) { tetpil=avma; p1=idealmul(nf,p1,(GEN)I[i]); }
2833: tetpil=avma; return gerepile(av,tetpil,p1);
2834: }
2835:
2836: GEN
2837: rnfdet(GEN nf, GEN order)
2838: {
2839: if (typ(order)!=t_VEC || lg(order)<3)
2840: err(talker,"not a pseudo-matrix in rnfdet");
2841: return rnfdet2(nf,(GEN)order[1],(GEN)order[2]);
2842: }
2843:
2844: GEN
2845: rnfdet0(GEN nf, GEN x, GEN y)
2846: {
2847: return y? rnfdet2(nf,x,y): rnfdet(nf,x);
2848: }
2849:
2850: /* given a pseudo-basis of an order in HNF [A,I] (or [A,I,D,d] it does
2851: * not matter), gives an nxn matrix (not in HNF) of a pseudo-basis and
2852: * an ideal vector [id,id,...,id,I] such that order=nf[7]^(n-1)xI.
2853: * Since it uses the approximation theorem, can be long.
2854: */
2855: GEN
2856: rnfsteinitz(GEN nf, GEN order)
2857: {
2858: long av=avma,tetpil,i,n;
2859: GEN Id,A,I,p1,a,b;
2860:
2861: nf = checknf(nf);
2862: Id = idmat(degpol(nf[1]));
2863: if (typ(order)==t_POL) order=rnfpseudobasis(nf,order);
2864: if (typ(order)!=t_VEC || lg(order)<3)
2865: err(talker,"not a pseudo-matrix in rnfsteinitz");
2866: A=dummycopy((GEN)order[1]);
2867: I=dummycopy((GEN)order[2]); n=lg(A)-1;
2868: if (typ(A) != t_MAT || typ(I) != t_VEC || lg(I) != n+1)
2869: err(typeer,"rnfsteinitz");
2870: for (i=1; i<n; i++)
2871: {
2872: a = (GEN)I[i];
2873: if (!gegal(a,Id))
2874: {
2875: GEN c1 = (GEN)A[i];
2876: GEN c2 = (GEN)A[i+1];
2877: b = (GEN)I[i+1];
2878: if (gegal(b,Id))
2879: {
2880: A[i] = (long)c2;
2881: A[i+1]= lneg(c1);
2882: I[i] = (long)b;
2883: I[i+1]= (long)a;
2884: }
2885: else
2886: {
2887: p1 = nfidealdet1(nf,a,b);
2888: A[i] = ladd(element_mulvec(nf,(GEN)p1[1], c1),
2889: element_mulvec(nf,(GEN)p1[2], c2));
2890: A[i+1]= ladd(element_mulvec(nf,(GEN)p1[3], c1),
2891: element_mulvec(nf,(GEN)p1[4], c2));
2892: I[i] =(long)Id;
2893: I[i+1]=(long)idealmul(nf,a,b); p1 = content((GEN)I[i+1]);
2894: if (!gcmp1(p1))
2895: {
2896: I[i+1] = ldiv((GEN)I[i+1],p1);
2897: A[i+1] = lmul(p1,(GEN)A[i+1]);
2898: }
2899: }
2900: }
2901: }
2902: tetpil=avma; p1=cgetg(lg(order),t_VEC);
2903: p1[1]=lcopy(A); p1[2]=lcopy(I);
2904: for (i=3; i<lg(order); i++) p1[i]=lcopy((GEN)order[i]);
2905: return gerepile(av,tetpil,p1);
2906: }
2907:
2908: /* Given bnf as output by buchinit and either an order as output by
2909: * rnfpseudobasis or a polynomial, and outputs a basis if it is free,
2910: * an n+1-generating set if it is not
2911: */
2912: GEN
2913: rnfbasis(GEN bnf, GEN order)
2914: {
2915: ulong av = avma;
2916: long j,N,n;
2917: GEN nf,A,I,classe,p1,p2,id;
2918:
2919: bnf = checkbnf(bnf);
2920: nf=(GEN)bnf[7]; N=degpol(nf[1]); id=idmat(N);
2921: if (typ(order)==t_POL) order=rnfpseudobasis(nf,order);
2922: if (typ(order)!=t_VEC || lg(order)<3)
2923: err(talker,"not a pseudo-matrix in rnfbasis");
2924: A=(GEN)order[1]; I=(GEN)order[2]; n=lg(A)-1;
2925: j=1; while (j<n && gegal((GEN)I[j],id)) j++;
2926: if (j<n) order=rnfsteinitz(nf,order);
2927: A=(GEN)order[1]; I=(GEN)order[2]; classe=(GEN)I[n];
2928: p1=isprincipalgen(bnf,classe);
2929: if (gcmp0((GEN)p1[1]))
2930: {
2931: p2=cgetg(n+1,t_MAT);
2932: p2[n]=(long)element_mulvec(nf,(GEN)p1[2],(GEN)A[n]);
2933: }
2934: else
2935: {
2936: p1=ideal_two_elt(nf,classe);
2937: p2=cgetg(n+2,t_MAT);
2938: p2[n]=lmul((GEN)p1[1],(GEN)A[n]);
2939: p2[n+1]=(long)element_mulvec(nf,(GEN)p1[2],(GEN)A[n]);
2940: }
2941: for (j=1; j<n; j++) p2[j]=A[j];
2942: return gerepilecopy(av,p2);
2943: }
2944:
2945: /* Given bnf as output by buchinit and either an order as output by
2946: * rnfpseudobasis or a polynomial, and outputs a basis (not pseudo)
2947: * in Hermite Normal Form if it exists, zero if not
2948: */
2949: GEN
2950: rnfhermitebasis(GEN bnf, GEN order)
2951: {
2952: ulong av = avma;
2953: long j,N,n;
2954: GEN nf,A,I,p1,id;
2955:
2956: bnf = checkbnf(bnf); nf=(GEN)bnf[7];
2957: N=degpol(nf[1]); id=idmat(N);
2958: if (typ(order)==t_POL)
2959: {
2960: order=rnfpseudobasis(nf,order);
2961: A=(GEN)order[1];
2962: }
2963: else
2964: {
2965: if (typ(order)!=t_VEC || lg(order)<3)
2966: err(talker,"not a pseudo-matrix in rnfbasis");
2967: A=gcopy((GEN)order[1]);
2968: }
2969: I=(GEN)order[2]; n=lg(A)-1;
2970: for (j=1; j<=n; j++)
2971: {
2972: if (!gegal((GEN)I[j],id))
2973: {
2974: p1=isprincipalgen(bnf,(GEN)I[j]);
2975: if (gcmp0((GEN)p1[1]))
2976: A[j]=(long)element_mulvec(nf,(GEN)p1[2],(GEN)A[j]);
2977: else { avma=av; return gzero; }
2978: }
2979: }
2980: return gerepilecopy(av,A);
2981: }
2982:
2983: long
2984: rnfisfree(GEN bnf, GEN order)
2985: {
2986: long av=avma,n,N,j;
2987: GEN nf,p1,id,I;
2988:
2989: bnf = checkbnf(bnf);
2990: if (gcmp1(gmael3(bnf,8,1,1))) return 1;
2991:
2992: nf=(GEN)bnf[7]; N=degpol(nf[1]); id=idmat(N);
2993: if (typ(order)==t_POL) order=rnfpseudobasis(nf,order);
2994: if (typ(order)!=t_VEC || lg(order)<3)
2995: err(talker,"not a pseudo-matrix in rnfisfree");
2996:
2997: I=(GEN)order[2]; n=lg(I)-1;
2998: j=1; while (j<=n && gegal((GEN)I[j],id)) j++;
2999: if (j>n) { avma=av; return 1; }
3000:
3001: p1=(GEN)I[j];
3002: for (j++; j<=n; j++)
3003: if (!gegal((GEN)I[j],id)) p1=idealmul(nf,p1,(GEN)I[j]);
3004: j = gcmp0(isprincipal(bnf,p1));
3005: avma=av; return j;
3006: }
3007:
3008: /**********************************************************************/
3009: /** **/
3010: /** COMPOSITUM OF TWO NUMBER FIELDS **/
3011: /** **/
3012: /**********************************************************************/
3013: extern GEN ZY_ZXY_resultant_all(GEN A, GEN B0, long *lambda, GEN *LPRS);
3014: extern GEN squff2(GEN x, long klim, long hint);
3015: extern GEN to_polmod(GEN x, GEN mod);
3016:
3017: /* modular version. TODO: check that compositum2 is not slower */
3018: GEN
3019: polcompositum0(GEN A, GEN B, long flall)
3020: {
3021: ulong av = avma;
3022: long v,k;
3023: GEN C, LPRS;
3024:
3025: if (typ(A)!=t_POL || typ(B)!=t_POL) err(typeer,"polcompositum0");
3026: if (degpol(A)<=0 || degpol(B)<=0) err(constpoler,"compositum");
3027: v = varn(A);
3028: if (varn(B) != v) err(talker,"not the same variable in compositum");
3029: C = content(A); if (!gcmp1(C)) A = gdiv(A, C);
3030: C = content(B); if (!gcmp1(C)) B = gdiv(B, C);
3031: check_pol_int(A,"compositum");
3032: check_pol_int(B,"compositum");
3033: if (!ZX_is_squarefree(A)) err(talker,"compositum: %Z not separable", A);
3034: if (!ZX_is_squarefree(B)) err(talker,"compositum: %Z not separable", B);
3035:
3036: k = 1; C = ZY_ZXY_resultant_all(A, B, &k, flall? &LPRS: NULL);
3037: C = squff2(C,0,0); /* C = Res_Y (A, B(X + kY)) guaranteed squarefree */
3038: if (flall)
3039: {
3040: long i,l = lg(C);
3041: GEN w,a,b; /* a,b,c root of A,B,C = compositum, c = b - k a */
3042: for (i=1; i<l; i++)
3043: { /* invmod possibly very costly */
3044: a = gmul((GEN)LPRS[1], ZX_invmod((GEN)LPRS[2], (GEN)C[i]));
3045: a = gneg_i(gmod(a, (GEN)C[i]));
3046: b = gadd(polx[v], gmulsg(k,a));
3047: w = cgetg(5,t_VEC); /* [C, a, b, n ] */
3048: w[1] = C[i];
3049: w[2] = (long)to_polmod(a, (GEN)w[1]);
3050: w[3] = (long)to_polmod(b, (GEN)w[1]);
3051: w[4] = lstoi(-k); C[i] = (long)w;
3052: }
3053: }
3054: settyp(C, t_VEC); return gerepilecopy(av, C);
3055: }
3056:
3057: GEN
3058: compositum(GEN pol1,GEN pol2)
3059: {
3060: return polcompositum0(pol1,pol2,0);
3061: }
3062:
3063: GEN
3064: compositum2(GEN pol1,GEN pol2)
3065: {
3066: return polcompositum0(pol1,pol2,1);
3067: }
3068:
3069: extern int isrational(GEN x);
3070: extern GEN nfgcd(GEN P, GEN Q, GEN nf, GEN den);
3071:
3072: int
3073: nfissquarefree(GEN nf, GEN x)
3074: {
3075: ulong av = avma;
3076: GEN g, y = derivpol(x);
3077: if (isrational(x))
3078: g = modulargcd(x, y);
3079: else
3080: g = nfgcd(x, y, nf, NULL);
3081: avma = av; return (degpol(g) == 0);
3082: }
3083:
3084: GEN
3085: rnfequation0(GEN nf, GEN B, long flall)
3086: {
3087: ulong av = avma;
3088: long v,vpol,k,lA,lB;
3089: GEN cC,A,C,LPRS;
3090:
3091: if (typ(nf)==t_POL) A=nf; else { nf=checknf(nf); A=(GEN)nf[1]; }
3092: B = fix_relative_pol(nf,B,1);
3093: v = varn(A); lA = lgef(A);
3094: vpol= varn(B); lB = lgef(B);
3095: if (lA<=3 || lB<=3) err(constpoler,"rnfequation");
3096:
3097: check_pol_int(A,"rnfequation");
3098: B = lift_intern(B); B = gdiv(B, content(B));
3099: for (k=2; k<lB; k++)
3100: if (lgef(B[k]) >= lA) B[k] = lres((GEN)B[k],A);
3101:
3102: if (!nfissquarefree(A,B))
3103: err(talker,"not k separable relative equation in rnfequation");
3104:
3105: k = 0; C = ZY_ZXY_resultant_all(A, B, &k, flall? &LPRS: NULL);
3106: if (gsigne(leadingcoeff(C)) < 0) C = gneg_i(C);
3107: C = primitive_part(C, &cC);
3108: if (flall)
3109: {
3110: GEN w,a,b; /* a,b,c root of A,B,C = compositum, c = b - k a */
3111: /* invmod possibly very costly */
3112: a = gmul((GEN)LPRS[1], ZX_invmod((GEN)LPRS[2], C));
3113: a = gneg_i(gmod(a, C));
3114: b = gadd(polx[v], gmulsg(k,a));
3115: w = cgetg(4,t_VEC); /* [C, a, n ] */
3116: w[1] = (long)C;
3117: w[2] = (long)to_polmod(a, (GEN)w[1]);
3118: w[3] = lstoi(-k); C = w;
3119: }
3120: return gerepilecopy(av, C);
3121: }
3122:
3123: GEN
3124: rnfequation(GEN nf,GEN pol2)
3125: {
3126: return rnfequation0(nf,pol2,0);
3127: }
3128:
3129: GEN
3130: rnfequation2(GEN nf,GEN pol2)
3131: {
3132: return rnfequation0(nf,pol2,1);
3133: }
3134:
3135: static GEN
3136: nftau(long r1, GEN x)
3137: {
3138: long i, ru = lg(x);
3139: GEN s;
3140:
3141: s = r1 ? (GEN)x[1] : gmul2n(greal((GEN)x[1]),1);
3142: for (i=2; i<=r1; i++) s=gadd(s,(GEN)x[i]);
3143: for ( ; i<ru; i++) s=gadd(s,gmul2n(greal((GEN)x[i]),1));
3144: return s;
3145: }
3146:
3147: static GEN
3148: nftocomplex(GEN nf, GEN x)
3149: {
3150: long ru,vnf,k;
3151: GEN p2,p3,ronf;
3152:
3153: p2 = (typ(x)==t_POLMOD)? (GEN)x[2]: gmul((GEN)nf[7],x);
3154: vnf=varn(nf[1]);
3155: ronf=(GEN)nf[6]; ru=lg(ronf); p3=cgetg(ru,t_COL);
3156: for (k=1; k<ru; k++) p3[k]=lsubst(p2,vnf,(GEN)ronf[k]);
3157: return p3;
3158: }
3159:
3160: static GEN
3161: rnfscal(GEN mth, GEN xth, GEN yth)
3162: {
3163: long n,ru,i,j,kk;
3164: GEN x,y,m,res,p1,p2;
3165:
3166: n=lg(mth)-1; ru=lg(gcoeff(mth,1,1));
3167: res=cgetg(ru,t_COL);
3168: for (kk=1; kk<ru; kk++)
3169: {
3170: m=cgetg(n+1,t_MAT);
3171: for (j=1; j<=n; j++)
3172: {
3173: p1=cgetg(n+1,t_COL); m[j]=(long)p1;
3174: for (i=1; i<=n; i++) { p2=gcoeff(mth,i,j); p1[i]=p2[kk]; }
3175: }
3176: x=cgetg(n+1,t_VEC);
3177: for (j=1; j<=n; j++) x[j]=(long)gconj((GEN)((GEN)xth[j])[kk]);
3178: y=cgetg(n+1,t_COL);
3179: for (j=1; j<=n; j++) y[j]=((GEN)yth[j])[kk];
3180: res[kk]=(long)gmul(x,gmul(m,y));
3181: }
3182: return res;
3183: }
3184:
3185: static GEN
3186: rnfdiv(GEN x, GEN y)
3187: {
3188: long i, ru = lg(x);
3189: GEN z;
3190:
3191: z=cgetg(ru,t_COL);
3192: for (i=1; i<ru; i++) z[i]=(long)gdiv((GEN)x[i],(GEN)y[i]);
3193: return z;
3194: }
3195:
3196: static GEN
3197: rnfmul(GEN x, GEN y)
3198: {
3199: long i, ru = lg(x);
3200: GEN z;
3201:
3202: z=cgetg(ru,t_COL);
3203: for (i=1; i<ru; i++) z[i]=(long)gmul((GEN)x[i],(GEN)y[i]);
3204: return z;
3205: }
3206:
3207: static GEN
3208: rnfvecmul(GEN x, GEN v)
3209: {
3210: long i, lx = lg(v);
3211: GEN y;
3212:
3213: y=cgetg(lx,typ(v));
3214: for (i=1; i<lx; i++) y[i]=(long)rnfmul(x,(GEN)v[i]);
3215: return y;
3216: }
3217:
3218: static GEN
3219: allonge(GEN v, long N)
3220: {
3221: long r,r2,i;
3222: GEN y;
3223:
3224: r=lg(v)-1; r2=N-r;
3225: y=cgetg(N+1,t_COL);
3226: for (i=1; i<=r; i++) y[i]=v[i];
3227: for ( ; i<=N; i++) y[i]=(long)gconj((GEN)v[i-r2]);
3228: return y;
3229: }
3230:
3231: static GEN
3232: findmin(GEN nf, GEN ideal, GEN muf,long prec)
3233: {
3234: long av=avma,N,tetpil,i;
3235: GEN m,y;
3236:
3237: m = qf_base_change(gmael(nf,5,3), ideal, 0); /* nf[5][3] = T2 */
3238: m = lllgramintern(m,4,1,prec);
3239: if (!m)
3240: {
3241: m = lllint(ideal);
3242: m = qf_base_change(gmael(nf,5,3), gmul(ideal,m), 0);
3243: m = lllgramintern(m,4,1,prec);
3244: if (!m) err(precer,"rnflllgram");
3245: }
3246: ideal=gmul(ideal,m);
3247: N=lg(ideal)-1; y=cgetg(N+1,t_MAT);
3248: for (i=1; i<=N; i++)
3249: y[i] = (long) allonge(nftocomplex(nf,(GEN)ideal[i]),N);
3250: m=ground(greal(gauss(y,allonge(muf,N))));
3251: tetpil=avma; return gerepile(av,tetpil,gmul(ideal,m));
3252: }
3253:
3254: #define swap(x,y) { long _t=x; x=y; y=_t; }
3255:
3256: /* given a base field nf (e.g main variable y), a polynomial pol with
3257: * coefficients in nf (e.g main variable x), and an order as output
3258: * by rnfpseudobasis, outputs a reduced order.
3259: */
3260: GEN
3261: rnflllgram(GEN nf, GEN pol, GEN order,long prec)
3262: {
3263: long av=avma,tetpil,i,j,k,l,kk,kmax,r1,ru,lx,vnf;
3264: GEN p1,p2,M,I,U,ronf,poll,unro,roorder,powreorder,mth,s,MC,MPOL,MCS;
3265: GEN B,mu,Bf,temp,ideal,x,xc,xpol,muf,mufc,muno,y,z,Ikk_inv;
3266:
3267: /* Initializations and verifications */
3268:
3269: nf=checknf(nf);
3270: if (typ(order)!=t_VEC || lg(order)<3)
3271: err(talker,"not a pseudo-matrix in rnflllgram");
3272: M=(GEN)order[1]; I=(GEN)order[2]; lx=lg(I);
3273: if (lx < 3) return gcopy(order);
3274: if (lx-1 != degpol(pol)) err(consister,"rnflllgram");
3275: U=idmat(lx-1); I = dummycopy(I);
3276:
3277: /* Compute the relative T2 matrix of powers of theta */
3278:
3279: vnf=varn(nf[1]); ronf=(GEN)nf[6]; ru=lg(ronf); poll=lift(pol);
3280: r1 = nf_get_r1(nf);
3281: unro=cgetg(lx,t_COL); for (i=1; i<lx; i++) unro[i]=un;
3282: roorder=cgetg(ru,t_VEC);
3283: for (i=1; i<ru; i++)
3284: roorder[i]=lroots(gsubst(poll,vnf,(GEN)ronf[i]),prec);
3285: powreorder=cgetg(lx,t_MAT);
3286: p1=cgetg(ru,t_COL); powreorder[1]=(long)p1;
3287: for (i=1; i<ru; i++) p1[i]=(long)unro;
3288: for (k=2; k<lx; k++)
3289: {
3290: p1=cgetg(ru,t_COL); powreorder[k]=(long)p1;
3291: for (i=1; i<ru; i++)
3292: {
3293: p2=cgetg(lx,t_COL); p1[i]=(long)p2;
3294: for (j=1; j<lx; j++)
3295: p2[j] = lmul(gmael(roorder,i,j),gmael3(powreorder,k-1,i,j));
3296: }
3297: }
3298: mth=cgetg(lx,t_MAT);
3299: for (l=1; l<lx; l++)
3300: {
3301: p1=cgetg(lx,t_COL); mth[l]=(long)p1;
3302: for (k=1; k<lx; k++)
3303: {
3304: p2=cgetg(ru,t_COL); p1[k]=(long)p2;
3305: for (i=1; i<ru; i++)
3306: {
3307: s=gzero;
3308: for (j=1; j<lx; j++)
3309: s = gadd(s,gmul(gconj(gmael3(powreorder,k,i,j)),
3310: gmael3(powreorder,l,i,j)));
3311: p2[i]=(long)s;
3312: }
3313: }
3314: }
3315:
3316: /* Transform the matrix M into a matrix with coefficients in K and also
3317: with coefficients polymod */
3318:
3319: MC=cgetg(lx,t_MAT); MPOL=cgetg(lx,t_MAT);
3320: for (j=1; j<lx; j++)
3321: {
3322: p1=cgetg(lx,t_COL); MC[j]=(long)p1;
3323: p2=cgetg(lx,t_COL); MPOL[j]=(long)p2;
3324: for (i=1; i<lx; i++)
3325: {
3326: p2[i]=(long)basistoalg(nf,gcoeff(M,i,j));
3327: p1[i]=(long)nftocomplex(nf,(GEN)p2[i]);
3328: }
3329: }
3330: MCS=cgetg(lx,t_MAT);
3331:
3332: /* Start LLL algorithm */
3333:
3334: mu=cgetg(lx,t_MAT); B=cgetg(lx,t_COL);
3335: for (j=1; j<lx; j++)
3336: {
3337: p1=cgetg(lx,t_COL); mu[j]=(long)p1; for (i=1; i<lx; i++) p1[i]=zero;
3338: B[j]=zero;
3339: }
3340: kk=2; if (DEBUGLEVEL) fprintferr("kk = %ld ",kk);
3341: kmax=1; B[1]=lreal(rnfscal(mth,(GEN)MC[1],(GEN)MC[1]));
3342: MCS[1]=lcopy((GEN)MC[1]);
3343: do
3344: {
3345: if (kk>kmax)
3346: {
3347: /* Incremental Gram-Schmidt */
3348: kmax=kk; MCS[kk]=lcopy((GEN)MC[kk]);
3349: for (j=1; j<kk; j++)
3350: {
3351: coeff(mu,kk,j) = (long) rnfdiv(rnfscal(mth,(GEN)MCS[j],(GEN)MC[kk]),
3352: (GEN) B[j]);
3353: MCS[kk] = lsub((GEN) MCS[kk], rnfvecmul(gcoeff(mu,kk,j),(GEN)MCS[j]));
3354: }
3355: B[kk] = lreal(rnfscal(mth,(GEN)MCS[kk],(GEN)MCS[kk]));
3356: if (gcmp0((GEN)B[kk])) err(lllger3);
3357: }
3358:
3359: /* RED(k,k-1) */
3360: l=kk-1; Ikk_inv=idealinv(nf, (GEN)I[kk]);
3361: ideal=idealmul(nf,(GEN)I[l],Ikk_inv);
3362: x=findmin(nf,ideal,gcoeff(mu,kk,l),2*prec-2);
3363: if (!gcmp0(x))
3364: {
3365: xpol=basistoalg(nf,x); xc=nftocomplex(nf,xpol);
3366: MC[kk]=lsub((GEN)MC[kk],rnfvecmul(xc,(GEN)MC[l]));
3367: U[kk]=lsub((GEN)U[kk],gmul(xpol,(GEN)U[l]));
3368: coeff(mu,kk,l)=lsub(gcoeff(mu,kk,l),xc);
3369: for (i=1; i<l; i++)
3370: coeff(mu,kk,i)=lsub(gcoeff(mu,kk,i),rnfmul(xc,gcoeff(mu,l,i)));
3371: }
3372: /* Test LLL condition */
3373: p1=nftau(r1,gadd((GEN) B[kk],
3374: gmul(gnorml2(gcoeff(mu,kk,kk-1)),(GEN)B[kk-1])));
3375: p2=gdivgs(gmulsg(9,nftau(r1,(GEN)B[kk-1])),10);
3376: if (gcmp(p1,p2)<=0)
3377: {
3378: /* Execute SWAP(k) */
3379: k=kk;
3380: swap(MC[k-1],MC[k]);
3381: swap(U[k-1],U[k]);
3382: swap(I[k-1],I[k]);
3383: for (j=1; j<=k-2; j++) swap(coeff(mu,k-1,j),coeff(mu,k,j));
3384: muf=gcoeff(mu,k,k-1);
3385: mufc=gconj(muf); muno=greal(rnfmul(muf,mufc));
3386: Bf=gadd((GEN)B[k],rnfmul(muno,(GEN)B[k-1]));
3387: p1=rnfdiv((GEN)B[k-1],Bf);
3388: coeff(mu,k,k-1)=(long)rnfmul(mufc,p1);
3389: temp=(GEN)MCS[k-1];
3390: MCS[k-1]=ladd((GEN)MCS[k],rnfvecmul(muf,(GEN)MCS[k-1]));
3391: MCS[k]=lsub(rnfvecmul(rnfdiv((GEN)B[k],Bf),temp),
3392: rnfvecmul(gcoeff(mu,k,k-1),(GEN)MCS[k]));
3393: B[k]=(long)rnfmul((GEN)B[k],p1); B[k-1]=(long)Bf;
3394: for (i=k+1; i<=kmax; i++)
3395: {
3396: temp=gcoeff(mu,i,k);
3397: coeff(mu,i,k)=lsub(gcoeff(mu,i,k-1),rnfmul(muf,gcoeff(mu,i,k)));
3398: coeff(mu,i,k-1) = ladd(temp, rnfmul(gcoeff(mu,k,k-1),gcoeff(mu,i,k)));
3399: }
3400: if (kk>2) { kk--; if (DEBUGLEVEL) fprintferr("%ld ",kk); }
3401: }
3402: else
3403: {
3404: for (l=kk-2; l; l--)
3405: {
3406: /* RED(k,l) */
3407: ideal=idealmul(nf,(GEN)I[l],Ikk_inv);
3408: x=findmin(nf,ideal,gcoeff(mu,kk,l),2*prec-2);
3409: if (!gcmp0(x))
3410: {
3411: xpol=basistoalg(nf,x); xc=nftocomplex(nf,xpol);
3412: MC[kk]=(long)gsub((GEN)MC[kk],rnfvecmul(xc,(GEN)MC[l]));
3413: U[kk]=(long)gsub((GEN)U[kk],gmul(xpol,(GEN)U[l]));
3414: coeff(mu,kk,l)=lsub(gcoeff(mu,kk,l),xc);
3415: for (i=1; i<l; i++)
3416: coeff(mu,kk,i) = lsub(gcoeff(mu,kk,i), rnfmul(xc,gcoeff(mu,l,i)));
3417: }
3418: }
3419: kk++; if (DEBUGLEVEL) fprintferr("%ld ",kk);
3420: }
3421: }
3422: while (kk<lx);
3423: if (DEBUGLEVEL) fprintferr("\n");
3424: p1=gmul(MPOL,U); tetpil=avma;
3425: y=cgetg(3,t_VEC); z=cgetg(3,t_VEC); y[1]=(long)z;
3426: z[2]=lcopy(I); z[1]=(long)algtobasis(nf,p1);
3427: y[2]=(long)algtobasis(nf,U);
3428: return gerepile(av,tetpil,y);
3429: }
3430:
3431: GEN
3432: rnfpolred(GEN nf, GEN pol, long prec)
3433: {
3434: ulong av = avma;
3435: long i,j,k,n,N, vpol = varn(pol);
3436: GEN id,id2,newid,newor,p1,p2,al,newpol,w,z;
3437: GEN bnf,zk,newideals,ideals,order,neworder;
3438:
3439: if (typ(pol)!=t_POL) err(typeer,"rnfpolred");
3440: if (typ(nf)!=t_VEC) err(idealer1);
3441: switch(lg(nf))
3442: {
3443: case 10: bnf = NULL; break;
3444: case 11: bnf = nf; nf = checknf((GEN)nf[7]); break;
3445: default: err(idealer1);
3446: return NULL; /* not reached */
3447: }
3448: if (degpol(pol) <= 1)
3449: {
3450: w=cgetg(2,t_VEC);
3451: w[1]=lpolx[vpol]; return w;
3452: }
3453: id=rnfpseudobasis(nf,pol); N=degpol(nf[1]);
3454: if (bnf && gcmp1(gmael3(bnf,8,1,1))) /* if bnf is principal */
3455: {
3456: ideals=(GEN)id[2]; n=lg(ideals)-1; order=(GEN)id[1];
3457: newideals=cgetg(n+1,t_VEC); neworder=cgetg(n+1,t_MAT);
3458: zk=idmat(N);
3459: for (j=1; j<=n; j++)
3460: {
3461: newideals[j]=(long)zk; p1=cgetg(n+1,t_COL); neworder[j]=(long)p1;
3462: p2=(GEN)order[j];
3463: al=(GEN)isprincipalgen(bnf,(GEN)ideals[j])[2];
3464: for (k=1; k<=n; k++)
3465: p1[k]=(long)element_mul(nf,(GEN)p2[k],al);
3466: }
3467: id=cgetg(3,t_VEC); id[1]=(long)neworder; id[2]=(long)newideals;
3468: }
3469: id2=rnflllgram(nf,pol,id,prec);
3470: z=(GEN)id2[1]; newid=(GEN)z[2]; newor=(GEN)z[1];
3471: n=lg(newor)-1; w=cgetg(n+1,t_VEC);
3472: for (j=1; j<=n; j++)
3473: {
3474: p1=(GEN)newid[j]; al=gmul(gcoeff(p1,1,1),(GEN)newor[j]);
3475: p1=basistoalg(nf,(GEN)al[n]);
3476: for (i=n-1; i; i--)
3477: p1=gadd(basistoalg(nf,(GEN)al[i]),gmul(polx[vpol],p1));
3478: newpol=gtopoly(gmodulcp(gtovec(caract2(lift(pol),lift(p1),vpol)),
3479: (GEN) nf[1]), vpol);
3480: p1 = ggcd(newpol, derivpol(newpol));
3481: if (degpol(p1)>0)
3482: {
3483: newpol=gdiv(newpol,p1);
3484: newpol=gdiv(newpol,leading_term(newpol));
3485: }
3486: w[j]=(long)newpol;
3487: if (DEBUGLEVEL>=4) outerr(newpol);
3488: }
3489: return gerepilecopy(av,w);
3490: }
3491:
3492: extern GEN vecpol_to_mat(GEN v, long n);
3493:
3494: /* given a relative polynomial pol over nf, compute a pseudo-basis for the
3495: * extension, then an absolute basis */
3496: GEN
3497: makebasis(GEN nf,GEN pol)
3498: {
3499: GEN elts,ids,polabs,plg,B,bs,p1,p2,a,den,vbs,vbspro,vpro,rnf;
3500: long av=avma,n,N,m,i,j, v = varn(pol);
3501:
3502: p1 = rnfequation2(nf,pol);
3503: polabs= (GEN)p1[1];
3504: plg = (GEN)p1[2];
3505: a = (GEN)p1[3];
3506: rnf = cgetg(12,t_VEC);
3507: for (i=2;i<=9;i++) rnf[i]=zero;
3508: rnf[1] =(long)pol;
3509: rnf[10]=(long)nf; p2=cgetg(4,t_VEC);
3510: rnf[11]=(long)p2; p2[1]=p2[2]=zero; p2[3]=(long)a;
3511: if (signe(a))
3512: pol = gsubst(pol,v,gsub(polx[v],
3513: gmul(a,gmodulcp(polx[varn(nf[1])],(GEN)nf[1]))));
3514: p1=rnfpseudobasis(nf,pol);
3515: elts= (GEN)p1[1];
3516: ids = (GEN)p1[2];
3517: if (DEBUGLEVEL>1) { fprintferr("relative basis computed\n"); flusherr(); }
3518: N=degpol(pol); n=degpol((GEN)nf[1]); m=n*N;
3519: den = denom(content(lift(plg)));
3520: vbs = cgetg(n+1,t_VEC);
3521: vbs[1] = un;
3522: vbs[2] = (long)plg; vbspro = gmul(den,plg);
3523: for(i=3;i<=n;i++)
3524: vbs[i] = ldiv(gmul((GEN)vbs[i-1],vbspro),den);
3525: bs = gmul(vbs, vecpol_to_mat((GEN)nf[7],n));
3526:
3527: vpro=cgetg(N+1,t_VEC);
3528: for (i=1;i<=N;i++)
3529: {
3530: p1=cgetg(3,t_POLMOD);
3531: p1[1]=(long)polabs;
3532: p1[2]=lpuigs(polx[v],i-1); vpro[i]=(long)p1;
3533: }
3534: vpro=gmul(vpro,elts); B = cgetg(m+1, t_MAT);
3535: for(i=1;i<=N;i++)
3536: for(j=1;j<=n;j++)
3537: {
3538: p1 = gmul(bs, element_mul(nf,(GEN)vpro[i],gmael(ids,i,j)));
3539: B[(i-1)*n+j] = (long)pol_to_vec(lift_intern(p1), m);
3540: }
3541: p1 = denom(B); B = gmul(B,p1);
3542: B = hnfmodid(B, p1); B = gdiv(B,p1);
3543: p1=cgetg(4,t_VEC);
3544: p1[1]=(long)polabs;
3545: p1[2]=(long)B;
3546: p1[3]=(long)rnf; return gerepilecopy(av, p1);
3547: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to base2.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