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