\\filename : inp.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

\\there are four major parts: the first three parts are concerned with the
\\calculation in the reflection group of 3E8 + H defined over eisenstein integers.
\\the fourth part is about 

\\1. The 26 simple roots in 3E8+H and some vectors fixed by diagram automorphisms.

\\2. Inner products, reflections and their matrices.

\\3. Generate more reflections by conjugating the simple reflections. Check that
\\   some reflections are in the group generated by a given set of reflections.

\\4. Checking some lattices are indefinite for the article: the diagrams for
\\   complex reflection groups 

h = 30;

\\111111111111111111111111111111111111111111111111111111111111111111111111
\\11111111                                                        11111111
\\11111111  The 26 simple roots in 3E8+H and some vectors fixed   11111111
\\11111111  by diagram automorphisms                              11111111
\\111111111111111111111111111111111111111111111111111111111111111111111111


r=vector(26);

r[ 1]=a = [ 0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0; 0, 0; 0, 0;  1, 0;-1,-1];

r[ 2]=c1= [-1, 1; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0];
r[ 3]=c2= [ 0, 0; 0, 0; 0, 0; 0, 0; -1, 1; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0];
r[ 4]=c3= [ 0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0; 0, 0; 0, 0; -1, 1; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0];

r[ 5]=e1= [ 0, 0;-1, 1; 0, 0; 0, 0;  0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0];
r[ 6]=e2= [ 0, 0; 0, 0; 0, 0; 0, 0;  0, 0;-1, 1; 0, 0; 0, 0;  0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0];
r[ 7]=e3= [ 0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0; 0, 0; 0, 0;  0, 0;-1, 1; 0, 0; 0, 0;  0, 0; 0, 0];

r[ 8]=a1= [ 0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0; 0, 0;-1, 1;  0, 0; 0, 0; 0, 0;-1, 1;  0, 1;-1,-1];
r[ 9]=a2= [ 0, 0; 0, 0; 0, 0;-1, 1;  0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0; 0, 0;-1, 1;  0, 1;-1,-1];
r[10]=a3= [ 0, 0; 0, 0; 0, 0;-1, 1;  0, 0; 0, 0; 0, 0;-1, 1;  0, 0; 0, 0; 0, 0; 0, 0;  0, 1;-1,-1];

r[11]=g1= [ 0, 0; 0, 0; 0, 0;-1, 1;  0, 0; 0, 0; 1,-1;-1, 1;  0, 0; 0, 0; 1,-1;-1, 1;  0, 2;-2,-2];
r[12]=g2= [ 0, 0; 0, 0; 1,-1;-1, 1;  0, 0; 0, 0; 0, 0;-1, 1;  0, 0; 0, 0; 1,-1;-1, 1;  0, 2;-2,-2];
r[13]=g3= [ 0, 0; 0, 0; 1,-1;-1, 1;  0, 0; 0, 0; 1,-1;-1, 1;  0, 0; 0, 0; 0, 0;-1, 1;  0, 2;-2,-2];

r[14]=f = [ 0, 0; 1, 0; 1, 0;-2, 0;  0, 0; 1, 0; 1, 0;-2, 0;  0, 0; 1, 0; 1, 0;-2, 0; -2, 1;-1,-2];

r[15]=f1= [ 0, 0; 1, 0; 1, 0; 1, 0;  0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0];
r[16]=f2= [ 0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 1, 0; 1, 0; 1, 0;  0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0];
r[17]=f3= [ 0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 1, 0; 1, 0; 1, 0;  0, 0; 0, 0];

r[18]=b1= [ 1, 0; 0, 0; 1, 0;-1, 0;  0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0; 0, 0; 0, 0; -1, 0; 0, 0]; 
r[19]=b2= [ 0, 0; 0, 0; 0, 0; 0, 0;  1, 0; 0, 0; 1, 0;-1, 0;  0, 0; 0, 0; 0, 0; 0, 0; -1, 0; 0, 0];
r[20]=b3= [ 0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0; 0, 0; 0, 0;  1, 0; 0, 0; 1, 0;-1, 0; -1, 0; 0, 0];

r[21]=z1= [ 0, 0; 1, 0; 1, 0;-2, 0;  1, 0; 0, 0; 1, 0;-1, 0;  1, 0; 0, 0; 1, 0;-1, 0; -1, 1;-1,-2];
r[22]=z2= [ 1, 0; 0, 0; 1, 0;-1, 0;  0, 0; 1, 0; 1, 0;-2, 0;  1, 0; 0, 0; 1, 0;-1, 0; -1, 1;-1,-2];
r[23]=z3= [ 1, 0; 0, 0; 1, 0;-1, 0;  1, 0; 0, 0; 1, 0;-1, 0;  0, 0; 1, 0; 1, 0;-2, 0; -1, 1;-1,-2];

r[24]=d1= [ 1, 0; 1, 0;-1, 0; 0, 0;  0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0];
r[25]=d2= [ 0, 0; 0, 0; 0, 0; 0, 0;  1, 0; 1, 0;-1, 0; 0, 0;  0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0];
r[26]=d3= [ 0, 0; 0, 0; 0, 0; 0, 0;  0, 0; 0, 0; 0, 0; 0, 0;  1, 0; 1, 0;-1, 0; 0, 0;  0, 0; 0, 0];

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

\\ the four special fixed vectors under L3(F3)

sp =[-1, 1; -1, 1; 2, -2; -5, 5; -1, 1; -1, 1; 2, -2; -5, 5; -1, 1; -1, 1; 2, -2; -5, 5; 1, 9; -10, -10];
sl=[4, 0; 4, 0; 5, 0; -6, 0; 4, 0; 4, 0; 5, 0; -6, 0; 4, 0; 4, 0; 5, 0; -6, 0; -8, 4; -4, -8];

wp=[0, 0; 0, 0; 1, -1; -2, 2; 0, 0; 0, 0; 1, -1; -2, 2; 0, 0; 0, 0; 1, -1; -2, 2; 0, 4; -4, -4];
wl=[1, 0; 1, 0; 2, 0; -3, 0; 1, 0; 1, 0; 2, 0; -3, 0; 1, 0; 1, 0; 2, 0; -3, 0; -3, 2; -2, -4];

\\222222222222222222222222222222222222222222222222222222222222222222222222
\\22222222                                                        22222222
\\22222222  Inner products, reflections and their matrices        22222222
\\22222222                                                        22222222
\\222222222222222222222222222222222222222222222222222222222222222222222222


\\negative defininte standard hermitian form on an eisenstein lattice

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

\\inner product of two vectors in the eisenstein lattice 3E_8+H.

ip(x,y)=
{
local(ipp, i);
ipp = 0; 
for( i= 1, 12, ipp= ipp - dot(bar(x[i,]),y[i,]) );
ipp = ipp + dot( bar(dot(x[13,],[1,2])) , y[14,]) + dot( dot (bar(x[14,]),[1,2]),y[13,]) ;
return(ipp);
}

\\given a set of vectors e, find the inner product matrix (using ipd and ip respectively)

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);
}

checkip(e)=
{
local(i,j,l,gram);
l = matsize(e)[2];
gram = matrix(l,l);
for( i = 1, l, for ( j = 1, l, gram[i,j] = ip(e[i],e[j])));
return(gram);
} 

\\ image of a under trireflection in a vector x of norm -3. 

phi(x,a)= a - sdot( dot([1, -1],ip(x,a))/(-3), x );

\\calculate the matrix of trireflection in a vector x of norm -3.

ma(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=phi(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

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

\\x is a collection of vectors. calculate the determinant of the gram matrix of the vectors.

gm(x)=
{
local(i,j,l);
l = matsize(x)[2];
gmm = matrix(l,l);
for( i = 1, l, for( j = 1, l, gmm[i,j] = ip(x[i],x[j]) ));
return(mdet(gmm));
}

X=[x1,y1; x2,y2; x3,y3; x4,y4; x5,y5; x6,y6; x7,y7; x8,y8; x9,y9; x10,y10; x11,y11; x12,y12; x13,y13; x14,y14];

U=[u1,v1; u2,v2; u3,v3; u4,v4; u5,v5; u6,v6; u7,v7; u8,v8; u9,v9; u10,v10; u11,v11; u12,v12; u13,v13; u14,v14];


\\given an eisenstein vector returns the matrix that corresponds to taking standard inner product with the vector in integer coords

multipd(x)=
{
local(i,j,l,mul, mulv, testv);
l = matsize(x)[1];
mul = mulv = vector(2);
mul[1] = mul[2] = matrix(l,2);
for( i = 1, l, for( j = 1, 2, testv = matrix(l,2);\
 testv[i,j] = 1; mul[1][i,j] = ipd(x,testv)[1]; mul[2][i,j] = ipd(x,testv)[2] ));
 mulv[1] = mtov(mul[1]);
 mulv[2] = mtov(mul[2]);
return(mulv);
}

\\given a set of vectors e from the lattice 3E_8+H finds a set of linearly independent vectors in its perp, 
\\whose simple multiples are in the lattice. uses the basis BB for 3E_8+H defined inside the program 
\\the 1st 12 vectors are the standard basis of 3 copies of E_8, consisting of the simple roots, i.e f_i, e-i,d_i,c_i

eperp(e)=
{
local(ht, l,i,j,Mb, cfeperpz,cfeperp,eperpp,Mbz, Mbj,BB);
ht = matsize(e)[2];
l = 14;
Mb = matrix( l, ht ); 
BB = vector(14);
for( i = 1, 14, BB[i] = matrix(14,2));
BB[1] = f1; BB[2] = e1; BB[3] = d1; BB[4] = c1; 
BB[5] = f2; BB[6] = e2; BB[7] = d2; BB[8] = c2; 
BB[9] = f3; BB[10]= e3; BB[11]= d3; BB[12]= c3; 
BB[13][13,1] = BB[14][14,1] = 1;
for( j = 1, ht, for( i = 1, l, Mb[i,j] = ip(BB[i],e[j]) ));
Mbz = matrix(2*l, 2*ht);
 for( j = 1, ht, \
 Mbj = matrix( l,2); for( i = 1, l, Mbj[i,] = Mb[i,j]);\
 Mbz[,2*j-1 ] = multipd(Mbj)[1]; Mbz[,2*j ] = multipd(Mbj)[2] );
cfeperpz = matkerint( Mbz~); 
cfeperp = eperpp = vector(l - ht);
for( j = 1,  l - ht, cfeperp[j] = vtom( cfeperpz[,2*j-1] ) );
\\for( j= 1, 2*(l-ht), print(cfeperpz[,j]) ); print( matrank(cfeperpz) );  
for( j = 1, l - ht, eperpp[j] = matrix(l,2);\
 for( i = 1, l, eperpp[j] = vad( eperpp[j] , sdot( cfeperp[j][i,] , BB[i]) ) ) );   
return( eperpp);
}

\\given a matrix with eisenstein integer entries, returns the integer matrix acting on the underlying Z lattice

intm(x)=
{
local(l1,l2,i,j,zx);
l1 = matsize(x)[1];l2 =matsize(x)[2];
zx = matrix(2*l1, 2*l2);
for( i = 1, l1, for( j = 1, l2, \
 zx[2*i -1, 2*j-1] = x[i,j][1];zx[2*i-1, 2*j] =-x[i,j][2] ;\
 zx[2*i   , 2*j-1] = x[i,j][2];zx[2*i  , 2*j] = x[i,j][1] - x[i,j][2]));
return(zx);
}

gal = matrix(28,28);
for( i = 1, 14, gal[2*i-1, 2*i-1] = 1; gal[2*i - 1, 2*i] = -1; gal[2*i, 2*i] = -1 );

aprime = [0,0;0,0;0,0;0,0;0,0;0,0;0,0;0,0;0,0;0,0;0,0;0,0; 1,0;0,-1];  

R = matrix(14,26); for(i= 1, 14, for( j = 1, 26, R[i,j] = r[j][i,]));
 

\\calculate teh kernel of the matrix x whose entries are eisenstein integeres

mker(x)=
{
local(intx,intkerx,k1,k2,j,kerx);
intx =intm(x); intkerx = matker(intx); 
k1 = matsize(intkerx)[1]/2; k2 = matsize(intkerx)[2]/2;
kerx = matrix(k1, k2);
for( j = 1, k2, tv1 = vtom(intkerx[,2*j -1]) ; tv2 = vtom(intkerx[,2*j]);\
print(sdot(w,tv1) == tv2 ); for(i = 1, k1, kerx[i,j] = tv1[i,])   );
return(kerx);
} 


\\333333333333333333333333333333333333333333333333333333333333333333333333
\\33333333   Generate more reflections by conjugating the simple  33333333
\\33333333   reflections. Check that some reflections are in the  33333333
\\33333333   group generated by a given set of reflections.       33333333
\\333333333333333333333333333333333333333333333333333333333333333333333333



\\ a hash function for the roots of 3E8+H

hf(x)=
{
local(hhf, i);
hff = 1;
for ( i = 2, 4, if ( divrem(x[i,1],2)[2] == 1  || divrem(x[i,2],2)[2] == 1  , hhf = hhf + 2^(12-i)  ) ;\
		if ( divrem(x[i+4,1],2)[2] == 1  || divrem(x[i,2],2)[2] == 1  , hhf = hhf + 2^(9-i) );\
		if ( divrem(x[i+8,1],2)[2] == 1  || divrem(x[i,2],2)[2] == 1  , hhf = hhf + 2^(6-i) );\
);
for ( i = 13,14,if ( divrem(x[i,1],2)[2] == 1 || divrem(x[i,2],2)[2] == 1  , hhf = hhf + 2^(14-i) ));
i=1;									
if( divrem(x[4,1]+x[8,1],2)[2] == 1  || divrem(x[4,1]+x[8,1],2)[2] == 1 , hhf = hhf + 2^11);
if( divrem(x[4,1]+x[3,1],2)[2] == 1  || divrem(x[4,1]+x[3,1],2)[2] == 1 , hhf = hhf + 2^12);
while ( i <= 14 &&  x[i,] == [0,0] , i++);			\\atmost +2^13
hhf=hhf+i;							\\atmost +14	
i=1;
while ( i <= 12 &&  x[13-i,] == [0,0] , i++);	
hhf=hhf+i*14;				    			\\atmost +196		
hhf = hhf + divrem(x[3,1]^2+x[3,2]^2-x[3,1]*x[3,2]        ,3)[2]*3^7\
	  + divrem(x[4,1]^2+x[4,2]^2-x[4,1]*x[4,2]        ,3)[2]*3^6\
	  + divrem(x[8,1]^2+x[8,2]^2-x[8,1]*x[8,2]        ,3)[2]*3^5 \
	  + divrem((x[8,1]+x[4,1])^2  + (x[8,2]+x[4,2])^2 - (x[8,1]+x[4,1])*(x[8,2]+x[4,2])   ,3)[2]*3^4 \
	  + divrem((x[12,1]-x[13,1])^2+(x[12,2]-x[13,2])^2-(x[12,1]-x[13,1])*(x[12,2]-x[13,2]),3)[2]*3^3 \
	  + divrem(x[10,1]^2+x[10,2]^2-x[10,1]*x[10,2]       ,3)[2]*3^2\
	  + divrem(x[13,1]^2+x[13,2]^2-x[13,1]*x[13,2]       ,3)[2]*3\
	  + divrem(x[14,1]^2+x[14,2]^2-x[14,1]*x[14,2],3)[2] ; 		\\atmost +3^8
hhf = divrem(hhf+1,2)[1];
return(hhf);			\\range of hash values : 1 through 2^14-1 +(3^8-1)/2 + 14^2+14
}

\\starts with a list of m vectors r[1], ..,r[m], and generates more by reflecting them in each other iteratively

genro()=
{
local(m,ttime,R,new,i,hh,ind,has, indmax,n,iold,jold,j,jj,ii,iii,k,flag,hashk,ij,newij);
m=26;ttime = 1024000;R=vector(ttime+1);new=matrix(14,2);has = matrix(5850,800);ind = vector(5850); 
for (i = 1 , m, R[i] = r[i];hh= hf(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 = phi(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 = hf(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;write(rolist, " i[",n,"]=",iold,"; j[", n,"]=",jold, "; ro[",n,"]=",newij,";")));
return();
}

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

\\trying a better hash function for a longer list of ro

hfn(x)=
{
local(hhf);
hhf = 1 + x[13,1] + 61*x[13,2] + 3721*x[14,1] +  226981*x[14,2] + x[1 ,1] + 5*x[1,2] + 25*x[5 ,1] + 125*x[5,2]+ 625*x[9,1] + 3125*x[9,2] + x[2 ,1] + 11*x[2,2] + 121*x[6 ,1] + 1331*x[6,2]+ 14641*x[10,1] + 161051*x[10,2] + x[3 ,1] + 13*x[3,2] +169*x[7 ,1] +2197*x[7,2]+ 28561*x[11,1] + 371293*x[11,2] +x[4 ,1] + 17*x[4,2] + 289*x[8 ,1] + 4913*x[8,2]+ 83521*x[12,1] + 1419857*x[12,2];
hhf = divrem(hhf, 19997)[2]+1;
return(hhf);
}

genron()=
{
local(m,ttime,R,new,i,hh,ind,has, indmax,n,iold,jold,j,jj,ii,iii,k,flag,hashk,ij,newij);
m=26;ttime = 2^20;R=vector(ttime+1);new=matrix(14,2);has = matrix(19998,200);ind = vector(19998); 
for (i = 1 , m, R[i] = r[i];hh= hfn(R[i]);ind[hh]++;has[hh,ind[hh]]=i;write(rolistn, " ro[",i,"]=",R[i],";" )); 
indmax=1; i=1; j=2; n = m;
while ( n < ttime, new = phi(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 = hfn(newij);\
while( k <= ind[hh] , if (newij == R[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, "  ",n));\
has[hh,ind[hh]]=n ; R[n]= new;write(rolistn, "ro[",n,"]=",newij,";")));
return();
}

\\begin??????????????????????????????????????????????????????????????????????????????????????????????????????????

MM=[[0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [1, 0], [0, 0], [0, 0], [0, 0], [1, 0], [-1, 0], [0, 1]; [0, 0], [1, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0]; [0, 0], [0, 0], [-1, 0], [2, 0], [-1, 0], [0, 0], [-1, 0], [2, 0], [-1, 0], [0, 0], [-1, 0], [2, 0], [-4, 0], [0, 4]; [1, 0], [0, 0], [0, 0], [-2, 0], [1, 0], [0, 0], [1, 0], [-2, 0], [1, 0], [0, 0], [1, 0], [-2, 0], [4, 0], [0, -4]; [0, 0],[0, 0], [0, 0], [1, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [1, 0], [-1, 0], [0, 1]; [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [1, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0]; [-1, 0], [0, 0], [-1, 0], [2, 0], [0, 0], [0, 0], [-1, 0], [2, 0], [-1, 0], [0, 0], [-1, 0], [2, 0], [-4, 0], [0, 4]; [1, 0], [0, 0], [1, 0], [-2, 0], [1, 0], [0, 0], [0, 0], [-2, 0], [1, 0], [0, 0], [1, 0], [-2, 0], [4, 0], [0, -4]; [0, 0], [0, 0], [0, 0], [1, 0], [0, 0], [0, 0], [0, 0], [1, 0], [0, 0], [0, 0], [0, 0], [0, 0], [-1, 0], [0, 1]; [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0,0], [0, 0], [1, 0], [0, 0], [0, 0], [0, 0], [0, 0]; [-1, 0], [0, 0], [-1, 0], [2, 0], [-1, 0], [0, 0], [-1, 0], [2, 0], [0, 0], [0, 0], [-1, 0], [2, 0], [-4, 0], [0, 4]; [1, 0], [0, 0], [1, 0], [-2, 0], [1, 0], [0, 0], [1, 0], [-2, 0], [1, 0], [0, 0], [0, 0], [-2, 0], [4, 0], [0, -4]; [2/3, -2/3], [0, 0], [2/3, -2/3], [-5/3, 5/3], [2/3, -2/3], [0, 0], [2/3, -2/3], [-5/3, 5/3], [2/3, -2/3], [0, 0], [2/3, -2/3], [-5/3, 5/3], [4, -3], [-3, -6]; [2/3, 4/3], [0, 0], [2/3, 4/3], [-5/3, -10/3], [2/3, 4/3], [0, 0], [2/3, 4/3], [-5/3, -10/3], [2/3, 4/3], [0, 0], [2/3, 4/3], [-5/3, -10/3], [3, 6], [7, 3]];

\\a norm zero vector that has satisfies <r,z>=\theta for r = bi,ci,di,ei,fi; 
\\Z has the form (\rho,\rho,\rho, * , h) , where h = 1 = (2+ \omega^2) + 2 + (2+ \omega^2) + 1 = 6 - 2\omega
\\and \rho is the weyl vector of E8
  
Z=[1, 1; 1, 1; 3, 4; -5, -7; 1, 1; 1, 1; 3, 4; -5, -7; 1, 1; 1, 1; 3, 4; -5, -7; -9, -6; 6, -2];

\\end?????????????????????????????????????????????????????????????????????????????????????????????????????????????


\\show that the 16 reflections of the M666 diagram generates the 26 simple reflections.

YgiveD()=
{
local();
print( phi(   phi( b2       ,         a)     ,    phi(b1       , phi(b1,c1))   )); print(); print(          a3," =a3"      ); print();
print( phi(   phi(phi(b1,c1), phi(a,b2))     ,    phi(phi(a,b1), phi(b3, a))   )); print(); print(sdot( W, z1)," = w^2*z1" ); print();
print( phi(   phi( z1       ,        c2)     ,    phi(b2       , phi( a,b1))   )); print(); print(sdot(-w, g3)," =  -w*g3" ); print();
print( phi(   phi(phi(b1,c1),        z2)     ,    phi(c3       ,        b3 )   )); print(); print(sdot(-w, f )," =  -w*f " ); print();
return();
}  

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

\\check that 14 reflections are enough to generate Aut(L)
\\we check that the clockwise and anticlockwise deflation relation holds and that reflection in f2 can be written
\\as a product of reflections in f1, e1, d1,c1,b1,a,b2,c2,d2,e2 

generatedby14()=
{
local(acl,cl,f1tod2,f1tod3,e2toe1,e3toe1,lhs);
acl= pr([f1,e1,d1,c1,b1,a,b2,c2,d2,e2]);
cl= pr([f2,e2,d2,c2,b2,a,b1,c1,d1,e1]);
print( md( acl , md( ma(f2), minv(acl))) == ma(a3),"  ",md( cl , md( ma(f1), minv(cl))) == ma(a3));
f1tod2=pr([f1,e1,d1,c1,b1,a,b2,c2,d2]);
e2toe1=pr([e2,d2,c2,b2,a,b1,c1,d1,e1]);
lhs=md( ma(e2), md(minv(f1tod2),e2toe1 ));
print( md( lhs, md( ma(f1), minv(lhs))) == ma(f2));
f1tod3=pr([f1,e1,d1,c1,b1,a,b3,c3,d3]);
e3toe1=pr([e3,d3,c3,b3,a,b1,c1,d1,e1]);
lhs=md( ma(e3), md(minv(f1tod3),e3toe1 ));
print( md( lhs, md( ma(f1), minv(lhs))) == ma(f3));
return();
}

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

\\ Let zi = wp + theta*wl be Leech cusp (because its stabilizer contains L3(F3) and looking at the list of autom. gp. of Eisenstin Niemier Lattices.
\\ Calculations to check that the roots (zi + omega^2*l) and (zi + x + l) where x is incident on l belong to the gp. gen. by 26 ref.  
\\this calculation is needed for the alternative verification of Daniel Allcock's theorem in the article: "on the Y555 reflection group"
checkht1roots()=
{
local(pe,i,zi,ss,q1,q2,q3,q1p,q2p);
pe = [2,1];
print(wp == sdot( bar(pe), f) + a + e1 + e2 + e3, "    wp == sdot( bar(pe), f) + a + e1 + e2 + e3  ");
print(wl == sdot( pe , a) + f + b1 + b2 + b3,     "    wp == sdot( pe, f) + f + b1 + b2 + b3  ");
print();
print(ip(wp, f) == pe,"    ip(wp, f) == pe");  print(bar(ip(wl, a)) == pe,"    bar(ip(wl, a)) == pe");  
print();
zi = wp + sdot(t,wl);
print("if ip(W*l,zi) or ip(x+l, zi) not equal to theta print PROBLEM"); 
print();print();
for(i=14,26,if(ip(sdot(W,r[i]),zi)!=t, print("PROBLEM")); for(j=1,13,if(ip(r[j],r[i]) !=0 && ip(r[j]+r[i],zi)!= t,print("PROBLEM"))));
ss = [sdot( -W, a1 + b1),sdot(-w, f2),sdot(-w, f3),sdot(-W, a + b1),sdot(-W, f + e1)];
print(ss == [sdot( -W, a1 + b1),sdot(-w, f2),sdot(-w, f3),sdot(-W, a + b1),sdot(-W, f + e1)],"    ss == [sdot( -W, a1 + b1),sdot(-w, f2),sdot(-w, f3),sdot(-W, a + b1),sdot(-W, f + e1)]" );
print();
print("check the inner product matrix of ss[i] and then then ip(ss[i], zi)   ");
print();
for( i = 1, 5, print(checkip(ss)[i,] ));
print();
print(ip(ss[1],zi), ip(ss[2],zi),ip(ss[3],zi),ip(ss[4],zi),ip(ss[5],zi));
print();
print("check the calculation from q_1 to q_2 ");
print();
q1 = zi + e1 + f; print( q1 == zi + e1 + f ,"    q1 == zi + e1 + f");
print();
q2 = zi+sdot(w,ss[5]+ss[4])+ss[3]+ss[2]-sdot(t,ss[1]); print(q2 == zi+sdot(w,ss[5]+ss[4])+ss[3]+ss[2]-sdot(t,ss[1]),"    q2 == zi+sdot(w,ss[5]+ss[4])+ss[3]+ss[2]-sdot(t,ss[1])");
print();
print(q1 == zi - sdot(w, ss[5]),"    q1 == zi - sdot(w, ss[5])");    
print();
print(phi(ss[5],phi(ss[5],q1))== zi + sdot(w, ss[5]),"    phi(ss[5],phi(ss[5],q1))== zi + sdot(w, ss[5])");
print();
print(phi(ss[4],phi(ss[4],zi + sdot(w, ss[5])))==zi + sdot(w,ss[5]+ss[4]),"    phi(ss[4],phi(ss[4],zi + sdot(w, ss[5])))==zi + sdot(w,ss[5]+ss[4])");
print();
print(phi(ss[2],phi(ss[3],zi + sdot(w, ss[4] + ss[5])))==zi + sdot(w,ss[5]+ss[4]) + ss[3] + ss[2],"    phi(ss[2],phi(ss[3],zi + sdot(w, ss[4] + ss[5])))==zi + sdot(w,ss[5]+ss[4]) + ss[3] + ss[2]");
print();
print(phi(ss[1],zi+sdot(w,ss[5]+ss[4])+ss[3]+ss[2])==q2,"    phi(ss[1],zi+sdot(w,ss[5]+ss[4])+ss[3]+ss[2])==q2");
print();
print(q2 == e2 + e3 + sdot(w, z1 - f) - sdot(W, wl),"    q2 == e2 + e3 + sdot(w, z1 - f) - sdot(W, wl)");
print();
print("check the calculation from q3 to q4");
print();
q3 = phi(e2, phi(e2, phi(e3, phi(e3, q2)))); print(q3 == phi(e2, phi(e2, phi(e3, phi(e3, q2)))),"    q3 == phi(e2, phi(e2, phi(e3, phi(e3, q2))))"); 
print();
print(q3 == sdot(w, z1 - f) - sdot( W, wl),"    q3 == sdot(w, z1 - f) - sdot( W, wl)");
print();
print(q3 == f + sdot( t, z1 + e1) - sdot( W, f1 + d1),"    q3 == f + sdot( t, z1 + e1) - sdot( W, f1 + d1)");
print();
print(phi(z1,q3)==f + sdot(t, e1) - sdot(W, f1 + d1),"    phi(z1,q3)==f + sdot(t, e1) - sdot(W, f1 + d1)");
print();
print(phi(e1, f + sdot(t, e1) - sdot(W, f1 + d1))== f - sdot(W, e1 + f1 + d1),"    phi(e1, f + sdot(t, e1) - sdot(W, f1 + d1))== f - sdot(W, e1 + f1 + d1)");
print();
print(phi(f,phi(f,f - sdot(W, e1 + f1 + d1))) == - sdot(W, e1 + f1 + d1),"    phi(f,phi(f,f - sdot(W, e1 + f1 + d1))) == - sdot(W, e1 + f1 + d1)");           
print();
print();
print();
print();
print("check the root zi + omega^2 * d1");
print();
q1p=zi + sdot(W, d1); print(q1p==zi + sdot(W, d1),"    q1p==zi + sdot(W, d1)");
print();
print(q1p==c1+g2+g3+d1+sdot(t,f+f1+z1)+sdot([1,3],e1),"    q1p==c1+g2+g3+d1+sdot(t,f+f1+z1)+sdot([1,3],e1)");
print();
q2p=d1+sdot(t,f+f1+z1)+sdot([1,3],e1); print( q2p==d1+sdot(t,f+f1+z1)+sdot([1,3],e1),"    q2p==d1+sdot(t,f+f1+z1)+sdot([1,3],e1)");
print();
print(q1p==phi(g3,phi(g2,phi(c1,q2p))),"    q1p==phi(g3,phi(g2,phi(c1,q2p)))");    
print();
print(phi(e1,phi(f,phi(f1,phi(e1,phi(z1,phi(e1,phi(f,phi(f1,q2p))))))))==d1,"    phi(e1,phi(f,phi(f1,phi(e1,phi(z1,phi(e1,phi(f,phi(f1,q2p))))))))==d1");
return();
}


\\444444444444444444444444444444444444444444444444444444444444444444444444
\\44444444                                                        44444444
\\44444444 Checking some lattices are indefinite for the article: 44444444
\\44444444   the diagrams for complex reflection groups           44444444
\\444444444444444444444444444444444444444444444444444444444444444444444444


\\the roots of Eisenstein E_8 upto units.

re8= vector(40);

re8[ 1] = [ 0 , 0 ;  1 , 0 ;  1 , 0 ;  1 , 0 ];re8[10] = [ 1 , 0 ;  0 , 0 ;  1 , 0 ; -1 , 0 ];
re8[ 2] = [ 0 , 0 ;  1 , 0 ;  1 , 0 ;  0 , 1 ];re8[11] = [ 1 , 0 ;  0 , 0 ;  1 , 0 ;  0 ,-1 ];
re8[ 3] = [ 0 , 0 ;  1 , 0 ;  1 , 0 ; -1 ,-1 ];re8[12] = [ 1 , 0 ;  0 , 0 ;  1 , 0 ;  1 , 1 ];
re8[ 4] = [ 0 , 0 ;  1 , 0 ;  0 , 1 ;  1 , 0 ];re8[13] = [ 1 , 0 ;  0 , 0 ;  0 , 1 ; -1 , 0 ];
re8[ 5] = [ 0 , 0 ;  1 , 0 ;  0 , 1 ;  0 , 1 ];re8[14] = [ 1 , 0 ;  0 , 0 ;  0 , 1 ;  0 ,-1 ];
re8[ 6] = [ 0 , 0 ;  1 , 0 ;  0 , 1 ; -1 ,-1 ];re8[15] = [ 1 , 0 ;  0 , 0 ;  0 , 1 ;  1 , 1 ];
re8[ 7] = [ 0 , 0 ;  1 , 0 ; -1 ,-1 ;  1 , 0 ];re8[16] = [ 1 , 0 ;  0 , 0 ; -1 ,-1 ; -1 , 0 ];
re8[ 8] = [ 0 , 0 ;  1 , 0 ; -1 ,-1 ;  0 , 1 ];re8[17] = [ 1 , 0 ;  0 , 0 ; -1 ,-1 ;  0 ,-1 ];
re8[ 9] = [ 0 , 0 ;  1 , 0 ; -1 ,-1 ; -1 ,-1 ];re8[18] = [ 1 , 0 ;  0 , 0 ; -1 ,-1 ;  1 , 1 ];

re8[19] = [ 1 , 0 ;  1 , 0 ; -1 , 0 ;  0 , 0 ];re8[28] = [ 1 , 0 ; -1 , 0 ;  0 , 0 ;  1 , 0 ];
re8[20] = [ 1 , 0 ;  1 , 0 ;  0 ,-1 ;  0 , 0 ];re8[29] = [ 1 , 0 ; -1 , 0 ;  0 , 0 ;  0 , 1 ];
re8[21] = [ 1 , 0 ;  1 , 0 ;  1 , 1 ;  0 , 0 ];re8[30] = [ 1 , 0 ; -1 , 0 ;  0 , 0 ; -1 ,-1 ];
re8[22] = [ 1 , 0 ;  0 , 1 ; -1 , 0 ;  0 , 0 ];re8[31] = [ 1 , 0 ;  0 ,-1 ;  0 , 0 ;  1 , 0 ];
re8[23] = [ 1 , 0 ;  0 , 1 ;  0 ,-1 ;  0 , 0 ];re8[32] = [ 1 , 0 ;  0 ,-1 ;  0 , 0 ;  0 , 1 ];
re8[24] = [ 1 , 0 ;  0 , 1 ;  1 , 1 ;  0 , 0 ];re8[33] = [ 1 , 0 ;  0 ,-1 ;  0 , 0 ; -1 ,-1 ];
re8[25] = [ 1 , 0 ; -1 ,-1 ; -1 , 0 ;  0 , 0 ];re8[34] = [ 1 , 0 ;  1 , 1 ;  0 , 0 ;  1 , 0 ];
re8[26] = [ 1 , 0 ; -1 ,-1 ;  0 ,-1 ;  0 , 0 ];re8[35] = [ 1 , 0 ;  1 , 1 ;  0 , 0 ;  0 , 1 ];
re8[27] = [ 1 , 0 ; -1 ,-1 ;  1 , 1 ;  0 , 0 ];re8[36] = [ 1 , 0 ;  1 , 1 ;  0 , 0 ; -1 ,-1 ];

re8[37] = [ 1 , 2 ;  0 , 0 ;  0 , 0 ;  0 , 0 ];
re8[38] = [ 0 , 0 ;  1 , 2 ;  0 , 0 ;  0 , 0 ];
re8[39] = [ 0 , 0 ;  0 , 0 ;  1 , 2 ;  0 , 0 ];
re8[40] = [ 0 , 0 ;  0 , 0 ;  0 , 0 ;  1 , 2 ];

\\check when a lattice with diagram given by triangle or square or pentagon have positive determinant
\\if is does not have positive determinant find a null vector with one of the coefficients equal to 1.
checkpolygon(n)=
{
local(i,G);
G = matrix(n,n); 
for(i = 1, n, for( j = 1, n, G[i,j] = [0,0] ); G[i,i] = [-3,0]);
for(i = 1, n-1, G[i,i+1] = [2,1] ; G[i+1, i] = bar([2,1]) );
for(i = 1, 6,\
 G[1,n] = dot(uni[i], [2,1]); G[n,1] = bar(G[1,n]);\
 if( (-1)^n*mdet(G)[1] > 0, print(uni[i],"  gives poitive definite ", mdet(G) ));\
 if( (-1)^n*mdet(G)[1]== 0, print(uni[i],"  gives singular ", mdet(G) )   );\
 if( (-1)^n*mdet(G)[1] < 0, print(uni[i],"  gives negative definite ", mdet(G) )));return();
}

\\Given a row vector c and a matrix Hermitian G, with entries in Eisenstein lattice calculate the c*.G.c~

nm(c,G)=
{
local(n,i,nmm);
n = matsize(G)[1]; nmm = 0 ; 
for( i = 1, n, nmm = nmm +  rea(dot( dot(bar(c[i]),c[i]),G[i,i] ))      );
for( i = 1, n, for( j = i+1, n, nmm = nmm + 2 * rea(dot( dot(bar(c[i]),c[j]),G[i,j] ))  ));
return(nmm); 
}


\\if G does not have positive determinant find a vector c with the 1st entry equal to 1 such that c*.G.c~ = 0


null3()=
{
local(i,rg,G,i2,i3,nmm,cf);
rg=3;
G = matrix(rg,rg); 
for(i = 1, rg, for( j = 1, rg, G[i,j] = [0,0] ); G[i,i] = [-3,0]);
for(i = 1, rg-1, G[i,i+1] = [2,1] ; G[i+1, i] = bar([2,1]) );
for(i = 1, 6, G[1,rg] = dot(uni[i], [2,1]); G[rg,1] = bar(G[1,rg]);\
 if( (-1)^rg*mdet(G)[1] > 0, print(uni[i],"  gives positive definite, the det is = ",(-1)^rg*mdet(G)[1]); print());\
 if( (-1)^rg*mdet(G)[1] <= 0, print(uni[i],"  gives non-positive definite, the det is = ",(-1)^rg*mdet(G)[1]); print();\
 for(i2=1,6, for(i3=1,6, cf =[[1,0],uni[i2],uni[i3]]; nmm = nm(cf, G); if(nmm == 0, print(cf); break(rg-1) ) )); print()  )); 
 return(); 
 }
 
 null4()=
{
local(i,rg,G,i2,i3,i4,nmm,cf);
rg=4;
G = matrix(rg,rg); 
for(i = 1, rg, for( j = 1, rg, G[i,j] = [0,0] ); G[i,i] = [-3,0]);
for(i = 1, rg-1, G[i,i+1] = [2,1] ; G[i+1, i] = bar([2,1]) );
for(i = 1, 6, G[1,rg] = dot(uni[i], [2,1]); G[rg,1] = bar(G[1,rg]);\
 if( (-1)^rg*mdet(G)[1] > 0, print(uni[i],"  gives positive definite, the det is = ",(-1)^rg*mdet(G)[1]); print());\
 if( (-1)^rg*mdet(G)[1] <= 0, print(uni[i],"  gives non-positive definite, the det is = ",(-1)^rg*mdet(G)[1]); print();\
 for(i2=1,6,for(i3=1,6,for(i4=1,6, cf =[[1,0],uni[i2],uni[i3],uni[i4]]; nmm = nm(cf, G); if(nmm == 0, print(cf);break(rg-1) );   ))); print()  )); 
 return(); 
 }
 
 null5()=
{
local(i,rg,G,i2,i3,i4,i5,nmm,cf);
rg=5;
G = matrix(rg,rg); 
for(i = 1, rg, for( j = 1, rg, G[i,j] = [0,0] ); G[i,i] = [-3,0]);
for(i = 1, rg-1, G[i,i+1] = [2,1] ; G[i+1, i] = bar([2,1]) );
for(i = 1, 6, G[1,rg] = dot(uni[i], [2,1]); G[rg,1] = bar(G[1,rg]);\
 if( (-1)^rg*mdet(G)[1] > 0, print(uni[i],"  gives positive definite, the det is = ",(-1)^rg*mdet(G)[1]); print());\
 if( (-1)^rg*mdet(G)[1] <= 0, print(uni[i],"  gives non-positive definite, the det is = ",(-1)^rg*mdet(G)[1]); print();\
 for(i2=1,6,for(i3=1,6,for(i4=1,6,for(i5=1,6, cf =[[1,0],uni[i2],uni[i3],uni[i4],uni[i5]]; nmm = nm(cf, G); if(nmm == 0, print(cf); break(rg-1))  )))); print()  )); 
 return(); 
 }
 
 
extriangle1()=
{
local(G, i, j);
G = matrix(5,5); cf = matrix(5,2); 
for(i = 1, 5, for( j = 1, 5, G[i,j] = [0,0] ); G[i,i] = [-3,0]);
G[1,2] = [2,1]; G[1,3] = [2,1]; G[1,4] = [2,1]; G[2,3] = [-2,-1]; G[3,5] = [2,1];
for( i = 1, 5, for( j = 1, i-1, G[i,j] = bar(G[j,i])  ));
print( mv( G, [2,1;0,-1;2,1;1,0;1,0]));
return(G);
}


extriangle2()=
{
local(G, i, j);
G = matrix(5,5); cf = matrix(5,2); 
for(i = 1, 5, for( j = 1, 5, G[i,j] = [0,0] ); G[i,i] = [-3,0]);
G[1,2] = [2,1]; G[1,3] = [2,1]; G[1,4] = [2,1]; G[2,3] = dot(w,[2,1]); G[3,5] = [2,1];
for( i = 1, 5, for( j = 1, i-1, G[i,j] = bar(G[j,i])  ));
print( mv( G,[2,1;1,1;1,-1;1,0;0,-1] ));
return(G);
}


\\the gram matrix for the diagrams  given below

\\                               x5                   c --<-- d means <c,d>=p=[2,1];
\\                        (-1) /    \
\\                 x1 -->-- x2 --<-- x3 -->-- x4   ; <x3, x5>=p ; (*): for <x2,x5> = (-1)p get G51; for <x2,x5>=W.p get G52  

\\ we can run mker(G51) to find the balanced numbering on G51

G51= [\
[-3, 0],[ 1,-1],[ 0, 0],[ 0, 0],[ 0, 0];\
[ 2, 1],[-3, 0],[ 2, 1],[ 0, 0],[-2,-1];\
[ 0, 0],[ 1,-1],[-3, 0],[ 1,-1],[ 2, 1];\
[ 0, 0],[ 0, 0],[ 2, 1],[-3, 0],[ 0, 0];\
[ 0, 0],[-1, 1],[ 1,-1],[ 0, 0],[-3, 0]];

G52= [\
[-3, 0],[ 1,-1],[ 0, 0],[ 0, 0],[ 0, 0];\
[ 2, 1],[-3, 0],[ 2, 1],[ 0, 0],[-1,-2];\
[ 0, 0],[ 1,-1],[-3, 0],[ 1,-1],[ 2, 1];\
[ 0, 0],[ 0, 0],[ 2, 1],[-3, 0],[ 0, 0];\
[ 0, 0],[ 1, 2],[ 1,-1],[ 0, 0],[-3, 0]];



checkbalance()=
{
local(pp,pb); 
pp=[2,1];pb=bar(pp);print();
print(" c --<-- d means <c,d>=p=[2,1]; " ); print();
print("checking the diagram: ");
print("                           x5    ");
print("                   (*)  /    |                                               ");   
print("            x1 -->-- x2 --<-- x3 -->-- x4   ; <x3, x5>=p ; (*): for <x2,x5> = (-1)p "); 
print();
print("at x1:  ", dot([1,0], pb) == dot(-w,pp));   
print("at x2:  ", dot( pp  , pb) ==  [1,0] + pb + dot([-1,0], W)    );   
print("at x3:  ", dot( pb  , pb) ==  dot(-w,pp) - w + W    );  
print("at x4:  ", pb             == pb);  
print("at x5:  ", dot(W,pb)      == dot(dot(-w,bar([-1,0])),pp) + dot(-w,pb)   );   
print();
print("checking the diagram: ");print();
print("                           x5   ");
print("                   (*)  /    |                                               ");   
print("            x1 -->-- x2 --<-- x3 -->-- x4   ; <x3, x5>=p ; (*): <x2,x5> =  W.p "); 
print();
print("at x1:  ", dot([1,0], pb) == dot(-w,pp));   
print("at x2:  ", dot( pp  , pb) ==  [1,0] + pp + dot(W, -W)    );   
print("at x3:  ", dot( pp  , pb) ==  dot(-w,pp) + dot(-w,-W) - W   );  
print("at x4:  ", dot( -W  , pb) ==   pp );  
print("at x5:  ", dot(-W,pb)     == dot(dot( -w, bar(W) ) ,pp) + dot(-w,pp)   );   
return();
}