\\filename : finite.gp

\\000000000000000000000000000000000000000000000000000000000000000000000000
\\00000000                                                        00000000
\\00000000                 preamble                               00000000
\\00000000                                                        00000000
\\000000000000000000000000000000000000000000000000000000000000000000000000

\\include the file ei-linear-alg.gp, which contains the basic linear algebra operations.

\r ei-linear-alg



\\111111111111111111111111111111111111111111111111111111111111111111111111
\\11111111                                                        11111111
\\11111111   positive definite inner products, reflections and    11111111
\\11111111   their matrices                                       11111111
\\111111111111111111111111111111111111111111111111111111111111111111111111


\\generate a random eisenstein vector of length n

ran(n,seed)=
{
local(i,vv);
vv = matrix(n,2);
for( i = 1 , n, vv[i,1] = random(2*seed-1)-seed; vv[i,2] = random(2*seed-1) - seed);
return(vv);
} 

\\change an eisenstein unit to a character ( for better viewing)

cha(u)=
{
local(uu); if ( u ==[ 0, 0] , uu = " 0"); 
if ( u ==[ 1, 0] , uu = " 1");if ( u ==[ 1, 1] , uu = "-W");if ( u ==[ 0, 1] , uu = " w");
if ( u ==[-1, 0] , uu = "-1");if ( u ==[-1,-1] , uu = " W");if ( u ==[ 0,-1] , uu = "-w");
return(uu);
}
\\change a vector similarly as above ( for better viewing)

chav(v)= 
{
local(i,l,vv); l = matsize(v)[1];vv = vector(l); for( i = 1, l, vv[i] = cha(v[i,]) ); return(vv);
}

\\positive defininte standard hermitian form on an eisenstein lattice

ipd(x,y)= 
{
local(ll,i,ipp); ll = matsize(x)[1]; ipp = 0; for( i = 1, ll, ipp = ipp + dot( bar(x[i,]),y[i,]) ); return(ipp);
}

\\ the norm of a vector 

n2d(x) = ipd(x,x)[1];

\\given a set of vectors e, calculate the gram matrix using the positive definite inner product.

checkipd(e)=
{
local(i,j,ll,gram);
ll = matsize(e)[2];
gram = matrix(ll,ll);
for( i = 1, ll, for ( j = 1, ll, gram[i,j] = ipd(e[i],e[j])  ));
return(gram);
}

\\the order 3 reflection and the order 2 reflection in the vector x applied to the vector a

phid3(x, a) = a-sdot(dot([1,-1],ipd(x,a))/ipd(x,x)[1],x);
phid2(x, a) = a-sdot(dot([2, 0],ipd(x,a))/ipd(x,x)[1],x);

\\ order k reflection in the vector x applied to a (k = 2 or 3)
\\the definition of phid below is made to work for E6, E8, G26, K12, K10, but not for G5 because in G5, we have order 3 reflection in the norm 6 vectors


phid(x,a)=
{
local(ph); ipp = ipd(x,x)[1]; if(ipp==3, ph = phid3(x,a)); if(ipp==2 || ipp==6, ph = phid2(x,a)); return(ph);
}

\\mad(x) is the matrix of the linear transformation phid(x, *).

mad(x)=
{
local (i,r,matr,n,unitvec);
n=matsize(x)[1];
matr=matrix(n,n);
unitvec=matrix(matsize(x)[1],matsize(x)[2]); 
for (r=1,n,\
  for ( i = 1, n, unitvec[i,] = [0,0]);\
  unitvec[r,] = [1,0];\
  ph=phid(x,unitvec);\
  for ( i = 1, n , matr[i,r] = ph[i,]) );
return(matr);
}

\\if x is a sequence of roots, pr(x) gives the product of the matrix of reflections in them
\\i.e. the matrix of phi_x1 o phi_x2 o .... o phi_xn

prd(x)=
{
local(lm,prr);
lm=matsize(x)[2];
prr= mad(x[1]);for(i=2,lm, prr = md(prr,mad(x[i]) ));
return(prr);
}


\\ the sine of the squared spherical distance between the hyperplane x-perp and the vector w.

distd(x,w) = n2( ipd( x, w) )[1]/(n2d(x)*n2d(w) )  ;




\\222222222222222222222222222222222222222222222222222222222222222222222222
\\22222222                                                        22222222
\\22222222   generate the projective roots of a lattice           22222222
\\22222222                                                        22222222
\\222222222222222222222222222222222222222222222222222222222222222222222222


\\*************************    Eisenstein E_8 lattice   **********************************


gre8()=
{
local(al,i,j,cou, re8, re8o);
al = vector(3);re8 = vector(40); re8o=vector(40);
cou = 0; al[1] = [1,0]; al[2]= [0,1]; al[3] = [-1,-1];
for( i = 1, 3, for( j = 1, 3, \
cou = cou+1;\
 re8[cou   ] = [ 0 , 0 ;  1 , 0 ;  1 , 0 ;  1 , 0 ]; \
re8[cou   ][3,] =dot( al[i],re8[cou   ][3,]); re8[cou   ][4,] =dot( al[j],re8[cou   ][4,]);\
re8[cou+ 9] = [ 1 , 0 ;  0 , 0 ;  1 , 0 ; -1 , 0 ];\
re8[cou+ 9][3,] =dot( al[i],re8[cou+ 9][3,]); re8[cou+ 9][4,] =dot( al[j],re8[cou+ 9][4,]);\
re8[cou+18] = [ 1 , 0 ;  1 , 0 ; -1 , 0 ;  0 , 0 ];\
re8[cou+18][2,] =dot( al[i],re8[cou+18][2,]); re8[cou+18][3,] =dot( al[j],re8[cou+18][3,]);\
re8[cou+27] = [ 1 , 0 ; -1 , 0 ;  0 , 0 ;  1 , 0 ];\
re8[cou+27][2,] =dot( al[i],re8[cou+27][2,]); re8[cou+27][4,] =dot( al[j],re8[cou+27][4,]);
));
for( i = 1, 4, re8[36+i] = [ 0 , 0 ;  0 , 0 ;  0 , 0 ;  0 , 0 ]; re8[36+i][i,] = [2,1]);
for(i=1,40, re8o[i] = 3);
return([re8,re8o]);
}


\\************************* A lattice with reflection group is G_26  *********************


\\ the roots of G26 were found by starting with the simple roots 
\\G26[1]=[ 0, 0; 1, 0; 1, 0; 1, 0];G26[2]=[ 0, 0; 0, 0; 0, 0;-2,-1];G26[3]=[ 1, 0; 2, 0; 1, 0; 0, 0] forming the diagram (2)==(3)--(3). 
\\and applying the algorithm genrofromx(G26,3,21); these projective are listed below: first the norm 3 roots, then the norm 6 roots.
\\it has 12 projective roots of norm 3 in which we have order 3 reflections and 9 projective roots of length 6 in which we have order 2 reflections 
\\the lattice is contained in Eisenstein E8 and contains Eisenstein E6.The short roots form the E6 root system.



gr26()=
{
local(rg26,rg26o); rg26=vector(21); rg26o  = vector(21);
rg26[ 1] = [ 0, 0; 1, 0; 1, 0; 1, 0]; rg26[ 2] = [ 0, 0; 0, 0; 0, 0;-2,-1]; rg26[ 3] = [ 0, 0; 1, 0; 1, 0;-1,-1]; rg26[ 4] = [ 0, 0; 1, 0; 1, 0; 0, 1];
rg26[ 5] = [-1, 0;-1, 0; 0, 0; 1, 0]; rg26[ 6] = [-1, 0;-1, 0; 0, 0;-1,-1]; rg26[ 7] = [-1, 0;-1, 0; 0, 0; 0, 1]; rg26[ 8] = [-1,-1; 0,-1; 1, 0; 0, 0];
rg26[ 9] = [ 0, 1; 0, 0; 0,-1;-1,-1]; rg26[10] = [-1, 0; 0, 0; 1, 0; 1, 1]; rg26[11] = [ 0, 1; 1, 1; 1, 0; 0, 0]; rg26[12] = [-1, 0; 0, 0; 1, 0;-1, 0];
rg26[13] = [ 1, 0; 2, 0; 1, 0; 0, 0]; rg26[14] = [ 1, 0; 1, 1; 0, 1;-1, 1]; rg26[15] = [ 1, 0; 1, 1; 0, 1; 2, 1];
rg26[16] = [ 1, 0; 1, 1; 0, 1;-1,-2]; rg26[17] = [ 1, 0; 0,-1;-1,-1;-2,-1]; rg26[18] = [ 1, 0;-1, 0;-2, 0; 0, 0];
rg26[19] = [ 0, 1; 1, 1; 1, 0;-2,-1]; rg26[20] = [-1,-1;-1, 0; 0, 1;-2,-1]; rg26[21] = [ 2, 0; 1, 0;-1, 0; 0, 0];
for( i = 1, 12, rg26o[i] = 3); for( i = 13, 21, rg26o[i] = 2);
return([rg26,rg26o]);
}


\\************ An Eisenstein lattice whose reflection group is G26   *********************
\\the form of the lattice given below is not modular

gr26var()=
{
local(i, rg26, rg26o, re6); rg26=vector(21); rg26o = vector(21);
re6 = gre6;
for( i = 1, 12, rg26[i] = re6[1][i] );
rg26[13]=[ 0, 0; 1, 0; -1, 0];  rg26[16] = [ 1, 0; 0, 0; -1, 0]; rg26[19] = [ 1, 0; -1, 0; 0, 0];
rg26[14]=[ 0, 0; 1, 0;  0,-1];  rg26[17] = [ 1, 0; 0, 0;  0,-1]; rg26[20] = [ 1, 0;  0,-1; 0, 0];
rg26[15]=[ 0, 0; 1, 0;  1, 1];  rg26[18] = [ 1, 0; 0, 0;  1, 1]; rg26[21] = [ 1, 0;  1, 1; 0, 0];
for( i = 1, 12, rg26o[i] = 3); for( i = 13, 21, rg26o[i] = 2);
return([rg26,rg26o]);
}


\\*************************    Eisenstein E_6 lattice   **********************************

gre6()=
{
local(al,re6, re6o,cou,i,j);
al = vector(3);re6= vector(12); re6o = vector(12);
cou = 0; al[1] = [1,0]; al[2]= [0,1]; al[3] = [-1,-1];
for( i = 1, 3, for( j = 1, 3, \
cou = cou+1;\
re6[cou   ] = [ 1 , 0 ;  1 , 0 ;  1 , 0 ]; \
re6[cou   ][2,] =dot( al[i],re6[cou   ][2,]); re6[cou   ][3,] =dot( al[j],re6[cou   ][3,]);\
));
for( i = 1, 3, re6[9+i] = [   0 , 0 ;  0 , 0 ;  0 , 0 ]; re6[9+i][i,] = [2,1]);
for( i = 1, 12, re6o[i] = 3);
return([re6,re6o]);
}


\\*********** sublattice of Eisenstein E_6 lattice: reflection group is G_5  *************

grG5()=
{
local(rG5, rG5o);
rG5 = vector(8); rG5o = vector(8);
rG5[1]=[ 1, 0; 1, 0; 1, 0];
rG5[2]=[ 0, 1; 0, 1; 1, 0];
rG5[3]=[-1,-1;-1,-1; 1, 0];
rG5[4]=[ 0, 0; 0, 0; 2, 1];
rG5[5]=[-2,-1;-2,-1; 0, 0];
rG5[6]=[ 0,-1; 0,-1; 2, 0];
rG5[7]=[ 0, 1; 0, 1; 2, 2];
rG5[8]=[-1,-1;-1,-1; 2, 2];
for( i  = 1, 8, rG5o[i] = 3);
return([rG5,rG5o]);
}

\\************************** the projective roots for G_4  *******************************

grG4()=
{
local(rG4, rG4o);
rG4 = vector(4); rG4o = vector(4);
rG4[1]=[ 1, 0; 1, 0; 1, 0];
rG4[2]=[ 0, 1; 0, 1; 1, 0];
rG4[3]=[-1,-1;-1,-1; 1, 0];
rG4[4]=[ 0, 0; 0, 0; 2, 1];
for( i  = 1, 4, rG4o[i] = 3);
return([rG4,rG4o]);
}

\\*************************    Eisenstein K_12 lattice  **********************************


grk12()=
{
local(rk12,rk12o,i,i1,i2,i3,j,k,j1,j2,cou,alph,bet);
rk12 = vector(126); rk12o = vector(126);
for( i = 1, 126, rk12[i] = matrix(6,2));
alph = [-1,0; 0,-1; 1,1];
bet = [0,1; -1,-1];
cou = 0;
\\..........................1 - 45.............................................
for( i = 1, 5, for( j = i+1, 6, for( k = 1, 3,\
cou = cou +1;\
rk12[ cou ][i,] =[1,0] ; rk12[cou][j,] = alph[k,] \
)));
\\..........................46 - 75............................................
for( i1 = 2, 6, for( j1 = 2,5, for(j2 = j1+1, 6, if ( j1 != i1 && j2 != i1 , \
cou = cou +1;\
rk12[cou ] = [1,0; 0, 1; 0,1; 0,1; 0,1; 0,1];\
rk12[cou][i1,] = [1,0]; rk12[cou][j1,] = rk12[cou][j2,] = [-1,-1];  \
rk12[cou] = sdivid(rk12[cou],t)\
  ))));
\\..........................76 - 95............................................
for( i1 = 2, 4, for( i2 = i1+1, 5, for( i3 = i2+1, 6,  for( k = 1, 2,\
cou = cou +1;\
rk12[cou ] = [1,0; 1,0; 1,0; 1,0; 1,0; 1,0];\
rk12[cou][i1,] = rk12[cou][i2,] = rk12[cou][i3,] = bet[k,]; \
rk12[cou] = sdivid(rk12[cou],t)\
    ))));
\\..........................96 - 125...........................................
 for( i1 = 1, 5, for( i2 = i1+1, 6, \
 cou = cou+1; rk12[cou] = [1,0; 1,0; 1,0; 1,0; 1,0; 1,0];\
rk12[cou][i1, ] = [0,1]; rk12[cou][i2,] = [-1,-1];\
rk12[cou] = sdivid(rk12[cou],t);\
cou = cou+1; rk12[cou] = [1,0; 1,0; 1,0; 1,0; 1,0; 1,0];\
rk12[cou][i2, ] = [0,1]; rk12[cou][i1,] = [-1,-1];\
rk12[cou] = sdivid(rk12[cou],t)\
));
\\..........................126................................................
rk12[126] = sdivid( [1,0; 1,0; 1,0; 1,0; 1,0; 1,0],t);
for( i = 1, 126, rk12o[i] = 2);
return([rk12,rk12o]);
}


\\*************************    Eisenstein K_10 lattice   *********************************

grk10()=
{
local(rk12, rk10, rk10o, cou, i,sk12n_1);
rk12 = grk12(); rk12 = rk12[1]; cou = 0; rk10 = vector(45); rk10o = vector(45);
sk12n_1 = [-1/3, 1/3; 2/3, 1/3; -1/3, 1/3; -1/3, 1/3; 2/3, 1/3; 2/3, 1/3]; \\a particular root of K12 whose perp is our K10
for( i = 1, 126, if( ipd( sk12n_1, rk12[i]) == [0,0], cou = cou+1; rk10[cou] = rk12[i] ));
for( i = 1, 45, rk10o[i] = 2);
return([rk10, rk10o]);
}




\\333333333333333333333333333333333333333333333333333333333333333333333333
\\33333333                                                        33333333
\\33333333  Find the set of simple roots for Eisenstein lattices  33333333
\\33333333                                                        33333333
\\333333333333333333333333333333333333333333333333333333333333333333333333


\\Calculate the Weyl vector 
\\input: r = a two component vector. 
\\First component = the set of projective roots of a root system of rank k,
\\Second componend= the order of reflection in these roots.
\\w = a vector in the ambient vector space.
\\output: the value of the function alpha(w) = \sum o(r)^(-2) * (<r,w>/|<r,w>| ) * (r/|r|), where the sum is over all projective roots.

nuweyl(x, w)= 
{
local(i,l,wey,ipp,nor,fac,odr,xr,xo); xr = x[1]; xo = x[2];
l = matsize(xr[1])[1];  wey = matrix(l,2);
for( i = 1, matsize(xr)[2], ipp=ipd( xr[i], w); nor = xo[i]^2*sqrt(n2(ipp)[1]* ipd(xr[i], xr[i])[1] );  \
fac = divid( ipp, [nor,0]);wey = wey + sdot( fac, xr[i]) );
return(wey);
}

\\ inputs: the set of projective roots r and a vector w in the ambient vector space and a number k.
\\ output: the k projective roots whose mirrors are closet to w. 

sim(r,w,k)=
{
local(num, simm, sv, i, sortv ); 
num = matsize(r)[2];  simm = vector(k); sv = vector(num); 
for( i = 1, num, sv[i] = distd(r[i],w)  ); 
sortv = vecsort(sv, ,1); 
for( i = 1, k , simm[i] = r[sortv[i]] );
return(simm);
} 

\\*********************************************************************************************

\\given a set of vectors in re8 find all the root obtained by repeated reflection in each other:

\\first write a hash function for roots of re8, called hfk


hfk(v)=
{
local(hfv, l ); l = matsize(v)[1]; v = 3*v; hfv = 1; for( i = 1, l, hfv = hfv + 2^(i-1)*divrem(v[i,1],2)[2] ); return(hfv);
}

\\given a set of roots x of length m, finds a larger set of roots (upto units) that can be obtained by repeated reflection of order t in the vectors x.
\\uses the hash function hfk
 
genrofromx(x,m,ttime)=
{
local(R,new,i,hh,ind,has, indmax,n,iold,jold,j,jj,ii,iii,k,flag,hashk,ij,newij);
\\ttime = 102400;
R=vector(ttime+1);new=matrix(14,2);has = matrix(1000,20);ind = vector(1000); 
for (i = 1 , m, R[i] = x[i];hh= hfk(R[i]);ind[hh]++;has[hh,ind[hh]]=R[i];write(rolist, " i[",i,"]=",i,"; j[", i,"]=",i, "; ro[",i,"]=",R[i],";" )); 
indmax=1; i=1; j=2; n = m;
while ( n < ttime, new = phid(R[i], R[j]);iold = i;jold = j;jj=3;\
if ( i-j == 1, j = i+1; i=1 ; jj = 1);if ( i < j  && jj == 3, ii=i ; i = j; j = ii; jj = 2);if (i > j+1 && jj == 3,  iii = j ; j = i ; i = iii+1;);\
for( ij = 1, 6, k=1;  flag=0; newij = sdot(uni[ij], new);hh = hfk(newij);\
while( k <= ind[hh] , if (newij == has[hh,k] ,  flag = 1); if(flag == 1, break(2)) ;k++)); 
if(  flag == 0 , n++ ; ind[hh] = ind[hh]+1 ; if(ind[hh] > indmax , indmax = ind[hh];print(indmax));\
has[hh,ind[hh]]=newij ; R[n]= new; print(n,"  ",i,"  ",j); write(rolist, " i[",n,"]=",iold,"; j[", n,"]=",jold, "; ro[",n,"]=",newij,";")));
return(R);
}


\\**************************************************************************************************************
\\one example of a simple system of E8

nwle8 = [-28.84738711905324514193592991, -87.85209334941779751737483596; -97.87061183257435692510971985, -108.7138238063479995647043409; 82.79755399688380249557439455, 104.1487451314407189547653939; -35.44892618629266535371691375, 81.38555433286938252816684416];

sime8 = vector(5);

sime8[1] = [-1, 0; 0, 0; -1, 0; 0, 1]; sime8[1] = sdot(uni[2], sime8[1]);
sime8[2] = [1, 0; 0, 1; 1, 1; 0, 0];  sime8[2] = sdot(uni[2], sime8[2]);
sime8[3] = [0, -1; 1, 0; 0, 0; 1, 1];
sime8[4] = [-1, 0; 0, 0; 1, 1; -1, -1]; 
sime8[5] = -vad(sime8[1], vad(sdot([1,-1], sime8[2]), vad(sdot([2,0],sime8[3]), sdot([1,-1], sime8[4]) ))); ; 

\\**************************************************************************************************************

\\return the valency vector of a gram matrix.
\\input:  a gram matrix of size (k,k) . 
\\output: a vector [v1, ... , vk] where the j th row conains v_j nonzero entries off the diagonal. 

valency(dip)=
{
local(val, i, j ,k, cou);
k = matsize(dip)[1];
val = vector(k);
for( i = 1, k, val[i] =-1; for( j = 1, k, if(  dip[i,j] != [0,0], val[i]++ )) );
print(" the valencies of the vertices in the diagram are: ",  val);
return(val);
}


\\given the gram matrix dip of the simple roots, check if they form the k12 diagram. 

checkk12d(dip)=
{
local(val, i, j, k );
val  = valency(dip);
if(vecsort(val) != [3,3,3,3,4,4], print("valency wrong"); return(0) );
j = vecsort(val, ,1)[5]; k = vecsort(val, ,1)[6]; \\ j and k are the two vertices with valency 4.
if( dip[j,k] == 0,print(" the two nodes with valency 4 are not connected"); return(0));
for( i = 1, 6, if( dip[i, j] == [0,0] && dip[i, k] == [0,0], print(" the two nodes with valency 4 does not connect to everything"); return(0) ));
print("correct diagram");
return(1);
}

\\given the gram matrix dip of the simple roots, check if they form the G26 diagram. 

checkg26d(dip)=
{
local(dit, i,j, pt,pn, fg);
fg=0;
dit = matrix(3,3); for( i = 1, 3, for( j = 1, 3, dit[i,j] = n2( dip[i,j] )[1] ));
for(j=1,6,pt=matrix(3,3);pn=numtoperm(3,j);for(i=1,3, pt[i,pn[i]]=1);if( pt*dit*pt^(-1) == [9,3,0;3,9,3;0,3,4], fg =1; return(fg)) );
fg = 0; 
return(fg);
}


\\given the gram matrix dip of the simple roots, check if they form the G5 diagram. 

checkg5d(dip)=
{
local();
fg=0; dit = matrix(2,2); for( i = 1, 2, for( j = 1, 2, dit[i,j] = n2( dip[i,j] )[1] ));
if( dit == [ 9, 12;12, 36] || dit==[36, 12; 12, 9], fg=1; return(fg) );
return(fg);
}


\\input :  r = a two component vector
\\First component = the set of projective roots of a root system of rank k, defined over Eisenstein numbers.
\\Second componend= the order of reflection in these roots.
\\n = number of trials. kk = minimal number of reflections needed to generate the reflection group.

\\ starts with a random vector and implements the algorithm to find a simple system n times.

\\let epsilon be the ratio of the norm of (w_n - w_{n-1}) and norm of w_n.
\\convergence criteria: iterate the sequence w_n until either epsilon is less than .00000001 or n exceeds 100. 
\\output: in each of the n trials print the following:
\\1. number of iteration done and the value of epsilon reached.
\\2. the gram matrix of the inner product between the simple roots.
\\3. the valency of each of the simple roots and whether the diagram is correct.
\\4. The Weyl vector and its norm.
\\5. After the n trials are complete it prints the ratio of trials in which the correct diagrams were obtained.

simk(r, n, kk)=
{
local(l,j,wl,wln,buf,cou,simm, ipdsim, checkdiagram, epsilon, correctcount, xr, lphistar, singularcount); xr = r[1]; lphistar=matsize(xr)[2];
print(lphistar);
l = matsize(xr[1])[1];  correctcount = 0; singularcount = 0;
for( j = 1, n, wl = ran( l,13847); wln=nuweyl(r,wl); cou=0; epsilon= ipd(vad( wl, -wln), vad(wl, -wln))[1]/ipd(wl,wl)[1]; \
until( cou > 100 || epsilon < .00000001, epsilon= ipd(vad( wl, -wln), vad(wl, -wln))[1]/ipd(wl,wl)[1];\
wl = wln; wln = nuweyl( r, wl );  cou++ );print("number of iteration = ", cou , ", epsilon =  ",epsilon); print(); \
simm = sim( xr,wln, kk);  ipdsim = checkipd(simm); \
print( " The gram matrix of the inner product between the simple roots are given below "); print(); \
for( j = 1, kk, print(ipdsim[j,]) ); print();\
print("det = ", mdet(ipdsim)); print(); if( mdet(ipdsim) == [0,0], singularcount = singularcount + 1);
if(lphistar == 126, checkdiagram = checkk12d(checkipd(simm))                             );\
if(lphistar == 45 , checkdiagram = (  vecsort(valency(checkipd(simm))) == [2,3,3,3,3])   );\
if(lphistar == 40 , checkdiagram = (  vecsort(valency(checkipd(simm))) == [1,1,2,2])     );\
if(lphistar == 21 , checkdiagram = checkg26d(checkipd(simm))                             );\
if(lphistar == 12 , checkdiagram = (  vecsort(valency(checkipd(simm))) == [1,1,2])       );\
if(lphistar == 8  , checkdiagram = checkg5d(checkipd(simm))                              );\
if(lphistar == 4  , checkdiagram = (  vecsort(valency(checkipd(simm))) == [1,1])         );\
correctcount=correctcount+checkdiagram; print( "diagramindex[", j , "] = ", checkdiagram); print();\
print( "the norm of the Weyl vector is = ", ipd( wln, wln)[1]); print(); print(" the weyl vector is = ",wln); print(); print(); print(); print()  );
print(" the ratio of trials with correct diagram is = ", correctcount/n);
print(" singularcount =", singularcount);
return(wln);
}


\\find the diagram for the reflection group G5

\\for(j=1,3,wy=ran(3,441);for(i=1,50,wy=nuweyl(g5,wy));sg=sim(g5,wy); sg2=[sg[1],sg[2]];print(checkipd(sg2),"  ",n2d(wy),"  ", sg2,"  ",n2d(wy-nuweyl(g5,wy)));print() );

\\the following command produces the right diagram for G26 !
\\ r26=gr26;for(j=1,20,wy=ran(4,2429);for(i=1,100, wy=nuweyl(r26,wy));print(n2d(wy-nuweyl(r26,wy))); s26=sim(r26,wy,.123);for(k=1,4,print(checkipd(s26)[k,])) );



\\444444444444444444444444444444444444444444444444444444444444444444444444
\\44444444                                                        44444444
\\44444444 Working over number rings other than Z[i] and Z[omega] 44444444
\\44444444                                                        44444444
\\444444444444444444444444444444444444444444444444444444444444444444444444



\\the arithmetic geometric and harmonic progression (for testing weyl vectors)/

ap(r,b,d)=
{
local(i,app); app = vector(r); for(i=1,r, app[i] = b + (i-1)*d  ); return(app);
}

gp(r,b,d)=
{
local(i,gpp); gpp = vector(r); for(i=1,r, gpp[i] = b*d^(i-1)  ); return(gpp);
}

hp(r,b,d)=
{
local(i,hpp); hpp = vector(r); for(i=1,r, hpp[i] = (b + (i-1)*d)^(-1)  ); return(hpp);
}

ee(x) = exp(2*Pi*I*x);

zet(m) = ee(1/m);

n2da(x) = conj(x)~ * x; 

\\the input e is a set of column vectos in C^n. The output is the matrix of inner product of these vectors.

checkipda(e)=
{
local(i,j,ll,gram);
ll = matsize(e)[2];
gram = matrix(ll,ll);
for( i = 1, ll, for ( j = 1, ll, gram[i,j] = conj(e[i])~*e[j]  ));
return(gram);
}


\\suppose g is the gram matrix g(i,j) = <r_i,r_j>. then mag(g,k,m) is the matrix of the order m reflection in the vector r_k 
\\where r_k is represented by the kth unit column vector and the matrix acts by left multiplication.
\\this is useful when we have to compute in the reflection group only given the gram matrix of the roots, but not
\\given the co-ordinates of the roots themselves.

mag(g,k,m)=
{
local(l,magg,xii,j);
l = matsize(g)[1]; xii = exp(2*Pi*I/m); magg = matid(l);
for( j = 1, l, magg[k,j] = - ( 1 - xii) * g[k,j]/g[k,k] ); magg[k,k] = xii; 
return(magg);
}

\\ image of the element a under the zeta(m)--reflection in x.

phida(x,m,a)= a - (1 - zet(m)) * ((conj(x)~ * a)/(conj(x)~ * x)) * x;

\\ the matrix of the above reflection when the vectors are represented as columns and the matrix acts on the left.

mada(x,m)= 
{
local(j,matr,n,unitvec);
n=matsize(x)[1];
matr=matrix(n,n);
for (j=1,n,  unitvec=vector(matsize(x)[1])~;  unitvec[j] = 1; matr[,j]=phida(x,m,unitvec) );
return(matr);
}

\\if x is a sequence of roots, pr(x) gives the product of the matrix of reflections of order m in them
\\i.e. the matrix of phida_x1 o phida_x2 o .... o phida_xn

\\the order 2 reflection in a root

phida2(x,a)= a - 2 * ((conj(x)~ * a)/(conj(x)~ * x)) * x;

\\ the matrix of the above reflection when the vectors are represented as columns and the matrix acts on the left.

mada2(x)= 
{
local(j,matr,n,unitvec);
n=matsize(x)[1];
matr=matrix(n,n);
for (j=1,n,  unitvec=vector(matsize(x)[1])~;  unitvec[j] = 1; matr[,j]=phida2(x,unitvec) );
return(matr);
}

prda2(x)=
{
local(lm,prr);
lm=matsize(x)[2];
prr= mada2(x[1]);for(i=2,lm, prr = prr*mada2(x[i]) );
return(prr);
}

\\calculate the order of a finite order matrix

orda(x)=
{
local(n, y, mid); 
mid=matid(matsize(x)[1]); n=1;y=x; while(norml2(y - mid) >10^(-20),n++;y=x*y;if(divrem(n,50)[2]==0,print(n)));return(n);
}

\\check if two column vectors a & b are proportional. rertuns the maximum value of |a[i]*b[j]- a[j]*b[i]|

checkprop(a,b)=
{
local(i,j,l,maxd);
l = matsize(a)[1]; maxd=0;
for(i = 1, l, for(j=i+1,l, if( abs(a[i]*b[j] - a[j]*b[i]) > maxd, maxd = abs(a[i]*b[j] - a[j]*b[i])) ));
return(maxd);
}

\\generate a random vector of length k starting with a random seed.

rana(k,seed)=
{
local(i,ranaa); ranaa = vector(k)~;
for( i = 1, k, ranaa[i] = random(2*seed)/seed - 1 );
return(ranaa);
}

\\the iterating function alpha 
\\input: rg = the list of projective roots (each projective root is a column vector), rgo[i]= th order of reflections in rg[i], 
\\ w = argument of the function

alpha(w,rg,rgo)=
{
local(alp,j,ipp);alp=0;for(j=1,matsize(rg)[2], ipp=conj(rg[j])~ * w; alp=alp + rgo[j]^(-2)*abs(n2da(rg[j]))^(-1/2)*(ipp/abs(ipp))* rg[j] );return(alp);
}

\\the function used to order the projective roots.

orderingfn(w,rg) = abs(conj(rg)~* w)/sqrt(abs(n2da(rg)));


\\find a simple systeim for the projective root system rm using the starting vector [we], the number of iteration has to be specified.
\\Also we have to input r, the number of projective roots closest to the weyl vector to be returned.
\\rm has to be two component vector. the first component rg is the list of the projective roots. 
\\the second component is the list of the orders of reflection in these roots.
\\returns the projective roots that are 


simka(rm,we,iter,r)=
{
local(rg,rgo,nn,l,j,seq, simr, ipsimr);
rg= rm[1]; rgo = rm[2]; nn=matsize(rg)[2]; l = matsize(rg[1])[1]; 
print("_____________________________________________________________________________________________________");
print(); print(" number of projective roots = ",nn); print();
print("starting random vector= ", we); print();
for( j =1, iter, we = alpha( we, rg, rgo) ); 
print("the weyl vector is = ",we); print();
print("level of convergence of iteration after ", iter," steps= ",norml2( we - alpha( we, rg,rgo))); print(); 
seq = vector(nn); for( j = 1, nn, seq[j]=orderingfn(we,rg[j])     ); seq=vecsort(seq, ,1);
simr = vector(r); for(j=1,r,simr[j] = rg[seq[j]] );
print("norm of the weyl vector is= ",norml2(we));
print();
ipsimr = checkipda(simr); for( i = 1, r, print( ipsimr[i,] )); print();
return( ipsimr );
}

\\the code below is to calculate 10 simple systems for G31 and to check how many times we get the right diagram.

\\cou=0;for(j=1,10,ski=simka(r31,rana(4,4291),600,5);vl=valency1(ski);if([vecsort(vl[1]),vl[2]]==[[1,1,2,2,2],1],cou++); print("count=",[cou,j]));

\\Given sk[1], ..., sk[cou], the following code checks if the set is closed under reflection. 
\\If not, increase cou by 1, and add a new vector to it, call it sk[cou+1]. print i,j such that reflecting sk[j] in sk[i] gives sk[cou+1].

\\ 
\\for(i=1,cou,for(j=1,cou,tsk=phida2(sk[i],sk[j]);fg=1;for(k=1,cou,for(l=1,4,if(tsk==sk[k]*I^l,fg=0)));if(fg==1,cou++;sk[cou]=tsk;print([i,j]);break(3) )));


\\******************************************************************************************
\\******************************************************************************************
\\*********************************A list of reflection groups******************************
\\******************************************************************************************
\\******************************************************************************************


\\generate the projective roots, together with the order of reflections in these.
\\the programs below return a two component vector. the first component is a list of projective roots.
\\the second component is the list of the order of reflections in these roots.

\\***********A lattice with reflection group G(de,e,r)  ***********************************

\\If d = 1, then all the roots have norm 2, otherwise there are norm 1 roots too. The indicator srt takes care of this.

rG(de,e,r)=
{
local(rg,j,k,cou,d,srt,rgo);
d=de/e;
srt=1; if(d ==1 , srt=0);
rg = vector(r*srt + de*r * (r-1)/2); rgo = vector(r*srt + de*r * (r-1)/2);
cou=0; if( srt == 1, for(j=1,r, cou++; rg[cou] = vector(r)~; rg[cou][j] = 1 ;  rgo[cou]=d ));
for( j = 1, r, for( k = j+1, r, for( t = 1, de, cou++; rg[cou] = vector(r)~; rg[cou][j] = exp(-2*Pi*I*t/de); rg[cou][k] = -1; rgo[cou]=2 ))); 
return([rg,rgo]);
}


\\***********A lattice over Z[zeta(12)] with reflection group G6  ************************

gr6()=
{
local(rg,rgo);
rg=vector(10); rgo=vector(10);
rg[1] = [2,0]~;                                rg[7] = [ 1, 1 + zet(12)]~;  
rg[2] = zet(3)*(1+I)*[-zet(12),1]~;            rg[8] = [-1, 1 + zet(12)]~;     
rg[3] = zet(3)*(1+I)*[ zet(12),1]~;            rg[9] = [-1-zet(-12), I]~;
rg[4]=  -(1+I)*[ zet(3),1]~;                   rg[10]= [ 1+zet(-12), I]~;
rg[5]=  -(1+I)*[-zet(3),1]~; 
rg[6]=   zet(-12)*[0,2]~;
for(i=1,6, rgo[i] = 2); for( i = 7,10, rgo[i] = 3);
return([rg,rgo]);
}

\\***********A lattice over Z[zeta(12)] with reflection group G6  ************************

gr9()=
{
local(rg,rgo);
rg = vector(18); rgo = vector(18);
rg[1] = [2,0]~; rg[3] = -sqrt(2)*[-1,1]~;
rg[2] = [0,2]~; rg[4] = -sqrt(2)*[ 1,1]~;
rg[5] = [-1-sqrt(2)*I,1]~;  rg[7] = [-1+sqrt(2)*I,1]~;
rg[6] = [ 1+sqrt(2)*I,1]~;  rg[8] = [ 1-sqrt(2)*I,1]~;
rg[9] = [1,-1-sqrt(2)*I]~;  rg[11]= [1,-1+sqrt(2)*I]~;
rg[10]= [1, 1+sqrt(2)*I]~;  rg[12]= [1, 1-sqrt(2)*I]~;
\\
rg[13]=[ 1, 1 + sqrt(2)]~;
rg[14]=[-1, 1 + sqrt(2)]~;
rg[15]=[ 1 + sqrt(2), 1]~;
rg[16]=[ 1 + sqrt(2),-1]~;
rg[17]=(1+zet(8))*[I,1]~;
rg[18]=(1+zet(8))*[-I,1]~;
for(i=1,12, rgo[i] = 2); for( i = 13, 18, rgo[i] = 4 );
return([rg,rgo]);
}

\\****projective roots of D4 lattice over Z[i] with reflction group G8 *****************

gr8()=
{
local(rg,rgo);
rg = vector(6); rgo=vector(6);
rg[1]=[1+I,0]~; rg[2]=[0,1+I]~; rg[3]=[1, 1]~; rg[4]=[1, I]~; rg[5]=[1,-1]~; rg[6]=[1,-I]~;
for(i=1,6,rgo[i] = 4);
return([rg,rgo]);
}

\\***projective roots of E8 lattice over Z[i] with reflction group G31 *****************

gr31()=
{
local(i,j,k,rg,rgo);
cou=0;rg=vector(60); rgo = vector(60);
for(i=1,4,cou++; rg[cou]=vector(4)~; rg[cou][i] = 2);                                                                \\4
cou++; rg[cou]=[1,1,1,1]~; cou++; rg[cou]=[1,-1,-1,-1]~;                                                             \\2
for(i=2,4,for(j=1,2, cou++;rg[cou]=[1,(-1)^j,(-1)^j,(-1)^j]~; rg[cou][i]=-(-1)^j  ));                                \\6
for(i=2,4,for(j=1,2, cou++;rg[cou]=[1, I, I, I]~; rg[cou][i]=(-1)^j  ));                                             \\6
for(i=2,4,for(j=1,2, cou++;rg[cou]=[1,-I,-I,-I]~; rg[cou][i]=(-1)^j  ));                                             \\6
for(j=1,2,cou++;rg[cou]=[1,(-1)^j, I,-I]~; cou++;rg[cou]=[1, I,(-1)^j,-I]~; cou++;rg[cou]=[1, I,-I,(-1)^j]~;\
cou++;rg[cou]=[1,(-1)^j,-I, I]~; cou++;rg[cou]=[1,-I,(-1)^j, I]~; cou++;rg[cou]=[1,-I, I,(-1)^j]~ );                 \\12
for(i=1,4,for(j=i+1,4,for(k=1,4,cou++; rg[cou]=vector(4)~; rg[cou][i] = 1+I; rg[cou][j] = (1+I)*I^k      )));        \\24
for(i=1,60, rgo[i] = 2);
return([rg,rgo]);
}

\\ five roots whose reflections generate G31. These were used while trying to calculate a presentation via deflation.

\\[0, 1 + I, -1 - I, 0]~[0, 0, 1 + I, 1 + I]~[1, 1, -I, -I]~[1, I, -I, 1]~[1, -1, 1, 1]~


\\a couple fixed point for alpha on the root system of G31, the nearest 5 roots generate G31. 
\\(these two weyl vectors have the same norm and the five closest mirrors give the same diagram.)
\\[-0.2117010544367148638425277545 - 5.910387625145599791457966197*I, 13.29879980815553035740468713 + 3.708294334055234902549231112*I, 18.59564199754155073632956414 - 0.6660674846533818796417315325*I, -12.99940930580244331099516586 + 4.650256004307780906595591115*I]~
\\[-9.308359631921902527848006695 - 7.133076409094987976745668103*I, 3.855648802988283043012298563 + 2.954616990090719657152465455*I, 22.45941930672612022021204367 - 2.971492142807650976558311516*I, 9.302996079870094957154826764 - 1.230832346036018797887651855*I]~
\\a fixed point for alpha for G31 such that nearest 5 roots do not generate G31. The value of S is lower here.
\\[1.840082709253807332055378089 + 5.04870979 E-29*I, -14.95334225766278600622077818 - 1.262177448 E-29*I, 12.35107233434782302888483231 + 3.786532346 E-29*I,-19.30713671433962102001691417 - 1.893266173 E-29*I]~


\\***projective roots of a lattice over Z[i] with reflction group G29 *****************

\\the roots of G29 are a subset of the roots of G31. They generate the same lattice E8.

gr29()=
{
local(r31,cou,rg,i); r31 = gr31[1];rg = vector(40); rgo=vector(40); cou = 0;
for(i =  1,  4, cou++; rg[cou] = r31[i] );
for(i = 13, 36, cou++; rg[cou] = r31[i] );
for(i =  1, 12, cou++; rg[cou] = r31[36+2*i] );
for(i=1,40, rgo[i] = 2);
return([rg,rgo]); 
}

\\a couple fixed point for alpha on the root system of G31, the nearest 5 roots generate G31. 
\\(these two weyl vectors have the same norm and the five closest mirrors give the same diagram.)
\\[6.241106617310397943250779365 + 2.943557280614459357138164380*I, -15.15388335281785986606034507 + 1.030646213016061707943684032*I, 0.2927849907313463133307530012 + 4.304900693337637546366781507*I, 5.785083302957902146059677030 - 3.761482746302544819741485565*I]~

\\reflection in the four roots below generate the reflection group G29
\\sg = vector(40); sg[1] = [2,0,0,0]~; sg[2] = [1+I, 1+I, 0, 0]~; sg[3] = [0, 1 + I,-1 - I, 0]~; sg[4] = [1,-I,1,-I]~;


\\555555555555555555555555555555555555555555555555555555555555555555555555
\\55555555                                                        55555555
\\55555555                                                        55555555
\\55555555                                                        55555555
\\555555555555555555555555555555555555555555555555555555555555555555555555



nuweylvar(x, w)= 
{
local(i,l,wey,ipp,nor,fac,odr); 
l = matsize(x[1])[1];  wey = matrix(l,2);
for( i = 1, matsize(x)[2], odr=ord(mad(x[i])); ipp=ipd( x[i], w); nor = odr^2*sqrt(  n2(ipp)[1] * ipd(x[i],x[i])[1]  );\
fac = divid( ipp, [nor,0]);  wey = wey + sdot( fac, x[i]) );
return(wey);
}


simkvar(r, n)=
{
local(l,j,wl,wln,buf,cou,simm, ipdsim, checkdiagram, epsilon, correctcount);
l = matsize(r[1])[1];  correctcount = 0;
for( j = 1, n, wl = ran( l,19246); wln=nuweylvar(r,wl); cou=0; epsilon= ipd(vad( wl, -wln), vad(wl, -wln))[1]/ipd(wl,wl)[1]; \
until( cou > 100 || epsilon < .00000001, epsilon= ipd(vad( wl, -wln), vad(wl, -wln))[1]/ipd(wl,wl)[1];\
wl = wln; wln = nuweylvar( r, wl );  cou++ );print("number of iteration = ", cou , ", epsilon =  ",epsilon); print(); \
simm = sim( r,wln);  ipdsim = checkipd(simm); \
print( " The gram matrix of the inner product between the simple roots are given below "); print(); \
for( j = 1, l, print(ipdsim[j,]) ); print();\
if(l==6, checkdiagram = checkk12d(checkipd(simm))                           );\
if(l==4, checkdiagram = (  vecsort(valency(checkipd(simm))) == [1,1,2,2])   );\
correctcount=correctcount+checkdiagram; print( "diagramindex[", j , "] = ", checkdiagram); print();\
print( "the norm of the Weyl vector is = ", ipd( wln, wln)[1]); print(); print(" the weyl vector is = ",wln); print(); print(); print(); print()  );
print(" the ratio of trials with correct diagram is = ", correctcount/n);
return();
}


alphavar(w,rg,rgo)=
{
local(alp,j,ipp);alp=0; for(j=1,matsize(rg)[2], ipp = conj(rg[j])~ * w ; alp = alp  + (ipp/abs(ipp)) * rg[j]/(rgo[j]^(.5)*sqrt(n2da(rg[j]))  ) );return(alp);
}





\\return the valency vector of a gram matrix.
\\input:  a gram matrix of size (k,k) . 
\\output: a two component vectror. 
\\The first component is a vector [v1, ... , vk] where the j th row conains v_j entries off the diagonal of norm 1. 
\\the second component is 1 if all the entries of the gram matrix is non-zero. Otherwise the second component is 0.

valency1(dip)=
{
local(val, i, j ,k, cou, fg);
k = matsize(dip)[1]; fg=1;
val = vector(k);
for( i = 1, k, val[i] =0; for( j = 1, k, if(dip[i,j] == 0, fg=0); if( (dip[i,j]/2)^4 == 1, val[i]++ )) );
print(" the valencies of the vertices in the diagram are: ",  val); print();
return([val,fg]);
}


simkavar(rm,we,iter,r)=
{
local(rg,rgo,nn,l,j,seq, simr, ipsimr);
rg= rm[1]; rgo = rm[2]; nn=matsize(rg)[2]; l = matsize(rg[1])[1]; 
print("_____________________________________________________________________________________________________");
print(); print("number of projective roots = ", nn); print();
print("starting random vector= ", we); print();
for( j =1, iter, we = alphavar( we, rg, rgo) ); 
print(" the Weyl vector is = ", we); print();
print("level of convergence of iteration after ", iter," steps= ",norml2( we - alphavar( we, rg,rgo))); print(); 
seq = vector(nn); for( j = 1, nn, seq[j]=orderingfn(we,rg[j])     ); seq=vecsort(seq, ,1);
simr = vector(r); for(j=1,r,simr[j] = rg[seq[j]] );
print("norm of the weyl vector is= ",norml2(we));
print();
ipsimr = checkipda(simr); for( i = 1, r, print( ipsimr[i,] )); print();
return( ipsimr );
}


\\**** Generating roots for the exceptional groups and  and relations between them ******

\\deflate a triangle (uses order 3 reflections phid and corresponding matrix mad)

dfei(x,y,z)={
local();
return(  prd([x,y,z,y])==prd([y,z,y,x])  );
}

\\check that the deflation relations for triangles and rectangles hold for G34 simple system. return the gram matrix of the simple system.
\\the simple system was obtained by running simk, modifying them by units to get the inner products equal to omega or omega^2,
\\ and then extending with sk12[7]. The Weyl vector is given below:
\\[72.26206777522818184439900094, -111.1708622481712411958990000; 169.2116163028080392408813988, -76.13361925217387438184393212;\
\\      -5.343435244496078676988554141, 157.3087827610078576343031621; -106.7540760457478428366891011, 47.84610068565485662634808191;\
\\     -155.0857800763446379532700784, 69.93324003225768525937687522; -36.13757479856714877693688810, -24.89698893127501116347622583];

checkrelG34()=
{
local(sk12);
sk12 = vector(7);
sk12[1] = [2/3, 1/3; -1/3, -2/3; -1/3, 1/3; -1/3, -2/3; 2/3, 1/3; -1/3, 1/3];
sk12[2] = -[1/3, -1/3; 1/3, -1/3; 1/3, -1/3; -2/3, -1/3; 1/3, -1/3; 1/3, 2/3];
sk12[3] = sdot(-W,[-1/3, -2/3; -1/3, 1/3; -1/3, -2/3; -1/3, 1/3; -1/3, 1/3; -1/3, -2/3]);
sk12[4] = sdot(W,[-1/3, 1/3; -1/3, 1/3; -1/3, 1/3; -1/3, 1/3; -1/3, 1/3; -1/3, 1/3]);
sk12[5] = sdot(-W,[-1/3, 1/3; 2/3, 1/3; -1/3, 1/3; -1/3, 1/3; 2/3, 1/3; 2/3, 1/3]);
sk12[6] = sdot(-w,[-1/3, -2/3; 2/3, 1/3; 2/3, 1/3; 2/3, 1/3; -1/3, -2/3; -1/3, -2/3]);
sk12[7] = -(sk12[1] + sk12[2] + sk12[3]+ sk12[4] + sk12[5] + sk12[6] );
print("checking relations"); print();
print("check deflation 1-2-3:  ",dfei(sk12[1],sk12[2],sk12[3]));
print("check deflation 2-3-4:  ",dfei(sk12[2],sk12[3],sk12[4]));
print("check deflation 3-4-5:  ",dfei(sk12[3],sk12[4],sk12[5]));
print("check deflation 4-5-6:  ",dfei(sk12[4],sk12[5],sk12[6]));
print("order of 1.3.5.6.5.3 is equal to: ", ord( prd([sk12[1],sk12[3],sk12[5],sk12[6],sk12[5],sk12[3]]) ) ); 
print();
return(checkipd(sk12));
}

checkrelG31()=
{
local(s31);
s31 = vector(5);
s31[1]=[1, 1, -I, -I]~;
s31[2]=-[0, 0, 1 + I, 1 + I]~;
s31[3]=[0, 1 + I, -1 - I, 0]~;
s31[4]=-I*[1, I, -I, 1]~;
s31[5]=-I*[1, -1, 1, 1]~;
print("checking relations"); print();
print( " the order of phi_i * phi_j "); print();
print("    i    j    order"); print("    ---------------");
for( i = 1, 5, for( j = i+1, 5, print("    ",i,"    ",j,"    ",orda(prda2([s31[i],s31[j]])) )    ));
print();
print( "order of 2.4.1.4 is equal to: ", orda(prda2([  s31[2],s31[4],s31[1],s31[4]   ]))  );print();
print( "order of 1.5.2.5 is equal to: ", orda(prda2([  s31[1],s31[5],s31[2],s31[5]   ]))  );print();
print( "order of 1.3.4.3 is equal to: ", orda(prda2([  s31[1],s31[3],s31[4],s31[3]   ]))  );print();
print( "order of 2.4.5.4 is equal to: ", orda(prda2([  s31[2],s31[4],s31[5],s31[4]   ]))  );print();
print( "check if 3.4.5==4.5.3 :       ", prda2([  s31[3],s31[4],s31[5]   ]) ==  prda2([  s31[4],s31[5],s31[3]   ])  );print();
print( "order of 1.2.3.2 is equal to: ", orda(prda2([  s31[1],s31[2],s31[3],s31[2]   ]))  );print();
return(checkipda(s31));
}


\\Returns the gram matrix of Circ_k,nu

circ(k,nu)=
{
local(gm, i);
gm= Mid(k);
for( i = 1, k, gm[i,i] = [3,0] );
for( i = 1, k-1, gm[i,i+1] = -[2,1] ; gm[i+1,i] = -bar([2,1]) );
gm[1,k] = -dot(uni[nu], [2,1]); gm[k,1] = bar(gm[1,k]);
return(gm);
}

\\calculate x^* Gm * y, where x, y are vectors with entries in Z[omega] and Gm is a matrix with entries in Z[omega].

ipcirc(x,Gm,y)=
{
local(ii,jj,i,j,ipp);
ii= matsize(Gm)[1]; jj = matsize(Gm)[2];
ipp = [0,0];
for( i = 1, ii, for( j = 1, jj, ipp = ipp + dot( bar( x[i] ), dot( Gm[i,j] , y[j] ))   ));
return(ipp);
} 


\\checks that certain vectors have zero norm in Circ_k, nu, for k = 3, 4, 5.

checkcircnull(k)=
{
local(nvec,nvecknu,i,j);
nvec = vector(6);
nvec[1] = [[1,0],   -w,    W,     W,     W];
nvec[2] = [[1,0],   -w,    W,[-1,0],[-1,0]];
nvec[3] = [[1,0],   -w,    W,[-1,0],     w];
nvec[4] = [[1,0],   -w,    W,[-1,0],     w];
nvec[5] = [[1,0],[1,0],[1,0], [1,0], [1,0]];
nvec[6] = [[1,0],   -w,   -w,    -w,    -w];
for( i = 1, 6, nvecknu = vector(k); for( j = 1, k, nvecknu[j] = nvec[i][j]); print( i,"  ", ipcirc(nvecknu, circ(k,i), nvecknu), "  "nvecknu); print() );
return(); 
}