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