Return to cl.gp CVS log

Up to [local] / OpenXM_contrib / pari-2.2 / examples

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