/*******************************************************************/
/* */
/* RAY CLASS FIELDS */
/* */
/*******************************************************************/
/* $Id: buch3.c,v 1.2 1999/09/23 17:50:56 karim Exp $ */
#include "pari.h"
#include "parinf.h"
GEN compute_class_number(GEN mit,GEN *met,GEN *u1,GEN *u2);
GEN logunitmatrix(GEN nf,GEN funits,GEN racunit,GEN bid);
GEN vconcat(GEN Q1, GEN Q2);
GEN ideleaddone_aux(GEN nf,GEN x,GEN ideal);
static GEN
get_full_rank(GEN nf, GEN v, GEN _0, GEN _1, GEN vecsign, GEN gen,
long ngen, long rankmax)
{
GEN v1,p1,alpha, bas=(GEN)nf[7], rac=(GEN)nf[6];
long rankinit=lg(v)-1, N=lgef(nf[1])-3, va=varn(nf[1]);
long limr,i,k,kk,r,rr;
for (r=1,rr=3; ; r++,rr+=2)
{
p1 = gpuigs(stoi(rr),N);
limr=(cmpis(p1,BIGINT)>1)? BIGINT: p1[2]; /* min(BIGINT,rr^N) */
limr = (limr-1)>>1;
for (k=rr; k<=limr; k++)
{
long av1=avma;
alpha = gzero;
for (kk=k,i=1; i<=N; i++,kk/=rr)
{
long lambda = (kk+r)%rr - r;
if (lambda)
alpha = gadd(alpha,gmulsg(lambda,(GEN)bas[i]));
}
for (i=1; i<=rankmax; i++)
vecsign[i] = (gsigne(gsubst(alpha,va,(GEN)rac[i])) > 0)? (long)_0
: (long)_1;
v1 = concatsp(v, vecsign);
if (rank(v1) == rankinit) avma=av1;
else
{
v=v1; rankinit++; ngen++;
gen[ngen] = (long) alpha;
if (rankinit == rankmax) return ginv(v); /* full rank */
}
}
}
}
GEN
buchnarrow(GEN bnf)
{
GEN nf,_0mod2,_1mod2,cyc,gen,v,matsign,arch;
GEN dataunit,p1,p2,p3,h,vecsign,clh,basecl,met,u1,u2;
long R1,R,i,j,ngen,sizeh,t,lo,c;
long av=avma,tetpil;
if (typ(bnf)!=t_VEC || lg(bnf)!=11)
err(talker,"not a big number field vector in buchnarrow");
nf=checknf(bnf); R1=itos(gmael(nf,2,1));
if (!R1) return gcopy(gmael(bnf,8,1));
_1mod2=gmodulss(1,2); _0mod2=gmodulss(0,2);
v=cgetg(R1+1,t_COL); vecsign=cgetg(R1+1,t_COL);
for (i=1; i<=R1; i++) v[i]=(long)_1mod2;
cyc=gmael3(bnf,8,1,2); gen=gmael3(bnf,8,1,3); ngen=lg(gen)-1;
matsign=signunits(bnf); R=lg(matsign); dataunit=cgetg(R+1,t_MAT);
for (j=1; j<R; j++)
{
p1=cgetg(R1+1,t_COL); dataunit[j]=(long)p1;
for (i=1; i<=R1; i++)
p1[i] = (signe(gcoeff(matsign,i,j)) > 0)? (long)_0mod2: (long)_1mod2;
}
dataunit[R]=(long)v; v=image(dataunit); t=lg(v)-1;
sizeh=ngen+R1-t; p1 = cgetg(sizeh+1,t_COL);
for (i=1; i<=ngen; i++) p1[i]=gen[i];
gen = p1;
if (t<R1)
v = get_full_rank(nf,v,_0mod2,_1mod2,vecsign,gen,ngen,R1);
h=cgetg(sizeh+1,t_MAT); arch = cgetg(R1+1,t_VEC);
for (i=1; i<=R1; i++) arch[i]=un;
for (j=1; j<=ngen; j++)
{
p1 = cgetg(sizeh+1,t_COL); h[j]=(long)p1;
p2 = idealpow(nf, (GEN)gen[j], (GEN)cyc[j]);
p2 = (GEN)isprincipalall(bnf,p2,nf_GEN | nf_FORCE)[2];
p2 = gmul(v,zsigne(nf,p2,arch));
for (i=1; i<=ngen; i++) p1[i] = (i==j)? cyc[j]: zero;
for ( ; i<=sizeh; i++) p1[i] = llift((GEN)p2[i+t-ngen]);
}
for ( ; j<=sizeh; j++)
{
p1 = cgetg(sizeh+1,t_COL); h[j]=(long)p1;
for (i=1; i<=sizeh; i++) p1[i]=(i==j)?deux:zero;
}
clh=compute_class_number(h,&met,&u1,&u2);
u1=reducemodmatrix(u1,h); lo=lg(met)-1; c=0;
for (c=1; c<=lo; c++)
if (gcmp1(gcoeff(met,c,c))) break;
basecl=cgetg(c,t_VEC);
for (j=1; j<c; j++)
{
p1=gcoeff(u1,1,j);
p3=idealpow(nf,(GEN)gen[1],p1);
if (signe(p1)<0) p3=numer(p3);
for (i=2; i<=lo; i++)
{
p1=gcoeff(u1,i,j);
if (signe(p1))
{
p3 = idealmul(nf,p3, idealpow(nf,(GEN)gen[i],p1));
p1 = content(p3); if (!gcmp1(p1)) p3 = gdiv(p3,p1);
}
}
basecl[j]=(long)p3;
}
tetpil=avma; v=cgetg(4,t_VEC);
v[1]=lcopy(clh); setlg(met,c);
v[2]=(long)mattodiagonal(met);
v[3]=lcopy(basecl); return gerepile(av,tetpil,v);
}
GEN idealaddtoone_i(GEN nf, GEN x, GEN y);
/* given two coprime ideals x (integral) and id, compute alpha in x,
* alpha = 1 mod (id), with x/alpha nearly reduced.
*/
static GEN
findalpha(GEN nf,GEN x,GEN id,long prec)
{
GEN p1,idprod,y;
GEN alp = idealaddtoone_i(nf,x,id);
idprod = idealmullll(nf,x,id);
y = lllgram(qf_base_change(gmael(nf,5,3),idprod,1), 2*prec-2);
y = gmul(idprod, (GEN)y[1]); /* small vector in idprod */
p1 = ground(element_div(nf,alp,y));
alp = gsub(alp,element_mul(nf,p1,y));
return gcmp0(alp)? y: alp;
}
static int
too_big(GEN nf, GEN bet)
{
GEN x = gnorm(basistoalg(nf,bet));
switch (typ(x))
{
case t_INT: return absi_cmp(x, gun);
case t_FRAC: return absi_cmp((GEN)x[1], (GEN)x[2]);
}
err(bugparier, "wrong type in too_big");
return 0; /* not reached */
}
static GEN
idealmodidele(GEN nf, GEN x, GEN ideal, GEN sarch, GEN arch, long prec)
{
long av = avma,i,l;
GEN p1,p2,alp,bet,b;
nf=checknf(nf); alp=findalpha(nf,x,ideal,prec);
p1=idealdiv(nf,alp,x);
bet = element_div(nf,findalpha(nf,p1,ideal,prec),alp);
if (too_big(nf,bet) > 0) { avma=av; return x; }
p1=(GEN)sarch[2]; l=lg(p1);
if (l > 1)
{
b=bet; p2=lift_intern(gmul((GEN)sarch[3],zsigne(nf,bet,arch)));
for (i=1; i<l; i++)
if (signe(p2[i])) bet = element_mul(nf,bet,(GEN)p1[i]);
if (b != bet && too_big(nf,bet) > 0) { avma=av; return x; }
}
return idealmul(nf,bet,x);
}
static GEN
idealmulmodidele(GEN nf,GEN x,GEN y, GEN ideal,GEN sarch,GEN arch,long prec)
{
return idealmodidele(nf,idealmul(nf,x,y),ideal,sarch,arch,prec);
}
/* assume n > 0 */
static GEN
idealpowmodidele(GEN nf,GEN x,GEN n, GEN ideal,GEN sarch,GEN arch,long prec)
{
long i,m,av=avma;
GEN y;
ulong j;
if (cmpis(n, 16) < 0)
{
if (gcmp1(n)) return x;
x = idealpow(nf,x,n);
x = idealmodidele(nf,x,ideal,sarch,arch,prec);
return gerepileupto(av,x);
}
i = lgefint(n)-1; m=n[i]; j=HIGHBIT;
while ((m&j)==0) j>>=1;
y = x;
for (j>>=1; j; j>>=1)
{
y = idealmul(nf,y,y);
if (m&j) y = idealmul(nf,y,x);
y = idealmodidele(nf,y,ideal,sarch,arch,prec);
}
for (i--; i>=2; i--)
for (m=n[i],j=HIGHBIT; j; j>>=1)
{
y = idealmul(nf,y,y);
if (m&j) y = idealmul(nf,y,x);
y = idealmodidele(nf,y,ideal,sarch,arch,prec);
}
return gerepileupto(av,y);
}
static GEN
buchrayall(GEN bnf,GEN module,long flag,long prec)
{
GEN nf,cyc,gen,genplus,fa2,sarch,hmatu,u,clg;
GEN dataunit,p1,p2,h,clh,basecl,met,u1,u2,u1old;
GEN racunit,bigres,bid,resbid2,resbid3,x,y,funits,hmat,vecel;
long RU,R3,i,j,ngen,lh,lo,c,av=avma,N;
bnf = checkbnf(bnf); nf=checknf(bnf); bigres=(GEN)bnf[8];
funits = check_units(bnf, "buchrayall");
N=lgef(nf[1])-3;
cyc=gmael(bigres,1,2);
gen=gmael(bigres,1,3); ngen=lg(cyc)-1;
bid = zidealstarinitall(nf,module,1);
resbid2=gmael(bid,2,2);
resbid3=gmael(bid,2,3);
R3=lg(resbid2)-1; lh=ngen+R3;
x = idealhermite(nf,module);
if (R3 || flag & (nf_INIT|nf_GEN))
{
vecel=cgetg(ngen+1,t_VEC);
for (j=1; j<=ngen; j++)
vecel[j]=(long)idealcoprime(nf,(GEN)gen[j],x);
}
if (flag & nf_GEN)
{
genplus=cgetg(lh+1,t_VEC);
for (j=1; j<=ngen; j++)
genplus[j] = (long) idealmul(nf,(GEN)vecel[j],(GEN)gen[j]);
for ( ; j<=lh; j++)
genplus[j] = resbid3[j-ngen];
}
if (!R3)
{
if (!(flag & nf_GEN)) clg=cgetg(3,t_VEC);
else
{ clg=cgetg(4,t_VEC); clg[3]=(long)genplus; }
clg[1]=mael(bigres,1,1);
clg[2]=(long)cyc;
if (!(flag & nf_INIT)) return gerepileupto(av,gcopy(clg));
y = cgetg(7,t_VEC);
y[1]=lcopy(bnf);
y[2]=lcopy(bid);
y[3]=lcopy(vecel);
y[4]=(long)idmat(ngen);
y[5]=lcopy(clg); u=cgetg(3,t_VEC);
y[6]=(long)u;
u[1]=lgetg(1,t_MAT);
u[2]=(long)idmat(lg(funits));
return gerepileupto(av,y);
}
fa2=(GEN)bid[4]; sarch=(GEN)fa2[lg(fa2)-1];
RU=lg(funits); dataunit=cgetg(RU+R3+1,t_MAT);
racunit=gmael(bigres,4,2);
dataunit[1] = (long)zideallog(nf,racunit,bid);
for (j=2; j<=RU; j++)
dataunit[j] = (long)zideallog(nf,(GEN)funits[j-1],bid);
for ( ; j<=RU+R3; j++)
{
p1=cgetg(R3+1,t_COL); dataunit[j]=(long)p1;
for (i=1; i<=R3; i++)
p1[i] = (i==(j-RU))? resbid2[i]: zero;
}
h=cgetg(lh+1,t_MAT);
for (j=1; j<=ngen; j++)
{
p1=cgetg(lh+1,t_COL); h[j]=(long)p1;
p2 = idealpow(nf, (GEN)gen[j], (GEN)cyc[j]);
p2 = (GEN)isprincipalall(bnf,p2,nf_GEN | nf_FORCE)[2];
p2 = element_mul(nf,p2,element_pow(nf,(GEN)vecel[j],(GEN)cyc[j]));
p2 = zideallog(nf,p2,bid);
for (i=1; i<=ngen; i++) p1[i] = (i==j)? cyc[j]: zero;
for ( ; i<=lh; i++) p1[i] = lnegi((GEN)p2[i-ngen]);
}
hmatu=hnfall(dataunit); hmat=(GEN)hmatu[1];
for ( ; j<=lh; j++)
{
p1=cgetg(lh+1,t_COL); h[j]=(long)p1;
for (i=1; i<=ngen; i++) p1[i]=zero;
for ( ; i<=lh; i++) p1[i]=coeff(hmat,i-ngen,j-ngen);
}
clh = compute_class_number(h,&met,&u1,&u2);
u1old=u1; lo=lg(met)-1;
for (c=1; c<=lo; c++)
if (gcmp1(gcoeff(met,c,c))) break;
if (flag & nf_GEN)
{
GEN Id=idmat(N), arch=(GEN)module[2];
u1 = reducemodmatrix(u1,h);
basecl=cgetg(c,t_VEC);
for (j=1; j<c; j++)
{
GEN *op, minus = Id, plus = Id;
long av1 = avma, s;
for (i=1; i<=lo; i++)
{
p1 = gcoeff(u1,i,j);
if (!(s = signe(p1))) continue;
if (s > 0) op = + else { op = − p1 = negi(p1); }
p1 = idealpowmodidele(nf,(GEN)genplus[i],p1,x,sarch,arch,prec);
*op = *op==Id? p1
: idealmulmodidele(nf,*op,p1,x,sarch,arch,prec);
}
if (minus == Id) p1 = plus;
else
{
p2 = ideleaddone_aux(nf,minus,module);
p1 = idealdivexact(nf,idealmul(nf,p2,plus),minus);
p1 = idealmodidele(nf,p1,x,sarch,arch,prec);
}
basecl[j]=lpileupto(av1,p1);
}
clg=cgetg(4,t_VEC); clg[3]=lcopy(basecl);
} else clg=cgetg(3,t_VEC);
clg[1]=licopy(clh); setlg(met,c);
clg[2]=(long)mattodiagonal(met);
if (!(flag & nf_INIT)) return gerepileupto(av,clg);
u2 = cgetg(R3+1,t_MAT);
u1 = cgetg(RU+1,t_MAT); u = (GEN)hmatu[2];
for (j=1; j<=RU; j++) { u1[j]=u[j]; setlg(u[j],RU+1); }
u += RU;
for (j=1; j<=R3; j++) { u2[j]=u[j]; setlg(u[j],RU+1); }
p1=lllint(u1); p2=ginv(hmat);
y=cgetg(7,t_VEC);
y[1]=lcopy(bnf);
y[2]=lcopy(bid);
y[3]=lcopy(vecel);
y[4]=linv(u1old);
y[5]=lcopy(clg); u=cgetg(3,t_VEC);
y[6]=(long)u;
u[1]=lmul(u2,p2);
u[2]=lmul(u1,p1);
return gerepileupto(av,y);
}
GEN
buchrayinitgen(GEN bignf, GEN ideal,long prec)
{
return buchrayall(bignf,ideal, nf_INIT | nf_GEN,prec);
}
GEN
buchrayinit(GEN bignf, GEN ideal,long prec)
{
return buchrayall(bignf,ideal, nf_INIT,prec);
}
GEN
buchray(GEN bignf, GEN ideal,long prec)
{
return buchrayall(bignf,ideal, nf_GEN,prec);
}
GEN
bnrclass0(GEN bignf, GEN ideal, long flag, long prec)
{
switch(flag)
{
case 0: flag = nf_GEN; break;
case 1: flag = nf_INIT; break;
case 2: flag = nf_INIT | nf_GEN; break;
default: err(flagerr,"bnrclass");
}
return buchrayall(bignf,ideal,flag,prec);
}
GEN
bnrinit0(GEN bignf, GEN ideal, long flag, long prec)
{
switch(flag)
{
case 0: flag = nf_INIT; break;
case 1: flag = nf_INIT | nf_GEN; break;
default: err(flagerr,"bnrinit");
}
return buchrayall(bignf,ideal,flag,prec);
}
GEN
rayclassno(GEN bnf,GEN ideal)
{
GEN nf,clno,dataunit,p1;
GEN racunit,bigres,bid,resbid,resbid2,funits,hmat;
long RU,R3,i,j,av=avma;
bnf = checkbnf(bnf); nf=(GEN)bnf[7]; bigres=(GEN)bnf[8];
funits = check_units(bnf,"rayclassno");
clno = gmael(bigres,1,1);
bid = zidealstarinitall(nf,ideal,0);
resbid=(GEN)bid[2]; resbid2=(GEN)resbid[2];
R3=lg(resbid2)-1; if (!R3) { avma=av; return icopy(clno); }
RU=lg(funits); dataunit=cgetg(RU+R3+1,t_MAT);
racunit=gmael(bigres,4,2); dataunit[1]=(long)zideallog(nf,racunit,bid);
for (j=2; j<=RU; j++)
dataunit[j]=(long)zideallog(nf,(GEN)funits[j-1],bid);
for ( ; j<=RU+R3; j++)
{
p1=cgetg(R3+1,t_COL); dataunit[j]=(long)p1;
for (i=1; i<=R3; i++)
p1[i] = (i==(j-RU))?resbid2[i]:zero;
}
hmat=hnfmod(dataunit,(GEN)resbid[1]);
for (i=lg(hmat)-1 ; i; i--) clno = mulii(clno,gcoeff(hmat,i,i));
avma=av; return icopy(clno);
}
GEN
isprincipalrayall(GEN bnr, GEN x, long flag)
{
long av=avma,tetpil,i,j,c,N,ngen,ngzk;
GEN bnf,nf,bid,vecel,vecep,matu,ep,p1,p2,p3,p4,beta,idep,y,rayclass;
GEN divray,genray,alpha,alphaall,racunit,res,funit,pol;
checkbnr(bnr); bnf=(GEN)bnr[1]; bid=(GEN)bnr[2];
vecel=(GEN)bnr[3]; matu=(GEN)bnr[4];
rayclass=(GEN)bnr[5]; nf=(GEN)bnf[7]; ngen=lg(vecel)-1;
idep=isprincipalall(bnf,x,nf_GEN | nf_FORCE);
if (lg(idep[1]) != ngen+1)
err(talker,"incorrect generator length in isprincipalray");
pol=(GEN)nf[1]; N=lgef(pol)-3;
p2=cgetg(N+1,t_COL); p2[1]=un;
for (i=2; i<=N; i++) p2[i]=zero;
ep=(GEN)idep[1];
for (i=1; i<=ngen; i++)
if (typ(vecel[i]) != t_INT)
p2=element_mul(nf,p2,element_pow(nf,(GEN)vecel[i],(GEN)ep[i]));
beta=element_div(nf,(GEN)idep[2],p2);
p3=zideallog(nf,beta,bid); ngzk=lg(p3)-1;
vecep=cgetg(ngen+ngzk+1,t_COL);
for (i=1; i<=ngen; i++) vecep[i]=ep[i];
for ( ; i<=ngen+ngzk; i++) vecep[i]=p3[i-ngen];
p1=gmul(matu,vecep);
divray=(GEN)rayclass[2]; c=lg(divray);
tetpil=avma; y=cgetg(c,t_COL);
for (i=1; i<c; i++)
y[i] = lmodii((GEN)p1[i],(GEN)divray[i]);
if (!(flag & nf_GEN)) return gerepile(av,tetpil,y);
if (lg(rayclass)<=3)
err(talker,"please apply bnrinit(,,1) and not bnrinit(,,0)");
genray=(GEN)rayclass[3]; p1=idmat(N);
for (i=1; i<c; i++)
p1=idealmul(nf,idealpow(nf,(GEN)genray[i],(GEN)y[i]),p1);
alphaall = isprincipalall(bnf,idealdiv(nf,x,p1),nf_GEN | nf_FORCE);
if (!gcmp0((GEN)alphaall[1])) err(bugparier,"isprincipalray (bug1)");
res=(GEN)bnf[8];
funit = check_units(bnf,"isprincipalrayall");
racunit=(GEN)res[4];
alpha = basistoalg(nf,(GEN)alphaall[2]);
p3=zideallog(nf,(GEN)alphaall[2],bid);
if (lg(p3)>1)
{
p4=(GEN)bnr[6]; p3=gmul((GEN)p4[1],p3);
if (!gcmp1(denom(p3))) err(bugparier,"isprincipalray (bug2)");
x=lllreducemodmatrix(p3,(GEN)p4[2]);
p3=gpui(gmodulcp((GEN)racunit[2],pol),(GEN)x[1],0);
for (j=1; j<lg(funit); j++)
p3=gmul(p3,gpui(gmodulcp((GEN)funit[j],pol),(GEN)x[j+1],0));
alpha = gdiv(alpha,p3);
}
tetpil=avma; p1=cgetg(4,t_VEC);
p1[1]=lcopy(y); p1[2]=(long)algtobasis(nf,alpha);
p1[3]=lmin((GEN)idep[3],(GEN)alphaall[3]);
return gerepile(av,tetpil,p1);
}
GEN
isprincipalray(GEN bignfray, GEN x)
{
return isprincipalrayall(bignfray,x,nf_REGULAR);
}
GEN
isprincipalraygen(GEN bignfray, GEN x)
{
return isprincipalrayall(bignfray,x,nf_GEN);
}
/* DK = |dK|; give c[N][R1] */
static GEN
zimmertbound(long N,long R1,GEN DK)
{
long av,tetpil,i,R2;
GEN w,p1,minkowski;
if (N<2) return gun;
av=avma;
if (N<21)
{
double **c=(double**)gpmalloc(sizeof(double*)*21);
for (i=1; i<=20; i++) c[i]=(double*)gpmalloc(sizeof(double)*21);
c[2][2] = -0.6931; c[3][3] = -1.71733859;
c[2][0] = -0.45158; c[3][1] = -1.37420604;
c[4][4] = -2.91799837; c[5][5] = -4.22701425;
c[4][2] = -2.50091538; c[5][3] = -3.75471588;
c[4][0] = -2.11943331; c[5][1] = -3.31196660;
c[6][6] = -5.61209925; c[7][7] = -7.05406203;
c[6][4] = -5.09730381; c[7][5] = -6.50550021;
c[6][2] = -4.60693851; c[7][3] = -5.97735406;
c[6][0] = -4.14303665; c[7][1] = -5.47145968;
c[8][8] = -8.54052636; c[9][9] = -10.0630022;
c[8][6] = -7.96438858; c[9][7] = -9.46382812;
c[8][4] = -7.40555445; c[9][5] = -8.87952524;
c[8][2] = -6.86558259; c[9][3] = -8.31139202;
c[8][0] = -6.34608077; c[9][1] = -7.76081149;
c[10][10]= -11.6153797; c[11][11]= -13.1930961;
c[10][8] = -10.9966020; c[11][9] = -12.5573772;
c[10][6] = -10.3907654; c[11][7] = -11.9330458;
c[10][4] = -9.79895170; c[11][5] = -11.3210061;
c[10][2] = -9.22232770; c[11][3] = -10.7222412;
c[10][0] = -8.66213267; c[11][1] = -10.1378082;
c[12][12]= -14.7926394; c[13][13]= -16.4112395;
c[12][10]= -14.1420915; c[13][11]= -15.7475710;
c[12][8] = -13.5016616; c[13][9] = -15.0929680;
c[12][6] = -12.8721114; c[13][7] = -14.4480777;
c[12][4] = -12.2542699; c[13][5] = -13.8136054;
c[12][2] = -11.6490374; c[13][3] = -13.1903162;
c[12][0] = -11.0573775; c[13][1] = -12.5790381;
c[14][14]= -18.0466672; c[15][15]= -19.6970961;
c[14][12]= -17.3712806; c[15][13]= -19.0111606;
c[14][10]= -16.7040780; c[15][11]= -18.3326615;
c[14][8] = -16.0456127; c[15][9] = -17.6620757;
c[14][6] = -15.3964878; c[15][7] = -16.9999233;
c[14][4] = -14.7573587; c[15][5] = -16.3467686;
c[14][2] = -14.1289364; c[15][3] = -15.7032228;
c[14][0] = -13.5119848; c[15][1] = -15.0699480;
c[16][16]= -21.3610081; c[17][17]= -23.0371259;
c[16][14]= -20.6655103; c[17][15]= -22.3329066;
c[16][12]= -19.9768082; c[17][13]= -21.6349299;
c[16][10]= -19.2953176; c[17][11]= -20.9435607;
c[16][8] = -18.6214885; c[17][9] = -20.2591899;
c[16][6] = -17.9558093; c[17][7] = -19.5822454;
c[16][4] = -17.2988108; c[17][5] = -18.9131878;
c[16][2] = -16.6510652; c[17][3] = -18.2525157;
c[16][0] = -16.0131906; c[17][1] = -17.6007672;
c[18][18]= -24.7243611; c[19][19]= -26.4217792;
c[18][16]= -24.0121449; c[19][17]= -25.7021950;
c[18][14]= -23.3056902; c[19][15]= -24.9879497;
c[18][12]= -22.6053167; c[19][13]= -24.2793271;
c[18][10]= -21.9113705; c[19][11]= -23.5766321;
c[18][8] = -21.2242247; c[19][9] = -22.8801952;
c[18][6] = -20.5442836; c[19][7] = -22.1903709;
c[18][4] = -19.8719830; c[19][5] = -21.5075437;
c[18][2] = -19.2077941; c[19][3] = -20.8321263;
c[18][0] = -18.5522234; c[19][1] = -20.1645647;
c[20][20]= -28.1285704;
c[20][18]= -27.4021674;
c[20][16]= -26.6807314;
c[20][14]= -25.9645140;
c[20][12]= -25.2537867;
c[20][10]= -24.5488420;
c[20][8] = -23.8499943;
c[20][6] = -23.1575823;
c[20][4] = -22.4719720;
c[20][2] = -21.7935548;
c[20][0] = -21.1227537;
w=gexp(dbltor(c[N][R1]),6);
for (i=1; i<=20; i++) free(c[i]); free(c);
p1=gmul(gsqrt(DK,MEDDEFAULTPREC),w);
tetpil=avma; return gerepile(av,tetpil,ground(p1));
}
R2=(N-R1)>>1; p1=gdiv(mpfact(N),gpuigs(stoi(N),N));
minkowski=ground(gmul(gmul(p1,gpuigs(gdivsg(4,mppi(MEDDEFAULTPREC)),R2)),gsqrt(DK,MEDDEFAULTPREC)));
if (cmpis(minkowski,500000)>0)
err(talker,"The field has degree more than 20 and the Minkowski bound is larger than 500 000: it is unrealistic to certify it");
tetpil=avma; return gerepile(av,tetpil,gcopy(minkowski));
}
/* all primes up to Minkowski bound factor on factorbase ? */
static void
testprime(GEN bnf,GEN minkowski)
{
long av = avma, pp,i,nbideal,k,pmax;
GEN f,p1,vectpp,fb,dK, nf=checknf(bnf);
byteptr delta = diffptr;
if (DEBUGLEVEL>=2)
fprintferr("PHASE 1: check primes to Zimmert bound = %Z\n\n",minkowski);
f=(GEN)nf[4]; dK=(GEN)nf[3];
if (!gcmp1(f))
{
GEN different = gmael(nf,5,5);
if (DEBUGLEVEL>=2)
fprintferr("**** Testing Different = %Z\n",different);
p1 = isprincipalall(bnf,different,nf_FORCE);
if (DEBUGLEVEL>=2) fprintferr(" is %Z\n",p1);
}
fb=(GEN)bnf[5];
p1 = gmael(fb, lg(fb)-1, 1); /* largest p in factorbase */
pp = 0; pmax = is_bigint(p1)? VERYBIGINT: itos(p1);
while (cmpsi(pp,minkowski)<1)
{
pp += *delta++; if (!*delta) err(primer1);
if (DEBUGLEVEL>=2) fprintferr("*** p = %ld\n",pp);
vectpp=primedec(bnf,stoi(pp)); nbideal=lg(vectpp)-1;
/* loop through all P | p if ramified, all but one otherwise */
if (!smodis(dK,pp)) nbideal++;
for (i=1; i<nbideal; i++)
{
GEN P = (GEN)vectpp[i]; /* non inert */
if (DEBUGLEVEL>=2)
fprintferr(" Testing P = %Z\n",P);
if (cmpii(idealnorm(bnf,P),minkowski) < 1)
{
if (pp <= pmax && (k = tablesearch(fb, P, cmp_prime_ideal)))
{
if (DEBUGLEVEL>=2) fprintferr(" #%ld in factor base\n",k);
}
else
{
p1=isprincipalall(bnf,P,nf_FORCE);
if (DEBUGLEVEL>=2) fprintferr(" is %Z\n",p1);
}
}
else if (DEBUGLEVEL>=2)
fprintferr(" Norm(P) > Zimmert bound\n");
}
}
avma=av;
if (DEBUGLEVEL>=2) { fprintferr("End of PHASE 1.\n\n"); flusherr(); }
}
/* rend constante d'Hermite^n si connue, sinon une borne sup */
static GEN
hermiteconstant(long n)
{
long av,tetpil;
GEN h,h1;
switch(n)
{
case 1: return gun;
case 2: h=cgetg(3,t_FRAC); h[1]=lstoi(4); h[2]=lstoi(3); return h;
case 3: return gdeux;
case 4: return stoi(4);
case 5: return stoi(8);
case 6: h=cgetg(3,t_FRAC); h[1]=lstoi(64); h[2]=lstoi(3); return h;
case 7: return stoi(64);
case 8: return stoi(256);
}
av = avma;
h = gpuigs(gdiv(gdeux,mppi(DEFAULTPREC)),n);
h1 = gsqr(ggamma(gdivgs(stoi(n+4),2),DEFAULTPREC));
tetpil=avma; return gerepile(av,tetpil,gmul(h,h1));
}
/* 1 primitif, 0 s'il est peut etre imprimitif... */
static long
isprimitive(GEN nf)
{
long N,first,i,l,ep;
GEN d,fa;
N = lgef(nf[1])-3; fa = (GEN)factor(stoi(N))[1]; /* primes | N */
first = itos((GEN)fa[1]); if (first==N) return 1;
d=absi((GEN)nf[3]); fa=(GEN)factor(d)[2]; /* expo. primes | disc(nf) */
if (mod2(d))
{ i=1; ep=1; }
else
{ i=2; ep=itos((GEN)fa[1])>>1; }
l=lg(fa);
for ( ; i < l; i++)
{
if (ep >= first) return 0;
ep = itos((GEN)fa[i]);
}
return 1;
}
static GEN
regulatorbound(GEN bnf)
{
long N,R1,R2,R;
GEN nf,dKa,bound,p1,c1;
nf=(GEN)bnf[7]; N=lgef(nf[1])-3;
bound=dbltor(0.2);
if (!isprimitive(nf))
{
if (DEBUGLEVEL>=2)
{ fprintferr("Default bound for regulator: 0.2\n"); flusherr(); }
return bound;
}
dKa=absi((GEN)nf[3]);
R1=itos(gmael(nf,2,1));
R2=itos(gmael(nf,2,2)); R=R1+R2-1;
if (!R2 && N<12) c1=gpuigs(stoi(4),N>>1); else c1=gpuigs(stoi(N),N);
if (cmpii(dKa,c1)<=0)
{
if (DEBUGLEVEL>=2)
{ fprintferr("Default bound for regulator: 0.2\n"); flusherr(); }
return bound;
}
p1 = gsqr(glog(gdiv(dKa,c1),DEFAULTPREC));
p1 = gdivgs(gmul2n(gpuigs(gdivgs(gmulgs(p1,3),N*(N*N-1)-6*R2),R),R2),N);
p1 = gsqrt(gdiv(p1, hermiteconstant(R)), DEFAULTPREC);
if (gcmp(p1,bound) > 0) bound = p1;
if (DEBUGLEVEL>=2)
{ fprintferr("Mahler bound for regulator: "); outerr(p1); flusherr(); }
return bound;
}
#define NBMAX 5000
/* should use smallvectors */
static GEN
minimforunits(GEN nf, long borne, long stockmax)
{
long av = avma,av1,n1,n,i,j,k,s,norme,normax,*x,fl1,cmpt;
GEN u,r,S,S1,a,base,p1;
double p;
double **q,*v,*y,*z;
double eps=0.000001;
if (DEBUGLEVEL>=2)
{
fprintferr("Searching minimum of T2-form on units:\n");
if (DEBUGLEVEL>2) fprintferr(" borne = %ld\n",borne);
flusherr();
}
a=gmael(nf,5,3); n1=lg(a);
n=n1-1;
x=(long*)gpmalloc(n1*sizeof(long));
y=(double*)gpmalloc(n1*sizeof(double));
z=(double*)gpmalloc(n1*sizeof(double));
v=(double*)gpmalloc(n1*sizeof(double));
q=(double**)gpmalloc(n1*sizeof(double*));
for (j=1; j<=n; j++) q[j]=(double*)gpmalloc(n1*sizeof(double));
u=lllgram(a,BIGDEFAULTPREC); base=gmul((GEN)nf[7],u);
a=gmul(qf_base_change(a,u,1), realun(BIGDEFAULTPREC));
r=sqred1(a);
for (j=1; j<=n; j++)
{
v[j]=rtodbl(gcoeff(r,j,j));
for (i=1; i<j; i++)
q[i][j]=rtodbl(gcoeff(r,i,j));
}
normax=0;
if (stockmax) S=cgetg(stockmax+1,t_MAT);
s=0; k=n; cmpt=0; y[n]=z[n]=0;
x[n]=(long)(sqrt(borne/v[n]+eps));
for(;;)
{
do
{
if (k>1)
{
k--; z[k]=0;
for (j=k+1; j<=n; j++) z[k]=z[k]+q[k][j]*x[j];
p=x[k+1]+z[k+1];
y[k]=y[k+1]+p*p*v[k+1];
x[k]=(long)floor(sqrt((borne-y[k]+eps)/v[k])-z[k]);
}
while (v[k]*(x[k]+z[k])*(x[k]+z[k]) > borne-y[k]+eps)
{
k++; x[k]--;
}
}
while (k>1);
if (!x[1] && y[1]<=eps) break;
cmpt++;
if (cmpt>NBMAX)
{
free(x); free(y); free(z); free(v);
for (j=1; j<=n; j++) free(q[j]); free(q);
av=avma; return NULL;
}
if (DEBUGLEVEL>8){ fprintferr("."); flusherr(); }
norme=(long)(y[1]+v[1]*(x[1]+z[1])*(x[1]+z[1])+eps);
if (norme>normax) normax=norme;
av1=avma; p1=gzero;
for (i=1; i<=n; i++) p1=gadd(p1,gmulsg(x[i],(GEN)base[i]));
fl1=gcmp1(gabs(subres(p1,(GEN)nf[1]),0)); avma=av1;
if (fl1)
{
if (DEBUGLEVEL>=2) { fprintferr("*"); flusherr(); }
s++; cmpt=0;
if (s<=stockmax)
{
p1=cgetg(n+1,t_COL);
for (i=1; i<=n; i++) p1[i]=lstoi(x[i]);
S[s]=lmul(u,p1);
}
}
x[k]--;
}
free(x); free(y); free(z); free(v);
for (j=1; j<=n; j++) free(q[j]); free(q);
if (DEBUGLEVEL>=2){ fprintferr("\n"); flusherr(); }
if (stockmax)
{
av1=avma;
k=(s<stockmax)? s:stockmax;
S1=cgetg(k+1,t_MAT);
for (j=1; j<=k; j++) S1[j]=lcopy((GEN)S[j]);
S=gerepile(av,av1,S1);
}
else { avma=av; S=cgetg(1,t_MAT); }
u=cgetg(4,t_VEC);
u[1]=lstoi(s<<1);
u[2]=lstoi(normax);
u[3]=(long)S;
return u;
}
#undef NBMAX
static GEN
compute_M0(GEN M_star,long N) /* On connait M_star; on calcule M0 */
{
long av1,tetpil,m1,m2,n1,n2,n3,k,kk,lr,lr1,lr2,i,j,l,vx,vy,vz,vM,PREC,prec;
GEN eps,pol,p1,p2,p3,p4,p5,p6,u,v,w,r,r1,r2,M0,M0_pro,S,P,M_formel;
GEN f1,f2,f3,g1,g2,g3,pg1,pg2,pg3,pf1,pf2,pf3,p7,p8,p9,p10,p11;
GEN x,y,z;
PREC=gprecision(M_star);
if (N==2) return gmul2n(gsqr(gach(gmul2n(M_star,-1),PREC)), -1);
vM = fetch_var(); M_formel=polx[vM];
vz = fetch_var(); z=polx[vz];
vy = fetch_var(); y=polx[vy];
vx = fetch_var(); x=polx[vx];
vx=varn(x); vy=varn(y); vz=varn(z); vM=varn(M_formel);
PREC = PREC>>1; if (!PREC) PREC=DEFAULTPREC;
eps=dbltor(0.0000001); M0=gzero; m1=(N-(N%3))/3;
for (n1=1; n1<=m1; n1++)
{
m2 = (N-n1)>>1;
for (n2=n1; n2<=m2; n2++)
{
av1=avma; n3=N-n1-n2; prec=PREC;
if (!(N%3) && n1==n2 && n1==n3)
{
p1=gdivgs(M_star,m1); p2=gaddsg(1,p1); p3=gsubgs(p1,3);
p4=gsqrt(gmul(p2,p3),prec); p5=gsubgs(p1,1);
u=gun; v=gmul2n(gadd(p5,p4),-1); w=gmul2n(gsub(p5,p4),-1);
M0_pro=gmul2n(gmulsg(m1,gadd(gsqr(glog(v,prec)),gsqr(glog(w,prec)))),-2);
if (DEBUGLEVEL>2)
{
fprintferr("[%ld,%ld,%ld]: ",n1,n2,n3);
outerr(M0_pro); flusherr();
}
if (gcmp0(M0) || gcmp(M0_pro,M0)<0)
{
M0=M0_pro; tetpil=avma; M0=gerepile(av1,tetpil,gcopy(M0));
}
else avma=av1;
}
else if (n1==n2 || n1==n3 || n2==n3)
{
if (n1==n2) k=n1; else if (n2==n3) k=n3;
kk=N-2*k;
p2=gsub(M_star,gmulgs(x,k));
p3=gmul(gpuigs(stoi(kk),kk),gpuigs(gsubgs(gmul(M_star,p2),kk*kk),k));
pol=gsub(p3,gmul(gmul(gpuigs(stoi(k),k),gpuigs(x,k)),gpuigs(p2,N-k)));
prec=gprecision(pol); if (!prec) prec=5;
r=roots(pol,prec); lr=lg(r)-1;
for (i=1; i<=lr; i++)
{
if (gcmp(gabs(gimag((GEN)r[i]),prec),eps) < 0 &&
gsigne(S=greal((GEN)r[i])) > 0)
{
p4=gsub(M_star,gmulsg(k,S));
P=gdiv(gmul(gmulsg(k,S),p4),gsubgs(gmul(M_star,p4),kk*kk));
p5=gsub(gsqr(S),gmul2n(P,2));
if (gsigne(p5)>=0)
{
p6=gsqrt(p5,prec);
u=gmul2n(gadd(S,p6),-1); v=gmul2n(gsub(S,p6),-1);
if ((gsigne(u)>0)&&(gsigne(v)>0))
{
w=gpui(P,gdivgs(stoi(-k),kk),prec);
p6=gmulsg(k,gadd(gsqr(glog(u,prec)),gsqr(glog(v,prec))));
M0_pro=gmul2n(gadd(p6,gmulsg(kk,gsqr(glog(w,prec)))),-2);
if (DEBUGLEVEL>2)
{
fprintferr("[%ld,%ld,%ld] : ",n1,n2,n3);
outerr(M0_pro); flusherr();
}
if (gcmp0(M0) || gcmp(M0_pro,M0)<0) M0=M0_pro;
}
}
}
}
tetpil=avma; M0=gerepile(av1,tetpil,gcopy(M0));
}
else
{
f1=gadd(gmulsg(n1,x),gadd(gmulsg(n2,y),gmulsg(n3,z)));
f1=gsub(f1,M_formel);
f2=gmulsg(n1,gmul(y,z));
f2=gadd(f2,gmulsg(n2,gmul(x,z)));
f2=gadd(f2,gmulsg(n3,gmul(x,y)));
f2=gsub(f2,gmul(M_formel,gmul(x,gmul(y,z))));
f3=gsub(gmul(gpuigs(x,n1),gmul(gpuigs(y,n2),gpuigs(z,n3))),gun);
g1=subres(f1,f2); g1=gdiv(g1,content(g1));
g2=subres(f1,f3); g2=gdiv(g2,content(g2));
g3=subres(g1,g2); g3=gdiv(g3,content(g3));
pf1=gsubst(f1,vM,M_star); pg1=gsubst(g1,vM,M_star);
pf2=gsubst(f2,vM,M_star); pg2=gsubst(g2,vM,M_star);
pf3=gsubst(f3,vM,M_star); pg3=gsubst(g3,vM,M_star);
prec=gprecision(pg3); if (!prec) prec=5;
r=roots(pg3,prec); lr=lg(r)-1;
for (i=1; i<=lr; i++)
{
if (gcmp(gabs(gimag((GEN)r[i]),prec),eps) < 0 &&
gsigne(w=greal((GEN)r[i])) > 0)
{
p1=gsubst(pg1,vz,w); p2=gsubst(pg2,vz,w);
p3=gsubst(pf1,vz,w); p4=gsubst(pf2,vz,w); p5=gsubst(pf3,vz,w);
prec=gprecision(p1); if (!prec) prec=5;
r1=roots(p1,prec); lr1=lg(r1)-1;
for (j=1; j<=lr1; j++)
{
if (gcmp(gabs(gimag((GEN)r1[j]),prec),eps) < 0 &&
gsigne(v=greal((GEN)r1[j])) > 0)
{
p6=gsubst(p2,vy,v);
if (gcmp(gabs(p6,prec),eps)<0)
{
p7=gsubst(p3,vy,v); p8=gsubst(p4,vy,v); p9=gsubst(p5,vy,v);
prec=gprecision(p7); if (!prec) prec=5;
r2=roots(p7,prec); lr2=lg(r2)-1;
for (l=1; l<=lr2; l++)
{
if (gcmp(gabs(gimag((GEN)r2[l]),prec),eps) < 0 &&
gsigne(u=greal((GEN)r2[l])) > 0)
{
p10=gsubst(p8,vx,u);
if (gcmp(gabs(p10,prec),eps)<0)
{
p11=gsubst(p9,vx,u);
if (gcmp(gabs(p11,prec),eps)<0)
{
M0_pro=gmulsg(n1,gsqr(glog(u,prec)));
M0_pro=gadd(M0_pro,gmulsg(n2,gsqr(glog(v,prec))));
M0_pro=gadd(M0_pro,gmulsg(n3,gsqr(glog(w,prec))));
M0_pro=gmul2n(M0_pro,-2);
if (DEBUGLEVEL>2)
{
fprintferr("[ %ld,%ld,%ld ] : ",n1,n2,n3);
outerr(M0_pro); flusherr();
}
if (gcmp0(M0) || gcmp(M0_pro,M0) < 0) M0=M0_pro;
}
}
}
}
}
}
}
}
}
tetpil=av1; M0_pro=gerepile(av1,tetpil,gcopy(M0));
}
}
}
for (i=1;i<=4;i++) delete_var();
return M0;
}
static GEN
lowerboundforregulator(GEN bnf,GEN units)
{
long N,R1,R2,RU,i,nbrootsofone,nbmin;
GEN rootsofone,nf,M0,M,m,col,T2,bound,minunit,newminunit;
GEN vecminim,colalg,p1,pol,y;
rootsofone=gmael(bnf,8,4); nbrootsofone=itos((GEN)rootsofone[1]);
nf=(GEN)bnf[7]; T2=gmael(nf,5,3); N=lgef(nf[1])-3;
R1=itos(gmael(nf,2,1)); R2=itos(gmael(nf,2,2)); RU=R1+R2-1;
if (RU==0) return gun;
units=algtobasis(bnf,units); minunit=qfeval(T2,(GEN)units[1]);
for (i=2; i<=RU; i++)
{
newminunit=qfeval(T2,(GEN)units[i]);
if (gcmp(newminunit,minunit)<0) minunit=newminunit;
}
if (gcmpgs(minunit,1000000000)>0) return regulatorbound(bnf);
vecminim=minimforunits(nf,itos(gceil(minunit)),10000);
if (!vecminim) return regulatorbound(bnf);
m=(GEN)vecminim[3]; nbmin=lg(m)-1;
if (nbmin==10000) return regulatorbound(bnf);
bound=gaddgs(minunit,1);
for (i=1; i<=nbmin; i++)
{
col=(GEN)m[i]; colalg=basistoalg(nf,col);
if (!gcmp1(lift_intern(gpuigs(colalg,nbrootsofone))))
{
newminunit=qfeval(T2,col);
if (gcmp(newminunit,bound)<0) bound=newminunit;
}
}
if (gcmp(bound,minunit)>0) err(talker,"bug in lowerboundforregulator");
if (DEBUGLEVEL>=2)
{
fprintferr("M* = "); outerr(bound); flusherr();
if (DEBUGLEVEL>2)
{
p1=polx[0]; pol=gaddgs(gsub(gpuigs(p1,N),gmul(bound,p1)),N-1);
fprintferr("pol = "); outerr(pol); flusherr();
p1=roots(pol,DEFAULTPREC);
if (N&1) y=greal((GEN)p1[3]); else y=greal((GEN)p1[2]);
fprintferr("y = "); outerr(y); flusherr();
M0=gmul2n(gmulsg(N*(N-1),gsqr(glog(y,DEFAULTPREC))),-2);
fprintferr("old fashion M0 = "); outerr(M0); flusherr();
}
}
M0=compute_M0(bound,N);
if (DEBUGLEVEL>=2) { fprintferr("M0 = "); outerr(M0); flusherr(); }
M=gmul2n(gdivgs(gdiv(gpuigs(M0,RU),hermiteconstant(RU)),N),R2);
if (gcmp(M,dbltor(0.04))<0) return regulatorbound(bnf);
M = gsqrt(M,DEFAULTPREC);
if (DEBUGLEVEL>=2)
{ fprintferr("(lower bound for regulator) M = "); outerr(M); flusherr(); }
return M;
}
/* Calcule une matrice carree de rang lg(beta) associee a une famille
* d'ideaux premiers P_i tels que : 1<=i<=nombre de beta; N(P_i) congru a 1
* mod pp v_(P_i)(beta[j])=0 pour tout 1<=j<=nbalphapro et 1<=i<=lg(beta)
*/
static void
primecertify(GEN bnf,GEN beta,long pp,GEN big)
{
long i,j,qq,nbcol,sizeofmat,nbqq,ra,N;
GEN nf,mat,mat1,qgen,decqq,newcol,eltgen,qrhall,Q;
nbcol=0; nf=(GEN)bnf[7]; N=lgef(nf[1])-3;
sizeofmat=lg(beta)-1; mat=cgetg(1,t_MAT); qq=1;
for(;;)
{
qq += 2*pp; qgen=stoi(qq);
if (smodis(big,qq)==0 || !isprime(qgen)) continue;
decqq=primedec(bnf,qgen); nbqq=lg(decqq)-1;
for (i=1; i<=nbqq; i++)
{
Q=(GEN)decqq[i];
if (!gcmp1((GEN)Q[4])) continue;
qrhall=nfmodprinit(nf,Q); nbcol++;
newcol=cgetg(sizeofmat+1,t_COL);
if (DEBUGLEVEL>=2)
fprintferr(" prime ideal Q: %Z\n",Q);
eltgen = gscalcol_i(lift(gener(qgen)), N);
for (j=1; j<=sizeofmat; j++)
newcol[j]=(long)nfshanks(nf,(GEN)beta[j],eltgen,Q,qrhall);
if (DEBUGLEVEL>=2)
{
fprintferr(" generator of (Zk/Q)^*: "); outerr(eltgen);
fprintferr(" column #%ld of the matrix log(b_j/Q): ",nbcol);
outerr(newcol);
}
mat1=concatsp(mat,newcol); ra=rank(mat1);
if (DEBUGLEVEL>=2)
{
fprintferr(" new rank of the matrix: %ld\n\n",ra); flusherr();
}
if (ra!=nbcol) nbcol--;
else
{
if (nbcol==sizeofmat) return;
mat=mat1;
}
}
}
}
static void
check_prime(long p, GEN bnf, GEN h, GEN cyc, long R, GEN alpha,
GEN funits, GEN rootsofone, GEN big)
{
long av = avma, nbpro,nbalpha,i, nbgen = lg(cyc)-1;
GEN p1,beta, nf = (GEN)bnf[7], nbrootsofone = (GEN)rootsofone[1];
if (DEBUGLEVEL>=2)
{ fprintferr("***** Testing prime p = %ld\n",p); flusherr(); }
if (smodis(h,p)) nbpro=0;
else
{
if (DEBUGLEVEL>=2) { fprintferr(" p divides cl(k)\n"); flusherr(); }
for (nbpro=nbgen; nbpro; nbpro--)
if (!smodis((GEN)cyc[nbpro],p)) break;
}
nbalpha=nbpro+R;
if (smodis(nbrootsofone,p)) p1 = (GEN) alpha[nbpro];
else
{
if (DEBUGLEVEL>=2) { fprintferr(" p divides w(k)\n"); flusherr(); }
nbalpha++; nbpro++; p1 = algtobasis(nf,(GEN)rootsofone[2]);
}
if (DEBUGLEVEL>=2)
{ fprintferr(" t+r+e = %ld\n",nbalpha); flusherr(); }
beta=cgetg(nbalpha+1,t_VEC);
if (nbpro)
{
for (i=1; i<nbpro; i++) beta[i]=alpha[i];
beta[nbpro]=(long) p1;
}
for (i=1; i<=R; i++)
beta[i+nbpro]=(long)algtobasis(nf,(GEN)funits[i]);
if (DEBUGLEVEL>=2)
{ fprintferr(" Beta list = "); outerr(beta); flusherr(); }
primecertify(bnf,beta,p,big); avma=av;
}
long
certifybuchall(GEN bnf)
{
long av = avma, nbgen,i,j,pp,N,R1,R2,R,nfa,nbf1,bound;
GEN big,nf,reg,rootsofone,funits,gen,p1,gbound,DK,alpha,factfd1,f1,h,cyc;
byteptr delta = diffptr;
bnf = checkbnf(bnf); nf = (GEN)bnf[7];
N=lgef(nf[1])-3; if (N==1) return 1;
R1=itos(gmael(nf,2,1)); R2=itos(gmael(nf,2,2)); R=R1+R2-1;
funits = check_units(bnf,"bnfcertify");
DK=absi((GEN)nf[3]);
testprime(bnf,zimmertbound(N,R1,DK));
p1=gmael(bnf,8,1); reg=gmael(bnf,8,2);
h=(GEN)p1[1];
cyc=(GEN)p1[2]; nbgen=lg(cyc)-1;
gen=(GEN)p1[3]; rootsofone=gmael(bnf,8,4);
if (DEBUGLEVEL>1)
{
fprintferr("Class number = "); outerr(h);
fprintferr("Cyclic components = "); outerr(cyc);
fprintferr("Generators = "); outerr(gen);
fprintferr("Regulator = "); outerr(reg);
fprintferr("Roots of one = "); outerr(rootsofone);
fprintferr("Fundamental units = "); outerr(funits);
}
alpha=cgetg(nbgen+1,t_VEC);
for (i=1; i<=nbgen; i++)
{
p1=idealpow(nf,(GEN)gen[i],(GEN)cyc[i]);
alpha[i]=isprincipalall(bnf,p1,nf_GEN | nf_FORCE)[2];
}
gbound = ground(gdiv(reg,lowerboundforregulator(bnf,funits)));
if (is_bigint(gbound))
err(talker,"sorry, too many primes to check");
bound = gbound[2];
if (DEBUGLEVEL>=2)
{
fprintferr("\nPHASE 2: are all primes good ?\n\n");
fprintferr(" Testing primes <= B (= %ld)\n\n",bound); flusherr();
}
for (big=gun,i=1; i<=nbgen; i++)
big = mulii(big,idealnorm(nf,(GEN)gen[i]));
for (pp = *delta++; pp <= bound; pp += *delta++)
check_prime(pp,bnf,h,cyc,R,alpha,funits,rootsofone,big);
nfa = 0;
if (nbgen)
{
factfd1=factor((GEN)cyc[1]);
nbf1=lg(factfd1[1]); f1=(GEN)factfd1[1];
for (i=1; i<nbf1; i++)
if (cmpis((GEN)f1[i],bound) > 0) nfa++;
}
if (DEBUGLEVEL>=2 && nfa)
{ fprintferr(" Testing primes > B (# = %ld)\n\n",nfa); flusherr(); }
for (j=1; j<=nfa; j++)
{
pp = itos((GEN)f1[nbf1-j]);
check_prime(pp,bnf,gzero,cyc,R,alpha,funits,rootsofone,big);
}
avma=av; return 1;
}
/*******************************************************************/
/* */
/* CORPS DE CLASSES DE RAYON : CONDUCTEURS ET DISCRIMINANTS */
/* */
/*******************************************************************/
/* Si s est la surjection de Cl_m sur Cl_n et H ssgroupe de Cl_m,
* retourne le ssgroupe s(H) de Cl_n
*/
static GEN
imageofgroup0(GEN gen,GEN bnr,GEN subgroup)
{
long j,l;
GEN E,Delta = diagonal(gmael(bnr,5,2));
if (gcmp0(subgroup)) return Delta;
l=lg(gen); E=cgetg(l,t_MAT);
for (j=1; j<l; j++)
E[j] = (long)isprincipalray(bnr,(GEN)gen[j]);
E = concatsp(gmul(E,subgroup),Delta);
return hnf(E);
}
static GEN
imageofgroup(GEN gen,GEN bnr,GEN subgroup)
{
long av = avma;
return gerepileupto(av,imageofgroup0(gen,bnr,subgroup));
}
/* retourne le cardinal de Cl_n / s(H) */
static GEN
cardofimagofquotientgroup(GEN H,GEN bnr2,GEN subgroup)
{
return dethnf_i(imageofgroup0(H,bnr2,subgroup));
}
static GEN
args_to_bnr(GEN arg0, GEN arg1, GEN arg2, GEN *subgroup, long prec)
{
GEN bnr,bnf;
if (typ(arg0)!=t_VEC)
err(talker,"neither bnf nor bnr in conductor or discray");
if (!arg1) arg1 = gzero;
if (!arg2) arg2 = gzero;
switch(lg(arg0))
{
case 7: /* bnr */
bnr=arg0; (void)checkbnf((GEN)bnr[1]);
*subgroup=arg1; break;
case 11: /* bnf */
bnf = checkbnf(arg0);
bnr=buchrayall(bnf,arg1,nf_INIT | nf_GEN,prec);
*subgroup=arg2; break;
default: err(talker,"neither bnf nor bnr in conductor or discray");
}
if (!gcmp0(*subgroup))
{
long tx=typ(*subgroup);
if (!is_matvec_t(tx))
err(talker,"bad subgroup in conductor or discray");
}
return bnr;
}
GEN
bnrconductor(GEN arg0,GEN arg1,GEN arg2,long all,long prec)
{
GEN sub=arg1, bnr=args_to_bnr(arg0,arg1,arg2,&sub,prec);
return conductor(bnr,sub,all,prec);
}
long
bnrisconductor(GEN arg0,GEN arg1,GEN arg2,long prec)
{
GEN sub=arg1, bnr=args_to_bnr(arg0,arg1,arg2,&sub,prec);
return itos(conductor(bnr,sub,-1,prec));
}
/* Given a number field bnf=bnr[1], a ray class group bnr (from buchrayinit),
* and a subgroup (HNF form) of the ray class group, compute the conductor
* of the subgroup (copy of discrayrelall) if all=0. If all > 0, compute
* furthermore the corresponding subgroup and output
* [[ideal,arch],[hm,cyc,gen],subgroup] if all = 1, and output
* [[ideal,arch],newbnr,subgroup] if all = 2. If all<0, answer only 1 is
* module is the conductor, 0 otherwise.
*/
GEN
conductor(GEN bnr,GEN subgroup,long all,long prec)
{
long r1,j,av=avma,tetpil,k,ep,trivial=0;
GEN bnf,nf,cl,cyc,gen,bid,ideal,arch,p1,p2,clhray,clhss;
GEN fa,arch2,bnr2,fa2,ex;
checkbnrgen(bnr); bnf=(GEN)bnr[1]; bid=(GEN)bnr[2];
cl=(GEN)bnr[5]; cyc=(GEN)cl[2]; gen=(GEN)cl[3];
nf=(GEN)bnf[7]; r1=itos(gmael(nf,2,1));
p1=(GEN)bid[1]; ideal=(GEN)p1[1]; arch=(GEN)p1[2];
if (gcmp0(subgroup)) { trivial=1; clhray=(GEN)cl[1]; }
else
{
p1 = gauss(subgroup,diagonal(cyc));
if (!gcmp1(denom(p1)))
err(talker,"incorrect subgroup in conductor");
if (gcmp1(det(p1))) trivial=1;
clhray = absi(det(subgroup));
}
fa=(GEN)bid[3]; fa2=(GEN)fa[1]; ex=(GEN)fa[2];
p2=cgetg(3,t_VEC); p2[2]=(long)arch;
for (k=1; k<lg(fa2); k++)
{
GEN pr=(GEN)fa2[k], prinv=idealinv(nf,pr);
ep = (all>=0)? itos((GEN)ex[k]): 1;
for (j=1; j<=ep; j++)
{
p2[1]=(long)idealmul(nf,ideal,prinv);
if (trivial) clhss=rayclassno(bnf,p2);
else
{
bnr2=buchrayall(bnf,p2,nf_INIT,prec);
clhss=cardofimagofquotientgroup(gen,bnr2,subgroup);
}
if (!egalii(clhss,clhray)) break;
if (all<0) { avma=av; return gzero; }
ideal = (GEN)p2[1];
}
}
p2[1]=(long)ideal; arch2=dummycopy(arch); p2[2]=(long)arch2;
for (k=1; k<=r1; k++)
if (signe(arch[k]))
{
arch2[k]=zero;
if (trivial) clhss=rayclassno(bnf,p2);
else
{
bnr2=buchrayall(bnf,p2,nf_INIT,prec);
clhss=cardofimagofquotientgroup(gen,bnr2,subgroup);
}
if (!egalii(clhss,clhray)) arch2[k]=un;
else if (all<0) { avma=av; return gzero; }
}
if (all<0) {avma=av; return gun;}
if (!all) { tetpil=avma; return gerepile(av,tetpil,gcopy(p2)); }
bnr2=buchrayall(bnf,p2,nf_INIT | nf_GEN,prec);
tetpil=avma; p1=cgetg(4,t_VEC);
p1[3]=(long)imageofgroup(gen,bnr2,subgroup);
if (all==1) bnr2=(GEN)bnr2[5];
p1[2]=lcopy(bnr2);
p1[1]=lcopy(p2); return gerepile(av,tetpil,p1);
}
/* etant donne un bnr et un polynome relatif, trouve le groupe des normes
correspondant a l'extension relative en supposant qu'elle est abelienne
et que le module donnant bnr est multiple du conducteur (tout ceci n'etant
bien sur pas verifie). */
GEN
rnfnormgroup(GEN bnr, GEN polrel)
{
long av=avma,i,j,reldeg,sizemat,prime,nfac,k;
GEN bnf,polreldisc,nf,raycl,group,detgroup,fa,pr,famo,ep,fac,col,p1;
byteptr d = diffptr;
checkbnr(bnr); bnf=(GEN)bnr[1]; raycl=(GEN)bnr[5];
if (typ(polrel)!=t_POL) err(typeer,"rnfnormgroup");
reldeg=lgef(polrel)-3; detgroup=(GEN)raycl[1];
group = diagonal((GEN)raycl[2]);
k = cmpis(detgroup,reldeg);
if (k<0) err(talker,"not an Abelian extension in rnfnormgroup?");
if (!k) return group;
polreldisc=discsr(polrel); nf=(GEN)bnf[7];
sizemat=lg(group)-1; prime = *d++;
/* tant que nffactormod est bugge pour p=2 on commence a prime = 3 */
for(;;)
{
prime += *d++; if (!*d) err(primer1);
fa=primedec(nf,stoi(prime));
for (i=1; i<lg(fa); i++)
{
pr = (GEN)fa[i];
if (element_val(nf,polreldisc,pr)==0)
{
famo=nffactormod(nf,polrel,pr);
ep=(GEN)famo[2]; fac=(GEN)famo[1];
nfac=lg(ep)-1; k=lgef((GEN)fac[1])-3;
for (j=1; j<=nfac; j++)
{
if (!gcmp1((GEN)ep[j])) err(bugparier,"rnfnormgroup");
if (lgef(fac[j])-3 != k)
err(talker,"non Galois extension in rnfnormgroup");
}
col=gmulsg(k,isprincipalrayall(bnr,pr,nf_REGULAR));
p1=cgetg(sizemat+2,t_MAT);
for (j=1; j<=sizemat; j++) p1[j]=group[j];
p1[sizemat+1]=(long)col;
group=hnf(p1); detgroup=dethnf(group);
k = cmpis(detgroup,reldeg);
if (k<0) err(talker,"not an Abelian extension in rnfnormgroup?");
if (!k) { cgiv(detgroup); return gerepileupto(av,group); }
}
}
}
}
/* Etant donne un bnf et un polynome relatif polrel definissant une extension
abelienne (ce qui n'est pas verifie), calcule le conducteur et le groupe de
congruence correspondant a l'extension definie par polrel sous la forme
[[ideal,arch],[hm,cyc,gen],group] ou [ideal,arch] est le conducteur, et
[hm,cyc,gen] est le groupe de classes de rayon correspondant. */
GEN
rnfconductor(GEN bnf, GEN polrel, long prec)
{
long av=avma,tetpil,R1,i,v;
GEN nf,module,arch,bnr,group,p1,pol2;
bnf = checkbnf(bnf); nf=(GEN)bnf[7];
module=cgetg(3,t_VEC); R1=itos(gmael(nf,2,1));
v=varn(polrel);
p1=unifpol((GEN)bnf[7],polrel,0);
p1=denom(gtovec(p1));
pol2=gsubst(polrel,v,gdiv(polx[v],p1));
pol2=gmul(pol2,gpuigs(p1,degree(pol2)));
module[1]=rnfdiscf(nf,pol2)[1]; arch=cgetg(R1+1,t_VEC);
module[2]=(long)arch; for (i=1; i<=R1; i++) arch[i]=un;
bnr=buchrayall(bnf,module,nf_INIT | nf_GEN,prec);
group=rnfnormgroup(bnr,pol2); tetpil=avma;
return gerepile(av,tetpil,conductor(bnr,group,1,prec));
}
/* Etant donnes un corps de nombres bnf=bnr[1], un groupe de classes de rayon
* bnr calcule par buchrayinit(), et une matrice HNF subgroup definissant un
* sous-groupe du groupe de classes de rayon, calcule [n,r1,dk] associe a ce
* sous groupe (cf. discrayall()). si flcond=0 le calcul est arrete si le
* module n'est pas le conducteur et le programme retourne gzero. Si
* flrel=0, seul la norme de l'ideal dk est calculee au lieu de dk lui-meme
*/
static GEN
discrayrelall(GEN bnr,GEN subgroup,long flag,long prec)
{
long r1,degk,j,av=avma,tetpil,k,ep,nbdezero;
long trivial=0, flrel = flag & nf_REL, flcond = flag & nf_COND;
GEN bnf,nf,cyc,gen,bid,ideal,arch,p1,p2,clhray,clhss;
GEN fa,dlk,arch2,bnr2,idealrel,fa2,ex,som;
checkbnrgen(bnr); bnf=(GEN)bnr[1];
cyc=gmael(bnr,5,2); gen=gmael(bnr,5,3);
nf=(GEN)bnf[7]; r1=itos(gmael(nf,2,1));
if (gcmp0(subgroup)) { trivial=1; clhray=gmael(bnr,5,1); }
else
{
p1=gauss(subgroup,diagonal(cyc));
if (!gcmp1(denom(p1)))
err(talker,"incorrect subgroup in discray");
if (gcmp1(det(p1))) trivial=1;
clhray = det(subgroup);
}
bid=(GEN)bnr[2]; ideal=gmael(bid,1,1); arch=gmael(bid,1,2);
fa=(GEN)bid[3]; fa2=(GEN)fa[1]; ex=(GEN)fa[2];
degk=lgef(nf[1])-3;
idealrel=flrel?idmat(degk):gun;
p2=cgetg(3,t_VEC); p2[2]=(long)arch;
for (k=1; k<lg(fa2); k++)
{
GEN pr=(GEN)fa2[k], prinv=idealinv(nf,pr);
ep=itos((GEN)ex[k]); p2[1]=(long)ideal; som=gzero;
for (j=1; j<=ep; j++)
{
p2[1]=(long)idealmul(nf,(GEN)p2[1],prinv);
if (trivial) clhss=rayclassno(bnf,p2);
else
{
bnr2=buchrayall(bnf,p2,nf_INIT,prec);
clhss=cardofimagofquotientgroup(gen,bnr2,subgroup);
}
if (j==1 && egalii(clhss,clhray) && flcond) { avma=av; return gzero; }
if (is_pm1(clhss)) { som = addis(som, ep-j+1); break; }
som = addii(som, clhss);
}
idealrel = flrel
? idealmul(nf,idealrel, idealpow(nf,pr, som))
: mulii(idealrel, powgi((GEN)pr[1], mulii(som,(GEN)pr[4])));
}
dlk = flrel
? idealdiv(nf,idealpow(nf,ideal,clhray),idealrel)
: divii(powgi(dethnf(ideal),clhray),idealrel);
p2[1]=(long)ideal; arch2=dummycopy(arch); p2[2]=(long)arch2; nbdezero=0;
for (k=1; k<=r1; k++)
{
if (signe(arch[k]))
{
arch2[k]=zero;
if (trivial) clhss=rayclassno(bnf,p2);
else
{
bnr2=buchrayall(bnf,p2,nf_INIT,prec);
clhss=cardofimagofquotientgroup(gen,bnr2,subgroup);
}
arch2[k]=un;
if (egalii(clhss,clhray))
{
if (flcond) { avma=av; return gzero; }
nbdezero++;
}
}
else nbdezero++;
}
tetpil=avma; p1=cgetg(4,t_VEC);
p1[1]=lcopy(clhray);
p1[2]=lstoi(nbdezero);
p1[3]=lcopy(dlk); return gerepile(av,tetpil,p1);
}
static GEN
discrayabsall(GEN bnr, GEN subgroup,long flag,long prec)
{
long av=avma,tetpil,degk,clhray,n,R1;
GEN p1,p2,p3,p4,nf,dkabs,bnf;
p1=discrayrelall(bnr,subgroup,flag,prec);
if (flag & nf_REL) return p1;
if (p1==gzero) { avma=av; return gzero; }
bnf=(GEN)bnr[1]; nf=(GEN)bnf[7]; degk=lgef(nf[1])-3;
dkabs=absi((GEN)nf[3]); p2=(GEN)p1[3];
clhray=itos((GEN)p1[1]); p3=gpuigs(dkabs,clhray);
n = degk * clhray; R1 = itos((GEN)p1[2]) * clhray;
if (((n-R1)&3)==2) p2=negi(p2);
tetpil=avma; p4=cgetg(4,t_VEC);
p4[1]=lstoi(n);
p4[2]=lstoi(R1);
p4[3]=lmulii(p2,p3); return gerepile(av,tetpil,p4);
}
GEN
bnrdisc0(GEN arg0, GEN arg1, GEN arg2, long flag, long prec)
{
GEN subgroup, bnr = args_to_bnr(arg0,arg1,arg2,&subgroup,prec);
return discrayabsall(bnr,subgroup,flag,prec);
}
GEN
discrayrel(GEN bnr, GEN subgroup,long prec)
{
return discrayrelall(bnr,subgroup,nf_REL,prec);
}
GEN
discrayrelcond(GEN bnr, GEN subgroup,long prec)
{
return discrayrelall(bnr,subgroup,nf_REL | nf_COND,prec);
}
GEN
discrayabs(GEN bnr, GEN subgroup,long prec)
{
return discrayabsall(bnr,subgroup,nf_REGULAR,prec);
}
GEN
discrayabscond(GEN bnr, GEN subgroup,long prec)
{
return discrayabsall(bnr,subgroup,nf_COND,prec);
}
/* Etant donnes un corps de nombres bnf=bnr[1], un groupe de classes de rayon
* bnr calcule par buchrayinit(), et un vecteur chi representant un caractere
* sur les generateurs bnr[2][3], calcule le conducteur de ce caractere.
*/
GEN
bnrconductorofchar(GEN bnr,GEN chi,long prec)
{
long nbgen,i,av=avma,tetpil;
GEN p1,m,u,d1,cl,cyclic;
checkbnrgen(bnr); cl=(GEN)bnr[5];
cyclic=(GEN)cl[2]; nbgen=lg(cyclic)-1;
if (lg(chi)-1 != nbgen)
err(talker,"incorrect character length in conductorofchar");
if (!nbgen) return conductor(bnr,gzero,0,prec);
d1=(GEN)cyclic[1]; m=cgetg(nbgen+2,t_MAT);
for (i=1; i<=nbgen; i++)
{
p1=cgetg(2,t_COL); m[i]=(long)p1;
p1[1]=ldiv(gmul((GEN)chi[i],d1),(GEN)cyclic[i]);
if (typ(p1[1])!=t_INT) err(typeer,"conductorofchar");
}
p1=cgetg(2,t_COL); m[i]=(long)p1;
p1[1]=(long)d1; u=(GEN)hnfall(m)[2];
tetpil=avma; setlg(u,nbgen+1);
for (i=1; i<=nbgen; i++) setlg(u[i],nbgen+1);
return gerepile(av,tetpil,conductor(bnr,u,0,prec));
}
/* Etant donne la liste des zidealstarinit et des matrices d'unites
* correspondantes, sort la liste des nombres de classes
*/
GEN
rayclassnolist(GEN bnf,GEN listes)
{
long av=avma,tetpil,i,j,k,kk,nc,nq,lx,ly;
GEN h,modulist,unitlist,classlist,sous,sousu,sousclass,p2,m,bid,q,cyclic;
if (typ(listes)!=t_VEC || lg(listes)!=3) err(typeer,"rayclassnolist");
bnf = checkbnf(bnf); h=gmael3(bnf,8,1,1);
modulist=(GEN)listes[1]; unitlist=(GEN)listes[2];
lx=lg(modulist); classlist=cgetg(lx,t_VEC);
for (i=1; i<lx; i++)
{
sous=(GEN)modulist[i]; sousu=(GEN)unitlist[i]; ly=lg(sous);
sousclass=cgetg(ly,t_VEC); classlist[i]=(long)sousclass;
for (j=1; j<ly; j++)
{
bid=(GEN)sous[j]; q=(GEN)sousu[j]; nq=lg(q)-1;
cyclic=gmael(bid,2,2); nc=lg(cyclic)-1;
if (lg(q[1]) != nc+1) err(bugparier,"rayclassnolist");
m=cgetg(nq+nc+1,t_MAT);
for (k=1; k<=nq; k++) m[k]=q[k];
for ( ; k<=nq+nc; k++)
{
p2=cgetg(nc+1,t_COL); m[k]=(long)p2;
for (kk=1; kk<=nc; kk++) p2[kk]=(kk==k-nq)?cyclic[kk]:zero;
}
sousclass[j] = lmul(h,dethnf(hnf(m)));
}
}
tetpil=avma; return gerepile(av,tetpil,gcopy(classlist));
}
static GEN
rayclassnolistes(GEN sous, GEN sousclass, GEN fac)
{
long i;
for (i=1; i<lg(sous); i++)
if (gegal(gmael(sous,i,3),fac)) return (GEN)sousclass[i];
err(bugparier,"discrayabslist");
return NULL; /* not reached */
}
static GEN
rayclassnolistessimp(GEN sous, GEN fac)
{
long i;
for (i=1; i<lg(sous); i++)
if (gegal(gmael(sous,i,1),fac)) return gmael(sous,i,2);
err(bugparier,"discrayabslistlong");
return NULL; /* not reached */
}
static GEN
factormul(GEN fa1,GEN fa2)
{
GEN pr,prnew,ex,exnew,v,p,y=cgetg(3,t_MAT);
long i,c,lx;
pr=concatsp((GEN)fa1[1],(GEN)fa2[1]); y[1]=(long)pr;
ex=concatsp((GEN)fa1[2],(GEN)fa2[2]); y[2]=(long)ex;
v=sindexsort(pr); lx=lg(pr);
prnew=cgetg(lx,t_COL); exnew=cgetg(lx,t_COL);
for (i=1; i<lx; i++) prnew[i]=pr[v[i]];
for (i=1; i<lx; i++) exnew[i]=ex[v[i]];
p=gzero; c=0;
for (i=1; i<lx; i++)
{
if (gegal((GEN)prnew[i],p))
ex[c]=laddii((GEN)ex[c],(GEN)exnew[i]);
else
{
c++; p=(GEN)prnew[i]; pr[c]=(long)p; ex[c]=exnew[i];
}
}
setlg(pr,c+1); setlg(ex,c+1); return y;
}
static GEN
factordivexact(GEN fa1,GEN fa2)
{
long i,j,k,c,lx1,lx2;
GEN listpr,listex,y,listpr1,listex1,listpr2,listex2,p1;
listpr1=(GEN)fa1[1]; listex1=(GEN)fa1[2]; lx1=lg(listpr1);
listpr2=(GEN)fa2[1]; listex2=(GEN)fa2[2]; lx2=lg(listpr1);
y=cgetg(3,t_MAT);
listpr=cgetg(lx1,t_COL); y[1]=(long)listpr;
listex=cgetg(lx1,t_COL); y[2]=(long)listex;
for (c=0,i=1; i<lx1; i++)
{
j=isinvector(listpr2,(GEN)listpr1[i],lx2-1);
if (!j) { c++; listpr[c]=listpr1[i]; listex[c]=listex1[i]; }
else
{
p1=subii((GEN)listex1[i],(GEN)listex2[j]); k=signe(p1);
if (k<0) err(talker,"factordivexact is not exact!");
if (k>0) { c++; listpr[c]=listpr1[i]; listex[c]=(long)p1; }
}
}
setlg(listpr,c+1); setlg(listex,c+1); return y;
}
static GEN
factorpow(GEN fa,long n)
{
GEN y=cgetg(3,t_MAT);
if (!n) { y[1]=lgetg(1,t_COL); y[2]=lgetg(1,t_COL); return y; }
y[1]=fa[1]; y[2]=lmulsg(n,(GEN)fa[2]); return y;
}
/* Etant donne la liste des zidealstarinit et des matrices d'unites
* correspondantes, sort la liste des discrayabs. On ne garde pas les modules
* qui ne sont pas des conducteurs
*/
GEN
discrayabslist(GEN bnf,GEN listes)
{
long av=avma,tetpil,ii,jj,i,j,k,clhss,ep,clhray,lx,ly,r1,som,degk,nbdezero;
long n,R1,normps,normi,lfa2;
GEN classlist,modulist,disclist,nf,dkabs,sous,sousclass,sousdisc;
GEN bid,module,ideal,arch,fa,fa2,ex,idealrel,p1,p2,pr;
GEN dlk,arch2,p3,fac,normp,fad,fad1,fad2,no1,no2;
classlist=rayclassnolist(bnf,listes); lx=lg(classlist);
modulist=(GEN)listes[1];
disclist=cgetg(lx,t_VEC); nf=(GEN)bnf[7]; r1=itos(gmael(nf,2,1));
degk=lgef(nf[1])-3; dkabs=gabs((GEN)nf[3],0);
for (ii=1; ii<lx; ii++)
{
sous=(GEN)modulist[ii]; sousclass=(GEN)classlist[ii];
ly=lg(sous); sousdisc=cgetg(ly,t_VEC); disclist[ii]=(long)sousdisc;
for (jj=1; jj<ly; jj++)
{
bid=(GEN)sous[jj]; clhray=itos((GEN)sousclass[jj]);
module=(GEN)bid[1]; ideal=(GEN)module[1]; arch=(GEN)module[2];
fa=(GEN)bid[3]; fa2=(GEN)fa[1]; ex=(GEN)fa[2]; fac=gcopy(fa);
idealrel=cgetg(3,t_MAT);
idealrel[1]=lgetg(1,t_COL);
idealrel[2]=lgetg(1,t_COL); lfa2=lg(fa2);
for (k=1; k<lfa2; k++)
{
pr=(GEN)fa2[k]; ep=itos((GEN)ex[k]); som=0;
normp=gpui((GEN)pr[1],(GEN)pr[4],0);
normps=itos(normp); normi=ii;
normp=cgetg(3,t_MAT);
no1=cgetg(2,t_COL); no1[1]=pr[1]; normp[1]=(long)no1;
no2=cgetg(2,t_COL); no2[1]=pr[4]; normp[2]=(long)no2;
for (j=1; j<=ep; j++)
{
if (j<ep) { coeff(fac,k,2)=lstoi(ep-j); fad=fac; }
else
{
fad=cgetg(3,t_MAT);
fad1=cgetg(lfa2-1,t_COL); fad[1]=(long)fad1;
fad2=cgetg(lfa2-1,t_COL); fad[2]=(long)fad2;
for (i=1; i<k; i++)
{ fad1[i]=coeff(fac,i,1); fad2[i]=coeff(fac,i,2); }
for ( ; i<lfa2-1; i++)
{ fad1[i]=coeff(fac,i+1,1); fad2[i]=coeff(fac,i+1,2); }
}
normi /= normps;
clhss = itos(rayclassnolistes((GEN)modulist[normi],
(GEN)classlist[normi],fad));
if (j==1 && clhss==clhray)
{
clhray=0; coeff(fac,k,2)=ex[k]; goto LLDISCRAY;
}
if (clhss==1) { som += ep-j+1; break; }
som += clhss;
}
coeff(fac,k,2)=ex[k];
idealrel=factormul(idealrel,factorpow(normp,som));
}
dlk=factordivexact(factorpow(factor(stoi(ii)),clhray),idealrel);
p2=cgetg(3,t_VEC);
p2[1]=(long)ideal; arch2=dummycopy(arch);
p2[2]=(long)arch2; nbdezero=0;
for (k=1; k<=r1; k++)
{
if (signe(arch[k]))
{
arch2[k]=zero;
clhss=itos(rayclassno(bnf,p2));
arch2[k]=un;
if (clhss==clhray) { clhray=0; break; }
}
else nbdezero++;
}
LLDISCRAY:
if (!clhray) sousdisc[jj]=lgetg(1,t_VEC);
else
{
p3=factorpow(factor(dkabs),clhray); n=clhray*degk; R1=nbdezero*clhray;
if (((n-R1)&3) == 2)
{
normp=cgetg(3,t_MAT);
no1=cgetg(2,t_COL); normp[1]=(long)no1; no1[1]=lstoi(-1);
no2=cgetg(2,t_COL); normp[2]=(long)no2; no2[1]=un;
dlk=factormul(normp,dlk);
}
p1=cgetg(4,t_VEC); p1[1]=lstoi(n);
p1[2]=lstoi(R1); p1[3]=(long)factormul(dlk,p3);
sousdisc[jj]=(long)p1;
}
}
}
tetpil=avma; return gerepile(av,tetpil,gcopy(disclist));
}
#define SHLGVINT 15
#define LGVINT 32768 /* must be 1<<SHLGVINT */
/* Attention: bound est le nombre de vraies composantes voulues. */
static GEN
bigcgetvec(long bound)
{
long nbcext,nbfinal,i;
GEN vext;
nbcext = ((bound-1)>>SHLGVINT)+1;
nbfinal = bound-((nbcext-1)<<SHLGVINT);
vext = cgetg(nbcext+1,t_VEC);
for (i=1; i<nbcext; i++) vext[i]=lgetg(LGVINT+1,t_VEC);
vext[nbcext]=lgetg(nbfinal+1,t_VEC); return vext;
}
static GEN
getcompobig(GEN vext,long i)
{
long cext;
if (i<=LGVINT) return gmael(vext,1,i);
cext=((i-1)>>SHLGVINT)+1;
return gmael(vext,cext,i-((cext-1)<<SHLGVINT));
}
static void
putcompobig(GEN vext,long i,GEN x)
{
long cext;
if (i<=LGVINT) { mael(vext,1,i)=(long)x; return; }
cext=((i-1)>>SHLGVINT)+1; mael(vext,cext,i-((cext-1)<<SHLGVINT))=(long)x;
return;
}
static GEN
zsimp(GEN bid, GEN matunit)
{
GEN y=cgetg(5,t_VEC);
y[1]=lcopy((GEN)bid[3]); y[2]=lcopy(gmael(bid,2,2));
y[3]=lcopy((GEN)bid[5]); y[4]=lcopy(matunit); return y;
}
static GEN
zsimpjoin(GEN bidsimp, GEN bidp, GEN faussefa, GEN matunit)
{
long i,j,lx,lx1,lx2,llx,llx1,llx2,nbgen,av=avma,tetpil,c;
GEN y,U1,U2,cyclic1,cyclic2,U,cyc,u1u2,p1,p2,met;
y=cgetg(5,t_VEC); y[1]=(long)vconcat((GEN)bidsimp[1],faussefa);
U1=(GEN)bidsimp[3]; U2=(GEN)bidp[5]; cyclic1=(GEN)bidsimp[2];
cyclic2=gmael(bidp,2,2); lx1=lg(U1); lx2=lg(U2); lx=lx1+lx2-1;
llx1=lg(cyclic1); llx2=lg(cyclic2);
llx=llx1+llx2-1; nbgen=llx-1; U=cgetg(lx,t_MAT);
if (nbgen)
{
cyc=diagonal(concatsp(cyclic1,cyclic2));
u1u2=smithclean(smith2(cyc)); met=(GEN)u1u2[3]; c=lg(met)-1;
for (j=1; j<lx1; j++)
{
p1=cgetg(llx,t_COL); p2=(GEN)U1[j]; U[j]=(long)p1;
for (i=1; i<llx1; i++) p1[i]=p2[i];
for ( ; i<llx; i++) p1[i]=zero;
}
for ( ; j<lx; j++)
{
p1=cgetg(llx,t_COL); p2=(GEN)U2[j-lx1+1]; U[j]=(long)p1;
for (i=1; i<llx1; i++) p1[i]=zero;
for ( ; i<llx; i++) p1[i]=p2[i-llx1+1];
}
y[3]=lmul((GEN)u1u2[1],U);
}
else
{
met=cgetg(1,t_MAT); for (j=1; j<lx; j++) U[j]=lgetg(1,t_COL);
y[3]=(long)U; c=0;
}
cyc=cgetg(c+1,t_VEC); for (i=1; i<=c; i++) cyc[i]=coeff(met,i,i);
y[2]=(long)cyc;
y[4]=(long)vconcat((GEN)bidsimp[4],matunit);
tetpil=avma; return gerepile(av,tetpil,gcopy(y));
}
static GEN
rayclassnointern(GEN sous, GEN clh)
{
long lx,nc,nq,k,kk,j;
GEN bidsimp,qm,sousray,cyclic,m,p2,p1;
lx=lg(sous); sousray=cgetg(lx,t_VEC);
for (j=1; j<lx; j++)
{
bidsimp=(GEN)sous[j]; qm=gmul((GEN)bidsimp[3],(GEN)bidsimp[4]);
nq=lg(qm)-1; cyclic=(GEN)bidsimp[2]; nc=lg(cyclic)-1;
if (lg(qm[1]) != nc+1) err(bugparier,"rayclassnolist");
m=cgetg(nq+nc+1,t_MAT);
for (k=1; k<=nq; k++) m[k]=qm[k];
for ( ; k<=nq+nc; k++)
{
p2=cgetg(nc+1,t_COL); m[k]=(long)p2;
for (kk=1; kk<=nc; kk++) p2[kk]=(kk==k-nq)?cyclic[kk]:zero;
}
p1=cgetg(3,t_VEC); p1[2]=lmul(clh,dethnf(hnf(m)));
p1[1]=bidsimp[1]; sousray[j]=(long)p1;
}
return sousray;
}
static GEN
rayclassnointernarch(GEN sous, GEN clh, GEN matarchunit)
{
long lx,nc,nq,k,kk,j,lm,lh,r1,jj,i,nba,nbarch,ii;
GEN bidsimp,qm,sousray,cyclic,m,p2,p1,p1all,p3,mm,mj,qmk,matk;
if (!matarchunit) return rayclassnointern(sous,clh);
lm=lg(matarchunit); if (!lm) err(talker,"no units in rayclassnointernarch");
r1=lg(matarchunit[1])-1; if (r1>15) err(talker,"r1>15 in rayclassnointernarch");
lx=lg(sous); sousray=cgetg(lx,t_VEC);
for (j=1; j<lx; j++)
{
bidsimp=(GEN)sous[j]; qm=gmul((GEN)bidsimp[3],(GEN)bidsimp[4]);
nq=lg(qm)-1; cyclic=(GEN)bidsimp[2]; nc=lg(cyclic)-1;
if (lm != nq+1) err(bugparier,"rayclassnointernarch (1)");
if (lg(qm[1]) != nc+1) err(bugparier,"rayclassnointernarch (2)");
m=cgetg(nq+nc+r1+1,t_MAT);
for (k=1; k<=nq; k++)
{
p2=cgetg(nc+r1+1,t_COL); m[k]=(long)p2; qmk=(GEN)qm[k];
matk=(GEN)matarchunit[k];
for (kk=1; kk<=nc; kk++) p2[kk]=qmk[kk];
for ( ; kk<=nc+r1; kk++) p2[kk]=matk[kk-nc];
}
for ( ; k<=nq+nc; k++)
{
p2=cgetg(nc+r1+1,t_COL); m[k]=(long)p2;
for (kk=1; kk<=nc+r1; kk++) p2[kk]=(kk==k-nq)?cyclic[kk]:zero;
}
for ( ; k<=nq+nc+r1; k++)
{
p2=cgetg(nc+r1+1,t_COL); m[k]=(long)p2;
for (kk=1; kk<=nc+r1; kk++) p2[kk]=(kk==k-nq)?deux:zero;
}
m=hnf(m);
nbarch=(1<<r1); p1all=cgetg(nbarch+1,t_VEC); lh=lg(m);
if (lh!=nc+r1+1) err(bugparier,"rayclassnointernarch (3)");
for (k=0; k<=nbarch-1; k++)
{
p2=cgetg(r1+1,t_COL); kk=k; nba=0;
for (jj=1; jj<=r1; jj++)
{
if (kk%2) { nba++; p2[jj]=un; } else p2[jj]=zero;
kk>>=1;
}
mm=cgetg(lh,t_MAT);
for (jj=1; jj<lh; jj++)
{
p3=cgetg(nc+nba+1,t_COL); mm[jj]=(long)p3; mj=(GEN)m[jj];
for (i=1; i<=nc; i++) p3[i]=mj[i];
for (ii=1; ii<=r1; ii++)
if (signe(p2[ii])) p3[i++]=mj[nc+ii];
}
p1all[k+1]=lmul(clh,dethnf(hnf(mm)));
}
p1=cgetg(3,t_VEC); p1[2]=(long)p1all; p1[1]=bidsimp[1];
sousray[j]=(long)p1;
}
return sousray;
}
GEN
decodemodule(GEN nf, GEN fa)
{
long n,k,j,fauxpr,av=avma;
GEN fa2,ex,id,pr;
nf=checknf(nf);
if (typ(fa)!=t_MAT || lg(fa)!=3)
err(talker,"incorrect factorisation in decodemodule");
n=lgef(nf[1])-3; id=idmat(n);
fa2=(GEN)fa[1]; ex=(GEN)fa[2];
for (k=1; k<lg(fa2); k++)
{
fauxpr=itos((GEN)fa2[k]);
j=(fauxpr%n)+1; fauxpr /= n*n;
pr = (GEN)primedec(nf,stoi(fauxpr))[j];
id = idealmul(nf,id,idealpow(nf,pr,(GEN)ex[k]));
}
return gerepileupto(av,id);
}
/* Fait tout a partir du debut, et bound peut aller jusqu'a 2^30. Pour
* l'instant ne s'occupe pas des sous-groupes. le format de sortie est le
* suivant: un vecteur vext dont les composantes vint ont exactement 2^LGVINT
* vraies composantes sauf le dernier qui peut etre plus court. vext[i][j]
* (ou j<=2^LGVINT) doit etre compris comme la composante d'indice
* I=(i-1)*2^LGVINT+j d'un grand vecteur unique V. Une telle composante
* d'indice I est un vecteur indexe par tous les ideaux de norme egal a I. Si
* m_0 est un tel ideal, la composante correspondante est la suivante:
*
* + si arch = NULL, on parcourt toutes les parties archimediennes possibles.
* L'ordre des arch est l'ordre lexicographique inverse, donc [0,0,..,0],
* [1,0,..,0], [0,1,..,0],... Dans ce cas: [m_0,V] ou V est un vecteur de
* 2^r1 composantes, donnant pour chaque arch, le triplet [N,R1,D], ou N, R1,
* D sont comme dans discrayabs sauf que D est sous forme factorisee.
*
* + sinon [m_0,arch,N,R1,D], ou N, R1, D sont comme ci-dessus.
*
* Si ramip est different de 0 et -1, ne prend que les modules sans facteur
* carres ailleurs qu'au dessus de ramip. Si ramip est strictement negatif,
* archsquare.
*/
static GEN
discrayabslistarchintern(GEN bnf, GEN arch, long bound, long ramip)
{
byteptr ptdif=diffptr;
long degk,lim,av0,av,av1,tetpil,i,j,k,p2s,lfa,lp1,sqbou,cex, allarch;
long ffs,fs,resp,flbou,tdi,nba, k2,karch,kka,nbarch,jjj,jj,square;
long ii2,ii,ly,clhray,lfa2,ep,som,clhss,normps,normi,nbdezero,r1,R1,n,lp4,c;
ulong q;
GEN nf,p,z,pol,p1,p2,p3,fa,pr,normp,ideal,bidp,z2,matarchunit;
GEN bigres,funits,racunit,embunit,sous,clh,sousray,raylist;
GEN clhrayall,discall,faall,Id,idealrel,idealrelinit;
GEN sousdisc,module,fa2,ex,fac,no1,no2,fad,fad1,fad2,fadkabs,pz;
GEN arch2,dlk,disclist,bidsimp,p4,faussefa,pex,fauxpr,gprime;
GEN *gptr[2];
/* ce qui suit recopie d'assez pres ideallistzstarall */
if (DEBUGLEVEL>2) timer2();
if (bound <= 0) err(talker,"non-positive bound in discrayabslist");
av0=avma; bnf = checkbnf(bnf); flbou=0;
nf=(GEN)bnf[7]; bigres=(GEN)bnf[8]; pol=(GEN)nf[1]; degk=lgef(pol)-3;
r1=itos(gmael(nf,2,1)); fadkabs=factor(absi((GEN)nf[3]));
clh=gmael(bigres,1,1);
funits = check_units(bnf,"discrayabslistarchintern");
racunit=gmael(bigres,4,2);
if (ramip >= 0) square = 0;
else
{
square = 1; ramip = -ramip;
if (ramip==1) ramip=0;
}
nba = 0; allarch = (arch==NULL);
if (allarch)
{ arch=cgetg(r1+1,t_VEC); for (i=1; i<=r1; i++) arch[i]=un; nba=r1; }
else if (gcmp0(arch))
{ arch=cgetg(r1+1,t_VEC); for (i=1; i<=r1; i++) arch[i]=zero; }
else
{
if (lg(arch)!=r1+1)
err(talker,"incorrect archimedean argument in discrayabslistlong");
for (i=1; i<=r1; i++) if (signe(arch[i])) nba++;
if (nba) module=cgetg(3,t_VEC);
}
p1=cgetg(3,t_VEC); p1[1]=(long)idmat(degk); p1[2]=(long)arch;
bidp=zidealstarinitall(nf,p1,0);
matarchunit = allarch? logunitmatrix(nf,funits,racunit,bidp): (GEN)NULL;
p=cgeti(3); p[1]=evalsigne(1)|evallgef(3);
sqbou=itos(racine(stoi(bound)))+1;
av=avma; lim=stack_lim(av,1);
z=bigcgetvec(bound); for (i=2; i<=bound; i++) putcompobig(z,i,cgetg(1,t_VEC));
if (allarch) bidp=zidealstarinitall(nf,idmat(degk),0);
embunit=logunitmatrix(nf,funits,racunit,bidp);
p1=cgetg(2,t_VEC); putcompobig(z,1,p1); p1[1]=(long)zsimp(bidp,embunit);
if (DEBUGLEVEL>=2) fprintferr("Starting zidealstarunits computations\n");
for (p[2]=0; p[2]<=bound; )
{
if (!*ptdif) err(primer1);
p[2] += *ptdif++;
if (!flbou && p[2]>sqbou)
{
flbou=1;
if (DEBUGLEVEL>=2)
{ fprintferr("\nStarting rayclassno computations\n"); flusherr(); }
tetpil=avma; z=gerepile(av,tetpil,gcopy(z));
av1=avma; raylist=bigcgetvec(bound);
/* maintenant on suit rayclassnolist */
for (i=1; i<=bound; i++)
{
sous=getcompobig(z,i);
sousray=rayclassnointernarch(sous,clh,matarchunit);
putcompobig(raylist,i,sousray);
}
tetpil=avma; raylist=gerepile(av1,tetpil,gcopy(raylist));
z2=bigcgetvec(sqbou);
for (i=1; i<=sqbou; i++)
{ p1=gcopy(getcompobig(z,i)); putcompobig(z2,i,p1); }
z = z2;
}
fa=primedec(nf,p); lfa=lg(fa)-1;
for (j=1; j<=lfa; j++)
{
pr=(GEN)fa[j]; normp=powgi(p,(GEN)pr[4]); cex=0;
if (DEBUGLEVEL>=2) { fprintferr("%ld ",p[2]); flusherr(); }
if (gcmpgs(normp,bound)<=0)
{
fauxpr=stoi(p[2]*degk*degk+(itos((GEN)pr[4])-1)*degk+j-1);
q=p2s=itos(normp); ideal=pr;
while (q <= (ulong)bound)
{
bidp=zidealstarinitall(nf,ideal,0);
faussefa=cgetg(3,t_MAT); p1=cgetg(2,t_COL);
faussefa[1]=(long)p1; p1[1]=(long)fauxpr;
pex=cgetg(2,t_COL); faussefa[2]=(long)pex;
cex++; pex[1]=lstoi(cex);
embunit=logunitmatrix(nf,funits,racunit,bidp);
for (i=q; i<=bound; i+=q)
{
p1=getcompobig(z,i/q); lp1=lg(p1);
if (lp1>1)
{
p2=cgetg(lp1,t_VEC); c=0;
for (k=1; k<lp1; k++)
{
p3=(GEN)p1[k]; if (i>q) { p4=gmael(p3,1,1); lp4=lg(p4)-1; }
if (i==q || !isinvector(p4,fauxpr,lp4))
{
c++;
p2[c]=(long)zsimpjoin(p3,bidp,faussefa,embunit);
}
}
setlg(p2,c+1);
if (c)
{
if (p[2]<=sqbou)
{
pz=getcompobig(z,i);
if (lg(pz)>1) putcompobig(z,i,concatsp(pz,p2));
else putcompobig(z,i,p2);
}
else
{
sousray=rayclassnointernarch(p2,clh,matarchunit);
pz=getcompobig(raylist,i);
if (lg(pz)>1) putcompobig(raylist,i,concatsp(pz,sousray));
else putcompobig(raylist,i,sousray);
}
}
}
}
if (ramip && ramip % p[2]) q=bound+1;
else
{
pz=mulss(q,p2s);
q=(gcmpgs(pz,bound)>0)?bound+1:pz[2];
if (q <= (ulong)bound) ideal=idealmul(nf,ideal,pr);
}
}
}
}
if (low_stack(lim, stack_lim(av,1)))
{
if(DEBUGMEM>1) err(warnmem,"[1]: discrayabslistarch");
if (!flbou)
{
if (DEBUGLEVEL>2)
fprintferr("avma = %ld, t(z) = %ld ",avma-bot,taille2(z));
tetpil=avma; z=gerepile(av,tetpil,gcopy(z));
}
else
{
if (DEBUGLEVEL>2)
fprintferr("avma = %ld, t(r) = %ld ",avma-bot,taille2(raylist));
gptr[0]=&z; gptr[1]=&raylist; gerepilemany(av,gptr,2);
}
if (DEBUGLEVEL>2) { fprintferr("avma = %ld ",avma-bot); flusherr(); }
}
}
if (!flbou) /* maintenant on suit rayclassnolist */
{
if (DEBUGLEVEL>=2)
{ fprintferr("\nStarting rayclassno computations\n"); flusherr(); }
tetpil=avma; z=gerepile(av,tetpil,gcopy(z));
av1=avma; raylist=bigcgetvec(bound);
for (i=1; i<=bound; i++)
{
sous=getcompobig(z,i);
sousray=rayclassnointernarch(sous,clh,matarchunit);
putcompobig(raylist,i,sousray);
}
}
if (DEBUGLEVEL>2)
fprintferr("avma = %ld, t(r) = %ld ",avma-bot,taille2(raylist));
tetpil=avma; raylist=gerepile(av,tetpil,gcopy(raylist));
if (DEBUGLEVEL>2)
{ fprintferr("avma = %ld ",avma-bot); msgtimer("zidealstarlist"); }
/* maintenant on suit discrayabslist */
if (DEBUGLEVEL>=2)
{ fprintferr("Starting discrayabs computations\n"); flusherr(); }
if (allarch) nbarch = 1<<r1;
else
{
nbarch = 1;
clhrayall = cgetg(2,t_VEC);
discall = cgetg(2,t_VEC);
faall = cgetg(2,t_VEC);
Id = idmat(degk);
}
idealrelinit=cgetg(3,t_MAT);
idealrelinit[1]=lgetg(1,t_COL);
idealrelinit[2]=lgetg(1,t_COL);
av1=avma; lim=stack_lim(av1,1);
if (square) bound = sqbou-1;
disclist=bigcgetvec(bound);
for (ii=1; ii<=bound; ii++) putcompobig(disclist,ii,cgetg(1,t_VEC));
for (ii2=1; ii2<=bound; ii2++)
{
ii = square? ii2*ii2: ii2;
if (DEBUGLEVEL>=2 && ii%100==0) { fprintferr("%ld ",ii); flusherr(); }
sous=getcompobig(raylist,ii); ly=lg(sous); sousdisc=cgetg(ly,t_VEC);
putcompobig(disclist, square? ii2: ii,sousdisc);
for (jj=1; jj<ly; jj++)
{
bidsimp=(GEN)sous[jj]; fa=(GEN)bidsimp[1]; fac=gcopy(fa);
fa2=(GEN)fa[1]; ex=(GEN)fa[2]; lfa2=lg(fa2);
if (allarch)
{
clhrayall = (GEN)bidsimp[2];
discall=cgetg(nbarch+1,t_VEC);
}
else
clhray = itos((GEN)bidsimp[2]);
for (karch=0; karch<nbarch; karch++)
{
if (!allarch) ideal = Id;
else
{
kka=karch; nba=0;
for (jjj=1; jjj<=r1; jjj++)
{
if (kka%2) nba++;
kka>>=1;
}
clhray=itos((GEN)clhrayall[karch+1]);
for (k2=1,k=1; k<=r1; k++,k2<<=1)
{
if (karch&k2 && clhray==itos((GEN)clhrayall[karch-k2+1]))
{ clhray=0; goto LDISCRAY; }
}
}
idealrel = idealrelinit;
for (k=1; k<lfa2; k++)
{
fauxpr=(GEN)fa2[k]; ep=itos((GEN)ex[k]); ffs=itos(fauxpr);
/* Hack for NeXTgcc 2.5.8 -- splitting "resp=fs%degk+1" */
fs=ffs/degk; resp=fs%degk; resp++; gprime=stoi((long)(fs/degk));
if (!allarch && nba)
{
p1=(GEN)primedec(nf,gprime)[ffs%degk+1];
ideal = idealmul(nf,ideal,idealpow(nf,p1,(GEN)ex[k]));
}
som=0; clhss=0;
normp=gpuigs(gprime,resp); normps=itos(normp); normi=ii;
normp=cgetg(3,t_MAT);
no1=cgetg(2,t_COL); no1[1]=(long)gprime; normp[1]=(long)no1;
no2=cgetg(2,t_COL); no2[1]=lstoi(resp); normp[2]=(long)no2;
for (j=1; j<=ep; j++)
{
if (clhss==1) som++;
else
{
if (j<ep) { coeff(fac,k,2)=lstoi(ep-j); fad=fac; }
else
{
fad=cgetg(3,t_MAT);
fad1=cgetg(lfa2-1,t_COL); fad[1]=(long)fad1;
fad2=cgetg(lfa2-1,t_COL); fad[2]=(long)fad2;
for (i=1; i<k; i++)
{ fad1[i]=coeff(fac,i,1); fad2[i]=coeff(fac,i,2); }
for ( ; i<lfa2-1; i++)
{ fad1[i]=coeff(fac,i+1,1); fad2[i]=coeff(fac,i+1,2); }
}
normi /= normps;
/* Hack for NeXTgcc 2.5.8 -- using dlk as temporary variable */
dlk=rayclassnolistessimp(getcompobig(raylist,normi),fad);
if (allarch) dlk = (GEN)dlk[karch+1];
clhss = itos(dlk);
if (j==1 && clhss==clhray)
{
clhray=0; coeff(fac,k,2)=ex[k]; goto LDISCRAY;
}
som += clhss;
}
}
coeff(fac,k,2)=ex[k];
idealrel=factormul(idealrel,factorpow(normp,som));
}
dlk=factordivexact(factorpow(factor(stoi(ii)),clhray),idealrel);
if (!allarch && nba)
{
module[1]=(long)ideal; arch2=gcopy(arch); module[2]=(long)arch2;
nbdezero=0;
for (k=1; k<=r1; k++)
{
if (signe(arch[k]))
{
arch2[k]=zero;
clhss=itos(rayclassno(bnf,module));
arch2[k]=un;
if (clhss==clhray) { clhray=0; goto LDISCRAY; }
}
else nbdezero++;
}
}
else nbdezero = r1-nba;
LDISCRAY:
p1=cgetg(4,t_VEC); discall[karch+1]=(long)p1;
if (!clhray) p1[1]=p1[2]=p1[3]=zero;
else
{
p3=factorpow(fadkabs,clhray); n=clhray*degk; R1=nbdezero*clhray;
if (((n-R1)&3)==2)
{
normp=cgetg(3,t_MAT);
no1=cgetg(2,t_COL); normp[1]=(long)no1; no1[1]=lstoi(-1);
no2=cgetg(2,t_COL); normp[2]=(long)no2; no2[1]=un;
dlk=factormul(normp,dlk);
}
p1[1]=lstoi(n); p1[2]=lstoi(R1); p1[3]=(long)factormul(dlk,p3);
}
}
if (allarch)
{ p1 = cgetg(3,t_VEC); p1[1]=(long)fa; p1[2]=(long)discall; }
else
{ faall[1]=(long)fa; p1 = concatsp(faall, p1); }
sousdisc[jj]=(long)p1;
if (low_stack(lim, stack_lim(av1,1)))
{
if(DEBUGMEM>1) err(warnmem,"[2]: discrayabslistarch");
if (DEBUGLEVEL>2)
fprintferr("avma = %ld, t(d) = %ld ",avma-bot,taille2(disclist));
tetpil=avma; disclist=gerepile(av1,tetpil,gcopy(disclist));
if (DEBUGLEVEL>2) { fprintferr("avma = %ld ",avma-bot); flusherr(); }
}
}
}
if (DEBUGLEVEL>2) msgtimer("discrayabs");
tdi=taille2(disclist);
if (DEBUGLEVEL>2)
{ fprintferr("avma = %ld, t(d) = %ld ",avma-bot,tdi); flusherr(); }
if (tdi<avma-bot)
{
tetpil=avma; disclist=gerepile(av0,tetpil,gcopy(disclist));
if (DEBUGLEVEL>2) { fprintferr("avma = %ld ",avma-bot); flusherr(); }
}
return disclist;
}
GEN
discrayabslistarch(GEN bnf, GEN arch, long bound)
{ return discrayabslistarchintern(bnf,arch,bound, 0); }
GEN
discrayabslistlong(GEN bnf, long bound)
{ return discrayabslistarchintern(bnf,gzero,bound, 0); }
GEN
discrayabslistarchsquare(GEN bnf, GEN arch, long bound)
{ return discrayabslistarchintern(bnf,arch,bound, -1); }
static GEN
computehuv(GEN bnr,GEN id, GEN arch,long prec)
{
long i,nbgen, av = avma;
GEN bnf,bnrnew,listgen,P,U,DC;
GEN newmod=cgetg(3,t_VEC);
newmod[1]=(long)id;
newmod[2]=(long)arch;
bnf=(GEN)bnr[1];
bnrnew=buchrayall(bnf,newmod,nf_INIT,prec);
listgen=gmael(bnr,5,3); nbgen=lg(listgen);
P=cgetg(nbgen,t_MAT);
for (i=1; i<nbgen; i++)
P[i] = (long)isprincipalray(bnrnew,(GEN)listgen[i]);
DC=diagonal(gmael(bnrnew,5,2));
U=(GEN)hnfall(concatsp(P,DC))[2];
setlg(U,nbgen); for (i=1; i<nbgen; i++) setlg(U[i], nbgen);
return gerepileupto(av, hnf(concatsp(U,diagonal(gmael(bnr,5,2)))));
}
/* 0 if hinv*list[i] has a denominator for all i, 1 otherwise */
static int
hnflistdivise(GEN list,GEN h)
{
long av = avma, i, I = lg(list);
GEN hinv = ginv(h);
for (i=1; i<I; i++)
if (gcmp1(denom(gmul(hinv,(GEN)list[i])))) break;
avma = av; return i < I;
}
static GEN
subgroupcond(GEN bnr, long indexbound, long prec)
{
long av=avma,tetpil,i,j,lgi,lp;
GEN li,p1,lidet,perm,nf,bid,ideal,arch,primelist,listkernels;
checkbnrgen(bnr); bid=(GEN)bnr[2];
ideal=gmael(bid,1,1);
arch =gmael(bid,1,2); nf=gmael(bnr,1,7);
primelist=gmael(bid,3,1); lp=lg(primelist)-1;
listkernels=cgetg(lp+lg(arch),t_VEC);
for (i=1; i<=lp; )
{
p1=idealdiv(nf,ideal,(GEN)primelist[i]);
listkernels[i++]=(long)computehuv(bnr,p1,arch,prec);
}
for (j=1; j<lg(arch); j++)
if (signe((GEN)arch[j]))
{
p1=dummycopy(arch); p1[j]=zero;
listkernels[i++]=(long)computehuv(bnr,ideal,p1,prec);
}
setlg(listkernels,i);
li=subgrouplist(gmael(bnr,5,2),indexbound);
lgi=lg(li);
for (i=1,j=1; j<lgi; j++)
if (!hnflistdivise(listkernels, (GEN)li[j])) li[i++] = li[j];
/* sort by increasing index */
lgi = i; setlg(li,i); lidet=cgetg(lgi,t_VEC);
for (i=1; i<lgi; i++) lidet[i]=(long)dethnf_i((GEN)li[i]);
perm = sindexsort(lidet); p1=li; li=cgetg(lgi,t_VEC);
for (i=1; i<lgi; i++) li[i] = p1[perm[lgi-i]];
tetpil=avma; return gerepile(av,tetpil,gcopy(li));
}
GEN
subgrouplist0(GEN bnr, long indexbound, long all, long prec)
{
if (typ(bnr)!=t_VEC) err(typeer,"subgrouplist");
if (lg(bnr)!=1 && typ(bnr[1])!=t_INT)
{
if (!all) return subgroupcond(bnr,indexbound,prec);
checkbnr(bnr); bnr = gmael(bnr,5,2);
}
return subgrouplist(bnr,indexbound);
}
GEN
bnrdisclist0(GEN bnf, GEN borne, GEN arch, long all)
{
if (typ(borne)==t_INT)
{
if (!arch) arch = gzero;
if (all==1) { arch = NULL; all = 0; }
return discrayabslistarchintern(bnf,arch,itos(borne),all);
}
return discrayabslist(bnf,borne);
}
![[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