Annotation of OpenXM_contrib/pari-2.2/examples/cl.gp, Revision 1.2
1.1 noro 1: rgcd(a,b)=
1.2 ! noro 2: { local(r);
1.1 noro 3:
1.2 ! noro 4: a = abs(a);
! 5: b = abs(b);
! 6: while (b > 0.01, r=a%b; a=b; b=r);
1.1 noro 7: a
8: }
9:
10: f(a,b)=
1.2 ! noro 11: { local(j,n,l,mv,u,vreg,cp,p);
1.1 noro 12:
1.2 ! noro 13: n = abs( norm(a + b*t) );
! 14: mv = vectorv(li);
! 15: forprime (p=2, plim,
! 16: if ((l = valuation(n,p)),
! 17: n /= p^l;
! 18: j = ind[p]; cp = v[j][2];
! 19: while((a+b*cp)%p,
! 20: j++; cp = v[j][2]
1.1 noro 21: );
22: mv[j]=l
23: )
24: );
25: if (n!=1, return);
26:
1.2 ! noro 27: /* found a relation */
! 28: vreg = vectorv(lireg,j,
! 29: u = a+b*re[j];
! 30: if (j<=r1, abs(u), norm(u))
! 31: );
! 32: mreg = concat(mreg, log(vreg));
! 33: m = concat(m,mv);
! 34: areg = concat(areg, a+b*t);
1.1 noro 35: print1("(" res++ ": " a "," b ")")
36: }
37:
1.2 ! noro 38: global(km,clh,R,nf,areg);
! 39:
! 40: clareg(pol, plim=19, lima=50, extra=5)=
! 41: { local(e,t,r,coreg,lireg,r1,ind,fa,li,co,res,a,b,m,mh,ms,mhs,mreg,mregh);
! 42:
! 43: nf=nfinit(pol); pol=nf.pol;
! 44: t=Mod(x,pol,1);
! 45: r=nf.roots; r1=nf.sign[1];
! 46: if (nf[4] > 1, /* index: power basis <==> nf[4] = 1 */
1.1 noro 47: error("sorry, the case f>1 is not implemented")
48: );
49:
1.2 ! noro 50: print("discriminant = " nf.disc ", signature = " nf.sign);
! 51:
! 52: lireg = (poldegree(pol) + r1) / 2; /* r1 + r2 */
1.1 noro 53: re=vector(lireg,j,
54: if (j<=r1, real(r[j]) , r[j])
55: );
56: ind=vector(plim); v=[];
1.2 ! noro 57: forprime(p=2,plim,
! 58: w = factormod(pol,p);
! 59: e = w[,2];
1.1 noro 60: find=0;
1.2 ! noro 61: for(l=1,#e,
! 62: fa = lift(w[l,1]);
! 63: if (poldegree(fa) == 1,
1.1 noro 64: if (!find,
1.2 ! noro 65: find=1; ind[p]=#v+1
1.1 noro 66: );
1.2 ! noro 67: v = concat(v, [[p,-polcoeff(fa,0),e[l]]])
1.1 noro 68: )
69: )
70: );
1.2 ! noro 71: li=#v; co=li+extra;
! 72: res=0; print("need " co " relations");
1.1 noro 73: areg=[]~; mreg = m = [;];
74: a=1; b=1; f(0,1);
75: while (res<co,
76: if (gcd(a,b)==1,
77: f(a,b); f(-a,b)
78: );
79: a++;
80: if (a*b>lima, b++; a=1)
81: );
82: print(" ");
83: mh=mathnf(m); ms=matsize(mh);
84: if (ms[1]!=ms[2],
85: print("not enough relations for class group: matrix size = ",ms);
86: return
87: );
88:
1.2 ! noro 89: mhs = matsnf(mh,4);
! 90: clh = prod(i=1,#mhs, mhs[i]);
! 91: print("class number = " clh ", class group = " mhs);
1.1 noro 92: areg=Mat(areg); km=matkerint(m); mregh=mreg*km;
93: if (lireg==1,
1.2 ! noro 94: R = 1
1.1 noro 95: ,
1.2 ! noro 96: coreg = #mregh;
! 97: if (coreg < lireg-1,
1.1 noro 98: print("not enough relations for regulator: matsize = " matsize(mregh));
1.2 ! noro 99: R = "(not given)";
1.1 noro 100: ,
1.2 ! noro 101: mreg1 = vecextract(mregh, Str(".." lireg-1), "..");
! 102: R = 0;
1.1 noro 103: for(j=lireg-1,coreg,
104: a = matdet(vecextract(mreg1, Str(j-lireg+2 ".." j)));
1.2 ! noro 105: R = rgcd(a,R)
! 106: )
1.1 noro 107: )
108: );
1.2 ! noro 109: print("regulator = " R)
1.1 noro 110: }
111:
112: check(lim=200) =
1.2 ! noro 113: { local(r1,r2,pol,z,Res,fa);
! 114:
! 115: r1=nf.sign[1];
! 116: r2=nf.sign[2]; pol=nf.pol;
! 117: z = 2^r1 * (2*Pi)^r2 / sqrt(abs(nf.disc)) / nfrootsof1(nf)[1];
! 118: Res = 1.; \\ Res (Zeta_K,s=1) ~ z * h * R
1.1 noro 119: forprime (q=2,lim,
1.2 ! noro 120: fa = factormod(pol,q,1)[,1];
! 121: Res *= (q-1)/q / prod(i=1, #fa, 1 - q^(-fa[i]))
1.1 noro 122: );
1.2 ! noro 123: z * clh * R / Res
1.1 noro 124: }
125:
1.2 ! noro 126: fu() = vector(#km, k, factorback(concat(areg, km[,k])))
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>