Annotation of OpenXM_contrib/pari/src/basemath/base4.c, Revision 1.1.1.1
1.1 maekawa 1: /*******************************************************************/
2: /* */
3: /* BASIC NF OPERATIONS */
4: /* (continued) */
5: /* */
6: /*******************************************************************/
7: /* $Id: base4.c,v 1.1.1.1 1999/09/16 13:47:22 karim Exp $ */
8: #include "pari.h"
9: #include "parinf.h"
10:
11: #define principalideal_aux(nf,x) (principalideal0((nf),(x),0))
12:
13: GEN element_muli(GEN nf, GEN x, GEN y);
14:
15: static GEN nfbezout(GEN nf, GEN a, GEN b, GEN ida, GEN idb, GEN *u, GEN *v, GEN *w, GEN *di);
16:
17: /*******************************************************************/
18: /* */
19: /* IDEAL OPERATIONS */
20: /* */
21: /*******************************************************************/
22:
23: /* A valid ideal is either principal (valid nf_element), or prime, or a matrix
24: * on the integer basis (preferably HNF).
25: * A prime ideal is of the form [p,a,e,f,b], where the ideal is p.Z_K+a.Z_K,
26: * p is a rational prime, a belongs to Z_K, e=e(P/p), f=f(P/p), and b
27: * Lenstra constant (p.P^(-1)= p Z_K + b Z_K).
28: *
29: * An idele is a couple[I,V] where I is a valid ideal and V a row vector
30: * with r1+r2 components (real or complex). For instance, if M=(a), V
31: * contains the complex logarithms of the first r1+r2 conjugates of a
32: * (depends on the chosen generator a). All subroutines work with either
33: * ideles or ideals (an omitted V is assumed to be 0).
34: *
35: * All the output ideals will be in HNF form.
36: */
37:
38: /* types and conversions */
39:
40: static long
41: idealtyp(GEN *ideal, GEN *arch)
42: {
43: GEN x = *ideal;
44: long t,lx,tx = typ(x);
45:
46: if (tx==t_VEC && lg(x)==3)
47: { *arch = (GEN)x[2]; x = (GEN)x[1]; tx = typ(x); }
48: else
49: *arch = NULL;
50: switch(tx)
51: {
52: case t_MAT: lx = lg(x);
53: if (lx>2) t = id_MAT;
54: else
55: {
56: t = id_PRINCIPAL;
57: x = (lx==2)? (GEN)x[1]: gzero;
58: }
59: break;
60:
61: case t_VEC: if (lg(x)!=6) err(idealer2);
62: t = id_PRIME; break;
63:
64: case t_POL: case t_POLMOD: case t_COL:
65: t = id_PRINCIPAL; break;
66: default:
67: if (tx!=t_INT && !is_frac_t(tx)) err(idealer2);
68: t = id_PRINCIPAL;
69: }
70: *ideal = x; return t;
71: }
72:
73: /* Assume ideal in HNF form */
74: long
75: ideal_is_zk(GEN ideal,long N)
76: {
77: long i,j, lx = lg(ideal);
78:
79: if (typ(ideal) != t_MAT || lx==1) return 0;
80: N++; if (lx != N || lg(ideal[1]) != N) return 0;
81: for (i=1; i<N; i++)
82: {
83: if (!gcmp1(gcoeff(ideal,i,i))) return 0;
84: for (j=i+1; j<N; j++)
85: if (!gcmp0(gcoeff(ideal,i,j))) return 0;
86: }
87: return 1;
88: }
89:
90: static GEN
91: prime_to_ideal_aux(GEN nf, GEN vp)
92: {
93: GEN m,el;
94: long i, N = lgef(nf[1])-3;
95:
96: m = cgetg(N+1,t_MAT); el = (GEN)vp[2];
97: for (i=1; i<=N; i++) m[i] = (long) element_mulid(nf,el,i);
98: return hnfmodid(m,(GEN)vp[1]);
99: }
100:
101: GEN
102: prime_to_ideal(GEN nf, GEN vp)
103: {
104: long av=avma;
105: if (typ(vp) == t_INT) return gscalmat(vp, lgef(nf[1])-3);
106: return gerepileupto(av, prime_to_ideal_aux(nf,vp));
107: }
108:
109: /* x = ideal in matrix form. Put it in hnf. */
110: static GEN
111: idealmat_to_hnf(GEN nf, GEN x)
112: {
113: long rx,i,j,N;
114: GEN m,dx;
115:
116: N=lgef(nf[1])-3; rx=lg(x)-1;
117: if (!rx) return gscalmat(gzero,N);
118:
119: dx=denom(x); if (gcmp1(dx)) dx = NULL; else x=gmul(dx,x);
120: if (rx >= N) m = x;
121: else
122: {
123: m=cgetg(rx*N + 1,t_MAT);
124: for (i=1; i<=rx; i++)
125: for (j=1; j<=N; j++)
126: m[(i-1)*N + j] = (long) element_mulid(nf,(GEN)x[i],j);
127: }
128: x = hnfmod(m,detint(m));
129: return dx? gdiv(x,dx): x;
130: }
131:
132: int
133: ishnfall(GEN x)
134: {
135: long i,j, lx = lg(x);
136: for (i=2; i<lx; i++)
137: {
138: if (gsigne(gcoeff(x,i,i)) <= 0) return 0;
139: for (j=1; j<i; j++)
140: if (!gcmp0(gcoeff(x,i,j))) return 0;
141: }
142: return (gsigne(gcoeff(x,1,1)) > 0);
143: }
144:
145: GEN
146: idealhermite_aux(GEN nf, GEN x)
147: {
148: long N,tx,lx;
149: GEN z;
150:
151: tx = idealtyp(&x,&z);
152: if (tx == id_PRIME) return prime_to_ideal(nf,x);
153: if (tx == id_PRINCIPAL)
154: {
155: x = principalideal(nf,x);
156: return idealmat_to_hnf(nf,x);
157: }
158: N=lgef(nf[1])-3; lx = lg(x);
159: if (lg(x[1]) != N+1) err(idealer2);
160:
161: if (lx == N+1 && ishnfall(x)) return x;
162: if (lx <= N) return idealmat_to_hnf(nf,x);
163: z=denom(x); if (gcmp1(z)) z=NULL; else x = gmul(z,x);
164: x = hnfmod(x,detint(x));
165: return z? gdiv(x,z): x;
166: }
167:
168: GEN
169: idealhermite(GEN nf, GEN x)
170: {
171: long av=avma;
172: GEN p1;
173: nf = checknf(nf); p1 = idealhermite_aux(nf,x);
174: if (p1==x || p1==(GEN)x[1]) return gcopy(p1);
175: return gerepileupto(av,p1);
176: }
177:
178: static GEN
179: principalideal0(GEN nf, GEN x, long copy)
180: {
181: GEN z = cgetg(2,t_MAT);
182: switch(typ(x))
183: {
184: case t_INT: case t_FRAC: case t_FRACN:
185: if (copy) x = gcopy(x);
186: x = gscalcol_i(x, lgef(nf[1])-3); break;
187:
188: case t_POLMOD:
189: if (!gegal((GEN)nf[1],(GEN)x[1]))
190: err(talker,"incompatible number fields in principalideal");
191: x=(GEN)x[2]; /* fall through */
192: case t_POL:
193: x = copy? algtobasis(nf,x): algtobasis_intern(nf,x);
194: break;
195:
196: case t_MAT:
197: if (lg(x)!=2) err(typeer,"principalideal");
198: x = (GEN)x[1];
199: case t_COL:
200: if (lg(x)==lgef(nf[1])-2)
201: {
202: if (copy) x = gcopy(x);
203: break;
204: }
205: default: err(typeer,"principalideal");
206: }
207: z[1]=(long)x; return z;
208: }
209:
210: GEN
211: principalideal(GEN nf, GEN x)
212: {
213: nf = checknf(nf); return principalideal0(nf,x,1);
214: }
215:
216: /* for internal use */
217: GEN
218: get_arch(GEN nf,GEN x,long prec)
219: {
220: long i,R1,RU;
221: GEN v,p1,p2;
222:
223: R1=itos(gmael(nf,2,1)); RU = R1+itos(gmael(nf,2,2));
224: if (typ(x)!=t_COL) x = algtobasis_intern(nf,x);
225: if (isnfscalar(x)) /* rational number */
226: {
227: v = cgetg(RU+1,t_VEC);
228: p1=glog((GEN)x[1],prec); if (RU!=R1) p2=gmul2n(p1,1);
229: for (i=1; i<=R1; i++) v[i]=(long)p1;
230: for ( ; i<=RU; i++) v[i]=(long)p2;
231: }
232: else
233: {
234: x = gmul(gmael(nf,5,1),x); v = cgetg(RU+1,t_VEC);
235: for (i=1; i<=R1; i++) v[i] = llog((GEN)x[i],prec);
236: for ( ; i<=RU; i++) v[i] = lmul2n(glog((GEN)x[i],prec),1);
237: }
238: return v;
239: }
240:
241: GEN
242: get_arch_real(GEN nf,GEN x,GEN *emb,long prec)
243: {
244: long i,R1,RU;
245: GEN v,p1,p2;
246:
247: R1=itos(gmael(nf,2,1)); RU = R1+itos(gmael(nf,2,2));
248: if (typ(x)!=t_COL) x = algtobasis_intern(nf,x);
249: if (isnfscalar(x)) /* rational number */
250: {
251: GEN u = (GEN)x[1];
252: v = cgetg(RU+1,t_COL);
253: i = signe(u);
254: if (!i) err(talker,"0 in get_arch_real");
255: p1= (i > 0)? glog(u,prec): gzero;
256: if (RU != R1) p2 = gmul2n(p1,1);
257: for (i=1; i<=R1; i++) v[i]=(long)p1;
258: for ( ; i<=RU; i++) v[i]=(long)p2;
259: }
260: else
261: {
262: x = gmul(gmael(nf,5,1),x); v = cgetg(RU+1,t_COL);
263: for (i=1; i<=R1; i++) v[i] = llog(gabs((GEN)x[i],prec),prec);
264: for ( ; i<=RU; i++) v[i] = llog(gnorm((GEN)x[i]),prec);
265: }
266: *emb = x; return v;
267: }
268:
269: GEN
270: principalidele(GEN nf, GEN x, long prec)
271: {
272: GEN p1,y = cgetg(3,t_VEC);
273: long av;
274:
275: nf = checknf(nf);
276: p1 = principalideal0(nf,x,1);
277: y[1] = (long)p1;
278: av =avma; p1 = get_arch(nf,(GEN)p1[1],prec);
279: y[2] = lpileupto(av,p1); return y;
280: }
281:
282: /* GP functions */
283:
284: GEN
285: ideal_two_elt0(GEN nf, GEN x, GEN a)
286: {
287: if (!a) return ideal_two_elt(nf,x);
288: return ideal_two_elt2(nf,x,a);
289: }
290:
291: GEN
292: idealpow0(GEN nf, GEN x, GEN n, long flag, long prec)
293: {
294: if (flag) return idealpowred(nf,x,n,prec);
295: return idealpow(nf,x,n);
296: }
297:
298: GEN
299: idealmul0(GEN nf, GEN x, GEN y, long flag, long prec)
300: {
301: if (flag) return idealmulred(nf,x,y,prec);
302: return idealmul(nf,x,y);
303: }
304:
305: GEN
306: idealpowred(GEN nf, GEN x, GEN n, long prec)
307: {
308: long av=avma, tetpil;
309: x = idealpow(nf,x,n); tetpil=avma;
310: return gerepile(av,tetpil, ideallllred(nf,x,NULL,prec));
311: }
312:
313: GEN
314: idealmulred(GEN nf, GEN x, GEN y, long prec)
315: {
316: long av=avma,tetpil;
317: x = idealmul(nf,x,y); tetpil=avma;
318: return gerepile(av,tetpil,ideallllred(nf,x,NULL,prec));
319: }
320:
321: GEN
322: idealinv0(GEN nf, GEN ix, long flag)
323: {
324: switch(flag)
325: {
326: case 0: return idealinv(nf,ix);
327: case 1: return oldidealinv(nf,ix);
328: default: err(flagerr,"idealinv");
329: }
330: return NULL; /* not reached */
331: }
332:
333: GEN
334: idealdiv0(GEN nf, GEN x, GEN y, long flag)
335: {
336: switch(flag)
337: {
338: case 0: return idealdiv(nf,x,y);
339: case 1: return idealdivexact(nf,x,y);
340: default: err(flagerr,"idealdiv");
341: }
342: return NULL; /* not reached */
343: }
344:
345: GEN
346: idealaddtoone0(GEN nf, GEN arg1, GEN arg2)
347: {
348: if (!arg2) return idealaddmultoone(nf,arg1);
349: return idealaddtoone(nf,arg1,arg2);
350: }
351:
352: static GEN
353: two_to_hnf(GEN nf, GEN a, GEN b)
354: {
355: a = principalideal_aux(nf,a);
356: b = principalideal_aux(nf,b);
357: a = concatsp(a,b);
358: if (lgef(nf[1])==5) /* quadratic field: a has to be turned into idealmat */
359: a = idealmul(nf,idmat(2),a);
360: return idealmat_to_hnf(nf, a);
361: }
362:
363: GEN
364: idealhnf0(GEN nf, GEN a, GEN b)
365: {
366: long av;
367: if (!b) return idealhermite(nf,a);
368:
369: /* HNF of aZ_K+bZ_K */
370: av = avma; nf=checknf(nf);
371: return gerepileupto(av, two_to_hnf(nf,a,b));
372: }
373:
374: GEN
375: idealhermite2(GEN nf, GEN a, GEN b)
376: {
377: return idealhnf0(nf,a,b);
378: }
379:
380: static GEN
381: check_elt(GEN a, GEN pol, GEN pnorm, GEN idz)
382: {
383: GEN x,norme, p2,p1;
384:
385: if (!signe(a)) return NULL;
386: p1 = content(a);
387: if (gcmp1(p1)) { x=a; p1=NULL; }
388: else { x=gdiv(a,p1); p2=gpowgs(p1, lgef(pol)-3); }
389:
390: norme = resultantducos(pol,x); if (p1) norme = gmul(norme,p2);
391: if (gcmp1(mppgcd(divii(norme,pnorm),pnorm))) return a;
392:
393: if (p1)
394: {
395: idz=gdiv(idz,p1);
396: if (typ(idz) == t_FRAC) /* should be exceedingly rare */
397: {
398: x = gmul(x,(GEN)idz[2]);
399: p2 = gdiv(p2, gpowgs((GEN)idz[2], lgef(pol)-3));
400: idz = (GEN)idz[1];
401: }
402: }
403: x = gadd(x,idz);
404:
405: norme = resultantducos(pol,x); if (p1) norme = gmul(norme,p2);
406: if (gcmp1(mppgcd(divii(norme,pnorm),pnorm))) return a;
407: return NULL;
408: }
409:
410: static void
411: setprec(GEN x, long prec)
412: {
413: long i,j, n=lg(x);
414: for (i=1;i<n;i++)
415: {
416: GEN p2,p1 = (GEN)x[i];
417: for (j=1;j<n;j++)
418: {
419: p2 = (GEN)p1[j];
420: if (typ(p2) == t_REAL) setlg(p2, prec);
421: }
422: }
423: }
424:
425: /* find a basis of x whose elements have small norm
426: * M a bound for the size of coeffs of x */
427: GEN
428: ideal_better_basis(GEN nf, GEN x, GEN M)
429: {
430: GEN a,b;
431: long nfprec = (long)nfnewprec(nf,0);
432: long prec = DEFAULTPREC + (expi(M) >> TWOPOTBITS_IN_LONG);
433:
434: if (typ(nf[5]) != t_VEC) return x;
435: if ((prec<<1) < nfprec) prec = (prec+nfprec) >> 1;
436: a = qf_base_change(gmael(nf,5,3),x,1);
437: setprec(a,prec);
438: b = lllgramintern(a,4,1,prec);
439: if (!b)
440: {
441: if (DEBUGLEVEL)
442: err(warner, "precision too low in ideal_better_basis (1)");
443: if (nfprec > prec)
444: {
445: setprec(a,nfprec);
446: b = lllgramintern(a,4,1,nfprec);
447: }
448: }
449: if (!b)
450: {
451: if (DEBUGLEVEL)
452: err(warner, "precision too low in ideal_better_basis (2)");
453: b = lllint(x); /* better than nothing */
454: }
455: return gmul(x, b);
456: }
457:
458: static GEN
459: mat_ideal_two_elt(GEN nf, GEN x)
460: {
461: GEN y,a,beta,pnorm,con,idz, pol = (GEN)nf[1];
462: long av,tetpil,i,z, N = lgef(pol)-3;
463:
464: y=cgetg(3,t_VEC); av=avma;
465: if (lg(x[1])!=N+1) err(typeer,"ideal_two_elt");
466: if (N == 2)
467: {
468: y[1] = lcopy(gcoeff(x,1,1));
469: y[2] = lcopy((GEN)x[2]); return y;
470: }
471:
472: con=content(x); if (!gcmp1(con)) x = gdiv(x,con);
473: if (lg(x) != N+1) x = idealhermite_aux(nf,x);
474: idz=gcoeff(x,1,1);
475: if (gcmp1(idz))
476: {
477: y[1]=lpileupto(av,gcopy(con));
478: y[2]=(long)gscalcol(con,N); return y;
479: }
480: pnorm = dethnf(x);
481: beta = gmul((GEN)nf[7], x);
482: for (i=2; i<=N; i++)
483: {
484: a = check_elt((GEN)beta[i], pol, pnorm, idz);
485: if (a) break;
486: }
487: if (i>N)
488: {
489: x = ideal_better_basis(nf,x,pnorm);
490: beta = gmul((GEN)nf[7], x);
491: for (i=1; i<=N; i++)
492: {
493: a = check_elt((GEN)beta[i], pol, pnorm, idz);
494: if (a) break;
495: }
496: }
497: if (i>N)
498: {
499: long c=0, av1=avma;
500:
501: if (DEBUGLEVEL>3) fprintferr("ideal_two_elt, hard case: ");
502: for(;;)
503: {
504: if (DEBUGLEVEL>3) fprintferr("%d ", ++c);
505: a = gzero;
506: for (i=1; i<=N; i++)
507: {
508: z = mymyrand() >> (BITS_IN_RANDOM-5); /* in [0,15] */
509: if (z >= 9) z -= 7;
510: a = gadd(a,gmulsg(z,(GEN)beta[i]));
511: }
512: a = check_elt(a, pol, pnorm, idz);
513: if (a) break;
514: avma=av1;
515: }
516: if (DEBUGLEVEL>3) fprintferr("\n");
517: }
518: a = centermod(algtobasis_intern(nf,a), idz);
519: tetpil=avma; y[1]=lmul(idz,con); y[2]=lmul(a,con);
520: gerepilemanyvec(av,tetpil,y+1,2); return y;
521: }
522:
523: /* Etant donne un ideal ix, ressort un vecteur [a,alpha] a deux composantes
524: * tel que a soit rationnel et ix=aZ_K+alpha Z_K, alpha etant un vecteur
525: * colonne sur la base d'entiers. On peut avoir a=0 ou alpha=0, mais on ne
526: * cherche pas a determiner si ix est principal.
527: */
528: GEN
529: ideal_two_elt(GEN nf, GEN x)
530: {
531: GEN z;
532: long N, tx = idealtyp(&x,&z);
533:
534: nf=checknf(nf);
535: if (tx==id_MAT) return mat_ideal_two_elt(nf,x);
536:
537: N=lgef(nf[1])-3; z=cgetg(3,t_VEC);
538: if (tx == id_PRINCIPAL)
539: {
540: switch(typ(x))
541: {
542: case t_INT: case t_FRAC: case t_FRACN:
543: z[1]=lcopy(x);
544: z[2]=(long)zerocol(N); return z;
545:
546: case t_POLMOD:
547: if (!gegal((GEN)nf[1],(GEN)x[1]))
548: err(talker,"incompatible number fields in ideal_two_elt");
549: x=(GEN)x[2]; /* fall through */
550: case t_POL:
551: z[1]=zero; z[2]=(long)algtobasis(nf,x); return z;
552: case t_COL:
553: if (lg(x)==N+1) { z[1]=zero; z[2]=lcopy(x); return z; }
554: }
555: }
556: else if (tx == id_PRIME)
557: {
558: z[1]=lcopy((GEN)x[1]);
559: z[2]=lcopy((GEN)x[2]); return z;
560: }
561: err(typeer,"ideal_two_elt");
562: return NULL; /* not reached */
563: }
564:
565: /* factorization */
566:
567: GEN
568: idealfactor(GEN nf, GEN x)
569: {
570: long av,tx, tetpil,i,j,k,lf,lff,N,ls,v,vd;
571: GEN d,f,f1,f2,ff,ff1,ff2,y1,y2,y,p1,p2,denx;
572:
573: tx = idealtyp(&x,&y);
574: if (tx == id_PRIME)
575: {
576: y=cgetg(3,t_MAT);
577: y[1]=lgetg(2,t_COL); mael(y,1,1)=lcopy(x);
578: y[2]=lgetg(2,t_COL); mael(y,2,1)=un; return y;
579: }
580: nf=checknf(nf); av=avma;
581: if (tx == id_PRINCIPAL) x = principalideal_aux(nf,x);
582:
583: N=lgef(nf[1])-3; if (lg(x) != N+1) x = idealmat_to_hnf(nf,x);
584: if (lg(x)==1) err(talker,"zero ideal in idealfactor");
585: denx=denom(x);
586: if (gcmp1(denx)) lff=1;
587: else
588: {
589: ff=factor(denx); ff1=(GEN)ff[1]; ff2=(GEN)ff[2];
590: lff=lg(ff1); x=gmul(denx,x);
591: }
592: for (d=gun,i=1; i<=N; i++) d=mulii(d,gcoeff(x,i,i));
593: f=factor(absi(d)); f1=(GEN)f[1]; f2=(GEN)f[2]; lf=lg(f1);
594: y1=cgetg((lf+lff-2)*N+1,t_COL); y2=new_chunk((lf+lff-2)*N+1);
595: k=0;
596: for (i=1; i<lf; i++)
597: {
598: p1=primedec(nf,(GEN)f1[i]); ls=itos((GEN)f2[i]);
599: vd=ggval(denx,(GEN)f1[i]);
600: for (j=1; j<lg(p1); j++)
601: {
602: p2=(GEN)p1[j];
603: if (ls)
604: {
605: v = idealval(nf,x,p2);
606: ls -= v*itos((GEN)p2[4]);
607: v -= vd*itos((GEN)p2[3]);
608: }
609: else v = - vd*itos((GEN)p2[3]);
610: if (v) { y1[++k]=(long)p2; y2[k]=v; }
611: }
612: }
613: for (i=1; i<lff; i++)
614: if (!divise(d,(GEN)ff1[i]))
615: {
616: p1=primedec(nf,(GEN)ff1[i]);
617: for (j=1; j<lg(p1); j++)
618: {
619: p2=(GEN)p1[j]; y1[++k]=(long)p2;
620: y2[k] = -itos((GEN)ff2[i])*itos((GEN)p2[3]);
621: }
622: }
623: tetpil=avma; y=cgetg(3,t_MAT);
624: p1=cgetg(k+1,t_COL); y[1]=(long)p1;
625: p2=cgetg(k+1,t_COL); y[2]=(long)p2;
626: for (i=1; i<=k; i++) { p1[i]=lcopy((GEN)y1[i]); p2[i]=lstoi(y2[i]); }
627: return gerepile(av,tetpil,y);
628: }
629:
630: /* vp prime ideal in primedec format. Return valuation(ix) at vp */
631: long
632: idealval(GEN nf, GEN ix, GEN vp)
633: {
634: long N,v,vd,w,av=avma,av1,lim,i,j,k, tx = typ(ix);
635: GEN mul,mat,a,x,y,r,bp,p,denx;
636:
637: nf=checknf(nf); checkprimeid(vp);
638: if (is_extscalar_t(tx) || tx==t_COL) return element_val(nf,ix,vp);
639: p=(GEN)vp[1]; N=lgef(nf[1])-3;
640: tx = idealtyp(&ix,&a);
641: denx=denom(ix); if (!gcmp1(denx)) ix=gmul(denx,ix);
642: if (tx != id_MAT)
643: ix = idealhermite_aux(nf,ix);
644: else
645: {
646: checkid(ix,N);
647: if (lg(ix) != N+1) ix=idealmat_to_hnf(nf,ix);
648: }
649: v = ggval(dethnf_i(ix), p);
650: vd = ggval(denx,p) * itos((GEN)vp[3]); /* v_p * e */
651: if (!v) return -vd;
652:
653: mul = cgetg(N+1,t_MAT); bp=(GEN)vp[5];
654: mat = cgetg(N+1,t_MAT);
655: for (j=1; j<=N; j++)
656: {
657: mul[j] = (long)element_mulid(nf,bp,j);
658: x = (GEN)ix[j];
659: y = cgetg(N+1, t_COL); mat[j] = (long)y;
660: for (i=1; i<=N; i++)
661: { /* compute (x.b)_i, ix in HNF ==> x[j+1..N] = 0 */
662: a = mulii((GEN)x[1], gcoeff(mul,i,1));
663: for (k=2; k<=j; k++) a = addii(a, mulii((GEN)x[k], gcoeff(mul,i,k)));
664: /* is it divisible by p ? */
665: y[i] = ldvmdii(a,p,&r);
666: if (signe(r)) { avma=av; return -vd; }
667: }
668: }
669: av1 = avma; lim=stack_lim(av1,3);
670: y = cgetg(N+1,t_COL);
671: for (w=1; w<v; w++)
672: for (j=1; j<=N; j++)
673: {
674: x = (GEN)mat[j];
675: for (i=1; i<=N; i++)
676: { /* compute (x.b)_i */
677: a = mulii((GEN)x[1], gcoeff(mul,i,1));
678: for (k=2; k<=N; k++) a = addii(a, mulii((GEN)x[k], gcoeff(mul,i,k)));
679: /* is it divisible by p ? */
680: y[i] = ldvmdii(a,p,&r);
681: if (signe(r)) { avma=av; return w - vd; }
682: }
683: r=x; mat[j]=(long)y; y=r;
684: if (low_stack(lim,stack_lim(av1,3)))
685: {
686: GEN *gptr[2]; gptr[0]=&y; gptr[1]=&mat;
687: if(DEBUGMEM>1) err(warnmem,"idealval");
688: gerepilemany(av1,gptr,2);
689: }
690: }
691: avma=av; return w - vd;
692: }
693:
694: /* gcd and generalized Bezout */
695:
696: GEN
697: idealadd(GEN nf, GEN x, GEN y)
698: {
699: long av=avma,N,tx,ty;
700: GEN z,p1,dx,dy,dz;
701:
702: tx = idealtyp(&x,&z);
703: ty = idealtyp(&y,&z);
704: nf=checknf(nf); N=lgef(nf[1])-3;
705: z = cgetg(N+1, t_MAT);
706: if (tx != id_MAT || lg(x)!=N+1) x = idealhermite_aux(nf,x);
707: if (ty != id_MAT || lg(y)!=N+1) y = idealhermite_aux(nf,y);
708: if (lg(x) == 1) return gerepileupto(av,y);
709: if (lg(y) == 1) return gerepileupto(av,x); /* check for 0 ideal */
710: dx=denom(x);
711: dy=denom(y); dz=mulii(dx,dy);
712: if (gcmp1(dz)) dz = NULL; else { x=gmul(x,dz); y=gmul(y,dz); }
713: p1=mppgcd(detint(x),detint(y));
714: if (gcmp1(p1))
715: {
716: long i,j;
717: if (!dz) { avma=av; return idmat(N); }
718: avma = (long)dz; dz = gerepileupto((long)z, ginv(dz));
719: for (i=1; i<=N; i++)
720: {
721: z[i]=lgetg(N+1,t_COL);
722: for (j=1; j<=N; j++)
723: coeff(z,j,i) = (i==j)? (long)dz: zero;
724: }
725: return z;
726: }
727: z=hnfmod(concatsp(x,y),p1); if (dz) z=gdiv(z,dz);
728: return gerepileupto(av,z);
729: }
730:
731: static GEN
732: get_p1(GEN nf, GEN x, GEN y,long fl)
733: {
734: GEN u,v,v1,v2,v3,v4;
735: long i,j,N;
736:
737: switch(fl)
738: {
739: case 1:
740: v1 = gcoeff(x,1,1);
741: v2 = gcoeff(y,1,1);
742: if (typ(v1)!=t_INT || typ(v2)!=t_INT)
743: err(talker,"ideals don't sum to Z_K in idealaddtoone");
744: if (gcmp1(bezout(v1,v2,&u,&v)))
745: return gmul(u,(GEN)x[1]);
746: default:
747: v=hnfperm(concatsp(x,y));
748: v1=(GEN)v[1]; v2=(GEN)v[2]; v3=(GEN)v[3];
749: j=0; N = lgef(nf[1])-3;
750: for (i=1; i<=N; i++)
751: {
752: if (!gcmp1(gcoeff(v1,i,i)))
753: err(talker,"ideals don't sum to Z_K in idealaddtoone");
754: if (gcmp1((GEN)v3[i])) j=i;
755: }
756: v4=(GEN)v2[N+j]; setlg(v4,N+1);
757: return gmul(x,v4);
758: }
759: }
760:
761: GEN
762: idealaddtoone_i(GEN nf, GEN x, GEN y)
763: {
764: long t, fl = 1;
765: GEN p1,xh,yh;
766:
767: if (DEBUGLEVEL>4)
768: {
769: fprintferr(" entering idealaddtoone:\n");
770: fprintferr(" x = %Z\n",x);
771: fprintferr(" y = %Z\n",y);
772: }
773: t = idealtyp(&x,&p1);
774: if (t != id_MAT || lg(x) != lg(x[1])) xh = idealhermite_aux(nf,x);
775: else
776: { xh=x; fl = isnfscalar((GEN)x[1]); }
777: t = idealtyp(&y,&p1);
778: if (t != id_MAT || lg(y) != lg(y[1])) yh = idealhermite_aux(nf,y);
779: else
780: { yh=y; if (fl) fl = isnfscalar((GEN)y[1]); }
781:
782: p1 = get_p1(nf,xh,yh,fl);
783: p1 = element_reduce(nf,p1, idealmullll(nf,x,y));
784: if (DEBUGLEVEL>4 && !gcmp0(p1))
785: fprintferr(" leaving idealaddtoone: %Z\n",p1);
786: return p1;
787: }
788:
789: /* ideal should be an idele (not mandatory). For internal use. */
790: GEN
791: ideleaddone_aux(GEN nf,GEN x,GEN ideal)
792: {
793: long i,nba,R1;
794: GEN p1,p2,p3,arch;
795:
796: idealtyp(&ideal,&arch);
797: if (!arch) return idealaddtoone_i(nf,x,ideal);
798:
799: R1=itos(gmael(nf,2,1));
800: if (typ(arch)!=t_VEC && lg(arch)!=R1+1)
801: err(talker,"incorrect idele in idealaddtoone");
802: for (nba=0,i=1; i<lg(arch); i++)
803: if (signe(arch[i])) nba++;
804: if (!nba) return idealaddtoone_i(nf,x,ideal);
805:
806: p3 = idealaddtoone_i(nf,x,ideal);
807: if (gcmp0(p3)) p3=(GEN)idealhermite_aux(nf,x)[1];
808: p1=idealmullll(nf,x,ideal);
809:
810: p2=zarchstar(nf,p1,arch,nba);
811: p1=lift_intern(gmul((GEN)p2[3],zsigne(nf,p3,arch)));
812: p2=(GEN)p2[2]; nba=0;
813: for (i=1; i<lg(p1); i++)
814: if (signe(p1[i])) { nba=1; p3=element_mul(nf,p3,(GEN)p2[i]); }
815: if (gcmp0(p3)) return gcopy((GEN)x[1]); /* can happen if ideal = Z_K */
816: return nba? p3: gcopy(p3);
817: }
818:
819: static GEN
820: unnf_minus_x(GEN x)
821: {
822: long i, N = lg(x);
823: GEN y = cgetg(N,t_COL);
824:
825: y[1] = lsub(gun,(GEN)x[1]);
826: for (i=2;i<N; i++) y[i] = lneg((GEN)x[i]);
827: return y;
828: }
829:
830: static GEN
831: addone(GEN f(GEN,GEN,GEN), GEN nf, GEN x, GEN y)
832: {
833: GEN z = cgetg(3,t_VEC);
834: long av=avma;
835:
836: nf=checknf(nf); x = gerepileupto(av, f(nf,x,y));
837: z[1]=(long)x; z[2]=(long)unnf_minus_x(x); return z;
838: }
839:
840: GEN
841: idealaddtoone(GEN nf, GEN x, GEN y)
842: {
843: return addone(idealaddtoone_i,nf,x,y);
844: }
845:
846: GEN
847: ideleaddone(GEN nf,GEN x,GEN idele)
848: {
849: return addone(ideleaddone_aux,nf,x,idele);
850: }
851:
852: GEN
853: nfmodprinit(GEN nf, GEN pr)
854: {
855: long av;
856: GEN p,e,p1,prhall;
857:
858: nf = checknf(nf); checkprimeid(pr);
859: prhall = cgetg(3,t_VEC);
860: prhall[1] = (long) prime_to_ideal(nf,pr);
861:
862: av = avma; p = (GEN)pr[1]; e = (GEN)pr[3];
863: p1 = cgetg(2,t_MAT);
864: p1[1] = ldiv(element_pow(nf,(GEN)pr[5],e), gpuigs(p,itos(e)-1));
865: p1 = hnfmodid(idealhermite_aux(nf,p1), p);
866: p1 = idealaddtoone_i(nf,pr,p1);
867:
868: /* p1 = 1 mod pr, p1 = 0 mod q^{e_q} for all other primes q | p */
869: prhall[2] = lpileupto(av, unnf_minus_x(p1)); return prhall;
870: }
871:
872: /* given an element x in Z_K and an integral ideal y with x, y coprime,
873: outputs an element inverse of x modulo y */
874: GEN
875: element_invmodideal(GEN nf, GEN x, GEN y)
876: {
877: long av=avma,N,i, fl = 1;
878: GEN v,p1,xh,yh;
879:
880: nf=checknf(nf); N=lgef(nf[1])-3;
881: if (ideal_is_zk(y,N)) return zerocol(N);
882: if (DEBUGLEVEL>4)
883: {
884: fprintferr(" entree dans element_invmodideal() :\n");
885: fprintferr(" x = "); outerr(x);
886: fprintferr(" y = "); outerr(y);
887: }
888: i = lg(y);
889: if (typ(y)!=t_MAT || i==1 || i != lg(y[1])) yh=idealhermite_aux(nf,y);
890: else
891: { yh=y; fl = isnfscalar((GEN)y[1]); }
892: switch (typ(x))
893: {
894: case t_POL: case t_POLMOD: case t_COL:
895: xh = idealhermite_aux(nf,x); break;
896: default: err(typeer,"element_invmodideal");
897: }
898: p1 = get_p1(nf,xh,yh,fl);
899: p1 = element_div(nf,p1,x);
900: v = gerepileupto(av, nfreducemodideal(nf,p1,y));
901: if (DEBUGLEVEL>2)
902: { fprintferr(" sortie de element_invmodideal : v = "); outerr(v); }
903: return v;
904: }
905:
906: GEN
907: idealaddmultoone(GEN nf, GEN list)
908: {
909: long av=avma,tetpil,N,i,i1,j,k;
910: GEN z,v,v1,v2,v3,p1;
911:
912: nf=checknf(nf); N=lgef(nf[1])-3;
913: if (DEBUGLEVEL>4)
914: {
915: fprintferr(" entree dans idealaddmultoone() :\n");
916: fprintferr(" list = "); outerr(list);
917: }
918: if (typ(list)!=t_VEC && typ(list)!=t_COL)
919: err(talker,"not a list in idealaddmultoone");
920: k=lg(list); z=cgetg(1,t_MAT); list = dummycopy(list);
921: if (k==1) err(talker,"ideals don't sum to Z_K in idealaddmultoone");
922: for (i=1; i<k; i++)
923: {
924: p1=(GEN)list[i];
925: if (typ(p1)!=t_MAT || lg(p1)!=lg(p1[1]))
926: list[i] = (long)idealhermite_aux(nf,p1);
927: z = concatsp(z,(GEN)list[i]);
928: }
929: v=hnfperm(z); v1=(GEN)v[1]; v2=(GEN)v[2]; v3=(GEN)v[3]; j=0;
930: for (i=1; i<=N; i++)
931: {
932: if (!gcmp1(gcoeff(v1,i,i)))
933: err(talker,"ideals don't sum to Z_K in idealaddmultoone");
934: if (gcmp1((GEN)v3[i])) j=i;
935: }
936:
937: v=(GEN)v2[(k-2)*N+j]; z=cgetg(k,t_VEC);
938: for (i=1; i<k; i++)
939: {
940: p1=cgetg(N+1,t_COL); z[i]=(long)p1;
941: for (i1=1; i1<=N; i1++) p1[i1]=v[(i-1)*N+i1];
942: }
943: tetpil=avma; v=cgetg(k,typ(list));
944: for (i=1; i<k; i++) v[i]=lmul((GEN)list[i],(GEN)z[i]);
945: if (DEBUGLEVEL>2)
946: { fprintferr(" sortie de idealaddmultoone v = "); outerr(v); }
947: return gerepile(av,tetpil,v);
948: }
949:
950: /* multiplication */
951:
952: /* x integral ideal (without archimedean component) in HNF form
953: * [a,alpha,n] corresponds to the ideal aZ_K+alpha Z_K of norm n (a is a
954: * rational integer). Multiply them
955: */
956: static GEN
957: idealmulspec(GEN nf, GEN x, GEN a, GEN alpha, GEN n)
958: {
959: long i, N=lg(x)-1;
960: GEN m;
961:
962: if (isnfscalar(alpha))
963: return gmul(mppgcd(a,(GEN)alpha[1]),x);
964: m = cgetg((N<<1)+1,t_MAT);
965: for (i=1; i<=N; i++) n = mulii(n,gcoeff(x,i,i));
966: for (i=1; i<=N; i++) m[i]=(long)element_muli(nf,alpha,(GEN)x[i]);
967: for (i=1; i<=N; i++) m[i+N]=lmul(a,(GEN)x[i]);
968: return hnfmod(m,n);
969: }
970:
971: /* x ideal (matrix form,maximal rank), vp prime ideal (primedec). Output the
972: * product. Can be used for arbitrary vp of the form [p,a,e,f,b], IF vp
973: * =pZ_K+aZ_K, p is an integer, and norm(vp) = p^f; e and b are not used. For
974: * internal use.
975: */
976: GEN
977: idealmulprime(GEN nf, GEN x, GEN vp)
978: {
979: GEN dx, denx = denom(x);
980:
981: if (gcmp1(denx)) denx = NULL; else x=gmul(denx,x);
982: dx = powgi((GEN)vp[1], (GEN)vp[4]);
983: x = idealmulspec(nf,x, (GEN)vp[1], (GEN)vp[2], dx);
984: return denx? gdiv(x,denx): x;
985: }
986:
987: /* Assume ix and iy are integral in HNF form (or ideles of the same form).
988: * Ideal in iy can be of the form [a,b,N], with iy = (a,b) and N = norm y
989: * For internal use. */
990: GEN
991: idealmulh(GEN nf, GEN ix, GEN iy)
992: {
993: long N,i, f = 0;
994: GEN res,x,y,dy;
995:
996: if (typ(ix)==t_VEC) {f=1; x=(GEN)ix[1];} else x=ix;
997: if (typ(iy)==t_VEC && lg(iy)==3) {f+=2; y=(GEN)iy[1];} else y=iy;
998: if (f) res = cgetg(3,t_VEC);
999:
1000: if (typ(y)==t_VEC)
1001: y = idealmulspec(nf,x,(GEN)y[1],(GEN)y[2],(GEN)y[3]);
1002: else
1003: {
1004:
1005: N=lg(x)-1; dy=gcoeff(y,1,1);
1006: for (i=2; i<=N; i++) dy=mulii(dy,gcoeff(y,i,i));
1007: y = ideal_two_elt(nf,y);
1008: y = idealmulspec(nf,x,(GEN)y[1],(GEN)y[2],dy);
1009: }
1010:
1011: if (!f) return y;
1012: res[1]=(long)y;
1013: if (f==3) y = gadd((GEN)ix[2],(GEN)iy[2]);
1014: else
1015: {
1016: y = (f==2)? (GEN)iy[2]: (GEN)ix[2];
1017: y = gcopy(y);
1018: }
1019: res[2]=(long)y; return res;
1020: }
1021:
1022: /* x and y are ideals in matrix form */
1023: static GEN
1024: idealmat_mul(GEN nf, GEN x, GEN y)
1025: {
1026: long i,j, rx=lg(x)-1, ry=lg(y)-1;
1027: GEN dx,dy,m;
1028:
1029: dx=denom(x); if (!gcmp1(dx)) x=gmul(dx,x);
1030: dy=denom(y); if (!gcmp1(dy)) y=gmul(dy,y);
1031: dx = mulii(dx,dy);
1032: if (rx<=2 || ry<=2)
1033: {
1034: m=cgetg(rx*ry+1,t_MAT);
1035: for (i=1; i<=rx; i++)
1036: for (j=1; j<=ry; j++)
1037: m[(i-1)*ry+j]=(long)element_muli(nf,(GEN)x[i],(GEN)y[j]);
1038: y=hnfmod(m, detint(m));
1039: }
1040: else
1041: {
1042: x=idealmat_to_hnf(nf,x);
1043: y=idealmat_to_hnf(nf,y); y=idealmulh(nf,x,y);
1044: }
1045: return gcmp1(dx)? y: gdiv(y,dx);
1046: }
1047:
1048: /* y is principal */
1049: static GEN
1050: add_arch(GEN nf, GEN ax, GEN y)
1051: {
1052: long tetpil, av=avma, prec=precision(ax);
1053:
1054: y = get_arch(nf,y,prec); tetpil=avma;
1055: return gerepile(av,tetpil,gadd(ax,y));
1056: }
1057:
1058: /* output the ideal product ix.iy (don't reduce) */
1059: GEN
1060: idealmul(GEN nf, GEN x, GEN y)
1061: {
1062: long tx,ty,av,f;
1063: GEN res,ax,ay,p1;
1064:
1065: tx = idealtyp(&x,&ax);
1066: ty = idealtyp(&y,&ay);
1067: if (tx>ty) {
1068: res=ax; ax=ay; ay=res;
1069: res=x; x=y; y=res;
1070: f=tx; tx=ty; ty=f;
1071: }
1072: f = (ax||ay); if (f) res = cgetg(3,t_VEC); /* product is an idele */
1073: nf=checknf(nf); av=avma;
1074: switch(tx)
1075: {
1076: case id_PRINCIPAL:
1077: switch(ty)
1078: {
1079: case id_PRINCIPAL:
1080: p1 = idealhermite_aux(nf, element_mul(nf,x,y));
1081: break;
1082: case id_PRIME:
1083: p1 = gmul((GEN)y[1],x);
1084: p1 = two_to_hnf(nf,p1, element_mul(nf,(GEN)y[2],x));
1085: break;
1086: default: /* id_MAT */
1087: p1 = idealmat_mul(nf,y, principalideal_aux(nf,x));
1088: }break;
1089:
1090: case id_PRIME:
1091: p1 = (ty==id_PRIME)? prime_to_ideal_aux(nf,y)
1092: : idealmat_to_hnf(nf,y);
1093: p1 = idealmulprime(nf,p1,x); break;
1094:
1095: default: /* id_MAT */
1096: p1 = idealmat_mul(nf,x,y);
1097: }
1098: p1 = gerepileupto(av,p1);
1099: if (!f) return p1;
1100:
1101: if (ax && ay) ax = gadd(ax,ay);
1102: else
1103: {
1104: if (ax)
1105: ax = (ty==id_PRINCIPAL)? add_arch(nf,ax,y): gcopy(ax);
1106: else
1107: ax = (tx==id_PRINCIPAL)? add_arch(nf,ay,x): gcopy(ay);
1108: }
1109: res[1]=(long)p1; res[2]=(long)ax; return res;
1110: }
1111:
1112: /* different of nf */
1113: GEN
1114: differente(GEN nf, GEN premiers)
1115: {
1116: long av=avma,i,j,vi,ei,v,nb_p,vpc;
1117: GEN ideal,diff,liste_id,p1,pcon,pr,pol,a,alpha;
1118:
1119: pol=(GEN)nf[1];
1120: if (DEBUGLEVEL>1) fprintferr("Computing different\n");
1121: if (gcmp1((GEN)nf[4]))
1122: {
1123: p1 = derivpol(pol);
1124: return gerepileupto(av,idealhermite_aux(nf,p1));
1125: }
1126:
1127: ideal = gmul((GEN)nf[3],ginv(gmael(nf,5,4)));
1128: pcon = content(ideal);
1129: if (!gcmp1(pcon)) ideal=gdiv(ideal,pcon);
1130:
1131: ideal=hnfmodid(ideal,divii((GEN)nf[3],pcon));
1132: if (DEBUGLEVEL>1) msgtimer("hnf(D*delta^-1)");
1133:
1134: if (!premiers)
1135: {
1136: premiers=factor(absi((GEN)nf[3]));
1137: if (DEBUGLEVEL>1) msgtimer("factor(D)");
1138: }
1139: diff=idmat(lgef(nf[1])-3); nb_p=lg(premiers[1]);
1140:
1141: for (i=1; i<nb_p; i++)
1142: {
1143: pr=gcoeff(premiers,i,1); liste_id = primedec(nf,pr);
1144: vi=itos(gcoeff(premiers,i,2)); vpc=ggval(pcon,pr);
1145: for (j=1; j<lg(liste_id); j++)
1146: {
1147: p1=(GEN)liste_id[j]; ei=itos((GEN)p1[3]);
1148: if (ei>1)
1149: {
1150: if (DEBUGLEVEL>1) fprintferr("treating %Z\n",p1);
1151: if (signe(ressi(ei,pr)))
1152: v = ei-1;
1153: else
1154: v = ei*(vi-vpc)-idealval(nf,ideal,p1);
1155: a = gpuigs(pr, 1+(v-1)/ei);
1156: alpha = element_pow(nf,(GEN)p1[2],stoi(v));
1157: v *= itos((GEN)p1[4]);
1158: diff = idealmulspec(nf,diff,a,alpha,gpuigs(pr,v));
1159: }
1160: }
1161: }
1162: return gerepileupto(av,diff);
1163: }
1164:
1165: /* norm of an ideal */
1166: GEN
1167: idealnorm(GEN nf, GEN x)
1168: {
1169: long av = avma,tetpil;
1170: GEN y;
1171:
1172: nf = checknf(nf);
1173: switch(idealtyp(&x,&y))
1174: {
1175: case id_PRIME:
1176: return powgi((GEN)x[1],(GEN)x[4]);
1177: case id_PRINCIPAL:
1178: x = gnorm(basistoalg(nf,x)); break;
1179: default:
1180: if (lg(x) != lgef(nf[1])-2) x = idealhermite_aux(nf,x);
1181: x = det(x);
1182: }
1183: tetpil=avma; return gerepile(av,tetpil,gabs(x,0));
1184: }
1185:
1186: /* inverse */
1187:
1188: /* P.M & M.H. */
1189: static GEN
1190: hnfideal_inv(GEN nf, GEN x)
1191: {
1192: long N = lgef(nf[1])-3;
1193: GEN denx = denom(x), detx,dual,p1;
1194:
1195: if (gcmp1(denx)) denx = NULL; else x = gmul(x,denx);
1196: detx = dethnf(x);
1197: if (gcmp0(detx)) err(talker, "cannot invert zero ideal");
1198: x = idealmulh(nf,x, gmael(nf,5,7));
1199: dual = gauss(x, gmul(detx, gmael(nf,5,6)));
1200: dual = gdiv(gtrans(dual), detx);
1201:
1202: /* nf[5][4] is a dense symmetric matrix. We computed
1203: * nf[5][6] = fieldd * ginv(nf[5][4]) in initalg.
1204: * x is upper triangular (HNF), and easily inverted.
1205: * The factor fieldd cancels while solving the linear equations.
1206: */
1207: p1 = denom(dual); dual = gmul(dual, p1);
1208: dual = hnfmod(dual, gdiv(gpowgs(p1,N),detx));
1209: if (denx) p1 = gdiv(p1,denx);
1210: return gdiv(dual,p1);
1211: }
1212:
1213: /* Calcule le dual de mat_id pour la forme trace */
1214: GEN
1215: oldidealinv(GEN nf, GEN x)
1216: {
1217: GEN res,dual,di,ax;
1218: long av,tetpil, tx = idealtyp(&x,&ax);
1219:
1220: if (tx!=id_MAT) return idealinv(nf,x);
1221: if (ax) res = cgetg(3,t_VEC);
1222: nf=checknf(nf); av=avma;
1223: if (lg(x)!=lgef(nf[1])) x = idealmat_to_hnf(nf,x);
1224:
1225: dual = ginv(gmul(gtrans(x), gmael(nf,5,4)));
1226: di=denom(dual); dual=gmul(di,dual);
1227: dual = idealmat_mul(nf,gmael(nf,5,5), dual);
1228: tetpil=avma; dual = gerepile(av,tetpil,gdiv(dual,di));
1229: if (!ax) return dual;
1230: res[1]=(long)dual; res[2]=lneg(ax); return res;
1231: }
1232:
1233: /* Calcule le dual de mat_id pour la forme trace */
1234: GEN
1235: idealinv(GEN nf, GEN x)
1236: {
1237: GEN res,ax,p1;
1238: long av=avma, tx = idealtyp(&x,&ax);
1239:
1240: if (ax) res = cgetg(3,t_VEC);
1241: nf=checknf(nf); av=avma;
1242: switch (tx)
1243: {
1244: case id_MAT:
1245: if (lg(x) != lg(x[1])) x = idealmat_to_hnf(nf,x);
1246: x = hnfideal_inv(nf,x); break;
1247: case id_PRINCIPAL: tx = typ(x);
1248: if (is_const_t(tx)) x = ginv(x);
1249: else
1250: {
1251: switch(tx)
1252: {
1253: case t_COL: x = gmul((GEN)nf[7],x); break;
1254: case t_POLMOD: x = (GEN)x[2]; break;
1255: }
1256: x = ginvmod(x,(GEN)nf[1]);
1257: }
1258: x = idealhermite_aux(nf,x); break;
1259: case id_PRIME:
1260: p1=cgetg(6,t_VEC); p1[1]=x[1]; p1[2]=x[5];
1261: p1[3]=p1[5]=zero; p1[4]=lsubsi(lgef(nf[1])-3, (GEN)x[4]);
1262: x = gdiv(prime_to_ideal_aux(nf,p1), (GEN)x[1]);
1263: }
1264: x = gerepileupto(av,x); if (!ax) return x;
1265: res[1]=(long)x; res[2]=lneg(ax); return res;
1266: }
1267:
1268: static GEN
1269: idealpowprime(GEN nf, GEN vp, GEN n)
1270: {
1271: GEN n1, x = dummycopy(vp);
1272: long s = signe(n);
1273:
1274: if (s < 0) n=negi(n);
1275: n1 = gceil(gdiv(n,(GEN)x[3]));
1276: x[1]=(long)powgi((GEN)x[1],n1);
1277: if (s < 0)
1278: x[2]=ldiv(element_pow(nf,(GEN)x[5],n), powgi((GEN)vp[1],subii(n,n1)));
1279: else
1280: x[2]=(long)element_pow(nf,(GEN)x[2],n);
1281: x = prime_to_ideal_aux(nf,x);
1282: if (s<0) x = gdiv(x, powgi((GEN)vp[1],n1));
1283: return x;
1284: }
1285:
1286: /* raise the ideal x to the power n (in Z) */
1287: GEN
1288: idealpow(GEN nf, GEN x, GEN n)
1289: {
1290: long tx,N,av,s,i;
1291: GEN res,ax,m,denx,denz,dx,n1,a,alpha;
1292:
1293: if (typ(n) != t_INT) err(talker,"non-integral exponent in idealpow");
1294: tx = idealtyp(&x,&ax);
1295: if (ax) res = cgetg(3,t_VEC);
1296: nf = checknf(nf);
1297: av=avma; N=lgef(nf[1])-3; s=signe(n);
1298: if (!s) x = idmat(N);
1299: else
1300: switch(tx)
1301: {
1302: case id_PRINCIPAL: tx = typ(x);
1303: if (!is_const_t(tx))
1304: switch(tx)
1305: {
1306: case t_COL: x = gmul((GEN)nf[7],x);
1307: case t_POL: x = gmodulcp(x,(GEN)nf[1]);
1308: }
1309: x = gpui(x,n,0);
1310: x = idealhermite_aux(nf,x); break;
1311: case id_PRIME:
1312: x = idealpowprime(nf,x,n); break;
1313: default:
1314: n1 = (s<0)? negi(n): n;
1315:
1316: denx=denom(x); if (gcmp1(denx)) denx=NULL; else x = gmul(x,denx);
1317: a=ideal_two_elt(nf,x); alpha=(GEN)a[2]; a=(GEN)a[1];
1318: dx=gcoeff(x,1,1); for (i=2; i<=N; i++) dx=mulii(dx,gcoeff(x,i,i));
1319:
1320: m = cgetg(N+1,t_MAT); a = gpui(a,n1,0);
1321: alpha = element_pow(nf,alpha,n1);
1322: for (i=1; i<=N; i++) m[i]=(long)element_mulid(nf,alpha,i);
1323: m = concatsp(m, gscalmat(a,N));
1324: x = hnfmod(m, gpui(dx,n1,0));
1325: if (s<0) x = hnfideal_inv(nf,x);
1326: if (denx) { denz=gpui(denx,negi(n),0); x=gmul(denz,x); }
1327: }
1328: x = gerepileupto(av, x);
1329: if (!ax) return x;
1330: res[1]=(long)x; res[2]=lmul(n,ax); return res;
1331: }
1332:
1333: /* Return ideal^e in number field nf. e is a C integer. */
1334: GEN
1335: idealpows(GEN nf, GEN ideal, long e)
1336: {
1337: long court[] = {evaltyp(t_INT) | m_evallg(3),0,0};
1338: affsi(e,court); return idealpow(nf,ideal,court);
1339: }
1340:
1341: /* compute vp^n (vp prime, n integer), reducing along the way if n is big */
1342: GEN
1343: idealpowred_prime(GEN nf, GEN vp, GEN n, long prec)
1344: {
1345: long av=avma,tetpil,i,m,RU, s = signe(n);
1346: GEN x = cgetg(3,t_VEC);
1347: ulong j;
1348:
1349: RU = itos(gmael(nf,2,1)) + itos(gmael(nf,2,2));
1350: x[2] =(long)zerocol(RU); settyp(x[2],t_VEC);
1351: if (absi_cmp(n,stoi(16)) < 0)
1352: {
1353: x[1] = s? (long)idealpowprime(nf,vp,n):
1354: (long)idmat(lgef(nf[1])-3);
1355: tetpil=avma;
1356: return gerepile(av,tetpil,ideallllred(nf,x,NULL,prec));
1357: }
1358:
1359: i = lgefint(n)-1; m=n[i]; j=HIGHBIT;
1360: while ((m&j)==0) j>>=1;
1361: x[1] = (long)prime_to_ideal_aux(nf,vp);
1362: for (j>>=1; j; j>>=1)
1363: {
1364: x = idealmul(nf,x,x);
1365: if (m&j) x[1] = (long)idealmulprime(nf,(GEN)x[1],vp);
1366: x = ideallllred(nf,x,NULL,prec);
1367: }
1368: for (i--; i>=2; i--)
1369: for (m=n[i],j=HIGHBIT; j; j>>=1)
1370: {
1371: x = idealmul(nf,x,x);
1372: if (m&j) x[1] = (long)idealmulprime(nf,(GEN)x[1],vp);
1373: x = ideallllred(nf,x,NULL,prec);
1374: }
1375: if (s < 0) x = idealinv(nf,x);
1376: return gerepileupto(av,x);
1377: }
1378:
1379: long
1380: isideal(GEN nf,GEN x)
1381: {
1382: long N,av,i,j,k,tx=typ(x),lx;
1383: GEN p1,minv;
1384:
1385: nf=checknf(nf); lx=lg(x);
1386: if (tx==t_VEC && lx==3) { x=(GEN)x[1]; tx=typ(x); lx=lg(x); }
1387: if (is_scalar_t(tx))
1388: return (tx==t_INT || tx==t_FRAC || tx==t_FRACN || tx==t_POL ||
1389: (tx==t_POLMOD && gegal((GEN)nf[1],(GEN)x[1])));
1390: if (typ(x)==t_VEC) return (lx==6);
1391: if (typ(x)!=t_MAT) return 0;
1392: if (lx == 1) return 1;
1393: N=lgef(nf[1])-2; if (lg(x[1]) != N) return 0;
1394:
1395: av=avma;
1396: if (lx != N) x = idealmat_to_hnf(nf,x);
1397: x = gdiv(x,content(x)); minv=ginv(x);
1398:
1399: for (i=1; i<N; i++)
1400: for (j=1; j<N; j++)
1401: {
1402: p1=gmul(minv, element_mulid(nf,(GEN)x[i],j));
1403: for (k=1; k<N; k++)
1404: if (typ(p1[k])!=t_INT) { avma=av; return 0; }
1405: }
1406: avma=av; return 1;
1407: }
1408:
1409: GEN
1410: idealdiv(GEN nf, GEN x, GEN y)
1411: {
1412: long av=avma,tetpil;
1413: GEN z=idealinv(nf,y);
1414:
1415: tetpil=avma; return gerepile(av,tetpil,idealmul(nf,x,z));
1416: }
1417:
1418: GEN
1419: idealdivexact(GEN nf, GEN x, GEN y)
1420: /* This routine computes the quotient x/y of two ideals in the number field
1421: * nf. It assumes that the quotient is an integral ideal.
1422: *
1423: * The idea is to find an ideal z dividing y
1424: * such that gcd(N(x)/N(z), N(z)) = 1. Then
1425: *
1426: * x + (N(x)/N(z)) x
1427: * --------------- = -----
1428: * y + (N(y)/N(z)) y
1429: *
1430: * When x and y are integral ideals, this identity can be checked by looking
1431: * at the exponent of a prime ideal p on both sides of the equation.
1432: *
1433: * Specifically, if a prime ideal p divides N(z), then it divides neither
1434: * N(x)/N(z) nor N(y)/N(z) (since N(x)/N(z) is the product of the integers
1435: * N(x/y) and N(y/z)). Both the numerator and the denominator on the left
1436: * will be coprime to p. So will x/y, since x/y is assumed integral and its
1437: * norm N(x/y) is coprime to p
1438: *
1439: * If instead p does not divide N(z), then the power of p dividing N(x)/N(z)
1440: * is the same as the power of p dividing N(x), which is at least as large
1441: * as the power of p dividing x. Likewise for N(y)/N(z). So the power of p
1442: * dividing the left side equals the power of dividing the right side.
1443: *
1444: * Peter Montgomery
1445: * July, 1994.
1446: */
1447: {
1448: long av = avma, tetpil,N;
1449: GEN x1,y1,detx1,dety1,detq,gcancel,gtemp, cy = content(y);
1450:
1451: nf=checknf(nf); N=lgef(nf[1])-3;
1452: if (gcmp0(cy)) err(talker, "cannot invert zero ideal");
1453:
1454: x1 = gdiv(x,cy); detx1 = idealnorm(nf,x1);
1455: if (gcmp0(detx1)) { avma = av; return gcopy(x); } /* numerator is zero */
1456:
1457: y1 = gdiv(y,cy); dety1 = idealnorm(nf,y1);
1458: detq = gdiv(detx1,dety1);
1459: if (!gcmp1(denom(x1)) || typ(detq) != t_INT)
1460: err(talker, "quotient not integral in idealdivexact");
1461: gcancel = dety1;
1462: /* Find a norm gcancel such that
1463: * (1) gcancel divides dety1;
1464: * (2) gcd(detx1/gcancel, gcancel) = 1.
1465: */
1466: do
1467: {
1468: gtemp = ggcd(gcancel, gdiv(detx1,gcancel));
1469: gcancel = gdiv(gcancel,gtemp);
1470: }
1471: while (!gcmp1(gtemp));
1472: /* x1 + (detx1/gcancel)
1473: * Replace x1/y1 by: --------------------
1474: * y1 + (dety1/gcancel)
1475: */
1476:
1477: x1 = idealadd(nf, x1, gscalmat(gdiv(detx1, gcancel), N));
1478: /* y1 reduced to unit ideal ? */
1479: if (gegal(gcancel,dety1)) return gerepileupto(av, x1);
1480:
1481: y1 = idealadd(nf,y1, gscalmat(gdiv(dety1,gcancel), N));
1482: y1 = hnfideal_inv(nf,y1); tetpil = avma;
1483: return gerepile(av, tetpil, idealmat_mul(nf,x1,y1));
1484: }
1485:
1486: GEN
1487: idealintersect(GEN nf, GEN x, GEN y)
1488: {
1489: long av=avma,lz,i,N;
1490: GEN z,dx,dy;
1491:
1492: nf=checknf(nf); N=lgef(nf[1])-3;
1493: if (idealtyp(&x,&z)!=t_MAT || lg(x)!=N+1) x=idealhermite_aux(nf,x);
1494: if (idealtyp(&y,&z)!=t_MAT || lg(y)!=N+1) y=idealhermite_aux(nf,y);
1495: dx=denom(x); if (!gcmp1(dx)) y=gmul(y,dx);
1496: dy=denom(y); if (!gcmp1(dy)) x=gmul(x,dy);
1497: dx = mulii(dx,dy);
1498: z=kerint(concatsp(x,y)); lz=lg(z);
1499: for (i=1; i<lz; i++) setlg(z[i], N+1);
1500: z=gmul(x,z); z = hnfmod(z,detint(z));
1501: if (!gcmp1(dx)) z = gdiv(z,dx);
1502: return gerepileupto(av,z);
1503: }
1504:
1505: static GEN
1506: computet2twist(GEN nf, GEN vdir)
1507: {
1508: long j, ru = lg(nf[6]);
1509: GEN p1,MC, mat = (GEN)nf[5];
1510:
1511: if (!vdir) return (GEN)mat[3];
1512: MC=(GEN)mat[2]; p1=cgetg(ru,t_MAT);
1513: for (j=1; j<ru; j++)
1514: {
1515: GEN v = (GEN)vdir[j];
1516: if (gcmp0(v))
1517: p1[j]=MC[j];
1518: else if (typ(v) == t_INT)
1519: p1[j]=lmul2n((GEN)MC[j],itos(v)<<1);
1520: else
1521: p1[j]=lmul((GEN)MC[j],gpui(stoi(4),v,0));
1522: }
1523: return mulmat_real(p1,(GEN)mat[1]);
1524: }
1525:
1526: GEN
1527: ideallllredall(GEN nf, GEN x, GEN vdir, long prec, long precint)
1528: {
1529: long tx,N,av,tetpil,i,j;
1530: GEN iax,ix,res,ax,p1,p2,y,alpha,beta,pol;
1531:
1532: nf=checknf(nf);
1533: if (vdir)
1534: {
1535: if (gcmp0(vdir)) vdir = NULL;
1536: else if (typ(vdir)!=t_VEC || lg(vdir) != lg(nf[6])) err(idealer5);
1537: }
1538: pol = (GEN)nf[1]; N=lgef(pol)-3;
1539: tx = idealtyp(&x,&ax); ix=x; iax=ax;
1540: if (ax) res = cgetg(3,t_VEC);
1541: av = avma;
1542: if (tx == id_PRINCIPAL)
1543: {
1544: if (ax)
1545: {
1546: res = cgetg(3,t_VEC);
1547: ax = get_arch(nf,x,prec); av=avma;
1548: }
1549: x = idealhermite_aux(nf,x);
1550: }
1551: else if (tx == id_PRIME)
1552: x = prime_to_ideal_aux(nf,x);
1553: else if (lg(x) != N+1) /* id_MAT */
1554: x = idealhermite_aux(nf,x);
1555:
1556: if (DEBUGLEVEL>=6) msgtimer("entering idealllred");
1557: p1=content(x); if (!gcmp1(p1)) x=gdiv(x,p1);
1558:
1559: for (i=1; ; i++)
1560: {
1561: p1=computet2twist(nf,vdir);
1562: if (DEBUGLEVEL>=6) msgtimer("twisted T2");
1563: y = qf_base_change(p1,x,1);
1564: j = 1 + (gexpo(y)>>TWOPOTBITS_IN_LONG);
1565: if (j<0) j=0;
1566: p1 = lllgramintern(y,100,1,j+precint);
1567: if (p1) break;
1568:
1569: if (i == MAXITERPOL) err(accurer,"ideallllredall");
1570: precint=(precint<<1)-2; prec=max(prec,precint);
1571: if (DEBUGLEVEL) err(warnprec,"ideallllredall",precint);
1572: nf=nfnewprec(nf,(j>>1)+precint);
1573: }
1574: y = gmul(x,(GEN)p1[1]);
1575: if (DEBUGLEVEL>=6) msgtimer("lllgram");
1576:
1577: i=2; while (i<=N && gcmp0((GEN)y[i])) i++;
1578: if (i>N)
1579: {
1580: if (x!=ix) x = gerepileupto(av,x); else { avma=av; x = gcopy(x); }
1581: if (!ax) return x;
1582: if (ax==iax) ax = gcopy(ax);
1583: res[1]=(long)x; res[2]=(long)ax; return res;
1584: }
1585: alpha = gmul((GEN)nf[7], y);
1586: /* beta = norm(alpha) / alpha */
1587: beta = gmul(subres(pol,alpha), ginvmod(alpha,pol));
1588: beta = algtobasis_intern(nf,beta);
1589: if (DEBUGLEVEL>=6) msgtimer("alpha/beta");
1590:
1591: p2 = cgetg(N+1,t_MAT);
1592: for (i=1; i<=N; i++)
1593: p2[i] = (long)element_muli(nf,beta,(GEN)x[i]);
1594: p1=content(p2); if (!gcmp1(p1)) p2=gdiv(p2,p1);
1595: if (DEBUGLEVEL>=6) msgtimer("new ideal");
1596: if (ax) y = gclone(gneg_i(get_arch(nf,y,prec)));
1597:
1598: p1 = detint(p2); tetpil = avma;
1599: p2 = gerepile(av,tetpil,hnfmod(p2,p1));
1600: if (DEBUGLEVEL>=6) msgtimer("final hnf");
1601: if (!ax) return p2;
1602: res[1]=(long)p2; res[2]=ladd(ax,y);
1603: gunclone(y); return res;
1604: }
1605:
1606: GEN
1607: ideallllred(GEN nf, GEN ix, GEN vdir, long prec)
1608: {
1609: return ideallllredall(nf,ix,vdir,prec,prec);
1610: }
1611:
1612: GEN
1613: minideal(GEN nf, GEN x, GEN vdir, long prec)
1614: {
1615: long N,av=avma,tetpil,j,RU,tx;
1616: GEN p1,p2,p3,y;
1617:
1618: if (!gcmp0(vdir) && typ(vdir)!=t_VEC) err(idealer5);
1619: nf=checknf(nf); N=lgef(nf[1])-3;
1620: tx = idealtyp(&x,&y); if (tx!=id_MAT) x = idealhermite_aux(nf,x);
1621: RU = N - itos(gmael(nf,2,2)); p1=(GEN)nf[5];
1622: if (gcmp0(vdir)) p1=(GEN)p1[3];
1623: else
1624: {
1625: p3=(GEN)p1[2]; p2=cgetg(RU+1,t_MAT);
1626: for (j=1; j<=RU; j++)
1627: p2[j] = lmul2n((GEN)p3[j], itos((GEN)vdir[j])<<1);
1628: p1=greal(gmul(p2,(GEN)p1[1]));
1629: }
1630: y = gmul(x,(GEN)lllgram(qf_base_change(p1,x,0),prec)[1]);
1631: tetpil = avma; return gerepile(av,tetpil,principalidele(nf,y,prec));
1632: }
1633: static GEN
1634: appr_reduce(GEN s, GEN y, long N)
1635: {
1636: GEN p1,u,z = cgetg(N+2,t_MAT);
1637: long i;
1638:
1639: s=gmod(s,gcoeff(y,1,1)); y=gmul(y,lllint(y));
1640: for (i=1; i<=N; i++) z[i]=y[i]; z[N+1]=(long)s;
1641: u=(GEN)ker(z)[1]; p1 = denom(u);
1642: if (!gcmp1(p1)) u=gmul(u,p1);
1643: p1=(GEN)u[N+1]; setlg(u,N+1);
1644: for (i=1; i<=N; i++) u[i]=lround(gdiv((GEN)u[i],p1));
1645: return gadd(s, gmul(y,u));
1646: }
1647:
1648: /* Given a fractional ideal x (if fl=0) or a prime ideal factorization
1649: * with possibly zero or negative exponents (if fl=1), gives a b such that
1650: * v_p(b)=v_p(x) for all prime ideals p dividing x (or in the factorization)
1651: * and v_p(b)>=0 for all other p, using the (standard) proof given in GTM 138.
1652: * Certainly not the most efficient, but sure.
1653: */
1654: GEN
1655: idealappr0(GEN nf, GEN x, long fl)
1656: {
1657: long av=avma,tetpil,i,j,k,l,N,r,r2;
1658: GEN fact,fact2,list,ep,ep1,ep2,y,z,v,p1,p2,p3,p4,s,pr,alpha,beta,den;
1659:
1660: if (DEBUGLEVEL>4)
1661: {
1662: fprintferr(" entree dans idealappr0() :\n");
1663: fprintferr(" x = "); outerr(x);
1664: }
1665: if (fl)
1666: {
1667: nf=checknf(nf); N=lgef(nf[1])-3;
1668: if (typ(x)!=t_MAT || lg(x)!=3)
1669: err(talker,"not a prime ideal factorization in idealappr0");
1670: fact=x; list=(GEN)fact[1]; ep=(GEN)fact[2]; r=lg(list);
1671: if (r==1) return gscalcol_i(gun,N);
1672: for (i=1; i<r; i++)
1673: if (signe(ep[i])<0)
1674: {
1675: ep1=cgetg(r,t_COL);
1676: for (i=1; i<r; i++)
1677: ep1[i] = (signe(ep[i])>=0)? zero: lnegi((GEN)ep[i]);
1678: fact[2]=(long)ep1; beta=idealappr0(nf,fact,1);
1679: fact2=idealfactor(nf,beta);
1680: p1=(GEN)fact2[1]; r2=lg(p1);
1681: ep2=(GEN)fact2[2]; l=r+r2-1;
1682: z=cgetg(l,t_VEC); for (i=1; i<r; i++) z[i]=list[i];
1683: ep1=cgetg(l,t_VEC);
1684: for (i=1; i<r; i++)
1685: ep1[i] = (signe(ep[i])<=0)? zero: licopy((GEN)ep[i]);
1686: j=r-1;
1687: for (i=1; i<r2; i++)
1688: {
1689: p3=(GEN)p1[i]; k=1;
1690: while (k<r &&
1691: ( !gegal((GEN)p3[1],gmael(list,k,1))
1692: || !element_val(nf,(GEN)p3[2],(GEN)list[k]) )) k++;
1693: if (k==r) { j++; z[j]=(long)p3; ep1[j]=ep2[i]; }
1694: }
1695: fact=cgetg(3,t_MAT);
1696: fact[1]=(long)z; setlg(z,j+1);
1697: fact[2]=(long)ep1; setlg(ep1,j+1);
1698: alpha=idealappr0(nf,fact,1); tetpil=avma;
1699: if (DEBUGLEVEL>2)
1700: {
1701: fprintferr(" alpha = "); outerr(alpha);
1702: fprintferr(" beta = "); outerr(beta);
1703: }
1704: return gerepile(av,tetpil,element_div(nf,alpha,beta));
1705: }
1706: y=idmat(N);
1707: for (i=1; i<r; i++)
1708: {
1709: pr=(GEN)list[i];
1710: if (signe(ep[i]))
1711: {
1712: p4=addsi(1,(GEN)ep[i]); p1=powgi((GEN)pr[1],p4);
1713: if (cmpis((GEN)pr[4],N))
1714: {
1715: p2=cgetg(3,t_MAT);
1716: p2[1]=(long)gscalcol_i(p1, N);
1717: p2[2]=(long)element_pow(nf,(GEN)pr[2],p4);
1718: y=idealmat_mul(nf,y,p2);
1719: }
1720: else y=gmul(p1,y);
1721: }
1722: else y=idealmulprime(nf,y,pr);
1723: }
1724: }
1725: else
1726: {
1727: den=denom(x); if (gcmp1(den)) den=NULL; else x=gmul(den,x);
1728: x=idealhermite_aux(nf,x); N=lgef(nf[1])-3;
1729: fact=idealfactor(nf,x);
1730: list=(GEN)fact[1]; ep=(GEN)fact[2]; r=lg(list);
1731: if (r==1) { avma=av; return gscalcol_i(gun,N); }
1732: if (den)
1733: {
1734: fact2=idealfactor(nf,den);
1735: p1=(GEN)fact2[1]; r2=lg(p1);
1736: l=r+r2-1;
1737: z=cgetg(l,t_COL); for (i=1; i<r; i++) z[i]=list[i];
1738: ep1=cgetg(l,t_COL); for (i=1; i<r; i++) ep1[i]=ep[i];
1739: j=r-1;
1740: for (i=1; i<r2; i++)
1741: {
1742: p3=(GEN)p1[i]; k=1;
1743: while (k<r && !gegal((GEN)list[k],p3)) k++;
1744: if (k==r){ j++; z[j]=(long)p3; ep1[j]=zero; }
1745: }
1746: fact=cgetg(3,t_MAT);
1747: fact[1]=(long)z; setlg(z,j+1);
1748: fact[2]=(long)ep1; setlg(ep1,j+1);
1749: alpha=idealappr0(nf,fact,1);
1750: if (DEBUGLEVEL>2) { fprintferr(" alpha = "); outerr(alpha); }
1751: tetpil=avma; return gerepile(av,tetpil,gdiv(alpha,den));
1752: }
1753: y=x; for (i=1; i<r; i++) y=idealmulprime(nf,y,(GEN)list[i]);
1754: }
1755:
1756: z=cgetg(r,t_VEC);
1757: for (i=1; i<r; i++)
1758: {
1759: pr=(GEN)list[i]; p4=addsi(1, (GEN)ep[i]); p1=powgi((GEN)pr[1],p4);
1760: if (cmpis((GEN)pr[4],N))
1761: {
1762: p2=cgetg(3,t_MAT);
1763: p2[1]=(long)gscalcol_i(p1,N);
1764: p2[2]=(long)element_pow(nf,(GEN)pr[5],p4);
1765: z[i]=ldiv(idealmat_mul(nf,y,p2),p1);
1766: }
1767: else z[i]=ldiv(y,p1);
1768: }
1769: v=idealaddmultoone(nf,z);
1770: s=cgetg(N+1,t_COL); for (i=1; i<=N; i++) s[i]=zero;
1771: for (i=1; i<r; i++)
1772: {
1773: pr=(GEN)list[i];
1774: if (signe(ep[i]))
1775: s=gadd(s,element_mul(nf,(GEN)v[i],element_pow(nf,(GEN)pr[2],(GEN)ep[i])));
1776: else s=gadd(s,(GEN)v[i]);
1777: }
1778: p3 = appr_reduce(s,y,N);
1779: if (DEBUGLEVEL>2)
1780: { fprintferr(" sortie de idealappr0 p3 = "); outerr(p3); }
1781: return gerepileupto(av,p3);
1782: }
1783:
1784: /* Given a prime ideal factorization x with possibly zero or negative exponents,
1785: * and a vector y of elements of nf, gives a b such that
1786: * v_p(b-y_p)>=v_p(x) for all prime ideals p in the ideal factorization
1787: * and v_p(b)>=0 for all other p, using the (standard) proof given in GTM 138.
1788: * Certainly not the most efficient, but sure.
1789: */
1790: GEN
1791: idealchinese(GEN nf, GEN x, GEN y)
1792: {
1793: long ty=typ(y),av=avma,i,j,k,l,N,r,r2;
1794: GEN fact,fact2,list,ep,ep1,ep2,z,t,v,p1,p2,p3,p4,s,pr,den;
1795:
1796: if (DEBUGLEVEL>4)
1797: {
1798: fprintferr(" entree dans idealchinese() :\n");
1799: fprintferr(" x = "); outerr(x);
1800: fprintferr(" y = "); outerr(y);
1801: }
1802: nf=checknf(nf); N=lgef(nf[1])-3;
1803: if (typ(x)!=t_MAT ||(lg(x)!=3))
1804: err(talker,"not a prime ideal factorization in idealchinese");
1805: fact=x; list=(GEN)fact[1]; ep=(GEN)fact[2]; r=lg(list);
1806: if (!is_vec_t(ty) || lg(y)!=r)
1807: err(talker,"not a suitable vector of elements in idealchinese");
1808: if (r==1) return gscalcol_i(gun,N);
1809:
1810: den=denom(y);
1811: if (!gcmp1(den))
1812: {
1813: fact2=idealfactor(nf,den);
1814: p1=(GEN)fact2[1]; r2=lg(p1);
1815: ep2=(GEN)fact2[2]; l=r+r2-1;
1816: z=cgetg(l,t_VEC); for (i=1; i<r; i++) z[i]=list[i];
1817: ep1=cgetg(l,t_VEC); for (i=1; i<r; i++) ep1[i]=ep[i];
1818: j=r-1;
1819: for (i=1; i<r2; i++)
1820: {
1821: p3=(GEN)p1[i]; k=1;
1822: while (k<r && !gegal((GEN)list[k],p3)) k++;
1823: if (k==r) { j++; z[j]=(long)p3; ep1[j]=ep2[i]; }
1824: else ep1[k]=ladd((GEN)ep1[k],(GEN)ep2[i]);
1825: }
1826: r=j+1; setlg(z,r); setlg(ep1,r); list=z; ep=ep1;
1827: }
1828: for (i=1; i<r; i++)
1829: if (signe(ep[i])<0) ep[i] = zero;
1830: t=idmat(N);
1831: for (i=1; i<r; i++)
1832: {
1833: pr=(GEN)list[i]; p4=(GEN)ep[i];
1834: if (signe(p4))
1835: {
1836: if (cmpis((GEN)pr[4],N))
1837: {
1838: p2=cgetg(3,t_MAT);
1839: p2[1]=(long)gscalcol_i(gpui((GEN)pr[1],p4,0), N);
1840: p2[2]=(long)element_pow(nf,(GEN)pr[2],p4);
1841: t=idealmat_mul(nf,t,p2);
1842: }
1843: else t=gmul(gpui((GEN)pr[1],p4,0),t);
1844: }
1845: }
1846: z=cgetg(r,t_VEC);
1847: for (i=1; i<r; i++)
1848: {
1849: pr=(GEN)list[i]; p4=(GEN)ep[i];
1850: if (cmpis((GEN)pr[4],N))
1851: {
1852: p2=cgetg(3,t_MAT); p1=gpui((GEN)pr[1],p4,0);
1853: p2[1]=(long)gscalcol_i(p1,N);
1854: p2[2]=(long)element_pow(nf,(GEN)pr[5],p4);
1855: z[i]=ldiv(idealmat_mul(nf,t,p2),p1);
1856: }
1857: else z[i]=ldiv(t,gpui((GEN)pr[1],p4,0));
1858: }
1859: v=idealaddmultoone(nf,z);
1860: s=cgetg(N+1,t_COL); for (i=1; i<=N; i++) s[i]=zero;
1861: for (i=1; i<r; i++)
1862: s = gadd(s,element_mul(nf,(GEN)v[i],(GEN)y[i]));
1863:
1864: p3 = appr_reduce(s,t,N);
1865: if (DEBUGLEVEL>2)
1866: { fprintferr(" sortie de idealchinese() : p3 = "); outerr(p3); }
1867: return gerepileupto(av,p3);
1868: }
1869:
1870: GEN
1871: idealappr(GEN nf, GEN x) { return idealappr0(nf,x,0); }
1872:
1873: GEN
1874: idealapprfact(GEN nf, GEN x) { return idealappr0(nf,x,1); }
1875:
1876: /* Given an integral ideal x and a in x, gives a b such that
1877: * x=aZ_K+bZ_K using a different algorithm than ideal_two_elt
1878: */
1879: GEN
1880: ideal_two_elt2(GEN nf, GEN x, GEN a)
1881: {
1882: long ta=typ(a), av=avma,tetpil,i,r;
1883: GEN con,ep,b,list,fact;
1884:
1885: nf = checknf(nf);
1886: if (!is_extscalar_t(ta) && typ(a)!=t_COL)
1887: err(typeer,"ideal_two_elt2");
1888: x = idealhermite_aux(nf,x);
1889: if (gcmp0(x))
1890: {
1891: if (!gcmp0(a)) err(talker,"element not in ideal in ideal_two_elt2");
1892: avma=av; return gcopy(a);
1893: }
1894: con = content(x);
1895: if (gcmp1(con)) con = NULL; else { x = gdiv(x,con); a = gdiv(a,con); }
1896: a = principalideal(nf,a);
1897: if (!gcmp1(denom(gauss(x,a))))
1898: err(talker,"element does not belong to ideal in ideal_two_elt2");
1899:
1900: fact=idealfactor(nf,a); list=(GEN)fact[1];
1901: r=lg(list); ep = (GEN)fact[2];
1902: for (i=1; i<r; i++) ep[i]=lstoi(idealval(nf,x,(GEN)list[i]));
1903: b = centermod(idealappr0(nf,fact,1), gcoeff(x,1,1));
1904: tetpil=avma; b = con? gmul(b,con): gcopy(b);
1905: return gerepile(av,tetpil,b);
1906: }
1907:
1908: /* Given 2 integral ideals x and y in a number field nf gives a beta
1909: * belonging to nf such that beta.x is an integral ideal coprime to y
1910: */
1911: GEN
1912: idealcoprime(GEN nf, GEN x, GEN y)
1913: {
1914: long av=avma,tetpil,i,r;
1915: GEN fact,list,p2,ep;
1916:
1917: if (DEBUGLEVEL>4)
1918: {
1919: fprintferr(" entree dans idealcoprime() :\n");
1920: fprintferr(" x = "); outerr(x);
1921: fprintferr(" y = "); outerr(y);
1922: }
1923: fact=idealfactor(nf,y); list=(GEN)fact[1];
1924: r=lg(list); ep = (GEN)fact[2];
1925: for (i=1; i<r; i++) ep[i]=lstoi(-idealval(nf,x,(GEN)list[i]));
1926: tetpil=avma; p2=idealappr0(nf,fact,1);
1927: if (DEBUGLEVEL>4)
1928: { fprintferr(" sortie de idealcoprime() : p2 = "); outerr(p2); }
1929: return gerepile(av,tetpil,p2);
1930: }
1931:
1932: /* returns the first index i<=n such that x=v[i] if it exits, 0 otherwise */
1933: long
1934: isinvector(GEN v, GEN x, long n)
1935: {
1936: long i;
1937:
1938: for (i=1; i<=n; i++)
1939: if (gegal((GEN)v[i],x)) return i;
1940: return 0;
1941: }
1942:
1943: /* Given an integral ideal x and three algebraic integers a, b and c in a
1944: * number field nf gives a beta belonging to nf such that beta.x^(-1) is an
1945: * integral ideal coprime to abc.Z_k
1946: */
1947: static GEN
1948: idealcoprimeinvabc(GEN nf, GEN x, GEN a, GEN b, GEN c)
1949: {
1950: long av=avma,tetpil,i,j,r,ra,rb,rc;
1951: GEN facta,factb,factc,fact,factall,p1,p2,ep;
1952:
1953: if (DEBUGLEVEL>4)
1954: {
1955: fprintferr(" entree dans idealcoprimeinvabc() :\n");
1956: fprintferr(" x = "); outerr(x); fprintferr(" a = "); outerr(a);
1957: fprintferr(" b = "); outerr(b); fprintferr(" c = "); outerr(c);
1958: flusherr();
1959: }
1960: facta=(GEN)idealfactor(nf,a)[1]; factb=(GEN)idealfactor(nf,b)[1];
1961: factc=(GEN)idealfactor(nf,c)[1]; ra=lg(facta); rb=lg(factb); rc=lg(factc);
1962: factall=cgetg(ra+rb+rc-2,t_COL);
1963: for (i=1; i<ra; i++) factall[i]=facta[i]; j=ra-1;
1964: for (i=1; i<rb; i++)
1965: if (!isinvector(factall,(GEN)factb[i],j)) factall[++j]=factb[i];
1966: for (i=1; i<rc; i++)
1967: if (!isinvector(factall,(GEN)factc[i],j)) factall[++j]=factc[i];
1968: r=j+1; fact=cgetg(3,t_MAT); p1=cgetg(r,t_COL); ep=cgetg(r,t_COL);
1969: for (i=1; i<r; i++) p1[i]=factall[i];
1970: for (i=1; i<r; i++) ep[i]=lstoi(idealval(nf,x,(GEN)p1[i]));
1971: fact[1]=(long)p1; fact[2]=(long)ep; tetpil=avma; p2=idealappr0(nf,fact,1);
1972: if (DEBUGLEVEL>2)
1973: { fprintferr(" sortie de idealcoprimeinvabc() : p2 = "); outerr(p2); }
1974: return gerepile(av,tetpil,p2);
1975: }
1976:
1977: /* Solve the equation ((b+aX).Z_k/((a,b).J),M)=Z_k. */
1978: static GEN
1979: findX(GEN nf, GEN a, GEN b, GEN J, GEN M)
1980: {
1981: long N,i,k,r,v;
1982: GEN p1,p2,abJ,fact,list,ve,ep,int0,int1,int2,pr;
1983:
1984: if (DEBUGLEVEL>4)
1985: {
1986: fprintferr(" entree dans findX() :\n");
1987: fprintferr(" a = "); outerr(a); fprintferr(" b = "); outerr(b);
1988: fprintferr(" J = "); outerr(J); fprintferr(" M = "); outerr(M);
1989: }
1990: N=lgef(nf[1])-3;
1991: p1=cgetg(3,t_MAT); p1[1]=(long)a; p1[2]=(long)b;
1992: if (N==2) p1=idealmul(nf,p1,idmat(2));
1993: abJ=idealmul(nf,p1,J);
1994: fact=idealfactor(nf,M); list=(GEN)fact[1]; r=lg(list);
1995: ve=cgetg(r,t_VEC); ep=cgetg(r,t_VEC);
1996: int0=cgetg(N+1,t_COL); int1=cgetg(N+1,t_COL); int2=cgetg(N+1,t_COL);
1997: for (i=2; i<=N; i++) int0[i]=int1[i]=int2[i]=zero;
1998: int0[1]=zero; int1[1]=un; int2[1]=deux;
1999: for (i=1; i<r; i++)
2000: {
2001: pr=(GEN)list[i]; v=element_val(nf,a,pr);
2002: if (v)
2003: {
2004: ep[i]=un;
2005: ve[i] = (element_val(nf,b,pr)<=v)? (long)int0: (long)int1;
2006: }
2007: else
2008: {
2009: v=idealval(nf,abJ,pr);
2010: p1=element_div(nf,idealaddtoone_i(nf,a,pr),a);
2011: ep[i]=lstoi(v+1); k=1;
2012: while (k<=v)
2013: {
2014: p1=element_mul(nf,p1,gsub(int2,element_mul(nf,a,p1)));
2015: k<<=1;
2016: }
2017: p1=element_mul(nf,p1,gsub(element_pow(nf,(GEN)pr[2],stoi(v)),b));
2018: ve[i]=lmod(p1,gpuigs((GEN)pr[1],v+1));
2019: }
2020: }
2021: fact[2]=(long)ep; p2=idealchinese(nf,fact,ve);
2022: if (DEBUGLEVEL>2) { fprintferr(" sortie de findX() : p2 = "); outerr(p2); }
2023: return p2;
2024: }
2025:
2026: /* A usage interne. Given a, b, c, d in nf, gives an algebraic integer y in
2027: * nf such that [c,d]-y.[a,b] is "small"
2028: */
2029: static GEN
2030: nfreducemat(GEN nf, GEN a, GEN b, GEN c, GEN d)
2031: {
2032: long av=avma,tetpil,i,i1,i2,j,j1,j2,k,N;
2033: GEN p1,p2,X,M,y,mult,s;
2034:
2035: mult=(GEN)nf[9]; N=lgef(nf[1])-3; X=cgetg(N+1,t_COL);
2036: for (j=1; j<=N; j++)
2037: {
2038: s=gzero;
2039: for (i=1; i<=N; i++) for (k=1; k<=N; k++)
2040: {
2041: p1=gcoeff(mult,k,j+(i-1)*N);
2042: if (!gcmp0(p1))
2043: s=gadd(s,gmul(p1,gadd(gmul((GEN)a[i],(GEN)c[k]),
2044: gmul((GEN)b[i],(GEN)d[k]))));
2045: }
2046: X[j]=(long)s;
2047: }
2048: M=cgetg(N+1,t_MAT);
2049: for (j2=1; j2<=N; j2++)
2050: {
2051: p1=cgetg(N+1,t_COL); M[j2]=(long)p1;
2052: for (j1=1; j1<=N; j1++)
2053: {
2054: s=gzero;
2055: for (i1=1; i1<=N; i1++)
2056: for (i2=1; i2<=N; i2++)
2057: for (k=1; k<=N; k++)
2058: {
2059: p2=gmul(gcoeff(mult,k,j1+(i1-1)*N),gcoeff(mult,k,j2+(i2-1)*N));
2060: if (!gcmp0(p2))
2061: s=gadd(s,gmul(p2,gadd(gmul((GEN)a[i1],(GEN)a[i2]),
2062: gmul((GEN)b[i1],(GEN)b[i2]))));
2063: }
2064: p1[j1]=(long)s;
2065: }
2066: }
2067: y=gauss(M,X); tetpil=avma;
2068: return gerepile(av,tetpil,ground(y));
2069: }
2070:
2071: /* Given 3 algebraic integers a,b,c in a number field nf given by their
2072: * components on the integral basis, gives a three-component vector [d,e,U]
2073: * whose first two components are algebraic integers d,e and the third a
2074: * unimodular 3x3-matrix U such that [a,b,c]*U=[0,d,e]
2075: */
2076: GEN
2077: threetotwo2(GEN nf, GEN a, GEN b, GEN c)
2078: {
2079: long av=avma,tetpil,i,N;
2080: GEN y,p1,p2,p3,M,X,Y,J,e,b1,c1,u,v,U,int0,Z,pk;
2081:
2082: if (DEBUGLEVEL>2)
2083: {
2084: fprintferr(" On entre dans threetotwo2() : \n");
2085: fprintferr(" a = "); outerr(a);
2086: fprintferr(" b = "); outerr(b);
2087: fprintferr(" c = "); outerr(c);
2088: }
2089: if (gcmp0(a))
2090: {
2091: y=cgetg(4,t_VEC); y[1]=lcopy(b); y[2]=lcopy(c); y[3]=(long)idmat(3);
2092: return y;
2093: }
2094: if (gcmp0(b))
2095: {
2096: y=cgetg(4,t_VEC); y[1]=lcopy(a); y[2]=lcopy(c);
2097: e = idmat(3); i=e[1]; e[1]=e[2]; e[2]=i;
2098: y[3]=(long)e; return y;
2099: }
2100: if (gcmp0(c))
2101: {
2102: y=cgetg(4,t_VEC); y[1]=lcopy(a); y[2]=lcopy(b);
2103: e = idmat(3); i=e[1]; e[1]=e[3]; e[3]=e[2]; e[2]=i;
2104: y[3]=(long)e; return y;
2105: }
2106:
2107: N=lgef(nf[1])-3;
2108: p1=cgetg(4,t_MAT); p1[1]=(long)a; p1[2]=(long)b;
2109: p1[3]=(long)c; p1=idealhermite_aux(nf,p1);
2110: if (DEBUGLEVEL>2)
2111: { fprintferr(" ideal a.Z_k+b.Z_k+c.Z_k = "); outerr(p1); }
2112: J=idealdiv(nf,e=idealcoprimeinvabc(nf,p1,a,b,c),p1);
2113: if (DEBUGLEVEL>2)
2114: { fprintferr(" ideal J = "); outerr(J); fprintferr(" e = "); outerr(e); }
2115: p1=cgetg(3,t_MAT); p1[1]=(long)a; p1[2]=(long)b; M=idealmul(nf,p1,J);
2116: if (DEBUGLEVEL>2)
2117: { fprintferr(" ideal M=(a.Z_k+b.Z_k).J = "); outerr(M); }
2118: X=findX(nf,a,b,J,M);
2119: if (DEBUGLEVEL>2){ fprintferr(" X = "); outerr(X); }
2120: p1=gadd(b,element_mul(nf,a,X));
2121: p2=cgetg(3,t_MAT); p2[1]=(long)element_mul(nf,a,p1);
2122: p2[2]=(long)element_mul(nf,c,p1);
2123: if (N==2) p2=idealhermite_aux(nf,p2);
2124: p3=cgetg(3,t_MAT); p3[1]=(long)a; p3[2]=(long)b;
2125: p3=idealhermite_aux(nf,p3);
2126: if (DEBUGLEVEL>2)
2127: { fprintferr(" ideal a.Z_k+b.Z_k = "); outerr(p3); }
2128: Y=findX(nf,a,c,J,idealdiv(nf,p2,p3));
2129: if (DEBUGLEVEL>2) { fprintferr(" Y = "); outerr(Y); }
2130: b1=element_div(nf,p1,e);
2131: if (DEBUGLEVEL>2) { fprintferr(" b1 = "); outerr(b1); }
2132: p2=gadd(c,element_mul(nf,a,Y));
2133: c1=element_div(nf,p2,e);
2134: if (DEBUGLEVEL>2) { fprintferr(" c1 = "); outerr(c1); }
2135: p1=idealhermite_aux(nf,b1);
2136: p2=idealhermite_aux(nf,c1);
2137: p3=idealaddtoone(nf,p1,p2);
2138: u=element_div(nf,(GEN)p3[1],b1); v=element_div(nf,(GEN)p3[2],c1);
2139: if (DEBUGLEVEL>2)
2140: { fprintferr(" u = "); outerr(u); fprintferr(" v = "); outerr(v); }
2141: U=cgetg(4,t_MAT);
2142: p1=cgetg(4,t_COL); p2=cgetg(4,t_COL); p3=cgetg(4,t_COL);
2143: U[1]=(long)p1; U[2]=(long)p2; U[3]=(long)p3;
2144: p1[1]=lsub(element_mul(nf,X,c1),element_mul(nf,Y,b1));
2145: p1[2]=(long)c1; p1[3]=lneg(b1);
2146: int0 = zerocol(N);
2147: p2[1]=(long)gscalcol_i(gun,N);
2148: p2[2]=p2[3]=(long)int0;
2149: Z=gadd(element_mul(nf,X,u),element_mul(nf,Y,v));
2150: pk=nfreducemat(nf,c1,(GEN)p1[3],u,v);
2151: p3[1]=(long)int0; p3[2]=lsub(u,element_mul(nf,pk,c1));
2152: p3[3]=(long)gadd(v,element_mul(nf,pk,b1));
2153: e=gadd(e,element_mul(nf,a,gsub(element_mul(nf,pk,(GEN)p1[1]),Z)));
2154: tetpil=avma;
2155: y=cgetg(4,t_VEC); y[1]=lcopy(a); y[2]=lcopy(e); y[3]=lcopy(U);
2156: if (DEBUGLEVEL>2)
2157: { fprintferr(" sortie de threetotwo2() : y = "); outerr(y); }
2158: return gerepile(av,tetpil,y);
2159: }
2160:
2161: /* Given 3 algebraic integers a,b,c in a number field nf given by their
2162: * components on the integral basis, gives a three-component vector [d,e,U]
2163: * whose first two components are algebraic integers d,e and the third a
2164: * unimodular 3x3-matrix U such that [a,b,c]*U=[0,d,e] Naive method which may
2165: * not work, but fast and small coefficients.
2166: */
2167: GEN
2168: threetotwo(GEN nf, GEN a, GEN b, GEN c)
2169: {
2170: long av=avma,tetpil,i,N;
2171: GEN pol,p1,p2,p3,p4,y,id,hu,h,V,U,r,vd,q1,q1a,q2,q2a,na,nb,nc,nr;
2172:
2173: nf=checknf(nf); pol=(GEN)nf[1]; N=lgef(pol)-3; id=idmat(N);
2174: na=gnorml2(a); nb=gnorml2(b); nc=gnorml2(c);
2175: U=gmul(idmat(3),gmodulsg(1,pol));
2176: if (gcmp(nc,nb)<0)
2177: {
2178: p1=c; c=b; b=p1; p1=nc; nc=nb; nb=p1;
2179: p1=(GEN)U[3]; U[3]=U[2]; U[2]=(long)p1;
2180: }
2181: if (gcmp(nc,na)<0)
2182: {
2183: p1=a; a=c; c=p1; p1=na; na=nc; nc=p1;
2184: p1=(GEN)U[1]; U[1]=U[3]; U[3]=(long)p1;
2185: }
2186: while (!gcmp0(gmin(na,nb)))
2187: {
2188: p1=cgetg(2*N+1,t_MAT);
2189: for (i=1; i<=N; i++)
2190: {
2191: p1[i]=(long)element_mul(nf,a,(GEN)id[i]);
2192: p1[i+N]=(long)element_mul(nf,b,(GEN)id[i]);
2193: }
2194: hu=hnfall(p1); h=(GEN)hu[1]; V=(GEN)hu[2];
2195: p2=(GEN)ker(concatsp(h,c))[1]; p3=(GEN)p2[N+1];
2196: p4=cgetg(N+1,t_COL);
2197: for (i=1; i<=N; i++) p4[i]=(long)ground(gdiv((GEN)p2[i],p3));
2198: r=gadd(c,gmul(h,p4));
2199: vd=cgetg(N+1,t_MAT); for (i=1; i<=N; i++) vd[i]=V[N+i];
2200: p2=gmul(vd,p4);
2201: q1=cgetg(N+1,t_COL); q2=cgetg(N+1,t_COL);
2202: for (i=1; i<=N; i++) { q1[i]=p2[i]; q2[i]=p2[i+N]; }
2203: q1a=basistoalg(nf,q1); q2a=basistoalg(nf,q2);
2204: U[3]=(long)gadd((GEN)U[3],gadd(gmul(q1a,(GEN)U[1]),gmul(q2a,(GEN)U[2])));
2205: nr=gnorml2(r);
2206: if (gcmp(nr,gmax(na,nb))>=0) err(talker,"threetotwo does not work");
2207: if (gcmp(na,nb)>=0)
2208: {
2209: c=a; nc=na; a=r; na=nr; p1=(GEN)U[1]; U[1]=U[3]; U[3]=(long)p1;
2210: }
2211: else
2212: {
2213: c=b; nc=nb; b=r; nb=nr; p1=(GEN)U[2]; U[2]=U[3]; U[3]=(long)p1;
2214: }
2215: }
2216: if (!gcmp0(na))
2217: {
2218: p1=a; a=b; b=p1; p1=(GEN)U[1]; U[1]=U[2]; U[2]=(long)p1;
2219: }
2220: tetpil=avma; y=cgetg(4,t_VEC); y[1]=lcopy(b); y[2]=lcopy(c);
2221: y[3]=(long)algtobasis(nf,U); return gerepile(av,tetpil,y);
2222: }
2223:
2224: /* Given 2 algebraic integers a,b in a number field nf given by their
2225: * components on the integral basis, gives a three-components vector [d,e,U]
2226: * whose first two component are algebraic integers d,e and the third a
2227: * unimodular 2x2-matrix U such that [a,b]*U=[d,e], with d and e hopefully
2228: * smaller than a and b.
2229: */
2230: GEN
2231: twototwo(GEN nf, GEN a, GEN b)
2232: {
2233: long av=avma,tetpil;
2234: GEN pol,p1,y,U,r,qr,qa,na,nb,nr;
2235:
2236: nf=checknf(nf);
2237: pol=(GEN)nf[1];
2238: na=gnorml2(a); nb=gnorml2(b);
2239: U=gmul(idmat(2),gmodulsg(1,pol));
2240: if (gcmp(na,nb)>0)
2241: {
2242: p1=a; a=b; b=p1; p1=na; na=nb; nb=p1;
2243: p1=(GEN)U[2]; U[2]=U[1]; U[1]=(long)p1;
2244: }
2245:
2246: while (!gcmp0(na))
2247: {
2248: qr=nfdivres(nf,b,a); r=(GEN)qr[2]; nr=gnorml2(r);
2249: if (gcmp(nr,na)<0)
2250: {
2251: b=a; a=r; nb=na; na=nr; qa=basistoalg(nf,(GEN)qr[1]);
2252: p1=gsub((GEN)U[2],gmul(qa,(GEN)U[1])); U[2]=U[1]; U[1]=(long)p1;
2253: }
2254: else
2255: {
2256: if (gcmp(nr,nb)<0)
2257: {
2258: qa=basistoalg(nf,(GEN)qr[1]);
2259: U[2]=lsub((GEN)U[2],gmul(qa,(GEN)U[1]));
2260: }
2261: break;
2262: }
2263: }
2264: tetpil=avma; y=cgetg(4,t_VEC); y[1]=lcopy(a); y[2]=lcopy(b);
2265: y[3]=(long)algtobasis(nf,U); return gerepile(av,tetpil,y);
2266: }
2267:
2268: GEN
2269: elt_mul_get_table(GEN nf, GEN x)
2270: {
2271: long i,lx = lg(x);
2272: GEN mul=cgetg(lx,t_MAT);
2273:
2274: /* assume w_1 = 1 */
2275: mul[1]=(long)x;
2276: for (i=2; i<lx; i++) mul[i] = (long)element_mulid(nf,x,i);
2277: return mul;
2278: }
2279:
2280: GEN
2281: elt_mul_table(GEN mul, GEN z)
2282: {
2283: long av = avma, i, lx = lg(mul);
2284: GEN p1 = gmul((GEN)z[1], (GEN)mul[1]);
2285:
2286: for (i=2; i<lx; i++)
2287: if (!gcmp0((GEN)z[i])) p1 = gadd(p1, gmul((GEN)z[i], (GEN)mul[i]));
2288: return gerepileupto(av, p1);
2289: }
2290:
2291: GEN
2292: element_mulvec(GEN nf, GEN x, GEN v)
2293: {
2294: long lv=lg(v),i;
2295: GEN mul = elt_mul_get_table(nf,x), y=cgetg(lv,t_COL);
2296:
2297: for (i=1; i<lv; i++)
2298: y[i] = (long)elt_mul_table(mul,(GEN)v[i]);
2299: return y;
2300: }
2301:
2302: static GEN
2303: element_mulvecrow(GEN nf, GEN x, GEN m, long i, long lim)
2304: {
2305: long lv,j;
2306: GEN y, mul = elt_mul_get_table(nf,x);
2307:
2308: lv=min(lg(m),lim+1); y=cgetg(lv,t_VEC);
2309: for (j=1; j<lv; j++)
2310: y[j] = (long)elt_mul_table(mul,gcoeff(m,i,j));
2311: return y;
2312: }
2313:
2314: /* Given an element x and an ideal in matrix form (not necessarily HNF),
2315: * gives an r such that x-r is in ideal and r is small. No checks
2316: */
2317: GEN
2318: element_reduce(GEN nf, GEN x, GEN ideal)
2319: {
2320: long tx=typ(x),av=avma,tetpil,N,i;
2321: GEN p1,u;
2322:
2323: if (is_extscalar_t(tx))
2324: x = algtobasis_intern(checknf(nf), x);
2325: N = lg(x); p1=cgetg(N+1,t_MAT);
2326: for (i=1; i<N; i++) p1[i]=ideal[i];
2327: p1[N]=(long)x; u=(GEN)ker(p1)[1];
2328: p1=(GEN)u[N]; setlg(u,N);
2329: for (i=1; i<N; i++) u[i]=lround(gdiv((GEN)u[i],p1));
2330: u=gmul(ideal,u); tetpil=avma;
2331: return gerepile(av,tetpil,gadd(u,x));
2332: }
2333:
2334: /* A torsion-free module M over Z_K will be given by a row vector [A,I] with
2335: * two components. I=[\a_1,...,\a_k] is a row vector of k fractional ideals
2336: * given in HNF. A is an nxk matrix (same k and n the rank of the module)
2337: * such that if A_j is the j-th column of A then M=\a_1A_1+...\a_kA_k. We say
2338: * that [A,I] is a pseudo-basis if k=n
2339: */
2340:
2341: /* Given a torsion-free module x as above outputs a pseudo-basis for x in
2342: * Hermite Normal Form
2343: */
2344: GEN
2345: nfhermite(GEN nf, GEN x)
2346: {
2347: long av=avma,tetpil,i,j,def,k,m;
2348: GEN p1,p2,y,A,I,J;
2349:
2350: nf=checknf(nf);
2351: if (typ(x)!=t_VEC || lg(x)!=3) err(talker,"not a module in nfhermite");
2352: A=(GEN)x[1]; I=(GEN)x[2]; k=lg(A)-1;
2353: if (typ(A)!=t_MAT) err(talker,"not a matrix in nfhermite");
2354: if (typ(I)!=t_VEC || lg(I)!=k+1)
2355: err(talker,"not a correct ideal list in nfhermite");
2356: if (!k) err(talker,"not a matrix of maximal rank in nfhermite");
2357: m=lg(A[1])-1;
2358: if (k<m) err(talker,"not a matrix of maximal rank in nfhermite");
2359:
2360: p1 = cgetg(k+1,t_MAT); for (j=1; j<=k; j++) p1[j]=A[j];
2361: A = p1; I = dummycopy(I);
2362: for (j=1; j<=k; j++)
2363: if (typ(I[j])!=t_MAT) I[j]=(long)idealhermite_aux(nf,(GEN)I[j]);
2364:
2365: J=cgetg(k+1,t_VEC); def=k+1;
2366: for (i=m; i>=1; i--)
2367: {
2368: GEN den,p4,p5,p6,u,v,newid, invnewid = NULL;
2369:
2370: def--; j=def; while (j>=1 && gcmp0(gcoeff(A,i,j))) j--;
2371: if (!j) err(talker,"not a matrix of maximal rank in nfhermite");
2372: if (j==def) j--;
2373: else
2374: {
2375: p1=(GEN)A[j]; A[j]=A[def]; A[def]=(long)p1;
2376: p1=(GEN)I[j]; I[j]=I[def]; I[def]=(long)p1;
2377: }
2378: p1=gcoeff(A,i,def); p2=element_inv(nf,p1);
2379: A[def]=(long)element_mulvec(nf,p2,(GEN)A[def]);
2380: I[def]=(long)idealmul(nf,p1,(GEN)I[def]);
2381: for ( ; j; j--)
2382: {
2383: p1=gcoeff(A,i,j);
2384: if (!gcmp0(p1))
2385: {
2386: p2=idealmul(nf,p1,(GEN)I[j]);
2387: newid = idealadd(nf,p2,(GEN)I[def]);
2388: invnewid = hnfideal_inv(nf,newid);
2389: p4 = idealmul(nf, p2, invnewid);
2390: p5 = idealmul(nf,(GEN)I[def],invnewid);
2391: y = idealaddtoone(nf,p4,p5);
2392: u=element_div(nf,(GEN)y[1],p1); v=(GEN)y[2];
2393: p6=gsub((GEN)A[j],element_mulvec(nf,p1,(GEN)A[def]));
2394: A[def]=ladd(element_mulvec(nf,u,(GEN)A[j]),
2395: element_mulvec(nf,v,(GEN)A[def]));
2396: A[j]=(long)p6;
2397: I[j]=(long)idealmul(nf,idealmul(nf,(GEN)I[j],(GEN)I[def]),invnewid);
2398: I[def]=(long)newid; den=denom((GEN)I[j]);
2399: if (!gcmp1(den))
2400: {
2401: I[j]=lmul(den,(GEN)I[j]);
2402: A[j]=ldiv((GEN)A[j],den);
2403: }
2404: }
2405: }
2406: if (!invnewid) invnewid = hnfideal_inv(nf,(GEN)I[def]);
2407: p1=(GEN)I[def]; J[def]=(long)invnewid;
2408: for (j=def+1; j<=k; j++)
2409: {
2410: p2=gsub(element_reduce(nf,gcoeff(A,i,j),idealmul(nf,p1,(GEN)J[j])),
2411: gcoeff(A,i,j));
2412: A[j]=ladd((GEN)A[j],element_mulvec(nf,p2,(GEN)A[def]));
2413: }
2414: }
2415: tetpil=avma; y=cgetg(3,t_VEC);
2416: p1=cgetg(m+1,t_MAT); y[1]=(long)p1;
2417: p2=cgetg(m+1,t_VEC); y[2]=(long)p2;
2418: for (j=1; j<=m; j++) p1[j]=lcopy((GEN)A[j+k-m]);
2419: for (j=1; j<=m; j++) p2[j]=lcopy((GEN)I[j+k-m]);
2420: return gerepile(av,tetpil,y);
2421: }
2422:
2423: /* A torsion module M over Z_K will be given by a row vector [A,I,J] with
2424: * three components. I=[b_1,...,b_n] is a row vector of k fractional ideals
2425: * given in HNF, J=[a_1,...,a_n] is a row vector of n fractional ideals in
2426: * HNF. A is an nxn matrix (same n) such that if A_j is the j-th column of A
2427: * and e_n is the canonical basis of K^n, then
2428: * M=(b_1e_1+...+b_ne_n)/(a_1A_1+...a_nA_n)
2429: */
2430:
2431: /* We input a torsion module x=[A,I,J] as above, and output the
2432: * smith normal form as K=[\c_1,...,\c_n] such that x=Z_K/\c_1+...+Z_K/\c_n.
2433: */
2434: GEN
2435: nfsmith(GEN nf, GEN x)
2436: {
2437: long av,tetpil,i,j,k,l,lim,c,fl,n,m,N;
2438: GEN p1,p2,p3,p4,z,b,u,v,w,d,dinv,unnf,A,I,J;
2439:
2440: nf=checknf(nf); N=lgef(nf[1])-3;
2441: if (typ(x)!=t_VEC || lg(x)!=4) err(talker,"not a module in nfsmith");
2442: A=(GEN)x[1]; I=(GEN)x[2]; J=(GEN)x[3];
2443: if (typ(A)!=t_MAT) err(talker,"not a matrix in nfsmith");
2444: n=lg(A)-1;
2445: if (typ(I)!=t_VEC || lg(I)!=n+1 || typ(J)!=t_VEC || lg(J)!=n+1)
2446: err(talker,"not a correct ideal list in nfsmith");
2447: if (!n) err(talker,"not a matrix of maximal rank in nfsmith");
2448: m=lg(A[1])-1;
2449: if (n<m) err(talker,"not a matrix of maximal rank in nfsmith");
2450: if (n>m) err(impl,"nfsmith for non square matrices");
2451:
2452: av=avma; lim=stack_lim(av,1);
2453: p1 = cgetg(n+1,t_MAT); for (j=1; j<=n; j++) p1[j]=A[j];
2454: A = p1; I = dummycopy(I); J=dummycopy(J);
2455: for (j=1; j<=n; j++)
2456: if (typ(I[j])!=t_MAT) I[j]=(long)idealhermite_aux(nf,(GEN)I[j]);
2457: for (j=1; j<=n; j++)
2458: if (typ(J[j])!=t_MAT) J[j]=(long)idealhermite_aux(nf,(GEN)J[j]);
2459: for (i=n; i>=2; i--)
2460: {
2461: do
2462: {
2463: c=0;
2464: for (j=i-1; j>=1; j--)
2465: {
2466: p1=gcoeff(A,i,j);
2467: if (!gcmp0(p1))
2468: {
2469: p2=gcoeff(A,i,i);
2470: d=nfbezout(nf,p2,p1,(GEN)J[i],(GEN)J[j],&u,&v,&w,&dinv);
2471: if (!gcmp0(u))
2472: {
2473: if (!gcmp0(v))
2474: b=gadd(element_mulvec(nf,u,(GEN)A[i]),
2475: element_mulvec(nf,v,(GEN)A[j]));
2476: else b=element_mulvec(nf,u,(GEN)A[i]);
2477: }
2478: else b=element_mulvec(nf,v,(GEN)A[j]);
2479: A[j]=lsub(element_mulvec(nf,p2,(GEN)A[j]),
2480: element_mulvec(nf,p1,(GEN)A[i]));
2481: A[i]=(long)b; J[j]=(long)w; J[i]=(long)d;
2482: }
2483: }
2484: for (j=i-1; j>=1; j--)
2485: {
2486: p1=gcoeff(A,j,i);
2487: if (!gcmp0(p1))
2488: {
2489: p2=gcoeff(A,i,i);
2490: d=nfbezout(nf,p2,p1,(GEN)I[i],(GEN)I[j],&u,&v,&w,&dinv);
2491: if (gcmp0(u))
2492: b=element_mulvecrow(nf,v,A,j,i);
2493: else
2494: {
2495: if (gcmp0(v))
2496: b=element_mulvecrow(nf,u,A,i,i);
2497: else
2498: b=gadd(element_mulvecrow(nf,u,A,i,i),
2499: element_mulvecrow(nf,v,A,j,i));
2500: }
2501: p3=gsub(element_mulvecrow(nf,p2,A,j,i),
2502: element_mulvecrow(nf,p1,A,i,i));
2503: for (k=1; k<=i; k++) { coeff(A,j,k)=p3[k]; coeff(A,i,k)=b[k]; }
2504: I[j]=(long)w; I[i]=(long)d; c++;
2505: }
2506: }
2507: if (!c)
2508: {
2509: b=gcoeff(A,i,i); if (gcmp0(b)) break;
2510:
2511: b=idealmul(nf,b,idealmul(nf,(GEN)J[i],(GEN)I[i]));
2512: fl=1;
2513: for (k=1; k<i && fl; k++)
2514: for (l=1; l<i && fl; l++)
2515: {
2516: p3=gcoeff(A,k,l);
2517: if (!gcmp0(p3))
2518: fl=gegal(idealadd(nf,b,idealmul(nf,p3,idealmul(nf,(GEN)J[l],(GEN)I[k]))),b);
2519: }
2520: if (!fl)
2521: {
2522: k--; l--;
2523: b=idealdiv(nf,(GEN)I[k],(GEN)I[i]);
2524: p4=gauss(idealdiv(nf,(GEN)J[i],idealmul(nf,p3,(GEN)J[l])),b);
2525: l=1; while (l<=N && gcmp1(denom((GEN)p4[l]))) l++;
2526: if (l>N) err(talker,"bug2 in nfsmith");
2527: p3=element_mulvecrow(nf,(GEN)b[l],A,k,i);
2528: for (l=1; l<=i; l++)
2529: coeff(A,i,l) = ladd(gcoeff(A,i,l),(GEN)p3[l]);
2530: }
2531: }
2532: if (low_stack(lim, stack_lim(av,1)))
2533: {
2534: GEN *gptr[3];
2535: if(DEBUGMEM>1) err(warnmem,"nfsmith");
2536: gptr[0]=&A; gptr[1]=&I; gptr[2]=&J; gerepilemany(av,gptr,3);
2537: }
2538: }
2539: while (c || !fl);
2540: }
2541: unnf=gscalcol_i(gun,N);
2542: p1=gcoeff(A,1,1); coeff(A,1,1)=(long)unnf;
2543: J[1]=(long)idealmul(nf,p1,(GEN)J[1]);
2544: for (i=2; i<=n; i++)
2545: if (!gegal(gcoeff(A,i,i),unnf)) err(talker,"bug in nfsmith");
2546: tetpil=avma; z=cgetg(n+1,t_VEC);
2547: for (i=1; i<=n; i++) z[i]=(long)idealmul(nf,(GEN)I[i],(GEN)J[i]);
2548: return gerepile(av,tetpil,z);
2549: }
2550:
2551: /*******************************************************************/
2552: /* */
2553: /* ALGEBRE LINEAIRE DANS LES CORPS DE NOMBRES */
2554: /* */
2555: /*******************************************************************/
2556:
2557: #define trivlift(x) ((typ(x)==t_POLMOD)? (GEN)x[2]: lift_intern(x))
2558:
2559: GEN
2560: element_mulmodpr2(GEN nf, GEN x, GEN y, GEN prhall)
2561: {
2562: long av=avma;
2563: GEN p1;
2564:
2565: nf=checknf(nf); checkprhall(prhall);
2566: p1 = element_mul(nf,x,y);
2567: return gerepileupto(av,nfreducemodpr(nf,p1,prhall));
2568: }
2569:
2570: /* On ne peut PAS definir ca comme les autres par
2571: * #define element_divmodpr() nfreducemodpr(element_div())
2572: * car le element_div ne marche pas en general
2573: */
2574: GEN
2575: element_divmodpr(GEN nf, GEN x, GEN y, GEN prhall)
2576: {
2577: long av=avma;
2578: GEN p1;
2579:
2580: nf=checknf(nf); checkprhall(prhall);
2581: p1=lift_intern(gdiv(gmodulcp(gmul((GEN)nf[7],trivlift(x)), (GEN)nf[1]),
2582: gmodulcp(gmul((GEN)nf[7],trivlift(y)), (GEN)nf[1])));
2583: p1=algtobasis_intern(nf,p1);
2584: return gerepileupto(av,nfreducemodpr(nf,p1,prhall));
2585: }
2586:
2587: GEN
2588: element_invmodpr(GEN nf, GEN y, GEN prhall)
2589: {
2590: long av=avma;
2591: GEN p1;
2592:
2593: p1=ginvmod(gmul((GEN)nf[7],trivlift(y)), (GEN)nf[1]);
2594: p1=algtobasis_intern(nf,p1);
2595: return gerepileupto(av,nfreducemodpr(nf,p1,prhall));
2596: }
2597:
2598: GEN
2599: element_powmodpr(GEN nf,GEN x,GEN k,GEN prhall)
2600: {
2601: long av=avma,N,s;
2602: GEN y,z;
2603:
2604: nf=checknf(nf); checkprhall(prhall);
2605: N=lgef(nf[1])-3;
2606: s=signe(k); k=(s>=0)?k:negi(k);
2607: z=x; y = gscalcol_i(gun,N);
2608: for(;;)
2609: {
2610: if (mpodd(k)) y=element_mulmodpr(nf,z,y,prhall);
2611: k=shifti(k,-1);
2612: if (signe(k)) z=element_sqrmodpr(nf,z,prhall);
2613: else
2614: {
2615: cgiv(k); if (s<0) y = element_invmodpr(nf,y,prhall);
2616: return gerepileupto(av,y);
2617: }
2618: }
2619: }
2620:
2621: /* x est une matrice dont les coefficients sont des vecteurs dans la base
2622: * d'entiers modulo un ideal premier prhall, sous forme reduite modulo prhall.
2623: */
2624: GEN
2625: nfkermodpr(GEN nf, GEN x, GEN prhall)
2626: {
2627: long i,j,k,r,t,n,m,av0,av,av1,av2,N,lim;
2628: GEN c,d,y,unnf,munnf,zeromodp,zeronf,p,pp,prh;
2629:
2630: nf=checknf(nf); checkprhall(prhall);
2631: if (typ(x)!=t_MAT) err(typeer,"nfkermodpr");
2632: n=lg(x)-1; if (!n) return cgetg(1,t_MAT);
2633: prh=(GEN)prhall[1]; av0=avma;
2634: N=lgef(nf[1])-3; pp=gcoeff(prh,1,1);
2635:
2636: zeromodp=gmodulsg(0,pp);
2637: unnf=cgetg(N+1,t_COL); unnf[1]=(long)gmodulsg(1,pp);
2638: zeronf=cgetg(N+1,t_COL); zeronf[1]=(long)zeromodp;
2639:
2640: av=avma; munnf=cgetg(N+1,t_COL); munnf[1]=(long)gmodulsg(-1,pp);
2641: for (i=2; i<=N; i++)
2642: zeronf[i] = munnf[i] = unnf[i] = (long)zeromodp;
2643:
2644: m=lg(x[1])-1; x=dummycopy(x); r=0;
2645: c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0;
2646: d=new_chunk(n+1); av1=avma; lim=stack_lim(av1,1);
2647: for (k=1; k<=n; k++)
2648: {
2649: j=1;
2650: while (j<=m && (c[j] || gcmp0(gcoeff(x,j,k)))) j++;
2651: if (j>m) { r++; d[k]=0; }
2652: else
2653: {
2654: p=element_divmodpr(nf,munnf,gcoeff(x,j,k),prhall);
2655: c[j]=k; d[k]=j; coeff(x,j,k)=(long)munnf;
2656: for (i=k+1; i<=n; i++)
2657: coeff(x,j,i)=(long)element_mulmodpr(nf,p,gcoeff(x,j,i),prhall);
2658: for (t=1; t<=m; t++)
2659: if (t!=j)
2660: {
2661: p=gcoeff(x,t,k); coeff(x,t,k)=(long)zeronf;
2662: for (i=k+1; i<=n; i++)
2663: coeff(x,t,i)=ladd(gcoeff(x,t,i),
2664: element_mulmodpr(nf,p,gcoeff(x,j,i),prhall));
2665: }
2666: if (low_stack(lim, stack_lim(av1,1)))
2667: {
2668: if (DEBUGMEM>1) err(warnmem,"nfkermodpr, k = %ld / %ld",k,n);
2669: av2=avma; x=gerepile(av1,av2,gcopy(x));
2670: }
2671: }
2672: }
2673: if (!r) { avma=av0; return cgetg(1,t_MAT); }
2674: av1=avma; y=cgetg(r+1,t_MAT);
2675: for (j=k=1; j<=r; j++,k++)
2676: {
2677: p=cgetg(n+1,t_COL); y[j]=(long)p; while (d[k]) k++;
2678: for (i=1; i<k; i++) p[i]=d[i]? lcopy(gcoeff(x,d[i],k)): (long)zeronf;
2679: p[k]=(long)unnf; for (i=k+1; i<=n; i++) p[i]=(long)zeronf;
2680: }
2681: return gerepile(av,av1,y);
2682: }
2683:
2684: /* a.x=b ou b est un vecteur */
2685: GEN
2686: nfsolvemodpr(GEN nf, GEN a, GEN b, GEN prhall)
2687: {
2688: long nbli,nbco,i,j,k,av=avma,tetpil;
2689: GEN aa,x,p,m,u;
2690:
2691: nf=checknf(nf); checkprhall(prhall);
2692: if (typ(a)!=t_MAT) err(typeer,"nfsolvemodpr");
2693: nbco=lg(a)-1; nbli=lg(a[1])-1;
2694: if (typ(b)!=t_COL) err(typeer,"nfsolvemodpr");
2695: if (lg(b)!=nbco+1) err(mattype1,"nfsolvemodpr");
2696: x=cgetg(nbli+1,t_COL);
2697: for (j=1; j<=nbco; j++) x[j]=b[j];
2698: aa=cgetg(nbco+1,t_MAT);
2699: for (j=1; j<=nbco; j++)
2700: {
2701: aa[j]=lgetg(nbli+1,t_COL);
2702: for (i=1; i<=nbli; i++) coeff(aa,i,j)=coeff(a,i,j);
2703: }
2704: for (i=1; i<nbli; i++)
2705: {
2706: p=gcoeff(aa,i,i); k=i;
2707: if (gcmp0(p))
2708: {
2709: k=i+1; while (k<=nbli && gcmp0(gcoeff(aa,k,i))) k++;
2710: if (k>nbco) err(matinv1);
2711: for (j=i; j<=nbco; j++)
2712: {
2713: u=gcoeff(aa,i,j); coeff(aa,i,j)=coeff(aa,k,j);
2714: coeff(aa,k,j)=(long)u;
2715: }
2716: u=(GEN)x[i]; x[i]=x[k]; x[k]=(long)u;
2717: p=gcoeff(aa,i,i);
2718: }
2719: for (k=i+1; k<=nbli; k++)
2720: {
2721: m=gcoeff(aa,k,i);
2722: if (!gcmp0(m))
2723: {
2724: m=element_divmodpr(nf,m,p,prhall);
2725: for (j=i+1; j<=nbco; j++)
2726: coeff(aa,k,j)=lsub(gcoeff(aa,k,j),
2727: element_mulmodpr(nf,m,gcoeff(aa,i,j),prhall));
2728: x[k]=lsub((GEN)x[k],element_mulmodpr(nf,m,(GEN)x[i],prhall));
2729: }
2730: }
2731: }
2732: /* Resolution systeme triangularise */
2733: p=gcoeff(aa,nbli,nbco); if (gcmp0(p)) err(matinv1);
2734:
2735: x[nbli]=(long)element_divmodpr(nf,(GEN)x[nbli],p,prhall);
2736: for (i=nbli-1; i>0; i--)
2737: {
2738: m=(GEN)x[i];
2739: for (j=i+1; j<=nbco; j++)
2740: m=gsub(m,element_mulmodpr(nf,gcoeff(aa,i,j),(GEN)x[j],prhall));
2741: x[i]=(long)element_divmodpr(nf,m,gcoeff(aa,i,i),prhall);
2742: }
2743: tetpil=avma; return gerepile(av,tetpil,gcopy(x));
2744: }
2745:
2746: GEN
2747: nfsuppl(GEN nf, GEN x, long n, GEN prhall)
2748: {
2749: long av=avma,av2,k,s,t,N,lx=lg(x);
2750: GEN y,p1,p2,p,unmodp,zeromodp,unnf,zeronf,prh;
2751: stackzone *zone;
2752:
2753: k=lx-1; if (k>n) err(suppler2);
2754: if (k && lg(x[1])!=n+1) err(talker,"incorrect dimension in nfsupl");
2755: N=lgef(nf[1])-3; prh=(GEN)prhall[1]; p=gcoeff(prh,1,1);
2756:
2757: zone = switch_stack(NULL, 2*(3 + 2*lg(p) + N+1) + (n+3)*(n+1));
2758: switch_stack(zone,1);
2759: unmodp=gmodulsg(1,p); zeromodp=gmodulsg(0,p);
2760: unnf=gscalcol_proto(unmodp,zeromodp,N);
2761: zeronf=gscalcol_proto(zeromodp,zeromodp,N);
2762: y = idmat_intern(n,unnf,zeronf);
2763: switch_stack(zone,0); av2=avma;
2764:
2765: for (s=1; s<=k; s++)
2766: {
2767: p1=nfsolvemodpr(nf,y,(GEN)x[s],prhall); t=s;
2768: while (t<=n && gcmp0((GEN)p1[t])) t++;
2769: avma=av2; if (t>n) err(suppler2);
2770: p2=(GEN)y[s]; y[s]=x[s]; if (s!=t) y[t]=(long)p2;
2771: }
2772: avma=av; y=gcopy(y);
2773: free(zone); return y;
2774: }
2775:
2776: /* Given two fractional ideals a and b, gives x in a, y in b, z in b^-1,
2777: t in a^-1 such that xt-yz=1. In the present version, z is in Z. */
2778: GEN
2779: nfidealdet1(GEN nf, GEN a, GEN b)
2780: {
2781: long av=avma,tetpil;
2782: GEN x,p1,p2,res,z,da,db;
2783:
2784: da=denom(a); if (gcmp1(da)) da = NULL; else a=gmul(da,a);
2785: db=denom(b); if (gcmp1(db)) db = NULL; else b=gmul(db,b);
2786: a = idealinv(nf,a); x=idealcoprime(nf,a,b);
2787: p1=idealmul(nf,x,a); p2=idealaddtoone(nf,p1,b);
2788:
2789: tetpil=avma; res=cgetg(5,t_VEC);
2790: res[1] = da? ldiv(x,da): lcopy(x);
2791: res[2] = db? ldiv((GEN)p2[2],db): lcopy((GEN)p2[2]);
2792: z = db? gneg_i(db): negi(gun);
2793: res[3] = (long) gscalcol_i(z, lgef(nf[1])-3);
2794: res[4] = (long) element_div(nf,(GEN)p2[1],(GEN)res[1]);
2795: return gerepile(av,tetpil,res);
2796: }
2797:
2798: /* Given a pseudo basis pseudo, outputs a multiple of its ideal determinant */
2799: GEN
2800: nfdetint(GEN nf,GEN pseudo)
2801: {
2802: GEN pass,c,v,det1,piv,pivprec,vi,p1,x,I,unnf,zeronf,id,idprod;
2803: long i,j,k,rg,n,n1,m,m1,av=avma,av1,tetpil,lim,cm=0,N;
2804:
2805: nf=checknf(nf); N=lgef(nf[1])-3;
2806: if (typ(pseudo)!=t_VEC || lg(pseudo)!=3)
2807: err(talker,"not a module in nfdetint");
2808: x=(GEN)pseudo[1]; I=(GEN)pseudo[2];
2809: if (typ(x)!=t_MAT) err(talker,"not a matrix in nfdetint");
2810: n1=lg(x); n=n1-1; if (!n) return gun;
2811:
2812: m1=lg(x[1]); m=m1-1;
2813: if (typ(I)!=t_VEC || lg(I)!=n1)
2814: err(talker,"not a correct ideal list in nfdetint");
2815:
2816: unnf=gscalcol_i(gun,N); zeronf=zerocol(N);
2817: id=idmat(N); c=new_chunk(m1); for (k=1; k<=m; k++) c[k]=0;
2818: piv = pivprec = unnf;
2819:
2820: av1=avma; lim=stack_lim(av1,1);
2821: det1=idprod=gzero; /* dummy for gerepilemany */
2822: pass=cgetg(m1,t_MAT); v=cgetg(m1,t_COL);
2823: for (j=1; j<=m; j++)
2824: {
2825: v[j] = zero; /* dummy */
2826: p1=cgetg(m1,t_COL); pass[j]=(long)p1;
2827: for (i=1; i<=m; i++) p1[i]=(long)zeronf;
2828: }
2829: for (rg=0,k=1; k<=n; k++)
2830: {
2831: long t = 0;
2832: for (i=1; i<=m; i++)
2833: if (!c[i])
2834: {
2835: vi=element_mul(nf,piv,gcoeff(x,i,k));
2836: for (j=1; j<=m; j++)
2837: if (c[j]) vi=gadd(vi,element_mul(nf,gcoeff(pass,i,j),gcoeff(x,j,k)));
2838: v[i]=(long)vi; if (!t && !gcmp0(vi)) t=i;
2839: }
2840: if (t)
2841: {
2842: pivprec = piv;
2843: if (rg == m-1)
2844: {
2845: if (!cm)
2846: {
2847: cm=1; idprod = id;
2848: for (i=1; i<=m; i++)
2849: if (i!=t)
2850: idprod = (idprod==id)? (GEN)I[c[i]]
2851: : idealmul(nf,idprod,(GEN)I[c[i]]);
2852: }
2853: p1 = idealmul(nf,(GEN)v[t],(GEN)I[k]); c[t]=0;
2854: det1 = (typ(det1)==t_INT)? p1: idealadd(nf,p1,det1);
2855: }
2856: else
2857: {
2858: rg++; piv=(GEN)v[t]; c[t]=k;
2859: for (i=1; i<=m; i++)
2860: if (!c[i])
2861: {
2862: for (j=1; j<=m; j++)
2863: if (c[j] && j!=t)
2864: {
2865: p1=gsub(element_mul(nf,piv,gcoeff(pass,i,j)),
2866: element_mul(nf,(GEN)v[i],gcoeff(pass,t,j)));
2867: coeff(pass,i,j) = rg>1? (long) element_div(nf,p1,pivprec)
2868: : (long) p1;
2869: }
2870: coeff(pass,i,t)=lneg((GEN)v[i]);
2871: }
2872: }
2873: }
2874: if (low_stack(lim, stack_lim(av1,1)))
2875: {
2876: GEN *gptr[6];
2877: if(DEBUGMEM>1) err(warnmem,"nfdetint");
2878: gptr[0]=&det1; gptr[1]=ϖ gptr[2]=&pivprec; gptr[3]=&pass;
2879: gptr[4]=&v; gptr[5]=&idprod; gerepilemany(av1,gptr,6);
2880: }
2881: }
2882: if (!cm) { avma=av; return gscalmat(gzero,N); }
2883: tetpil=avma; return gerepile(av,tetpil,idealmul(nf,idprod,det1));
2884: }
2885:
2886: /* clean in place (destroy x) */
2887: static void
2888: nfcleanmod(GEN nf, GEN x, long lim, GEN detmat)
2889: {
2890: long lx=lg(x),i;
2891:
2892: if (lim<=0 || lim>=lx) lim=lx-1;
2893: for (i=1; i<=lim; i++)
2894: x[i]=(long)element_reduce(nf,(GEN)x[i],detmat);
2895: }
2896:
2897: static GEN
2898: zero_nfbezout(GEN nf,GEN b, GEN ida,GEN idb,GEN *u,GEN *v,GEN *w,GEN *di)
2899: {
2900: long av, tetpil, j, N=lgef(nf[1])-3;
2901: GEN pab,d;
2902:
2903: d=idealmulelt(nf,b,idb); *di=idealinv(nf,d);
2904: av=avma; pab=idealmul(nf,ida,idb); tetpil=avma;
2905: *w=gerepile(av,tetpil, idealmul(nf,pab,*di));
2906:
2907: *u=cgetg(N+1,t_COL); for (j=1; j<=N; j++) (*u)[j]=zero;
2908: *v=element_inv(nf,b); return d;
2909: }
2910:
2911: /* a usage interne
2912: * Given elements a,b, ideals ida, idb, outputs d=a.ida+b.idb and gives
2913: * di=d^-1, w=ida.idb.di, u, v such that au+bv=1 and u in ida.di, v in
2914: * idb.di. We assume ida, idb non-zero, but a and b can be zero. Error if a
2915: * and b are both zero.
2916: */
2917: static GEN
2918: nfbezout(GEN nf,GEN a,GEN b, GEN ida,GEN idb, GEN *u,GEN *v,GEN *w,GEN *di)
2919: {
2920: GEN pab,pu,pv,pw,uv,d,dinv,pa,pb,pa1,pb1, *gptr[5];
2921: long av,tetpil;
2922:
2923: if (gcmp0(a))
2924: {
2925: if (gcmp0(b)) err(talker,"both elements zero in nfbezout");
2926: return zero_nfbezout(nf,b,ida,idb,u,v,w,di);
2927: }
2928: if (gcmp0(b))
2929: return zero_nfbezout(nf,a,idb,ida,v,u,w,di);
2930:
2931: av = avma;
2932: pa=idealmulelt(nf,a,ida);
2933: pb=idealmulelt(nf,b,idb);
2934:
2935: d=idealadd(nf,pa,pb); dinv=idealinv(nf,d);
2936: pa1=idealmullll(nf,pa,dinv);
2937: pb1=idealmullll(nf,pb,dinv);
2938: uv=idealaddtoone(nf,pa1,pb1);
2939: pab=idealmul(nf,ida,idb); tetpil=avma;
2940:
2941: pu=element_div(nf,(GEN)uv[1],a);
2942: pv=element_div(nf,(GEN)uv[2],b);
2943: d=gcopy(d); dinv=gcopy(dinv);
2944: pw=idealmul(nf,pab,dinv);
2945:
2946: *u=pu; *v=pv; *w=pw; *di=dinv;
2947: gptr[0]=u; gptr[1]=v; gptr[2]=w; gptr[3]=di;
2948: gptr[4]=&d; gerepilemanysp(av,tetpil,gptr,5);
2949: return d;
2950: }
2951:
2952: /* A usage interne. Pas de verifs ni gestion de pile */
2953: GEN
2954: idealoplll(GEN op(GEN,GEN,GEN), GEN nf, GEN x, GEN y)
2955: {
2956: GEN z = op(nf,x,y), den = denom(z);
2957:
2958: if (gcmp1(den)) den = NULL; else z=gmul(den,z);
2959: z=gmul(z,lllintpartial(z));
2960: return den? gdiv(z,den): z;
2961: }
2962:
2963: /* A usage interne. Pas de verifs ni gestion de pile */
2964: GEN
2965: idealmulelt(GEN nf, GEN elt, GEN x)
2966: {
2967: long lx=lg(x),j;
2968: GEN z=cgetg(lx,t_MAT);
2969: for (j=1; j<lx; j++) z[j]=(long)element_mul(nf,elt,(GEN)x[j]);
2970: return z;
2971: }
2972:
2973: GEN
2974: nfhermitemod(GEN nf, GEN pseudo, GEN detmat)
2975: {
2976: long av0=avma,li,co,av,tetpil,i,j,jm1,def,ldef,lim,N;
2977: GEN b,q,w,p1,p2,d,u,v,den,x,I,J,dinv,unnf,wh;
2978:
2979: nf=checknf(nf); N=lgef(nf[1])-3;
2980: if (typ(pseudo)!=t_VEC || lg(pseudo)!=3)
2981: err(talker,"not a module in nfhermitemod");
2982: x=(GEN)pseudo[1]; I=(GEN)pseudo[2];
2983: if (typ(x)!=t_MAT) err(talker,"not a matrix in nfhermitemod");
2984: co=lg(x);
2985: if (typ(I)!=t_VEC || lg(I)!=co)
2986: err(talker,"not a correct ideal list in nfhermitemod");
2987: if (co==1) return cgetg(1,t_MAT);
2988:
2989: li=lg(x[1]); x=dummycopy(x); I=dummycopy(I);
2990: unnf=gscalcol_i(gun,N);
2991: for (j=1; j<co; j++)
2992: if (typ(I[j])!=t_MAT) I[j]=(long)idealhermite_aux(nf,(GEN)I[j]);
2993:
2994: den=denom(detmat); if (!gcmp1(den)) detmat=gmul(den,detmat);
2995: detmat=gmul(detmat,lllintpartial(detmat));
2996:
2997: av=avma; lim=stack_lim(av,1);
2998: def=co; ldef=(li>co)?li-co+1:1;
2999: for (i=li-1; i>=ldef; i--)
3000: {
3001: def--; j=def-1; while (j && gcmp0(gcoeff(x,i,j))) j--;
3002: while (j)
3003: {
3004: jm1=j-1; if (!jm1) jm1=def;
3005: d=nfbezout(nf,gcoeff(x,i,j),gcoeff(x,i,jm1),(GEN)I[j],(GEN)I[jm1],
3006: &u,&v,&w,&dinv);
3007: if (!gcmp0(u))
3008: {
3009: p1=element_mulvec(nf,u,(GEN)x[j]);
3010: if (!gcmp0(v)) p1=gadd(p1, element_mulvec(nf,v,(GEN)x[jm1]));
3011: }
3012: else p1=element_mulvec(nf,v,(GEN)x[jm1]);
3013: x[j]=lsub(element_mulvec(nf,gcoeff(x,i,j),(GEN)x[jm1]),
3014: element_mulvec(nf,gcoeff(x,i,jm1),(GEN)x[j]));
3015: nfcleanmod(nf,(GEN)x[j],i,idealdivlll(nf,detmat,w));
3016: nfcleanmod(nf,p1,i,idealmullll(nf,detmat,dinv));
3017: x[jm1]=(long)p1; I[j]=(long)w; I[jm1]=(long)d;
3018: j--; while (j && gcmp0(gcoeff(x,i,j))) j--;
3019: }
3020: if (low_stack(lim, stack_lim(av,1)))
3021: {
3022: GEN *gptr[2];
3023: if(DEBUGMEM>1) err(warnmem,"[1]: nfhermitemod");
3024: gptr[0]=&x; gptr[1]=&I; gerepilemany(av,gptr,2);
3025: }
3026: }
3027: b=detmat; wh=cgetg(li,t_MAT); def--;
3028: for (i=li-1; i>=1; i--)
3029: {
3030: d = nfbezout(nf,gcoeff(x,i,i+def),unnf,(GEN)I[i+def],b,&u,&v,&w,&dinv);
3031: p1 = element_mulvec(nf,u,(GEN)x[i+def]);
3032: nfcleanmod(nf,p1,i,idealmullll(nf,b,dinv));
3033: wh[i]=(long)p1; coeff(wh,i,i)=(long)unnf; I[i+def]=(long)d;
3034: if (i>1) b=idealmul(nf,b,dinv);
3035: }
3036: J=cgetg(li,t_VEC); J[1]=zero;
3037: for (j=2; j<li; j++) J[j]=(long)idealinv(nf,(GEN)I[j+def]);
3038: for (i=li-2; i>=1; i--)
3039: {
3040: for (j=i+1; j<li; j++)
3041: {
3042: q=idealmul(nf,(GEN)I[i+def],(GEN)J[j]);
3043: p1=gsub(element_reduce(nf,gcoeff(wh,i,j),q),gcoeff(wh,i,j));
3044: wh[j]=(long)gadd((GEN)wh[j],element_mulvec(nf,p1,(GEN)wh[i]));
3045: }
3046: if (low_stack(lim, stack_lim(av,1)))
3047: {
3048: GEN *gptr[3];
3049: if(DEBUGMEM>1) err(warnmem,"[2]: nfhermitemod");
3050: gptr[0]=&wh; gptr[1]=&I; gptr[2]=&J; gerepilemany(av,gptr,3);
3051: }
3052: }
3053: tetpil=avma; p1=cgetg(3,t_VEC); p1[1]=lcopy(wh);
3054: p2=cgetg(li,t_VEC); p1[2]=(long)p2;
3055: for (j=1; j<li; j++) p2[j]=lcopy((GEN)I[j+def]);
3056: return gerepile(av0,tetpil,p1);
3057: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to base4.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