Annotation of OpenXM_contrib/pari/src/basemath/base3.c, Revision 1.1.1.1
1.1 maekawa 1: /*******************************************************************/
2: /* */
3: /* BASIC NF OPERATIONS */
4: /* */
5: /*******************************************************************/
6: /* $Id: base3.c,v 1.1.1.1 1999/09/16 13:47:20 karim Exp $ */
7: #include "pari.h"
8: int ker_trivial_mod_p(GEN x, GEN p);
9: long ideal_is_zk(GEN ideal,long N);
10: GEN idealaddtoone_i(GEN nf, GEN x, GEN y);
11:
12: /*******************************************************************/
13: /* */
14: /* OPERATIONS OVER NUMBER FIELD ELEMENTS. */
15: /* These are always represented as column vectors over the */
16: /* integral basis nf[7] */
17: /* */
18: /*******************************************************************/
19:
20: int
21: isnfscalar(GEN x)
22: {
23: long lx=lg(x),i;
24:
25: for (i=2; i<lx; i++)
26: if (!gcmp0((GEN)x[i])) return 0;
27: return 1;
28: }
29:
30: static GEN
31: checknfelt_mod(GEN nf, GEN x)
32: {
33: if (gegal((GEN)x[1],(GEN)nf[1])) return (GEN) x[2];
34: err(talker,"not the same polynomial in number field operation");
35: return NULL; /* not reached */
36: }
37:
38: static GEN
39: scal_mul(GEN nf, GEN x, GEN y, long ty)
40: {
41: long av=avma, tetpil;
42: GEN p1;
43:
44: if (!is_extscalar_t(ty))
45: {
46: if (ty!=t_COL) err(typeer,"nfmul");
47: y = gmul((GEN)nf[7],y);
48: }
49: p1 = gmul(x,y); tetpil=avma;
50: return gerepile(av,tetpil,algtobasis(nf,p1));
51: }
52:
53: /* product of x and y in nf */
54: GEN
55: element_mul(GEN nf, GEN x, GEN y)
56: {
57: long av,i,j,k,N,tx,ty;
58: GEN s,v,c,p1,tab;
59:
60: if (x == y) return element_sqr(nf,x);
61:
62: tx=typ(x); ty=typ(y);
63: nf=checknf(nf); tab = (GEN)nf[9]; N=lgef(nf[1])-3;
64: if (tx==t_POLMOD) x=checknfelt_mod(nf,x);
65: if (ty==t_POLMOD) y=checknfelt_mod(nf,y);
66: if (is_extscalar_t(tx)) return scal_mul(nf,x,y,ty);
67: if (is_extscalar_t(ty)) return scal_mul(nf,y,x,tx);
68: if (isnfscalar(x)) return gmul((GEN)x[1],y);
69: if (isnfscalar(y)) return gmul((GEN)y[1],x);
70:
71: v=cgetg(N+1,t_COL); av=avma;
72: for (k=1; k<=N; k++)
73: {
74: if (k == 1)
75: s = gmul((GEN)x[1],(GEN)y[1]);
76: else
77: s = gadd(gmul((GEN)x[1],(GEN)y[k]),
78: gmul((GEN)x[k],(GEN)y[1]));
79: for (i=2; i<=N; i++)
80: {
81: c=gcoeff(tab,k,(i-1)*N+i);
82: if (signe(c))
83: {
84: p1 = gmul((GEN)x[i],(GEN)y[i]);
85: if (!gcmp1(c)) p1 = gmul(p1,c);
86: s = gadd(s, p1);
87: }
88: for (j=i+1; j<=N; j++)
89: {
90: c=gcoeff(tab,k,(i-1)*N+j);
91: if (signe(c))
92: {
93: p1 = gadd(gmul((GEN)x[i],(GEN)y[j]),
94: gmul((GEN)x[j],(GEN)y[i]));
95: if (!gcmp1(c)) p1 = gmul(p1,c);
96: s = gadd(s, p1);
97: }
98: }
99: }
100: v[k]=(long)gerepileupto(av,s); av=avma;
101: }
102: return v;
103: }
104:
105: /* inverse of x in nf */
106: GEN
107: element_inv(GEN nf, GEN x)
108: {
109: long av=avma,tetpil,flx,i,N,tx=typ(x);
110: GEN p1,p;
111:
112: nf=checknf(nf); N=lgef(nf[1])-3;
113: if (is_extscalar_t(tx))
114: {
115: if (tx==t_POLMOD) checknfelt_mod(nf,x);
116: else if (tx==t_POL) x=gmodulcp(x,(GEN)nf[1]);
117: p1=ginv(x); tetpil=avma;
118: return gerepile(av,tetpil,algtobasis(nf,p1));
119: }
120: if (isnfscalar(x))
121: {
122: p1=cgetg(N+1,t_COL); p1[1]=linv((GEN)x[1]);
123: for (i=2; i<=N; i++) p1[i]=lcopy((GEN)x[i]);
124: return p1;
125: }
126: flx=1;
127: for (i=1; i<=N; i++)
128: if (typ(x[i])==t_INTMOD)
129: {
130: p=gmael(x,i,1); x=lift(x);
131: flx=0; break;
132: }
133: p1 = ginvmod(gmul((GEN)nf[7],x), (GEN)nf[1]);
134: p1 = algtobasis_intern(nf,p1);
135:
136: if (flx) return gerepileupto(av,p1);
137: tetpil=avma; return gerepile(av,tetpil,Fp_vec(p1, p));
138: }
139:
140: /* quotient of x and y in nf */
141: GEN
142: element_div(GEN nf, GEN x, GEN y)
143: {
144: long av=avma,tetpil,flx,i,N,tx=typ(x),ty=typ(y);
145: GEN p1,p;
146:
147: nf=checknf(nf); N=lgef(nf[1])-3;
148: if (tx==t_POLMOD) checknfelt_mod(nf,x);
149: else if (tx==t_POL) x=gmodulcp(x,(GEN)nf[1]);
150:
151: if (ty==t_POLMOD) checknfelt_mod(nf,y);
152: else if (ty==t_POL) y=gmodulcp(y,(GEN)nf[1]);
153:
154: if (is_extscalar_t(tx))
155: {
156: if (is_extscalar_t(ty)) p1=gdiv(x,y);
157: else
158: {
159: if (ty!=t_COL) err(typeer,"nfdiv");
160: p1=gdiv(x,gmodulcp(gmul((GEN)nf[7],y),(GEN)nf[1]));
161: }
162: tetpil=avma; return gerepile(av,tetpil,algtobasis(nf,p1));
163: }
164: if (is_extscalar_t(ty))
165: {
166: if (tx!=t_COL) err(typeer,"nfdiv");
167: p1=gdiv(gmodulcp(gmul((GEN)nf[7],x),(GEN)nf[1]),y);
168: tetpil=avma; return gerepile(av,tetpil,algtobasis(nf,p1));
169: }
170:
171: if (isnfscalar(y)) return gdiv(x,(GEN)y[1]);
172: if (isnfscalar(x))
173: {
174: p1=element_inv(nf,y); tetpil=avma;
175: return gerepile(av,tetpil,gmul((GEN)x[1],p1));
176: }
177:
178: flx=1;
179: for (i=1; i<=N; i++)
180: if (typ(x[i])==t_INTMOD)
181: {
182: p=gmael(x,i,1); x=lift(x);
183: flx=0; break;
184: }
185: for (i=1; i<=N; i++)
186: if (typ(y[i])==t_INTMOD)
187: {
188: p=gmael(y,i,1); y=lift(y);
189: flx=0; break;
190: }
191:
192: p1 = gmul(gmul((GEN)nf[7],x), ginvmod(gmul((GEN)nf[7],y), (GEN)nf[1]));
193: p1 = algtobasis_intern(nf, gres(p1, (GEN)nf[1]));
194: if (flx) return gerepileupto(av,p1);
195: tetpil=avma; return gerepile(av,tetpil, Fp_vec(p1,p));
196: }
197:
198: /* product of INTEGERS (i.e vectors with integral coeffs) x and y in nf */
199: GEN
200: element_muli(GEN nf, GEN x, GEN y)
201: {
202: long av,i,j,k,N=lgef(nf[1])-3;
203: GEN p1,s,v,c,tab = (GEN)nf[9];
204:
205: v=cgetg(N+1,t_COL); av=avma;
206: for (k=1; k<=N; k++)
207: {
208: if (k == 1)
209: s = mulii((GEN)x[1],(GEN)y[1]);
210: else
211: s = addii(mulii((GEN)x[1],(GEN)y[k]),
212: mulii((GEN)x[k],(GEN)y[1]));
213: for (i=2; i<=N; i++)
214: {
215: c=gcoeff(tab,k,(i-1)*N+i);
216: if (signe(c))
217: {
218: p1 = mulii((GEN)x[i],(GEN)y[i]);
219: if (!gcmp1(c)) p1 = mulii(p1,c);
220: s = addii(s,p1);
221: }
222: for (j=i+1; j<=N; j++)
223: {
224: c=gcoeff(tab,k,(i-1)*N+j);
225: if (signe(c))
226: {
227: p1 = addii(mulii((GEN)x[i],(GEN)y[j]),
228: mulii((GEN)x[j],(GEN)y[i]));
229: if (!gcmp1(c)) p1 = mulii(p1,c);
230: s = addii(s,p1);
231: }
232: }
233: }
234: v[k]=(long) gerepileuptoint(av,s); av=avma;
235: }
236: return v;
237: }
238:
239: /* product of INTEGERS (i.e vectors with integral coeffs) x and y in nf */
240: GEN
241: element_sqri(GEN nf, GEN x)
242: {
243: long av,i,j,k,N=lgef(nf[1])-3;
244: GEN p1,s,v,c,tab = (GEN)nf[9];
245:
246: v=cgetg(N+1,t_COL); av=avma;
247: for (k=1; k<=N; k++)
248: {
249: if (k == 1)
250: s = sqri((GEN)x[1]);
251: else
252: s = shifti(mulii((GEN)x[1],(GEN)x[k]), 1);
253: for (i=2; i<=N; i++)
254: {
255: c=gcoeff(tab,k,(i-1)*N+i);
256: if (signe(c))
257: {
258: p1 = sqri((GEN)x[i]);
259: if (!gcmp1(c)) p1 = mulii(p1,c);
260: s = addii(s,p1);
261: }
262: for (j=i+1; j<=N; j++)
263: {
264: c=gcoeff(tab,k,(i-1)*N+j);
265: if (signe(c))
266: {
267: p1 = shifti(mulii((GEN)x[i],(GEN)x[j]),1);
268: if (!gcmp1(c)) p1 = mulii(p1,c);
269: s = addii(s,p1);
270: }
271: }
272: }
273: v[k]=(long) gerepileuptoint(av,s); av=avma;
274: }
275: return v;
276: }
277:
278: /* square of x in nf */
279: GEN
280: element_sqr(GEN nf, GEN x)
281: {
282: long av,i,j,k,N=lgef(nf[1])-3;
283: GEN p1,s,v,c, tab = (GEN)nf[9];
284:
285: if (isnfscalar(x))
286: {
287: s=cgetg(N+1,t_COL); s[1]=lsqr((GEN)x[1]);
288: for (i=2; i<=N; i++) s[i]=lcopy((GEN)x[i]);
289: return s;
290: }
291: v=cgetg(N+1,t_COL); av=avma;
292: for (k=1; k<=N; k++)
293: {
294: if (k == 1)
295: s = gsqr((GEN)x[1]);
296: else
297: s = gmul2n(gmul((GEN)x[1],(GEN)x[k]), 1);
298: for (i=2; i<=N; i++)
299: {
300: c=gcoeff(tab,k,(i-1)*N+i);
301: if (signe(c))
302: {
303: p1 = gsqr((GEN)x[i]);
304: if (!gcmp1(c)) p1 = gmul(p1,c);
305: s = gadd(s,p1);
306: }
307: for (j=i+1; j<=N; j++)
308: {
309: c=gcoeff(tab,k,(i-1)*N+j);
310: if (signe(c))
311: {
312: p1 = gmul((GEN)x[i],(GEN)x[j]);
313: p1 = gcmp1(c)? gmul2n(p1,1): gmul(p1,shifti(c,1));
314: s = gadd(s,p1);
315: }
316: }
317: }
318: v[k]=(long)gerepileupto(av,s); av=avma;
319: }
320: return v;
321: }
322:
323: /* Compute x^n in nf, left-shift binary powering */
324: GEN
325: element_pow(GEN nf, GEN x, GEN n)
326: {
327: long s,av=avma,N,i,j,m;
328: GEN y,p1;
329:
330: if (typ(n)!=t_INT) err(talker,"not an integer exponent in nfpow");
331: nf=checknf(nf); N=lgef(nf[1])-3;
332: s=signe(n); if (!s) return gscalcol_i(gun,N);
333: if (typ(x)!=t_COL) x=algtobasis(nf,x);
334:
335: if (isnfscalar(x))
336: {
337: y = gscalcol_i(gun,N);
338: y[1] = (long)powgi((GEN)x[1],n); return y;
339: }
340: p1 = n+2; m = *p1;
341: y=x; j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
342: for (i=lgefint(n)-2;;)
343: {
344: for (; j; m<<=1,j--)
345: {
346: y = element_sqr(nf, y);
347: if (m<0) y=element_mul(nf, y, x);
348: }
349: if (--i == 0) break;
350: m = *++p1, j = BITS_IN_LONG;
351: }
352: if (s<0) y = element_inv(nf, y);
353: return av==avma? gcopy(y): gerepileupto(av,y);
354: }
355:
356: /* x in Z^n, compute lift(x^n mod p) */
357: GEN
358: element_pow_mod_p(GEN nf, GEN x, GEN n, GEN p)
359: {
360: long s,av=avma,N,i,j,m;
361: GEN y,p1;
362:
363: if (typ(n)!=t_INT) err(talker,"not an integer exponent in nfpow");
364: nf=checknf(nf); N=lgef(nf[1])-3;
365: s=signe(n); if (!s) return gscalcol_i(gun,N);
366: if (typ(x)!=t_COL) x=algtobasis(nf,x);
367:
368: if (isnfscalar(x))
369: {
370: y = gscalcol_i(gun,N);
371: y[1] = (long)powmodulo((GEN)x[1],n,p); return y;
372: }
373: p1 = n+2; m = *p1;
374: y=x; j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
375: for (i=lgefint(n)-2;;)
376: {
377: for (; j; m<<=1,j--)
378: {
379: y = element_sqri(nf, y);
380: if (m<0) y=element_muli(nf, y, x);
381: y = Fp_vec_red(y, p);
382: }
383: if (--i == 0) break;
384: m = *++p1, j = BITS_IN_LONG;
385: }
386: if (s<0) y = Fp_vec_red(element_inv(nf,y), p);
387: return av==avma? gcopy(y): gerepileupto(av,y);
388: }
389:
390: /* x = I-th vector of the Z-basis of Z_K, in Z^n, compute lift(x^n mod p) */
391: GEN
392: element_powid_mod_p(GEN nf, long I, GEN n, GEN p)
393: {
394: long s,av=avma,N,i,j,m;
395: GEN y,p1;
396:
397: if (typ(n)!=t_INT) err(talker,"not an integer exponent in nfpow");
398: nf=checknf(nf); N=lgef(nf[1])-3;
399: s=signe(n);
400: if (!s || I == 1) return gscalcol_i(gun,N);
401: p1 = n+2; m = *p1;
402: y = zerocol(N); y[I] = un;
403: j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
404: for (i=lgefint(n)-2;;)
405: {
406: for (; j; m<<=1,j--)
407: {
408: y = element_sqri(nf, y);
409: if (m<0) y=element_mulid(nf, y, I);
410: y = Fp_vec_red(y, p);
411: }
412: if (--i == 0) break;
413: m = *++p1, j = BITS_IN_LONG;
414: }
415: if (s<0) y = Fp_vec_red(element_inv(nf,y), p);
416: return av==avma? gcopy(y): gerepileupto(av,y);
417: }
418:
419: /* Outputs x.w_i, where w_i is the i-th elt of the integral basis */
420: GEN
421: element_mulid(GEN nf, GEN x, long i)
422: {
423: long av,j,k,N;
424: GEN s,v,c,p1, tab;
425:
426: if (i==1) return gcopy(x);
427: N = lgef(nf[1])-3;
428: tab = (GEN)nf[9]; tab += (i-1)*N;
429: v=cgetg(N+1,t_COL); av=avma;
430: for (k=1; k<=N; k++)
431: {
432: s = gzero;
433: for (j=1; j<=N; j++)
434: if (signe(c = gcoeff(tab,k,j)) && !gcmp0(p1 = (GEN)x[j]))
435: {
436: if (!gcmp1(c)) p1 = gmul(p1,c);
437: s = gadd(s,p1);
438: }
439: v[k]=lpileupto(av,s); av=avma;
440: }
441: return v;
442: }
443:
444: /* valuation of integer x, with resp. to prime ideal P above p.
445: * p.P^(-1) = b Z_K, v <= val_p(norm(x)), and N = deg(nf)
446: */
447: long
448: int_elt_val(GEN nf, GEN x, GEN p, GEN b, long v)
449: {
450: long i,k,w, N = lgef(nf[1])-3;
451: GEN r,a,y, mul = cgetg(N+1,t_MAT);
452:
453: for (i=1; i<=N; i++) mul[i] = (long)element_mulid(nf,b,i);
454: y = cgetg(N+1, t_COL); /* will hold the new x */
455: x = dummycopy(x);
456: for(w=0; w<=v; w++)
457: {
458: for (i=1; i<=N; i++)
459: { /* compute (x.b)_i */
460: a = mulii((GEN)x[1], gcoeff(mul,i,1));
461: for (k=2; k<=N; k++) a = addii(a, mulii((GEN)x[k], gcoeff(mul,i,k)));
462: /* is it divisible by p ? */
463: y[i] = ldvmdii(a,p,&r);
464: if (signe(r)) return w;
465: }
466: r=x; x=y; y=r;
467: }
468: return w;
469: }
470:
471: long
472: element_val(GEN nf, GEN x, GEN vp)
473: {
474: long av = avma,N,w,vcx,e;
475: GEN cx,p;
476:
477: if (gcmp0(x)) return VERYBIGINT;
478: nf=checknf(nf); N=lgef(nf[1])-3;
479: checkprimeid(vp); p=(GEN)vp[1]; e=itos((GEN)vp[3]);
480: switch(typ(x))
481: {
482: case t_INT: case t_FRAC: case t_FRACN:
483: return ggval(x,p)*e;
484: case t_POLMOD: x = (GEN)x[2]; /* fall through */
485: case t_POL:
486: x = algtobasis_intern(nf,x); break;
487: case t_COL:
488: if (lg(x)==N+1) break;
489: default: err(typeer,"element_val");
490: }
491: if (isnfscalar(x)) return ggval((GEN)x[1],p)*e;
492:
493: cx = content(x);
494: if (gcmp1(cx)) vcx=0; else { x = gdiv(x,cx); vcx = ggval(cx,p); }
495: w = int_elt_val(nf,x,p,(GEN)vp[5],VERYBIGINT);
496: avma=av; return w + vcx*e;
497: }
498:
499: /* d = a multiple of norm(x), x integral */
500: long
501: element_val2(GEN nf, GEN x, GEN d, GEN vp)
502: {
503: GEN p = (GEN)vp[1];
504: long av, v = ggval(d,p);
505:
506: if (!v) return 0;
507: av=avma;
508: v = int_elt_val(nf,x,p,(GEN)vp[5],v);
509: avma=av; return v;
510: }
511:
512: /* polegal without comparing variables */
513: long
514: polegal_spec(GEN x, GEN y)
515: {
516: long i = lgef(x);
517:
518: if (i != lgef(y)) return 0;
519: for (i--; i > 1; i--)
520: if (!gegal((GEN)x[i],(GEN)y[i])) return 0;
521: return 1;
522: }
523:
524: GEN
525: basistoalg(GEN nf, GEN x)
526: {
527: long tx=typ(x),lx=lg(x),i;
528: GEN z;
529:
530: nf=checknf(nf);
531: switch(tx)
532: {
533: case t_COL:
534: for (i=1; i<lx; i++)
535: {
536: long t = typ(x[i]);
537: if (is_matvec_t(t)) break;
538: }
539: if (i==lx)
540: {
541: z = cgetg(3,t_POLMOD);
542: z[1] = lcopy((GEN)nf[1]);
543: z[2] = lmul((GEN)nf[7],x); return z;
544: }
545: /* fall through */
546:
547: case t_VEC: case t_MAT: z=cgetg(lx,tx);
548: for (i=1; i<lx; i++) z[i]=(long)basistoalg(nf,(GEN)x[i]);
549: return z;
550:
551: case t_POLMOD:
552: if (!polegal_spec((GEN)nf[1],(GEN)x[1]))
553: err(talker,"not the same number field in basistoalg");
554: return gcopy(x);
555: default: z=cgetg(3,t_POLMOD);
556: z[1]=lcopy((GEN)nf[1]);
557: z[2]=lmul(x,polun[varn(nf[1])]); return z;
558: }
559: }
560:
561: /* return the (N-dimensional) vector of coeffs of p */
562: GEN
563: pol_to_vec(GEN x, long N)
564: {
565: long i,l=lgef(x)-1;
566: GEN z = cgetg(N+1,t_COL);
567: x++;
568: for (i=1; i<l ; i++) z[i]=x[i];
569: for ( ; i<=N; i++) z[i]=zero;
570: return z;
571: }
572:
573: /* valid for scalars and polynomial, degree less than N.
574: * No garbage collecting. No check (SEGV for vectors).
575: */
576: GEN
577: algtobasis_intern(GEN nf,GEN x)
578: {
579: long i,tx=typ(x),N=lgef(nf[1])-3;
580: GEN z;
581:
582: if (tx==t_POLMOD) { x=(GEN)x[2]; tx=typ(x); }
583: if (tx==t_POL)
584: {
585: if (lgef(x)-3 >= N) x=gres(x,(GEN)nf[1]);
586: return gmul((GEN)nf[8], pol_to_vec(x, N));
587: }
588: z = cgetg(N+1,t_COL);
589: z[1]=lcopy(x); for (i=2; i<=N; i++) z[i]=zero;
590: return z;
591: }
592:
593: GEN
594: algtobasis(GEN nf, GEN x)
595: {
596: long tx=typ(x),lx=lg(x),av=avma,i,N;
597: GEN z;
598:
599: nf=checknf(nf);
600: switch(tx)
601: {
602: case t_VEC: case t_COL: case t_MAT:
603: z=cgetg(lx,tx);
604: for (i=1; i<lx; i++) z[i]=(long)algtobasis(nf,(GEN)x[i]);
605: return z;
606: case t_POLMOD:
607: if (!polegal_spec((GEN)nf[1],(GEN)x[1]))
608: err(talker,"not the same number field in algtobasis");
609: x = (GEN)x[2]; /* fall through */
610: case t_POL:
611: return gerepileupto(av,algtobasis_intern(nf,x));
612:
613: default: N=lgef(nf[1])-3; return gscalcol(x,N);
614: }
615: }
616:
617: /* Given a and b in nf, gives an algebraic integer y in nf such that a-b.y
618: * is "small"
619: */
620: GEN
621: nfdiveuc(GEN nf, GEN a, GEN b)
622: {
623: long av=avma, tetpil;
624: a = element_div(nf,a,b); tetpil=avma;
625: return gerepile(av,tetpil,ground(a));
626: }
627:
628: /* Given a and b in nf, gives a "small" algebraic integer r in nf
629: * of the form a-b.y
630: */
631: GEN
632: nfmod(GEN nf, GEN a, GEN b)
633: {
634: long av=avma,tetpil;
635: GEN p1=gneg_i(element_mul(nf,b,ground(element_div(nf,a,b))));
636: tetpil=avma; return gerepile(av,tetpil,gadd(a,p1));
637: }
638:
639: /* Given a and b in nf, gives a two-component vector [y,r] in nf such
640: * that r=a-b.y is "small".
641: */
642: GEN
643: nfdivres(GEN nf, GEN a, GEN b)
644: {
645: long av=avma,tetpil;
646: GEN p1,z, y = ground(element_div(nf,a,b));
647:
648: p1=gneg_i(element_mul(nf,b,y)); tetpil=avma;
649: z=cgetg(3,t_VEC); z[1]=lcopy(y); z[2]=ladd(a,p1);
650: return gerepile(av,tetpil,z);
651: }
652:
653: /*************************************************************************/
654: /** **/
655: /** (Z_K/I)^* **/
656: /** **/
657: /*************************************************************************/
658:
659: /* return (column) vector of R1 signatures of x (coeff modulo 2)
660: * if arch = NULL, assume arch = [0,..0]
661: */
662: GEN
663: zsigne(GEN nf,GEN x,GEN arch)
664: {
665: GEN _0,_1, vecsign, rac = (GEN)nf[6];
666: long av, i,j,l;
667:
668: if (!arch) return cgetg(1,t_COL);
669: switch(typ(x))
670: {
671: case t_COL: x = gmul((GEN)nf[7],x); break;
672: case t_POLMOD: x = (GEN)x[2];
673: }
674: if (gcmp0(x)) err(talker,"zero element in zsigne");
675:
676: l = lg(arch); vecsign = cgetg(l,t_COL);
677: _0 = gmodulss(0,2);
678: _1 = gmodulss(1,2); av = avma;
679: for (j=1,i=1; i<l; i++)
680: if (signe(arch[i]))
681: vecsign[j++] = (gsigne(poleval(x,(GEN)rac[i])) > 0)? (long)_0
682: : (long)_1;
683: avma = av; setlg(vecsign,j); return vecsign;
684: }
685:
686: /* For internal use. Reduce x modulo ideal. We want a non-zero result */
687: GEN
688: nfreducemodideal(GEN nf,GEN x,GEN ideal)
689: {
690: long N = lg(x)-1, do_copy = 1, i;
691: GEN p1,q;
692:
693: ideal=idealhermite(nf,ideal);
694: for (i=N; i>=1; i--)
695: {
696: p1=gcoeff(ideal,i,i); q=gdivround((GEN)x[i],p1);
697: if (signe(q)) { x=gsub(x,gmul(q,(GEN)ideal[i])); do_copy=0; }
698: }
699: if (gcmp0(x)) return (GEN) ideal[1];
700: return do_copy? gcopy(x) : x;
701: }
702:
703: /* a usage interne...Reduction de la colonne x modulo la matrice y inversible
704: utilisant LLL */
705: GEN
706: lllreducemodmatrix(GEN x,GEN y)
707: {
708: long av = avma,tetpil;
709: GEN z = gmul(y,lllint(y));
710:
711: z = gneg(gmul(z, ground(gauss(z,x))));
712: tetpil=avma; return gerepile(av,tetpil,gadd(x,z));
713: }
714:
715: /* Reduce column x modulo y in HNF */
716: static GEN
717: colreducemodmat(GEN x,GEN y)
718: {
719: long av = avma, i;
720: GEN q;
721:
722: for (i=lg(x)-1; i>=1; i--)
723: {
724: q = gdivround((GEN)x[i], gcoeff(y,i,i));
725: if (signe(q)) { q = gneg(q); x = gadd(x, gmul(q,(GEN)y[i])); }
726: }
727: return avma==av? gcopy(x): gerepileupto(av,x);
728: }
729:
730: /* for internal use...Reduce matrix x modulo matrix y */
731: GEN
732: reducemodmatrix(GEN x, GEN y)
733: {
734: long i, lx = lg(x);
735: GEN z = cgetg(lx,t_MAT);
736:
737: if (DEBUGLEVEL>=8)
738: {
739: fprintferr("entering reducemodmatrix; avma-bot = %ld\n",avma-bot);
740: flusherr();
741: }
742: y = hnfmod(y,detint(y));
743: for (i=1; i<lx; i++)
744: {
745: if(DEBUGLEVEL>=8) { fprintferr("%ld ",i); flusherr(); }
746: z[i]=(long)colreducemodmat((GEN)x[i],y);
747: }
748: if(DEBUGLEVEL>=8) { fprintferr("\n"); flusherr(); }
749: return z;
750: }
751:
752: /* compute elt = x mod idele, with sign(elt) = sign(x) at arch */
753: GEN
754: nfreducemodidele(GEN nf,GEN x,GEN idele,GEN sarch)
755: {
756: GEN p1,p2,arch,gen;
757: long nba,i;
758:
759: if (gcmp0(x)) return gcopy(x);
760: if (!sarch || typ(idele)!=t_VEC || lg(idele)!=3)
761: return nfreducemodideal(nf,x,idele);
762:
763: arch =(GEN)idele[2];
764: nba = lg(sarch[1]);
765: gen =(GEN)sarch[2];
766: p1 = nfreducemodideal(nf,x,(GEN)idele[1]);
767: p2 = gadd(zsigne(nf,p1,arch), zsigne(nf,x,arch));
768: p2 = lift_intern(gmul((GEN)sarch[3],p2));
769: for (i=1; i<nba; i++)
770: if (signe(p2[i])) p1 = element_mul(nf,p1,(GEN)gen[i]);
771: return (gcmp(gnorml2(p1),gnorml2(x)) > 0)? x: p1;
772: }
773:
774: GEN
775: element_powmodideal(GEN nf,GEN x,GEN k,GEN ideal)
776: {
777: GEN y = gscalcol_i(gun,lgef(nf[1])-3);
778: for(;;)
779: {
780: if (mpodd(k)) y=element_mulmodideal(nf,x,y,ideal);
781: k=shifti(k,-1); if (!signe(k)) return y;
782: x = element_sqrmodideal(nf,x,ideal);
783: }
784: }
785:
786: GEN
787: element_powmodidele(GEN nf,GEN x,GEN k,GEN idele,GEN sarch)
788: {
789: GEN y = gscalcol_i(gun,lgef(nf[1])-3);
790: for(;;)
791: {
792: if (mpodd(k)) y=element_mulmodidele(nf,x,y,idele,sarch);
793: k=shifti(k,-1); if (!signe(k)) return y;
794: x = element_sqrmodidele(nf,x,idele,sarch);
795: }
796: }
797:
798: /* given 2 integral ideals x, y in HNF s.t x|y|x^2, compute the quotient
799: (1+x)/(1+y) in the form [[cyc],[gen],ux^-1]. */
800: static GEN
801: zidealij(GEN x, GEN y)
802: {
803: GEN p1,p2,p3,p4,d,z,x1;
804: long j,N,c;
805:
806: if(DEBUGLEVEL>=6)
807: {fprintferr("entree dans zidealij; avma-bot = %ld\n",avma-bot); flusherr();}
808: x1 = ginv(x);
809: p1 = gmul(x1,y);
810: if(DEBUGLEVEL>=6)
811: {fprintferr("p1 = y/x; avma-bot = %ld\n",avma-bot); flusherr();}
812: p2 = smith2(p1);
813: if(DEBUGLEVEL>=6)
814: {fprintferr("p2 = smith2(p1); avma-bot = %ld\n",avma-bot); flusherr();}
815: p3 = ginv((GEN)p2[1]);
816: if(DEBUGLEVEL>=6)
817: {fprintferr("p3 = 1/p2[1]; avma-bot = %ld\n",avma-bot); flusherr();}
818: p3 = reducemodmatrix(p3,p1);
819: p3 = gmul(x,p3); N=lg(p3)-1;
820: if(DEBUGLEVEL>=6)
821: {fprintferr("p3 = x.p3 mod p1; avma-bot = %ld\n",avma-bot); flusherr();}
822: for (j=1; j<=N; j++) coeff(p3,1,j)=laddsi(1,gcoeff(p3,1,j));
823: p4 = smithclean(p2); d=(GEN)p4[3];
824:
825: z=cgetg(4,t_VEC); c=lg(d); p1=cgetg(c,t_VEC);
826: /* transform p3 in a vector (gen) */
827: p3[0] = evaltyp(t_VEC) | evallg(c);
828: for (j=1; j<c; j++) p1[j] = coeff(d,j,j);
829: z[1]=(long)p1;
830: z[2]=(long)p3;
831: z[3] = lmul((GEN)p4[1],x1); return z;
832: }
833:
834: #if 0
835: /* un element g generateur d'un p^k divisant x etant donne, calcule
836: un element congru a g modulo p^k et a 1 modulo x/p^k et de plus
837: positif aux places de arch */
838: GEN
839: zconvert(GEN nf,GEN uv,GEN x,GEN arch,GEN structarch,GEN g)
840: {
841: long i,nba;
842: GEN p1,p2,generator;
843:
844: p1=nfreducemodideal(nf,gadd((GEN)uv[1],element_mul(nf,g,(GEN)uv[2])),x);
845: nba=lg(structarch[1])-1; generator=(GEN)structarch[2];
846: p2 = zsigne(nf,p1,arch);
847: for (i=1; i<=nba; i++)
848: if (signe(p2[i])) p1 = element_mul(nf,p1,(GEN)generator[i]);
849: return p1;
850: }
851: #endif
852:
853: static GEN
854: get_multab(GEN nf, GEN x)
855: {
856: long lx = lg(x), i;
857: GEN multab = cgetg(lx, t_MAT);
858: for (i=1; i<lx; i++)
859: multab[i]=(long)element_mulid(nf,x,i);
860: return multab;
861: }
862:
863: /* return mat * y mod prh */
864: static GEN
865: mul_matvec_mod_pr(GEN mat, GEN y, GEN prh)
866: {
867: long av,i,j, lx = lg(mat);
868: GEN v, res = cgetg(lx,t_COL), p = gcoeff(prh,1,1);
869:
870: av = avma; (void)new_chunk((lx-1) * lgefint(p));
871: v = zerocol(lx-1);
872: for (i=lx-1; i; i--)
873: {
874: GEN p1 = (GEN)v[i], t = (GEN)prh[i];
875: for (j=1; j<lx; j++)
876: p1 = addii(p1, mulii(gcoeff(mat,i,j), (GEN)y[j]));
877: p1 = modii(p1, p);
878: if (p1 != gzero && is_pm1(t[i]))
879: {
880: for (j=1; j<i; j++)
881: v[j] = lsubii((GEN)v[j], mulii(p1, (GEN)t[j]));
882: p1 = gzero;
883: }
884: if (p1 == gzero) /* intended */
885: res[i] = zero;
886: else
887: av = res[i] = (long)icopy_av(p1, (GEN)av);
888: }
889: avma = av; return res;
890: }
891:
892: /* smallest integer n such that g0^n=x modulo p prime */
893: static GEN
894: Fp_shanks(GEN x,GEN g0,GEN p)
895: {
896: long av=avma,av1,lim,lbaby,i,k,c;
897: GEN p1,smalltable,giant,perm,v,g0inv;
898:
899: x = modii(x,p);
900: if (is_pm1(x) || egalii(p,gdeux)) { avma = av; return gzero; }
901: p1 = addsi(-1, p);
902: if (egalii(p1,x)) { avma = av; return shifti(p,-1); }
903: p1 = racine(p1);
904: if (cmpis(p1,LGBITS) >= 0) err(talker,"module too large in Fp_shanks");
905: lbaby=itos(p1)+1; smalltable=cgetg(lbaby+1,t_VEC);
906: g0inv = mpinvmod(g0, p); p1 = x;
907:
908: c = 3 * lgefint(p);
909: for (i=1;;i++)
910: {
911: av1 = avma;
912: if (is_pm1(p1)) { avma=av; return stoi(i-1); }
913: smalltable[i]=(long)p1; if (i==lbaby) break;
914: new_chunk(c); p1 = mulii(p1,g0inv);
915: avma = av1; p1 = resii(p1, p);
916: }
917: giant = resii(mulii(x, mpinvmod(p1,p)), p);
918: p1=cgetg(lbaby+1,t_VEC);
919: perm = gen_sort(smalltable, cmp_IND | cmp_C, cmpii);
920: for (i=1; i<=lbaby; i++) p1[i]=smalltable[perm[i]];
921: smalltable=p1; p1=giant;
922:
923: av1 = avma; lim=stack_lim(av1,2);
924: for (k=1;;k++)
925: {
926: i=tablesearch(smalltable,p1,cmpii);
927: if (i)
928: {
929: v=addis(mulss(lbaby-1,k),perm[i]);
930: return gerepileuptoint(av,addsi(-1,v));
931: }
932: p1 = resii(mulii(p1,giant), p);
933:
934: if (low_stack(lim, stack_lim(av1,2)))
935: {
936: if(DEBUGMEM>1) err(warnmem,"Fp_shanks, k = %ld", k);
937: p1 = gerepileuptoint(av1, p1);
938: }
939: }
940: }
941:
942: /* smallest integer n such that g0^n=x modulo pr, assume g0 reduced */
943: /* TODO: should be done in F_(p^f), not in Z_k mod pr (done for f=1) */
944: GEN
945: nfshanks(GEN nf,GEN x,GEN g0,GEN pr,GEN prhall)
946: {
947: long av=avma,av1,lim,lbaby,i,k, f = itos((GEN)pr[4]);
948: GEN p1,smalltable,giant,perm,v,g0inv,prh = (GEN)prhall[1];
949: GEN multab, p = (GEN)pr[1];
950:
951: x = lift_intern(nfreducemodpr(nf,x,prhall));
952: if (f == 1) return gerepileuptoint(av, Fp_shanks((GEN)x[1],(GEN)g0[1],p));
953: p1 = addsi(-1, gpowgs(p,f));
954: if (isnfscalar(x))
955: {
956: x = (GEN)x[1];
957: if (gcmp1(x) || egalii((GEN)pr[1], gdeux)) { avma = av; return gzero; }
958: if (egalii(x, p1)) return gerepileuptoint(av,shifti(p1,-1));
959: v = divii(p1, addsi(-1,p));
960: g0 = lift_intern((GEN)element_powmodpr(nf,g0,v,prhall)[1]);
961: return gerepileuptoint(av, mulii(v, Fp_shanks(x,g0,p)));
962: }
963: p1 = racine(p1);
964: if (cmpis(p1,LGBITS) >= 0) err(talker,"module too large in nfshanks");
965: lbaby=itos(p1)+1; smalltable=cgetg(lbaby+1,t_VEC);
966: g0inv = lift_intern(element_invmodpr(nf,g0,prhall));
967: p1 = x;
968:
969: multab = get_multab(nf, g0inv);
970: for (i=lg(multab)-1; i; i--)
971: multab[i]=(long)Fp_vec_red((GEN)multab[i], p);
972:
973: for (i=1;;i++)
974: {
975: if (isnfscalar(p1) && gcmp1((GEN)p1[1])) { avma=av; return stoi(i-1); }
976:
977: smalltable[i]=(long)p1; if (i==lbaby) break;
978: p1 = mul_matvec_mod_pr(multab,p1,prh);
979: }
980: giant=lift_intern(element_divmodpr(nf,x,p1,prhall));
981: p1=cgetg(lbaby+1,t_VEC);
982: perm = gen_sort(smalltable, cmp_IND | cmp_C, cmp_vecint);
983: for (i=1; i<=lbaby; i++) p1[i]=smalltable[perm[i]];
984: smalltable=p1; p1=giant;
985:
986: multab = get_multab(nf, giant);
987: for (i=lg(multab)-1; i; i--)
988: multab[i]=(long)Fp_vec_red((GEN)multab[i], p);
989:
990: av1 = avma; lim=stack_lim(av1,2);
991: for (k=1;;k++)
992: {
993: i=tablesearch(smalltable,p1,cmp_vecint);
994: if (i)
995: {
996: v=addis(mulss(lbaby-1,k),perm[i]);
997: return gerepileuptoint(av,addsi(-1,v));
998: }
999: p1 = mul_matvec_mod_pr(multab,p1,prh);
1000:
1001: if (low_stack(lim, stack_lim(av1,2)))
1002: {
1003: if(DEBUGMEM>1) err(warnmem,"nfshanks");
1004: p1 = gerepileupto(av1, p1);
1005: }
1006: }
1007: }
1008:
1009: GEN
1010: znlog(GEN x, GEN g0)
1011: {
1012: long av=avma;
1013: GEN p = (GEN)g0[1];
1014: if (typ(g0) != t_INTMOD)
1015: err(talker,"not an element of (Z/pZ)* in znlog");
1016: switch(typ(x))
1017: {
1018: case t_INT: break;
1019: default: x = gmul(x, gmodulsg(1,p));
1020: if (typ(x) != t_INTMOD)
1021: err(talker,"not an element of (Z/pZ)* in znlog");
1022: case t_INTMOD: x = (GEN)x[2]; break;
1023: }
1024: return gerepileuptoint(av, Fp_shanks(x,(GEN)g0[2],p));
1025: }
1026:
1027: GEN
1028: dethnf(GEN mat)
1029: {
1030: long av,i,l = lg(mat);
1031: GEN s;
1032:
1033: if (l<3) return l<2? gun: icopy(gcoeff(mat,1,1));
1034: av = avma; s = gcoeff(mat,1,1);
1035: for (i=2; i<l; i++) s = gmul(s,gcoeff(mat,i,i));
1036: return av==avma? gcopy(s): gerepileupto(av,s);
1037: }
1038:
1039: GEN
1040: dethnf_i(GEN mat)
1041: {
1042: long av,i,l = lg(mat);
1043: GEN s;
1044:
1045: if (l<3) return l<2? gun: icopy(gcoeff(mat,1,1));
1046: av = avma; s = gcoeff(mat,1,1);
1047: for (i=2; i<l; i++) s = mulii(s,gcoeff(mat,i,i));
1048: return gerepileuptoint(av,s);
1049: }
1050:
1051: static GEN
1052: makeprimetoideal(GEN nf,GEN id,GEN uv,GEN x)
1053: {
1054: GEN p1 = gadd((GEN)uv[1], element_mul(nf,x,(GEN)uv[2]));
1055: return nfreducemodideal(nf,p1,id);
1056: }
1057:
1058: static GEN
1059: makeprimetoidealvec(GEN nf,GEN ideal,GEN uv,GEN listgen)
1060: {
1061: long i, lx = lg(listgen);
1062: GEN y = cgetg(lx,t_VEC);
1063:
1064: for (i=1; i<lx; i++)
1065: y[i] = (long)makeprimetoideal(nf,ideal,uv,(GEN)listgen[i]);
1066: return y;
1067: }
1068:
1069: /* Given an ideal pr^ep, and an integral ideal x (in HNF form) compute a list
1070: * of vectors, each with 5 components as follows :
1071: * [[clh],[gen1],[gen2],[signat2],U.X^-1]. Each component corresponds to
1072: * d_i,g_i,g'_i,s_i. Generators g_i are not necessarily prime to x, the
1073: * generators g'_i are. signat2 is the (horizontal) vector made of the
1074: * signatures (column vectors) of the g'_i. If x = NULL, the original ideal
1075: * was a prime power
1076: */
1077: static GEN
1078: zprimestar(GEN nf,GEN pr,GEN ep,GEN x,GEN arch)
1079: {
1080: long av=avma,av1,N,f,nbp,j,n,m,tetpil,i,e,a,b;
1081: GEN prh,p,pf_1,list,v,p1,p2,p3,p4,prk,uv,g0,newgen,pra,prb;
1082: GEN *gptr[2];
1083:
1084: if(DEBUGLEVEL>=4)
1085: { fprintferr("treating pr = %Z ^ %Z\n",pr,ep); flusherr(); }
1086: prh=prime_to_ideal(nf,pr); N=lg(prh)-1;
1087: f=itos((GEN)pr[4]); p=(GEN)pr[1];
1088:
1089: pf_1 = addis(gpowgs(p,f), -1);
1090: v = zerocol(N);
1091: if (f==1) v[1]=gener(p)[2];
1092: else
1093: {
1094: GEN prhall = cgetg(3,t_VEC);
1095: long psim;
1096: if (is_bigint(p)) err(talker,"prime too big in zprimestar");
1097: psim = itos(p);
1098: list = (GEN)factor(pf_1)[1]; nbp=lg(list)-1;
1099: prhall[1]=(long)prh; prhall[2]=zero;
1100: for (n=psim; n>=0; n++)
1101: {
1102: m=n;
1103: for (i=1; i<=N; i++)
1104: if (!gcmp1(gcoeff(prh,i,i))) { v[i]=lstoi(m%psim); m/=psim; }
1105: for (j=1; j<=nbp; j++)
1106: {
1107: p1 = divii(pf_1,(GEN)list[j]);
1108: p1 = lift_intern(element_powmodpr(nf,v,p1,prhall));
1109: if (isnfscalar(p1) && gcmp1((GEN)p1[1])) break;
1110: }
1111: if (j>nbp) break;
1112: }
1113: if (n < 0) err(talker,"prime too big in zprimestar");
1114: }
1115: /* v generates (Z_K / pr)^* */
1116: if(DEBUGLEVEL>=4) {fprintferr("v calcule\n");flusherr();}
1117: e = itos(ep); prk=(e==1)? pr: idealpow(nf,pr,ep);
1118: if(DEBUGLEVEL>=4) {fprintferr("prk calcule\n");flusherr();}
1119: g0 = v;
1120: if (x)
1121: {
1122: uv = idealaddtoone(nf,prk,idealdivexact(nf,x,prk));
1123: g0 = makeprimetoideal(nf,x,uv,v);
1124: }
1125: if(DEBUGLEVEL>=4) {fprintferr("g0 calcule\n");flusherr();}
1126:
1127: list=cgetg(2,t_VEC);
1128: p1=cgetg(6,t_VEC); list[1]=(long)p1; p1[5]=un;
1129: p2=cgetg(2,t_VEC); p1[1]=(long)p2; p2[1]=(long)pf_1;
1130: p2=cgetg(2,t_VEC); p1[2]=(long)p2; p2[1]=(long)v;
1131: p2=cgetg(2,t_VEC); p1[3]=(long)p2; p2[1]=(long)g0;
1132: p2=cgetg(2,t_VEC); p1[4]=(long)p2; p2[1]=(long)zsigne(nf,g0,arch);
1133: if (e==1)
1134: {
1135: tetpil=avma; return gerepile(av,tetpil,gcopy(list));
1136: }
1137:
1138: a=1; b=2; av1=avma;
1139: pra = prh; prb = (e==2)? prk: idealpow(nf,pr,gdeux);
1140: for(;;)
1141: {
1142: if(DEBUGLEVEL>=4)
1143: {fprintferr("on traite a = %ld, b = %ld\n",a,b); flusherr();}
1144: p1 = zidealij(pra,prb);
1145: newgen = dummycopy((GEN)p1[2]);
1146: p3 = cgetg(lg(newgen),t_VEC);
1147: if(DEBUGLEVEL>=4) {fprintferr("zidealij fait\n"); flusherr();}
1148: for (i=1; i<lg(newgen); i++)
1149: {
1150: if (x) newgen[i]=(long)makeprimetoideal(nf,x,uv,(GEN)newgen[i]);
1151: p3[i]=(long)zsigne(nf,(GEN)newgen[i],arch);
1152: }
1153: p2=cgetg(2,t_VEC); p4=cgetg(6,t_VEC); p2[1]=(long)p4;
1154: p4[1] = p1[1];
1155: p4[2] = p1[2];
1156: p4[3] = (long)newgen;
1157: p4[4] = (long)p3;
1158: p4[5] = p1[3];
1159:
1160: a=b; b=min(e,b<<1); tetpil = avma;
1161: list = concat(list,p2);
1162: if (a==b) return gerepile(av,tetpil,list);
1163:
1164: pra = gcopy(prb);
1165: gptr[0]=&pra; gptr[1]=&list;
1166: gerepilemanysp(av1,tetpil,gptr,2);
1167: prb = (b==(a<<1))? idealpow(nf,pra,gdeux): prk;
1168: }
1169: }
1170:
1171: /* x ideal, compute elements in 1+x whose sign matrix is invertible */
1172: GEN
1173: zarchstar(GEN nf,GEN x,GEN arch,long nba)
1174: {
1175: long av,N,i,j,k,r,rr,limr,zk,lgmat;
1176: GEN p1,y,bas,genarch,mat,lambda;
1177:
1178: if (!nba)
1179: {
1180: y=cgetg(4,t_VEC);
1181: y[1]=lgetg(1,t_VEC);
1182: y[2]=lgetg(1,t_VEC);
1183: y[3]=lgetg(1,t_MAT); return y;
1184: }
1185: N=lgef(nf[1])-3; y=cgetg(4,t_VEC);
1186: p1=cgetg(nba+1,t_VEC); for (i=1; i<=nba; i++) p1[i]=deux;
1187: y[1]=(long)p1; av = avma;
1188: if (N==1)
1189: {
1190: p1 = scalarpol(subsi(1, shifti(gcoeff(x,1,1),1)), varn(nf[1]));
1191: p1 = gerepileupto(av, p1);
1192: y[2]=lgetg(2,t_VEC); mael(y,2,1) = (long)p1;
1193: y[3]=(long)gscalmat(gun,1); return y;
1194: }
1195: zk = ideal_is_zk(x,N);
1196: x = gmul(x,lllintpartial(x));
1197: genarch = cgetg(nba+1,t_VEC);
1198: mat = cgetg(nba+1,t_MAT); setlg(mat,2);
1199: for (i=1; i<=nba; i++) mat[i] = lgetg(nba+1, t_MAT);
1200: lambda = new_chunk(N+1);
1201:
1202: bas = dummycopy(gmael(nf,5,1)); r = lg(arch);
1203: for (k=1,i=1; i<r; i++)
1204: if (signe(arch[i]))
1205: {
1206: if (k < i)
1207: for (j=1; j<=N; j++) coeff(bas,k,j) = coeff(bas,i,j);
1208: k++;
1209: }
1210: r = nba+1;
1211: for (i=1; i<=N; i++) setlg(bas[i], r);
1212: bas = gmul(bas, x);
1213:
1214: for (lgmat=1,r=1, rr=3; ; r<<=1, rr=(r<<1)+1)
1215: {
1216: if (DEBUGLEVEL>5) fprintferr("zarchstar: r = %ld\n",r);
1217: p1 = gpuigs(stoi(rr),N);
1218: limr = (cmpis(p1,BIGINT) > 0)? BIGINT: p1[2];
1219: limr = (limr-1)>>1;
1220: k = zk? rr: 1; /* to avoid 0 */
1221: for ( ; k<=limr; k++)
1222: {
1223: long av1=avma, kk = k;
1224: GEN alpha = gzero;
1225: for (i=1; i<=N; i++)
1226: {
1227: lambda[i] = (kk+r)%rr - r; kk/=rr;
1228: if (lambda[i])
1229: alpha = gadd(alpha, gmulsg(lambda[i],(GEN)bas[i]));
1230: }
1231: for (i=1; i<=nba; i++)
1232: alpha[i] = ladd((GEN)alpha[i], gun);
1233: p1 = (GEN)mat[lgmat];
1234: for (i=1; i<=nba; i++)
1235: p1[i] = (gsigne((GEN)alpha[i]) > 0)? zero: un;
1236:
1237: if (ker_trivial_mod_p(mat, gdeux)) avma = av1;
1238: else
1239: { /* new vector indep. of previous ones */
1240: avma = av1; alpha = gzero;
1241: for (i=1; i<=N; i++)
1242: if (lambda[i])
1243: alpha = gadd(alpha, gmulsg(lambda[i],(GEN)x[i]));
1244: alpha[1] = ladd((GEN)alpha[1], gun);
1245: genarch[lgmat++] = (long)alpha;
1246: if (lgmat > nba)
1247: {
1248: GEN *gptr[2];
1249: mat = ginv(mat); gptr[0]=&genarch; gptr[1]=&mat;
1250: gerepilemany(av,gptr,2);
1251: y[2]=(long)genarch;
1252: y[3]=(long)mat; return y;
1253: }
1254: setlg(mat,lgmat+1);
1255: }
1256: }
1257: }
1258: }
1259:
1260: /* Retourne la decomposition de a sur les nbgen generateurs successifs
1261: * contenus dans list_set et si index !=0 on ne fait ce calcul que pour
1262: * l'ideal premier correspondant a cet index en mettant a 0 les autres
1263: * composantes
1264: */
1265: static GEN
1266: zinternallog(GEN nf,GEN list_set,long nbgen,GEN arch,GEN fa,GEN a,long index)
1267: {
1268: GEN prlist,ep,y,ainit,list,pr,prk,cyc,gen,psigne,p1,p2,p3;
1269: long av,nbp,cp,i,j,k;
1270:
1271: y = cgetg(nbgen+1,t_COL); cp=0; av=avma;
1272: prlist=(GEN)fa[1]; ep=(GEN)fa[2]; nbp=lg(ep)-1;
1273: i=typ(a); if (is_extscalar_t(i)) a = algtobasis(nf,a);
1274: if (DEBUGLEVEL>3)
1275: {
1276: fprintferr("entering zinternallog, "); flusherr();
1277: if (DEBUGLEVEL>5) fprintferr("with a = %Z\n",a);
1278: }
1279: ainit = a; psigne = zsigne(nf,ainit,arch);
1280: for (k=1; k<=nbp; k++)
1281: {
1282: list=(GEN)list_set[k];
1283: if (index && index!=k)
1284: {
1285: for (j=1; j<lg(list); j++)
1286: {
1287: cyc = gmael(list,j,1);
1288: for (i=1; i<lg(cyc); i++) y[++cp]=zero;
1289: }
1290: continue;
1291: }
1292: pr=(GEN)prlist[k]; prk=idealpow(nf,pr,(GEN)ep[k]);
1293: for (j=1; j<lg(list); j++)
1294: {
1295: p1 = (GEN)list[j]; cyc=(GEN)p1[1]; gen=(GEN)p1[2];
1296: if (j==1)
1297: {
1298: if (DEBUGLEVEL>5) { fprintferr("do nfshanks\n"); flusherr(); }
1299: a=ainit; p3=nfmodprinit(nf,pr);
1300: p3 = nfshanks(nf,a,(GEN)gen[1],pr,p3);
1301: }
1302: else
1303: {
1304: p3 = (GEN)a[1]; a[1] = laddsi(-1,(GEN)a[1]);
1305: p2 = gmul((GEN)p1[5],a); a[1] = (long)p3;
1306: if (lg(p2)!=lg(cyc)) err(bugparier,"zinternallog");
1307: p3 = (GEN)p2[1];
1308: }
1309: for(i=1;;)
1310: {
1311: p3 = modii(negi(p3), (GEN)cyc[i]);
1312: y[++cp] = lnegi(p3);
1313: if (signe(p3))
1314: {
1315: if (mpodd((GEN)y[cp])) psigne = gadd(psigne,gmael(p1,4,i));
1316: if (DEBUGLEVEL>5) fprintferr("do element_powmodideal\n");
1317: p3 = element_powmodideal(nf,(GEN)gen[i],p3,prk);
1318: a = element_mulmodideal(nf,a,p3,prk);
1319: }
1320: i++; if (i==lg(cyc)) break;
1321: p3 = (GEN)p2[i];
1322: }
1323: }
1324: }
1325: p1=lift_intern(gmul(gmael(list_set,nbp+1,3), psigne));
1326: avma=av; for (i=1; i<lg(p1); i++) y[++cp] = p1[i];
1327: if (DEBUGLEVEL>3) { fprintferr("leaving\n"); flusherr(); }
1328: for (i=1; i<=nbgen; i++) y[i] = licopy((GEN)y[i]);
1329: return y;
1330: }
1331:
1332: GEN
1333: compute_gen(GEN nf, GEN u1, GEN met, GEN gen, GEN module, long nbp, GEN sarch)
1334: {
1335: long i,j,s,nba, c = lg(met), lgen = lg(gen);
1336: GEN basecl=cgetg(c,t_VEC), unnf = gscalcol_i(gun, lgef(nf[1])-3);
1337: GEN arch,x,p1,p2,p3,p4,p5,generator;
1338:
1339: if (sarch)
1340: {
1341: arch = (GEN)module[2];
1342: x = (GEN)module[1];
1343: generator = (GEN)sarch[2];
1344: nba = lg(generator)-1;
1345: }
1346: else
1347: {
1348: x = (typ(module) == t_MAT)? module: (GEN)module[1];
1349: nba = 0;
1350: }
1351: for (j=1; j<c; j++)
1352: {
1353: GEN *op, plus = unnf, minus = unnf;
1354:
1355: for (i=1; i<lgen; i++)
1356: {
1357: p1 = gcoeff(u1,i,j);
1358: if (!(s = signe(p1))) continue;
1359:
1360: if (s > 0) op = + else { op = − p1 = negi(p1); }
1361: p1 = element_powmodidele(nf,(GEN)gen[i],p1,module,sarch);
1362: *op = *op==unnf? p1: element_mulmodidele(nf,*op,p1,module,sarch);
1363: }
1364: if (nbp)
1365: {
1366: p5=idealaddtoone_i(nf,minus,x);
1367: p4=element_div(nf,p5,minus); /* = minus^(-1) mod x */
1368: p3=element_mulmodideal(nf,plus,p4,x);
1369: }
1370: else p3=unnf;
1371:
1372: if (nba)
1373: { /* correct so that sign(p3) == sign(plus/minus) */
1374: p4=gadd(gadd(zsigne(nf,p3,arch),
1375: zsigne(nf,plus,arch)),
1376: zsigne(nf,minus,arch));
1377: p2=lift_intern(gmul((GEN)sarch[3],p4));
1378: for (i=1; i<=nba; i++)
1379: if (signe(p2[i])) p3=element_mul(nf,p3,(GEN)generator[i]);
1380: }
1381: basecl[j]=(long)p3;
1382: }
1383: return basecl;
1384: }
1385:
1386: /* Compute [[ideal,arch], [h,[cyc],[gen]], idealfact, [liste], U]
1387: gen not computed unless add_gen != 0 */
1388: GEN
1389: zidealstarinitall(GEN nf, GEN ideal,long add_gen)
1390: {
1391: long av=avma,tetpil,i,j,k,nba,nbp,R1,nbgen,cp;
1392: GEN p1,p2,p3,p4,y,h,met,u1,u1u2,u1u2clean,fa,fa2,ep,x,arch,gen,list,sarch;
1393:
1394: nf = checknf(nf); R1=itos(gmael(nf,2,1));
1395: if (typ(ideal)==t_VEC && lg(ideal)==3)
1396: {
1397: arch=(GEN)ideal[2]; ideal = (GEN)ideal[1];
1398: i = typ(arch);
1399: if (!is_vec_t(i) || lg(arch) != R1+1)
1400: err(talker,"incorrect archimedean component in zidealstarinit");
1401: for (nba=0,i=1; i<=R1; i++)
1402: if (signe(arch[i])) nba++;
1403: }
1404: else
1405: {
1406: arch=cgetg(R1+1,t_VEC);
1407: for (nba=0,i=1; i<=R1; i++) arch[i]=zero;
1408: }
1409: x = idealhermite(nf,ideal);
1410: if (!gcmp1(denom(x)))
1411: err(talker,"zidealstarinit needs an integral ideal. x =\n%Z\n",x);
1412: p1=cgetg(3,t_VEC); ideal=p1;
1413: p1[1]=(long)x;
1414: p1[2]=(long)arch;
1415:
1416: fa=idealfactor(nf,x); list=(GEN)fa[1]; ep=(GEN)fa[2];
1417: nbp=lg(list)-1; fa2=cgetg(nbp+2,t_VEC);
1418:
1419: /* compute them */
1420: gen = cgetg(1,t_VEC);
1421: p2 = (nbp==1)? (GEN)NULL: x;
1422: for (i=1; i<=nbp; i++)
1423: {
1424: p1 = zprimestar(nf,(GEN)list[i],(GEN)ep[i],p2,arch);
1425: fa2[i]=(long)p1;
1426: for (j=1; j<lg(p1); j++)
1427: gen = concatsp(gen,gmael(p1,j,3));
1428: }
1429: sarch = zarchstar(nf,x,arch,nba);
1430: fa2[i]=(long)sarch;
1431: gen = concatsp(gen,(GEN)sarch[2]);
1432: nbgen = lg(gen)-1;
1433: h=cgetg(nbgen+1,t_MAT); cp=0;
1434: for (i=1; i<=nbp; i++)
1435: {
1436: list=(GEN)fa2[i];
1437: for (j=1; j<lg(list); j++)
1438: {
1439: p1=(GEN)list[j]; p2=(GEN)p1[1]; p3=(GEN)p1[3];
1440: for (k=1; k<lg(p3); k++)
1441: {
1442: if (DEBUGLEVEL>=6)
1443: { fprintferr("entering element_powmodidele\n"); flusherr(); }
1444: p4 = element_powmodidele(nf,(GEN)p3[k],(GEN)p2[k],ideal,sarch);
1445: h[++cp] = lneg(zinternallog(nf,fa2,nbgen,arch,fa,p4,i));
1446: coeff(h,cp,cp) = p2[k];
1447: }
1448: }
1449: }
1450: for (j=1; j<=nba; j++)
1451: { h[++cp]=(long)zerocol(nbgen); coeff(h,cp,cp)=deux; }
1452: if (cp!=nbgen) err(talker,"bug in zidealstarinit");
1453: u1u2 = smith2(h); u1u2clean = smithclean(u1u2);
1454: met=(GEN)u1u2clean[3];
1455: if (add_gen)
1456: {
1457: u1 = reducemodmatrix(ginv((GEN)u1u2[1]),h);
1458: p1 = cgetg(4,t_VEC);
1459: p1[3] = (long)compute_gen(nf,u1,met,gen,ideal, nbp,sarch);
1460: }
1461: else p1 = cgetg(3,t_VEC);
1462: y = cgetg(6,t_VEC);
1463: y[1]=(long)ideal;
1464: y[2]=(long)p1;
1465: p1[1]=(long)dethnf(met);
1466: p1[2]=(long)mattodiagonal(met);
1467: y[3]=(long)fa;
1468: y[4]=(long)fa2;
1469: y[5] = u1u2clean[1];
1470: tetpil=avma; return gerepile(av,tetpil,gcopy(y));
1471: }
1472:
1473: GEN
1474: zidealstarinitgen(GEN nf, GEN ideal)
1475: {
1476: return zidealstarinitall(nf,ideal,1);
1477: }
1478:
1479: GEN
1480: zidealstarinit(GEN nf, GEN ideal)
1481: {
1482: return zidealstarinitall(nf,ideal,0);
1483: }
1484:
1485: GEN
1486: zidealstar(GEN nf, GEN ideal)
1487: {
1488: long av = avma,tetpil;
1489: GEN y = zidealstarinitall(nf,ideal,1);
1490: tetpil=avma; return gerepile(av,tetpil,gcopy((GEN)y[2]));
1491: }
1492:
1493: GEN
1494: idealstar0(GEN nf, GEN ideal,long flag)
1495: {
1496: switch(flag)
1497: {
1498: case 0: return zidealstar(nf,ideal);
1499: case 1: return zidealstarinit(nf,ideal);
1500: case 2: return zidealstarinitgen(nf,ideal);
1501: default: err(flagerr,"idealstar");
1502: }
1503: return NULL; /* not reached */
1504: }
1505:
1506: /* x is not integral, but check that v_p(x)=0 for all prime divisors of the
1507: * ideal. We need x = a/b with integral a and b, prime to ideal
1508: * denmat = den * id.
1509: */
1510: static GEN
1511: rat_zinternallog(GEN nf, GEN x, GEN bigideal, GEN denmat)
1512: {
1513: long nbp,i,v,k, N = lgef(nf[1])-3;
1514: GEN den,fa,list,ep,pr,p1,p2,p3,x1,dinv,ideal;
1515:
1516: ideal = (GEN)bigideal[1];
1517: if (lg(ideal) == 3) ideal = (GEN)ideal[1];
1518: fa=(GEN)bigideal[3]; list=(GEN)fa[1]; ep=(GEN)fa[2];
1519: den=gmael(denmat,1,1); k=1; nbp=lg(list)-1;
1520: for (i=1; i<=nbp; i++)
1521: {
1522: pr=(GEN)list[i];
1523: v = (ggval(den,(GEN)pr[1])*itos((GEN)pr[3])) / itos((GEN)ep[i]) + 1;
1524: if (v>k) k=v;
1525: }
1526: p3=idealpow(nf,ideal,stoi(k));
1527: p1=idealadd(nf,denmat,p3); dinv=idealinv(nf,p1);
1528: p2=idealmullll(nf,denmat,dinv);
1529: p3=idealmullll(nf,p3,dinv);
1530: x1=idealaddtoone_i(nf,p2,p3);
1531: if (gcmp0(x1)) x1 = (GEN)denmat[1];
1532: /* x = a/b is suitable, with a=x1*x, b=x1 */
1533: p1=element_mul(nf,x1,x);
1534: /* x1 is prime to ideal. Check that x1 * x also */
1535: p2=idealadd(nf,p1,ideal);
1536: if (! ideal_is_zk(p2,N))
1537: err(talker,"element is not coprime to ideal in zideallog");
1538: p1=zideallog(nf,p1,bigideal);
1539: p2=zideallog(nf,x1,bigideal);
1540: return gsub(p1,p2);
1541: }
1542:
1543: /* Given x (not necessarily integral), and bid as output by zidealstarinit,
1544: * compute the vector of components on the generators bid[2]
1545: */
1546: GEN
1547: zideallog(GEN nf, GEN x, GEN bid)
1548: {
1549: long av,l,i,N,c;
1550: GEN fa,fa2,ideal,arch,den,p1,cyc,y;
1551:
1552: nf=checknf(nf); checkbid(bid);
1553: cyc=gmael(bid,2,2); c=lg(cyc);
1554: y=cgetg(c,t_COL); av=avma;
1555: N = lgef(nf[1])-3; ideal = (GEN) bid[1];
1556: if (typ(ideal)==t_VEC && lg(ideal)==3)
1557: arch = (GEN)ideal[2];
1558: else
1559: arch = NULL;
1560: switch(typ(x))
1561: {
1562: case t_INT: case t_FRAC: case t_FRACN:
1563: x = gscalcol_i(x,N); break;
1564: case t_POLMOD: case t_POL:
1565: x = algtobasis(nf,x); break;
1566: case t_COL: break;
1567: default: err(talker,"not an element in zideallog");
1568: }
1569: if (lg(x) != N+1) err(talker,"not an element in zideallog");
1570:
1571: den=denom(x);
1572: if (!gcmp1(den))
1573: p1 = rat_zinternallog(nf,x,bid, gscalmat(den,N));
1574: else
1575: {
1576: l=lg(bid[5])-1; fa=(GEN)bid[3]; fa2=(GEN)bid[4];
1577: p1 = zinternallog(nf,fa2,l,arch,fa,x,0);
1578: p1 = gmul((GEN)bid[5],p1); /* apply smith */
1579: }
1580: if (lg(p1)!=c) err(bugparier,"zideallog");
1581: for (i=1; i<c; i++)
1582: y[i] = lmodii((GEN)p1[i],(GEN)cyc[i]);
1583: avma=av; /* following line does a gerepile ! */
1584: for (i=1; i<c; i++)
1585: y[i] = (long)icopy((GEN)y[i]);
1586: return y;
1587: }
1588:
1589: /* Etant donnes bid1, bid2 resultats de zidealstarinit pour deux modules m1
1590: * et m2 premiers entre eux sans partie archimedienne, calcule le
1591: * zidealstarinit [[ideal,arch],[h,[cyc],[gen]],idealfact,[liste],U] du
1592: * produit
1593: */
1594: static GEN
1595: zidealstarinitjoinall(GEN nf, GEN bid1, GEN bid2, long add_gen)
1596: {
1597: long av=avma,i,j,nbp,nbgen,lx1,lx2,llx1,llx2,lx,llx;
1598: GEN module1,module2,struct1,struct2,fact1,fact2,liste1,liste2,U1,U2;
1599: GEN module,liste,fact,U,cyc,ex1,ex2,uv;
1600: GEN p1,p2,y,met,u1;
1601: GEN fa1,fa2,gen,u1u2,u1u2clean;
1602:
1603: nf=checknf(nf); checkbid(bid1); checkbid(bid2);
1604: module1=(GEN)bid1[1]; struct1=(GEN)bid1[2]; fact1=(GEN)bid1[3];
1605: module2=(GEN)bid2[1]; struct2=(GEN)bid2[2]; fact2=(GEN)bid2[3];
1606: module=cgetg(3,t_VEC);
1607: module[1]=(long)idealmul(nf,(GEN)module1[1],(GEN)module2[1]);
1608: module[2]=ladd((GEN)module1[2],(GEN)module2[2]);
1609: if (gcmpgs(vecmax((GEN)module[2]),1)>=0)
1610: err(talker,"nontrivial Archimedian components in zidealstarinitjoin");
1611:
1612: fa1=(GEN)fact1[1]; ex1=(GEN)fact1[2];
1613: fa2=(GEN)fact2[1]; ex2=(GEN)fact2[2];
1614: fact=cgetg(3,t_MAT);
1615: fact[1]=lconcat(fa1,fa2); lx1=lg(fa1);
1616: fact[2]=lconcat(ex1,ex2); lx2=lg(fa2);
1617: nbp=lx1+lx2-2;
1618: for (i=1; i<lx1; i++)
1619: if (isinvector(fa2,(GEN)fa1[i],lx2-1))
1620: err(talker,"noncoprime ideals in zidealstarinitjoin");
1621:
1622: liste1=(GEN)bid1[4]; lx1=lg(liste1);
1623: liste2=(GEN)bid2[4]; lx2=lg(liste2);
1624: lx=lx1+lx2-2; liste = cgetg(lx,t_VEC);
1625: for (i=1; i<lx1-1; i++) liste[i]=liste1[i];
1626: for ( ; i<lx; i++) liste[i]=liste2[i-lx1+2];
1627: U1=(GEN)bid1[5]; lx1=lg(U1);
1628: U2=(GEN)bid2[5]; lx2=lg(U2);
1629: lx=lx1+lx2-1; U = cgetg(lx,t_MAT);
1630:
1631: llx1=lg(struct1[2]);
1632: llx2=lg(struct2[2]);
1633: llx=llx1+llx2-1; nbgen=llx-1;
1634: cyc = diagonal(concatsp((GEN)struct1[2],(GEN)struct2[2]));
1635: u1u2=smith2(cyc); u1u2clean=smithclean(u1u2);
1636: met=(GEN)u1u2clean[3];
1637: if (nbgen)
1638: {
1639: for (j=1; j<lx1; j++)
1640: {
1641: p1=cgetg(llx,t_COL); U[j]=(long)p1;
1642: p2=(GEN)U1[j];
1643: for (i=1; i<llx1; i++) p1[i]=p2[i];
1644: for ( ; i<llx; i++) p1[i]=zero;
1645: }
1646: for ( ; j<lx; j++)
1647: {
1648: p1=cgetg(llx,t_COL); U[j]=(long)p1;
1649: p2=((GEN)U2[j-(lx1-1)]) + 1-llx1;
1650: for (i=1; i<llx1; i++) p1[i]=zero;
1651: for ( ; i<llx; i++) p1[i]=p2[i];
1652: }
1653: U = gmul((GEN)u1u2clean[1],U);
1654: }
1655: else
1656: for (j=1; j<lx; j++) U[j]=lgetg(1,t_COL);
1657:
1658: if (add_gen)
1659: {
1660: if (lg(struct1)<=3 || lg(struct2)<=3)
1661: err(talker,"please apply idealstar(,,2) and not idealstar(,,1)");
1662:
1663: uv = idealaddtoone(nf,(GEN)module1[1],(GEN)module2[1]);
1664: p1 = makeprimetoidealvec(nf,(GEN)module[1],uv,(GEN)struct1[3]);
1665: p2=(GEN)uv[1]; uv[1]=uv[2]; uv[2]=(long)p2;
1666: p2 = makeprimetoidealvec(nf,(GEN)module[1],uv,(GEN)struct2[3]);
1667: gen=concatsp(p1,p2);
1668:
1669: u1 = reducemodmatrix(ginv((GEN)u1u2[1]),cyc);
1670: p1 = cgetg(4,t_VEC);
1671: p1[3] = (long)compute_gen(nf,u1,met,gen,module, nbp,NULL);
1672: }
1673: else p1 = cgetg(3,t_VEC);
1674:
1675: y=cgetg(6,t_VEC);
1676: y[1]=(long)module;
1677: y[2]=(long)p1;
1678: p1[1]=(long)dethnf(met);
1679: p1[2]=(long)mattodiagonal(met);
1680: y[3]=(long)fact;
1681: y[4]=(long)liste;
1682: y[5]=(long)U;
1683: return gerepileupto(av,gcopy(y));
1684: }
1685:
1686: GEN
1687: zidealstarinitjoin(GEN nf, GEN bid1, GEN bid2)
1688: {
1689: return zidealstarinitjoinall(nf,bid1,bid2,0);
1690: }
1691:
1692: GEN
1693: zidealstarinitjoingen(GEN nf, GEN bid1, GEN bid2)
1694: {
1695: return zidealstarinitjoinall(nf,bid1,bid2,1);
1696: }
1697:
1698: /* Etant donnes bid1 resultat de zidealstarinit pour un module m1 sans partie
1699: * archimedienne et une partie archimedienne arch, calcule le zidealstarinit
1700: * [[ideal,arch],[h,[cyc],[gen]],idealfact,[liste],U] du produit
1701: */
1702: static GEN
1703: zidealstarinitjoinarchall(GEN nf, GEN bid1, GEN arch, long nba, long add_gen)
1704: {
1705: long av=avma,i,nbp,lx1,lx;
1706: GEN module1,struct1,fact1,liste1,U1;
1707: GEN module,liste,h,p1,y,met,u1,x,gen,sarch,u1u2,u1u2clean;
1708:
1709: nf=checknf(nf); checkbid(bid1);
1710: module1=(GEN)bid1[1]; struct1=(GEN)bid1[2]; fact1=(GEN)bid1[3];
1711: nbp=lg((GEN)fact1[1])-1;
1712: x=(GEN)module1[1];
1713: sarch = zarchstar(nf,x,arch,nba);
1714: module=cgetg(3,t_VEC);
1715: module[1]=(long)x;
1716: module[2]=(long)arch;
1717: if (gcmpgs(vecmax((GEN)module1[2]),1)>=0)
1718: err(talker,"nontrivial Archimedian components in zidealstarinitjoinarchall");
1719: liste1=(GEN)bid1[4]; lx=lg(liste1);
1720: liste=cgetg(lx,t_VEC); for (i=1; i<lx-1; i++) liste[i]=liste1[i];
1721: liste[lx-1]=(long)sarch;
1722: U1=(GEN)bid1[5]; lx1=lg(U1);
1723: lx=lx1+nba;
1724: h = diagonal(concatsp((GEN)struct1[2],(GEN)sarch[1]));
1725: u1u2 = smith2(h); u1u2clean = smithclean(u1u2);
1726: met=(GEN)u1u2clean[3];
1727:
1728: if (add_gen)
1729: {
1730: if (lg(struct1)<=3)
1731: err(talker,"please apply idealstar(,,2) and not idealstar(,,1)");
1732: gen=concatsp((GEN)struct1[3],(GEN)sarch[2]);
1733:
1734: u1 = reducemodmatrix(ginv((GEN)u1u2[1]),h);
1735: p1 = cgetg(4,t_VEC);
1736: p1[3] = (long)compute_gen(nf,u1,met,gen,module, nbp,sarch);
1737: }
1738: else p1 = cgetg(3,t_VEC);
1739:
1740: y=cgetg(6,t_VEC);
1741: y[1]=(long)module;
1742: y[2]=(long)p1;
1743: p1[1]=(long)dethnf(met);
1744: p1[2]=(long)mattodiagonal(met);
1745: y[3]=(long)fact1;
1746: y[4]=(long)liste;
1747: y[5] = u1u2clean[1];
1748: return gerepileupto(av,gcopy(y));
1749: }
1750:
1751: GEN
1752: zidealstarinitjoinarch(GEN nf, GEN bid1, GEN arch, long nba)
1753: {
1754: return zidealstarinitjoinarchall(nf,bid1,arch,nba,0);
1755: }
1756:
1757: GEN
1758: zidealstarinitjoinarchgen(GEN nf, GEN bid1, GEN arch, long nba)
1759: {
1760: return zidealstarinitjoinarchall(nf,bid1,arch,nba,1);
1761: }
1762:
1763: /* calcule la matrice des zinternallog des unites */
1764: GEN
1765: logunitmatrix(GEN nf,GEN funits,GEN racunit,GEN bid)
1766: {
1767: long R,j,sizeh;
1768: GEN m,fa2,fa,arch;
1769:
1770: R=lg(funits)-1; m=cgetg(R+2,t_MAT);
1771: fa2=(GEN)bid[4]; sizeh=lg(bid[5])-1; arch=gmael(bid,1,2);
1772: fa=(GEN)bid[3];
1773: m[1]=(long)zinternallog(nf,fa2,sizeh,arch,fa,racunit,0);
1774: for (j=2; j<=R+1; j++)
1775: m[j]=(long)zinternallog(nf,fa2,sizeh,arch,fa,(GEN)funits[j-1],0);
1776: return m;
1777: }
1778:
1779: /* calcule la matrice des zinternallog des unites */
1780: static GEN
1781: logunitmatrixarch(GEN nf,GEN funits,GEN racunit,GEN bid)
1782: {
1783: long R,j;
1784: GEN m,liste,structarch,arch;
1785:
1786: R=lg(funits)-1; m=cgetg(R+2,t_MAT); arch=gmael(bid,1,2);
1787: liste=(GEN)bid[4]; structarch=(GEN)liste[lg(liste)-1];
1788: m[1]=(long)zsigne(nf,racunit,arch);
1789: for (j=2; j<=R+1; j++)
1790: m[j]=(long)zsigne(nf,(GEN)funits[j-1],arch);
1791: return lift_intern(gmul((GEN)structarch[3],m));
1792: }
1793:
1794: /* concatenation verticale de Q1 et Q2. Ne cree pas le resultat. */
1795: GEN
1796: vconcat(GEN Q1, GEN Q2)
1797: {
1798: long lc,lr,lx1,lx2,i,j;
1799: GEN m,p1,p2,p3;
1800:
1801: lc=lg(Q1); if (lc==1) return Q1;
1802: lx1=lg(Q1[1]); lx2=lg(Q2[1]); lr=lx1+lx2-1;
1803: m=cgetg(lc,t_MAT);
1804: for (j=1; j<lc; j++)
1805: {
1806: p1=cgetg(lr,t_COL); m[j]=(long)p1; p2=(GEN)Q1[j]; p3=(GEN)Q2[j];
1807: for (i=1; i<lx1; i++) p1[i]=p2[i];
1808: for ( ; i<lr; i++) p1[i]=p3[i-lx1+1];
1809: }
1810: return m;
1811: }
1812:
1813: static void
1814: init_units(GEN bnf, GEN *funits, GEN *racunit)
1815: {
1816: GEN p1;
1817: bnf = checkbnf(bnf); p1=(GEN)bnf[8];
1818: if (lg(p1)==5) *funits=(GEN)buchfu(bnf)[1];
1819: else
1820: {
1821: if (lg(p1)!=7) err(talker,"incorrect big number field");
1822: *funits=(GEN)p1[5];
1823: }
1824: *racunit=gmael(p1,4,2);
1825: }
1826:
1827: /* flag &1 : generateurs, sinon non
1828: * flag &2 : unites, sinon pas.
1829: * flag &4 : ideaux, sinon zidealstar.
1830: */
1831: static GEN
1832: ideallistzstarall(GEN bnf,long bound,long flag)
1833: {
1834: byteptr ptdif=diffptr;
1835: long lim,av0=avma,av,i,j,k,l,q2,lp1,q;
1836: long do_gen = flag & 1, do_units = flag & 2, big_id = !(flag & 4);
1837: GEN y,nf,p,z,z2,p1,p2,p3,fa,pr,ideal,lu,lu2,funits,racunit,embunit;
1838:
1839: nf=checknf(bnf);
1840: z = cgetg(bound+1,t_VEC);
1841: for (i=2; i<=bound; i++) z[i] = lgetg(1,t_VEC);
1842: if (do_units)
1843: {
1844: init_units(bnf,&funits,&racunit);
1845: lu = cgetg(bound+1,t_VEC);
1846: for (i=2; i<=bound; i++) lu[i]=lgetg(1,t_VEC);
1847: }
1848:
1849: ideal = idmat(lgef(nf[1])-3);
1850: if (big_id) ideal = zidealstarinitall(nf,ideal,do_gen);
1851: p2 = cgetg(2,t_VEC); z[1] = (long)p2;
1852: p2[1] = (long)ideal;
1853: if (do_units)
1854: {
1855: p2 = cgetg(2,t_VEC); lu[1] = (long)p2;
1856: p2[1] = (long)logunitmatrix(nf,funits,racunit,ideal);
1857: }
1858:
1859: p=cgeti(3); p[1]=evalsigne(1) | evallgefint(3);
1860: av=avma; lim=stack_lim(av,1);
1861: for (p[2]=0; p[2]<=bound; )
1862: {
1863: if (!*ptdif) err(primer1);
1864: p[2] += *ptdif++;
1865: if (DEBUGLEVEL>1) { fprintferr("%ld ",p[2]); flusherr(); }
1866: fa = primedec(nf,p);
1867: for (j=1; j<lg(fa); j++)
1868: {
1869: pr = (GEN)fa[j]; p1 = powgi(p,(GEN)pr[4]);
1870: if (is_bigint(p1) || (q = itos(p1)) > bound) continue;
1871:
1872: q2=q; ideal=pr; z2=dummycopy(z);
1873: if (do_units) lu2=dummycopy(lu);
1874: for (l=2; ;l++)
1875: {
1876: if (big_id) ideal = zidealstarinitall(nf,ideal,do_gen);
1877: if (do_units) embunit = logunitmatrix(nf,funits,racunit,ideal);
1878: for (i=q; i<=bound; i+=q)
1879: {
1880: p1 = (GEN)z[i/q];
1881: if ((lp1 = lg(p1)) == 1) continue;
1882:
1883: p2 = cgetg(lp1,t_VEC);
1884: for (k=1; k<lp1; k++)
1885: if (big_id)
1886: p2[k] = (long)zidealstarinitjoinall(nf,(GEN)p1[k],ideal,do_gen);
1887: else
1888: p2[k] = (long)idealmul(nf,(GEN)p1[k],ideal);
1889: z2[i] = (long)concatsp((GEN)z2[i],p2);
1890: if (do_units)
1891: {
1892: p1 = (GEN)lu[i/q];
1893: p2 = cgetg(lp1,t_VEC);
1894: for (k=1; k<lp1; k++)
1895: p2[k] = (long)vconcat((GEN)p1[k],embunit);
1896: lu2[i] = (long)concatsp((GEN)lu2[i],p2);
1897: }
1898: }
1899: q *= q2; if ((ulong)q > (ulong)bound) break;
1900: ideal = idealpows(nf,pr,l);
1901: }
1902: z = z2; if (do_units) lu = lu2;
1903: }
1904: if (low_stack(lim, stack_lim(av,1)))
1905: {
1906: GEN *gptr[2]; gptr[0]=&z; gptr[1]=&lu;
1907: if(DEBUGMEM>1) err(warnmem,"ideallistzstarall");
1908: gerepilemany(av,gptr,do_units?2:1);
1909: }
1910: }
1911: if (!do_units) return gerepileupto(av0, gcopy(z));
1912: y = cgetg(3,t_VEC);
1913: y[1] = lcopy(z);
1914: lu2 = cgetg(lg(z),t_VEC);
1915: for (i=1; i<lg(z); i++)
1916: {
1917: p1=(GEN)z[i]; p2=(GEN)lu[i]; lp1=lg(p1);
1918: p3=cgetg(lp1,t_VEC); lu2[i]=(long)p3;
1919: for (j=1; j<lp1; j++) p3[j] = lmul(gmael(p1,j,5),(GEN)p2[j]);
1920: }
1921: y[2]=(long)lu2; return gerepileupto(av0, y);
1922: }
1923:
1924: GEN
1925: ideallist0(GEN bnf,long bound, long flag)
1926: {
1927: if (flag<0 || flag>4) err(flagerr,"ideallist");
1928: return ideallistzstarall(bnf,bound,flag);
1929: }
1930:
1931: GEN
1932: ideallistzstar(GEN nf,long bound)
1933: {
1934: return ideallistzstarall(nf,bound,0);
1935: }
1936:
1937: GEN
1938: ideallistzstargen(GEN nf,long bound)
1939: {
1940: return ideallistzstarall(nf,bound,1);
1941: }
1942:
1943: GEN
1944: ideallistunit(GEN nf,long bound)
1945: {
1946: return ideallistzstarall(nf,bound,2);
1947: }
1948:
1949: GEN
1950: ideallistunitgen(GEN nf,long bound)
1951: {
1952: return ideallistzstarall(nf,bound,3);
1953: }
1954:
1955: GEN
1956: ideallist(GEN bnf,long bound)
1957: {
1958: return ideallistzstarall(bnf,bound,4);
1959: }
1960:
1961: static GEN
1962: ideallist_arch(GEN nf,GEN list,GEN arch,long flun)
1963: {
1964: long nba,i,j,lx,ly;
1965: GEN p1,z,p2;
1966:
1967: nba=0; for (i=1; i<lg(arch); i++) if (signe(arch[i])) nba++;
1968: lx=lg(list); z=cgetg(lx,t_VEC);
1969: for (i=1; i<lx; i++)
1970: {
1971: p2=(GEN)list[i]; checkbid(p2); ly=lg(p2);
1972: p1=cgetg(ly,t_VEC); z[i]=(long)p1;
1973: for (j=1; j<ly; j++)
1974: p1[j]=(long)zidealstarinitjoinarchall(nf,(GEN)p2[j],arch,nba,flun);
1975: }
1976: return z;
1977: }
1978:
1979: static GEN
1980: ideallistarchall(GEN bnf,GEN list,GEN arch,long flag)
1981: {
1982: long av,tetpil,i,j,lp1;
1983: long do_units = flag & 2;
1984: GEN nf = checknf(bnf), p1,p2,p3,racunit,funits,lu2,lu,embunit,z,y;
1985:
1986: if (typ(list) != t_VEC || (do_units && lg(list) != 3))
1987: err(typeer, "ideallistarch");
1988: if (lg(list) == 1) return cgetg(1,t_VEC);
1989: if (do_units)
1990: {
1991: y = cgetg(3,t_VEC);
1992: z = (GEN)list[1];
1993: lu= (GEN)list[2]; if (typ(lu) != t_VEC) err(typeer, "ideallistarch");
1994: }
1995: else z = list;
1996: if (typ(z) != t_VEC) err(typeer, "ideallistarch");
1997: z = ideallist_arch(nf,z,arch, flag & 1);
1998: if (!do_units) return z;
1999: y[1]=(long)z; av=avma;
2000: init_units(bnf,&funits,&racunit);
2001: lu2=cgetg(lg(z),t_VEC);
2002: for (i=1; i<lg(z); i++)
2003: {
2004: p1=(GEN)z[i]; p2=(GEN)lu[i]; lp1=lg(p1);
2005: p3=cgetg(lp1,t_VEC); lu2[i]=(long)p3;
2006: for (j=1; j<lp1; j++)
2007: {
2008: embunit = logunitmatrixarch(nf,funits,racunit,(GEN)p1[j]);
2009: p3[j] = lmul(gmael(p1,j,5), vconcat((GEN)p2[j],embunit));
2010: }
2011: }
2012: tetpil=avma; y[2]=lpile(av,tetpil,gcopy(lu2)); return y;
2013: }
2014:
2015: GEN
2016: ideallistarch(GEN nf, GEN list, GEN arch)
2017: {
2018: return ideallistarchall(nf,list,arch,0);
2019: }
2020:
2021: GEN
2022: ideallistarchgen(GEN nf, GEN list, GEN arch)
2023: {
2024: return ideallistarchall(nf,list,arch,1);
2025: }
2026:
2027: GEN
2028: ideallistunitarch(GEN bnf,GEN list,GEN arch)
2029: {
2030: return ideallistarchall(bnf,list,arch,2);
2031: }
2032:
2033: GEN
2034: ideallistunitarchgen(GEN bnf,GEN list,GEN arch)
2035: {
2036: return ideallistarchall(bnf,list,arch,3);
2037: }
2038:
2039: GEN
2040: ideallistarch0(GEN nf, GEN list, GEN arch,long flag)
2041: {
2042: if (!arch) arch=cgetg(1,t_VEC);
2043: if (flag<0 || flag>3) err(flagerr,"ideallistarch");
2044: return ideallistarchall(nf,list,arch,flag);
2045: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to base3.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari / src / basemath