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