\\filename : isom.gp

\\*****************************************************C12***********************************************************************

\\generate the C4 Code
C4 = vector(9);
C4[1] = [ 0,0; 0,0; 0,0; 0,0];
C4[2] = [ 1,0; 1,0;-1,0; 0,0];C4[3] = [-1,0;-1,0; 1,0; 0,0];
C4[4] = [ 1,0;-1,0; 0,0; 1,0];C4[5] = [-1,0; 1,0; 0,0;-1,0];
C4[6] = [ 1,0; 0,0; 1,0;-1,0];C4[7] = [-1,0; 0,0;-1,0; 1,0];
C4[8] = [ 0,0; 1,0; 1,0; 1,0];C4[9] = [ 0,0;-1,0;-1,0;-1,0];

\\check if a vector is in E_8, returns 1 if its in E_8, zero otherwise,
checkE8(v)=
{
local(i, flag);
flag = 0; for( i = 1, 9, if( modt(v) == C4[i] , flag = 1));
return(flag);
} 


\\generate the C12 with wt distribution 
\\x^12 + y^12 + z^12 + 22(x^6y^6 + x^6z^6 + y^6z^6) +220( x^6y^3z^3 + x^3y^6z^3 + x^3y^3z^6)
\\basis taken from page 85 of SPLAG
\\ note that is is not the same subspace used in the minimog construction or in wilson.

C12()=
{
local(gen,i,j,i1,i2,i3,i4,i5,i6, word,countt);
gen = vector(6);
for( i = 1, 6, gen[i] = vector(12));
gen[1]= [1, 0, 0, 0, 0, 0,   0, 1, 1, 1, 1, 1]; 
gen[2]= [0, 1, 0, 0, 0, 0,  -1, 0, 1,-1,-1, 1];
gen[3]= [0, 0, 1, 0, 0, 0,  -1, 1, 0, 1,-1,-1];
gen[4]= [0, 0, 0, 1, 0, 0,  -1,-1, 1, 0, 1,-1];
gen[5]= [0, 0, 0, 0, 1, 0,  -1,-1,-1, 1, 0, 1];
gen[6]= [0, 0, 0, 0, 0, 1,  -1, 1,-1,-1, 1, 0];
countt = 0;word = vector(3^6);
for( i1 = -1,1, for( i2 = -1,1, for( i3 = -1,1,for( i4 = -1,1, for( i5 = -1,1, for( i6 = -1,1,\
countt = countt+1 ; \
word[countt] = i1*gen[1]+i2*gen[2]+i3*gen[3]+i4*gen[4]+i5*gen[5]+i6*gen[6];\
for( j = 1, 12, word[countt][j] = divrem(word[countt][j],3)[2]);\
write(outttt, "wordc12[",countt,"]=",word[countt],";")\
))))));
return();
}

\\the C12 coordinates that are used in wilson

shuf(v)=
{
local(sv);
cyc= [7,10,6,8,9,12,3,2,4,11,5];
for( i = 1, 10,sv[cyc[i+1]] = v[ cyc[i]  ]); 
return(sv);
}

Wt( v, eps)=
{
local(l,i,cou);
l = matsize(v)[2];
cou = 0;
for( i = 1, l, if( v[i] == eps, cou = cou+1 ) ) ;
return(cou);
}

checkwt()=
{
local(wt3,wt6,wt9,wt12,i,j);
cou = matrix(13,13);
for( i = 1, 729, \
cou[Wt(wordc12[i],1)+1, Wt(wordc12[i], -1)+1] = cou[Wt(wordc12[i],1)+1, Wt(wordc12[i], -1)+1] +1);
for( i = 1, 13, for( j = 1, 13, if( cou[i,j] != 0, print( " +",i-1," -",j-1, " ", cou[i,j]) )));
return();
}

monomial()=
{
local(monom, cou, i,j);
monom = vector(243);
cou = 0;
for( i = 1, 729, if( wordc12[i][1] == 0, cou = cou+1;\
monom[cou] = matrix(12,2);\
for( j = 1, 12, \
if( wordc12[i][j]== 0,   monom[cou][j,] = [ 1, 0]) ;\
if( wordc12[i][j]== 1,   monom[cou][j,] = [ 0, 1]) ;\
if( wordc12[i][j]==-1,   monom[cou][j,] = [-1,-1]) )));
return(monom);
}

\\****************************************************LEECH*********************************************************************

\\return 1 if the vector v is in the leech lattice and 0 otherwise
\\ the leech lattice vectors are those of the form (m + theta*y_i + 3* z_i) , where
\\ m is 0, 1, or -1, y is in c12 , and the sum of the z_i 's is equal to m mod theta.

checkleech(v)=
{
local(flag, mtv, i,v1,y,c12word,zz,sumz,fff);
flag = 0;
mtv = modt(v) ; m = mtv[1,1];
for( i = 1, 12, if( mtv[i,1] != m , return(flag)));
v1 = sdivid(vad( v, -mtv), t); 
y = modt(v1); c12word = y[,1]; 
for( i = 1, 729, if( c12word == wordc12[i]~, flag = 1 ));
if( flag == 0, return(flag));
z= sdivid( vad( v1 , -y), -t);
sumz = [0,0]; for( i = 1, 12, sumz = sumz + z[i,]);
fff= divid(( sumz - [m,0]), t );
if( frac(fff[1]) != 0 || frac(fff[2]) != 0, flag = 0);
return(flag);  
}

\\generate the norm 2 vectors of Leech lattice upto units

Leech21()=
{
local(i,j,k,word,cou);
cou  = 0;
for( i = 1, 12, for( j = i+1 , 12, for( k = 1, 3, cou = cou+1;\
word = matrix(12,2); word[i,] = [3,0]; word[j,] = dot( uni[2*k], [-3,0]) ;\
write( outttt, "leech2[",cou,"]=",word, ";") )));
return(cou); 
} 
 
Leech22()=
{
local(i,cword, word, k2,k3,k4,k5,k6, rr,  ind, j, sumz ,cou);
cword = vector(22);cou = 0; rr = 0; ind = vector(6);
for( i = 1, 729, if( Wt(wordc12[i] ,1) == 6 && Wt(wordc12[i], -1) == 0, rr = rr+1;cword[rr] = wordc12[i] ));
for( i = 1, 22, rr = 0; \
for( j = 1, 12, if( cword[i][j] != 0, rr = rr+1; ind[rr] = j));for( k2 = 1, 3, for( k3 = 1, 3, for( k4 = 1, 3, for( k5 = 1, 3,for( k6 = 1, 3, \
sumz = uni[2*k2] + uni[2*k3] +uni[2*k4] + uni[2*k5] + uni[2*k6] -[5,0];\
if( divrem(sumz[1],3)[2] == 0 || divrem( sumz[2] ,3)[2] == 0,\
cou = cou+1; word = matrix(12,2); \
word[ind[1],] = (cword[i][ind[1]])*t;                  word[ind[2],] = (cword[i][ind[2]])*dot( t, uni[2*k2]);\
word[ind[3],] = (cword[i][ind[3]])*dot( t, uni[2*k3]); word[ind[4],] = (cword[i][ind[4]])*dot( t, uni[2*k4]);\
word[ind[5],] = (cword[i][ind[5]])*dot( t, uni[2*k5]); word[ind[6],] = (cword[i][ind[6]])*dot( t, uni[2*k6]);\
write( outttt, "leech2[",cou+198,"]=",word, ";") )   )))))  );
return(cou);
}

Leech23()=
{
local(i,cword, word, k2,k3,k4,k5,k6, rr,  ind, j, sumz ,cou);
cword = vector(110);cou = 0; rr = 0; ind = vector(6);
for( i = 1, 729, if( Wt(wordc12[i] ,1) == 3 && Wt(wordc12[i], -1) == 3,\
positn= 1; while(  wordc12[i][positn] == 0, positn = positn+1 ); if(  wordc12[i][positn] == 1,\
 rr = rr+1;cword[rr] = wordc12[i] ))); 
for( i = 1, 110, rr = 0; \
for( j = 1, 12, if( cword[i][j] != 0, rr = rr+1; ind[rr] = j));for( k2 = 1, 3, for( k3 = 1, 3, for( k4 = 1, 3, for( k5 = 1, 3,for( k6 = 1, 3, \
sumz = cword[i][ind[2]]*uni[2*k2] + cword[i][ind[3]]*uni[2*k3] + cword[i][ind[4]]*uni[2*k4]+\
cword[i][ind[5]]*uni[2*k5] +cword[i][ind[6]]*uni[2*k6] +[1,0] ;\
if( divrem(sumz[1],3)[2] == 0 || divrem( sumz[2] ,3)[2] == 0,\
cou = cou+1; word = matrix(12,2); \
word[ind[1],] = (cword[i][ind[1]])*t;                  word[ind[2],] = (cword[i][ind[2]])*dot( t, uni[2*k2]);\
word[ind[3],] = (cword[i][ind[3]])*dot( t, uni[2*k3]); word[ind[4],] = (cword[i][ind[4]])*dot( t, uni[2*k4]);\
word[ind[5],] = (cword[i][ind[5]])*dot( t, uni[2*k5]); word[ind[6],] = (cword[i][ind[6]])*dot( t, uni[2*k6]);\
write( outttt, "leech2[",cou +1980,"]=",word, ";") )   )))))  );
return(cou);
}

Leech24()=
{
local(mult, cword, i,j,k,r, word,cou);
mult = monomial();
word = matrix(12,2);cou= 0;
for( i = 1, 12, for( j = i+1, 12, cword = [1,0;1,0;1,0;1,0;1,0;1,0;1,0;1,0;1,0;1,0;1,0;1,0]; cword[i,] = [-2,0]; cword[j,] = [-2,0];\
 for( k = 1, 243, cou = cou+1; for( r = 1, 12, word[r,] = dot( mult[k][r,], cword[r,]) );\
 write( outttt, "leech2[",cou+ 10890,"]=",word, ";")   )));
return(cou);
}

Leech25()=
{
local(mult, cword, i,k,r, word,cou);
mult = monomial();
word = matrix(12,2);cou= 0;
for( i = 1, 12, cword = [1,0;1,0;1,0;1,0;1,0;1,0;1,0;1,0;1,0;1,0;1,0;1,0]; cword[i,] = [1,3]; \
 for( k = 1, 243, cou = cou+1; for( r = 1, 12, word[r,] = dot( mult[k][r,], cword[r,]) );\
 write( outttt, "leech2[",cou +26928,"]=",word, ";")   ));
return(cou);
}

Leech26()=
{
local(mult, cword, i,k,r, word,cou);
mult = monomial();
word = matrix(12,2);cou= 0;
for( i = 1, 12, cword = [1,0;1,0;1,0;1,0;1,0;1,0;1,0;1,0;1,0;1,0;1,0;1,0]; cword[i,] = [-2,-3]; \
 for( k = 1, 243, cou = cou+1; for( r = 1, 12, word[r,] = dot( mult[k][r,], cword[r,]) );\
 write( outttt, "leech2[",cou +29844,"]=",word, ";")   ));
return(cou);
}

\\***************************************************Leech+H*********************************************************************

\\inner product of two vectors in the Eisenstein Lattice Leech+H ( in Wilson's coordinates)

ipw(x,y)=
{
local(ipp, i);
ipp = 0; 
for( i= 1, 12, ipp= ipp - (1/3)*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);
}

phiw(x,a)= a - sdot( dot([1, -1],ipw(x,a))/(-3), x );

\\given a set of vectors e from Leech+H finds the inner product matrix

checkipw(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] = ipw(e[i],e[j])));
return(gram);
} 

maw(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=phiw(x,unitvec);\
  for ( i = 1, n , matr[i,r] = ph[i,]) );
return(matr);
}

\\*********************************************************************************************************************************
\\given an eisenstein integer 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);
}

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

imbyt(a,b) = rea(divid(1/3*ipd(a,b), t)); 

ee = [-2, 0; 1, 0; 1, 0; 1, 0; -2, 0; 1, 0; 1, 0; 1, 0; 1, 0; 1, 0; 1, 0; 1, 0];

\\generate the vectors in the 1st shell of the leech lattice that are at a distance -6 from the vector ee.
\\these vectors are stored in lis

list6(ee)=
{
local(cou,i,j,tv,tvv);
cou = 0;
for( i = 1, 32760, for ( j = 1, 6 , tv = sdot( uni[j], leech2[i]); tvv = vad( tv, -ee); if( ipd(tvv,tvv) == [-18, 0], \
cou = cou+1; write(Lis, "list[",cou,"] =", tv,";") ) )); 
return(cou);
}

\\input a vector ee and the vectors "lis" in the 1st shell of leech lattice that are at a distance 6 from ee, 
\\ output is a maximal simplex with vectors from "lis" containing ee, such that distance between any two vectors is "-6"

dist6(ee,lis)=
{
local(ind, d6, cou, ocou, oldlist, newlist, tv, i);
ind = 1; d6 = vector(24);  \\ might need bigger than 17 for other bigger simplices that may lie in 1st shell of Leech
d6[1] = ee;
cou = matsize(lis)[2];
ocou = cou;
oldlist = lis;
while( cou > 0,ind = ind+1; d6[ind] = oldlist[1];\
newlist = vector( cou); cou = 0;\
for( i = 1, ocou, tv = vad( oldlist[i], - d6[ind]);\
 if( ipd(tv,tv) == [-18,0],  cou = cou+1 ; newlist[cou] = oldlist[i] ));\
print(ind, " ",cou); oldlist = vector(cou); for( i = 1, cou, oldlist[i] = newlist[i] ); ocou = cou);
print(ind);
return(d6);
}

checkimip(e)=
{
local();
ll =matsize(e)[2];
matimm = matrix(ll,ll);
for( i = 1, ll, for( j = 1, ll, matimm[i,j] = imbyt(e[i], e[j] ) ));
return(matimm);
}

\\**********************************Some output**********************************************************************************
\\the entries of the simp are 24 vectors from the 1st shell such that the difference between any two of them is again in the 1st shell
simp = vector(24);

simp[1] =[-2, 0; 1, 0; 1, 0; 1, 0; -2, 0; 1, 0; 1, 0; 1, 0; 1, 0; 1, 0; 1, 0; 1, 0];
simp[2] =[-3, 0; 3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0];
simp[3] =[-3, 0; 0, 0; 3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0];
simp[4] =[-3, 0; 0, 0; 0, 0; 3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0];
simp[5] =[-3, 0; 0, 0; 0, 0; 0, 0; 0, 3; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0];
simp[6] =[-3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0];
simp[7] =[-3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0];
simp[8] =[-3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 3, 0; 0, 0; 0, 0; 0, 0; 0, 0];
simp[9] =[-3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 3, 0; 0, 0; 0, 0; 0, 0];
simp[10] =[-3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 3, 0; 0, 0; 0, 0];
simp[11] =[-3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 3, 0; 0, 0];
simp[12] =[-3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 3, 0];
simp[13] =[-2, 1; 0, -1; 0, -1; 1, 1; 1, 1; 1, 1; 0, -1; 0, -1; 0, -1; 1, 1; 0, -1; 1, 1];
simp[14] =[-2, 1; 0, -1; 1, 1; 0, -1; 1, 1; 0, -1; 1, 1; 1, 1; 0, -1; 0, -1; 0, -1; 1, 1];
simp[15] =[-2, 1; 0, -1; 1, 1; 1, 1; 1, 1; 0, -1; 0, -1; 0, -1; 1, 1; 0, -1; 1, 1; 0, -1];
simp[16] =[-2, 1; 1, 1; 0, -1; 0, -1; 1, 1; 0, -1; 1, 1; 0, -1; 1, 1; 1, 1; 0, -1; 0, -1];
simp[17] =[-2, 1; 1, 1; 0, -1; 0, -1; 1, 1; 1, 1; 0, -1; 1, 1; 0, -1; 0, -1; 1, 1; 0, -1];
simp[18] =[-2, 1; 1, 1; 1, 1; 1, 1; 1, 1; 1, 1; 1, 1; 1, 1; 1, 1; 1, 1; 1, 1; 1, 1];
simp[19] =[-3, -1; 0, -1; 0, -1; 0, -1; 1, 1; 1, 1; 1, 1; 0, -1; 1, 1; 0, -1; 1, 1; 1, 1];
simp[20] =[-3, -1; 0, -1; 0, -1; 1, 1; 1, 1; 0, -1; 1, 1; 1, 1; 0, -1; 1, 1; 1, 1; 0, -1];
simp[21] =[-3, -1; 0, -1; 1, 1; 0, -1; 1, 1; 1, 1; 0, -1; 1, 1; 1, 1; 1, 1; 0, -1; 0, -1];
simp[22] =[-3, -1; 1, 1; 0, -1; 1, 1; 1, 1; 0, -1; 0, -1; 1, 1; 1, 1; 0, -1; 0, -1; 1, 1];
simp[23] =[-3, -1; 1, 1; 1, 1; 0, -1; 1, 1; 0, -1; 0, -1; 0, -1; 0, -1; 1, 1; 1, 1; 1, 1];
simp[24] =[-3, -1; 1, 1; 1, 1; 1, 1; 1, 1; 1, 1; 1, 1; 0, -1; 0, -1; 0, -1; 0, -1; 0, -1];

GG = matrix(24,24);

 GG[1, ] =[0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, -1, 0, 0, 0, 0, 0, 0];
 GG[2, ] =[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, -1, -1, -1];
 GG[3, ] =[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, -1, 0, -1, -1];
 GG[4, ] =[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, -1, 0, -1, 0, -1];
 GG[5, ] =[-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0];
 GG[6, ] =[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, -1, 0, -1, 0, 0, -1];
 GG[7, ] =[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, -1, -1, 0, 0, 0, -1];
 GG[8, ] =[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, -1, -1, -1, 0, 0];
 GG[9, ] =[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, -1, 0, -1, -1, 0, 0];
 GG[10,] =[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, -1, -1, 0, -1, 0];
 GG[11,] =[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, -1, -1, 0, 0, -1, 0];
 GG[12,] =[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, -1, 0, 0, -1, -1, 0];
 GG[13,] =[-1, -1, -1, 0, -1, 0, -1, -1, -1, 0, -1, 0, 0, 0, 0, 0, 0, -1, -1, -1, -1, -1, -1, -1];
 GG[14,] =[-1, -1, 0, -1, -1, -1, 0, 0, -1, -1, -1, 0, 0, 0, 0, 0, 0, -1, -1, -1, -1, -1, -1, -1];
 GG[15,] =[-1, -1, 0, 0, -1, -1, -1, -1, 0, -1, 0, -1, 0, 0, 0, 0, 0, -1, -1, -1, -1, -1, -1, -1];
 GG[16,] =[-1, 0, -1, -1, -1, -1, 0, -1, 0, 0, -1, -1, 0, 0, 0, 0, 0, -1, -1, -1, -1, -1, -1, -1];
 GG[17,] =[-1, 0, -1, -1, -1, 0, -1, 0, -1, -1, 0, -1, 0, 0, 0, 0, 0, -1, -1, -1, -1, -1, -1, -1];
 GG[18,] =[1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0];
 GG[19,] =[0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0];
 GG[20,] =[0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0];
 GG[21,] =[0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0];
 GG[22,] =[0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0];
 GG[23,] =[0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0];
 GG[24,] =[0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0];

\\start with the 1st two vectors (below ) from the 1st shell, and successively find a vector in the 1st shell such that the gram matrix has 
\\minimum determinant. iterate this to make a list of 12 vectors , this makes a basis for the leech lattice with 12 vectors from the 1st shell.

findbasis()=
{
local(bold, b, i,j,testb,detold, detnew,innerp);
bold = vector(1);
bold[1] = [3, 0;-3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0];
for( i = 2, 12,  b =testb = vector(i); for( j = 1, i - 1,\
 b[j] = bold[j]; testb[j] = bold[j]); detold = 3^16;\
 for( j = 1, 32760, testb[i] = leech2[j]; innerp = checkipd(testb); detnew = abs(mdet(innerp)[1])/3^i;if( divrem(j,5000)[2] == 0, print(j));\
 if( detnew != 0 && detnew < detold ,b[i] = testb[i]; detold = detnew)); print(b); print(checkipd(b));  bold = b; print(mdet(checkipd(b))[1]/3^i) ) ;
return(b);
}

\\ a basis for leech consisting of 1st shell vectors found by starting with the 1st one on the list.

bb = vector(12);

 bb[1] =[3, 0; -3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0];
 bb[2] =[1, 2; -1, -2; 2, 1; 0, 0; 2, 1; 0, 0; 0, 0; -2, -1; -2, -1; 0, 0; 0, 0; 0, 0];
 bb[3] =[0, 0; 0, 0; 3, 0; 0, 0; 0, 0; 0, 0; 0, 0; -3, 0; 0, 0; 0, 0; 0, 0; 0, 0];
 bb[4] =[0, 0; 0, 0; 1, 2; 0, 0; 2, 1; 0, 0; 0, 0; -1, -2; -2, -1; 0, 0; 2, 1; -2, -1];
 bb[5] =[3, 0; 0, 0; 0, -3; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0];
 bb[6] =[0, 0; 1, 2; 1, -1; 0, 0; 0, 0; 0, 0; -2, -1; 1, 2; -2, -1; 0, 0; 1, -1; 0, 0];
 bb[7] =[0, 0; 0, 0; 0, 0; 0, 0; 1, 2; -2, -1; -2, -1; 0, 0; -2, -1; 0, 0; -2, -1; 1, -1];
 bb[8] =[-2, 0; 2, 2; 1, 0; 0, 1; -1, -1; 1, 0; 0, 1; 1, 0; 0, 1; 1, 0; -1, -1; 1, 0];
 bb[9] =[3, 0; 0, -3; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0];
 bb[10] =[3, 0; 0, 0; 0, 0; 0, 0; 0, -3; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0];
 bb[11] =[3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, -3; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0];
 bb[12] =[0, 0; 0, 0; 1, 2; 0, 0; -2, -1; -2, -1; 0, 0; -2, -1; -2, -1; 1, -1; 0, 0; 0, 0];

\\given a set of vectors e from the lattice Leech+H finds a set of linearly independent vectors in its perp, 
\\whose simple multiples are in the lattice. uses the basis BB for Leech+H defined inside the program 

eperpw(e)=
{
local(ht, l,i,j,Mb, cfeperpz,cfeperp,eperpp,Mbz, Mbj,BB);
ht = matsize(e)[2];
l = 14;
Mb = matrix( l, h ); 
BB = vector(14);
for( i = 1, 14, BB[i] = matrix(14,2));
for( i = 1, 12, for( j = 1, 12, BB[i][j,] = bb[i][j,] )); BB[13][13,1] = BB[14][14,1] = 1;
for( j = 1, ht, for( i = 1, l, Mb[i,j] = ipw(BB[i],e[j]) ));
Mbz = matrix(2*l, 2*h);
 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);
}

\\ calculate a lx2 matrix which corresponds in the Z form to taking inner product with a eisenstein vector
zmat(G)=
{
local(l1,l2,Gz,i,testv1,testv2);
l1 = matsize(G)[1]; l2 = matsize(G)[2];
Gz = matrix(2*l1, 2*l2);
for( i = 1, l1, testv1 = testv2 = matrix(l2,2); testv1[i,] = [1,0]; testv2[i,] = [0,1]; \
Gz[,2*i-1] = mtov( mv( G, testv1)); Gz[,2*i] = mtov( mv(G, testv2) ) );
return(Gz);
}

\\****************************************************************************************************************************
\\find the three orthogonal copies of E_8 inside Lorentzian leech lattice.

\\ check if the submatrix of GG corresponding to row & column given by the entries of v, form an E_8 
checke8(v)=
{
local(c124,c134,c234,ccc);
c124 = GG[v[1],v[2]] - GG[v[1],v[4]] + GG[v[2],v[4]]; if( c124^2 != 1, return(0)); 
c134 = GG[v[1],v[3]] - GG[v[1],v[4]] + GG[v[3],v[4]]; if( c134^2 != 1, return(0));
c234 = GG[v[2],v[3]] - GG[v[2],v[4]] + GG[v[3],v[4]]; ccc = c234 - c134; if( ccc^2 != 1, return(0));
return(1);
}

\\find the indices for which the corresponding entries of simp make an E_8
finde8()=
{
local(i1,i2,i3,i4,cou);
cou = 0;
for( i1 =1, 24, for( i2 = 1, 24, for( i3 = 1, 24, for( i4 = 1, 24, \
if( i1 != i2 && i1 != i3 && i1 != i4 && i2 != i3 && i2 != i4 && i3 != i4 ,\
if( checke8([i1,i2,i3,i4]) == 1, cou = cou +1;write( outttt, "inde8[",cou,"] = [", i1,",",i2,",",i3,",",i4, "];" ) ))))); print(i1));
return();
}

\\check if the submatrix of GG corresponding to row & column given by the entries of v1 & v2, form an E_8 + E_8

checke8_2(v1,v2)=
{
local();
 if( -GG[v1[1],v2[2]] + GG[v1[1],v2[1]] != -GG[v2[1], v2[4]] + GG[v2[2],v2[4]],return(0)); 
 if( -GG[v1[1],v2[3]] + GG[v1[1],v2[1]] != -GG[v2[1], v2[3]]  ,return(0));                 
 if( -GG[v1[1],v2[4]] + GG[v1[1],v2[1]] != -GG[v2[1], v2[4]]  ,return(0));                 
 
 if( -GG[v1[2],v2[1]] + GG[v1[1],v2[1]] !=  GG[v1[1], v1[4]] - GG[v1[2], v1[4]],return(0)); 
 if( -GG[v1[2],v2[2]] + GG[v1[1],v2[1]] != -GG[v2[1], v2[4]] + GG[v2[2], v2[4]] + GG[v1[1], v1[4]] - GG[v1[2],v1[4]],return(0)); 
 if( -GG[v1[2],v2[3]] + GG[v1[1],v2[1]] != -GG[v2[1], v2[3]] + GG[v1[1], v1[4]] - GG[v1[2],v1[4]],return(0)); 
 if( -GG[v1[2],v2[4]] + GG[v1[1],v2[1]] != -GG[v2[1], v2[4]] + GG[v1[1], v1[4]] - GG[v1[2],v1[4]],return(0)); 
 
 if( -GG[v1[3],v2[1]] + GG[v1[1],v2[1]] !=  GG[v1[1], v1[3]],return(0)); 
 if( -GG[v1[3],v2[2]] + GG[v1[1],v2[1]] != -GG[v2[1], v2[4]] + GG[v2[2], v2[4]] + GG[v1[1], v1[3]],return(0)); 
 if( -GG[v1[3],v2[3]] + GG[v1[1],v2[1]] != -GG[v2[1], v2[3]] + GG[v1[1], v1[3]],return(0)); 
 if( -GG[v1[3],v2[4]] + GG[v1[1],v2[1]] != -GG[v2[1], v2[4]] + GG[v1[1], v1[3]],return(0)); 
 
 if( -GG[v1[4],v2[1]] + GG[v1[1],v2[1]] !=  GG[v1[1], v1[4]],return(0)); 
 if( -GG[v1[4],v2[2]] + GG[v1[1],v2[1]] != -GG[v2[1], v2[4]] + GG[v2[2], v2[4]] + GG[v1[1], v1[4]],return(0)); 
 if( -GG[v1[4],v2[3]] + GG[v1[1],v2[1]] != -GG[v2[1], v2[3]] + GG[v1[1], v1[4]],return(0)); 
 if( -GG[v1[4],v2[4]] + GG[v1[1],v2[1]] != -GG[v2[1], v2[4]] + GG[v1[1], v1[4]],return(0)); 
return(1);
}

\\find indices for which the corresponding entries of the simp make E_8+E_8

finde8_2()=
{
local(chk);
for( i = 1, 5484, for( j = 1, i-1, if( setintersect( Set(inde8[i]), Set(inde8[j])) == [], \
chk = checke8_2(inde8[i], inde8[j] ); print(chk);\
if( chk == 1, write(outttt,  i, " ", j, " ", inde8[i], " ", inde8[j] )))); if( divrem(i,100)[2] == 0, print(i)));
return();
}

\\given the 8 indices corresponding to 8 entries of simp, calculate the 8 vectors that form E_8+E_8

getbeta(v1,v2)=
{
local(beta1,beta2, i,jj, chk, gen);
beta1 = beta2  = vector(4);
beta1[1] = -1/2; beta2[1] = beta1[1] + GG[ v1[1], v2[1]] ; 
beta1[2] = beta1[1] + GG[v1[1],v1[4]] - GG[v1[2],v1[4]] ; beta1[3] = beta1[1] + GG[v1[1],v1[3]]; beta1[4] = beta1[1] + GG[v1[1],v1[4]];
beta2[2] = beta2[1] + GG[v2[1],v2[4]] - GG[v2[2],v2[4]] ; beta2[3] = beta2[1] + GG[v2[1],v2[3]]; beta2[4] = beta2[1] + GG[v2[1],v2[4]];
gen = vector(8);
for( i = 1, 8, gen[i] = matrix(14,2));
for( i = 1, 4, for( jj = 1, 12, gen[i  ][jj,] = simp[v1[i]][jj,] ) ; gen[i  ][13,] = [1,0]; gen[i  ][14,] = [beta1[i],0] + 1/2* t   );
for( i = 1, 4, for( jj = 1, 12, gen[i+4][jj,] = simp[v2[i]][jj,] ) ; gen[i+4][13,] = [1,0]; gen[i+4][14,] = [beta2[i],0] + 1/2* t   );
for( i = 1, 8, print(gen[i]));
chk = checkipw(gen);
return(chk) 
}
\\*****************************************************************************************************************

\\ uptill here we can find 8 vectors forming an E_8+E_8 . calculate the eperp by lattice reduction algorithm and the one can easily find a hyperbolic cell in the perp , hence getting coordinates for L1 = E_8+ E_8 + H\\next calculate the perp again by lattice reduction and get a basis for the complementary lattice K isomorphic to E_8, It is a positive definite lattice. 
\\By grid search its easy to list its 40 roots ( upto unity) and then find a simple system out of it.

\\one set of resulting 14 vectors is given below

e8s = vector(14);
e8s[1] = [-3, 0; 3, 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];
e8s[2] =[-2, 1; 0, -1; 0, -1; 1, 1; 1, 1; 1, 1; 0, -1; 0, -1; 0, -1; 1, 1; 0, -1; 1, 1; 1, 0; 0, 1];
e8s[3] =[-3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 1, 0; 0, 1];
e8s[4] =[-3, -1; 1, 1; 0, -1; 1, 1; 1, 1; 0, -1; 0, -1; 1, 1; 1, 1; 0, -1; 0, -1; 1, 1; 1, 0; -1, 1];

e8s[5] =[-3, 0; 0, 0; 0, 0; 3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 1, 0; 0, 1];
e8s[6] =[-2, 1; 1, 1; 0, -1; 0, -1; 1, 1; 0, -1; 1, 1; 0, -1; 1, 1; 1, 1; 0, -1; 0, -1; 1, 0; 0, 1];
e8s[7] =[-3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 3, 0; 1, 0; 0, 1];
e8s[8] =[-3, -1; 1, 1; 1, 1; 1, 1; 1, 1; 1, 1; 1, 1; 0, -1; 0, -1; 0, -1; 0, -1; 0, -1; 1, 0; -1, 1];

e8s[9] =[-14, -3; 3, 1; 0, -2; 3, 1; 5, 2; 1, 0; 3, 1; 0, -2; 1, 0; 2, -1; -4, -4; 3, 1; 5, 1; -1, 5];
e8s[10] =[9, -1; -2, 1; 1, 1; -2, 1; -4, 0; 0, 2; -2, 1; 1, 1; 0, 2; -1, 0; 2, 3; -2, 1; -3, 1; -1, -5];
e8s[11] =[-1, 2; -1, -1; -1, -1; -1, -1; -1, -1; -1, -1; -1, -1; -1, -1; -1, -1; -1, -1; -1, -1; -1, -1; 0, -1; 1, 1];
e8s[12] =[13, -4; -4, 1; 2, 4; -4, 1; -6, 0; 1, 2; -4, 1; 2, 4; 1, 2; 0, 3; 3, 3; -4, 1; -5, 2; -2, -8];

e8s[13] =[-11, 20; 4, -4; -5, -4; 4, -4; 7, -4; -2, -4; 4, -4; -5, -4; -2, -4; -2, -4; -8, -4; 4, -4; 4, -8; 9, 14];
e8s[14]=[-36, -39; 8, 10; 4, -4; 8, 10; 9, 15; 5, 1; 8, 10; 4, -4; 5, 1; 6, 3; 1, -7; 8, 10; 14, 14; -18, -1];

\\**************************************************************************************************************************************
\\another basis for 3E8 + H in Leech+ H that have smaller coordinates and all the simple roots of the E8's have ht 1.
\\this is found in isom2
e8s1 = vector(14);
e8s1[1]=[-3, 0; 0, 0; 3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 1, 0; 0, 1];
e8s1[2]=[2, -1; 0, 1; -1, -1; 0, 1; -1, -1; 0, 1; -1, -1; -1, -1; 0, 1; 0, 1; 0, 1; -1, -1; -1, 0; -1, -1];
e8s1[3]=[-3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 1, 0; 0, 1];
e8s1[4]=[3, 1; 0, 1; 0, 1; 0, 1; -1, -1; -1, -1; -1, -1; 0, 1; -1, -1; 0, 1; -1, -1; -1, -1; -1, 0; 0, -1];
e8s1[5]=[-3, 0; 3, 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];
e8s1[6]=[3, 1; -1, -1; 0, 1; -1, -1; -1, -1; 0, 1; 0, 1; -1, -1; -1, -1; 0, 1; 0, 1; -1, -1; -1, 0; 0, -1];
e8s1[7]=[-3, 0; 0, 0; 0, 0; 3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 1, 0; 0, 1];
e8s1[8]=[2, -1; 0, 1; 0, 1; -1, -1; -1, -1; -1, -1; 0, 1; 0, 1; 0, 1; -1, -1; 0, 1; -1, -1; -1, 0; -1, -1];
e8s1[9]=[-3, 0; 0, 0; 0, 0; 0, 0; 0, 3; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 1, 0; 0, 1];
e8s1[10]=[2, 0; -1, 0; -1, 0; -1, 0; 2, 0; -1, 0; -1, 0; -1, 0; -1, 0; -1, 0; -1, 0; -1, 0; -1, 0; 0, -1];
e8s1[11]=[-3, -2; 1, 0; 1, 0; 1, 0; 0, 1; -1, -1; 1, 0; -1, -1; 0, 1; 1, 0; 1, 0; -1, -1; 1, 1; 0, 1];
e8s1[12]=[-2, -2; 1, 1; 1, 1; 1, 1; -1, 0; 0, 2; 1, 1; -1, 0; 1, 1; 0, -1; 0, -1; -1, 0; 1, 1; -1, 0];
e8s1[13]=[12, -2; -3, 1; -3, 1; -3, 1; -1, -4; -2, 0; -3, 1; -1, -1; -3, -2; -2, 0; -2, 0; -1, -4; -5, 0; -2, -6];
e8s1[14]=[12, 16; -4, -4; -4, -4; -4, -4; 4, 0; 0, -2; -4, -4; 2, -1; 0, -2; -2, -3; -2, -3; 2, -1; -4, -6; 6, 0];

\\This matrix e8s1 is called E1 in the paper. This has the same inner product "checkipw" as f1,w*e1,d1,w*c1, (0^12,0,1), (0^12,1,0)

\\writing down the change of basis matrix between the two coordinate systems of L.

e8s2 = vector(14);
e8s2[ 1] = f1; e8s2[ 2] = sdot(w,e1); e8s2[ 3] = d1; e8s2[ 4] =sdot(w, c1);
e8s2[ 5] = f2; e8s2[ 6] = sdot(w,e2); e8s2[ 7] = d2; e8s2[ 8] =sdot(w, c2);
e8s2[ 9] = f3; e8s2[10] = sdot(w,e3); e8s2[11] = d3; e8s2[12] =sdot(w, c3);
e8s2[13] = [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];
e8s2[14] = [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,0];

\\then the ip matrix of the e8s2 and the ipw matrix of e81s are the same. So a change of basis can be obtained by sending 
\\e8s1[i] to e8s2[i]. Left multiplication by C sends Leech + H --> 3E8 + H
\\ If the matrix of an automorphism is A in the co-ordinate Leech+ H, it becomes C*A*C^{-1} in the coordinate 3E8+ H

EE1 = EE2 = matrix(14,14);
for( j = 1, 14,for( i = 1, 14, EE1[i,j] = e8s1[ j][i,]));
for( j = 1, 14,for( i = 1, 14, EE2[i,j] = e8s2[ j][i,]));
C = md( EE2, minv(EE1));

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

\\translation T_{l,z} in the coordinate system Leech+H
\\inner product in the negative definite leech lattice is (1/3)*ipd ; 
\\there was a mistake in the earlier formula for Trans in the last row first 12 entries, but it does not matter
\\because the generators are obtained by hitting rr[1] and rr[2] by Trans, and
\\both rr[1] and rr[2] has first 12 coordinates zero.
 
Trans(l,z)=
{
local(T,i,j);
T = matrix(14,14);
for( i = 1, 14, for( j = 1, 14, T[i,j] = [0,0]); T[i,i] = [1,0]);
for( i = 1, 12, T[i,13] = l[i,]; T[14,i] = -divid(bar(l[i,]) ,3*t));
T[14,13] = -divid( z - ipd(l,l)/6  ,t);
return(T);
}

\\generators of the whole reflection group of L in coord system Leech+H
rr = vector(2);

rr[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; 1,0; -1,-1];
rr[2] = [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]; 

gen = vector(50);

for( i = 1, 12, gen[2*i-1] = mv(Trans(bb[i], [0,0]), rr[1] ); gen[2*i] = mv( Trans( sdot( w,bb[i] ),[0,0]),rr[1]) ;\
                gen[24+2*i-1] = mv(Trans(bb[i], [0,0]), rr[2] ); gen[24+2*i] = mv( Trans( sdot( w,bb[i] ),[0,0]),rr[2]) );
gen[49] = rr[1];
gen[50] = rr[2];

\\for( i = 1, 50, print(i," ", checkleech( vecextract(gen[i]~,2^12 -1 )~), " ",gen[i] ));

\\program to write down the 50 generators 
genge()=
{
local(tge,i);
tge=vector(50);
tge[49] = [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];
tge[50] = [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]; 
for( i = 1, 12, tge[2*i-1] = mv(Trans(bb[i], [0,0]), tge[49] );\
	        tge[2*i] = mv( Trans( sdot( w,bb[i] ),[0,0]),tge[49]) ;\
                tge[24+2*i-1] = mv(Trans(bb[i], [0,0]), tge[50] );\
                tge[24+2*i] = mv( Trans( sdot( w,bb[i] ),[0,0]),tge[50]) );
return(tge);
}


\\**************************************************************************************************************************************
n0=[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]; 

valencye8(m)=
{
local(val,i,j,cou1,cou2,ind1,ind2,sim);
val=vector(4); cou1 = cou2 = 0; ind1 = ind2 = sim = vector(4);
for(i =1, 4, for(j = 1, 4, if( i != j && ipw(liste8[m[i]],liste8[m[j]]) !=[0,0], val[i] = val[i]+1 ) ));
for( i = 1, 4, if(val[i] == 1, cou1 = cou1+1 ;ind1[cou1]= m[i] ); if(val[i] ==2, cou2 = cou2+1;ind2[cou2] = m[i] ) );
if( cou1 ==2 && cou2 == 2, print( m, " ",val );\
 sim[1]=liste8[ind1[1]]; sim[4]=liste8[ind1[2]];sim[2]=liste8[ind2[1]]; sim[3]=liste8[ind2[2]];\
if(ipw(sim[1],sim[2]) == [0,0], sim[2]=liste8[ind2[2]]; sim[3]=liste8[ind2[1]] );\
sim[2] = sdot( divid(t,ipw(sim[1],sim[2])) , sim[2]);\
sim[3] = sdot( divid(-t,ipw(sim[2],sim[3])) , sim[3]);\
sim[4] = sdot( divid(t,ipw(sim[3],sim[4])) , sim[4]); for(i = 1, 4, print(sim[i]));ipp=checkipw(sim); for(i=1,4, print(ipp[i,]))  );
return(sim);
}

checkval(v)=
{
local(ll,i1,i2,i3,i4);
ll = matsize(v)[2];
for(i1 = 1, ll , for(i2= i1+1, ll, for(i3 = i2+1, ll, for(i4 = i3 +1,ll, valencye8([i1,i2,i3,i4] ) ))));
return();
}

e8s1n=vector(14);

lf1=e8s1n[1]=[-3, 0; 0, 0; 3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 1, 0; 0, 1];
le1=e8s1n[2]=[2, -1; 0, 1; -1, -1; 0, 1; -1, -1; 0, 1; -1, -1; -1, -1; 0, 1; 0, 1; 0, 1; -1, -1; -1, 0; -1, -1];
ld1=e8s1n[3]=[-3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 1, 0; 0, 1];
lc1=e8s1n[4]=[3, 1; 0, 1; 0, 1; 0, 1; -1, -1; -1, -1; -1, -1; 0, 1; -1, -1; 0, 1; -1, -1; -1, -1; -1, 0; 0, -1];
lf2=e8s1n[5]=[-3, 0; 3, 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];
le2=e8s1n[6]=[3, 1; -1, -1; 0, 1; -1, -1; -1, -1; 0, 1; 0, 1; -1, -1; -1, -1; 0, 1; 0, 1; -1, -1; -1, 0; 0, -1];
ld2=e8s1n[7]=[-3, 0; 0, 0; 0, 0; 3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 1, 0; 0, 1];
lc2=e8s1n[8]=[2, -1; 0, 1; 0, 1; -1, -1; -1, -1; -1, -1; 0, 1; 0, 1; 0, 1; -1, -1; 0, 1; -1, -1; -1, 0; -1, -1];
lf3=e8s1n[9]=[-2, 0; 1, 0; 1, 0; 1, 0; 0, 1; 2, 2; 1, 0; 0, 1; 1, 0; -1, -1; -1, -1; 0, 1; 1, 0; 0, 1];
le3=e8s1n[10]=[2, -1; 0, 1; 0, 1; 0, 1; -1, -1; 1, 0; 0, 1; 1, 0; -1, -1; 0, 1; 0, 1; 1, 0; -1, 0; -1, -1];
ld3=e8s1n[11]=[-3, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 0; 0, 3; 1, 0; 0, 1];
lc3=e8s1n[12]=[2, 0; -1, 0; -1, 0; -1, 0; -1, 0; -1, 0; -1, 0; -1, 0; -1, 0; -1, 0; -1, 0; 2, 0; -1, 0; 0, -1];
    e8s1n[13]=[-5, -18; 2, 5; 2, 5; 2, 5; -4, -1; -1, 2; 2, 5; -2, 0; -2, 0; 1, 3; 1, 3; -5, -3; 1, 6; -8, -4];
    e8s1n[14]=[12, -2; -3, 1; -3, 1; -3, 1; -1, -4; -2, 0; -3, 1; -1, -1; -3, -2; -2, 0; -2, 0; -1, -4; -5, 0; -2, -6];

\\e8s1n is called E_1' in the paper. it has the same gram matrix ( checkipw) as that of e8s2 : f_i, w*e_i, d_i, w*c_i, n_1, n_2 
\\that is e8s1, e8s1n has the same gram matrix (checkipw) in Leech+H as e8s2 (checkip) in 3E8+H 

e8s3 = vector(14);
for( i = 1, 4, e8s3[i] = e8s1n[4+i]; e8s3[4+i] =e8s1n[i] );
for(i = 9,14, e8s3[i] = e8s1n[i] );

e8s4 = vector(14);
for( i = 1, 4, e8s4[i] = e8s1n[i]; e8s4[4+i] =e8s1n[8+i]; e8s4[8+i] = e8s1n[4+i]);
for(i = 13,14, e8s4[i] = e8s1n[i] );

EE1n=EE3=EE4=matrix(14,14);

for( j = 1, 14,for( i = 1, 14, EE1n[i,j]= e8s1n[ j][i,]));
for( j = 1, 14,for( i = 1, 14, EE3[i,j] = e8s3[ j][i,]));
for( j = 1, 14,for( i = 1, 14, EE4[i,j] = e8s4[ j][i,]));

q12= md( EE3, minv(EE1n));
q23=md( EE4, minv(EE1n));
Cn = md( EE2, minv(EE1n)); \\left mult by Cn is also an isom from Leech+h -> 3E8 + H, it takes (0^12, 0,1) to a vector that is invariant under interchanging the three E_8

\\q12 swaps the 1st two hands of Y and fixes the third q23 fixes the 1st swaps the other two. they are in R(L).

n0=[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]; 
n1=[-2, 0; -1, -1; -1, -1; 1, 0; 0, 1; 0, -2; 1, 0; -1, -1; 1, 0; 1, 0; 1, 0; 0, 1; 1, 0; 1, 2];
\\  -2   ,  w^2  ,  w^2  ,  1  ,  w  , -2w  ,  1  , w^2   ,  1  ,  1  ,   1 ,   w ,  1  ,   t 
\\n0, n1 are fixed by q12 and q23 and they form a hyperbolic cell whose perp is the Leech lattice 
\\this shows that q12 and q23 are in R(L) meet 6.Suz and so R(L) is everything.

low=vector(3);
la = vad( e8s1n[14], sdot(W,e8s1n[13] ));
lb1=low[1]=vad(-e8s1n[14],-vad(e8s1n[ 1],vad(sdot([-1,-2],e8s1n[ 2]),vad(sdot([2,0],e8s1n[ 3]),sdot([-1,-2],e8s1n[ 4]) ))));
lb2=low[2]=vad(-e8s1n[14],-vad(e8s1n[ 5],vad(sdot([-1,-2],e8s1n[ 6]),vad(sdot([2,0],e8s1n[ 7]),sdot([-1,-2],e8s1n[ 8]) ))));
lb3=low[3]=vad(-e8s1n[14],-vad(e8s1n[ 9],vad(sdot([-1,-2],e8s1n[10]),vad(sdot([2,0],e8s1n[11]),sdot([-1,-2],e8s1n[12]) ))));

\\la, lb1, lb2, lb3, e8s1n[i], i = 1, ..., 12, have the same inner product as
\\a,   b1,  b2,  b3, e8s2[i] , i = 1, ..., 12. Note that e8s2[i] are fi, w*e_i, d_i, w*c_i. 

genq12()=
{
local(q1,q2,q3,q4,q5,q6,h3,cox3,qq,qqq);
q1=md(maw(lf1), md(maw(le1), md(maw(ld1), md(maw(lc1), md(maw(lb1),md(maw(la),md(maw(lb2),  md(maw(lc2), md(maw(ld2), md(maw(le2), md(maw(lf2),\
md(maw(le2),md(maw(ld2), md(maw(lc2), md(maw(lb2),  md(maw(la), md(maw(lb1),md(maw(lc1),md(maw(ld1), md(maw(le1), maw(lf1)  ))))))))))))))))))));
q2=md(maw(le1), md(maw(ld1), md(maw(lc1), md(maw(lb1), md(maw(la),md(maw(lb2),md(maw(lc2), md(maw(ld2), md(maw(le2),\
md(maw(ld2), md(maw(lc2), md(maw(lb2),  md(maw(la), md(maw(lb1),md(maw(lc1),md(maw(ld1), maw(le1) ))))))))))))))));
q3=md(maw(ld1), md(maw(lc1), md(maw(lb1), md(maw(la),md(maw(lb2),  md(maw(lc2),md(maw(ld2),\
md(maw(lc2), md(maw(lb2),  md(maw(la), md(maw(lb1),md(maw(lc1), maw(ld1)))))))))))));
q4=md(maw(lc1), md(maw(lb1), md(maw(la),md(maw(lb2),  md(maw(lc2), md(maw(lb2),md(maw(la), md(maw(lb1), maw(lc1) ))))))));
q5=md(maw(lb1), md(maw(la),md(maw(lb2),  md( maw(la), maw(lb1)))));
q6=maw(la);
h3=md( maw(lf3), md( maw(le3), md(maw(ld3), maw(lc3))));
cox3=md(h3, md(h3, md(h3,md(h3, h3)))); cox3 = md(cox3, md(cox3, cox3));
qq=md(q1,md(q2,md(q3,md(q4,md(q5,q6)))));
qqq=-md(qq, cox3);
return(qqq);
}

genq23()=
{
local(q1,q2,q3,q4,q5,q6,h1,cox1,qq,qqq);
q1=md(maw(lf2), md(maw(le2), md(maw(ld2), md(maw(lc2), md(maw(lb2),md(maw(la),md(maw(lb3),  md(maw(lc3), md(maw(ld3), md(maw(le3), md(maw(lf3),\
md(maw(le3),md(maw(ld3), md(maw(lc3), md(maw(lb3),  md(maw(la), md(maw(lb2),md(maw(lc2),md(maw(ld2), md(maw(le2), maw(lf2)  ))))))))))))))))))));
q2=md(maw(le2), md(maw(ld2), md(maw(lc2), md(maw(lb2), md(maw(la),md(maw(lb3),md(maw(lc3), md(maw(ld3), md(maw(le3),\
md(maw(ld3), md(maw(lc3), md(maw(lb3),  md(maw(la), md(maw(lb2),md(maw(lc2),md(maw(ld2), maw(le2) ))))))))))))))));
q3=md(maw(ld2), md(maw(lc2), md(maw(lb2), md(maw(la),md(maw(lb3),  md(maw(lc3),md(maw(ld3),\
md(maw(lc3), md(maw(lb3),  md(maw(la), md(maw(lb2),md(maw(lc2), maw(ld2)))))))))))));
q4=md(maw(lc2), md(maw(lb2), md(maw(la),md(maw(lb3),  md(maw(lc3), md(maw(lb3),md(maw(la), md(maw(lb2), maw(lc2) ))))))));
q5=md(maw(lb2), md(maw(la),md(maw(lb3),  md( maw(la), maw(lb2)))));
q6=maw(la);
h1=md( maw(lf1), md( maw(le1), md(maw(ld1), maw(lc1))));
cox1=md(h1, md(h1, md(h1,md(h1, h1)))); cox1 = md(cox1, md(cox1, cox1));
qq=md(q1,md(q2,md(q3,md(q4,md(q5,q6)))));
qqq=-md(qq, cox1);
return(qqq);
}

\\e8s1 in verbatim form for the tex file:
\\[-3, 0] [0, 0] [ 3,0] [ 0,0] [ 0,0] [ 0,0] [ 0,0] [ 0,0] [ 0,0] [ 0,0] [ 0,0] [ 0,0] [ 1,0] [ 0,1]
\\[ 2,-1] [0, 1] -[1,1] [ 0,1] -[1,1] [ 0,1] -[1,1] -[1,1] [ 0,1] [ 0,1] [ 0,1] -[1,1] [-1,0] -[1,1]
\\[-3, 0] [0, 0] [ 0,0] [ 0,0] [ 0,0] [ 0,0] [ 3,0] [ 0,0] [ 0,0] [ 0,0] [ 0,0] [ 0,0] [ 1,0] [ 0,1]
\\[ 3, 1] [0, 1] [ 0,1] [ 0,1] -[1,1] -[1,1] -[1,1] [ 0,1] -[1,1] [ 0,1] -[1,1] -[1,1] [-1,0] [0,-1]
\\[-3, 0] [3, 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]
\\[ 3, 1] -[1,1] [ 0,1] -[1,1] -[1,1] [ 0,1] [ 0,1] -[1,1] -[1,1] [ 0,1] [ 0,1] -[1,1] [-1,0] [0,-1]
\\[-3, 0] [0, 0] [ 0,0] [ 3,0] [ 0,0] [ 0,0] [ 0,0] [ 0,0] [ 0,0] [ 0,0] [ 0,0] [ 0,0] [ 1,0] [ 0,1]
\\[ 2,-1] [0, 1] [ 0,1] -[1,1] -[1,1] -[1,1] [ 0,1] [ 0,1] [ 0,1] -[1,1] [ 0,1] -[1,1] [-1,0] -[1,1]
\\[-3, 0] [0, 0] [ 0,0] [ 0,0] [ 0,3] [ 0,0] [ 0,0] [ 0,0] [ 0,0] [ 0,0] [ 0,0] [ 0,0] [ 1,0] [ 0,1]
\\[ 2, 0] [-1,0] [-1,0] [-1,0] [ 2,0] [-1,0] [-1,0] [-1,0] [-1,0] [-1,0] [-1,0] [-1,0] [-1,0] [0,-1]
\\[-3,-2] [1, 0] [ 1,0] [ 1,0] [ 0,1] -[1,1] [ 1,0] -[1,1] [ 0,1] [ 1,0] [ 1,0] -[1,1] [ 1,1] [ 0,1]
\\[-2,-2] [1, 1] [ 1,1] [ 1,1] [-1,0] [ 0,2] [ 1,1] [-1,0] [ 1,1] [0,-1] [0,-1] [-1,0] [ 1,1] [-1,0]
\\[12,-2] [-3,1] [-3,1] [-3,1] -[1,4] [-2,0] [-3,1] -[1,1] -[3,2] [-2,0] [-2,0] -[1,4] [-5,0] -[2,6]
\\[12,16] -[4,4] -[4,4] -[4,4] [ 4,0] [0,-2] -[4,4] [2,-1] [0,-2] -[2,3] -[2,3] [2,-1] -[4,6] [ 6,0]


\\********************************************************************************************************
\\some calculation to check the indecomposible form of a quadratic form over 2 - adic numbers
\\need the calculation for the gram matrix of the real leech lattice 

\\calculate v_p(n)= power of p that divides n maximum value infinity= 10^7 (for 0)
va(p,n)=
{
local(vanum,vaden,num, den);
num = numerator(n); den = denominator(n);
if( n==0, return(10^7));
vanum=0;while( divrem(num,p)[2]==0 && vanum < 10^7, vanum++;num=num/p );
vaden=0;while( divrem(den,p)[2]==0 && vaden < 10^7, vaden++;den=den/p );
return(vanum-vaden);
}

trun(M,p)=
{
local(MM,n,i,j);
n=matsize(M)[1];
MM=matrix(n-1,n-1);
for(i=1,p-1, for(j=1,p-1, MM[i,j]= M[i,j]);for(j=p+1,n, MM[i,j-1]= M[i,j])  );
for(i=p+1,n, for(j=1,p-1, MM[i-1,j]= M[i,j]);for(j=p+1,n, MM[i-1,j-1]= M[i,j]) );
return(MM);
}

\\given a symmetric integer matrix want to break it into indecomposible one and two dimensional forms over Z_(2)
 
decomp(M)=
{
local(n,vd,vo,pp,p,q,i,j,alpha,beta,size,de,cou,dee);
de=vector(matsize(M)[1]); cou=0;
while(matsize(M)[1] > 2, n=matsize(M)[1]; cou++;\
vd=va(2,M[1,1]); pp=1; for(i = 2, n, if( va(2,M[i,i]) < vd, vd=va(2,M[i,i]); pp=i ));vo=va(2,M[1,2]); p=1;q=2; for(i = 2,n,for(j=i+1,n, if( va(2,M[i,j]) < vo, vo=va(2,M[i,j]); p=i; q=j ))); \
if( vd <= vo, \
for(i=1,n,if(i != pp, alpha = M[i,pp]/M[pp,pp];\
for(j = 1, n, M[i,j] = M[i,j] - alpha*M[pp,j]);for(j=1,n, M[j,i] = M[j,i] - alpha*M[j,pp] )));\
de[cou]= M[pp,pp]; M=trun(M,pp) );\
if(vd > vo,\
for(i=1,n,if(i != p && i != q,alpha = (M[q,q]*M[i,p] - M[p,q]*M[i,q])/( M[p,p]*M[q,q] - M[p,q]^2);\
beta  =(-M[p,q]*M[i,p] + M[p,p]*M[i,q])/( M[p,p]*M[q,q] - M[p,q]^2);\
for(j=1,n, M[i,j] = M[i,j] - alpha*M[p,j] - beta*M[q,j]);for(j=1,n,M[j,i] = M[j,i] - alpha*M[j,p] - beta*M[j,q]) ));\
de[cou]= [M[p,p],M[p,q];M[q,p],M[q,q] ] ; M=trun(trun(M,q),p)   ));
cou++; de[cou]=M;dee=vector(cou);for(i=1,cou,dee[i]=de[i]);
return(dee);
}

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