Annotation of OpenXM_contrib/pari/src/basemath/buch3.c, Revision 1.1.1.1
1.1 maekawa 1: /*******************************************************************/
2: /* */
3: /* RAY CLASS FIELDS */
4: /* */
5: /*******************************************************************/
6: /* $Id: buch3.c,v 1.2 1999/09/23 17:50:56 karim Exp $ */
7: #include "pari.h"
8: #include "parinf.h"
9:
10: GEN compute_class_number(GEN mit,GEN *met,GEN *u1,GEN *u2);
11: GEN logunitmatrix(GEN nf,GEN funits,GEN racunit,GEN bid);
12: GEN vconcat(GEN Q1, GEN Q2);
13: GEN ideleaddone_aux(GEN nf,GEN x,GEN ideal);
14:
15: static GEN
16: get_full_rank(GEN nf, GEN v, GEN _0, GEN _1, GEN vecsign, GEN gen,
17: long ngen, long rankmax)
18: {
19: GEN v1,p1,alpha, bas=(GEN)nf[7], rac=(GEN)nf[6];
20: long rankinit=lg(v)-1, N=lgef(nf[1])-3, va=varn(nf[1]);
21: long limr,i,k,kk,r,rr;
22:
23: for (r=1,rr=3; ; r++,rr+=2)
24: {
25: p1 = gpuigs(stoi(rr),N);
26: limr=(cmpis(p1,BIGINT)>1)? BIGINT: p1[2]; /* min(BIGINT,rr^N) */
27: limr = (limr-1)>>1;
28: for (k=rr; k<=limr; k++)
29: {
30: long av1=avma;
31: alpha = gzero;
32: for (kk=k,i=1; i<=N; i++,kk/=rr)
33: {
34: long lambda = (kk+r)%rr - r;
35: if (lambda)
36: alpha = gadd(alpha,gmulsg(lambda,(GEN)bas[i]));
37: }
38: for (i=1; i<=rankmax; i++)
39: vecsign[i] = (gsigne(gsubst(alpha,va,(GEN)rac[i])) > 0)? (long)_0
40: : (long)_1;
41: v1 = concatsp(v, vecsign);
42: if (rank(v1) == rankinit) avma=av1;
43: else
44: {
45: v=v1; rankinit++; ngen++;
46: gen[ngen] = (long) alpha;
47: if (rankinit == rankmax) return ginv(v); /* full rank */
48: }
49: }
50: }
51: }
52:
53: GEN
54: buchnarrow(GEN bnf)
55: {
56: GEN nf,_0mod2,_1mod2,cyc,gen,v,matsign,arch;
57: GEN dataunit,p1,p2,p3,h,vecsign,clh,basecl,met,u1,u2;
58: long R1,R,i,j,ngen,sizeh,t,lo,c;
59: long av=avma,tetpil;
60:
61: if (typ(bnf)!=t_VEC || lg(bnf)!=11)
62: err(talker,"not a big number field vector in buchnarrow");
63: nf=checknf(bnf); R1=itos(gmael(nf,2,1));
64: if (!R1) return gcopy(gmael(bnf,8,1));
65:
66: _1mod2=gmodulss(1,2); _0mod2=gmodulss(0,2);
67: v=cgetg(R1+1,t_COL); vecsign=cgetg(R1+1,t_COL);
68: for (i=1; i<=R1; i++) v[i]=(long)_1mod2;
69: cyc=gmael3(bnf,8,1,2); gen=gmael3(bnf,8,1,3); ngen=lg(gen)-1;
70: matsign=signunits(bnf); R=lg(matsign); dataunit=cgetg(R+1,t_MAT);
71: for (j=1; j<R; j++)
72: {
73: p1=cgetg(R1+1,t_COL); dataunit[j]=(long)p1;
74: for (i=1; i<=R1; i++)
75: p1[i] = (signe(gcoeff(matsign,i,j)) > 0)? (long)_0mod2: (long)_1mod2;
76: }
77: dataunit[R]=(long)v; v=image(dataunit); t=lg(v)-1;
78: sizeh=ngen+R1-t; p1 = cgetg(sizeh+1,t_COL);
79: for (i=1; i<=ngen; i++) p1[i]=gen[i];
80: gen = p1;
81: if (t<R1)
82: v = get_full_rank(nf,v,_0mod2,_1mod2,vecsign,gen,ngen,R1);
83:
84: h=cgetg(sizeh+1,t_MAT); arch = cgetg(R1+1,t_VEC);
85: for (i=1; i<=R1; i++) arch[i]=un;
86: for (j=1; j<=ngen; j++)
87: {
88: p1 = cgetg(sizeh+1,t_COL); h[j]=(long)p1;
89: p2 = idealpow(nf, (GEN)gen[j], (GEN)cyc[j]);
90: p2 = (GEN)isprincipalall(bnf,p2,nf_GEN | nf_FORCE)[2];
91: p2 = gmul(v,zsigne(nf,p2,arch));
92: for (i=1; i<=ngen; i++) p1[i] = (i==j)? cyc[j]: zero;
93: for ( ; i<=sizeh; i++) p1[i] = llift((GEN)p2[i+t-ngen]);
94: }
95: for ( ; j<=sizeh; j++)
96: {
97: p1 = cgetg(sizeh+1,t_COL); h[j]=(long)p1;
98: for (i=1; i<=sizeh; i++) p1[i]=(i==j)?deux:zero;
99: }
100: clh=compute_class_number(h,&met,&u1,&u2);
101: u1=reducemodmatrix(u1,h); lo=lg(met)-1; c=0;
102: for (c=1; c<=lo; c++)
103: if (gcmp1(gcoeff(met,c,c))) break;
104: basecl=cgetg(c,t_VEC);
105: for (j=1; j<c; j++)
106: {
107: p1=gcoeff(u1,1,j);
108: p3=idealpow(nf,(GEN)gen[1],p1);
109: if (signe(p1)<0) p3=numer(p3);
110: for (i=2; i<=lo; i++)
111: {
112: p1=gcoeff(u1,i,j);
113: if (signe(p1))
114: {
115: p3 = idealmul(nf,p3, idealpow(nf,(GEN)gen[i],p1));
116: p1 = content(p3); if (!gcmp1(p1)) p3 = gdiv(p3,p1);
117: }
118: }
119: basecl[j]=(long)p3;
120: }
121: tetpil=avma; v=cgetg(4,t_VEC);
122: v[1]=lcopy(clh); setlg(met,c);
123: v[2]=(long)mattodiagonal(met);
124: v[3]=lcopy(basecl); return gerepile(av,tetpil,v);
125: }
126:
127: GEN idealaddtoone_i(GEN nf, GEN x, GEN y);
128:
129: /* given two coprime ideals x (integral) and id, compute alpha in x,
130: * alpha = 1 mod (id), with x/alpha nearly reduced.
131: */
132: static GEN
133: findalpha(GEN nf,GEN x,GEN id,long prec)
134: {
135: GEN p1,idprod,y;
136: GEN alp = idealaddtoone_i(nf,x,id);
137:
138: idprod = idealmullll(nf,x,id);
139: y = lllgram(qf_base_change(gmael(nf,5,3),idprod,1), 2*prec-2);
140: y = gmul(idprod, (GEN)y[1]); /* small vector in idprod */
141:
142: p1 = ground(element_div(nf,alp,y));
143: alp = gsub(alp,element_mul(nf,p1,y));
144: return gcmp0(alp)? y: alp;
145: }
146:
147: static int
148: too_big(GEN nf, GEN bet)
149: {
150: GEN x = gnorm(basistoalg(nf,bet));
151: switch (typ(x))
152: {
153: case t_INT: return absi_cmp(x, gun);
154: case t_FRAC: return absi_cmp((GEN)x[1], (GEN)x[2]);
155: }
156: err(bugparier, "wrong type in too_big");
157: return 0; /* not reached */
158: }
159:
160: static GEN
161: idealmodidele(GEN nf, GEN x, GEN ideal, GEN sarch, GEN arch, long prec)
162: {
163: long av = avma,i,l;
164: GEN p1,p2,alp,bet,b;
165:
166: nf=checknf(nf); alp=findalpha(nf,x,ideal,prec);
167: p1=idealdiv(nf,alp,x);
168: bet = element_div(nf,findalpha(nf,p1,ideal,prec),alp);
169: if (too_big(nf,bet) > 0) { avma=av; return x; }
170: p1=(GEN)sarch[2]; l=lg(p1);
171: if (l > 1)
172: {
173: b=bet; p2=lift_intern(gmul((GEN)sarch[3],zsigne(nf,bet,arch)));
174: for (i=1; i<l; i++)
175: if (signe(p2[i])) bet = element_mul(nf,bet,(GEN)p1[i]);
176: if (b != bet && too_big(nf,bet) > 0) { avma=av; return x; }
177: }
178: return idealmul(nf,bet,x);
179: }
180:
181: static GEN
182: idealmulmodidele(GEN nf,GEN x,GEN y, GEN ideal,GEN sarch,GEN arch,long prec)
183: {
184: return idealmodidele(nf,idealmul(nf,x,y),ideal,sarch,arch,prec);
185: }
186:
187: /* assume n > 0 */
188: static GEN
189: idealpowmodidele(GEN nf,GEN x,GEN n, GEN ideal,GEN sarch,GEN arch,long prec)
190: {
191: long i,m,av=avma;
192: GEN y;
193: ulong j;
194:
195: if (cmpis(n, 16) < 0)
196: {
197: if (gcmp1(n)) return x;
198: x = idealpow(nf,x,n);
199: x = idealmodidele(nf,x,ideal,sarch,arch,prec);
200: return gerepileupto(av,x);
201: }
202:
203: i = lgefint(n)-1; m=n[i]; j=HIGHBIT;
204: while ((m&j)==0) j>>=1;
205: y = x;
206: for (j>>=1; j; j>>=1)
207: {
208: y = idealmul(nf,y,y);
209: if (m&j) y = idealmul(nf,y,x);
210: y = idealmodidele(nf,y,ideal,sarch,arch,prec);
211: }
212: for (i--; i>=2; i--)
213: for (m=n[i],j=HIGHBIT; j; j>>=1)
214: {
215: y = idealmul(nf,y,y);
216: if (m&j) y = idealmul(nf,y,x);
217: y = idealmodidele(nf,y,ideal,sarch,arch,prec);
218: }
219: return gerepileupto(av,y);
220: }
221:
222: static GEN
223: buchrayall(GEN bnf,GEN module,long flag,long prec)
224: {
225: GEN nf,cyc,gen,genplus,fa2,sarch,hmatu,u,clg;
226: GEN dataunit,p1,p2,h,clh,basecl,met,u1,u2,u1old;
227: GEN racunit,bigres,bid,resbid2,resbid3,x,y,funits,hmat,vecel;
228: long RU,R3,i,j,ngen,lh,lo,c,av=avma,N;
229:
230: bnf = checkbnf(bnf); nf=checknf(bnf); bigres=(GEN)bnf[8];
231: funits = check_units(bnf, "buchrayall");
232: N=lgef(nf[1])-3;
233: cyc=gmael(bigres,1,2);
234: gen=gmael(bigres,1,3); ngen=lg(cyc)-1;
235:
236: bid = zidealstarinitall(nf,module,1);
237: resbid2=gmael(bid,2,2);
238: resbid3=gmael(bid,2,3);
239: R3=lg(resbid2)-1; lh=ngen+R3;
240:
241: x = idealhermite(nf,module);
242: if (R3 || flag & (nf_INIT|nf_GEN))
243: {
244: vecel=cgetg(ngen+1,t_VEC);
245: for (j=1; j<=ngen; j++)
246: vecel[j]=(long)idealcoprime(nf,(GEN)gen[j],x);
247: }
248: if (flag & nf_GEN)
249: {
250: genplus=cgetg(lh+1,t_VEC);
251: for (j=1; j<=ngen; j++)
252: genplus[j] = (long) idealmul(nf,(GEN)vecel[j],(GEN)gen[j]);
253: for ( ; j<=lh; j++)
254: genplus[j] = resbid3[j-ngen];
255: }
256: if (!R3)
257: {
258: if (!(flag & nf_GEN)) clg=cgetg(3,t_VEC);
259: else
260: { clg=cgetg(4,t_VEC); clg[3]=(long)genplus; }
261: clg[1]=mael(bigres,1,1);
262: clg[2]=(long)cyc;
263: if (!(flag & nf_INIT)) return gerepileupto(av,gcopy(clg));
264: y = cgetg(7,t_VEC);
265: y[1]=lcopy(bnf);
266: y[2]=lcopy(bid);
267: y[3]=lcopy(vecel);
268: y[4]=(long)idmat(ngen);
269: y[5]=lcopy(clg); u=cgetg(3,t_VEC);
270: y[6]=(long)u;
271: u[1]=lgetg(1,t_MAT);
272: u[2]=(long)idmat(lg(funits));
273: return gerepileupto(av,y);
274: }
275: fa2=(GEN)bid[4]; sarch=(GEN)fa2[lg(fa2)-1];
276:
277: RU=lg(funits); dataunit=cgetg(RU+R3+1,t_MAT);
278: racunit=gmael(bigres,4,2);
279: dataunit[1] = (long)zideallog(nf,racunit,bid);
280: for (j=2; j<=RU; j++)
281: dataunit[j] = (long)zideallog(nf,(GEN)funits[j-1],bid);
282: for ( ; j<=RU+R3; j++)
283: {
284: p1=cgetg(R3+1,t_COL); dataunit[j]=(long)p1;
285: for (i=1; i<=R3; i++)
286: p1[i] = (i==(j-RU))? resbid2[i]: zero;
287: }
288: h=cgetg(lh+1,t_MAT);
289: for (j=1; j<=ngen; j++)
290: {
291: p1=cgetg(lh+1,t_COL); h[j]=(long)p1;
292: p2 = idealpow(nf, (GEN)gen[j], (GEN)cyc[j]);
293: p2 = (GEN)isprincipalall(bnf,p2,nf_GEN | nf_FORCE)[2];
294: p2 = element_mul(nf,p2,element_pow(nf,(GEN)vecel[j],(GEN)cyc[j]));
295: p2 = zideallog(nf,p2,bid);
296: for (i=1; i<=ngen; i++) p1[i] = (i==j)? cyc[j]: zero;
297: for ( ; i<=lh; i++) p1[i] = lnegi((GEN)p2[i-ngen]);
298: }
299:
300: hmatu=hnfall(dataunit); hmat=(GEN)hmatu[1];
301: for ( ; j<=lh; j++)
302: {
303: p1=cgetg(lh+1,t_COL); h[j]=(long)p1;
304: for (i=1; i<=ngen; i++) p1[i]=zero;
305: for ( ; i<=lh; i++) p1[i]=coeff(hmat,i-ngen,j-ngen);
306: }
307: clh = compute_class_number(h,&met,&u1,&u2);
308: u1old=u1; lo=lg(met)-1;
309: for (c=1; c<=lo; c++)
310: if (gcmp1(gcoeff(met,c,c))) break;
311:
312: if (flag & nf_GEN)
313: {
314: GEN Id=idmat(N), arch=(GEN)module[2];
315: u1 = reducemodmatrix(u1,h);
316: basecl=cgetg(c,t_VEC);
317: for (j=1; j<c; j++)
318: {
319: GEN *op, minus = Id, plus = Id;
320: long av1 = avma, s;
321: for (i=1; i<=lo; i++)
322: {
323: p1 = gcoeff(u1,i,j);
324: if (!(s = signe(p1))) continue;
325:
326: if (s > 0) op = + else { op = − p1 = negi(p1); }
327: p1 = idealpowmodidele(nf,(GEN)genplus[i],p1,x,sarch,arch,prec);
328: *op = *op==Id? p1
329: : idealmulmodidele(nf,*op,p1,x,sarch,arch,prec);
330: }
331: if (minus == Id) p1 = plus;
332: else
333: {
334: p2 = ideleaddone_aux(nf,minus,module);
335: p1 = idealdivexact(nf,idealmul(nf,p2,plus),minus);
336: p1 = idealmodidele(nf,p1,x,sarch,arch,prec);
337: }
338: basecl[j]=lpileupto(av1,p1);
339: }
340: clg=cgetg(4,t_VEC); clg[3]=lcopy(basecl);
341: } else clg=cgetg(3,t_VEC);
342: clg[1]=licopy(clh); setlg(met,c);
343: clg[2]=(long)mattodiagonal(met);
344: if (!(flag & nf_INIT)) return gerepileupto(av,clg);
345:
346: u2 = cgetg(R3+1,t_MAT);
347: u1 = cgetg(RU+1,t_MAT); u = (GEN)hmatu[2];
348: for (j=1; j<=RU; j++) { u1[j]=u[j]; setlg(u[j],RU+1); }
349: u += RU;
350: for (j=1; j<=R3; j++) { u2[j]=u[j]; setlg(u[j],RU+1); }
351: p1=lllint(u1); p2=ginv(hmat);
352: y=cgetg(7,t_VEC);
353: y[1]=lcopy(bnf);
354: y[2]=lcopy(bid);
355: y[3]=lcopy(vecel);
356: y[4]=linv(u1old);
357: y[5]=lcopy(clg); u=cgetg(3,t_VEC);
358: y[6]=(long)u;
359: u[1]=lmul(u2,p2);
360: u[2]=lmul(u1,p1);
361: return gerepileupto(av,y);
362: }
363:
364: GEN
365: buchrayinitgen(GEN bignf, GEN ideal,long prec)
366: {
367: return buchrayall(bignf,ideal, nf_INIT | nf_GEN,prec);
368: }
369:
370: GEN
371: buchrayinit(GEN bignf, GEN ideal,long prec)
372: {
373: return buchrayall(bignf,ideal, nf_INIT,prec);
374: }
375:
376: GEN
377: buchray(GEN bignf, GEN ideal,long prec)
378: {
379: return buchrayall(bignf,ideal, nf_GEN,prec);
380: }
381:
382: GEN
383: bnrclass0(GEN bignf, GEN ideal, long flag, long prec)
384: {
385: switch(flag)
386: {
387: case 0: flag = nf_GEN; break;
388: case 1: flag = nf_INIT; break;
389: case 2: flag = nf_INIT | nf_GEN; break;
390: default: err(flagerr,"bnrclass");
391: }
392: return buchrayall(bignf,ideal,flag,prec);
393: }
394:
395: GEN
396: bnrinit0(GEN bignf, GEN ideal, long flag, long prec)
397: {
398: switch(flag)
399: {
400: case 0: flag = nf_INIT; break;
401: case 1: flag = nf_INIT | nf_GEN; break;
402: default: err(flagerr,"bnrinit");
403: }
404: return buchrayall(bignf,ideal,flag,prec);
405: }
406:
407: GEN
408: rayclassno(GEN bnf,GEN ideal)
409: {
410: GEN nf,clno,dataunit,p1;
411: GEN racunit,bigres,bid,resbid,resbid2,funits,hmat;
412: long RU,R3,i,j,av=avma;
413:
414: bnf = checkbnf(bnf); nf=(GEN)bnf[7]; bigres=(GEN)bnf[8];
415: funits = check_units(bnf,"rayclassno");
416: clno = gmael(bigres,1,1);
417: bid = zidealstarinitall(nf,ideal,0);
418: resbid=(GEN)bid[2]; resbid2=(GEN)resbid[2];
419: R3=lg(resbid2)-1; if (!R3) { avma=av; return icopy(clno); }
420:
421: RU=lg(funits); dataunit=cgetg(RU+R3+1,t_MAT);
422: racunit=gmael(bigres,4,2); dataunit[1]=(long)zideallog(nf,racunit,bid);
423: for (j=2; j<=RU; j++)
424: dataunit[j]=(long)zideallog(nf,(GEN)funits[j-1],bid);
425: for ( ; j<=RU+R3; j++)
426: {
427: p1=cgetg(R3+1,t_COL); dataunit[j]=(long)p1;
428: for (i=1; i<=R3; i++)
429: p1[i] = (i==(j-RU))?resbid2[i]:zero;
430: }
431: hmat=hnfmod(dataunit,(GEN)resbid[1]);
432: for (i=lg(hmat)-1 ; i; i--) clno = mulii(clno,gcoeff(hmat,i,i));
433: avma=av; return icopy(clno);
434: }
435:
436: GEN
437: isprincipalrayall(GEN bnr, GEN x, long flag)
438: {
439: long av=avma,tetpil,i,j,c,N,ngen,ngzk;
440: GEN bnf,nf,bid,vecel,vecep,matu,ep,p1,p2,p3,p4,beta,idep,y,rayclass;
441: GEN divray,genray,alpha,alphaall,racunit,res,funit,pol;
442:
443: checkbnr(bnr); bnf=(GEN)bnr[1]; bid=(GEN)bnr[2];
444: vecel=(GEN)bnr[3]; matu=(GEN)bnr[4];
445: rayclass=(GEN)bnr[5]; nf=(GEN)bnf[7]; ngen=lg(vecel)-1;
446: idep=isprincipalall(bnf,x,nf_GEN | nf_FORCE);
447: if (lg(idep[1]) != ngen+1)
448: err(talker,"incorrect generator length in isprincipalray");
449: pol=(GEN)nf[1]; N=lgef(pol)-3;
450: p2=cgetg(N+1,t_COL); p2[1]=un;
451: for (i=2; i<=N; i++) p2[i]=zero;
452: ep=(GEN)idep[1];
453: for (i=1; i<=ngen; i++)
454: if (typ(vecel[i]) != t_INT)
455: p2=element_mul(nf,p2,element_pow(nf,(GEN)vecel[i],(GEN)ep[i]));
456: beta=element_div(nf,(GEN)idep[2],p2);
457: p3=zideallog(nf,beta,bid); ngzk=lg(p3)-1;
458: vecep=cgetg(ngen+ngzk+1,t_COL);
459: for (i=1; i<=ngen; i++) vecep[i]=ep[i];
460: for ( ; i<=ngen+ngzk; i++) vecep[i]=p3[i-ngen];
461: p1=gmul(matu,vecep);
462: divray=(GEN)rayclass[2]; c=lg(divray);
463: tetpil=avma; y=cgetg(c,t_COL);
464: for (i=1; i<c; i++)
465: y[i] = lmodii((GEN)p1[i],(GEN)divray[i]);
466: if (!(flag & nf_GEN)) return gerepile(av,tetpil,y);
467:
468: if (lg(rayclass)<=3)
469: err(talker,"please apply bnrinit(,,1) and not bnrinit(,,0)");
470:
471: genray=(GEN)rayclass[3]; p1=idmat(N);
472: for (i=1; i<c; i++)
473: p1=idealmul(nf,idealpow(nf,(GEN)genray[i],(GEN)y[i]),p1);
474: alphaall = isprincipalall(bnf,idealdiv(nf,x,p1),nf_GEN | nf_FORCE);
475: if (!gcmp0((GEN)alphaall[1])) err(bugparier,"isprincipalray (bug1)");
476:
477: res=(GEN)bnf[8];
478: funit = check_units(bnf,"isprincipalrayall");
479: racunit=(GEN)res[4];
480: alpha = basistoalg(nf,(GEN)alphaall[2]);
481: p3=zideallog(nf,(GEN)alphaall[2],bid);
482: if (lg(p3)>1)
483: {
484: p4=(GEN)bnr[6]; p3=gmul((GEN)p4[1],p3);
485: if (!gcmp1(denom(p3))) err(bugparier,"isprincipalray (bug2)");
486:
487: x=lllreducemodmatrix(p3,(GEN)p4[2]);
488: p3=gpui(gmodulcp((GEN)racunit[2],pol),(GEN)x[1],0);
489: for (j=1; j<lg(funit); j++)
490: p3=gmul(p3,gpui(gmodulcp((GEN)funit[j],pol),(GEN)x[j+1],0));
491: alpha = gdiv(alpha,p3);
492: }
493: tetpil=avma; p1=cgetg(4,t_VEC);
494: p1[1]=lcopy(y); p1[2]=(long)algtobasis(nf,alpha);
495: p1[3]=lmin((GEN)idep[3],(GEN)alphaall[3]);
496: return gerepile(av,tetpil,p1);
497: }
498:
499: GEN
500: isprincipalray(GEN bignfray, GEN x)
501: {
502: return isprincipalrayall(bignfray,x,nf_REGULAR);
503: }
504:
505: GEN
506: isprincipalraygen(GEN bignfray, GEN x)
507: {
508: return isprincipalrayall(bignfray,x,nf_GEN);
509: }
510:
511: /* DK = |dK|; give c[N][R1] */
512: static GEN
513: zimmertbound(long N,long R1,GEN DK)
514: {
515: long av,tetpil,i,R2;
516: GEN w,p1,minkowski;
517:
518: if (N<2) return gun;
519: av=avma;
520: if (N<21)
521: {
522: double **c=(double**)gpmalloc(sizeof(double*)*21);
523: for (i=1; i<=20; i++) c[i]=(double*)gpmalloc(sizeof(double)*21);
524: c[2][2] = -0.6931; c[3][3] = -1.71733859;
525: c[2][0] = -0.45158; c[3][1] = -1.37420604;
526:
527: c[4][4] = -2.91799837; c[5][5] = -4.22701425;
528: c[4][2] = -2.50091538; c[5][3] = -3.75471588;
529: c[4][0] = -2.11943331; c[5][1] = -3.31196660;
530:
531: c[6][6] = -5.61209925; c[7][7] = -7.05406203;
532: c[6][4] = -5.09730381; c[7][5] = -6.50550021;
533: c[6][2] = -4.60693851; c[7][3] = -5.97735406;
534: c[6][0] = -4.14303665; c[7][1] = -5.47145968;
535:
536: c[8][8] = -8.54052636; c[9][9] = -10.0630022;
537: c[8][6] = -7.96438858; c[9][7] = -9.46382812;
538: c[8][4] = -7.40555445; c[9][5] = -8.87952524;
539: c[8][2] = -6.86558259; c[9][3] = -8.31139202;
540: c[8][0] = -6.34608077; c[9][1] = -7.76081149;
541:
542: c[10][10]= -11.6153797; c[11][11]= -13.1930961;
543: c[10][8] = -10.9966020; c[11][9] = -12.5573772;
544: c[10][6] = -10.3907654; c[11][7] = -11.9330458;
545: c[10][4] = -9.79895170; c[11][5] = -11.3210061;
546: c[10][2] = -9.22232770; c[11][3] = -10.7222412;
547: c[10][0] = -8.66213267; c[11][1] = -10.1378082;
548:
549: c[12][12]= -14.7926394; c[13][13]= -16.4112395;
550: c[12][10]= -14.1420915; c[13][11]= -15.7475710;
551: c[12][8] = -13.5016616; c[13][9] = -15.0929680;
552: c[12][6] = -12.8721114; c[13][7] = -14.4480777;
553: c[12][4] = -12.2542699; c[13][5] = -13.8136054;
554: c[12][2] = -11.6490374; c[13][3] = -13.1903162;
555: c[12][0] = -11.0573775; c[13][1] = -12.5790381;
556:
557: c[14][14]= -18.0466672; c[15][15]= -19.6970961;
558: c[14][12]= -17.3712806; c[15][13]= -19.0111606;
559: c[14][10]= -16.7040780; c[15][11]= -18.3326615;
560: c[14][8] = -16.0456127; c[15][9] = -17.6620757;
561: c[14][6] = -15.3964878; c[15][7] = -16.9999233;
562: c[14][4] = -14.7573587; c[15][5] = -16.3467686;
563: c[14][2] = -14.1289364; c[15][3] = -15.7032228;
564: c[14][0] = -13.5119848; c[15][1] = -15.0699480;
565:
566: c[16][16]= -21.3610081; c[17][17]= -23.0371259;
567: c[16][14]= -20.6655103; c[17][15]= -22.3329066;
568: c[16][12]= -19.9768082; c[17][13]= -21.6349299;
569: c[16][10]= -19.2953176; c[17][11]= -20.9435607;
570: c[16][8] = -18.6214885; c[17][9] = -20.2591899;
571: c[16][6] = -17.9558093; c[17][7] = -19.5822454;
572: c[16][4] = -17.2988108; c[17][5] = -18.9131878;
573: c[16][2] = -16.6510652; c[17][3] = -18.2525157;
574: c[16][0] = -16.0131906; c[17][1] = -17.6007672;
575:
576: c[18][18]= -24.7243611; c[19][19]= -26.4217792;
577: c[18][16]= -24.0121449; c[19][17]= -25.7021950;
578: c[18][14]= -23.3056902; c[19][15]= -24.9879497;
579: c[18][12]= -22.6053167; c[19][13]= -24.2793271;
580: c[18][10]= -21.9113705; c[19][11]= -23.5766321;
581: c[18][8] = -21.2242247; c[19][9] = -22.8801952;
582: c[18][6] = -20.5442836; c[19][7] = -22.1903709;
583: c[18][4] = -19.8719830; c[19][5] = -21.5075437;
584: c[18][2] = -19.2077941; c[19][3] = -20.8321263;
585: c[18][0] = -18.5522234; c[19][1] = -20.1645647;
586:
587: c[20][20]= -28.1285704;
588: c[20][18]= -27.4021674;
589: c[20][16]= -26.6807314;
590: c[20][14]= -25.9645140;
591: c[20][12]= -25.2537867;
592: c[20][10]= -24.5488420;
593: c[20][8] = -23.8499943;
594: c[20][6] = -23.1575823;
595: c[20][4] = -22.4719720;
596: c[20][2] = -21.7935548;
597: c[20][0] = -21.1227537;
598: w=gexp(dbltor(c[N][R1]),6);
599: for (i=1; i<=20; i++) free(c[i]); free(c);
600: p1=gmul(gsqrt(DK,MEDDEFAULTPREC),w);
601: tetpil=avma; return gerepile(av,tetpil,ground(p1));
602: }
603: R2=(N-R1)>>1; p1=gdiv(mpfact(N),gpuigs(stoi(N),N));
604: minkowski=ground(gmul(gmul(p1,gpuigs(gdivsg(4,mppi(MEDDEFAULTPREC)),R2)),gsqrt(DK,MEDDEFAULTPREC)));
605: if (cmpis(minkowski,500000)>0)
606: err(talker,"The field has degree more than 20 and the Minkowski bound is larger than 500 000: it is unrealistic to certify it");
607:
608: tetpil=avma; return gerepile(av,tetpil,gcopy(minkowski));
609: }
610:
611: /* all primes up to Minkowski bound factor on factorbase ? */
612: static void
613: testprime(GEN bnf,GEN minkowski)
614: {
615: long av = avma, pp,i,nbideal,k,pmax;
616: GEN f,p1,vectpp,fb,dK, nf=checknf(bnf);
617: byteptr delta = diffptr;
618:
619: if (DEBUGLEVEL>=2)
620: fprintferr("PHASE 1: check primes to Zimmert bound = %Z\n\n",minkowski);
621: f=(GEN)nf[4]; dK=(GEN)nf[3];
622: if (!gcmp1(f))
623: {
624: GEN different = gmael(nf,5,5);
625: if (DEBUGLEVEL>=2)
626: fprintferr("**** Testing Different = %Z\n",different);
627: p1 = isprincipalall(bnf,different,nf_FORCE);
628: if (DEBUGLEVEL>=2) fprintferr(" is %Z\n",p1);
629: }
630: fb=(GEN)bnf[5];
631: p1 = gmael(fb, lg(fb)-1, 1); /* largest p in factorbase */
632: pp = 0; pmax = is_bigint(p1)? VERYBIGINT: itos(p1);
633: while (cmpsi(pp,minkowski)<1)
634: {
635: pp += *delta++; if (!*delta) err(primer1);
636: if (DEBUGLEVEL>=2) fprintferr("*** p = %ld\n",pp);
637: vectpp=primedec(bnf,stoi(pp)); nbideal=lg(vectpp)-1;
638: /* loop through all P | p if ramified, all but one otherwise */
639: if (!smodis(dK,pp)) nbideal++;
640: for (i=1; i<nbideal; i++)
641: {
642: GEN P = (GEN)vectpp[i]; /* non inert */
643: if (DEBUGLEVEL>=2)
644: fprintferr(" Testing P = %Z\n",P);
645: if (cmpii(idealnorm(bnf,P),minkowski) < 1)
646: {
647: if (pp <= pmax && (k = tablesearch(fb, P, cmp_prime_ideal)))
648: {
649: if (DEBUGLEVEL>=2) fprintferr(" #%ld in factor base\n",k);
650: }
651: else
652: {
653: p1=isprincipalall(bnf,P,nf_FORCE);
654: if (DEBUGLEVEL>=2) fprintferr(" is %Z\n",p1);
655: }
656: }
657: else if (DEBUGLEVEL>=2)
658: fprintferr(" Norm(P) > Zimmert bound\n");
659: }
660: }
661: avma=av;
662: if (DEBUGLEVEL>=2) { fprintferr("End of PHASE 1.\n\n"); flusherr(); }
663: }
664:
665: /* rend constante d'Hermite^n si connue, sinon une borne sup */
666: static GEN
667: hermiteconstant(long n)
668: {
669: long av,tetpil;
670: GEN h,h1;
671:
672: switch(n)
673: {
674: case 1: return gun;
675: case 2: h=cgetg(3,t_FRAC); h[1]=lstoi(4); h[2]=lstoi(3); return h;
676: case 3: return gdeux;
677: case 4: return stoi(4);
678: case 5: return stoi(8);
679: case 6: h=cgetg(3,t_FRAC); h[1]=lstoi(64); h[2]=lstoi(3); return h;
680: case 7: return stoi(64);
681: case 8: return stoi(256);
682: }
683: av = avma;
684: h = gpuigs(gdiv(gdeux,mppi(DEFAULTPREC)),n);
685: h1 = gsqr(ggamma(gdivgs(stoi(n+4),2),DEFAULTPREC));
686: tetpil=avma; return gerepile(av,tetpil,gmul(h,h1));
687: }
688:
689: /* 1 primitif, 0 s'il est peut etre imprimitif... */
690: static long
691: isprimitive(GEN nf)
692: {
693: long N,first,i,l,ep;
694: GEN d,fa;
695:
696: N = lgef(nf[1])-3; fa = (GEN)factor(stoi(N))[1]; /* primes | N */
697: first = itos((GEN)fa[1]); if (first==N) return 1;
698:
699: d=absi((GEN)nf[3]); fa=(GEN)factor(d)[2]; /* expo. primes | disc(nf) */
700: if (mod2(d))
701: { i=1; ep=1; }
702: else
703: { i=2; ep=itos((GEN)fa[1])>>1; }
704: l=lg(fa);
705: for ( ; i < l; i++)
706: {
707: if (ep >= first) return 0;
708: ep = itos((GEN)fa[i]);
709: }
710: return 1;
711: }
712:
713: static GEN
714: regulatorbound(GEN bnf)
715: {
716: long N,R1,R2,R;
717: GEN nf,dKa,bound,p1,c1;
718:
719: nf=(GEN)bnf[7]; N=lgef(nf[1])-3;
720: bound=dbltor(0.2);
721: if (!isprimitive(nf))
722: {
723: if (DEBUGLEVEL>=2)
724: { fprintferr("Default bound for regulator: 0.2\n"); flusherr(); }
725: return bound;
726: }
727: dKa=absi((GEN)nf[3]);
728: R1=itos(gmael(nf,2,1));
729: R2=itos(gmael(nf,2,2)); R=R1+R2-1;
730: if (!R2 && N<12) c1=gpuigs(stoi(4),N>>1); else c1=gpuigs(stoi(N),N);
731: if (cmpii(dKa,c1)<=0)
732: {
733: if (DEBUGLEVEL>=2)
734: { fprintferr("Default bound for regulator: 0.2\n"); flusherr(); }
735: return bound;
736: }
737: p1 = gsqr(glog(gdiv(dKa,c1),DEFAULTPREC));
738: p1 = gdivgs(gmul2n(gpuigs(gdivgs(gmulgs(p1,3),N*(N*N-1)-6*R2),R),R2),N);
739: p1 = gsqrt(gdiv(p1, hermiteconstant(R)), DEFAULTPREC);
740: if (gcmp(p1,bound) > 0) bound = p1;
741: if (DEBUGLEVEL>=2)
742: { fprintferr("Mahler bound for regulator: "); outerr(p1); flusherr(); }
743: return bound;
744: }
745:
746: #define NBMAX 5000
747: /* should use smallvectors */
748: static GEN
749: minimforunits(GEN nf, long borne, long stockmax)
750: {
751: long av = avma,av1,n1,n,i,j,k,s,norme,normax,*x,fl1,cmpt;
752: GEN u,r,S,S1,a,base,p1;
753: double p;
754: double **q,*v,*y,*z;
755: double eps=0.000001;
756:
757: if (DEBUGLEVEL>=2)
758: {
759: fprintferr("Searching minimum of T2-form on units:\n");
760: if (DEBUGLEVEL>2) fprintferr(" borne = %ld\n",borne);
761: flusherr();
762: }
763: a=gmael(nf,5,3); n1=lg(a);
764: n=n1-1;
765: x=(long*)gpmalloc(n1*sizeof(long));
766: y=(double*)gpmalloc(n1*sizeof(double));
767: z=(double*)gpmalloc(n1*sizeof(double));
768: v=(double*)gpmalloc(n1*sizeof(double));
769: q=(double**)gpmalloc(n1*sizeof(double*));
770: for (j=1; j<=n; j++) q[j]=(double*)gpmalloc(n1*sizeof(double));
771: u=lllgram(a,BIGDEFAULTPREC); base=gmul((GEN)nf[7],u);
772: a=gmul(qf_base_change(a,u,1), realun(BIGDEFAULTPREC));
773: r=sqred1(a);
774: for (j=1; j<=n; j++)
775: {
776: v[j]=rtodbl(gcoeff(r,j,j));
777: for (i=1; i<j; i++)
778: q[i][j]=rtodbl(gcoeff(r,i,j));
779: }
780: normax=0;
781: if (stockmax) S=cgetg(stockmax+1,t_MAT);
782: s=0; k=n; cmpt=0; y[n]=z[n]=0;
783: x[n]=(long)(sqrt(borne/v[n]+eps));
784:
785: for(;;)
786: {
787: do
788: {
789: if (k>1)
790: {
791: k--; z[k]=0;
792: for (j=k+1; j<=n; j++) z[k]=z[k]+q[k][j]*x[j];
793: p=x[k+1]+z[k+1];
794: y[k]=y[k+1]+p*p*v[k+1];
795: x[k]=(long)floor(sqrt((borne-y[k]+eps)/v[k])-z[k]);
796: }
797: while (v[k]*(x[k]+z[k])*(x[k]+z[k]) > borne-y[k]+eps)
798: {
799: k++; x[k]--;
800: }
801: }
802: while (k>1);
803: if (!x[1] && y[1]<=eps) break;
804:
805: cmpt++;
806: if (cmpt>NBMAX)
807: {
808: free(x); free(y); free(z); free(v);
809: for (j=1; j<=n; j++) free(q[j]); free(q);
810: av=avma; return NULL;
811: }
812: if (DEBUGLEVEL>8){ fprintferr("."); flusherr(); }
813: norme=(long)(y[1]+v[1]*(x[1]+z[1])*(x[1]+z[1])+eps);
814: if (norme>normax) normax=norme;
815: av1=avma; p1=gzero;
816: for (i=1; i<=n; i++) p1=gadd(p1,gmulsg(x[i],(GEN)base[i]));
817: fl1=gcmp1(gabs(subres(p1,(GEN)nf[1]),0)); avma=av1;
818: if (fl1)
819: {
820: if (DEBUGLEVEL>=2) { fprintferr("*"); flusherr(); }
821: s++; cmpt=0;
822: if (s<=stockmax)
823: {
824: p1=cgetg(n+1,t_COL);
825: for (i=1; i<=n; i++) p1[i]=lstoi(x[i]);
826: S[s]=lmul(u,p1);
827: }
828: }
829: x[k]--;
830: }
831: free(x); free(y); free(z); free(v);
832: for (j=1; j<=n; j++) free(q[j]); free(q);
833: if (DEBUGLEVEL>=2){ fprintferr("\n"); flusherr(); }
834: if (stockmax)
835: {
836: av1=avma;
837: k=(s<stockmax)? s:stockmax;
838: S1=cgetg(k+1,t_MAT);
839: for (j=1; j<=k; j++) S1[j]=lcopy((GEN)S[j]);
840: S=gerepile(av,av1,S1);
841: }
842: else { avma=av; S=cgetg(1,t_MAT); }
843: u=cgetg(4,t_VEC);
844: u[1]=lstoi(s<<1);
845: u[2]=lstoi(normax);
846: u[3]=(long)S;
847: return u;
848: }
849:
850: #undef NBMAX
851:
852: static GEN
853: compute_M0(GEN M_star,long N) /* On connait M_star; on calcule M0 */
854: {
855: long av1,tetpil,m1,m2,n1,n2,n3,k,kk,lr,lr1,lr2,i,j,l,vx,vy,vz,vM,PREC,prec;
856: GEN eps,pol,p1,p2,p3,p4,p5,p6,u,v,w,r,r1,r2,M0,M0_pro,S,P,M_formel;
857: GEN f1,f2,f3,g1,g2,g3,pg1,pg2,pg3,pf1,pf2,pf3,p7,p8,p9,p10,p11;
858: GEN x,y,z;
859:
860: PREC=gprecision(M_star);
861: if (N==2) return gmul2n(gsqr(gach(gmul2n(M_star,-1),PREC)), -1);
862: vM = fetch_var(); M_formel=polx[vM];
863: vz = fetch_var(); z=polx[vz];
864: vy = fetch_var(); y=polx[vy];
865: vx = fetch_var(); x=polx[vx];
866: vx=varn(x); vy=varn(y); vz=varn(z); vM=varn(M_formel);
867:
868: PREC = PREC>>1; if (!PREC) PREC=DEFAULTPREC;
869: eps=dbltor(0.0000001); M0=gzero; m1=(N-(N%3))/3;
870: for (n1=1; n1<=m1; n1++)
871: {
872: m2 = (N-n1)>>1;
873: for (n2=n1; n2<=m2; n2++)
874: {
875: av1=avma; n3=N-n1-n2; prec=PREC;
876: if (!(N%3) && n1==n2 && n1==n3)
877: {
878: p1=gdivgs(M_star,m1); p2=gaddsg(1,p1); p3=gsubgs(p1,3);
879: p4=gsqrt(gmul(p2,p3),prec); p5=gsubgs(p1,1);
880: u=gun; v=gmul2n(gadd(p5,p4),-1); w=gmul2n(gsub(p5,p4),-1);
881: M0_pro=gmul2n(gmulsg(m1,gadd(gsqr(glog(v,prec)),gsqr(glog(w,prec)))),-2);
882: if (DEBUGLEVEL>2)
883: {
884: fprintferr("[%ld,%ld,%ld]: ",n1,n2,n3);
885: outerr(M0_pro); flusherr();
886: }
887: if (gcmp0(M0) || gcmp(M0_pro,M0)<0)
888: {
889: M0=M0_pro; tetpil=avma; M0=gerepile(av1,tetpil,gcopy(M0));
890: }
891: else avma=av1;
892: }
893: else if (n1==n2 || n1==n3 || n2==n3)
894: {
895: if (n1==n2) k=n1; else if (n2==n3) k=n3;
896: kk=N-2*k;
897: p2=gsub(M_star,gmulgs(x,k));
898: p3=gmul(gpuigs(stoi(kk),kk),gpuigs(gsubgs(gmul(M_star,p2),kk*kk),k));
899: pol=gsub(p3,gmul(gmul(gpuigs(stoi(k),k),gpuigs(x,k)),gpuigs(p2,N-k)));
900: prec=gprecision(pol); if (!prec) prec=5;
901: r=roots(pol,prec); lr=lg(r)-1;
902: for (i=1; i<=lr; i++)
903: {
904: if (gcmp(gabs(gimag((GEN)r[i]),prec),eps) < 0 &&
905: gsigne(S=greal((GEN)r[i])) > 0)
906: {
907: p4=gsub(M_star,gmulsg(k,S));
908: P=gdiv(gmul(gmulsg(k,S),p4),gsubgs(gmul(M_star,p4),kk*kk));
909: p5=gsub(gsqr(S),gmul2n(P,2));
910: if (gsigne(p5)>=0)
911: {
912: p6=gsqrt(p5,prec);
913: u=gmul2n(gadd(S,p6),-1); v=gmul2n(gsub(S,p6),-1);
914: if ((gsigne(u)>0)&&(gsigne(v)>0))
915: {
916: w=gpui(P,gdivgs(stoi(-k),kk),prec);
917: p6=gmulsg(k,gadd(gsqr(glog(u,prec)),gsqr(glog(v,prec))));
918: M0_pro=gmul2n(gadd(p6,gmulsg(kk,gsqr(glog(w,prec)))),-2);
919: if (DEBUGLEVEL>2)
920: {
921: fprintferr("[%ld,%ld,%ld] : ",n1,n2,n3);
922: outerr(M0_pro); flusherr();
923: }
924: if (gcmp0(M0) || gcmp(M0_pro,M0)<0) M0=M0_pro;
925: }
926: }
927: }
928: }
929: tetpil=avma; M0=gerepile(av1,tetpil,gcopy(M0));
930: }
931: else
932: {
933: f1=gadd(gmulsg(n1,x),gadd(gmulsg(n2,y),gmulsg(n3,z)));
934: f1=gsub(f1,M_formel);
935: f2=gmulsg(n1,gmul(y,z));
936: f2=gadd(f2,gmulsg(n2,gmul(x,z)));
937: f2=gadd(f2,gmulsg(n3,gmul(x,y)));
938: f2=gsub(f2,gmul(M_formel,gmul(x,gmul(y,z))));
939: f3=gsub(gmul(gpuigs(x,n1),gmul(gpuigs(y,n2),gpuigs(z,n3))),gun);
940: g1=subres(f1,f2); g1=gdiv(g1,content(g1));
941: g2=subres(f1,f3); g2=gdiv(g2,content(g2));
942: g3=subres(g1,g2); g3=gdiv(g3,content(g3));
943: pf1=gsubst(f1,vM,M_star); pg1=gsubst(g1,vM,M_star);
944: pf2=gsubst(f2,vM,M_star); pg2=gsubst(g2,vM,M_star);
945: pf3=gsubst(f3,vM,M_star); pg3=gsubst(g3,vM,M_star);
946: prec=gprecision(pg3); if (!prec) prec=5;
947: r=roots(pg3,prec); lr=lg(r)-1;
948: for (i=1; i<=lr; i++)
949: {
950: if (gcmp(gabs(gimag((GEN)r[i]),prec),eps) < 0 &&
951: gsigne(w=greal((GEN)r[i])) > 0)
952: {
953: p1=gsubst(pg1,vz,w); p2=gsubst(pg2,vz,w);
954: p3=gsubst(pf1,vz,w); p4=gsubst(pf2,vz,w); p5=gsubst(pf3,vz,w);
955: prec=gprecision(p1); if (!prec) prec=5;
956: r1=roots(p1,prec); lr1=lg(r1)-1;
957: for (j=1; j<=lr1; j++)
958: {
959: if (gcmp(gabs(gimag((GEN)r1[j]),prec),eps) < 0 &&
960: gsigne(v=greal((GEN)r1[j])) > 0)
961: {
962: p6=gsubst(p2,vy,v);
963: if (gcmp(gabs(p6,prec),eps)<0)
964: {
965: p7=gsubst(p3,vy,v); p8=gsubst(p4,vy,v); p9=gsubst(p5,vy,v);
966: prec=gprecision(p7); if (!prec) prec=5;
967: r2=roots(p7,prec); lr2=lg(r2)-1;
968: for (l=1; l<=lr2; l++)
969: {
970: if (gcmp(gabs(gimag((GEN)r2[l]),prec),eps) < 0 &&
971: gsigne(u=greal((GEN)r2[l])) > 0)
972: {
973: p10=gsubst(p8,vx,u);
974: if (gcmp(gabs(p10,prec),eps)<0)
975: {
976: p11=gsubst(p9,vx,u);
977: if (gcmp(gabs(p11,prec),eps)<0)
978: {
979: M0_pro=gmulsg(n1,gsqr(glog(u,prec)));
980: M0_pro=gadd(M0_pro,gmulsg(n2,gsqr(glog(v,prec))));
981: M0_pro=gadd(M0_pro,gmulsg(n3,gsqr(glog(w,prec))));
982: M0_pro=gmul2n(M0_pro,-2);
983: if (DEBUGLEVEL>2)
984: {
985: fprintferr("[ %ld,%ld,%ld ] : ",n1,n2,n3);
986: outerr(M0_pro); flusherr();
987: }
988: if (gcmp0(M0) || gcmp(M0_pro,M0) < 0) M0=M0_pro;
989: }
990: }
991: }
992: }
993: }
994: }
995: }
996: }
997: }
998: tetpil=av1; M0_pro=gerepile(av1,tetpil,gcopy(M0));
999: }
1000: }
1001: }
1002: for (i=1;i<=4;i++) delete_var();
1003: return M0;
1004: }
1005:
1006: static GEN
1007: lowerboundforregulator(GEN bnf,GEN units)
1008: {
1009: long N,R1,R2,RU,i,nbrootsofone,nbmin;
1010: GEN rootsofone,nf,M0,M,m,col,T2,bound,minunit,newminunit;
1011: GEN vecminim,colalg,p1,pol,y;
1012:
1013: rootsofone=gmael(bnf,8,4); nbrootsofone=itos((GEN)rootsofone[1]);
1014: nf=(GEN)bnf[7]; T2=gmael(nf,5,3); N=lgef(nf[1])-3;
1015: R1=itos(gmael(nf,2,1)); R2=itos(gmael(nf,2,2)); RU=R1+R2-1;
1016: if (RU==0) return gun;
1017:
1018: units=algtobasis(bnf,units); minunit=qfeval(T2,(GEN)units[1]);
1019: for (i=2; i<=RU; i++)
1020: {
1021: newminunit=qfeval(T2,(GEN)units[i]);
1022: if (gcmp(newminunit,minunit)<0) minunit=newminunit;
1023: }
1024: if (gcmpgs(minunit,1000000000)>0) return regulatorbound(bnf);
1025:
1026: vecminim=minimforunits(nf,itos(gceil(minunit)),10000);
1027: if (!vecminim) return regulatorbound(bnf);
1028: m=(GEN)vecminim[3]; nbmin=lg(m)-1;
1029: if (nbmin==10000) return regulatorbound(bnf);
1030: bound=gaddgs(minunit,1);
1031: for (i=1; i<=nbmin; i++)
1032: {
1033: col=(GEN)m[i]; colalg=basistoalg(nf,col);
1034: if (!gcmp1(lift_intern(gpuigs(colalg,nbrootsofone))))
1035: {
1036: newminunit=qfeval(T2,col);
1037: if (gcmp(newminunit,bound)<0) bound=newminunit;
1038: }
1039: }
1040: if (gcmp(bound,minunit)>0) err(talker,"bug in lowerboundforregulator");
1041: if (DEBUGLEVEL>=2)
1042: {
1043: fprintferr("M* = "); outerr(bound); flusherr();
1044: if (DEBUGLEVEL>2)
1045: {
1046: p1=polx[0]; pol=gaddgs(gsub(gpuigs(p1,N),gmul(bound,p1)),N-1);
1047: fprintferr("pol = "); outerr(pol); flusherr();
1048: p1=roots(pol,DEFAULTPREC);
1049: if (N&1) y=greal((GEN)p1[3]); else y=greal((GEN)p1[2]);
1050: fprintferr("y = "); outerr(y); flusherr();
1051: M0=gmul2n(gmulsg(N*(N-1),gsqr(glog(y,DEFAULTPREC))),-2);
1052: fprintferr("old fashion M0 = "); outerr(M0); flusherr();
1053: }
1054: }
1055: M0=compute_M0(bound,N);
1056: if (DEBUGLEVEL>=2) { fprintferr("M0 = "); outerr(M0); flusherr(); }
1057: M=gmul2n(gdivgs(gdiv(gpuigs(M0,RU),hermiteconstant(RU)),N),R2);
1058: if (gcmp(M,dbltor(0.04))<0) return regulatorbound(bnf);
1059: M = gsqrt(M,DEFAULTPREC);
1060: if (DEBUGLEVEL>=2)
1061: { fprintferr("(lower bound for regulator) M = "); outerr(M); flusherr(); }
1062: return M;
1063: }
1064:
1065: /* Calcule une matrice carree de rang lg(beta) associee a une famille
1066: * d'ideaux premiers P_i tels que : 1<=i<=nombre de beta; N(P_i) congru a 1
1067: * mod pp v_(P_i)(beta[j])=0 pour tout 1<=j<=nbalphapro et 1<=i<=lg(beta)
1068: */
1069: static void
1070: primecertify(GEN bnf,GEN beta,long pp,GEN big)
1071: {
1072: long i,j,qq,nbcol,sizeofmat,nbqq,ra,N;
1073: GEN nf,mat,mat1,qgen,decqq,newcol,eltgen,qrhall,Q;
1074:
1075: nbcol=0; nf=(GEN)bnf[7]; N=lgef(nf[1])-3;
1076: sizeofmat=lg(beta)-1; mat=cgetg(1,t_MAT); qq=1;
1077: for(;;)
1078: {
1079: qq += 2*pp; qgen=stoi(qq);
1080: if (smodis(big,qq)==0 || !isprime(qgen)) continue;
1081:
1082: decqq=primedec(bnf,qgen); nbqq=lg(decqq)-1;
1083: for (i=1; i<=nbqq; i++)
1084: {
1085: Q=(GEN)decqq[i];
1086: if (!gcmp1((GEN)Q[4])) continue;
1087:
1088: qrhall=nfmodprinit(nf,Q); nbcol++;
1089: newcol=cgetg(sizeofmat+1,t_COL);
1090: if (DEBUGLEVEL>=2)
1091: fprintferr(" prime ideal Q: %Z\n",Q);
1092: eltgen = gscalcol_i(lift(gener(qgen)), N);
1093: for (j=1; j<=sizeofmat; j++)
1094: newcol[j]=(long)nfshanks(nf,(GEN)beta[j],eltgen,Q,qrhall);
1095: if (DEBUGLEVEL>=2)
1096: {
1097: fprintferr(" generator of (Zk/Q)^*: "); outerr(eltgen);
1098: fprintferr(" column #%ld of the matrix log(b_j/Q): ",nbcol);
1099: outerr(newcol);
1100: }
1101: mat1=concatsp(mat,newcol); ra=rank(mat1);
1102: if (DEBUGLEVEL>=2)
1103: {
1104: fprintferr(" new rank of the matrix: %ld\n\n",ra); flusherr();
1105: }
1106: if (ra!=nbcol) nbcol--;
1107: else
1108: {
1109: if (nbcol==sizeofmat) return;
1110: mat=mat1;
1111: }
1112: }
1113: }
1114: }
1115:
1116: static void
1117: check_prime(long p, GEN bnf, GEN h, GEN cyc, long R, GEN alpha,
1118: GEN funits, GEN rootsofone, GEN big)
1119: {
1120: long av = avma, nbpro,nbalpha,i, nbgen = lg(cyc)-1;
1121: GEN p1,beta, nf = (GEN)bnf[7], nbrootsofone = (GEN)rootsofone[1];
1122:
1123: if (DEBUGLEVEL>=2)
1124: { fprintferr("***** Testing prime p = %ld\n",p); flusherr(); }
1125: if (smodis(h,p)) nbpro=0;
1126: else
1127: {
1128: if (DEBUGLEVEL>=2) { fprintferr(" p divides cl(k)\n"); flusherr(); }
1129: for (nbpro=nbgen; nbpro; nbpro--)
1130: if (!smodis((GEN)cyc[nbpro],p)) break;
1131: }
1132: nbalpha=nbpro+R;
1133: if (smodis(nbrootsofone,p)) p1 = (GEN) alpha[nbpro];
1134: else
1135: {
1136: if (DEBUGLEVEL>=2) { fprintferr(" p divides w(k)\n"); flusherr(); }
1137: nbalpha++; nbpro++; p1 = algtobasis(nf,(GEN)rootsofone[2]);
1138: }
1139: if (DEBUGLEVEL>=2)
1140: { fprintferr(" t+r+e = %ld\n",nbalpha); flusherr(); }
1141: beta=cgetg(nbalpha+1,t_VEC);
1142: if (nbpro)
1143: {
1144: for (i=1; i<nbpro; i++) beta[i]=alpha[i];
1145: beta[nbpro]=(long) p1;
1146: }
1147: for (i=1; i<=R; i++)
1148: beta[i+nbpro]=(long)algtobasis(nf,(GEN)funits[i]);
1149: if (DEBUGLEVEL>=2)
1150: { fprintferr(" Beta list = "); outerr(beta); flusherr(); }
1151: primecertify(bnf,beta,p,big); avma=av;
1152: }
1153:
1154: long
1155: certifybuchall(GEN bnf)
1156: {
1157: long av = avma, nbgen,i,j,pp,N,R1,R2,R,nfa,nbf1,bound;
1158: GEN big,nf,reg,rootsofone,funits,gen,p1,gbound,DK,alpha,factfd1,f1,h,cyc;
1159: byteptr delta = diffptr;
1160:
1161: bnf = checkbnf(bnf); nf = (GEN)bnf[7];
1162: N=lgef(nf[1])-3; if (N==1) return 1;
1163: R1=itos(gmael(nf,2,1)); R2=itos(gmael(nf,2,2)); R=R1+R2-1;
1164: funits = check_units(bnf,"bnfcertify");
1165: DK=absi((GEN)nf[3]);
1166: testprime(bnf,zimmertbound(N,R1,DK));
1167: p1=gmael(bnf,8,1); reg=gmael(bnf,8,2);
1168: h=(GEN)p1[1];
1169: cyc=(GEN)p1[2]; nbgen=lg(cyc)-1;
1170: gen=(GEN)p1[3]; rootsofone=gmael(bnf,8,4);
1171: if (DEBUGLEVEL>1)
1172: {
1173: fprintferr("Class number = "); outerr(h);
1174: fprintferr("Cyclic components = "); outerr(cyc);
1175: fprintferr("Generators = "); outerr(gen);
1176: fprintferr("Regulator = "); outerr(reg);
1177: fprintferr("Roots of one = "); outerr(rootsofone);
1178: fprintferr("Fundamental units = "); outerr(funits);
1179: }
1180: alpha=cgetg(nbgen+1,t_VEC);
1181: for (i=1; i<=nbgen; i++)
1182: {
1183: p1=idealpow(nf,(GEN)gen[i],(GEN)cyc[i]);
1184: alpha[i]=isprincipalall(bnf,p1,nf_GEN | nf_FORCE)[2];
1185: }
1186: gbound = ground(gdiv(reg,lowerboundforregulator(bnf,funits)));
1187: if (is_bigint(gbound))
1188: err(talker,"sorry, too many primes to check");
1189:
1190: bound = gbound[2];
1191: if (DEBUGLEVEL>=2)
1192: {
1193: fprintferr("\nPHASE 2: are all primes good ?\n\n");
1194: fprintferr(" Testing primes <= B (= %ld)\n\n",bound); flusherr();
1195: }
1196: for (big=gun,i=1; i<=nbgen; i++)
1197: big = mulii(big,idealnorm(nf,(GEN)gen[i]));
1198: for (pp = *delta++; pp <= bound; pp += *delta++)
1199: check_prime(pp,bnf,h,cyc,R,alpha,funits,rootsofone,big);
1200:
1201: nfa = 0;
1202: if (nbgen)
1203: {
1204: factfd1=factor((GEN)cyc[1]);
1205: nbf1=lg(factfd1[1]); f1=(GEN)factfd1[1];
1206: for (i=1; i<nbf1; i++)
1207: if (cmpis((GEN)f1[i],bound) > 0) nfa++;
1208: }
1209: if (DEBUGLEVEL>=2 && nfa)
1210: { fprintferr(" Testing primes > B (# = %ld)\n\n",nfa); flusherr(); }
1211: for (j=1; j<=nfa; j++)
1212: {
1213: pp = itos((GEN)f1[nbf1-j]);
1214: check_prime(pp,bnf,gzero,cyc,R,alpha,funits,rootsofone,big);
1215: }
1216: avma=av; return 1;
1217: }
1218:
1219: /*******************************************************************/
1220: /* */
1221: /* CORPS DE CLASSES DE RAYON : CONDUCTEURS ET DISCRIMINANTS */
1222: /* */
1223: /*******************************************************************/
1224:
1225: /* Si s est la surjection de Cl_m sur Cl_n et H ssgroupe de Cl_m,
1226: * retourne le ssgroupe s(H) de Cl_n
1227: */
1228: static GEN
1229: imageofgroup0(GEN gen,GEN bnr,GEN subgroup)
1230: {
1231: long j,l;
1232: GEN E,Delta = diagonal(gmael(bnr,5,2));
1233:
1234: if (gcmp0(subgroup)) return Delta;
1235:
1236: l=lg(gen); E=cgetg(l,t_MAT);
1237: for (j=1; j<l; j++)
1238: E[j] = (long)isprincipalray(bnr,(GEN)gen[j]);
1239: E = concatsp(gmul(E,subgroup),Delta);
1240: return hnf(E);
1241: }
1242:
1243: static GEN
1244: imageofgroup(GEN gen,GEN bnr,GEN subgroup)
1245: {
1246: long av = avma;
1247: return gerepileupto(av,imageofgroup0(gen,bnr,subgroup));
1248: }
1249:
1250: /* retourne le cardinal de Cl_n / s(H) */
1251: static GEN
1252: cardofimagofquotientgroup(GEN H,GEN bnr2,GEN subgroup)
1253: {
1254: return dethnf_i(imageofgroup0(H,bnr2,subgroup));
1255: }
1256:
1257: static GEN
1258: args_to_bnr(GEN arg0, GEN arg1, GEN arg2, GEN *subgroup, long prec)
1259: {
1260: GEN bnr,bnf;
1261:
1262: if (typ(arg0)!=t_VEC)
1263: err(talker,"neither bnf nor bnr in conductor or discray");
1264: if (!arg1) arg1 = gzero;
1265: if (!arg2) arg2 = gzero;
1266:
1267: switch(lg(arg0))
1268: {
1269: case 7: /* bnr */
1270: bnr=arg0; (void)checkbnf((GEN)bnr[1]);
1271: *subgroup=arg1; break;
1272:
1273: case 11: /* bnf */
1274: bnf = checkbnf(arg0);
1275: bnr=buchrayall(bnf,arg1,nf_INIT | nf_GEN,prec);
1276: *subgroup=arg2; break;
1277:
1278: default: err(talker,"neither bnf nor bnr in conductor or discray");
1279: }
1280: if (!gcmp0(*subgroup))
1281: {
1282: long tx=typ(*subgroup);
1283: if (!is_matvec_t(tx))
1284: err(talker,"bad subgroup in conductor or discray");
1285: }
1286: return bnr;
1287: }
1288:
1289: GEN
1290: bnrconductor(GEN arg0,GEN arg1,GEN arg2,long all,long prec)
1291: {
1292: GEN sub=arg1, bnr=args_to_bnr(arg0,arg1,arg2,&sub,prec);
1293: return conductor(bnr,sub,all,prec);
1294: }
1295:
1296: long
1297: bnrisconductor(GEN arg0,GEN arg1,GEN arg2,long prec)
1298: {
1299: GEN sub=arg1, bnr=args_to_bnr(arg0,arg1,arg2,&sub,prec);
1300: return itos(conductor(bnr,sub,-1,prec));
1301: }
1302:
1303: /* Given a number field bnf=bnr[1], a ray class group bnr (from buchrayinit),
1304: * and a subgroup (HNF form) of the ray class group, compute the conductor
1305: * of the subgroup (copy of discrayrelall) if all=0. If all > 0, compute
1306: * furthermore the corresponding subgroup and output
1307: * [[ideal,arch],[hm,cyc,gen],subgroup] if all = 1, and output
1308: * [[ideal,arch],newbnr,subgroup] if all = 2. If all<0, answer only 1 is
1309: * module is the conductor, 0 otherwise.
1310: */
1311: GEN
1312: conductor(GEN bnr,GEN subgroup,long all,long prec)
1313: {
1314: long r1,j,av=avma,tetpil,k,ep,trivial=0;
1315: GEN bnf,nf,cl,cyc,gen,bid,ideal,arch,p1,p2,clhray,clhss;
1316: GEN fa,arch2,bnr2,fa2,ex;
1317:
1318: checkbnrgen(bnr); bnf=(GEN)bnr[1]; bid=(GEN)bnr[2];
1319: cl=(GEN)bnr[5]; cyc=(GEN)cl[2]; gen=(GEN)cl[3];
1320: nf=(GEN)bnf[7]; r1=itos(gmael(nf,2,1));
1321: p1=(GEN)bid[1]; ideal=(GEN)p1[1]; arch=(GEN)p1[2];
1322: if (gcmp0(subgroup)) { trivial=1; clhray=(GEN)cl[1]; }
1323: else
1324: {
1325: p1 = gauss(subgroup,diagonal(cyc));
1326: if (!gcmp1(denom(p1)))
1327: err(talker,"incorrect subgroup in conductor");
1328: if (gcmp1(det(p1))) trivial=1;
1329: clhray = absi(det(subgroup));
1330: }
1331: fa=(GEN)bid[3]; fa2=(GEN)fa[1]; ex=(GEN)fa[2];
1332: p2=cgetg(3,t_VEC); p2[2]=(long)arch;
1333: for (k=1; k<lg(fa2); k++)
1334: {
1335: GEN pr=(GEN)fa2[k], prinv=idealinv(nf,pr);
1336: ep = (all>=0)? itos((GEN)ex[k]): 1;
1337: for (j=1; j<=ep; j++)
1338: {
1339: p2[1]=(long)idealmul(nf,ideal,prinv);
1340: if (trivial) clhss=rayclassno(bnf,p2);
1341: else
1342: {
1343: bnr2=buchrayall(bnf,p2,nf_INIT,prec);
1344: clhss=cardofimagofquotientgroup(gen,bnr2,subgroup);
1345: }
1346: if (!egalii(clhss,clhray)) break;
1347: if (all<0) { avma=av; return gzero; }
1348: ideal = (GEN)p2[1];
1349: }
1350: }
1351: p2[1]=(long)ideal; arch2=dummycopy(arch); p2[2]=(long)arch2;
1352: for (k=1; k<=r1; k++)
1353: if (signe(arch[k]))
1354: {
1355: arch2[k]=zero;
1356: if (trivial) clhss=rayclassno(bnf,p2);
1357: else
1358: {
1359: bnr2=buchrayall(bnf,p2,nf_INIT,prec);
1360: clhss=cardofimagofquotientgroup(gen,bnr2,subgroup);
1361: }
1362: if (!egalii(clhss,clhray)) arch2[k]=un;
1363: else if (all<0) { avma=av; return gzero; }
1364: }
1365: if (all<0) {avma=av; return gun;}
1366: if (!all) { tetpil=avma; return gerepile(av,tetpil,gcopy(p2)); }
1367:
1368: bnr2=buchrayall(bnf,p2,nf_INIT | nf_GEN,prec);
1369: tetpil=avma; p1=cgetg(4,t_VEC);
1370: p1[3]=(long)imageofgroup(gen,bnr2,subgroup);
1371: if (all==1) bnr2=(GEN)bnr2[5];
1372: p1[2]=lcopy(bnr2);
1373: p1[1]=lcopy(p2); return gerepile(av,tetpil,p1);
1374: }
1375:
1376: /* etant donne un bnr et un polynome relatif, trouve le groupe des normes
1377: correspondant a l'extension relative en supposant qu'elle est abelienne
1378: et que le module donnant bnr est multiple du conducteur (tout ceci n'etant
1379: bien sur pas verifie). */
1380: GEN
1381: rnfnormgroup(GEN bnr, GEN polrel)
1382: {
1383: long av=avma,i,j,reldeg,sizemat,prime,nfac,k;
1384: GEN bnf,polreldisc,nf,raycl,group,detgroup,fa,pr,famo,ep,fac,col,p1;
1385: byteptr d = diffptr;
1386:
1387: checkbnr(bnr); bnf=(GEN)bnr[1]; raycl=(GEN)bnr[5];
1388: if (typ(polrel)!=t_POL) err(typeer,"rnfnormgroup");
1389: reldeg=lgef(polrel)-3; detgroup=(GEN)raycl[1];
1390: group = diagonal((GEN)raycl[2]);
1391: k = cmpis(detgroup,reldeg);
1392: if (k<0) err(talker,"not an Abelian extension in rnfnormgroup?");
1393: if (!k) return group;
1394:
1395: polreldisc=discsr(polrel); nf=(GEN)bnf[7];
1396: sizemat=lg(group)-1; prime = *d++;
1397: /* tant que nffactormod est bugge pour p=2 on commence a prime = 3 */
1398: for(;;)
1399: {
1400: prime += *d++; if (!*d) err(primer1);
1401: fa=primedec(nf,stoi(prime));
1402: for (i=1; i<lg(fa); i++)
1403: {
1404: pr = (GEN)fa[i];
1405: if (element_val(nf,polreldisc,pr)==0)
1406: {
1407: famo=nffactormod(nf,polrel,pr);
1408: ep=(GEN)famo[2]; fac=(GEN)famo[1];
1409: nfac=lg(ep)-1; k=lgef((GEN)fac[1])-3;
1410: for (j=1; j<=nfac; j++)
1411: {
1412: if (!gcmp1((GEN)ep[j])) err(bugparier,"rnfnormgroup");
1413: if (lgef(fac[j])-3 != k)
1414: err(talker,"non Galois extension in rnfnormgroup");
1415: }
1416: col=gmulsg(k,isprincipalrayall(bnr,pr,nf_REGULAR));
1417: p1=cgetg(sizemat+2,t_MAT);
1418: for (j=1; j<=sizemat; j++) p1[j]=group[j];
1419: p1[sizemat+1]=(long)col;
1420: group=hnf(p1); detgroup=dethnf(group);
1421: k = cmpis(detgroup,reldeg);
1422: if (k<0) err(talker,"not an Abelian extension in rnfnormgroup?");
1423: if (!k) { cgiv(detgroup); return gerepileupto(av,group); }
1424: }
1425: }
1426: }
1427: }
1428:
1429: /* Etant donne un bnf et un polynome relatif polrel definissant une extension
1430: abelienne (ce qui n'est pas verifie), calcule le conducteur et le groupe de
1431: congruence correspondant a l'extension definie par polrel sous la forme
1432: [[ideal,arch],[hm,cyc,gen],group] ou [ideal,arch] est le conducteur, et
1433: [hm,cyc,gen] est le groupe de classes de rayon correspondant. */
1434: GEN
1435: rnfconductor(GEN bnf, GEN polrel, long prec)
1436: {
1437: long av=avma,tetpil,R1,i,v;
1438: GEN nf,module,arch,bnr,group,p1,pol2;
1439:
1440: bnf = checkbnf(bnf); nf=(GEN)bnf[7];
1441: module=cgetg(3,t_VEC); R1=itos(gmael(nf,2,1));
1442: v=varn(polrel);
1443: p1=unifpol((GEN)bnf[7],polrel,0);
1444: p1=denom(gtovec(p1));
1445: pol2=gsubst(polrel,v,gdiv(polx[v],p1));
1446: pol2=gmul(pol2,gpuigs(p1,degree(pol2)));
1447: module[1]=rnfdiscf(nf,pol2)[1]; arch=cgetg(R1+1,t_VEC);
1448: module[2]=(long)arch; for (i=1; i<=R1; i++) arch[i]=un;
1449: bnr=buchrayall(bnf,module,nf_INIT | nf_GEN,prec);
1450: group=rnfnormgroup(bnr,pol2); tetpil=avma;
1451: return gerepile(av,tetpil,conductor(bnr,group,1,prec));
1452: }
1453:
1454: /* Etant donnes un corps de nombres bnf=bnr[1], un groupe de classes de rayon
1455: * bnr calcule par buchrayinit(), et une matrice HNF subgroup definissant un
1456: * sous-groupe du groupe de classes de rayon, calcule [n,r1,dk] associe a ce
1457: * sous groupe (cf. discrayall()). si flcond=0 le calcul est arrete si le
1458: * module n'est pas le conducteur et le programme retourne gzero. Si
1459: * flrel=0, seul la norme de l'ideal dk est calculee au lieu de dk lui-meme
1460: */
1461: static GEN
1462: discrayrelall(GEN bnr,GEN subgroup,long flag,long prec)
1463: {
1464: long r1,degk,j,av=avma,tetpil,k,ep,nbdezero;
1465: long trivial=0, flrel = flag & nf_REL, flcond = flag & nf_COND;
1466: GEN bnf,nf,cyc,gen,bid,ideal,arch,p1,p2,clhray,clhss;
1467: GEN fa,dlk,arch2,bnr2,idealrel,fa2,ex,som;
1468:
1469: checkbnrgen(bnr); bnf=(GEN)bnr[1];
1470: cyc=gmael(bnr,5,2); gen=gmael(bnr,5,3);
1471: nf=(GEN)bnf[7]; r1=itos(gmael(nf,2,1));
1472:
1473: if (gcmp0(subgroup)) { trivial=1; clhray=gmael(bnr,5,1); }
1474: else
1475: {
1476: p1=gauss(subgroup,diagonal(cyc));
1477: if (!gcmp1(denom(p1)))
1478: err(talker,"incorrect subgroup in discray");
1479: if (gcmp1(det(p1))) trivial=1;
1480: clhray = det(subgroup);
1481: }
1482: bid=(GEN)bnr[2]; ideal=gmael(bid,1,1); arch=gmael(bid,1,2);
1483: fa=(GEN)bid[3]; fa2=(GEN)fa[1]; ex=(GEN)fa[2];
1484:
1485: degk=lgef(nf[1])-3;
1486: idealrel=flrel?idmat(degk):gun;
1487:
1488: p2=cgetg(3,t_VEC); p2[2]=(long)arch;
1489: for (k=1; k<lg(fa2); k++)
1490: {
1491: GEN pr=(GEN)fa2[k], prinv=idealinv(nf,pr);
1492:
1493: ep=itos((GEN)ex[k]); p2[1]=(long)ideal; som=gzero;
1494: for (j=1; j<=ep; j++)
1495: {
1496: p2[1]=(long)idealmul(nf,(GEN)p2[1],prinv);
1497: if (trivial) clhss=rayclassno(bnf,p2);
1498: else
1499: {
1500: bnr2=buchrayall(bnf,p2,nf_INIT,prec);
1501: clhss=cardofimagofquotientgroup(gen,bnr2,subgroup);
1502: }
1503: if (j==1 && egalii(clhss,clhray) && flcond) { avma=av; return gzero; }
1504: if (is_pm1(clhss)) { som = addis(som, ep-j+1); break; }
1505: som = addii(som, clhss);
1506: }
1507: idealrel = flrel
1508: ? idealmul(nf,idealrel, idealpow(nf,pr, som))
1509: : mulii(idealrel, powgi((GEN)pr[1], mulii(som,(GEN)pr[4])));
1510: }
1511: dlk = flrel
1512: ? idealdiv(nf,idealpow(nf,ideal,clhray),idealrel)
1513: : divii(powgi(dethnf(ideal),clhray),idealrel);
1514:
1515: p2[1]=(long)ideal; arch2=dummycopy(arch); p2[2]=(long)arch2; nbdezero=0;
1516: for (k=1; k<=r1; k++)
1517: {
1518: if (signe(arch[k]))
1519: {
1520: arch2[k]=zero;
1521: if (trivial) clhss=rayclassno(bnf,p2);
1522: else
1523: {
1524: bnr2=buchrayall(bnf,p2,nf_INIT,prec);
1525: clhss=cardofimagofquotientgroup(gen,bnr2,subgroup);
1526: }
1527: arch2[k]=un;
1528: if (egalii(clhss,clhray))
1529: {
1530: if (flcond) { avma=av; return gzero; }
1531: nbdezero++;
1532: }
1533: }
1534: else nbdezero++;
1535: }
1536: tetpil=avma; p1=cgetg(4,t_VEC);
1537: p1[1]=lcopy(clhray);
1538: p1[2]=lstoi(nbdezero);
1539: p1[3]=lcopy(dlk); return gerepile(av,tetpil,p1);
1540: }
1541:
1542: static GEN
1543: discrayabsall(GEN bnr, GEN subgroup,long flag,long prec)
1544: {
1545: long av=avma,tetpil,degk,clhray,n,R1;
1546: GEN p1,p2,p3,p4,nf,dkabs,bnf;
1547:
1548: p1=discrayrelall(bnr,subgroup,flag,prec);
1549: if (flag & nf_REL) return p1;
1550: if (p1==gzero) { avma=av; return gzero; }
1551:
1552: bnf=(GEN)bnr[1]; nf=(GEN)bnf[7]; degk=lgef(nf[1])-3;
1553: dkabs=absi((GEN)nf[3]); p2=(GEN)p1[3];
1554: clhray=itos((GEN)p1[1]); p3=gpuigs(dkabs,clhray);
1555: n = degk * clhray; R1 = itos((GEN)p1[2]) * clhray;
1556: if (((n-R1)&3)==2) p2=negi(p2);
1557: tetpil=avma; p4=cgetg(4,t_VEC);
1558: p4[1]=lstoi(n);
1559: p4[2]=lstoi(R1);
1560: p4[3]=lmulii(p2,p3); return gerepile(av,tetpil,p4);
1561: }
1562:
1563: GEN
1564: bnrdisc0(GEN arg0, GEN arg1, GEN arg2, long flag, long prec)
1565: {
1566: GEN subgroup, bnr = args_to_bnr(arg0,arg1,arg2,&subgroup,prec);
1567: return discrayabsall(bnr,subgroup,flag,prec);
1568: }
1569:
1570: GEN
1571: discrayrel(GEN bnr, GEN subgroup,long prec)
1572: {
1573: return discrayrelall(bnr,subgroup,nf_REL,prec);
1574: }
1575:
1576: GEN
1577: discrayrelcond(GEN bnr, GEN subgroup,long prec)
1578: {
1579: return discrayrelall(bnr,subgroup,nf_REL | nf_COND,prec);
1580: }
1581:
1582: GEN
1583: discrayabs(GEN bnr, GEN subgroup,long prec)
1584: {
1585: return discrayabsall(bnr,subgroup,nf_REGULAR,prec);
1586: }
1587:
1588: GEN
1589: discrayabscond(GEN bnr, GEN subgroup,long prec)
1590: {
1591: return discrayabsall(bnr,subgroup,nf_COND,prec);
1592: }
1593:
1594: /* Etant donnes un corps de nombres bnf=bnr[1], un groupe de classes de rayon
1595: * bnr calcule par buchrayinit(), et un vecteur chi representant un caractere
1596: * sur les generateurs bnr[2][3], calcule le conducteur de ce caractere.
1597: */
1598: GEN
1599: bnrconductorofchar(GEN bnr,GEN chi,long prec)
1600: {
1601: long nbgen,i,av=avma,tetpil;
1602: GEN p1,m,u,d1,cl,cyclic;
1603:
1604: checkbnrgen(bnr); cl=(GEN)bnr[5];
1605: cyclic=(GEN)cl[2]; nbgen=lg(cyclic)-1;
1606: if (lg(chi)-1 != nbgen)
1607: err(talker,"incorrect character length in conductorofchar");
1608: if (!nbgen) return conductor(bnr,gzero,0,prec);
1609:
1610: d1=(GEN)cyclic[1]; m=cgetg(nbgen+2,t_MAT);
1611: for (i=1; i<=nbgen; i++)
1612: {
1613: p1=cgetg(2,t_COL); m[i]=(long)p1;
1614: p1[1]=ldiv(gmul((GEN)chi[i],d1),(GEN)cyclic[i]);
1615: if (typ(p1[1])!=t_INT) err(typeer,"conductorofchar");
1616: }
1617: p1=cgetg(2,t_COL); m[i]=(long)p1;
1618: p1[1]=(long)d1; u=(GEN)hnfall(m)[2];
1619: tetpil=avma; setlg(u,nbgen+1);
1620: for (i=1; i<=nbgen; i++) setlg(u[i],nbgen+1);
1621: return gerepile(av,tetpil,conductor(bnr,u,0,prec));
1622: }
1623:
1624: /* Etant donne la liste des zidealstarinit et des matrices d'unites
1625: * correspondantes, sort la liste des nombres de classes
1626: */
1627: GEN
1628: rayclassnolist(GEN bnf,GEN listes)
1629: {
1630: long av=avma,tetpil,i,j,k,kk,nc,nq,lx,ly;
1631: GEN h,modulist,unitlist,classlist,sous,sousu,sousclass,p2,m,bid,q,cyclic;
1632:
1633: if (typ(listes)!=t_VEC || lg(listes)!=3) err(typeer,"rayclassnolist");
1634: bnf = checkbnf(bnf); h=gmael3(bnf,8,1,1);
1635: modulist=(GEN)listes[1]; unitlist=(GEN)listes[2];
1636: lx=lg(modulist); classlist=cgetg(lx,t_VEC);
1637: for (i=1; i<lx; i++)
1638: {
1639: sous=(GEN)modulist[i]; sousu=(GEN)unitlist[i]; ly=lg(sous);
1640: sousclass=cgetg(ly,t_VEC); classlist[i]=(long)sousclass;
1641: for (j=1; j<ly; j++)
1642: {
1643: bid=(GEN)sous[j]; q=(GEN)sousu[j]; nq=lg(q)-1;
1644: cyclic=gmael(bid,2,2); nc=lg(cyclic)-1;
1645: if (lg(q[1]) != nc+1) err(bugparier,"rayclassnolist");
1646: m=cgetg(nq+nc+1,t_MAT);
1647: for (k=1; k<=nq; k++) m[k]=q[k];
1648: for ( ; k<=nq+nc; k++)
1649: {
1650: p2=cgetg(nc+1,t_COL); m[k]=(long)p2;
1651: for (kk=1; kk<=nc; kk++) p2[kk]=(kk==k-nq)?cyclic[kk]:zero;
1652: }
1653: sousclass[j] = lmul(h,dethnf(hnf(m)));
1654: }
1655: }
1656: tetpil=avma; return gerepile(av,tetpil,gcopy(classlist));
1657: }
1658:
1659: static GEN
1660: rayclassnolistes(GEN sous, GEN sousclass, GEN fac)
1661: {
1662: long i;
1663: for (i=1; i<lg(sous); i++)
1664: if (gegal(gmael(sous,i,3),fac)) return (GEN)sousclass[i];
1665: err(bugparier,"discrayabslist");
1666: return NULL; /* not reached */
1667: }
1668:
1669: static GEN
1670: rayclassnolistessimp(GEN sous, GEN fac)
1671: {
1672: long i;
1673: for (i=1; i<lg(sous); i++)
1674: if (gegal(gmael(sous,i,1),fac)) return gmael(sous,i,2);
1675: err(bugparier,"discrayabslistlong");
1676: return NULL; /* not reached */
1677: }
1678:
1679: static GEN
1680: factormul(GEN fa1,GEN fa2)
1681: {
1682: GEN pr,prnew,ex,exnew,v,p,y=cgetg(3,t_MAT);
1683: long i,c,lx;
1684:
1685: pr=concatsp((GEN)fa1[1],(GEN)fa2[1]); y[1]=(long)pr;
1686: ex=concatsp((GEN)fa1[2],(GEN)fa2[2]); y[2]=(long)ex;
1687: v=sindexsort(pr); lx=lg(pr);
1688: prnew=cgetg(lx,t_COL); exnew=cgetg(lx,t_COL);
1689: for (i=1; i<lx; i++) prnew[i]=pr[v[i]];
1690: for (i=1; i<lx; i++) exnew[i]=ex[v[i]];
1691: p=gzero; c=0;
1692: for (i=1; i<lx; i++)
1693: {
1694: if (gegal((GEN)prnew[i],p))
1695: ex[c]=laddii((GEN)ex[c],(GEN)exnew[i]);
1696: else
1697: {
1698: c++; p=(GEN)prnew[i]; pr[c]=(long)p; ex[c]=exnew[i];
1699: }
1700: }
1701: setlg(pr,c+1); setlg(ex,c+1); return y;
1702: }
1703:
1704: static GEN
1705: factordivexact(GEN fa1,GEN fa2)
1706: {
1707: long i,j,k,c,lx1,lx2;
1708: GEN listpr,listex,y,listpr1,listex1,listpr2,listex2,p1;
1709:
1710: listpr1=(GEN)fa1[1]; listex1=(GEN)fa1[2]; lx1=lg(listpr1);
1711: listpr2=(GEN)fa2[1]; listex2=(GEN)fa2[2]; lx2=lg(listpr1);
1712: y=cgetg(3,t_MAT);
1713: listpr=cgetg(lx1,t_COL); y[1]=(long)listpr;
1714: listex=cgetg(lx1,t_COL); y[2]=(long)listex;
1715: for (c=0,i=1; i<lx1; i++)
1716: {
1717: j=isinvector(listpr2,(GEN)listpr1[i],lx2-1);
1718: if (!j) { c++; listpr[c]=listpr1[i]; listex[c]=listex1[i]; }
1719: else
1720: {
1721: p1=subii((GEN)listex1[i],(GEN)listex2[j]); k=signe(p1);
1722: if (k<0) err(talker,"factordivexact is not exact!");
1723: if (k>0) { c++; listpr[c]=listpr1[i]; listex[c]=(long)p1; }
1724: }
1725: }
1726: setlg(listpr,c+1); setlg(listex,c+1); return y;
1727: }
1728:
1729: static GEN
1730: factorpow(GEN fa,long n)
1731: {
1732: GEN y=cgetg(3,t_MAT);
1733:
1734: if (!n) { y[1]=lgetg(1,t_COL); y[2]=lgetg(1,t_COL); return y; }
1735: y[1]=fa[1]; y[2]=lmulsg(n,(GEN)fa[2]); return y;
1736: }
1737:
1738: /* Etant donne la liste des zidealstarinit et des matrices d'unites
1739: * correspondantes, sort la liste des discrayabs. On ne garde pas les modules
1740: * qui ne sont pas des conducteurs
1741: */
1742: GEN
1743: discrayabslist(GEN bnf,GEN listes)
1744: {
1745: long av=avma,tetpil,ii,jj,i,j,k,clhss,ep,clhray,lx,ly,r1,som,degk,nbdezero;
1746: long n,R1,normps,normi,lfa2;
1747: GEN classlist,modulist,disclist,nf,dkabs,sous,sousclass,sousdisc;
1748: GEN bid,module,ideal,arch,fa,fa2,ex,idealrel,p1,p2,pr;
1749: GEN dlk,arch2,p3,fac,normp,fad,fad1,fad2,no1,no2;
1750:
1751: classlist=rayclassnolist(bnf,listes); lx=lg(classlist);
1752: modulist=(GEN)listes[1];
1753: disclist=cgetg(lx,t_VEC); nf=(GEN)bnf[7]; r1=itos(gmael(nf,2,1));
1754: degk=lgef(nf[1])-3; dkabs=gabs((GEN)nf[3],0);
1755: for (ii=1; ii<lx; ii++)
1756: {
1757: sous=(GEN)modulist[ii]; sousclass=(GEN)classlist[ii];
1758: ly=lg(sous); sousdisc=cgetg(ly,t_VEC); disclist[ii]=(long)sousdisc;
1759: for (jj=1; jj<ly; jj++)
1760: {
1761: bid=(GEN)sous[jj]; clhray=itos((GEN)sousclass[jj]);
1762: module=(GEN)bid[1]; ideal=(GEN)module[1]; arch=(GEN)module[2];
1763: fa=(GEN)bid[3]; fa2=(GEN)fa[1]; ex=(GEN)fa[2]; fac=gcopy(fa);
1764: idealrel=cgetg(3,t_MAT);
1765: idealrel[1]=lgetg(1,t_COL);
1766: idealrel[2]=lgetg(1,t_COL); lfa2=lg(fa2);
1767: for (k=1; k<lfa2; k++)
1768: {
1769: pr=(GEN)fa2[k]; ep=itos((GEN)ex[k]); som=0;
1770: normp=gpui((GEN)pr[1],(GEN)pr[4],0);
1771: normps=itos(normp); normi=ii;
1772: normp=cgetg(3,t_MAT);
1773: no1=cgetg(2,t_COL); no1[1]=pr[1]; normp[1]=(long)no1;
1774: no2=cgetg(2,t_COL); no2[1]=pr[4]; normp[2]=(long)no2;
1775: for (j=1; j<=ep; j++)
1776: {
1777: if (j<ep) { coeff(fac,k,2)=lstoi(ep-j); fad=fac; }
1778: else
1779: {
1780: fad=cgetg(3,t_MAT);
1781: fad1=cgetg(lfa2-1,t_COL); fad[1]=(long)fad1;
1782: fad2=cgetg(lfa2-1,t_COL); fad[2]=(long)fad2;
1783: for (i=1; i<k; i++)
1784: { fad1[i]=coeff(fac,i,1); fad2[i]=coeff(fac,i,2); }
1785: for ( ; i<lfa2-1; i++)
1786: { fad1[i]=coeff(fac,i+1,1); fad2[i]=coeff(fac,i+1,2); }
1787: }
1788: normi /= normps;
1789: clhss = itos(rayclassnolistes((GEN)modulist[normi],
1790: (GEN)classlist[normi],fad));
1791: if (j==1 && clhss==clhray)
1792: {
1793: clhray=0; coeff(fac,k,2)=ex[k]; goto LLDISCRAY;
1794: }
1795: if (clhss==1) { som += ep-j+1; break; }
1796: som += clhss;
1797: }
1798: coeff(fac,k,2)=ex[k];
1799: idealrel=factormul(idealrel,factorpow(normp,som));
1800: }
1801: dlk=factordivexact(factorpow(factor(stoi(ii)),clhray),idealrel);
1802: p2=cgetg(3,t_VEC);
1803: p2[1]=(long)ideal; arch2=dummycopy(arch);
1804: p2[2]=(long)arch2; nbdezero=0;
1805: for (k=1; k<=r1; k++)
1806: {
1807: if (signe(arch[k]))
1808: {
1809: arch2[k]=zero;
1810: clhss=itos(rayclassno(bnf,p2));
1811: arch2[k]=un;
1812: if (clhss==clhray) { clhray=0; break; }
1813: }
1814: else nbdezero++;
1815: }
1816: LLDISCRAY:
1817: if (!clhray) sousdisc[jj]=lgetg(1,t_VEC);
1818: else
1819: {
1820: p3=factorpow(factor(dkabs),clhray); n=clhray*degk; R1=nbdezero*clhray;
1821: if (((n-R1)&3) == 2)
1822: {
1823: normp=cgetg(3,t_MAT);
1824: no1=cgetg(2,t_COL); normp[1]=(long)no1; no1[1]=lstoi(-1);
1825: no2=cgetg(2,t_COL); normp[2]=(long)no2; no2[1]=un;
1826: dlk=factormul(normp,dlk);
1827: }
1828: p1=cgetg(4,t_VEC); p1[1]=lstoi(n);
1829: p1[2]=lstoi(R1); p1[3]=(long)factormul(dlk,p3);
1830: sousdisc[jj]=(long)p1;
1831: }
1832: }
1833: }
1834: tetpil=avma; return gerepile(av,tetpil,gcopy(disclist));
1835: }
1836:
1837: #define SHLGVINT 15
1838: #define LGVINT 32768 /* must be 1<<SHLGVINT */
1839:
1840: /* Attention: bound est le nombre de vraies composantes voulues. */
1841: static GEN
1842: bigcgetvec(long bound)
1843: {
1844: long nbcext,nbfinal,i;
1845: GEN vext;
1846:
1847: nbcext = ((bound-1)>>SHLGVINT)+1;
1848: nbfinal = bound-((nbcext-1)<<SHLGVINT);
1849: vext = cgetg(nbcext+1,t_VEC);
1850: for (i=1; i<nbcext; i++) vext[i]=lgetg(LGVINT+1,t_VEC);
1851: vext[nbcext]=lgetg(nbfinal+1,t_VEC); return vext;
1852: }
1853:
1854: static GEN
1855: getcompobig(GEN vext,long i)
1856: {
1857: long cext;
1858:
1859: if (i<=LGVINT) return gmael(vext,1,i);
1860: cext=((i-1)>>SHLGVINT)+1;
1861: return gmael(vext,cext,i-((cext-1)<<SHLGVINT));
1862: }
1863:
1864: static void
1865: putcompobig(GEN vext,long i,GEN x)
1866: {
1867: long cext;
1868:
1869: if (i<=LGVINT) { mael(vext,1,i)=(long)x; return; }
1870: cext=((i-1)>>SHLGVINT)+1; mael(vext,cext,i-((cext-1)<<SHLGVINT))=(long)x;
1871: return;
1872: }
1873:
1874: static GEN
1875: zsimp(GEN bid, GEN matunit)
1876: {
1877: GEN y=cgetg(5,t_VEC);
1878: y[1]=lcopy((GEN)bid[3]); y[2]=lcopy(gmael(bid,2,2));
1879: y[3]=lcopy((GEN)bid[5]); y[4]=lcopy(matunit); return y;
1880: }
1881:
1882: static GEN
1883: zsimpjoin(GEN bidsimp, GEN bidp, GEN faussefa, GEN matunit)
1884: {
1885: long i,j,lx,lx1,lx2,llx,llx1,llx2,nbgen,av=avma,tetpil,c;
1886: GEN y,U1,U2,cyclic1,cyclic2,U,cyc,u1u2,p1,p2,met;
1887:
1888: y=cgetg(5,t_VEC); y[1]=(long)vconcat((GEN)bidsimp[1],faussefa);
1889: U1=(GEN)bidsimp[3]; U2=(GEN)bidp[5]; cyclic1=(GEN)bidsimp[2];
1890: cyclic2=gmael(bidp,2,2); lx1=lg(U1); lx2=lg(U2); lx=lx1+lx2-1;
1891: llx1=lg(cyclic1); llx2=lg(cyclic2);
1892: llx=llx1+llx2-1; nbgen=llx-1; U=cgetg(lx,t_MAT);
1893: if (nbgen)
1894: {
1895: cyc=diagonal(concatsp(cyclic1,cyclic2));
1896: u1u2=smithclean(smith2(cyc)); met=(GEN)u1u2[3]; c=lg(met)-1;
1897: for (j=1; j<lx1; j++)
1898: {
1899: p1=cgetg(llx,t_COL); p2=(GEN)U1[j]; U[j]=(long)p1;
1900: for (i=1; i<llx1; i++) p1[i]=p2[i];
1901: for ( ; i<llx; i++) p1[i]=zero;
1902: }
1903: for ( ; j<lx; j++)
1904: {
1905: p1=cgetg(llx,t_COL); p2=(GEN)U2[j-lx1+1]; U[j]=(long)p1;
1906: for (i=1; i<llx1; i++) p1[i]=zero;
1907: for ( ; i<llx; i++) p1[i]=p2[i-llx1+1];
1908: }
1909: y[3]=lmul((GEN)u1u2[1],U);
1910: }
1911: else
1912: {
1913: met=cgetg(1,t_MAT); for (j=1; j<lx; j++) U[j]=lgetg(1,t_COL);
1914: y[3]=(long)U; c=0;
1915: }
1916: cyc=cgetg(c+1,t_VEC); for (i=1; i<=c; i++) cyc[i]=coeff(met,i,i);
1917: y[2]=(long)cyc;
1918: y[4]=(long)vconcat((GEN)bidsimp[4],matunit);
1919: tetpil=avma; return gerepile(av,tetpil,gcopy(y));
1920: }
1921:
1922: static GEN
1923: rayclassnointern(GEN sous, GEN clh)
1924: {
1925: long lx,nc,nq,k,kk,j;
1926: GEN bidsimp,qm,sousray,cyclic,m,p2,p1;
1927:
1928: lx=lg(sous); sousray=cgetg(lx,t_VEC);
1929: for (j=1; j<lx; j++)
1930: {
1931: bidsimp=(GEN)sous[j]; qm=gmul((GEN)bidsimp[3],(GEN)bidsimp[4]);
1932: nq=lg(qm)-1; cyclic=(GEN)bidsimp[2]; nc=lg(cyclic)-1;
1933: if (lg(qm[1]) != nc+1) err(bugparier,"rayclassnolist");
1934: m=cgetg(nq+nc+1,t_MAT);
1935: for (k=1; k<=nq; k++) m[k]=qm[k];
1936: for ( ; k<=nq+nc; k++)
1937: {
1938: p2=cgetg(nc+1,t_COL); m[k]=(long)p2;
1939: for (kk=1; kk<=nc; kk++) p2[kk]=(kk==k-nq)?cyclic[kk]:zero;
1940: }
1941: p1=cgetg(3,t_VEC); p1[2]=lmul(clh,dethnf(hnf(m)));
1942: p1[1]=bidsimp[1]; sousray[j]=(long)p1;
1943: }
1944: return sousray;
1945: }
1946:
1947: static GEN
1948: rayclassnointernarch(GEN sous, GEN clh, GEN matarchunit)
1949: {
1950: long lx,nc,nq,k,kk,j,lm,lh,r1,jj,i,nba,nbarch,ii;
1951: GEN bidsimp,qm,sousray,cyclic,m,p2,p1,p1all,p3,mm,mj,qmk,matk;
1952:
1953: if (!matarchunit) return rayclassnointern(sous,clh);
1954:
1955: lm=lg(matarchunit); if (!lm) err(talker,"no units in rayclassnointernarch");
1956: r1=lg(matarchunit[1])-1; if (r1>15) err(talker,"r1>15 in rayclassnointernarch");
1957: lx=lg(sous); sousray=cgetg(lx,t_VEC);
1958: for (j=1; j<lx; j++)
1959: {
1960: bidsimp=(GEN)sous[j]; qm=gmul((GEN)bidsimp[3],(GEN)bidsimp[4]);
1961: nq=lg(qm)-1; cyclic=(GEN)bidsimp[2]; nc=lg(cyclic)-1;
1962: if (lm != nq+1) err(bugparier,"rayclassnointernarch (1)");
1963: if (lg(qm[1]) != nc+1) err(bugparier,"rayclassnointernarch (2)");
1964: m=cgetg(nq+nc+r1+1,t_MAT);
1965: for (k=1; k<=nq; k++)
1966: {
1967: p2=cgetg(nc+r1+1,t_COL); m[k]=(long)p2; qmk=(GEN)qm[k];
1968: matk=(GEN)matarchunit[k];
1969: for (kk=1; kk<=nc; kk++) p2[kk]=qmk[kk];
1970: for ( ; kk<=nc+r1; kk++) p2[kk]=matk[kk-nc];
1971: }
1972: for ( ; k<=nq+nc; k++)
1973: {
1974: p2=cgetg(nc+r1+1,t_COL); m[k]=(long)p2;
1975: for (kk=1; kk<=nc+r1; kk++) p2[kk]=(kk==k-nq)?cyclic[kk]:zero;
1976: }
1977: for ( ; k<=nq+nc+r1; k++)
1978: {
1979: p2=cgetg(nc+r1+1,t_COL); m[k]=(long)p2;
1980: for (kk=1; kk<=nc+r1; kk++) p2[kk]=(kk==k-nq)?deux:zero;
1981: }
1982: m=hnf(m);
1983: nbarch=(1<<r1); p1all=cgetg(nbarch+1,t_VEC); lh=lg(m);
1984: if (lh!=nc+r1+1) err(bugparier,"rayclassnointernarch (3)");
1985: for (k=0; k<=nbarch-1; k++)
1986: {
1987: p2=cgetg(r1+1,t_COL); kk=k; nba=0;
1988: for (jj=1; jj<=r1; jj++)
1989: {
1990: if (kk%2) { nba++; p2[jj]=un; } else p2[jj]=zero;
1991: kk>>=1;
1992: }
1993: mm=cgetg(lh,t_MAT);
1994: for (jj=1; jj<lh; jj++)
1995: {
1996: p3=cgetg(nc+nba+1,t_COL); mm[jj]=(long)p3; mj=(GEN)m[jj];
1997: for (i=1; i<=nc; i++) p3[i]=mj[i];
1998: for (ii=1; ii<=r1; ii++)
1999: if (signe(p2[ii])) p3[i++]=mj[nc+ii];
2000: }
2001: p1all[k+1]=lmul(clh,dethnf(hnf(mm)));
2002: }
2003: p1=cgetg(3,t_VEC); p1[2]=(long)p1all; p1[1]=bidsimp[1];
2004: sousray[j]=(long)p1;
2005: }
2006: return sousray;
2007: }
2008:
2009: GEN
2010: decodemodule(GEN nf, GEN fa)
2011: {
2012: long n,k,j,fauxpr,av=avma;
2013: GEN fa2,ex,id,pr;
2014:
2015: nf=checknf(nf);
2016: if (typ(fa)!=t_MAT || lg(fa)!=3)
2017: err(talker,"incorrect factorisation in decodemodule");
2018: n=lgef(nf[1])-3; id=idmat(n);
2019: fa2=(GEN)fa[1]; ex=(GEN)fa[2];
2020: for (k=1; k<lg(fa2); k++)
2021: {
2022: fauxpr=itos((GEN)fa2[k]);
2023: j=(fauxpr%n)+1; fauxpr /= n*n;
2024: pr = (GEN)primedec(nf,stoi(fauxpr))[j];
2025: id = idealmul(nf,id,idealpow(nf,pr,(GEN)ex[k]));
2026: }
2027: return gerepileupto(av,id);
2028: }
2029:
2030: /* Fait tout a partir du debut, et bound peut aller jusqu'a 2^30. Pour
2031: * l'instant ne s'occupe pas des sous-groupes. le format de sortie est le
2032: * suivant: un vecteur vext dont les composantes vint ont exactement 2^LGVINT
2033: * vraies composantes sauf le dernier qui peut etre plus court. vext[i][j]
2034: * (ou j<=2^LGVINT) doit etre compris comme la composante d'indice
2035: * I=(i-1)*2^LGVINT+j d'un grand vecteur unique V. Une telle composante
2036: * d'indice I est un vecteur indexe par tous les ideaux de norme egal a I. Si
2037: * m_0 est un tel ideal, la composante correspondante est la suivante:
2038: *
2039: * + si arch = NULL, on parcourt toutes les parties archimediennes possibles.
2040: * L'ordre des arch est l'ordre lexicographique inverse, donc [0,0,..,0],
2041: * [1,0,..,0], [0,1,..,0],... Dans ce cas: [m_0,V] ou V est un vecteur de
2042: * 2^r1 composantes, donnant pour chaque arch, le triplet [N,R1,D], ou N, R1,
2043: * D sont comme dans discrayabs sauf que D est sous forme factorisee.
2044: *
2045: * + sinon [m_0,arch,N,R1,D], ou N, R1, D sont comme ci-dessus.
2046: *
2047: * Si ramip est different de 0 et -1, ne prend que les modules sans facteur
2048: * carres ailleurs qu'au dessus de ramip. Si ramip est strictement negatif,
2049: * archsquare.
2050: */
2051: static GEN
2052: discrayabslistarchintern(GEN bnf, GEN arch, long bound, long ramip)
2053: {
2054: byteptr ptdif=diffptr;
2055: long degk,lim,av0,av,av1,tetpil,i,j,k,p2s,lfa,lp1,sqbou,cex, allarch;
2056: long ffs,fs,resp,flbou,tdi,nba, k2,karch,kka,nbarch,jjj,jj,square;
2057: long ii2,ii,ly,clhray,lfa2,ep,som,clhss,normps,normi,nbdezero,r1,R1,n,lp4,c;
2058: ulong q;
2059: GEN nf,p,z,pol,p1,p2,p3,fa,pr,normp,ideal,bidp,z2,matarchunit;
2060: GEN bigres,funits,racunit,embunit,sous,clh,sousray,raylist;
2061: GEN clhrayall,discall,faall,Id,idealrel,idealrelinit;
2062: GEN sousdisc,module,fa2,ex,fac,no1,no2,fad,fad1,fad2,fadkabs,pz;
2063: GEN arch2,dlk,disclist,bidsimp,p4,faussefa,pex,fauxpr,gprime;
2064: GEN *gptr[2];
2065:
2066: /* ce qui suit recopie d'assez pres ideallistzstarall */
2067: if (DEBUGLEVEL>2) timer2();
2068: if (bound <= 0) err(talker,"non-positive bound in discrayabslist");
2069: av0=avma; bnf = checkbnf(bnf); flbou=0;
2070: nf=(GEN)bnf[7]; bigres=(GEN)bnf[8]; pol=(GEN)nf[1]; degk=lgef(pol)-3;
2071: r1=itos(gmael(nf,2,1)); fadkabs=factor(absi((GEN)nf[3]));
2072: clh=gmael(bigres,1,1);
2073: funits = check_units(bnf,"discrayabslistarchintern");
2074: racunit=gmael(bigres,4,2);
2075:
2076: if (ramip >= 0) square = 0;
2077: else
2078: {
2079: square = 1; ramip = -ramip;
2080: if (ramip==1) ramip=0;
2081: }
2082: nba = 0; allarch = (arch==NULL);
2083: if (allarch)
2084: { arch=cgetg(r1+1,t_VEC); for (i=1; i<=r1; i++) arch[i]=un; nba=r1; }
2085: else if (gcmp0(arch))
2086: { arch=cgetg(r1+1,t_VEC); for (i=1; i<=r1; i++) arch[i]=zero; }
2087: else
2088: {
2089: if (lg(arch)!=r1+1)
2090: err(talker,"incorrect archimedean argument in discrayabslistlong");
2091: for (i=1; i<=r1; i++) if (signe(arch[i])) nba++;
2092: if (nba) module=cgetg(3,t_VEC);
2093: }
2094: p1=cgetg(3,t_VEC); p1[1]=(long)idmat(degk); p1[2]=(long)arch;
2095: bidp=zidealstarinitall(nf,p1,0);
2096: matarchunit = allarch? logunitmatrix(nf,funits,racunit,bidp): (GEN)NULL;
2097:
2098: p=cgeti(3); p[1]=evalsigne(1)|evallgef(3);
2099: sqbou=itos(racine(stoi(bound)))+1;
2100: av=avma; lim=stack_lim(av,1);
2101: z=bigcgetvec(bound); for (i=2; i<=bound; i++) putcompobig(z,i,cgetg(1,t_VEC));
2102: if (allarch) bidp=zidealstarinitall(nf,idmat(degk),0);
2103: embunit=logunitmatrix(nf,funits,racunit,bidp);
2104: p1=cgetg(2,t_VEC); putcompobig(z,1,p1); p1[1]=(long)zsimp(bidp,embunit);
2105: if (DEBUGLEVEL>=2) fprintferr("Starting zidealstarunits computations\n");
2106: for (p[2]=0; p[2]<=bound; )
2107: {
2108: if (!*ptdif) err(primer1);
2109: p[2] += *ptdif++;
2110: if (!flbou && p[2]>sqbou)
2111: {
2112: flbou=1;
2113: if (DEBUGLEVEL>=2)
2114: { fprintferr("\nStarting rayclassno computations\n"); flusherr(); }
2115: tetpil=avma; z=gerepile(av,tetpil,gcopy(z));
2116: av1=avma; raylist=bigcgetvec(bound);
2117: /* maintenant on suit rayclassnolist */
2118: for (i=1; i<=bound; i++)
2119: {
2120: sous=getcompobig(z,i);
2121: sousray=rayclassnointernarch(sous,clh,matarchunit);
2122: putcompobig(raylist,i,sousray);
2123: }
2124: tetpil=avma; raylist=gerepile(av1,tetpil,gcopy(raylist));
2125: z2=bigcgetvec(sqbou);
2126: for (i=1; i<=sqbou; i++)
2127: { p1=gcopy(getcompobig(z,i)); putcompobig(z2,i,p1); }
2128: z = z2;
2129: }
2130: fa=primedec(nf,p); lfa=lg(fa)-1;
2131: for (j=1; j<=lfa; j++)
2132: {
2133: pr=(GEN)fa[j]; normp=powgi(p,(GEN)pr[4]); cex=0;
2134: if (DEBUGLEVEL>=2) { fprintferr("%ld ",p[2]); flusherr(); }
2135: if (gcmpgs(normp,bound)<=0)
2136: {
2137: fauxpr=stoi(p[2]*degk*degk+(itos((GEN)pr[4])-1)*degk+j-1);
2138: q=p2s=itos(normp); ideal=pr;
2139: while (q <= (ulong)bound)
2140: {
2141: bidp=zidealstarinitall(nf,ideal,0);
2142: faussefa=cgetg(3,t_MAT); p1=cgetg(2,t_COL);
2143: faussefa[1]=(long)p1; p1[1]=(long)fauxpr;
2144: pex=cgetg(2,t_COL); faussefa[2]=(long)pex;
2145: cex++; pex[1]=lstoi(cex);
2146: embunit=logunitmatrix(nf,funits,racunit,bidp);
2147: for (i=q; i<=bound; i+=q)
2148: {
2149: p1=getcompobig(z,i/q); lp1=lg(p1);
2150: if (lp1>1)
2151: {
2152: p2=cgetg(lp1,t_VEC); c=0;
2153: for (k=1; k<lp1; k++)
2154: {
2155: p3=(GEN)p1[k]; if (i>q) { p4=gmael(p3,1,1); lp4=lg(p4)-1; }
2156: if (i==q || !isinvector(p4,fauxpr,lp4))
2157: {
2158: c++;
2159: p2[c]=(long)zsimpjoin(p3,bidp,faussefa,embunit);
2160: }
2161: }
2162: setlg(p2,c+1);
2163: if (c)
2164: {
2165: if (p[2]<=sqbou)
2166: {
2167: pz=getcompobig(z,i);
2168: if (lg(pz)>1) putcompobig(z,i,concatsp(pz,p2));
2169: else putcompobig(z,i,p2);
2170: }
2171: else
2172: {
2173: sousray=rayclassnointernarch(p2,clh,matarchunit);
2174: pz=getcompobig(raylist,i);
2175: if (lg(pz)>1) putcompobig(raylist,i,concatsp(pz,sousray));
2176: else putcompobig(raylist,i,sousray);
2177: }
2178: }
2179: }
2180: }
2181: if (ramip && ramip % p[2]) q=bound+1;
2182: else
2183: {
2184: pz=mulss(q,p2s);
2185: q=(gcmpgs(pz,bound)>0)?bound+1:pz[2];
2186: if (q <= (ulong)bound) ideal=idealmul(nf,ideal,pr);
2187: }
2188: }
2189: }
2190: }
2191: if (low_stack(lim, stack_lim(av,1)))
2192: {
2193: if(DEBUGMEM>1) err(warnmem,"[1]: discrayabslistarch");
2194: if (!flbou)
2195: {
2196: if (DEBUGLEVEL>2)
2197: fprintferr("avma = %ld, t(z) = %ld ",avma-bot,taille2(z));
2198: tetpil=avma; z=gerepile(av,tetpil,gcopy(z));
2199: }
2200: else
2201: {
2202: if (DEBUGLEVEL>2)
2203: fprintferr("avma = %ld, t(r) = %ld ",avma-bot,taille2(raylist));
2204: gptr[0]=&z; gptr[1]=&raylist; gerepilemany(av,gptr,2);
2205: }
2206: if (DEBUGLEVEL>2) { fprintferr("avma = %ld ",avma-bot); flusherr(); }
2207: }
2208: }
2209: if (!flbou) /* maintenant on suit rayclassnolist */
2210: {
2211: if (DEBUGLEVEL>=2)
2212: { fprintferr("\nStarting rayclassno computations\n"); flusherr(); }
2213: tetpil=avma; z=gerepile(av,tetpil,gcopy(z));
2214: av1=avma; raylist=bigcgetvec(bound);
2215: for (i=1; i<=bound; i++)
2216: {
2217: sous=getcompobig(z,i);
2218: sousray=rayclassnointernarch(sous,clh,matarchunit);
2219: putcompobig(raylist,i,sousray);
2220: }
2221: }
2222: if (DEBUGLEVEL>2)
2223: fprintferr("avma = %ld, t(r) = %ld ",avma-bot,taille2(raylist));
2224: tetpil=avma; raylist=gerepile(av,tetpil,gcopy(raylist));
2225: if (DEBUGLEVEL>2)
2226: { fprintferr("avma = %ld ",avma-bot); msgtimer("zidealstarlist"); }
2227: /* maintenant on suit discrayabslist */
2228: if (DEBUGLEVEL>=2)
2229: { fprintferr("Starting discrayabs computations\n"); flusherr(); }
2230:
2231: if (allarch) nbarch = 1<<r1;
2232: else
2233: {
2234: nbarch = 1;
2235: clhrayall = cgetg(2,t_VEC);
2236: discall = cgetg(2,t_VEC);
2237: faall = cgetg(2,t_VEC);
2238: Id = idmat(degk);
2239: }
2240: idealrelinit=cgetg(3,t_MAT);
2241: idealrelinit[1]=lgetg(1,t_COL);
2242: idealrelinit[2]=lgetg(1,t_COL);
2243: av1=avma; lim=stack_lim(av1,1);
2244: if (square) bound = sqbou-1;
2245: disclist=bigcgetvec(bound);
2246: for (ii=1; ii<=bound; ii++) putcompobig(disclist,ii,cgetg(1,t_VEC));
2247: for (ii2=1; ii2<=bound; ii2++)
2248: {
2249: ii = square? ii2*ii2: ii2;
2250: if (DEBUGLEVEL>=2 && ii%100==0) { fprintferr("%ld ",ii); flusherr(); }
2251: sous=getcompobig(raylist,ii); ly=lg(sous); sousdisc=cgetg(ly,t_VEC);
2252: putcompobig(disclist, square? ii2: ii,sousdisc);
2253: for (jj=1; jj<ly; jj++)
2254: {
2255: bidsimp=(GEN)sous[jj]; fa=(GEN)bidsimp[1]; fac=gcopy(fa);
2256: fa2=(GEN)fa[1]; ex=(GEN)fa[2]; lfa2=lg(fa2);
2257:
2258: if (allarch)
2259: {
2260: clhrayall = (GEN)bidsimp[2];
2261: discall=cgetg(nbarch+1,t_VEC);
2262: }
2263: else
2264: clhray = itos((GEN)bidsimp[2]);
2265: for (karch=0; karch<nbarch; karch++)
2266: {
2267: if (!allarch) ideal = Id;
2268: else
2269: {
2270: kka=karch; nba=0;
2271: for (jjj=1; jjj<=r1; jjj++)
2272: {
2273: if (kka%2) nba++;
2274: kka>>=1;
2275: }
2276: clhray=itos((GEN)clhrayall[karch+1]);
2277: for (k2=1,k=1; k<=r1; k++,k2<<=1)
2278: {
2279: if (karch&k2 && clhray==itos((GEN)clhrayall[karch-k2+1]))
2280: { clhray=0; goto LDISCRAY; }
2281: }
2282: }
2283: idealrel = idealrelinit;
2284: for (k=1; k<lfa2; k++)
2285: {
2286: fauxpr=(GEN)fa2[k]; ep=itos((GEN)ex[k]); ffs=itos(fauxpr);
2287: /* Hack for NeXTgcc 2.5.8 -- splitting "resp=fs%degk+1" */
2288: fs=ffs/degk; resp=fs%degk; resp++; gprime=stoi((long)(fs/degk));
2289: if (!allarch && nba)
2290: {
2291: p1=(GEN)primedec(nf,gprime)[ffs%degk+1];
2292: ideal = idealmul(nf,ideal,idealpow(nf,p1,(GEN)ex[k]));
2293: }
2294: som=0; clhss=0;
2295: normp=gpuigs(gprime,resp); normps=itos(normp); normi=ii;
2296: normp=cgetg(3,t_MAT);
2297: no1=cgetg(2,t_COL); no1[1]=(long)gprime; normp[1]=(long)no1;
2298: no2=cgetg(2,t_COL); no2[1]=lstoi(resp); normp[2]=(long)no2;
2299: for (j=1; j<=ep; j++)
2300: {
2301: if (clhss==1) som++;
2302: else
2303: {
2304: if (j<ep) { coeff(fac,k,2)=lstoi(ep-j); fad=fac; }
2305: else
2306: {
2307: fad=cgetg(3,t_MAT);
2308: fad1=cgetg(lfa2-1,t_COL); fad[1]=(long)fad1;
2309: fad2=cgetg(lfa2-1,t_COL); fad[2]=(long)fad2;
2310: for (i=1; i<k; i++)
2311: { fad1[i]=coeff(fac,i,1); fad2[i]=coeff(fac,i,2); }
2312: for ( ; i<lfa2-1; i++)
2313: { fad1[i]=coeff(fac,i+1,1); fad2[i]=coeff(fac,i+1,2); }
2314: }
2315: normi /= normps;
2316: /* Hack for NeXTgcc 2.5.8 -- using dlk as temporary variable */
2317: dlk=rayclassnolistessimp(getcompobig(raylist,normi),fad);
2318: if (allarch) dlk = (GEN)dlk[karch+1];
2319: clhss = itos(dlk);
2320: if (j==1 && clhss==clhray)
2321: {
2322: clhray=0; coeff(fac,k,2)=ex[k]; goto LDISCRAY;
2323: }
2324: som += clhss;
2325: }
2326: }
2327: coeff(fac,k,2)=ex[k];
2328: idealrel=factormul(idealrel,factorpow(normp,som));
2329: }
2330: dlk=factordivexact(factorpow(factor(stoi(ii)),clhray),idealrel);
2331: if (!allarch && nba)
2332: {
2333: module[1]=(long)ideal; arch2=gcopy(arch); module[2]=(long)arch2;
2334: nbdezero=0;
2335: for (k=1; k<=r1; k++)
2336: {
2337: if (signe(arch[k]))
2338: {
2339: arch2[k]=zero;
2340: clhss=itos(rayclassno(bnf,module));
2341: arch2[k]=un;
2342: if (clhss==clhray) { clhray=0; goto LDISCRAY; }
2343: }
2344: else nbdezero++;
2345: }
2346: }
2347: else nbdezero = r1-nba;
2348: LDISCRAY:
2349: p1=cgetg(4,t_VEC); discall[karch+1]=(long)p1;
2350: if (!clhray) p1[1]=p1[2]=p1[3]=zero;
2351: else
2352: {
2353: p3=factorpow(fadkabs,clhray); n=clhray*degk; R1=nbdezero*clhray;
2354: if (((n-R1)&3)==2)
2355: {
2356: normp=cgetg(3,t_MAT);
2357: no1=cgetg(2,t_COL); normp[1]=(long)no1; no1[1]=lstoi(-1);
2358: no2=cgetg(2,t_COL); normp[2]=(long)no2; no2[1]=un;
2359: dlk=factormul(normp,dlk);
2360: }
2361: p1[1]=lstoi(n); p1[2]=lstoi(R1); p1[3]=(long)factormul(dlk,p3);
2362: }
2363: }
2364: if (allarch)
2365: { p1 = cgetg(3,t_VEC); p1[1]=(long)fa; p1[2]=(long)discall; }
2366: else
2367: { faall[1]=(long)fa; p1 = concatsp(faall, p1); }
2368: sousdisc[jj]=(long)p1;
2369: if (low_stack(lim, stack_lim(av1,1)))
2370: {
2371: if(DEBUGMEM>1) err(warnmem,"[2]: discrayabslistarch");
2372: if (DEBUGLEVEL>2)
2373: fprintferr("avma = %ld, t(d) = %ld ",avma-bot,taille2(disclist));
2374: tetpil=avma; disclist=gerepile(av1,tetpil,gcopy(disclist));
2375: if (DEBUGLEVEL>2) { fprintferr("avma = %ld ",avma-bot); flusherr(); }
2376: }
2377: }
2378: }
2379: if (DEBUGLEVEL>2) msgtimer("discrayabs");
2380: tdi=taille2(disclist);
2381: if (DEBUGLEVEL>2)
2382: { fprintferr("avma = %ld, t(d) = %ld ",avma-bot,tdi); flusherr(); }
2383: if (tdi<avma-bot)
2384: {
2385: tetpil=avma; disclist=gerepile(av0,tetpil,gcopy(disclist));
2386: if (DEBUGLEVEL>2) { fprintferr("avma = %ld ",avma-bot); flusherr(); }
2387: }
2388: return disclist;
2389: }
2390:
2391: GEN
2392: discrayabslistarch(GEN bnf, GEN arch, long bound)
2393: { return discrayabslistarchintern(bnf,arch,bound, 0); }
2394:
2395: GEN
2396: discrayabslistlong(GEN bnf, long bound)
2397: { return discrayabslistarchintern(bnf,gzero,bound, 0); }
2398:
2399: GEN
2400: discrayabslistarchsquare(GEN bnf, GEN arch, long bound)
2401: { return discrayabslistarchintern(bnf,arch,bound, -1); }
2402:
2403: static GEN
2404: computehuv(GEN bnr,GEN id, GEN arch,long prec)
2405: {
2406: long i,nbgen, av = avma;
2407: GEN bnf,bnrnew,listgen,P,U,DC;
2408: GEN newmod=cgetg(3,t_VEC);
2409: newmod[1]=(long)id;
2410: newmod[2]=(long)arch;
2411:
2412: bnf=(GEN)bnr[1];
2413: bnrnew=buchrayall(bnf,newmod,nf_INIT,prec);
2414: listgen=gmael(bnr,5,3); nbgen=lg(listgen);
2415: P=cgetg(nbgen,t_MAT);
2416: for (i=1; i<nbgen; i++)
2417: P[i] = (long)isprincipalray(bnrnew,(GEN)listgen[i]);
2418: DC=diagonal(gmael(bnrnew,5,2));
2419: U=(GEN)hnfall(concatsp(P,DC))[2];
2420: setlg(U,nbgen); for (i=1; i<nbgen; i++) setlg(U[i], nbgen);
2421: return gerepileupto(av, hnf(concatsp(U,diagonal(gmael(bnr,5,2)))));
2422: }
2423:
2424: /* 0 if hinv*list[i] has a denominator for all i, 1 otherwise */
2425: static int
2426: hnflistdivise(GEN list,GEN h)
2427: {
2428: long av = avma, i, I = lg(list);
2429: GEN hinv = ginv(h);
2430:
2431: for (i=1; i<I; i++)
2432: if (gcmp1(denom(gmul(hinv,(GEN)list[i])))) break;
2433: avma = av; return i < I;
2434: }
2435:
2436: static GEN
2437: subgroupcond(GEN bnr, long indexbound, long prec)
2438: {
2439: long av=avma,tetpil,i,j,lgi,lp;
2440: GEN li,p1,lidet,perm,nf,bid,ideal,arch,primelist,listkernels;
2441:
2442: checkbnrgen(bnr); bid=(GEN)bnr[2];
2443: ideal=gmael(bid,1,1);
2444: arch =gmael(bid,1,2); nf=gmael(bnr,1,7);
2445: primelist=gmael(bid,3,1); lp=lg(primelist)-1;
2446: listkernels=cgetg(lp+lg(arch),t_VEC);
2447: for (i=1; i<=lp; )
2448: {
2449: p1=idealdiv(nf,ideal,(GEN)primelist[i]);
2450: listkernels[i++]=(long)computehuv(bnr,p1,arch,prec);
2451: }
2452: for (j=1; j<lg(arch); j++)
2453: if (signe((GEN)arch[j]))
2454: {
2455: p1=dummycopy(arch); p1[j]=zero;
2456: listkernels[i++]=(long)computehuv(bnr,ideal,p1,prec);
2457: }
2458: setlg(listkernels,i);
2459:
2460: li=subgrouplist(gmael(bnr,5,2),indexbound);
2461: lgi=lg(li);
2462: for (i=1,j=1; j<lgi; j++)
2463: if (!hnflistdivise(listkernels, (GEN)li[j])) li[i++] = li[j];
2464: /* sort by increasing index */
2465: lgi = i; setlg(li,i); lidet=cgetg(lgi,t_VEC);
2466: for (i=1; i<lgi; i++) lidet[i]=(long)dethnf_i((GEN)li[i]);
2467: perm = sindexsort(lidet); p1=li; li=cgetg(lgi,t_VEC);
2468: for (i=1; i<lgi; i++) li[i] = p1[perm[lgi-i]];
2469: tetpil=avma; return gerepile(av,tetpil,gcopy(li));
2470: }
2471:
2472: GEN
2473: subgrouplist0(GEN bnr, long indexbound, long all, long prec)
2474: {
2475: if (typ(bnr)!=t_VEC) err(typeer,"subgrouplist");
2476: if (lg(bnr)!=1 && typ(bnr[1])!=t_INT)
2477: {
2478: if (!all) return subgroupcond(bnr,indexbound,prec);
2479: checkbnr(bnr); bnr = gmael(bnr,5,2);
2480: }
2481: return subgrouplist(bnr,indexbound);
2482: }
2483:
2484: GEN
2485: bnrdisclist0(GEN bnf, GEN borne, GEN arch, long all)
2486: {
2487: if (typ(borne)==t_INT)
2488: {
2489: if (!arch) arch = gzero;
2490: if (all==1) { arch = NULL; all = 0; }
2491: return discrayabslistarchintern(bnf,arch,itos(borne),all);
2492: }
2493: return discrayabslist(bnf,borne);
2494: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to buch3.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