version 1.1, 2001/10/02 11:16:59 |
version 1.2, 2002/09/11 07:26:45 |
|
|
rgcd(a,b)= |
rgcd(a,b)= |
{ |
{ local(r); |
local(r); |
|
|
|
a=abs(a); b=abs(b); |
a = abs(a); |
while(b > 0.01, |
b = abs(b); |
r = a - b*truncate(a/b); |
while (b > 0.01, r=a%b; a=b; b=r); |
a=b; b=r |
|
); |
|
a |
a |
} |
} |
|
|
f(a,b)= |
f(a,b)= |
{ |
{ local(j,n,l,mv,u,vreg,cp,p); |
local(j,s,n,l,mv,vreg,cp); |
|
|
|
s=a+b*t; n=abs(norm(s)); |
n = abs( norm(a + b*t) ); |
mv=vectorv(li); |
mv = vectorv(li); |
forprime(k=2,plim, |
forprime (p=2, plim, |
l=0; while (!(n%k), l++; n/=k); |
if ((l = valuation(n,p)), |
if (l, |
n /= p^l; |
j=ind[k];cp=v[j][2]; |
j = ind[p]; cp = v[j][2]; |
while((a+b*cp)%k, |
while((a+b*cp)%p, |
j++; cp=v[j][2] |
j++; cp = v[j][2] |
); |
); |
mv[j]=l |
mv[j]=l |
) |
) |
); |
); |
if (n!=1, return); |
if (n!=1, return); |
|
|
vreg=vectorv(lireg,j, |
/* found a relation */ |
if (j<=r1, |
vreg = vectorv(lireg,j, |
log(abs(a+b*re[j])) |
u = a+b*re[j]; |
, |
if (j<=r1, abs(u), norm(u)) |
log(norm(a+b*re[j])) |
|
) |
|
); |
); |
mreg=concat(mreg,vreg); m=concat(m,mv); |
mreg = concat(mreg, log(vreg)); |
areg=concat(areg,a+b*t); |
m = concat(m,mv); |
|
areg = concat(areg, a+b*t); |
print1("(" res++ ": " a "," b ")") |
print1("(" res++ ": " a "," b ")") |
} |
} |
|
|
clareg(p, plim=19, lima=50, extra=5)= |
global(km,clh,R,nf,areg); |
{ |
|
vi=nfinit(p); p=vi.pol; |
clareg(pol, plim=19, lima=50, extra=5)= |
t=Mod(x,p,1); |
{ local(e,t,r,coreg,lireg,r1,ind,fa,li,co,res,a,b,m,mh,ms,mhs,mreg,mregh); |
r=vi.roots; r1=vi.sign[1]; |
|
findex=vi[4]; /* index: power basis <==> findex = 1 */ |
nf=nfinit(pol); pol=nf.pol; |
if (findex>1, |
t=Mod(x,pol,1); |
|
r=nf.roots; r1=nf.sign[1]; |
|
if (nf[4] > 1, /* index: power basis <==> nf[4] = 1 */ |
error("sorry, the case f>1 is not implemented") |
error("sorry, the case f>1 is not implemented") |
); |
); |
|
|
print("discriminant = " vi.disc ", signature = " vi.sign); |
print("discriminant = " nf.disc ", signature = " nf.sign); |
dep=length(p)-1; |
|
lireg=(dep+r1)/2; |
lireg = (poldegree(pol) + r1) / 2; /* r1 + r2 */ |
re=vector(lireg,j, |
re=vector(lireg,j, |
if (j<=r1, real(r[j]) , r[j]) |
if (j<=r1, real(r[j]) , r[j]) |
); |
); |
ind=vector(plim); v=[]; |
ind=vector(plim); v=[]; |
forprime(k=2,plim, |
forprime(p=2,plim, |
w=lift(factormod(p,k)); |
w = factormod(pol,p); |
|
e = w[,2]; |
find=0; |
find=0; |
for(l=1,length(w[,1]), |
for(l=1,#e, |
fa=w[l,1]; |
fa = lift(w[l,1]); |
if (length(fa)==2, |
if (poldegree(fa) == 1, |
if (!find, |
if (!find, |
find=1; ind[k]=length(v)+1 |
find=1; ind[p]=#v+1 |
); |
); |
v=concat(v,[[k,-polcoeff(fa,0),w[l,2]]]) |
v = concat(v, [[p,-polcoeff(fa,0),e[l]]]) |
) |
) |
) |
) |
); |
); |
li=length(v); co=li+extra; |
li=#v; co=li+extra; |
res=0; print("need " co); |
res=0; print("need " co " relations"); |
areg=[]~; mreg = m = [;]; |
areg=[]~; mreg = m = [;]; |
a=1; b=1; f(0,1); |
a=1; b=1; f(0,1); |
while (res<co, |
while (res<co, |
Line 88 clareg(p, plim=19, lima=50, extra=5)= |
|
Line 86 clareg(p, plim=19, lima=50, extra=5)= |
|
return |
return |
); |
); |
|
|
mhs=matsnf(mh); clh=mhs[1]; mh1=[clh]; |
mhs = matsnf(mh,4); |
j=1; |
clh = prod(i=1,#mhs, mhs[i]); |
while (j<length(mhs), |
print("class number = " clh ", class group = " 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; |
areg=Mat(areg); km=matkerint(m); mregh=mreg*km; |
if (lireg==1, |
if (lireg==1, |
a1=1 |
R = 1 |
, |
, |
coreg=length(mregh); |
coreg = #mregh; |
if (coreg<lireg-1, |
if (coreg < lireg-1, |
print("not enough relations for regulator: matsize = " matsize(mregh)); |
print("not enough relations for regulator: matsize = " matsize(mregh)); |
a1 = "(not given)"; |
R = "(not given)"; |
, |
, |
mreg1=vecextract(mregh, Str(1 ".." lireg-1), ".."); |
mreg1 = vecextract(mregh, Str(".." lireg-1), ".."); |
a1 = 0; |
R = 0; |
for(j=lireg-1,coreg, |
for(j=lireg-1,coreg, |
a = matdet(vecextract(mreg1, Str(j-lireg+2 ".." j))); |
a = matdet(vecextract(mreg1, Str(j-lireg+2 ".." j))); |
a1 = if (a1, a, rgcd(a1,a)) |
R = rgcd(a,R) |
); |
) |
) |
) |
); |
); |
print("regulator = " a1) |
print("regulator = " R) |
} |
} |
|
|
check(lim=200) = |
check(lim=200) = |
{ |
{ local(r1,r2,pol,z,Res,fa); |
r1=vi.sign[1]; r2=vi.sign[2]; pol=vi.pol; |
|
z = 2^(-r1) * (2*Pi)^(-r2) * sqrt(abs(vi.disc)); |
r1=nf.sign[1]; |
|
r2=nf.sign[2]; pol=nf.pol; |
|
z = 2^r1 * (2*Pi)^r2 / sqrt(abs(nf.disc)) / nfrootsof1(nf)[1]; |
|
Res = 1.; \\ Res (Zeta_K,s=1) ~ z * h * R |
forprime (q=2,lim, |
forprime (q=2,lim, |
z *= (q-1)/q; fa=factormod(pol,q,1)[,1]; |
fa = factormod(pol,q,1)[,1]; |
z /= prod(i=1, length(fa), 1 - 1/q^fa[i]) |
Res *= (q-1)/q / prod(i=1, #fa, 1 - q^(-fa[i])) |
); |
); |
clh*a1/nfrootsof1(vi)[1]/z |
z * clh * R / Res |
} |
} |
|
|
fu() = vector(length(km), k, factorback(concat(areg, km[,k]))) |
fu() = vector(#km, k, factorback(concat(areg, km[,k]))) |