\\filename : ht.gp
\\here we do some floating point calculations using eisenstein numbers as usual complex numbers.
\\ the prefix "a" infront of the codes denotes that the numbers are in floating point mode.
\\need to read in inp.gp
\\ the symbols uses are exclusive from those in inp, so compatiable with them.
\\try to use reflections in the 26 nodes to decrease the height of a vector,
\\the height of a vector is obtained by taking inner product with the weyl vector, which is not in the lattice, 
\\this makes floating point calculation necessary
\\need to read in inp.gp 


 
at = sqrt(-3); aw = (-1 + at)/2; aW = (-1 - at)/2;  xi = exp( -Pi * I / 6);
auni = vector(6);
auni[1] = -aW; auni[2] = aw; auni[3] = -1; auni[4]= aW; auni[5] = -aw; auni[6] = 1;

\\the 50 generators of the whole reflection group transferred to the 3E8+H coordinates.

ge = vector(50);

ge[1]=[-2, -4; 0, -4; 4, -8; -10, 12; 0, -4; -3, -4; 3, -8; -9, 12; 4, 2; -2, -4; 6, -6; -10, 10; 8, 22; -24, -8];
ge[2]=[1, -1; 2, 0; 6, 2; -11, -5; 2, 0; 1, -2; 6, 1; -10, -5; 1, 2; 1, -1; 6, 3; -10, -5; -7, 4; -8, -12];
ge[3]=[0, -2; 2, -1; 6, -1; -10, 1; 2, -1; -1, -3; 5, -2; -9, 1; 2, 2; 0, -3; 5, -1; -10, 0; -2, 10; -13, -10];
ge[4]=[-1, -1; -1, -1; -1, -3; 0, 6; -1, -1; 0, -1; 0, -3; 1, 6; 1, 0; 0, 0; 1, -3; 0, 5; 5, 6; -6, 1];
ge[5]=[-1, -2; 0, -2; 2, -4; -2, 6; 0, -2; -2, -2; 1, -4; -2, 6; 1, 0; -1, -2; 1, -3; -3, 5; 4, 9; -9, -1];
ge[6]=[-1, -2; 0, -2; 2, -4; -5, 9; 0, -2; -1, -3; 2, -5; -4, 8; 3, 1; -1, -2; 3, -5; -5, 7; 6, 13; -15, -4];
ge[7]=[-1, -2; 0, -2; 2, -4; -4, 8; 0, -2; -2, -3; 1, -5; -2, 8; 3, 1; -1, -3; 2, -5; -5, 6; 6, 12; -13, -3];
ge[8]=[-1, -2; 0, -2; 2, -4; -7, 5; 0, -2; 0, -2; 3, -4; -6, 6; 2, 2; 0, -1; 5, -2; -6, 4; 3, 12; -13, -6];
ge[9]=[-2, -3; 0, -4; 4, -7; -10, 8; 0, -3; -2, -3; 4, -6; -8, 9; 3, 2; -1, -3; 5, -4; -10, 7; 5, 18; -20, -8];
ge[10]=[0, -2; 2, 0; 5, 1; -7, -1; 1, -1; 0, -2; 4, 0; -7, -1; 1, 1; 0, -1; 4, 0; -7, -2; -3, 5; -8, -8];
ge[11]=[-2, -3; 0, -2; 3, -4; -5, 5; 0, -2; -1, -1; 2, -3; -3, 6; 2, 1; -1, -2; 3, -3; -5, 5; 4, 11; -11, -4];
ge[12]=[0, -2; 0, -2; 2, -3; -4, 7; 0, -2; -2, -3; 1, -5; -4, 7; 2, 1; -1, -2; 3, -3; -5, 5; 4, 11; -12, -3];
ge[13]=[-2, -3; 1, -2; 4, -4; -8, 6; 1, -2; 0, -2; 4, -4; -7, 6; 3, 2; 0, -3; 5, -3; -7, 5; 3, 14; -15, -7];
ge[14]=[0, -2; 0, -1; 2, -2; -5, 3; 0, -1; -1, -1; 2, -2; -4, 3; 1, 1; 0, 0; 3, -1; -5, 3; 1, 7; -9, -4];
ge[15]=[-1, -1; 0, -1; 0, -3; 0, 4; 0, -1; 0, 0; 1, -2; 0, 4; 1, 0; -1, 0; 0, -3; 0, 4; 4, 5; -5, 0];
ge[16]=[-2, -3; -1, -3; 1, -7; -3, 13; -1, -3; -3, -3; 0, -7; -2, 12; 3, 1; -3, -4; 3, -6; -4, 11; 10, 17; -16, -1];
ge[17]=[-2, -3; 0, -3; 4, -6; -10, 9; 0, -3; -2, -4; 4, -7; -8, 8; 3, 2; -1, -3; 5, -4; -10, 7; 5, 18; -20, -8];
ge[18]=[0, -2; 1, -1; 4, 0; -8, -2; 1, -1; 1, -1; 5, 1; -6, 0; 1, 1; 0, -1; 4, 0; -7, -2; -3, 5; -8, -8];
ge[19]=[-2, -4; 0, -4; 4, -8; -10, 12; 0, -4; -2, -4; 4, -8; -8, 12; 4, 2; -3, -4; 5, -6; -11, 10; 8, 22; -24, -8];
ge[20]=[1, -1; 2, 0; 6, 2; -11, -5; 2, 0; 1, -1; 6, 2; -10, -4; 1, 2; 1, -2; 6, 2; -10, -6; -7, 4; -8, -12];
ge[21]=[-2, -4; 0, -3; 4, -7; -10, 10; 0, -4; -2, -3; 4, -7; -8, 10; 4, 3; -1, -3; 5, -5; -10, 8; 6, 19; -22, -8];
ge[22]=[1, -1; 1, -1; 5, 1; -9, -3; 2, 0; 0, -2; 5, 1; -8, -2; 0, 1; 0, -1; 5, 1; -8, -3; -4, 5; -8, -10];
ge[23]=[-2, -3; 0, -2; 1, -5; -3, 9; 0, -3; -2, -3; 1, -6; -2, 9; 3, 0; -1, -2; 2, -5; -3, 9; 7, 13; -13, -1];
ge[24]=[0, -2; 0, -2; 3, -4; -8, 5; 1, -1; 0, -2; 4, -3; -7, 5; 3, 3; -1, -2; 5, -2; -9, 3; 2, 12; -15, -8];
ge[25]=[-3, -4; -1, -4; 1, -10; -5, 16; -1, -4; -4, -4; 0, -10; -4, 16; 4, 1; -3, -4; 3, -8; -5, 14; 13, 22; -22, -2];
ge[26]=[0, -1; 1, 0; 3, 0; -6, -1; 1, 0; 0, -2; 3, -1; -5, -1; 1, 1; 0, -1; 3, 1; -5, -1; -2, 4; -6, -6];
ge[27]=[-1, -2; 1, -1; 3, -3; -5, 5; 1, -1; -2, -3; 2, -4; -4, 5; 2, 1; -1, -3; 2, -3; -5, 4; 3, 10; -11, -4];
ge[28]=[-2, -1; -2, -1; -4, -5; 5, 10; -2, -1; -1, -1; -3, -5; 6, 10; 1, -1; -1, 0; -2, -5; 5, 9; 10, 6; -4, 7];
ge[29]=[-2, -2; -1, -2; -1, -6; 3, 10; -1, -2; -3, -2; -2, -6; 3, 10; 1, -1; -2, -2; -2, -5; 2, 9; 9, 9; -7, 5];
ge[30]=[-2, -2; -1, -2; -1, -6; 0, 13; -1, -2; -2, -3; -1, -7; 1, 12; 3, 0; -2, -2; 0, -7; 0, 11; 11, 13; -13, 2];
ge[31]=[-2, -2; -1, -2; -1, -6; 1, 12; -1, -2; -3, -3; -2, -7; 3, 12; 3, 0; -2, -3; -1, -7; 0, 10; 11, 12; -11, 3];
ge[32]=[-2, -2; -1, -2; -1, -6; -2, 9; -1, -2; -1, -2; 0, -6; -1, 10; 2, 1; -1, -1; 2, -4; -1, 8; 8, 12; -11, 0];
ge[33]=[-3, -3; -1, -4; 1, -9; -5, 12; -1, -3; -3, -3; 1, -8; -3, 13; 3, 1; -2, -3; 2, -6; -5, 11; 10, 18; -18, -2];
ge[34]=[-1, -2; 1, 0; 2, -1; -2, 3; 0, -1; -1, -2; 1, -2; -2, 3; 1, 0; -1, -1; 1, -2; -2, 2; 2, 5; -6, -2];
ge[35]=[-3, -3; -1, -2; 0, -6; 0, 9; -1, -2; -2, -1; -1, -5; 2, 10; 2, 0; -2, -2; 0, -5; 0, 9; 9, 11; -9, 2];
ge[36]=[-1, -2; -1, -2; -1, -5; 1, 11; -1, -2; -3, -3; -2, -7; 1, 11; 2, 0; -2, -2; 0, -5; 0, 9; 9, 11; -10, 3];
ge[37]=[-3, -3; 0, -2; 1, -6; -3, 10; 0, -2; -1, -2; 1, -6; -2, 10; 3, 1; -1, -3; 2, -5; -2, 9; 8, 14; -13, -1];
ge[38]=[-1, -2; -1, -1; -1, -4; 0, 7; -1, -1; -2, -1; -1, -4; 1, 7; 1, 0; -1, 0; 0, -3; 0, 7; 6, 7; -7, 2];
ge[39]=[-2, -1; -1, -1; -3, -5; 5, 8; -1, -1; -1, 0; -2, -4; 5, 8; 1, -1; -2, 0; -3, -5; 5, 8; 9, 5; -3, 6];
ge[40]=[-3, -3; -2, -3; -2, -9; 2, 17; -2, -3; -4, -3; -3, -9; 3, 16; 3, 0; -4, -4; 0, -8; 1, 15; 15, 17; -14, 5];
ge[41]=[-3, -3; -1, -3; 1, -8; -5, 13; -1, -3; -3, -4; 1, -9; -3, 12; 3, 1; -2, -3; 2, -6; -5, 11; 10, 18; -18, -2];
ge[42]=[-1, -2; 0, -1; 1, -2; -3, 2; 0, -1; 0, -1; 2, -1; -1, 4; 1, 0; -1, -1; 1, -2; -2, 2; 2, 5; -6, -2];
ge[43]=[-3, -4; -1, -4; 1, -10; -5, 16; -1, -4; -3, -4; 1, -10; -3, 16; 4, 1; -4, -4; 2, -8; -6, 14; 13, 22; -22, -2];
ge[44]=[0, -1; 1, 0; 3, 0; -6, -1; 1, 0; 0, -1; 3, 0; -5, 0; 1, 1; 0, -2; 3, 0; -5, -2; -2, 4; -6, -6];
ge[45]=[-3, -4; -1, -3; 1, -9; -5, 14; -1, -4; -3, -3; 1, -9; -3, 14; 4, 2; -2, -3; 2, -7; -5, 12; 11, 19; -20, -2];
ge[46]=[0, -1; 0, -1; 2, -1; -4, 1; 1, 0; -1, -2; 2, -1; -3, 2; 0, 0; -1, -1; 2, -1; -3, 1; 1, 5; -6, -4];
ge[47]=[-3, -3; -1, -2; -2, -7; 2, 13; -1, -3; -3, -3; -2, -8; 3, 13; 3, -1; -2, -2; -1, -7; 2, 13; 12, 13; -11, 5];
ge[48]=[-1, -2; -1, -2; 0, -6; -3, 9; 0, -1; -1, -2; 1, -5; -2, 9; 3, 2; -2, -2; 2, -4; -4, 7; 7, 12; -13, -2];
ge[49]=[0, 0; 1, 0; 1, 0; -2, 0; 1, 0; 0, 0; 1, 0; -1, 0; 0, 0; 0, 0; 2, 1; -2, -1; -1, 1; -2, -2];
ge[50]=[-1, 0; 0, 0; -2, -2; 3, 4; 0, 0; -1, 0; -2, -2; 4, 4; 0, -1; -1, 0; -1, -1; 3, 3; 4, 1; 0, 4];

\\given a matrix of eisenstein integers in discrete form returns the complex matrix

amat(a)=
{
local(l1,l2,ana,i,j);
l1 = matsize(a)[1];l2 = matsize(a)[2];ana= matrix(l1,l2);
for( i = 1, l1, for( j = 1, l2, ana[i,j] = a[i,j][1] + aw* a[i,j][2] ));
return(ana);
}

\\given a vector of eisenstein integers in discrete form returns the complex vector

avec(a)=
{
local(l1,ana,i);
ww = (-1 + sqrt(-3))/2;
l1 = matsize(a)[1];ana= vector(l1);
for( i = 1, l1, ana[i] = a[i,1] + aw* a[i,2] );
return(ana);
}

\\the inner product in 3E8+H in continuous coordinates in floating point 

aip(x,y)=-conj(vecextract(x,2^12-1))*vecextract(y,2^12-1)~ - conj(x[13])*sqrt(3)*I*y[14] + conj(x[14])*sqrt(3)*I*y[13];

\\The reflection in floating point


\\the negative definite inner product in the floating point

aipd(x,y)= sum( i = 1, matsize(x)[2], -conj(x[i])*y[i] );

\\check the inner product matrix of a set of vectors e.

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

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

aro = vector(26);
for( i = 1, 13, aro[i] = avec(r[i]); aro[13+i] = xi * avec(r[13+i]) );

rob = vector(14); for( i = 1, 26, rob = rob + aro[i]/26 );

height( x) = aip(x, rob);

h0 = abs(height(rob));

\\for(j=14,26,for(k=1,6, testv = sdot(uni[k],r[j]);for(i =1,23165, if(testv ==ro[i],print(i," ", j," ", k ))))); 

ahta(x) = abs(height(avec(x)));

\\ image of a under reflection in a vector x 

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

aphi(x,a)= a - (1- aw)*(aip(avec(x),a)/aip(avec(x),avec(x) ) ) * avec(x) ;


phipm(x,a,j)=
{
local(phip);
if( j == 1,phip = phi(x,a));
if( j == 0,phip = phi(x, phi(x,a)));
return(phip);
}

aphipm(x,a,j)=
{
local(phip);
if( j == 1,phip = aphi(x,a));
if( j == 0,phip = aphi(x, aphi(x,a)));
return(phip);
}


\\Decrease the height of a vector in 3E8+H using the height function h(x) = |<x,rob>|

decreaseht(x)=
{
local(abshtr, testr1,testr2,i,flag,cou,track);
track = vector(10000);
while( ahta(x) > h0 + .00000001, cou++; abshtr = ahta(x);flag = 0;\
for( i = 1, 26, testr1 = phi( r[i], x); testr2 = phi( r[i], testr1);
if( ahta(testr1) < abshtr -.00000001, x = testr1; track[cou] = [i,1];print(x); flag = 1;break());
if( ahta(testr2) < abshtr -.00000001, x = testr2; track[cou] = [i,2];print(x); flag = 1;break()) );\
if( flag==0,print(cou); return(x))   );
print(cou);return(track);
}

\\Decrease the height of a vector in 3E8+H using the height function h(x) = |<x,cusp>| where 
\\cusp is the norm zero vector obtained by reducing height of (wp + t * wl) as much as possible


decreasehtcusp(x)=
{
local(abshtr, testr1,testr2,i,flag,cou,track,cusp);
track = vector(1000);
cusp = [1, 2; 1, 0; 3, 1; -5, -3; 1, 2; 1, 0; 3, 1; -5, -3; 1, 2; 1, 0; 3, 1; -5, -3; -6, -1; 0, -5];
while( n2(ip(x,cusp))[1] > 0.00000001, cou++; abshtr = n2(ip(x,cusp))[1] ;flag = 0;\
for( i = 1, 26, testr1 = phi( r[i], x); testr2 = phi( r[i], testr1);
if( n2(ip(testr1,cusp))[1] < abshtr -.00000001, x = testr1; track[cou] = [i,1];print(n2(ip(testr1,cusp))[1],"  ",x); flag = 1;break());
if( n2(ip(testr2,cusp))[1] < abshtr -.00000001, x = testr2; track[cou] = [i,2];print(n2(ip(testr2,cusp))[1],"  ",x); flag = 1;break()) );\
if( flag==0,print(cou); return(x))   );
print(cou);return(track);
}

\\Decrease the height of a cusp in 3E8+H using the height function h(x) = |<x,rob>| 

decreasehtwrtz(x)=
{
local(abshtr, testr1,testr2,i,flag,cou,track);
track = vector(1000);
while( ahta(x) >  .00000001, cou++; abshtr = ahta(x);flag = 0;\
for( i = 1, 26, testr1 = phi( r[i], x); testr2 = phi( r[i], testr1);
if( ahta(testr1) < abshtr -.00000001, x = testr1; track[cou] = [i,1];print(x); flag = 1;break());
if( ahta(testr2) < abshtr -.00000001, x = testr2; track[cou] = [i,2];print(x); flag = 1;break()) );\
if( flag==0,print(cou); return(x))   );
print(cou); return(x);
}


decreasehtnew(x,m)=
{
local(absht, tr,i,flag,cou,track,trackk, ii,jj);
track = vector(1000);
while( ahta(x) > h0 + .00000001, cou++; absht = ahta(x);flag = 0; i = 0;
while( i <= 51 && flag ==0 ,i++; jj = divrem(i,2)[2]; ii = (i+jj)/2;  tr = phipm( r[ii], x,jj); 
if( ahta(tr) < absht -.00000001, x = tr; track[cou] = i;print(ahta(x)," ",x); flag = 1  ) );\
if( flag==0,print(cou-1); trackk = vector(cou-1);for(i=1,cou-1,trackk[i]=track[i]);write( path,"trac1[",m,"]= ",trackk,";"); return(x)));
print(cou);trackk = vector(cou); for( i = 1, cou, trackk[i] = track[i]); write( path, "trac1[",m,"]=", trackk,";"); return(x);
}

checktrack(m)=
{
local(xx,i,ni,nj,nii,kk,k);
xx = ge[m];for(i = 1,matsize(trac1[m])[2],ni = trac1[m][i]; nj=divrem(ni,2)[2]; nii=(ni+nj)/2; xx=phipm(r[nii],xx,nj);print(xx) );
if( perturb[m] == 0, for(k=1,26,for(kk=1,6,if(sdot(uni[kk],r[k])==xx,print("yes  ",k," ",kk)))); return( abs(ahta(xx) - h0) < .00000000001 ));
xx= phi( ge[perturb[m]],xx);
for(i = 1,matsize(trac2[m])[2],ni = trac2[m][i]; nj=divrem(ni,2)[2]; nii=(ni+nj)/2; xx=phipm(r[nii],xx,nj); print(xx));
for(k=1,26,for(kk=1,6,if(sdot(uni[kk],r[k])==xx,print("yes ",k," ",kk))));return(abs( ahta(xx) - h0) < .00000000001);
}

perturb = vector(50); 

perturb[ 1]=3;  perturb[ 2]=3;  perturb[ 5]=3;   perturb[ 7]=3;   perturb[12]=4;   perturb[18]=4;
perturb[19]=4;  perturb[20]=4;  perturb[22]=4;   perturb[25]=4;   perturb[26]=4;   perturb[27]=4;
perturb[35]=4;  perturb[36]=3;  perturb[39]=3;   perturb[40]=6;   perturb[43]=3;   perturb[46]=3; perturb[48]=4;

runtrac(i)=
{
local(tvec);
tvec = decreasehtnew(ge[i],i); write(path," ");\
if( perturb[i] != 0, tvec = phi( ge[perturb[i]], tvec); decreasehtnew(tvec,i) ); write(path," ")   ;
return();
} 

\\****************************************************************************************************************************************
htm(x) = abs(aip(( ( amat(x)*rob~)~), rob)); 

decreasehtmat(x)=
{
local(absht, tr,i,flag,cou,track,trackk, ii,jj,trij);
track = vector(1000);
while( htm(x) > h0 + .00000001, cou++; absht = htm(x);flag = 0; i = 0;
while( i <= 51 && flag ==0 ,i++; jj = divrem(i,2)[2]; ii = (i+jj)/2; if( jj == 1,trij=ma(r[ii]));if(jj==0,trij=minv(ma(r[ii])));tr=md(x,trij); 
if( htm(tr) < absht -.00000001, x = tr; track[cou] = i;print(i,"  ",htm(x)); flag = 1  ) );\
if( flag==0,print(cou-1); trackk = vector(cou-1);for(i=1,cou-1,trackk[i]=track[i]);write( path,"trac1[",m,"]= ",trackk,";"); return(x)));
print(cou);trackk = vector(cou); for( i = 1, cou, trackk[i] = track[i]); write( path, "trac1[",m,"]=", trackk,";"); return(x);
}

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];

htp(r) = n2( ip(wp, r))[1];
htl(r) = n2( ip(wl, r))[1];

\\****************************************************************************************************************************************
\\****************************************************************************************************************************************
\\****************************************************************************************************************************************
\\****************************************************************************************************************************************
\\****************************************************************************************************************************************
\\I think I wont need the codes below this any more. Need to check what they do and delete them
\\****************************************************************************************************************************************
\\****************************************************************************************************************************************
\\****************************************************************************************************************************************
\\****************************************************************************************************************************************
\\****************************************************************************************************************************************

decreasehtpair(x,m)=
{
local(htpx,htlx, tr,i,flag,cou,track,trackk, ii,jj);
track = vector(1000);
while( htp(x) > 0 && htl(x)> 3 , cou++; htpx = htp(x); htlx = htl(x);flag = 0; i = 0;
while( i <= 51 && flag ==0 ,i++; jj = divrem(i,2)[2]; ii = (i+jj)/2;  tr = phipm( r[ii], x,jj);
if( (htp(tr) < htpx && htl(tr) <= htlx) || (htp(tr) <= htpx && htl(tr) < htlx)   ,\
 x = tr; track[cou] = i;print(x); flag = 1  ) );\
if( flag==0,print(cou-1); trackk = vector(cou-1);for(i=1,cou-1,trackk[i]=track[i]);write( path,"trac1[",m,"]= ",trackk,";"); return(x)));
print(cou);trackk = vector(cou); for( i = 1, cou, trackk[i] = track[i]); write( path, "trac1[",m,"]=", trackk,";"); return(x);
}

Plus(X,a)=
{
local(cou, i,j,pl,mx,flag,pll,eps);
eps=10^(-12);
mx = matsize(X)[2];
pl = vector( mx + 2*mx^2 ); cou = 0;
for( i = 1, mx, cou = cou+1;pl[cou] = X[i]);
for(i=1,mx,for(j=1,mx,try = X[i] + a*X[j] ; flag = 0; for( k = 1, cou,if(norm(pl[k]-try)<eps,flag=1));if(flag==0,cou=cou+1;pl[cou]=try ));print(i));
for(i=1,mx,for(j=1,mx,try = X[i] + a^(-1)*X[j] ; flag = 0; for( k = 1, cou,if(norm(pl[k]-try)<eps,flag=1));if(flag==0,cou=cou+1;pl[cou]=try ));print(i));
pll=vector(cou); for(i=1, cou, pll[i]=pl[i]);
print(cou);
return(pll);
}


htofmat(seed, n,m)=
{
local(ind, i,j, matr, htt);
seed = divrem(seed, 400000)[2];ind = vector(m);
for( i = 1, n, vecr = rob~;for( j = 1, m, ind[j] = divrem(random( seed ),26)[2]+1; vecr = analm( ma(ro[ind[j]]))* vecr );\
htt = aip( vecr~, rob);\
if( norm( htt) < htm , print( ind, " " , norm(htt) )); print(i," ", norm(htt)) );
return();
}

\\generates a matrix for the C12 codeword wo, acting on Leech lattice as muliplicative golay code (Z/3)^5

multc12(wo)=
{
local(mt, i, j);
mt = Mid(12);
for( i = 1, 12, if( wo[i]==1, mt[i,i] = [0,1]); if( wo[i] == -1, mt[i,i] = [-1,-1]) );
return(mt);
}

\\check if the elements of the mul golay code (Z/3)^5 inside 6.Suz fixes the weyl vector rob.

checkrofix(m)=
{
local(suzd,suzl,suze,sr,j,i);
suzd = multc12(wordc12[m]);
for( j = 1,12, if(checkleech( mv( suzd, bb[j]) )!= 1, print(problem);return("problem" ) ) );
suzl= Mid(14);
for( j = 1, 12, suzl[j,j] = suzd[j,j] );  
suze = md( C, md( suzl, minv(C))); 
sr = vector(26);for( i = 1, 26, sr[i] = mv( suze, r[i] ) ); print( "ipcheck"," ",checkip(sr) == checkip(r));
if( norml2( rob - (amat( suze)* rob~)~ ) < 1, print("problem");return(problem) ); 
return();
}

\\mathewd is an element of order 5 in the M11 acting on Eisenstein Leech by permuting coords.

mathewd =matrix(12,12);
for( i = 1, 12, for( j = 1, 12, mathewd[i,j] = [0,0] ));
mathewd[1,1] =mathewd[3,2] =mathewd[4,3] =mathewd[5,4] =mathewd[6,5] =mathewd[2,6] =[1,0];
mathewd[7,7] =mathewd[9,8] =mathewd[10,9] =mathewd[11,10] =mathewd[12,11] =mathewd[8,12] =[1,0]; 


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

nfxi = nfinit(x^4 - x^2 +1 );

xidot(x, y) = nfeltmul( nfxi, x,y);

ipf(xx) = (aa + bb * xx^8)+ xx*(cc + dd*xx^8);
 
NN(aa,bb,cc,dd)=aa^4 - 2*bb*aa^3 + (3*bb^2 + (-cc^2 - 2*dd*cc + 2*dd^2))*aa^2 + (-2*bb^3 + (4*cc^2 - 4*dd*cc - 2*dd^2)*bb)*aa + (bb^4 + (-cc^2 + 4*dd*cc - dd^2)*bb^2 + (cc^4 - 2*dd*cc^3 + 3*dd^2*cc^2 - 2*dd^3*cc + dd^4));

N1(r)=NN( ip(r,sp)[1]  , ip(r,sp)[2] , ip(r,sl)[1] , ip(r,sl)[2] )/ 39^2 ;

robm = (avec(sp) -  xi* avec(sl))/26;

h0m = abs(aip(robm, robm));

N2(r) = abs((aip( avec(r), rob)/h0 )*(aip(avec(r), robm)/h0m))^2;

\\N1 and N2 compute the norm of the root r ( in discrete and floatig point).

mul13(r,eps)=
{
local(aa1,aa2,bb1,bb2);
aa1=ip(r,sp)[1];aa2=ip(r,sp)[2];bb1=ip(r,sl)[1];bb2=ip(r,sl)[2];
print( (bb1+bb2 + 4*aa1)/13 , "   ",(bb1 - 2*bb2 - 4*aa2)/13 );
if( eps==0, return( ((aa1+aa2)^2-3*bb2^2)/39 ) );
if( eps==1, return( ((bb1-2*bb2)^2-3*aa2^2)/39) );
}

\\**********************************************************************************************************************
\\Looking for a notion of positivity of roots
\\
\\ indicator function of the cone with arguement between (-a , a]
cone(x,a) = 
{
local(flag, ar);
flag = 0;
ar = arg( rea(x) + sqrt(-3) * ima(x) ) ;
if( ar <= a && ar > -a , flag = 1) ;
return( flag);
}

genpint()=
{
local(PP,wpv,wlv,i,j,pv,lv);
PP = matrix(28,28);
wpv=mtov(wp);wlv = mtov(wl);
for( j = 1, 13, pv= mtov(r[j]); lv=mtov(r[13+j]);for(i = 1, 28,PP[i,j]=pv[i]; PP[i,j+14]= lv[i]  )   );
for(i = 1,28, PP[i,14]=wpv[i]; PP[i,28] = wlv[i] );
return(PP);
}

pinv = genpint^(-1);

pcfint(x)=
{
local(iprob,i,cf,ux);
iprob = aip(rob,avec(x));
for( i = 1, 6, if( arg( auni[i]*iprob) <= Pi && arg( auni[i]*iprob) > -Pi,\
ux = sdot(uni[i], x); cf= pinv*mtov(ux); print(cf); print();  ));
return();
}

pcf(x)=
{
local(cf,i);
Pinv = minv(genPeis);
cf = sdot( t, mv( Pinv, x));
for( i = 1, 6, if( cone(dot( uni[i], cf[14,]), Pi/6) == 1, cf = sdot(uni[i], cf);print(i); return(cf) ));
return();
} 

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

\\ Trying to check if there is a root with absolute height less than one.

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

\\the hyperbolic distance between two points x,y or two hyperplanes determined by virtual points x,y are given by 
\\acosh(cd(x,y));
\\the hyperbolic distance between a point x and a hyperplane determined by virtual point y is given by 
\\asinh( cd(x,y));

cd(x, y)=
{
local(cdd);
if( matsize(x)[2] == 2, x = avec(x));
if( matsize(y)[2] == 2,  y = avec(y));
cdd =   abs( aip( x, y))/ (abs( aip(x,x)*aip(y,y) ) ^(.5))  ;
return(cdd);
}

\\distance of v from the 26 roots
mcd(v)=
{
local(dv,i,mdv);
dv = vector(26);
for( i = 1, 26, dv[i] = cd(v, r[i]) );
mdv = vecsort(dv)[1];
return(mdv);
}

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

\\move from \bar{\rho} = rob a little bit in random direction and see if one can decrease mcd

\\mcd0=mcd(rob);for(i=1,10,robp=rob+(avec(ran(14,5342))/10^8);print(mcd(robp));if(mcd(robp)>mcd0 + 10^(-10),print(robp," ",mcd(robp))));


\\bound on the absolute inner product of a root r of height one and a known root x

ipbd(x) = 3*cosh( asinh( sqrt(h0/3) ) + asinh( ahta(x)/sqrt(3*h0) ));

boundedE8()=
{
local(i,set1,set2,i1,i2,i3,i4,bdd,testv);
set1 = vector(13); set2 = vector(10);
for( i = 1, 6, set1[2*i - 1] = uni[i] ; set1[2*i] = dot( t, uni[i] ));set1[13] = [0,0];  \\zero and vectors in 1st 2 shell
set2[1] = [0,0];  set2[2] = [1,0];  set2[3] = [2,0];  set2[4] = [3,0];  set2[5] = [2,1]; \\zero and the vectors in 1st 4 shell 
set2[6] = [3,1];  set2[7] = [4,1];  set2[8] = [3,2];  set2[9] = [4,2];  set2[10]= [4,3]; \\ upto units.  
bdd = vector(2722);
cou = 0;
for( i1 = 1, 13, for( i2 = 1, 13, for( i3 = 1, 13, for( i4 = 1, 10, \
testv = matrix(4,2); testv[1,] = set1[i1] ; testv[2,] = set1[i2] ; testv[3,] = set1[i3] ; testv[4,] = set2[i4] ;\
if( checkE8(testv) == 1, cou = cou + 1 ; bdd[cou] = testv ) ))) );
print(cou);
return(bdd);
}

e8pairs()=
{
local(cou, i1, i2, k2, vv1,vv2,vv);
cou=0;
for( i1 = 1,782,for(i2= i1, 782,for(k2 =1,6,\
vv1= bdd[i1];vv2=sdot(uni[k2],bdd[i2]);\
if( n2( vv1[4,] - vv2[4,])[1] <=12 &&\
n2( vv1[3,] - vv2[3,])[1] <=12 &&\
n2( (vv1[1,] +vv1[3,]- vv1[4,])-(vv2[1,] +vv2[3,]- vv2[4,]) )[1] <= 36 &&\
n2((vv1[2,] -vv1[1,]- vv1[4,])-(vv2[2,] -vv2[1,]- vv2[4,]) )[1] <= 36 ,\
cou = cou+1; vv = [i1,i2,k2];write(e8bdd2,"pr[",cou,"] =", vv,";" ) ) )));
return(cou);
} 

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

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

phipm1(x,a,j)=
{
local(phip);
if( j == 1,phip = phi1(x,a));
if( j == 0,phip = phi1(x, phi1(x,a)));
return(phip);
}

r1 = vector(28);
for( i = 1, 26, r1[i] = r[i]); r1[27] = wp ; r1[28] = wl;

decreasehtnew1(x,m)=
{
local(absht, tr,i,flag,cou,track,trackk, ii,jj);
track = vector(1000);
while( ahta(x) > h0 + .00000001, cou++; absht = ahta(x);flag = 0; i = 0;
while( i <= 55 && flag ==0 ,i++; jj = divrem(i,2)[2]; ii = (i+jj)/2;  tr = phipm1( r1[ii], x,jj); 
if( ahta(tr) < absht -.00000001, x = tr; track[cou] = i;print(ahta(x)," ",x); flag = 1  ) );\
if( flag==0,print(cou-1); trackk = vector(cou-1);for(i=1,cou-1,trackk[i]=track[i]);write( path,"trac1[",m,"]= ",trackk,";"); return(x)));
print(cou);trackk = vector(cou); for( i = 1, cou, trackk[i] = track[i]); write( path, "trac1[",m,"]=", trackk,";"); return(x);
}

\\*********************************************************************
\\looking for the cusps
zz = vad( sp, sdot(-t, sl));
zzz = decreasehtnew(zz, 1); 

\\experimentally zzz is a null vector that is apparently not perpendicular to any root
\\and zzz has the minimal height that has been found among these. approx value of |<zzz, \rhobar>| = 4.27007073514
 
hzz= ahta(zzz);

\\see if any conjugate of zz has height lesser than zz  
tryref(z,k,n)=
{
local(j,ind,tv);
for( j = 1, n, \
ind = vector(k);\
for(i = 1, k, ind[i] = random(400000)+1); \
tv = z; for( i = 1, k, tv = phi(ro[ind[i]], tv) );\
if( ahta(tv) < hzz, print(ind); return(tv)  ) );
return();
}

\\for(i = 1, 420000, ttv = phi(ro[i], zz); if(  ahta(ttv) < hzz || ahta(phi(ro[i], ttv)) < hzz, print(i) ));


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

decreasehtwrt(z,x,m)=
{
local(absht, tr,i,flag,cou,track,trackk, ii,jj);
track = vector(1000);
while( n2(ip(z,x))[1] > h0 + .00000001, cou++; absht = n2(ip(z,x))[1];flag = 0; i = 0;
while( i <= 51 && flag ==0 ,i++; jj = divrem(i,2)[2]; ii = (i+jj)/2;  tr = phipm( r[ii], x,jj); 
if( n2(ip(z,tr))[1] < absht -.00000001, x = tr; track[cou] = i;print( n2(ip(z,x))," ",x); flag = 1  ) );\
if( flag==0,print(cou-1); trackk = vector(cou-1);for(i=1,cou-1,trackk[i]=track[i]);write( path,"trac1[",m,"]= ",trackk,";"); return(x)));
print(cou);trackk = vector(cou); for( i = 1, cou, trackk[i] = track[i]); write( path, "trac1[",m,"]=", trackk,";"); return(x);
}

\\*************************the null vector of minimal height of type leech found so far

\\zn=matrix(14,2);zn[14,1]=1;zn=mv(C,zn);zn=decreasehtnew(zn,1);

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


\\********define height of r as the norm of I(r)= Intersection of the mirror of r with F = <wp, wl> when you get stuck

Iu(x)= vad( sdot(ip(x,wp),wl), -sdot(ip(x,wl),wp));

Ih(x)=-ip(Iu(x),Iu(x))[1];

decreaseIh(x,m)=
{
local(absht, tr,i,flag,cou,track,trackk, ii,jj);
track = vector(1000);
while( Ih(x) > 0, cou++; absht = Ih(x);flag = 0; i = 0;
while( i <= 51 && flag ==0 ,i++; jj = divrem(i,2)[2]; ii = (i+jj)/2;  tr = phipm( r[ii], x,jj); 
if( Ih(tr) < absht -.00000001, x = tr; track[cou] = i;print(Ih(x)," ",x); flag = 1  ) );\
if( flag==0,print(cou-1); trackk = vector(cou-1);for(i=1,cou-1,trackk[i]=track[i]);write( path,"trac1[",m,"]= ",trackk,";"); return(x)));
print(cou);trackk = vector(cou); for( i = 1, cou, trackk[i] = track[i]); write( path, "trac1[",m,"]=", trackk,";"); return(x);
}


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

hexagon()=
{
local(rop, c, Gd, s,t,grid,i);
rop = vector(7);
rop[1]=rop[7]=rob; rop[4]=aphi(b1, aphi(a, aphi(b1, rob)));
rop[2] = aphi( a, rob); rop[3]=aphi(a , aphi(b1, rob));
rop[6] = aphi(b1, rob); rop[5]=aphi(b1, aphi(a , rob));
c = (rop[1] + rop[2] + rop[3] + rop[4] + rop[5] + rop[6])/6;
\\for( i = 1, 6, print("distance of rop[",i,"] and rop[",i+1,"]=", acosh( cd(rop[i], rop[i+1] )) ));
\\for( i = 1, 6, print("distance of rop[",i,"] and c=", acosh( cd(rop[i], c )) ));
print( "max possible height of a ropot that may meet Delta_1 or whose reflection in a may meet Delta_2 =",\
sinh( acosh(cd(rop[1], rop[2])))*sqrt(3/aip(rob,rob))   );
\\all values are between .05 and .3
Gd=100;
print("checking if the distance of a simple root from a point on Delta_1 is less than .05, using a grid length =",Gd);
for(s=0,Gd,for(t=0,Gd-s,grid=(s/Gd)*rop[1]+(t/Gd)*rop[2]+(1-(s+t)/Gd)*c;for(i=1,26,if(asinh(cd(grid,avec(r[i])))/2 <.05, print(s" ",t," ",i)) )));
print("checking if the distance of phi_a of a simple root from a point on Delta_2 is less than .05, using a grid length =",Gd);
for(s=0,Gd,for(t=0,Gd-s,grid=(s/Gd)*rop[2]+(t/Gd)*rop[3]+(1-(s+t)/Gd)*c;for(i=1,26,if(asinh(cd(grid,avec(phi(a,r[i]))))/2 <.05, print(s" ",t," ",i)) )));
return();
}


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

dist(x,y) = 2*abs(acosh( cd(x, y))); 

dist0=2*asinh( cd(avec(a), rob)); \\ distance between robar and the closest mirrors
\\****************************************************************************************************
\\Given a point x in the hyperbolic plane try to find the conjugate of \bar{\rho} that is closet to it, i.e. within dist0 with margin epsilon
\\try reflecting on both sides. i.e. when ro_i= M.rob, try all the vectors M.phi(r_i).rob and phi(r_i).M.rob, for 
\\$i = 1, 2, .., 26$ and see if the distance from x is reduced.
\\If we reach a M.rho that is dist0+epsilon of x then return [1,M]. Else return [0,M] where $M$ is the matrix
\\that took us as close as we could get to x.  
\\****************************************************************************************************

\\the function onestep takes an automorphism matrix M and a point x in the hyperbolic space and
\\returns [1,i,j] if dist( M.phi^j(r_i).rob, x) < dist( M.rob, x), returns [0,i,j] if dist( phi^j(r_i). M. rob, x) < dist( M.rob, rob)
\\and returns [0,0,0] otherwise. 

onestep(M,x)=
{
local(i,j,dnow);
dnow=dist( (M*rob~)~, x);
for( i=1,26, for(j=0,1, if( dist( (M*aphipm( r[i], rob, j)~)~, x) < dnow, return( [1,i,2*j-1] )))); 
for( i=1,26, for(j=0,1, if( dist( aphipm( r[i], (M*rob~)~, j), x) < dnow, return( [0,i,2*j-1] )))); 
return([0,0,0]);
}

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

decreaseht2(x)=
{
local(M,ind);
M=Mid(14); 
while( dist( (amat(M)*rob~)~, x) > .89, \
ind = onestep(amat(M),x); print(ind);  \
if(  ind[2] == 0 , return([0,M] )); \
if( ind[1] == 1 && ind[3] ==  1, M = md( M ,     ma( r[ind[2]])   )  );\
if( ind[1] == 1 && ind[3] == -1, M = md( M , minv(ma( r[ind[2]])) )  );\
if( ind[1] == 0 && ind[3] ==  1, M = md(ma( r[ind[2]]) ,  M       )  );\
if( ind[1] == 0 && ind[3] == -1, M = md(minv(ma( r[ind[2]])) , M  )  ));
return([1,M]);
}

deflate(II)=
{
local(ys,rop,rom,i,j, Gd, grid, Matr, Mat2,s,t);
ys=vector(12);
ys[ 1]=f1; ys[ 2]=e1; ys[ 3]=d1; ys[ 4]=c1; ys[ 5]=b1; ys[ 6]=a; ys[ 7]=b2; ys[ 8]=c2;ys[ 9]=d2; ys[10]=e2; ys[11]=f2; ys[12]=a3; 
rop=rom=vector(12);
for( i = 1, 12, rop[i] = rob; for( j = 1, i-1, rop[i] = aphi(ys[i - j], rop[i] ) ));
rom[1]=rob;for( i = 2, 12, rom[i] = rob; for( j = 2, i-1, rom[i] = aphi(ys[i - j], rom[i] ) ); rom[i] = aphi( ys[12], rom[i] ) );
Gd = 4;
for(i=II,II,print(i); print();for(s=0,Gd,for(t=0,Gd,grid= ( 1 - s/Gd)*((1-t/Gd)*rop[i]+ (t/Gd)*rop[i+1] ) + (s/Gd)*((1-t/Gd)*rom[i]+ (t/Gd)*rom[i+1] );\
Mat2=decreaseht2(grid); Matr = Mat2[2]; \
if( Mat2[1] == 0, print( "******Closest we could get is=",dist( (amat(Matr)*rob~)~, grid), "    for grid index = ", [i, s, t] ));\
if( Mat2[1] == 1, print(" could get it within distance " dist( (amat(Matr)*rob~)~, grid), "    for grid index = ", [i, s, t] );\
for( j=1,130, if( 2* asinh( cd( mv( Matr, ro[i]), grid)) < .2, print( "too close for gridindex =",[i,s,t], " root number=", j ))))   ))); 
return();
}



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

\\ Check some inner products (needed for a table used in a lemma in el2)

aipn(x,y) = aip(x,y)/h0;

ipxlzi()=
{
local(cp,ava,avf,pr1,pr2);
cp = 1 + 3^(-1/2); ava = avec(a); avf = avec(f);
pr1 = rob + (h0/3)*ava;  pr2= rob + (h0*xi/3)*avf;
print("column 1");
print(aipn(rob, ava) - 1);
print(aipn(rob, avf) - xi^(-1));
print();
print("column 2");
print(aipn( pr1, ava)             );
print(aipn( pr1, avf) - xi^(-1)*cp);
print();
print("column 3");
print(aipn( aphi(a, rob) , ava) -aW );                
print(aipn( aphi(a, rob) , avf) -(xi^(-1) - aW) );
print();
print("column 4");
print(aipn( aphi(a, pr2) , ava) - aW*cp);
print(aipn( aphi(a, pr2) , avf) + aW*cp);
print();
print("column 5");
print(aipn( aphi(a,aphi(f, rob)) , ava) -(aW -I) );
print(aipn( aphi(a,aphi(f, rob)) , avf) + aW     );
print();
print("column 6");
print(aipn( aphi(a,aphi(f, pr1)) , ava) + I*cp);
print(aipn( aphi(a,aphi(f, pr1)) , avf));
print();
print("column 7");
print(aipn( aphi(f,aphi(a,aphi(f, rob))) , ava) + I);
print(aipn( aphi(f,aphi(a,aphi(f, rob))) , avf) + aw );
print();
return()
}