rgcd(a,b)=
{
local(r);
a=abs(a); b=abs(b);
while(b > 0.01,
r = a - b*truncate(a/b);
a=b; b=r
);
a
}
f(a,b)=
{
local(j,s,n,l,mv,vreg,cp);
s=a+b*t; n=abs(norm(s));
mv=vectorv(li);
forprime(k=2,plim,
l=0; while (!(n%k), l++; n/=k);
if (l,
j=ind[k];cp=v[j][2];
while((a+b*cp)%k,
j++; cp=v[j][2]
);
mv[j]=l
)
);
if (n!=1, return);
vreg=vectorv(lireg,j,
if (j<=r1,
log(abs(a+b*re[j]))
,
log(norm(a+b*re[j]))
)
);
mreg=concat(mreg,vreg); m=concat(m,mv);
areg=concat(areg,a+b*t);
print1("(" res++ ": " a "," b ")")
}
clareg(p, plim=19, lima=50, extra=5)=
{
vi=nfinit(p); p=vi.pol;
t=Mod(x,p,1);
r=vi.roots; r1=vi.sign[1];
findex=vi[4]; /* index: power basis <==> findex = 1 */
if (findex>1,
error("sorry, the case f>1 is not implemented")
);
print("discriminant = " vi.disc ", signature = " vi.sign);
dep=length(p)-1;
lireg=(dep+r1)/2;
re=vector(lireg,j,
if (j<=r1, real(r[j]) , r[j])
);
ind=vector(plim); v=[];
forprime(k=2,plim,
w=lift(factormod(p,k));
find=0;
for(l=1,length(w[,1]),
fa=w[l,1];
if (length(fa)==2,
if (!find,
find=1; ind[k]=length(v)+1
);
v=concat(v,[[k,-polcoeff(fa,0),w[l,2]]])
)
)
);
li=length(v); co=li+extra;
res=0; print("need " co);
areg=[]~; mreg = m = [;];
a=1; b=1; f(0,1);
while (res<co,
if (gcd(a,b)==1,
f(a,b); f(-a,b)
);
a++;
if (a*b>lima, b++; a=1)
);
print(" ");
mh=mathnf(m); ms=matsize(mh);
if (ms[1]!=ms[2],
print("not enough relations for class group: matrix size = ",ms);
return
);
mhs=matsnf(mh); clh=mhs[1]; mh1=[clh];
j=1;
while (j<length(mhs),
j++;
if (mhs[j]>1,
mh1=concat(mh1,mhs[j]);
clh *= mhs[j]
)
);
print("class number = " clh ", class group = " mh1);
areg=Mat(areg); km=matkerint(m); mregh=mreg*km;
if (lireg==1,
a1=1
,
coreg=length(mregh);
if (coreg<lireg-1,
print("not enough relations for regulator: matsize = " matsize(mregh));
a1 = "(not given)";
,
mreg1=vecextract(mregh, Str(1 ".." lireg-1), "..");
a1 = 0;
for(j=lireg-1,coreg,
a = matdet(vecextract(mreg1, Str(j-lireg+2 ".." j)));
a1 = if (a1, a, rgcd(a1,a))
);
)
);
print("regulator = " a1)
}
check(lim=200) =
{
r1=vi.sign[1]; r2=vi.sign[2]; pol=vi.pol;
z = 2^(-r1) * (2*Pi)^(-r2) * sqrt(abs(vi.disc));
forprime (q=2,lim,
z *= (q-1)/q; fa=factormod(pol,q,1)[,1];
z /= prod(i=1, length(fa), 1 - 1/q^fa[i])
);
clh*a1/nfrootsof1(vi)[1]/z
}
fu() = vector(length(km), k, factorback(concat(areg, km[,k])))