Annotation of OpenXM_contrib/pari-2.2/src/basemath/buch3.c, Revision 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