\\filename: quat.gp

\\ quarternionic Lorentzian leech lattice

\\an element of the quaternion is represented as [a+b*I,c+d*I] = a+ b*I + (c + d*I)*J
\\the vectors of length n over quaternions are represented as matrices of size n by 2, the rows giving the entries. 

qi = [I,0]; qj = [0,1]; qk = [0,I]; p=[1-I,0]; P=[1+I,0];al=[(1 + I)/2, (1 + I)/2]; om = [ ( - 1 + I)/2, ( 1 + I)/2]; Om = -al; 

qunii = vector(8);
qunii[1] = [1,0] ; qunii[2] = [-1,0]; qunii[3] = qi; qunii[4] = -qi; qunii[5] = qj; qunii[6] = -qj; qunii[7] = qk ; qunii[8] = -qk;


quni = vector(24); quni[1] = [1,0] ; quni[2] = qi; quni[3] = qj; quni[4] = qk;
quni[5]=[(1 + I)/2, (1 + I)/2];quni[6]=[(1 + I)/2, (1 - I)/2];quni[7]=[(1 + I)/2, (-1 + I)/2];quni[8]=[(1 + I)/2, (-1 - I)/2];
quni[9]=[(1 - I)/2, (1 + I)/2];quni[10]=[(1 - I)/2, (1 - I)/2];quni[11]=[(1 - I)/2, (-1 + I)/2];quni[12]=[(1 - I)/2, (-1 - I)/2];
for( i = 1, 12, quni[25 - i] = -quni[i] );

qr= vector(14);
pp=p[1];
qr[ 1] = qa = [  0, 0  ;  0, 0  ;      0, 0  ;  0, 0  ;      0, 0  ;  0, 0  ;      1, 0  ; -1, 0  ];

qr[ 2] = qc1= [  1, 0  ; -1, 0  ;      0, 0  ;  0, 0  ;      0, 0  ;  0, 0  ;      0, 0  ;  0, 0  ];
qr[ 3] = qc2= [  0, 0  ;  0, 0  ;      1, 0  ; -1, 0  ;      0, 0  ;  0, 0  ;      0, 0  ;  0, 0  ];
qr[ 4] = qc3= [  0, 0  ;  0, 0  ;      0, 0  ;  0, 0  ;      1, 0  ; -1, 0  ;      0, 0  ;  0, 0  ];

qr[ 5] = qe1= [  0, 0  ;  0, 0  ;      1, 0  ;  1, 0  ;      1, 0  ;  1, 0  ;     -I, 0  ; -1, 0  ];
qr[ 6] = qe2= [  1, 0  ;  1, 0  ;      0, 0  ;  0, 0  ;      1, 0  ;  1, 0  ;     -I, 0  ; -1, 0  ];
qr[ 7] = qe3= [  1, 0  ;  1, 0  ;      1, 0  ;  1, 0  ;      0, 0  ;  0, 0  ;     -I, 0  ; -1, 0  ];

qr[ 8] = qf = [  0, 0  ; pp, 0  ;      0, 0  ; pp, 0  ;      0, 0  ; pp, 0  ;  -I*pp, 0  ;-pp, 0  ];

qr[ 9] = qb1= [  0, 0  ; pp, 0  ;      0, 0  ;  0, 0  ;      0, 0  ;  0, 0  ;     -1, 0  ;  0, 0  ];
qr[10] = qb2= [  0, 0  ;  0, 0  ;      0, 0  ; pp, 0  ;      0, 0  ;  0, 0  ;     -1, 0  ;  0, 0  ];
qr[11] = qb3= [  0, 0  ;  0, 0  ;      0, 0  ;  0, 0  ;      0, 0  ; pp, 0  ;     -1, 0  ;  0, 0  ];

qr[12] = qd1= [-pp, 0  ;  0, 0  ;      0, 0  ;  0, 0  ;      0, 0  ;  0, 0  ;      0, 0  ;  0, 0  ];
qr[13] = qd2= [  0, 0  ;  0, 0  ;    -pp, 0  ;  0, 0  ;      0, 0  ;  0, 0  ;      0, 0  ;  0, 0  ];
qr[14] = qd3= [  0, 0  ;  0, 0  ;      0, 0  ;  0, 0  ;    -pp, 0  ;  0, 0  ;      0, 0  ;  0, 0  ];


\\product of two quaternions

qdot(a,b) = [a[1]*b[1] - a[2]*conj(b[2]) , a[1]*b[2] + a[2]*conj(b[1]) ];

qn(a) = qdot(qconj(a),a);

\\conjugate of a quaternion

qconj(a) = [conj(a[1]), - a[2]];

qinv(a) = qconj(a)/qn(a)[1];

qreal(x) = real(x[1]);
qimag(x) = [ I*imag(x[1]), x[2]];  \\ x = qreal(x) + qimag(x)

\\ define the 6 roots of the form (0, 1, qeta), and of the form (0, xi, 1) where eta is an unit.
qreta = qrxi = qeta = vector(6);
qeta[1] = qi; qeta[2] = [-1,0];
qeta[3] = qdot(qj,al);qeta[5] = -qconj(al); qeta[4] = qdot(qi,al);qeta[5] = -qconj(al);  qeta[6] = -qdot(qconj(al),qk);
 
for( i = 1, 6, qreta[i] =qrxi[i] =  matrix(8,2); qreta[i][7,] = qrxi[i][8,] = [1,0]; qreta[i][8,] = qeta[i]; qrxi[i][7,] = qconj(qeta[i]) );

\\inner product of two element of 3*Lambda_2 + H ( Lambda_2 is the quaternionic E_8 lattice) 

qip(x,y)=-qdot(qconj(x[1,]),y[1,])-qdot(qconj(x[2,]),y[2,])-qdot(qconj(x[3,]),y[3,])\
-qdot(qconj(x[4,]),y[4,])-qdot(qconj(x[5,]),y[5,])-qdot(qconj(x[6,]),y[6,])\
+qdot(qdot(qconj(x[7,]),P),y[8,])+qdot(qdot(qconj(x[8,]),p),y[7,]);


qipd(x,y)=
{
local(i,qipdd);qipdd = [0,0];for( i = 1 , matsize(x)[1], qipdd = qipdd -qdot(qconj(x[i,]),y[i,]));return(qipdd);
}

qipw(x,y) = (-qdot(qconj(x[1,]),y[1,])-qdot(qconj(x[2,]),y[2,])-qdot(qconj(x[3,]),y[3,])\
-qdot(qconj(x[4,]),y[4,])-qdot(qconj(x[5,]),y[5,])-qdot(qconj(x[6,]),y[6,]) )/2 \
+qdot(qdot(qconj(x[7,]),P),y[8,])+qdot(qdot(qconj(x[8,]),p),y[7,]);


\\check the inner products of a set of vectors

checkqip(v)=
{
local(i,j,l,cck);
l = matsize(v)[2];
cck = matrix(l,l);
for( i = 1, l, for( j = 1, l, cck[i,j] = qip(v[i], v[j])  ));
return(cck);
}


checkqipw(v)=
{
local(i,j,l,cck);
l = matsize(v)[2];
cck = matrix(l,l);
for( i = 1, l, for( j = 1, l, cck[i,j] = qipw(v[i], v[j])  ));
return(cck);
}


\\scalar product of a vector and an quaternion, note that the scalar is on the right!

qsdot(x,f)=
{
local(i,sd);
sd=matrix(matsize(x)[1],matsize(x)[2]);
for ( i = 1, matsize(x)[1], sd[i,]=qdot(x[i,],f));
return(sd);
}

\\product of a matrix and an quaternion

qmdot( a, x) =
{
local( i,j,l1,l2,mdott );
l1 = matsize(x)[1]; l2 = matsize(x)[2];
mdott = matrix( l1,l2); 
for( i = 1, l1, for( j = 1, l2, mdott[i,j] = qdot( a, x[i,j]) ));
return(mdott);
}

\\add two vectors whose entries are quaternions
qvad(x,y)=
{
local(i,j,r,s,ad);
r=matsize(x)[1];
s=matsize(x)[2];
ad = matrix(r,s);
for( i = 1, r, for( j = 1, s, ad[i,j] = x[i,j] + y[i,j] ));
return(ad);
}

\\ multiply two matrices whose entries are quaternions

qmd(a,b)=
{
local(mdd,im,jm,km,i,j,k);
im=matsize(a)[1];
jm=matsize(a)[2];
km=matsize(b)[2];
mdd = matrix(im,km);
if( matsize(a)[2] == matsize(b)[1], for( i = 1, im, for ( k = 1, km,  for ( j = 1, jm, mdd[i,k] = mdd[i,k] + qdot(a[i,j], b[j,k]) ))));
return(mdd);
} 


\\multiply a matrix with a vector to get a vector, all entries in quaternions

qmv(a,b)=
{
local(mdd,im,jm,i,j);
im=matsize(a)[1];
jm=matsize(a)[2];
mdd = matrix(im,2);
 for( i = 1, im,   for ( j = 1, jm, mdd[i,] = mdd[i,] + qdot( a[i,j], b[j,]) ));
return(mdd);
}

\\generate the identity matrix of size n

Mid(n) =
{
local(mid, i,j);
mid= matrix(n,n);
for ( i = 1, n, for ( j = 1, n, mid[i,j] = [0,0]); mid[i,i]=[1,0]);
return(mid);
}

\\calculate the inverse of a matrix by sweep out

qminv(A) =
{
local(B,s1,s2,i,j,k,ii,jj,temp, fac);
s1=matsize(A)[1];
s2 = matsize(A)[2];
B = matrix(s1,s2);
for ( i = 1, s1, for ( j = 1, s2, B[i,j] = [0,0]); B[i,i]=[1,0]);
for( j = 1, s1, i = j; \
 while ( i < s1+1 && A[i,j] == [0,0] , i=i+1); if( i == s1+1 , return(singular));fac = A[i,j];\
 for( k = 1, s2,\
  temp = A[i,k]; A[i,k]=A[j,k];A[j,k]=qdot(qinv(fac),temp);temp = B[i,k]; B[i,k]=B[j,k];B[j,k]=qdot(qinv(fac),temp) );\
  for( ii = 1, s2 , if( ii != j,\
   fac = A[ii,j];\
   for( jj = 1, s2, A[ii,jj]=A[ii,jj]-qdot(fac,A[j,jj]);\
   B[ii, jj]=B[ii,jj]-qdot(fac,B[j,jj])\
    );\
   )));
return(B);
}

\\calculate the determinant of a matrix by upper triangularizing
qmdet(A) = 
{
local(B,s1,s2,i,j,k,ii,jj,temp, fac,de);
s1=matsize(A)[1];
s2 = matsize(A)[2];
de = [1,0];
for( j = 1, s1, i = j; \
if(A[i,i]==[0,0], while ( i < s1+1 && A[i,j] == [0,0] , i=i+1); if( i == s1+1 , return([0,0]));\
 for( k = 1, s2,temp = A[i,k]; A[i,k]=A[j,k];A[j,k]=temp ); de = qdot([-1,0],de) ); \
  for( ii = j+1, s2 ,\
   fac = qdot(qinv(A[j,j]),A[ii,j]);\
   for( jj = 1, s2, A[ii,jj]=A[ii,jj]-qdot(fac,A[j,jj])) ); de = qdot(de, A[j,j])\
   );
return(de);
}

\\calculate the order of a matrix 

qord(x)=
{
  local(n,y,mid,ll);
  ll= matsize(x)[1];
  mid = matrix(matsize(x)[1], matsize(x)[1]);
  for ( i = 1, ll, for ( j = 1, ll, mid[i,j] = [0,0]); mid[i,i]=[1,0]);
if( x == mid , return(1));
  n = 1;
  y=x;
  until(y == mid, n++; y = qmd(x,y);if(divrem(n,10)[2]==0,print(n)));
  return(n);
}

\\reflection of the vector a in the root x, wrt the unit un

qphi(x,a,un) = a - qsdot(x,qdot(([1,0] - un),qip( x,a)/(-2)));

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

qmaH(x,un)=
{
local(qmaa,l,qmah);
qmaa = qma(x,un);
l = matsize(qmaa)[1];
qmah = [qmaa[l-1,l-1],qmaa[l-1,l];qmaa[l,l-1],qmaa[l,l]];
return(qmah);
}

qTrans(l,z)=
{
local(T,i,j);
T = matrix(8,8);
for( i = 1, 8, for( j = 1, 8, T[i,j] = [0,0]); T[i,i] = [1,0]);
for( i = 1, 6, T[i,7] = l[i,]; T[8,i] = qdot(qinv(qconj(p)), qconj(l[i,]) ));
T[8,7] = qdot( qinv(qconj(p)),(z - qipd(l,l)/2) );
return(T);
}

qTransw(l,z)=
{
local(T,i,j);
T = matrix(8,8);
for( i = 1, 8, for( j = 1, 8, T[i,j] = [0,0]); T[i,i] = [1,0]);
for( i = 1, 6, T[i,7] = l[i,]; T[8,i] = qdot(qinv(qconj(p)), qconj(l[i,]) )/2);
T[8,7] = qdot( qinv(qconj(p)),(z - qipd(l,l)/4) );
return(T);
}

XX=[x11+ I*x12,x13 +I*x14 ; x21+ I*x22,x23 +I*x24  ;   x31+ I*x32,x33 +I*x34  ;  x41+ I*x42,x43 +I*x44;  \
  x51+ I*x52,x53 +I*x54  ; x61+ I*x62,x63 +I*x64];

zz = [I*zz2,zz3+I*zz4];

\\aa = qmd( qma( qreta[3], qi), qma(qreta[2],qj )); aainv = qminv(aa);  eps = aa[7,7];
\\qmd(qTrans(-XX, -zz),qmd(aa,qmd(qTrans(XX,zz), aainv))) == qTrans(qsdot(XX,qconj(eps)-[1,0]), qdot(eps,qdot(zz,qconj(eps))) - zz\
\\- qimag(qipd(qsdot(XX,qconj(eps)),XX  ))  )

\\ the above formula is to check that for every $l $ in Leech, a = matrix of phi( qreta[1],alpha) phi(qreta[2], i*alpha*i^-1), we have
\\T( l,z)^-1 * a * T(l,z) * a^-1 =  T(l(ebar - 1), ebar*z*e - z + Im<l ebar, - l>) where e = -i*alpha = (-+-+); a= (e,0;*,d) with d*pbar^-1=pbar^-1 *e
\\but the left hand side of the equation is in $R(L)$, so is the right hand side. So R(L) has a translation in each l in Leech, and this allows
\\us to get all translations from finitely many roots and then all height one roots from them.

qwl = [0, 0; 1 - I, 0; 0, 0; 1 - I, 0; 0, 0; 1 - I, 0; -2 - I, 0; -1 + I, 0];
qwp = [1, 0; 1, 0; 1, 0; 1, 0; 1, 0; 1, 0; -2*I, 0; -2, 0];
\\These are two vectors fixed under diagram automorphisms L2(F2),
\\qwl = point*p      + Sum of the lines on that point = perp to the lines, norm 2
\\qwp = line*conj(P) + Sum of the points on that line = perp to the points, norm 2


\\****************************************************************************************************************************************
\\reduce an quaternion modulo p= ( 1 - I)

\\cosets generated by 0,1, ( 1+ i+ j+ k)/2, and ( -1 + i + j+ k)/2
modpold(v)=
{
local(an,bn,cn,dn);
cn = real(v[2]) + floor(imag(v[2]));
dn = imag(v[2]) - floor(imag(v[2])); 
bn = imag(v[1]) + floor(cn);
cn = cn - floor(cn); 
an = real(v[1]) + floor(bn) ; 
bn = bn - floor(bn);
an = divrem(2*an, 4)[2]/2; if(an == 3/2, an = -1/2);
return([an + bn*I, cn + dn*I]);
}

\\cosets generated by 0, 1, ( - 1 - i - j - k)/ 2 and ( - 1 + i + j + k)/2 as in wilson's paper

modp(v)=
{
local(summ, modpv);
summ = real(v[1]) + imag(v[1]) + real(v[2]) + imag(v[2]);
if( frac(real(v[1])) == 0, modpv = [0,0] ; if( divrem(summ,2)[2] == 1, modpv = [1,0]));
if( frac(real(v[1])) != 0, modpv = -al ; if( divrem(summ,2)[2] == 1, modpv = [ (-1 + I)/2,( 1 + I)/2]) );
return(modpv);
}

\\reduce a vector over quaternions modulo p

modpvec(v)=
{
local(i, modpvv);
modpvv = matrix( matsize(v)[1], matsize(v)[2]);
for( i = 1, matsize(v)[1], modpvv[i,] = modp( v[i,]) );
return(modpvv);
}

\\change a lx2 matrix  v thought of as a Hurwitz column vector  to a 4*lx1 integer column vector and conversely 

qmtov(v)=
{
local(l,i,j,vect);
l = matsize(v)[1];
vect = vector(4*l)~;
for( i = 1, l , for(  j = 1, 2, vect[4*(i-1)+2*j -1 ] = real(v[i,j]) ;vect[4*(i-1)+2*j ] = imag(v[i,j])  ));
return(vect);
} 

qvtom(v)=
{
local(l,i,j,vect);
l = matsize(v)[1]/4;
vect = matrix(l,2);
for( i = 1, l , for(  j = 1, 2, vect[i,j] =v[4*(i-1)+2*j-1] + I*v[4*(i-1)+2*j]    ));
return(vect);
} 

\\generate the hexacode over GF(4)
genhex()=
{
local(wd,wdd,i,j,k,ii,cou);
wd = vector(6); wdd = vector(3); cou = 0; cword = matrix(6,2);
wd[1] = [  0   ,  0   ;   0  ,  0   ;  0   ,  0   ;  0   ,  0   ;  0   ,  0   ;  0   ,  0   ];
wd[2] = [  1   ,  0   ;   1  ,  0   ;  1   ,  0   ;  1   ,  0   ;  1   ,  0   ;  1   ,  0   ];
wd[3] = [  1   ,  0   ;   1  ,  0   ;om[1] ,om[2] ;Om[1] ,Om[2] ;Om[1] ,Om[2] ;om[1] ,om[2] ];
wd[4] = [  0   ,  0   ;   0  ,  0   ;Om[1] ,Om[2] ;om[1] ,om[2] ;om[1] ,om[2] ;Om[1] ,Om[2] ];
wd[5] = [om[1] ,om[2] ;om[1] ,om[2] ;  0   ,  0   ;Om[1] ,Om[2] ;Om[1] ,Om[2] ;  0   ,  0   ];
wd[6] = [Om[1] ,Om[2] ;Om[1] ,Om[2] ;om[1] ,om[2] ;  0   ,  0   ;  0   ,  0   ;om[1] ,om[2] ];
cou = cou + 1; print( cou, "  " , "hex[", cou, "]=",wd[1],";" );
wdd[1] = wdd[2] = wdd[3] = wd[2];
for(ii = 1,6, wdd[2][ii,] = modp( qdot(wdd[2][ii,],om));wdd[3][ii,] = modp( qdot(wdd[3][ii,],Om)) );
for( j = 1, 3, cou = cou + 1 ; print( cou, "  " , "hex[", cou, "]=",  wdd[j],";" ));
for( i =3, 6, wdd[1] = wdd[2] = wdd[3] = wd[i] ;\
for(ii = 1,6, wdd[2][ii,] = modp( qdot(wdd[2][ii,],om));wdd[3][ii,] = modp( qdot(wdd[3][ii,],Om)) );\
for( j = 1, 3, for(k = 0, 4, cword[1,] = wdd[j][1,] ; \
for( ii = 0, 4, cword[divrem(ii + k,5)[2] + 2,] = wdd[j][ii+2,] );\
cou = cou+1; print(cou, "  " , "hex[", cou, "]=",cword,";" )   )));
return();
}


L25()=
{
local(g,i,j,k,cosetg);
g =vector(60);
for( i = 1, 60, g[i] = matrix(6,6); for( j = 1, 6, for( k = 1, 6, g[i][j,k] = [0,0])));
g[1] = Mid(6);
g[2][1,1] = g[2][2,6] =  g[2][3,2] =  g[2][4,3] =  g[2][5,4] =  g[2][6,5] = [1,0]; 
for(j = 3,5, g[j] = qmd(g[j-1], g[2]));
g[6][1,1] = g[6][2,2] =  g[6][3,6] =  g[6][6,3] =  g[6][4,5] =  g[6][5,4] = [1,0]; 
for(j = 7, 10, g[j] = qmd( g[j-5],g[6]));
cosetg = vector(5);
for( i = 1, 5, cosetg[i] = matrix(6,6); for( j = 1, 6, for( k = 1, 6, cosetg[i][j,k] = [0,0])) );
cosetg[1][1,2] = cosetg[1][2,1] =cosetg[1][3,6] =cosetg[1][6,3] =cosetg[1][4,4] =cosetg[1][5,5] = [1,0];
for( i = 2, 5, cosetg[i] =qmd( g[i], cosetg[1]) );
for( i = 1, 5, for( j = 1, 10, g[10*i + j] = qmd( cosetg[i], g[j] ) ));
return(g);
}

writel25()=
{
local(tr,ee,i,j,k);
tr = ee= vector(6);for( i = 1, 6, ee[i] = matrix(6,2); ee[i][i,] = [1,0] );
for( i = 1, 60, for( j = 1, 6, for( k =1 , 6, if( qmv(g[i],ee[j]) == ee[k], tr[j] = k ))); print([1,2,3,4,5,6]);print(tr); print() );
}

\\basis of Leech from page 453 of Wilson

bb = vector(6);
bb[1] = [2+2*I,  0  ;  0  ,  0  ;  0  ,   0  ;  0  ,  0  ;   0  ,  0  ;  0  ,  0  ];    
bb[2] = [  2  ,  0  ;  2  ,  0  ;  0  ,   0  ;  0  ,  0  ;   0  ,  0  ;  0  ,  0  ];    
bb[3] = [  0  ,  0  ;  2  ,  0  ;  2  ,   0  ;  0  ,  0  ;   0  ,  0  ;  0  ,  0  ];    
bb[4] = [  I  ,1+I  ;  1  ,  0  ;  1  ,   0  ;  1  ,  0  ;   1  ,  0  ;  1  ,  0  ];    
bb[5] = [  0  ,  0  ;  0  ,  0  ;  1  ,   I  ;  1  ,  1  ;   1  ,  1  ;  1  ,  I  ];    
bb[6] = [  0  ,  0  ;  1  ,  1  ;  1  ,   1  ;  1  ,  I  ;   0  ,  0  ;  1  ,  I  ];    

\\an hash function for norm -4 vectors of Leech

qhf(x)=
{
local(hhf, i);
hff = 1;
x= qsdot(x,[2,0]); 
for( i = 1, 6, hhf = hhf + divrem( real(x[i,1]) + 2*i , 3)[2]*3^(i-1));
for( i = 1, 6, hhf = hhf + divrem( imag(x[i,1]) + 5*i , 2)[2]*2^(i-1));
for( i = 1, 5, hhf = hhf + divrem( real(x[i,2])+imag(x[i,2]) + 3*i , 5)[2]*5^(i-1));
return(hhf);
}



\\starts with a list of m norm -4 vectors r[1], ..,r[m] in Leech , and generates more by adding them in each other iteratively
\\uses the above hash function

genql2()=
{
local(m,ttime,R,new,i,hh,ind,has, indmax,n,iold,jold,j,jj,ii,iii,k,flag,hashk,ij,newij);
m=6;ttime = 1024000;R=vector(ttime+1);new=matrix(6,2);has = matrix(4000,800);ind = vector(4000); 
for (i = 1 , m, R[i] = bb[i];hh= qhf(R[i]);ind[hh]++;has[hh,ind[hh]]=R[i];write(ql2list, " i[",i,"]=",i,"; j[", i,"]=",i, ";ql2[",i,"]=",R[i],";" )); 
indmax=1; i=1; j=2; n = m; iu = 1;
while ( n < ttime, new = qvad(R[i], qsdot(R[j],quni[iu]) ); iu++  ;iold = i;jold = j;jj=3;\
if(iu==25, iu = 1;\
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;));\
if( qipd(new,new)[1] == -8,for( ij = 1, 24, k=1;  flag=0; newij = qsdot(new,quni[ij]);hh = qhf(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(ql2list, " i[",n,"]=",iold,"; j[", n,"]=",jold, "; ql2[",n,"]=",newij,";")) ) );
return();
}

\\list the pair of vectors in first shell of leech lattice that can give a E_8
liste8()=
{
local(ll,i,j,k,vv,cou);
ll = matsize(ql2)[2];cou = 0;
for( i = 1, ll, for( j = i + 1, ll, for( k = 1, 24, vv = qvad(ql2[i], -qsdot(ql2[j],quni[k]));\
if( qipd(vv,vv)[1] == -12, cou = cou + 1;indd = [i,j,k];write(e8ist,"ind[",cou,"]=", indd,";" ); break() ))));
return();
}

\\given an E8 generated by x1 & x2 find all the first shell vectors in Leech that may belong to a complementary E8
\\store the output in the variable "dis" to be used by the next code.

dist4e8(x1,x2)=
{
local(cou, i, j, vv, diff1,diff2,diss,dis);
cou = 0;
diss = vector(1000);
for(i = 1, 7403, for(j = 1, 24, vv = qsdot(ql2[i],quni[j]) ; diff1 = qvad(x1, - vv); diff2= qvad(x2, -vv); \
imcycle=(qimag(qipd(x1,x2)) -qimag(qipd(x1,vv)) + qimag(qipd(x2, vv)))/2 ; \
if( qipd(diff1,diff1)[1] == -8 && qipd(diff2, diff2)[1] == -8 && imcycle == [I,0] ,\
cou = cou+1; diss[cou] = vv  ) ));
dis = vector(cou); for( i = 1, cou , dis[i] = diss[i] );
return(dis);
}

\\look through the list "dis" found by the previous code, and find the pairs that make an E8

findE82(x1,x2)=
{
local(qgg,i,j,flag,xx,diff,ii,jj);
qgg = matrix(4,4);
xx = vector(4); xx[1] = x1; xx[2] = x2 ;
ll = matsize(dis)[2]; 
for( ii = 1, ll, for( jj = ii+1, ll,\
xx[3] = dis[ii]; xx[4]= dis[jj];\
imcycle= -qimag(qipd(xx[1],xx[2]))+ qimag(qipd(xx[3],xx[4]))- qimag(qipd(xx[2],xx[4]))+ qimag(qipd(xx[1],xx[3])) ;\
diff = qvad( xx[3], -xx[4]);\
if( qipd(diff, diff)[1] == -12 && imcycle == [0,0] , print( ii, "  ", jj) ) ));
return();
}

\\FINDINGS: 
\\had a big list of first shell vectors "ql2" compiled by iterating from the basis using a hash function 
\\Started with del[1] and del[2] which can make an E8 ( found by hand).
\\dis = dist4e8(del[1], del[2] ) found 20 first shell vectors from the list "ql2" that may be a part of a complementary E8
\\findE82(del[1], del[2] ) finds only 4 pairs out of the 20 vectors in "dis" that form complementary E8 pairs.
\\two out of these give orthogonal E8 We call them del[3], del[4] and del[5], del[6] respectively.
\\the Lorentzian roots now have the form (delta[i] ; 1, qinv(qconj(p))*( 1 + beta[i] ));
\\beta[1] is taken to be I and the others are then determined.

del = vector(6);
del[1] =[2, 0; 2, 0; 0, 0; 0, 0; 0, 0; 0, 0];
del[2] =[I, 1 + I; 1, 0; 1, 0; 1, 0; 1, 0; 1, 0];
del[3] =[1/2 + 1/2*I, 1/2 + 1/2*I; 3/2 + 1/2*I, -1/2 + 1/2*I; 1/2 - 1/2*I, 1/2 - 1/2*I; 1/2 - 1/2*I, -1/2 + 1/2*I; -1/2 - 1/2*I, -1/2 - 1/2*I; 1/2 - 1/2*I,-1/2 + 1/2*I];
del[4] =[1 + I, I; 1, 0; 0, I; 0, -I; 1, 0; 0, -1];
del[5] = [1 + I, I; 1, 0; I, 0; 0, -1; -I, 0; 1, 0];
del[6] = [1 + I, 0; 1, I; 1 - I, 0; 1, -I; 0, 0; 0, 0];

beta = vector(6);
beta[1] = [I,0]; 
beta[2] = beta[1]- qimag( qipd(del[1], del[2]))/2 + [ I,0] ; 
beta[3] = beta[1]- qimag( qipd(del[1], del[3]))/2;
beta[4] = beta[1]- qimag( qipd(del[1], del[4]))/2;
beta[5] = beta[1]- qimag( qipd(del[1], del[5]))/2;
beta[6] = beta[1]- qimag( qipd(del[1], del[6]))/2;


\\the 8 vectors rr form a basis for Leech + H where the first 6 make 3 copies of E_8 rr[i] = ( del[i]; 1, beta[i] )
\\rr[7] and rr[8] is a hyperbolic cell orthogonal to the three E8's found by using qeperp and fiddling with them
\\these vectors are to be used to make an isomorphism from 3E8 + H to Leech + H

rr = vector(8);
for( i = 1, 6, rr[i] = matrix(8,2); for( j = 1,6, rr[i][j,]=del[i][j,]);rr[i][7,] = [1,0]; rr[i][8,] = qdot(p, [1,0] + beta[i])/2 ); 
rr[7] = [-7/2 + 1/2*I, -3/2 - 5/2*I; -7/2 + 3/2*I, -3/2 + 1/2*I; -1/2 + 3/2*I, -1/2 + 1/2*I; -1/2 + 1/2*I, -1/2 + 3/2*I; -1/2 + 3/2*I, 1/2 + 3/2*I; -3/2 + 1/2*I, -1/2 + 1/2*I; -5/2 + 3/2*I, -3/2 + 1/2*I; -4 + 2*I, -3 - I];
rr[8] = [-5/2 - 3/2*I, -1/2 - 5/2*I; -7/2 + 1/2*I, -1/2 - 3/2*I; -1/2 + 1/2*I, -1/2 - 1/2*I; -3/2 + 1/2*I, -1/2 + 1/2*I; -1/2 + 3/2*I, -1/2 + 1/2*I; -3/2 + 1/2*I, 1/2 - 1/2*I; -5/2 + 1/2*I, -1/2 - 1/2*I; -4, -2 - 2*I];



\\the possible pairs of dis that work for the 2nd E8. we have taken the first two, dis[4] & dis[6]
\\[1/2 + 1/2*I, 1/2 + 1/2*I; 3/2 + 1/2*I, -1/2 + 1/2*I; 1/2 - 1/2*I, 1/2 - 1/2*I; 1/2 - 1/2*I, -1/2 + 1/2*I; -1/2 - 1/2*I, -1/2 - 1/2*I; 1/2 - 1/2*I, -1/2 + 1/2*I];
\\[1 + I, I; 1, 0; 0, I; 0, -I; 1, 0; 0, -1];
\\[1 + I, I; 1, 0; 0, -I; 0, I; 1, 0; 0, -1];
\\[3/2 + 1/2*I, 1/2 + 1/2*I; 1/2 + 1/2*I, 1/2 - 1/2*I; 1/2 - 1/2*I, 1/2 + 1/2*I; 1/2 - 1/2*I, -1/2 - 1/2*I; 1/2 - 1/2*I, 1/2 + 1/2*I; 1/2 + 1/2*I, 1/2 - 1/2*I];
\\[1 + I, I; 1, 0; I, 0; 0, -1; -I, 0; 1, 0];
\\[1 + I, 0; 1, I; 1 - I, 0; 1, -I; 0, 0; 0, 0];
\\[3/2 + 1/2*I, 1/2 + 1/2*I; 1/2 + 1/2*I, -1/2 + 1/2*I; 1/2 - 1/2*I, 1/2 + 1/2*I; 1/2 - 1/2*I, 1/2 + 1/2*I; 1/2 + 1/2*I, -1/2 + 1/2*I; 1/2 - 1/2*I, -1/2 - 1/2*I];
\\[1/2 + 1/2*I, 1/2 + 1/2*I; 3/2 + 1/2*I, 1/2 - 1/2*I; -1/2 - 1/2*I, -1/2 - 1/2*I; 1/2 - 1/2*I, -1/2 + 1/2*I; 1/2 - 1/2*I, 1/2 - 1/2*I; 1/2 - 1/2*I, 1/2 - 1/2*I];

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

qmultipd(x)=
{
local(i,j,l,mul, mulv, testv,ipp);
l = matsize(x)[1];
mul  = matrix(4,4*l); mulv = vector(4);
for( i = 1, l, \
testv = matrix(l,2);testv[i,1] =   1 ; ipp = qipd(x,testv); mul[,4*i-3]=[real(ipp[1]),imag(ipp[1]),real(ipp[2]),imag(ipp[2])]~; \
testv = matrix(l,2);testv[i,1] =   I ; ipp = qipd(x,testv); mul[,4*i-2]=[real(ipp[1]),imag(ipp[1]),real(ipp[2]),imag(ipp[2])]~; \
testv = matrix(l,2);testv[i,2] =   1 ; ipp = qipd(x,testv); mul[,4*i-1]=[real(ipp[1]),imag(ipp[1]),real(ipp[2]),imag(ipp[2])]~; \
testv = matrix(l,2);testv[i,2] =   I ; ipp = qipd(x,testv); mul[,4*i  ]=[real(ipp[1]),imag(ipp[1]),real(ipp[2]),imag(ipp[2])]~);
for( i = 1,4,mulv[i] = mul[i,]);
return(mulv);
}

\\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 

qeperp(e)=
{
local(ht, l,i,j,k,Mb, cfeperpz,cfeperp,eperpp,Mbz, Mbit,BB);
ht = matsize(e)[2];
l = 8;
Mb = matrix( ht,l ); 
BB = vector(8);
for( i = 1, 8, BB[i] = matrix(8,2));
for( i = 1, 6, for(j = 1, 6, BB[i][j,] = bb[i][j,])); BB[7][7,1] = BB[8][8,1] = 1;
for( i = 1, ht, for( j = 1, l, Mb[i,j] = qipw(e[i],BB[j]) ));
Mbz = matrix(4*ht, 4*l);
 for( i = 1, ht, \
 Mbit = matrix(l,2); for( j = 1, l, Mbit[j,] = qconj(Mb[i,j]));\
 for( k = 1, 4, Mbz[4*(i-1) + k, ] = qmultipd(Mbit)[k]));
print("Mbz"); print(matsize(Mbz), "  ",matrank(Mbz));
cfeperpz = matkerint( Mbz); 
cfeperp = eperpp = vector(l - ht);
for( j = 1,  l - ht, cfeperp[j] = qvtom( cfeperpz[,4*j-3] ) );  
for( j = 1, l - ht, eperpp[j] = matrix(l,2);\
 for( i = 1, l, eperpp[j] = qvad( eperpp[j] , qsdot( BB[i],cfeperp[j][i,] ) ) ) );   
return( eperpp);
}

\\check if the first 6 coordinates of a vector give a leech lattice vector

checkqleech(v)=
{
local();
print("1st condition  ", modp(v[4,] - v[5,]) , "  ", modp(v[4,] - v[6,] ) );
print("? zero mod 2  ", qdot(v[3,] + v[6,], Om) + qdot(v[4,] + v[5,], om) , "  ", qdot(v[2,] + v[3,],om) + qdot(v[4,] + v[6,],Om) );
print("? zero mod 2+2i  ", modp((qdot(-v[1,],[I,1+I])+v[2,] + v[3,] +v[4,] + v[5,] + v[6,])/2) );
return();
} 

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

\\generators for the diagram automorphisms 8.P^2(F_2) acting on 3E8 + H

\\some computation on P^2(F_2)

qp = qP = vector(14);
\\points:		        lines

qp[ 1] =  Mod([ 1, 1, 1],2)~; qp[ 8] = qp[ 1]~ ;  qP[ 1] =AA  ; qP[ 8] =FF  ;
qp[ 2] =  Mod([ 0, 1, 1],2)~; qp[ 9] = qp[ 2]~ ;  qP[ 2] =C1  ; qP[ 9] =B1  ;
qp[ 3] =  Mod([ 1, 0, 1],2)~; qp[10] = qp[ 3]~ ;  qP[ 3] =C2  ; qP[10] =B2  ;            
qp[ 4] =  Mod([ 1, 1, 0],2)~; qp[11] = qp[ 4]~ ;  qP[ 4] =C3  ; qP[11] =B3  ;       
qp[ 5] =  Mod([ 1, 0, 0],2)~; qp[12] = qp[ 5]~ ;  qP[ 5] =E1  ; qP[12] =D1  ;        
qp[ 6] =  Mod([ 0, 1, 0],2)~; qp[13] = qp[ 6]~ ;  qP[ 6] =E2  ; qP[13] =D2  ;             
qp[ 7] =  Mod([ 0, 0, 1],2)~; qp[14] = qp[ 7]~ ;  qP[ 7] =E3  ; qP[14] =D3  ;              

\\act by a linear transformation T on P^2(F_2);

\\T has in entered as a 3x3 matrix mod 2.

\\compute the action of T on the indices 1,2,...14.the vector trp(T) stores the answer

qtrp(T) = 
{
local(i,j,tp,pp,ll);
tp=vector(14);
for ( i = 1, 7, \
  pp = T*qp[i]; ll = qp[7+i]*T^-1; \
   for(j = 1, 7,\
    if( pp == qp[j] || pp == -qp[j], tp[i]=j);\
    if( ll == qp[7+j] || ll == -qp[7+j], tp[7+i]=7+j);\
     )\
   );
return(tp);
}

\\ The permutation of points and lines by T

qperm(T)=
{
local(i,j,tp);
tp=qtrp(T);
print("points  :");
for( i = 1, 7, print1(qP[i],"  "));
print(" ");
for( i = 1, 7, print1(qP[tp[i]],"  "));
print(" ");
print("lines  :");
for( i = 1, 7, print1(qP[7+i],"  "));
print(" ");
for( i = 1, 7, print1(qP[tp[7+i]],"  "));
return();
}

\\take the generator matrix [a,c_i, e_i, f] of H^(1,7) and transform by T and store the anwer in tgen 
\\check if the action of G = PGl_3(F_2) on the vectors of length 14 are lineear.
\\See if for the element T in G the action on the basis matches the action on the other 6 lines 

qgg = matrix (8,8);
for ( j = 1, 8, for( i = 1, 8,  qgg[i,j] = qr[j][i,]));
qgginv=qminv(qgg);

qtrg(T)=
{
local(i,j,tp,tg);
tp=qtrp(T);
tg = matrix(8,8);
for ( j = 1, 8, for( i = 1, 8,  tg[i,j] = qr[tp[j]][i,])); 
\\for ( i = 9, 14, newv = qr[tp[i]] - qmv(tg ,qmv( qgginv , qr[i] )) ; print( qip(newv, newv))    );
return(tg);
}

\\the two standard generators of P^2(F_3)

qTa = Mod( [1,0,0; 0,1,0; 0,1,1],2);
qTb = Mod( [0,1,0; 0,0,1; 1,0,0],2);

\\the 14x14 matrices for the two standard generators in its representation on C^(1,13)

qMa = qmd(qtrg(qTa) , qgginv );
qMb = qmd(qtrg(qTb) , qgginv );

\\check that the action of Ma and Mb on the vectors r[i] is same as the permutation action of a & b on L3F3
\\ for( i = 1, 14, print( qmv( Ma, qr[i]) == qr[trp(Ta)[i]], "  ", qmv( Mb, qr[i]) == qr[trp(Tb)[i]]) );

\\XXX = [x11 + I*x12, x13 + I*x14; x21 + I*x22, x23 + I*x24; x31 + I*x32, x33 + I*x34; x41 + I*x42, x43 + I*x44;\
\\ x51 + I*x52, x53 + I*x54; x61 + I*x62, x63 + I*x64; x71 + I*x72, x73 + I*x74; x81 + I*x82, x83 + I*x84];

presL2F2()=
{
local(qc,qd);
qc = qmd(qMa, qMb); qd = qmd( qc, qmd( qminv( qMa), qminv(qMb)));
print( "order of a    = ",qord(qMa));
print( "order of b    = ",qord(qMb));
print( "order of ab   = ",qord(qc));
print( "order of [a,b]= ",qord(qd));
return();
}

\\the automorphism of order 8 interchaning points and lines: (a, f), (c_i, b_i), (e_i, f_i), (line) to (point), (point) to (line)(i)

genqsig()=
{
local(ta, tb, i, j);
ta = tb = matrix(8,8);
for( j = 1, 7, linei = qsdot( qr[7+j], qi); for( i = 1, 8, ta[i,j] = qr[j][i,] ; tb[i,j] = linei[i,] ));
for( i = 1, 8, ta[i,8] = qr[8][i,] ; tb[i,8] = qr[1][i,] );
return( qmd( tb, qminv( ta)) );
}

\\check that it does what it is supposed to do
\\for( i =1, 7, print( qmv( qsig, qr[i] ) == qsdot(qr[7+i], qi), "  ", qmv(qsig, qr[7+i]) == qr[i]  ));

\\the height of a root defined as its inner product with the fixed vector wp + wl( 1 + i)/sqrt(2)

qht(x) = real(qn( qip( qwp, x) + qdot( p, qip( qwl, x))/sqrt(2) )[1]);

\\generators for the reflection group in coordinates Leech + H
\\use the Z- basis for Leech given by bb[i] multiplied by, 1, i, j , om for i = 1, ..., 6. Note 1,i,j,on is an additive basis for hurwitz int

genrefw()=
{
local(i,zbasis,rr1,rr2);
zbasis = vector(24); genw = vector(81);
for( i = 1, 6, zbasis[4*i - 3] = bb[i]; zbasis[4*i - 2] = qsdot(bb[i],qi); zbasis[4*i - 1] = qsdot(bb[i], qj) ; zbasis[4*i]=qsdot(bb[i],om));
rr1 = [0,0; 0,0; 0,0; 0,0; 0,0; 0,0; 1,0; I,0];
rr2 = [0,0; 0,0; 0,0; 0,0; 0,0; 0,0; 1,0;-1,0];
rr3 = [0,0; 0,0; 0,0; 0,0; 0,0; 0,0; 1,0;(-1+I)/2,(1-I)/2];
zeroo = [0,0; 0,0; 0,0; 0,0; 0,0; 0,0];
for( i = 1, 24,\
genw[3*i - 2] = qmv( qTransw( zbasis[i],[0,0]),rr1);\
genw[3*i - 1] = qmv( qTransw( zbasis[i],[0,0]),rr2 );\
genw[3*i    ] = qmv( qTransw( zbasis[i],[0,0]),rr3 ) );\
genw[73] = qmv(qTransw( zeroo,[I,1]),rr1 ); genw[74] = qmv( qTransw(zeroo,[I,1]),rr2);genw[75] = qmv( qTransw(zeroo,[I,1]),rr3);
genw[76] = qmv(qTransw( zeroo,[I,I]),rr1 ); genw[77] = qmv( qTransw(zeroo,[I,I]),rr2);genw[78] = qmv( qTransw(zeroo,[I,I]),rr3);
genw[79] = rr1; genw[80] = rr2;genw[81] = rr3; 
return(genw);
}

XXX= [x11 + I*x12, x13 + I*x14; x21 + I*x22, x23 + I*x24; x31 + I*x32, x33 + I*x34; x41 + I*x42, x43 + I*x44;\
 x51 + I*x52, x53 + I*x54; x61 + I*x62, x63 + I*x64; x71 + I*x72, x73 + I*x74; x81 + I*x82, x83 + I*x84];


\\the matrix that qCC takes 3E8 + H to Leech + H

rre = vector(8);
rre[1] = qd1; rre[2] = qc1;rre[3] = qd2; rre[4] = qc2;rre[5] = qd3; rre[6] = qc3;
rre[7]=[0,0;0,0;0,0;0,0;0,0;0,0;1,0;0,0];
rre[8]=[0,0;0,0;0,0;0,0;0,0;0,0;0,0;1,0];

qe8s1 = qe8s2 = matrix(8,8);

for( i = 1, 8, for( j = 1, 8,qe8s2[i,j] = rr[j][i,]  ; qe8s1[i,j] = rre[j][i,] ));

qCC = qmd( qe8s2, qminv(qe8s1));

genrefe()=
{
local(i,gene,qCCi);
gene = vector(81);
qCCi = qminv(qCC);
Genw = genrefw;
for( i = 1, 81, gene[i] = qmv(qCCi, Genw[i] ) );
return(gene);
}

qGe = genrefe;

\\qGe is a set of 81 vectors in 3E8 + H that generate the reflection group

qdecreasehtnew(x)=
{
local(absht, tr,i,flag,cou,track,trackk, ii,jj);
track = vector(1000);
while( qht(x) > 2 + .00000001, cou++; absht = qht(x);flag = 0; i = 1; jj = 1;
while( i <= 14 && flag ==0 , tr = qphi( qr[i], x, qunii[jj]); if( qht(tr) < absht -.00000001, x = tr; track[cou] = [i,jj]; \
print(qht(x)," ",x); flag = 1  )  ;if( jj < 8, jj++, i++; jj=1)  ); if( flag == 0, return(x)) );
return();
}

\\perturbing is needed in the following cases and always reflecting in qGe[2] works - either i or j reflection

qperturb = vector(81); 
qperturb[1 ] =[2,qi];qperturb[8 ] =[2,qj];qperturb[37] =[2,qi];qperturb[44] =[2,qi];qperturb[50] =[2,qi];qperturb[51] =[2,qj];
qperturb[54] =[2,qi];qperturb[58] =[2,qi];qperturb[60] =[2,qj];qperturb[67] =[2,qj];qperturb[74] =[2,qi];qperturb[81] =[2,qi];



checkq3e8(v)=
{
local( i, v71,v72,v73,v74,v81,v82,v83,v84);
print( modp( v[1,] - v[2,] )," ", modp(v[3,] - v[4,] )," ",modp(v[5,] - v[6,]) );
v71 = 2*real(v[7,1]); v72 = 2*imag( v[7,1]); v73 = 2*real( v[7,2]); v74 = 2*imag( v[7,2]); 
v81 = 2*real(v[8,1]); v82 = 2*imag( v[8,1]); v83 = 2*real( v[8,2]); v84 = 2*imag( v[8,2]); 
print( divrem(v71 - v72, 2)[2],"  ",divrem(v71 - v73, 2)[2],"  ",divrem(v71 - v74, 2)[2],"  ",\
divrem(v81 - v82, 2)[2],"  ",divrem(v81 - v83, 2)[2],"  ",divrem(v81 - v84, 2)[2] );
return();
}

\\  \begin{equation*}
\\  \begin{pmatrix}
\\  2           & 2          & 0         & 0         & 0           & 0       & 1              & 1         \\
\\  i+j+k      & 1          & 1         & 1         & 1           & 1       & 1              & 2+i+j     \\
\\  \frac{1+i+j+k}{2}    & \frac{3+i-j+k}{2}    &  \frac{1-i+j-k}{2}  & \frac{1 -i-j+k}{2}  & \frac{-1-i-j-k}{2}    & \frac{1-i-j+k}{2} & 
\\    1 & \frac{3+i+j+k }{2} \\
\\  i+i+k       & 1          & k         & -k        &  1          & -j      & 1              & \frac{3+i+j+k}{2}\\
\\  1+i+k          & 1              & i             & -j            & -i             & 1           & 1              & \frac{3+i+j+k}{2}\\
\\  1+i            & 1+k            & 1-i           & 1-k           & 0              & 0           & 1              & \frac{3+i+j+k}{2}\\
\\  \frac{-7+i-3j-5k}{2} & \frac{-7+3i-3j+k}{2} & \frac{-1+3i-j+k}{2} & \frac{-1+i-j+3k}{2} & \frac{-1+3i+j+3k}{}2 & \frac{-3+i-j+k}{2} & 
\\    \frac{-5+3i-3j+k}{2} & -4+2i-3j-k \\
\\  \frac{-5-3i-j-5k}{2} & \frac{-7+i-j-3k}{2}  & \frac{-1+i-j-k}{2}  & \frac{-3 +i -j+k}{2}& \frac{-1+3i-j+k}{2}  & \frac{-3+i+j-k}{2} &
\\    \frac{-5+i-j-k}{2}   & -4-2j-2k \\ 
\\  \end{pmatrix}
\\  \end{equation*}
\\  \end{tiny}


textform(z)= [real(z[1]) , imag(z[1])*ii , real(z[2])*jj   ,imag(z[2])*kk ] ;

\\ product of the reflections in the unit un in a sequence of vectors given by the enties of x
qpr(x,un)=
{
local(lm,prr);
lm=matsize(x)[2];
prr= qma(x[1],un);for(i=2,lm, prr = qmd(prr,qma(x[i],un) ));
return(prr);
}

\\qmv( qpr([qd1,qc1,qb1,qa,qb2,qc2],qi), qd2 ) == -qe3