program mod2matrices;
{ procedures to work with square matrices with entries in Z/2^n }
{ where 0 <= n <= 7 }
label 9999;
const

  (* to change size of matrix, change both these *)
  nsize =  3;  (* # of rows *)
  nsize2 = 6;  (* twice number of cols for row reduction calc *)

  (* arith. is done mod 2^exp2. conv. to have other forms *)
  exp2 =   1;
  exp21 =  1; (* power for computing inverses *)
  pow21 =  $1; (* mask for adjusting arith. operations *)


type
   colvec =  array[1..nsize] of integer;
   rowvec =  array[1..nsize] of integer;
   rowvec2 = array[1..nsize2] of integer;
   matrix =  array[1..nsize,1..nsize] of integer;
   doublemat = array[1..nsize,1..nsize2] of integer; (* for row reductions *)

 var
    tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp : matrix;
    tua,tub,tuc,tud,tue,tuf : matrix;  (* signed transpositions *)
    ua,ub,uc,ud,ue,uf : matrix;        (* transpositions *)
    e1,e2 : matrix; { generate abelian subgroup Z2 x Z2 in upper row }
    tau : matrix; { element of order two that switches e1 and e2 }
    A,B,c,d,g1,g2 : matrix; { generate lower GL(2,2) block }
    u,uu,uuu,t,tt,ttt,s,ss,sss,ssss,sssss,ssssss : matrix; { order 7 }
    ii,jj,kk,gg,count : integer;
    alpha,alphabar : integer;
    v1,v2,v3,v4,v5,v6,v7 : colvec;
    { roots of x^2 + x + 2 = 0, alpha == 0 mod 2}

procedure findalpha; { solve for alpha, alphabar }
var i :longint;
begin
for i:= 1 to pow21 do
  if ( ((i*i) and pow21) + i + 2) and pow21 = 0 then begin
     if (i and 1) <> 0 then alpha:= (-i-1) and pow21 else alpha:=i;
     alphabar:= (-alpha-1) and pow21;
     writeln('Alpha = ',alpha and pow21 : 6);
     writeln('Alphabar = ', alphabar and pow21 : 6);
     exit;
   end;

end;

procedure printvector ( var v : colvec);
var i : integer;
begin
  writeln;
  write('[');
  for i:=1 to nsize do write(v[i]:8);
  write(']');
  writeln;
end;

procedure printmatrix ( var A : matrix );
var i,j : integer;
begin
writeln;
for i:=1 to nsize do begin
   for j:=1 to nsize do write(A[i,j]:8);
   writeln;
end;
end;

procedure printdouble ( var A : doublemat );
var i,j: integer;
begin
writeln;
for i:=1 to nsize do begin
   for j:=1 to nsize2 do write(A[i,j]:4);
   writeln;
end;
end;

procedure smult(var a:matrix; n : integer);
var i,j : integer;
begin
   for i:=1 to nsize do for j:=1 to nsize do a[i,j]:= ( n*a[i,j]) and pow21;
end;

procedure switchrows(var A : doublemat; n,m : integer);
var temp : integer;
    k : integer;
begin
for k:=1 to nsize2 do begin
  temp:=A[n,k];
  A[n,k]:=A[m,k];
  A[m,k]:=temp;
end;
end;

procedure addrows (var A : doublemat; i,j : integer);
{ does mod 2 addition of rows }
var  k : integer;
begin
  for k:=1 to nsize2 do A[j,k]:=(A[j,k] +  A[i,k]) and 1;
end;

procedure mkscalar(var a : matrix ; n : integer );
var i : integer;
begin
  fillchar(a,sizeof(a),0);
  for i:=1 to nsize do a[i,i]:= n and pow21;
end;


procedure reduce( var a : matrix);
var i,j : integer;
begin
  for i:=1 to nsize do for j:=1 to nsize do a[i,j]:= a[i,j] and pow21;
end;

procedure addmat( var a,b,c : matrix);
var i,j : integer;
begin
  for i:=1 to nsize do for j:=1 to nsize do
     c[i,j]:= (a[i,j] + b[i,j]) and pow21;
end;

procedure transpose(var a : matrix);
begin
end;

function invert2( var a,a_inverse : matrix): boolean;
{ reduce a mod 2 and invert it }
label 1234;
var i,j,k,r : integer;
    W : doublemat;
begin
 invert2:=false;
 fillchar(w,sizeof(w),0);
 for k:=1 to nsize do w[k,nsize+k] := 1;
 (* put a into a double matrix with identity on other side *)
 for i:=1 to nsize do
   for k:= 1 to nsize do w[i,k]:=a[i,k] and 1;
 (* printmatrix(a);
 writeln;
 printdouble(w);
 readln; *)
 (* now reduce going down *)
 for j:=1 to nsize do begin  (* go across columns j *)
 (* writeln;
  printdouble(W);
  readln;  *)
     { look in col j for 1st non-zero entry at or below starting}
     for i:=j to nsize do begin
       if w[i,j] = 1 then begin
          if i <> j then begin
             switchrows(W,i,j);
           (*  printdouble(W);
             readln; *)
          end;
          k:=i;
          goto 1234; (* success *)
       end;
     end;  (* end of i loop *)
     writeln;
     writeln('Failure to invert.');
     printdouble(W);
     exit; { failed to find non-zero entry, so matrix is singular }
1234:
    for r:=k+1 to nsize do if w[r,j]= 1 then begin
         addrows(W,j,r);
      (*   printdouble(W);
         readln;   *)
    end;

end; (* bottom of j loop *)
invert2:=true; (* ok, got all 1's on diagonal *)
{ now start at bottom corner and add back }
{ not too efficient, since adds zero entries also }

for r:=nsize downto 2 do
   for k:= (r-1) downto 1 do if w[k,r] = 1 then addrows(W,r,k);
for i:= 1 to nsize do for j:=1 to nsize do a_inverse[i,j] := W[i,nsize+j];
(* printdouble(W);
readln;  *)
end; {invert 2 }

procedure matmult( var a,b,c : matrix);
{ c = a * b }
var i,j,k,sum : integer;
begin
  for i:= 1 to nsize do for j:=1 to nsize do
  begin
    sum:=0 ;
    for k:=1 to nsize do sum:=(sum + a[i,k]*b[k,j]) and pow21;
    c[i,j]:= sum;
   end;
end; { matmult }

procedure pwrsum(var P,A : matrix);
{ computes A = 1 + P +  p^2 + ..... }
  var sum,newsum,D : matrix;
     i : integer;
begin
  mkscalar(SUM,1); D:=SUM;
  for i:= 1 to exp21 do begin
      matmult(D,P,NEWSUM);
      D:=NEWSUM;
      addmat(SUM,D,NEWSUM);
      SUM:=NEWSUM;
      (* writeln;
         printmatrix(SUM);  *)
   end;
   A:=SUM;
end;

procedure invertid2( var a,a_2inv : matrix);
{ invert a matrix which is id mod 2 by
  (1 -2P)^-1 = 1 + 2p + 4p^2 +......
}
var i,n,j : integer;
    P : matrix;
begin
    for i:=1 to nsize do for j:=1 to nsize do
       if i<>j then p[i,j]:= ((-a[i,j]) and pow21 )
       else p[i,i] := (( 1 - a[i,i]) and pow21) ;
    pwrsum(P,a_2inv);    { P is 2* F }
 (*   writeln;
    printmatrix(a);
    writeln;
    printmatrix(a_2inv);    *)
end;

function invert( var a,a_inv: matrix) : boolean;
{ invert any nonsingular matrix in two steps }
{ 1) invert mod 2 with matrix c
  2) invert AC = 1-2P to get D
  3) A^-1 = CD
}

var b,c,d : matrix;
begin
  invert:=false;
  b:=a;
  if not invert2(b,c) then exit;
  { c is a mod 2 inverse to a }
  matmult(b,c,d);
  (* writeln;
  printmatrix(d); *)
  invertid2(d,b);
(*  writeln;
  printmatrix(b);
  writeln;  *)
  matmult(c,b,a_inv);
  invert:=true;
end;

procedure matconj(var A,B,C : matrix);
(* calculate C = A(B)Ainv   *)
var t,tt,ttt : matrix;
begin
  if not invert(A,t) then exit;
  matmult(A,B,tt);
  matmult(tt,t,C);
end;




procedure vectmat( var a: matrix; var z, w  : colvec );
{ w = a*z;}
var i,sum,k : integer;
begin
  for i:=1 to nsize do begin  (* across rows *)
    sum:=0;
    for k:= 1 to nsize do
        sum:= ( sum + a[i,k] * z[k] ) and pow21;
    w[i]:=  sum and pow21;
  end;
end;

procedure vectadd(var x,y,z : colvec);
{ z = x + y }
var i : integer;
begin
  for i:=1 to nsize do z[i]:= (x[i] + y[i]) and pow21;
end;


procedure detmatrix( var a: matrix; var n : integer);
begin

end;

procedure matpower( var A,B: matrix; n : integer);
{ dumb take the n-th power of the matrix }
var i : integer;
    t,tt,ttt : matrix;
begin
  t:=A;
  for i:= 2 to n do begin
     matmult(A,t,tt);
     t:=tt;
  end;
  B:=t;
end;

function comparx( var a,b : matrix) : boolean;
var i,j : integer;
begin
comparx:=false;
for i:=1 to nsize do for j:=1 to nsize do
   if a[i,j]<> b[i,j] then exit;
comparx:=true;
end;

function gorder(var a: matrix) : integer;
var t,tt,ttt : matrix;
    g : integer;
begin
 mkscalar(t,1);
 tt:=a;
 gorder:=1;
 g:=1;
 repeat
   matmult(a,tt,ttt);
   tt:=ttt; { get ready for next iteration }
   g:=g+1;
 until comparx( tt,t) or ( g > 32535);
 gorder := g;
end;



begin
findalpha;
(* generators for Z/2 x Z/2 given by upper row *)

(* want to choose a lift of e1,e2 which commute and have order 2 *)
(* and have trace -1 *)

(* ua .. uf are transpositions that generate the Sigma(4) rep
   of x_1 +...x+4 = 0 *)

(* tua.. tuf are the same, but tensored with the sign rep *)

  fillchar(ua,sizeof(ua),0);
  ua[1,1]:=(-1) and pow21;
  ua[1,2]:= 1;
  ua[2,2]:=1;
  ua[3,3]:=1;

  fillchar(ub,sizeof(ub),0);
  ub[1,2]:= (-1) and pow21;
  ub[1,3]:= (1) and pow21;
  ub[2,1]:= (-1) and pow21;
  ub[2,3]:= (1) and pow21;
  ub[3,3]:= (1) and pow21;

  fillchar(uc,sizeof(uc),0);
  uc[1,3]:= (-1) and pow21;
  uc[2,2]:= (1) and pow21;
  uc[2,1]:= (-1) and pow21;
  uc[2,3]:= (-1) and pow21;
  uc[3,1]:= (-1) and pow21;

  fillchar(ud,sizeof(ud),0);
  ud[1,1]:= (1) and pow21;
  ud[2,2]:= (-1)  and pow21;
  ud[2,1]:= (1) and pow21;
  ud[2,3]:= (1) and pow21;
  ud[3,3]:= (1) and pow21;

  fillchar(ue,sizeof(ue),0);
  ue[1,1]:= (1) and pow21;
  ue[2,1]:= (1)  and pow21;
  ue[3,1]:= (1) and pow21;
  ue[3,2]:= (-1) and pow21;
  ue[2,3]:= (-1) and pow21;

  fillchar(uf,sizeof(uf),0);
  uf[1,1]:= (1) and pow21;
  uf[2,2]:= (1)  and pow21;
  uf[3,3]:= (-1) and pow21;
  uf[3,2]:= 1;

(* tua.. tuf are the same, but tensored with the sign rep *)
  mkscalar(tmp,(-1));
  matmult(ua,tmp,tua);
  matmult(ub,tmp,tub);
  matmult(uc,tmp,tuc);
  matmult(ud,tmp,tud);
  matmult(ue,tmp,tue);
  matmult(uf,tmp,tuf);

  e1:= tua;
  e2:= tuf;
  writeln('order e1 = ',gorder(e1));
  writeln('e1 = ');
  printmatrix(e1);
  writeln;
  writeln('order e2 = ',gorder(e2));
  writeln('e2 = ');
  printmatrix(e2);
  writeln;
  readln;
  matmult(e1,e2,tmp1);
  matmult(e2,e1,tmp2);
  writeln(' e1e2 = ');
  printmatrix(tmp1);
  writeln('e2e1 = ');
  printmatrix(tmp2);
  if comparx(tmp1,tmp2) then writeln('e1 and e2 commute.')
  else writeln('e2 and e1 do not commute.');
readln;

(* make the element of order 2 which interchanges e1 and e2 *)
matmult(tub,tue,tmp2);
writeln('This element has order ',gorder(tmp2),' and switches e1 and e2 ');
writeln('conjugation.');
printmatrix(tmp2);
matmult(tmp2,e1,tmp3);
matmult(tmp3,tmp2,tmp4);
writeln('Checking conjugation of e1 by this element...');
printmatrix(tmp4);
readln;
if comparx(tmp4,e2) then writeln('Conjugation verified') else
writeln('Conjugation verification failed.');
readln;
tau:=tmp2;
mkscalar(s,0);
s[1,3]:= 1;
s[2,1]:= 1;
s[2,3]:= (-alphabar) and pow21 ;
s[3,2]:= 1;
s[3,3]:= alpha;
writeln('This matrix has order ',gorder(s));
printmatrix(s);
readln;

(* find an element of order 3 with the property *)
(* TSTT=SS *)

(* for example, in S_7, if S = (1234567), then
   T = (235)(467)
*)

(* so one brute force method is to factor S^7 -1
   = (S-1)( S^3 -alpha    S^2 + alphabar S -1)
          ( S^3 -alphabar S^2 + alpha S -1 )
 if U = (S-1)(S^3 -alphabar S^2 + alpha S -1 )
 then get a 3 dimensional subspace of the 7-dim space
 as UV
 So, if S =     0  0  1
                1  0  -alphabar
                0  1  alpha
        then Tv = [1,-alpha, alpha] in the v,Sv, S^2 v basis
        where v = [1,0,0]
        and other columns are SSTv, SSSSTv
*)
(* so figure these out *)
  v1[1]:= 1;
  v1[2]:= (2 + alpha) and pow21;
  v1[3]:= alpha and pow21;
  matpower(s,ss,2);
  writeln('Here''s the square of S ');
  printmatrix(ss);
  writeln;
  readln;
  writeln('Here''s the fourth power of S ');
  matpower(s,ssss,4);
  printmatrix(ssss);
  readln;
  (* find the other columns *)

  
  vectmat(ss,v1,v2);
  printvector(v2);
  vectmat(ssss,v2,v3);
   printvector(v3);
writeln('Here is T ');
mkscalar(t,0);
(* for ii:=1 to 3 do t[ii,1]:=v1[ii];
for ii:=1 to 3 do t[ii,2]:=v2[ii];
for ii:=1 to 3 do t[ii,3]:=v3[ii];
writeln; *)
t[1,1]:=1;
t[2,1]:= (2 + alpha) and pow21;
t[2,2]:= alpha;
t[2,3]:= ( -1 and pow21);
t[3,1]:= alpha;
t[3,2] := ( -1 and pow21);
t[3,3]:= (-1 - alpha) and pow21;

printmatrix(t);
writeln('Calculating order of t, for t = ');
printmatrix(t);
writeln;
gg:=gorder(t);
writeln('Order of t is ',gg);
matpower(t,tt,2);
printmatrix(t);
writeln( 'TSTT = ');
matmult(t,s,tmp1);
matmult(tmp1,tt, tmp2);
writeln;
printmatrix(tmp2);
readln;
if comparx(tmp2,ss) then writeln('S conjugated by T is S^2') else
writeln('Wrong conjugate of S by T.');
readln;

writeln;
writeln('Check relation on tau*S = ? ');
mkscalar(tau,0);
tau[1,1]:=1;
tau[3,2]:=1;
tau[3,3]:=1;
tau[2,3]:=1;

matmult(tau,s,c);
for ii:=1 to 5 do begin
  writeln('TAU * S to ',ii:2,'-th power = ');
  matpower(c,d,ii);
  printmatrix(d);
  writeln;
  readln;
end;

v1[1]:= 0;
v1[2]:= ( alphabar + 1) and pow21;
v1[3]:= 1  and pow21;

vectmat(T,v1,v2);
printvector(v1);
printvector(v2);

readln;
mkscalar(A,0);
A[1,1]:=3;
A[2,1]:= (2*alpha + 2 ) and pow21;A[2,2]:=(alphabar + 1) and pow21;
    A[2,3]:= (1 + alpha) and pow21;
A[3,1]:= alpha;
A[3,2]:= 1;
A[3,3]:= (-1) and pow21;
writeln('The transformation matrix is A = ');
printmatrix(A);
if invert(A,B)
then writeln('It''s inverse is Ainv = ') else
writeln('Oophs, A is not invertible.');
printmatrix(B);
readln;

{ now look for an element of order 3 that tau normalizes }
{ and which takes e1 to e1*e2, e2 to e1 }
{ can get elements of order 3 as products of pair of transpositions
  with common index }

 9999:
 end.