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