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