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

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