\\filename :ei-linear algebra.gp
\\contains the basic linear algebra operations to manipulate vectors and matrices over Eisenstein integers.
\\no reference to any inner product

\\A few special vectors

uni= vector(6); uni[1] = [ 1, 1]; uni[2] = [ 0, 1]; uni[3] = [-1, 0]; uni[4] = [-1,-1]; uni[5] = [ 0,-1]; uni[6] = [ 1, 0];

t=[1,2];w=[0,1];W =[-1,-1]; 

tt = sqrt(-3); ww = (-1+tt)/2; WW = ( -1 -tt)/2; 

\\ Eisenstein integers are written as, [a,b] = a + omega*b, where a ,b are integers.


bar(u)=[u[1]-u[2],-u[2]]; \\conjugate of a number u

dot(u,v)=[u[1]*v[1] - u[2]*v[2] , u[2]*v[1] + u[1]*v[2] - u[2]*v[2]]; \\product of two numbers u and v

n2(u) = dot( bar( u), u) ; \\norm of a number u

rea(u) = u[1] - u[2]/2; \\the real part of u

ima(u) = u[2]/2; \\the imaginary part of u

\\ if z is a complex number then om1(z) and om2(z) are defined so that z = om1(z) + omega * om2(z);

om2(z) = 2*imag(z)/sqrt(3);

om1(z) = real(z) + imag(z)/sqrt(3);

\\nearest integer of a

nint(a)=
{
local(fr, ninta); fr = frac(a); ninta = floor(a); if( frac(a) >= .5 , ninta = floor(a)+1 ); return(ninta);
}

err(a) = abs( a - nint(a)); \\the differernce between a and the nearest integer to a.

dig(z) = [ nint(om1(z)), nint(om2(z))]; \\write a complex number that is an (approx) eisentein integer in the [a,b] form

argue(x) = arg( rea(x) + sqrt(-3) * ima(x) ) ; \\the arguement of a complex number x = [a,b], -pi < arg <= pi

\\divide a number u by a number v

divid(u,v)=
{
local(denom,quo); denom = dot(bar(v),v)[1]; quo = dot( u, bar(v))/denom; return(quo);
}

\\ scalar product of a eisenstein integer and a vector in eisenstein lattice

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

\\product of a matrix by an eisenstin number

mdot( 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] = dot( a, x[i,j]) ));
return(mdott);
}

\\ scalar division of a vector by an eisenstein integer

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

\\add two vectors whose entries are eisentein integers
vad(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 eisenstein integers

md(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] + dot(a[i,j], b[j,k]) ))));
return(mdd);
} 


\\multiply a matrix with a vector to get a vector, all entries in eisenstein integer.

mv(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,] + dot(a[i,j], b[j,]) ));
return(mdd);
}


\\calculate the (p,q)th minor of the matrix A.

minor( A, p, q) = 
{
local( mino, i, j, s1,s2);
s1=matsize(A)[1]-1;
s2=matsize(A)[2]-1;
mino = matrix( s1,s2);
for ( i = 1, p-1, \
for( j = 1 ,q-1, mino[i,j] = A[i,j]);\
for( j = q, s2, mino[i,j] = A[i,j+1]));
for( i = p, s1, \
for( j = 1 ,q-1, mino[i,j] = A[i+1,j]);\
for( j = q, s2, mino[i,j] = A[i+1,j+1]));
return(mino);
}

\\calculate the inverse of a matrix by sweep out

minv(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]=divid(temp,fac);temp = B[i,k]; B[i,k]=B[j,k];B[j,k]=divid(temp,fac) );\
  for( ii = 1, s2 , if( ii != j,\
   fac = A[ii,j];\
   for( jj = 1, s2, A[ii,jj]=A[ii,jj]-dot(fac,A[j,jj]);\
   B[ii, jj]=B[ii,jj]-dot(fac,B[j,jj])\
    );\
   )));
return(B);
}

\\calculate the determinant of a matrix by upper triangularizing
mdet(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 = dot([-1,0],de) ); \
  for( ii = j+1, s2 ,\
   fac = divid(A[ii,j],A[j,j]);\
   for( jj = 1, s2, A[ii,jj]=A[ii,jj]-dot(fac,A[j,jj])) ); de = dot(de, A[j,j])\
   );
return(de);
}

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

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

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

\\reduce a vector v modulo theta

modt(v)=
{
local(l,red,test,i);
l = matsize(v)[1];
red = matrix(l,2);
for( i = 1, l, 
test = divid(v[i,]      ,t);if( frac(test[1])==0 && frac(test[2])==0, red[i,1] = 0, red[i,2] = 0 );\
test = divid(v[i,]-[1,0],t);if( frac(test[1])==0 && frac(test[2])==0, red[i,1] = 1, red[i,2] = 0 );\
test = divid(v[i,]+[1,0],t);if( frac(test[1])==0 && frac(test[2])==0, red[i,1] =-1, red[i,2] = 0 )\
);
return(red);
}

\\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 order of a matrix 

ord(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]);
  n = 1;
  y=x;
  while(y != mid, n++; y = md(x,y);if(divrem(n,50)[2]==0,print(n))  );
  
return(n);
}


\\discretize a vector a for better viewing

discretize(a) = 
{
local(da);
da = a;
for( j = 1, 6, if( a == dot(uni[j],[-3, 0]) , da = 3 );
if( a == dot( uni[j],[ 2, 1]) , da = 2 ));
if( a == [ 0, 0] , da = 0 );
return(da);
}