Return to cl.gp CVS log

Up to [local] / OpenXM_contrib / pari / examples

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])))