\\ -*-gp-script-*-	


tensor(a,b)=vector(length(a)*length(b),j,a[(j-1)%length(a)+1]*b[(j-1)\length(a)+1]);

rat(q)=qq=lindep([q,-1],-1);qq[2]/qq[1];

vec_to_mat(v)=matrix(length(v),length(v[1]),i,j,v[i][j]);
mat_to_vec(m)=vector(matsize(m)[1],i,m[i,]);

{
 pic(bnf)=
\\   bnf=bnfinit(x^2-j);
   rc=vector(bnf.sign[1],j,1);
   res=[bnf.clgp[1],length(idealfactor(bnf,2)~),idealfactor(bnf,2)[,2]~,bnrinit(bnf,[1,rc])[5][1]]
}


{
 torsionpart(u,v,lendec,lendecpq)=
   ld=lendec+1;
   concat(
     concat(0, \\u[1]*v[1]*ww,
        vector(ld-1,j,u[1]*v[j+1]-u[j+1]*v[1]+ww*v[j+1]*u[j+1]))%(2*ww),
     tensor(vector(lendec,j,u[j+1]),v)%(2*ww))     
}

{
  symwedge(u,v,t,s,level)=
\\ associates to the vectors u (length t) and v (length s) the 
\\ corresponding vector  u\odot{m-2}\tensor u\wedge v in SYM^(m-2)\tensor Wedge
    vw=vector((binomial(t,2)+(s-t)*t)*binomial(t+level-3,level-2),n,0);
    if(level>2, veclim=vector(level-2,j,[1,t]);,veclim=vector(1,j,[1,1]));
      lauf=0;
      forvec(a=veclim,
        if(a==vecsort(a),
          valfactor=prod(j=1,level-2,u[a[j]]);
         for(i1=1,t-1,
            for(i2=i1+1,t,lauf=lauf+1;vw[lauf]=valfactor*(u[i1]*v[i2]-u[i2]*v[i1]);
            )
          );
          for(i1=1,t,
            for(i2=t+1,s,
              vw[lauf++]=valfactor*u[i1]*v[i2]
            )
          );
        ,valfactor=1;
        );
      );vw;
}

prenull(m1)=
{ local(i,j,dims);
  dims=matsize(m1);
  for(i=1,dims[1],
    for(j=1,dims[2],
      if(abs(m1[i,j])<10^-50,
        m1[i,j]=0);
    );
  );
  m1;
}

{
  covol(mm)=
\\ gives the covolume of the image of the real matrix   mm
    dims=matsize(mm);
    mm=prenull(mm);
    imm=matimage(mm);
    imagemat=if(matrank(imm)>0,matsupplement(imm),matid(dims[1]));
    prerat=1/imagemat*mm;
    mmrat=matrix(dims[1],dims[2],j,k,rat(prerat[j,k]));
    hm=mathnf(mmrat);
    if(matsize(hm)[1]==matsize(hm)[2]&&hm<>[;],dh=matdet(hm),dh=0);
    abs(matdet(imagemat))*dh;
}


{
  makerationallattice(mm)=
\\ input a real matrix  mm  whose image has (presumably) a  Z-structure
\\ completes the image to obtain a square matrix  matsupplement(imd) 
\\ and multiplies its inverse to mm
    dims=matsize(mm);
    imd=matimage(mm~);
    if(length(imd)==0,imdinv=matid(dims[2]),
      if(length(imd)==length(imd~),
        if(matdet(imd)>epsi,imdinv=1/imd,imdinv=1/matsupplement(imd));
      ,imdinv=1/matsupplement(imd)));
    mm*imdinv~;
}

{
  rationallattice(mm,blocks,eps)=
\\ 
    ms=matsize(mm);
    len=ms[2]/blocks;
    blockvec=vector(blocks,j,0);  
    mx=matrix(ms[1],ms[2],j,k,if(abs(mm[j,k])<eps,0,mm[j,k]));
    for(jj=1,blocks,
      blockvec[jj]=makerationallattice(matrix(ms[1],len,j,k,mx[j,k+(jj-1)*len]))
    );
    result=matrix(ms[1],ms[2],j,k,rat(blockvec[(k-1)\len+1][j,(k-1)%len+1]));
    result;
}

{
  good_kernel_base(mm,n)=
\\ forms and improves a basis for the kernel of mm   n  times
    if(matsize(mm)[1]>1&&matrank(mm)>0,
\\print("in qflll");
\\write("bug",mm~/content(mm));
      better=qflll(mm~/content(mm),4)[1];
\\print("out qflll");
      if(length(better)>1,
        for(j=1,n,print1("improving ");better=better*qflll(better,1));
      )
    ,better=mm);better;
}

{
 idlist(bnf,list, lst, j)=lst=[];
\\ creates the list of ideals above list
    for(j=1,length(list),
      lst=concat(lst,idealprimedec(bnf,list[j])));lst
}

{ 
  takepolylog(m)=
\\ evaluates the m-log on  x_i  and finds the real lattice given by the
\\ linear combinations from the kernel and  turns it into a  Z-structure
\\ i.e. gives an integer matrix which is a  
    print(" taking ",m,"log ("level")");
\\  extra=factorial(m)/2/Pi;
  extra=1;
  extra=factorial(m)/(2*Pi)^(m-1);

    symdim=binomial(lendec+level-m-1,level-m);
    if(even(m),nhyperbolic=nminus; nspherical=nplus;
       ,nhyperbolic=nplus; nspherical=nminus;
    );

    polyl=matrix(triples,symdim*nhyperbolic,j,k,0);
\\    print(" polyl ",level," (",m,"): ");
    if(level>m,
      veclim=vector(level-m,j,[1,lendec]),
      veclim=vector(1,j,[1,1])
    );
\\    print(veclim);
    valfactor=vector(symdim,j,0);
    for(j=1,triples,lauf=0;
      if(even(m),
        vp=vector(nhyperbolic,l,extra*polylog(m,conjvec(m2[j])[2*l+r1],3));
        ,vp=vector(nhyperbolic,l,
        if(l>r1,
          extra*polylog(m,conjvec(m2[j])[2*(l-r1)+r1],3),
          extra*polylog(m,conjvec(m2[j])[l],3)));
      );
      u=freedec_x[j,];
      forvec(a=veclim,
        if(a==vecsort(a),
          valfactor[lauf++]=prod(k=1,level-m,u[a[k]])
        ,)
      );
    polyl[j,]=tensor(vp,valfactor);
  );


  if(m<level,
    ratlat=rationallattice(wedge_relations~*polyl,symdim,epsi);
    contvec=vector(length(ratlat),j,content(ratlat[,j]));
    integral_lattice=matrix(length(ratlat~),length(ratlat),j,k,
        if(contvec[k]<>0,ratlat[j,k]/contvec[k],0));
    ker_of_lattice=good_kernel_base(integral_lattice,2);
    wedge_relations*=ker_of_lattice;
\\   print1(" ",matsize(wedge_relations)," ",truncate(log(1+norml2(wedge_relations))));
    rk=matrank(matimage(ker_of_lattice)); 
print(" rank still equal to  "rk);
   ,

wedge_relations=wedge_relations*qflll(wedge_relations,1);
wedge_relations=wedge_relations*qflll(wedge_relations,1);
  lat=wedge_relations; ww2=2*ww;
  for(j=1,length(wedge_relations),
    tor=(torsionwedge*wedge_relations[,j])%(ww2);co=lift(gcd(ww2,content(tor)));
    if(tor<>0,wedge_relations[,j]=(ww2)/co*wedge_relations[,j],);
  );

  regulator_lattice=polyl~*wedge_relations;
  regul=covol(regulator_lattice);
  rk=matrank(matimage(regulator_lattice));
  );
}

{
  Belong(a,level, e)=
\\ checks if one of the associates of  a under the S_3 -action already 
\\ occurs in the list, for level=2  there are  6 possibilities, for others 2
    ind = setsearch(m1, a, 1);
    if(level==2,  ind == 0 || setsearch(m1, 1 - a) ||
      setsearch(m1, e = 1/a) || setsearch(m1, 1 - e) || 
      setsearch(m1, e = 1/(1-a)) || setsearch(m1, 1 - e),
     ind == 0   ||  setsearch(m1,  1/a)
    )
}

drop(v,r)=concat(vector(r-1,j,v[j]),vector(length(v)-r,j,v[j+r]));

{
 decompose(alpha,lst)= 
\\ associates to an element alpha its exponent vector in the S-unit base
   if(lst==idealp,
      bnfissunit(bnf,S_units_p,alpha),
      bnfissunit(bnf,S_units_pq,alpha)
   );
}

{
  issmooth(n, listp)=
  no = abs(n);
  no==1 ||  no == prod(j=1, length(listp), listp[j]^valuation(no,listp[j]))
}

{
  search_integral_S_units(veclim_add,candmax,pol,np)=
\\ creates a list of S-units and then checks which differences are also S-units
    candnum=vector(candmax,j,0);
    candcont=vector(candmax,j,0);
    m1 = listcreate(candmax+2*ww);
    lauf=0; maxlauf = 30*binomial(lendec,2); 

    triples=0;
    if(fu_number,fu_number = min(fu_number, r1+r2-1);
       forvec(a=veclim_mult, 
         if(lauf<candmax,
           elt=prod(j=1,fu_number, fundunits[j]^a[j]);
           tu=bnf.tu[2]; 
           for(i=1,2*ww,elt*=tu;
             if(issmooth(norm(1-elt), listpq),
               if(!Belong(elt,level),triples++; listinsert(m1, Str(elt), ind))
             );
           )
         );
       );
     );

    if(!fu_number,
      forvec(a=veclim_add, if (a==0, next);
        if(lauf<candmax, alpha = nfbasistoalg(bnf, a~);
          if(issmooth(norm(alpha),listp),lauf++;
            candnum[lauf]=alpha;candcont[lauf]=content(a);
            if (lauf == maxlauf, break)
          ,);
        ,)
      );
 if(level==3,
   if(setsearch(Set(listp),"3")<>0,
      triples++;listinsert(m1,eval(Str("Mod(-3,"bnf.pol")")),triples);
      triples++;listinsert(m1,eval(Str("Mod(3,"bnf.pol")")),triples);
   );
   if(setsearch(Set(listp),"2")<>0,
      triples++;listinsert(m1,eval(Str("Mod(-1,"bnf.pol")")),triples);
   );
 );
      for(j=1,lauf,cn=candnum[j];
        for(k=j+1,lauf,cd=candnum[k];
          if((triples<candmax),
            if((gcd(candcont[j],candcont[k])==1),
              cq=cn/cd;
              if(issmooth(norm(cd-cn),listpq),
                if(!Belong(cq,level),triples++; listinsert(m1, Str(cq), ind)));
              if(issmooth(norm(cd+cn),listpq),
                if(!Belong(-cq,level),triples++; listinsert(m1, Str(-cq), ind)));
            ,)
          ,)
        );
      );
    );
    print("found "triples"  S-units,  [wanted: >= "candmax" ]");
    m2=vector(triples,j,eval(m1[j]));
}


{
  newsplitprime(last)=forprime(j=last+1,1000,
      if (length(idealfactor(bnf,j)[,1])>1, erg=j;break()));erg;
}        

{
  numberfielddata(pol,listp_orig,extrap,listq_orig,extraq)=   
\\ creates the lists of primes and prime ideals  and collects further 
\\ necessary data for the number field 
    deg=length(pol)-1;
    bnf=bnfinit(pol);  classnumber=bnf.clgp[1];

    listp=listp_orig;
    if(length(listp)>0,bound = listp[length(listp)],bound = 1);
    for(j=1,extrap, 
        listp = concat(listp, newsplitprime(bound));
        bound=listp[length(listp)];
    );

    listq=listq_orig;
    listq=vecsort(eval(setunion(Set(listq),Set(listp))));
/*  throw away superfluous primes which occur in both p and q  */

    for(j=1,extraq,
        listq = concat(listq, newsplitprime(bound));
        bound=listq[length(listq)];
    );

    listpq=listq;
    r1=bnf.sign[1]; r2=bnf.sign[2]; d=bnf.disc; 
    fundunits=bnf.fu; ww=bnf.tu[1]/2; 
    nplus=r1+r2; nminus=r2;
       print(pol"  , d = "d" , signature = "bnf.sign);
    
    idealp=idlist(bnf,listp);
    S_units_p=bnfsunit(bnf,idealp);

    idealpq=idlist(bnf,listpq);
    S_units_pq=bnfsunit(bnf,idealpq);

} 

{ 
 getmax(deg)=
\\ heuristic bounds for the coefficients in additive search
   if(deg==1,amax=2500;bmax=30;,
     if(deg==2,amax=500;bmax=30;,
       if(deg==3,amax=50;bmax=15,
         amax=8;bmax=5;
\\       amax=18;bmax=8
       )
     )
   );[amax,bmax];
}


{
   constants(dummy)=
     epsi=10^-(prc-20);
     eps1=10^-(prc\2);
     if(deg<=2,bd1=20;bd2=10,if(deg==3,bd1=20;bd2=10,bd1=20;bd2=8));
     if(deg<=3,sz=50,sz=20);
  sz=20;
     veclim_mult=vector(fu_number,j,[round(-log(2^sz)/log(j+2)^2/2),round(log(2^sz)/log(j+2)^2/2)]);

     amax=getmax(deg)[1];
     bmax=getmax(deg)[2];

     degmax=6;

   veclim_add=vector(deg,j,if(j==1,[1,amax],if(j<degmax,[-bmax,bmax],[0,0])));
   lendec=length(idealp)+r1+r2-1;
   lendecpq=length(idealpq)+r1+r2-1;
}

{getktor(nf,nn)=
  longfactbase=[];
  for(j=1,50,idfac=idealfactor(bnf,prime(j));
    for(k=1,length(idfac~),
	longfactbase=concat(longfactbase,[idfac[k,1]]);
    )
  );
  ktor=0;lfb=length(longfactbase);
  for(j=1,lfb,no=idealnorm(nf,longfactbase[j]);
    ktor=gcd(ktor,no^10*(no^nn-1))
  );ktor;}



even(z)=z%2==0;

odd(z)=z%2==1;

{
  getknownfactor1(level)=
    wmf=getktor(bnf,level);
    if(even(level)==0,	
      knownfactor1=abs(Pi^((level-1)*(r1+r2))*(sqrt(abs(d))/2^r2/Pi^(deg/2))^(2*level-1)*wmf
         *gamma(level/2)^r1*gamma(level)^r2);
    , knownfactor1=abs(Pi^((level-1)*(r2))*(sqrt(abs(d))/2^r2/Pi^(deg/2))^(2*level-1)*wmf
         *(gamma(level/2)/gamma((1-level)/2))^r1*gamma(level)^r2););
}


{
  getknownfactor0(level)=
    wmf=getktor(bnf,level);
    deg=r1+2*r2;
    aa=2^(-r2)*Pi^(-deg/2)*sqrt(abs(d));
    if(level%2==1,
      zneg=aa^(2*level-1)/2^r1*(gamma(level/2)*gamma((level+1)/2))^r1*gamma(level)^(2*r2);
      knownfactor=wmf*zneg*Pi^((r1+r2)*(level-1))/(Pi)^(r1+r2)/2^(2*r2);
      , 
     zneg=aa^(2*level-1)*(gamma(level/2)/gamma((level+1)/2))^r1*gamma(level)^(2*r2);
     knownfactor=wmf*zneg*Pi^(r2*(level-1))/(Pi)^(r2)/2^(2*r1);
    );
     knownfactor;
}

getknownfactor(level)=
{\\ compute the "known factors" (Pi, discriminant of bnf, Gamma factors) which occur via the 
\\functional equation of the Dedekind zeta function as well as  K_{2n-1}^tor 
    wmf=getktor(bnf,level);
    deg=r1+2*r2;
    aa=2^(-r2)*Pi^(-deg/2)*sqrt(abs(d));
    if(level%2==1,
      zneg=aa^(2*level-1)/2^r1*(gamma(level/2)*gamma((level+1)/2))^r1*gamma(level)^(2*r2);
      knownfactor=wmf*zneg;
\\*Pi^((r1+r2)*(level-1));
    , 
      zneg=aa^(2*level-1)*(gamma(level/2)/gamma((1-level)/2))^r1*gamma(level)^(2*r2);
      knownfactor=wmf*zneg;
\\*Pi^(r2*(level-1));
\\  BEST:  1/(2*Pi)^(n*m)*Pi^(r1+r2)*2^r2*gamma(m)^n*sqrt(abs(d))^(2*m-1)
\\ n=deg;  m= level odd, for m even  interchange role of r2 and r1+r2;

    );
    knownfactor;
}

{
  getzetavalue(level)=
    default(realprecision,18);
    if(abs(d)<10^6,
      zetavalue=zetak(zetakinit(bnf.pol),level);
    ,bd=10000;
      dircoeffs=dirzetak(bnf,bd);
      zetavalue=0.+sum(j=1,bd,dircoeffs[j]/j^level)
    );zetavalue
}

{
  showme(level,fn1,fn2)= 
     indexappr=bestappr(zetavalue/regul*knownfactor,100);
     dyadicprimes=length(idealfactor(bnf,2)~);
     if(indexappr<10^10,index=round(indexappr),extrap=1);
     write(fn1,""2*level-2": {"bnf.pol","indexappr","length(p)","ww","[r1,r2]
      ","d","d%9"(9),"d%16"(16),"dyadicprimes",cl "classnumber"}, ord:"
      ,factor(index)"  d: "factor(abs(d))" ("triples") "concat(listp,[])
       " "eval(setminus(Set(listpq),Set(listp))));
     write(fn2,"{"2*level-2",{"r1","r2"},"d","ww","dyadicprimes","index"},");
}

{
  showmescreen(level)=
     indexappr=abs(zetavalue/regul*knownfactor);
     dyadicprimes=length(idealfactor(bnf,2)~);
     if(indexappr<10^10,index=round(indexappr),extrap=1);
     print("K"2*level-2" for "bnf.pol"  is conjecturally equal (up to powers of 2) to");
     print("  *** "rat(indexappr)" ***    ="factor(rat(indexappr)));
     print();
\\,"length(listp)","ww",",[r1,r2],
\\      ","d","d%9"(9),"d%16"(16),"dyadicprimes",cl "classnumber"}, ord:"
\\      ,factor(bestappr(indexappr,10^8))"  d: "factor(abs(d))" ("triples") "concat(listp,[])
\\       " "eval(setminus(Set(listpq),Set(listp))));
\\     print("{"2*level-2",{"r1","r2"},"d","ww","pic(bnf)","index"},");
}

{
  showmescreen_old(level)= 
     indexappr=zetavalue/regul*knownfactor;
     dyadicprimes=length(idealfactor(bnf,2)~);
     if(indexappr<10^10,index=round(indexappr),extrap=1);
     print("K"2*level-2": {"bnf.pol","indexappr","length(listp)","ww",",[r1,r2],
      ","d","d%9"(9),"d%16"(16),"dyadicprimes",cl "classnumber"}, ord:"
      ,factor(bestappr(indexappr,10^8))"  d: "factor(abs(d))" ("triples") "concat(listp,[])
       " "eval(setminus(Set(listpq),Set(listp))));
     print("{"2*level-2",{"r1","r2"},"d","ww","pic(bnf)","index"},");
}

{
 factnum(u, p,e,str)=
     if(u<0, str="-";u=-u, str="");
     u=factor(u);
     p=u[,1]; e=u[,2];
     for(j=1,length(p),str=Str(str "" p[j]"^{"e[j]"} "));
     str;
}

{
  showmetex(level,fn)= 
     k2order=bestappr(zetavalue/regul*knownfactor,100);
     write(fn,factnum(d) "\t&  "d%9 "," d%16" \t&    "
         k2order" \t&  "factnum(k2order)" \\cr");
}

{bloch_element(m)=
  if(r2==1,
    for(j=1,length(wedge_relations~),wr=wedge_relations[j,m];
      if(wr>0,print1("+",wr,"[",lift(m2[j]),"]")
        ,
        if(wr<0,print1(wr,"[",lift(m2[j]),"]")
         ,
        )
      );
    );
    print();
    print(" regulator * ",rat((polyl[,1]~*wedge_relations[,m])/regul));
    print();
    print(polyl[,1]~*wedge_relations[,m]);
    print();
  ,)
}

{bloch_dilog(m)=
  if(r2==1,
    for(j=1,length(wedge_relations~),wr=wedge_relations[j,m];
      if(wr>0,print1("+",wr," ",polyl[j,]," ")
        ,
        if(wr<0,print1(wr," ",polyl[j,]," ")
         ,
        )
      );
    );
    print();
    print(" regulator * ",rat((polyl[,1]~*wedge_relations[,m])/regul));
  ,)
}
 
{bloch_vec(m, fact=1)=
  if(r2==1,
    for(j=1,length(wedge_relations~),wr=wedge_relations[j,m];
      if(wr>0,print1("+",wr/fact," "vecwedge[j]" ")
        ,
        if(wr<0,print1("+",-wr/fact," "-vecwedge[j]" ")
         ,
        )
      );
    );
    print();
\\    print(" regulator * ",rat((polyl[,1]~*wedge_relations[,m])/regul));
  ,)
}
 
decompose_triples()=
{  print("decompose the S-units on the given base");
   fulldec_x=matrix(triples,lendec+1,j,k,decompose(m2[j],idealp)[k]);
   dropvec_x=vector(triples,j,drop(fulldec_x[j,],r1+r2));
   freedec_x=matrix(triples,lendec,j,k,dropvec_x[j][k]);
   torsdec_x=matrix(triples,1,j,k,fulldec_x[j,r1+r2]);

   fulldec_1minusx=matrix(triples,lendecpq+1,j,k,decompose(1-m2[j],idealpq)[k]);
   dropvec_1minusx=vector(triples,j,drop(fulldec_1minusx[j,],r1+r2));
   freedec_1minusx=matrix(triples,lendec,j,k,dropvec_1minusx[j][k]);
   torsdec_1minusx=matrix(triples,1,j,k,fulldec_1minusx[j,r1+r2]);
}


{liza(pol,listp_orig=[2,3],extrap=0,listq=[],extraq=0,candmax=70,prc=150,level=2,fu_number=0)=

  print();
  print("Approximating the order of K_"2*level-2" O_F, where");
  print1(" minpol(F) = ");
   default(realprecision,prc);
   setrand(2);

   numberfielddata(pol,listp_orig,extrap,listq,extraq);

   constants();
if(SLOW,print1("?");input(););
   search_integral_S_units(veclim_add,candmax,pol,listp);
if(SLOW,print1("?");input(););
   if(triples<5, return);

   decompose_triples();
if(SLOW,print1("?");input(););
  print(" taking tensor product of exponent vectors ");
   torsionwedge=vector(triples,j,
       torsionpart(fulldec_x[j,],fulldec_1minusx[j,],lendec,lendecpq));
   vecwedge=vector(triples,j,
       symwedge(dropvec_x[j],dropvec_1minusx[j],lendec,lendecpq,level));

   matwedge=vec_to_mat(vecwedge);
   print("size of matrix for LLLkerim: "matsize(matwedge));
   wedge_relations=qflll(matwedge~,4)[1];
if(SLOW,print1("?");input(););

   merk=vector(level,j,0);
   if(level>2 && length(wedge_relations~)>2,
     for(i=1,2, wedge_relations=wedge_relations*qflll(wedge_relations,1));
        print(" qflll(,4) finished.  ");
   );

   rk=matrank(wedge_relations);
   regul=0;lev=1;
   while(lev<level&&rk>0,
     lev++;
     if(r2>0||odd(lev), takepolylog(lev))
   );
   default(realprecision,50);
   print;print("regulator   "regul);
if(SLOW,print1("?");input(););
   if(abs(regul)<epsi, return);
   getknownfactor(level);
   zetavalue=getzetavalue(level); 

/*  the following gives an output for conjectural order of K_2 and 
    supplementary data: */

   if(rk==nhyperbolic, 
     if(triples>tripmin, 
       showmescreen(level);
\\       showmetex(level,filenametex);
      ,showmescreen(level);
\\       showmetex(level,filenametex_todo)
     );
   );
}


{
 k2(disc,vec_of_primes=[2,3])=liza(x^2-disc,vec_of_primes,0,[],0,150,120,2);
}


/*  default for bound on number of  exceptional  S-units  */
tripmin=150;





/*  COMMENTS+++COMMENTS+++COMMENTS+++COMMENTS+++COMMENTS+++COMMENTS+++

      typical call of  "Lichtenbaum-Zagier", short  "liza"  :

  pol=x^3-x-1;  \\ the polynomial defining the number field 
  vecp=[2];   \\ prime factors  of the  x_i:  the primes above  {2}
  extrap=0;   \\ let the program choose extrap more primes for the x_i
  vecq=[3];   \\ additional primes which are allowed to divide  1-x_i
  extraq=0;   \\ let the program choose extraq more primes for the 1-x_i

  bnd=150;    \\ default for bound on number of  exceptional  S-units
  prec=100;   \\ precision for (most of) the computations

  polylog_level=2;    \\ the program also works for higher polylogs 

  fu_number =0;       \\ take all fundamental units

\\ Then a typical call would look as follows:

liza(pol, vecp,extrap, vecq,extraq, bnd, prec, polylog_level, fu_number)



\\ or, if one wants to concentrate on imaginary quadratic fields and 
\\ the dilogarithm, 
  
k2(-23,vecp)

*/

/*  A typical in-  and  output, e.g. in the case    d=-404, is

INPUT:    liza(x^2+101,[2,3,5,7])

OUTPUT:    
regulator   1560.6998678495772023190002920566127302750087745273
K2: {x^2 + 101,1.00000000000000000,4,1,[0, 1],-404,1(9),12(16),1,cl 14}, ord:matrix(0,2,j,k,0)  d: [2, 2; 101, 1] (69) [2, 3, 5, 7] []
   
{x^2 + 101,                input polynomial for F
1.00000000000000000,       ratio  zeta_F(2)* known factors / regulator
                  =conjectural order of K_2 O_F (if enough primes)
4,                         number of primes                     
1,                         number of rootsof1
[0, 1],                    signature 
-404,                      discriminant
1(9),                      discriminant mod  9
12(16),                    discriminant mod  16
1,                         number of primes above 2
cl 14},                    class number 
ord:matrix(0,2,j,k,0)      factored conjectural order
d: [2, 2; 101, 1]          factored discriminant
(69)                       number of exceptional S-units found
[2, 3, 5, 7]               S  for building  x_i
[]                         S' - S  (extra primes for building  1-x_i )
   

COMMENTS:  Too many primes will make matrices too large to handle
conveniently  (e.g. qflll  for  300x500 is still possible on a large machine
but takes quite long)

For small discriminants  2  or  3  primes are often enough.
Try e.g. 
   for(j=2,50,if(isfundamental(-j),k2(-j)))

The procedure

    bloch_element(m)

gives as output the  m-th  element in the Bloch group (the ordering coming
from the relation matrix), together with its multiple of the regulator
and its actual dilog value.

Example:

INPUT:  bloch_element(6) 

OUTPUT:  
+1756[-1/10*x - 1/5]+6140[-1/14*x - 23/14]+1964[-1/14]+2336[-1/15]+5100[-1/2]+4900[-1/20*x + 13/20]-3080[-1/245*x - 12/245]+3276[-1/252*x - 5/252]-1140[-1/3*x - 2/3]+5968[-1/35*x + 12/35]+2336[-1/36*x + 13/36]-2564[-1/4]-5796[-1/42*x + 20/21]-3104[-1/49*x + 37/49]+1384[-1/5*x - 1/10]-2544[-1/5]+344[-1/525*x - 2/525]-1168[-1/567*x + 40/567]+2692[-1/6]+7688[-1/63*x + 40/63]-6024[-1/63*x - 5/63]+1964[-1/7]+1190[-2/105*x + 16/21]-264[-2/3*x - 1/3]

		               the element

 regulator * -2                the (integral) multiple of the regulator
-3121.39973569915440           the actual value

*/