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