[BACK]Return to cl.gp CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / pari-2.2 / examples

Diff for /OpenXM_contrib/pari-2.2/examples/Attic/cl.gp between version 1.1 and 1.2

version 1.1, 2001/10/02 11:16:59 version 1.2, 2002/09/11 07:26:45
Line 1 
Line 1 
 rgcd(a,b)=  rgcd(a,b)=
 {  { local(r);
   local(r);  
   
   a=abs(a); b=abs(b);    a = abs(a);
   while(b > 0.01,    b = abs(b);
     r = a - b*truncate(a/b);    while (b > 0.01, r=a%b; a=b; b=r);
     a=b; b=r  
   );  
   a    a
 }  }
   
 f(a,b)=  f(a,b)=
 {  { local(j,n,l,mv,u,vreg,cp,p);
   local(j,s,n,l,mv,vreg,cp);  
   
   s=a+b*t; n=abs(norm(s));    n = abs( norm(a + b*t) );
   mv=vectorv(li);    mv = vectorv(li);
   forprime(k=2,plim,    forprime (p=2, plim,
     l=0; while (!(n%k), l++; n/=k);      if ((l = valuation(n,p)),
     if (l,        n /= p^l;
       j=ind[k];cp=v[j][2];        j = ind[p]; cp = v[j][2];
       while((a+b*cp)%k,        while((a+b*cp)%p,
         j++; cp=v[j][2]          j++; cp = v[j][2]
       );        );
       mv[j]=l        mv[j]=l
     )      )
   );    );
   if (n!=1, return);    if (n!=1, return);
   
   vreg=vectorv(lireg,j,    /* found a relation */
     if (j<=r1,    vreg = vectorv(lireg,j,
       log(abs(a+b*re[j]))      u = a+b*re[j];
     ,      if (j<=r1, abs(u), norm(u))
       log(norm(a+b*re[j]))  
     )  
   );    );
   mreg=concat(mreg,vreg); m=concat(m,mv);    mreg = concat(mreg, log(vreg));
   areg=concat(areg,a+b*t);    m = concat(m,mv);
     areg = concat(areg, a+b*t);
   print1("(" res++ ": " a "," b ")")    print1("(" res++ ": " a "," b ")")
 }  }
   
 clareg(p, plim=19, lima=50, extra=5)=  global(km,clh,R,nf,areg);
 {  
   vi=nfinit(p); p=vi.pol;  clareg(pol, plim=19, lima=50, extra=5)=
   t=Mod(x,p,1);  { local(e,t,r,coreg,lireg,r1,ind,fa,li,co,res,a,b,m,mh,ms,mhs,mreg,mregh);
   r=vi.roots; r1=vi.sign[1];  
   findex=vi[4]; /* index: power basis <==> findex = 1 */    nf=nfinit(pol); pol=nf.pol;
   if (findex>1,    t=Mod(x,pol,1);
     r=nf.roots; r1=nf.sign[1];
     if (nf[4] > 1, /* index: power basis <==> nf[4] = 1 */
     error("sorry, the case f>1 is not implemented")      error("sorry, the case f>1 is not implemented")
   );    );
   
   print("discriminant = " vi.disc ", signature = " vi.sign);    print("discriminant = " nf.disc ", signature = " nf.sign);
   dep=length(p)-1;  
   lireg=(dep+r1)/2;    lireg = (poldegree(pol) + r1) / 2; /* r1 + r2 */
   re=vector(lireg,j,    re=vector(lireg,j,
     if (j<=r1, real(r[j]) , r[j])      if (j<=r1, real(r[j]) , r[j])
   );    );
   ind=vector(plim); v=[];    ind=vector(plim); v=[];
   forprime(k=2,plim,    forprime(p=2,plim,
     w=lift(factormod(p,k));      w = factormod(pol,p);
       e = w[,2];
     find=0;      find=0;
     for(l=1,length(w[,1]),      for(l=1,#e,
       fa=w[l,1];        fa = lift(w[l,1]);
       if (length(fa)==2,        if (poldegree(fa) == 1,
         if (!find,          if (!find,
           find=1; ind[k]=length(v)+1            find=1; ind[p]=#v+1
         );          );
         v=concat(v,[[k,-polcoeff(fa,0),w[l,2]]])          v = concat(v, [[p,-polcoeff(fa,0),e[l]]])
       )        )
     )      )
   );    );
   li=length(v); co=li+extra;    li=#v; co=li+extra;
   res=0; print("need " co);    res=0; print("need " co " relations");
   areg=[]~; mreg = m = [;];    areg=[]~; mreg = m = [;];
   a=1; b=1; f(0,1);    a=1; b=1; f(0,1);
   while (res<co,    while (res<co,
Line 88  clareg(p, plim=19, lima=50, extra=5)=
Line 86  clareg(p, plim=19, lima=50, extra=5)=
     return      return
   );    );
   
   mhs=matsnf(mh); clh=mhs[1]; mh1=[clh];    mhs = matsnf(mh,4);
   j=1;    clh = prod(i=1,#mhs, mhs[i]);
   while (j<length(mhs),    print("class number = " clh ", class group = " mhs);
     j++;  
     if (mhs[j]>1,  
       mh1=concat(mh1,mhs[j]);  
       clh *= mhs[j]  
     )  
   );  
   print("class number = " clh ", class group = " mh1);  
   areg=Mat(areg); km=matkerint(m); mregh=mreg*km;    areg=Mat(areg); km=matkerint(m); mregh=mreg*km;
   if (lireg==1,    if (lireg==1,
     a1=1      R = 1
   ,    ,
     coreg=length(mregh);      coreg = #mregh;
     if (coreg<lireg-1,      if (coreg < lireg-1,
       print("not enough relations for regulator: matsize = " matsize(mregh));        print("not enough relations for regulator: matsize = " matsize(mregh));
       a1 = "(not given)";        R = "(not given)";
     ,      ,
       mreg1=vecextract(mregh, Str(1 ".." lireg-1), "..");        mreg1 = vecextract(mregh, Str(".." lireg-1), "..");
       a1 = 0;        R = 0;
       for(j=lireg-1,coreg,        for(j=lireg-1,coreg,
         a = matdet(vecextract(mreg1, Str(j-lireg+2 ".." j)));          a = matdet(vecextract(mreg1, Str(j-lireg+2 ".." j)));
         a1 = if (a1, a, rgcd(a1,a))          R = rgcd(a,R)
       );        )
     )      )
   );    );
   print("regulator = " a1)    print("regulator = " R)
 }  }
   
 check(lim=200) =  check(lim=200) =
 {  { local(r1,r2,pol,z,Res,fa);
   r1=vi.sign[1]; r2=vi.sign[2]; pol=vi.pol;  
   z = 2^(-r1) * (2*Pi)^(-r2) * sqrt(abs(vi.disc));    r1=nf.sign[1];
     r2=nf.sign[2]; pol=nf.pol;
     z = 2^r1 * (2*Pi)^r2 / sqrt(abs(nf.disc)) / nfrootsof1(nf)[1];
     Res = 1.; \\ Res (Zeta_K,s=1) ~ z * h * R
   forprime (q=2,lim,    forprime (q=2,lim,
     z *= (q-1)/q; fa=factormod(pol,q,1)[,1];      fa = factormod(pol,q,1)[,1];
     z /= prod(i=1, length(fa), 1 - 1/q^fa[i])      Res *= (q-1)/q / prod(i=1, #fa, 1 - q^(-fa[i]))
   );    );
   clh*a1/nfrootsof1(vi)[1]/z    z * clh * R / Res
 }  }
   
 fu() = vector(length(km), k, factorback(concat(areg, km[,k])))  fu() = vector(#km, k, factorback(concat(areg, km[,k])))
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>