File: [local] / OpenXM_contrib / pari-2.2 / examples / Attic / cl.gp (download)
Revision 1.2, Wed Sep 11 07:26:45 2002 UTC (21 years, 11 months ago) by noro
Branch: MAIN
CVS Tags: RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2 Changes since 1.1: +65 -71
lines
Upgraded pari-2.2 to pari-2.2.4.
|
rgcd(a,b)=
{ local(r);
a = abs(a);
b = abs(b);
while (b > 0.01, r=a%b; a=b; b=r);
a
}
f(a,b)=
{ local(j,n,l,mv,u,vreg,cp,p);
n = abs( norm(a + b*t) );
mv = vectorv(li);
forprime (p=2, plim,
if ((l = valuation(n,p)),
n /= p^l;
j = ind[p]; cp = v[j][2];
while((a+b*cp)%p,
j++; cp = v[j][2]
);
mv[j]=l
)
);
if (n!=1, return);
/* found a relation */
vreg = vectorv(lireg,j,
u = a+b*re[j];
if (j<=r1, abs(u), norm(u))
);
mreg = concat(mreg, log(vreg));
m = concat(m,mv);
areg = concat(areg, a+b*t);
print1("(" res++ ": " a "," b ")")
}
global(km,clh,R,nf,areg);
clareg(pol, plim=19, lima=50, extra=5)=
{ local(e,t,r,coreg,lireg,r1,ind,fa,li,co,res,a,b,m,mh,ms,mhs,mreg,mregh);
nf=nfinit(pol); pol=nf.pol;
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")
);
print("discriminant = " nf.disc ", signature = " nf.sign);
lireg = (poldegree(pol) + r1) / 2; /* r1 + r2 */
re=vector(lireg,j,
if (j<=r1, real(r[j]) , r[j])
);
ind=vector(plim); v=[];
forprime(p=2,plim,
w = factormod(pol,p);
e = w[,2];
find=0;
for(l=1,#e,
fa = lift(w[l,1]);
if (poldegree(fa) == 1,
if (!find,
find=1; ind[p]=#v+1
);
v = concat(v, [[p,-polcoeff(fa,0),e[l]]])
)
)
);
li=#v; co=li+extra;
res=0; print("need " co " relations");
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,4);
clh = prod(i=1,#mhs, mhs[i]);
print("class number = " clh ", class group = " mhs);
areg=Mat(areg); km=matkerint(m); mregh=mreg*km;
if (lireg==1,
R = 1
,
coreg = #mregh;
if (coreg < lireg-1,
print("not enough relations for regulator: matsize = " matsize(mregh));
R = "(not given)";
,
mreg1 = vecextract(mregh, Str(".." lireg-1), "..");
R = 0;
for(j=lireg-1,coreg,
a = matdet(vecextract(mreg1, Str(j-lireg+2 ".." j)));
R = rgcd(a,R)
)
)
);
print("regulator = " R)
}
check(lim=200) =
{ local(r1,r2,pol,z,Res,fa);
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,
fa = factormod(pol,q,1)[,1];
Res *= (q-1)/q / prod(i=1, #fa, 1 - q^(-fa[i]))
);
z * clh * R / Res
}
fu() = vector(#km, k, factorback(concat(areg, km[,k])))