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