(* 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 *)


  initmatrix(ua);
  ua[1,1].r:=-1;
  ua[1,2].r:= 1;
  ua[2,2].r:=1;
  ua[3,3].r:=1;

  initmatrix(ub);
  ub[1,2].r:= (-1);
  ub[1,3].r:= (1) ;
  ub[2,1].r:= (-1);
  ub[2,3].r:= (1) ;
  ub[3,3].r:= (1) ;

  initmatrix(uc);
  uc[1,3].r:= (-1);
  uc[2,2].r:= (1) ;
  uc[2,1].r:= (-1);
  uc[2,3].r:= (-1);
  uc[3,1].r:= (-1);

  initmatrix(ud);
  ud[1,1].r:= (1) ;
  ud[2,2].r:= (-1);
  ud[2,1].r:= (1) ;
  ud[2,3].r:= (1) ;
  ud[3,3].r:= (1) ;

  initmatrix(ue);
  ue[1,1].r:= (1);
  ue[2,1].r:= (1);
  ue[3,1].r:= (1);
  ue[3,2].r:= (-1);
  ue[2,3].r:= (-1);

  initmatrix(uf);
  uf[1,1].r:= (1) ;
  uf[2,2].r:= (1) ;
  uf[3,3].r:= (-1);
  uf[3,2].r:= 1;

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

  mkscalar(tmp,minus1);
  matmult(ua,tmp,tua);
  matmult(ub,tmp,tub);
  matmult(uc,tmp,tuc);
  matmult(ud,tmp,tud);
  matmult(ue,tmp,tue);
  matmult(uf,tmp,tuf);

  writeln('Initialized order 2 matrices.');
  (* check gcd routines *)

  writeln;
  scal2quad(x1,7 );
  scal2quad(x2,84);
  writeln('The norm of alpha = ',normquad(alpha));
  writeln('The norm of alphabar = ',normquad(alphabar));
  dividequad(x2,x1,x3);
  writeln('84 divided by 7 = ');
  printquad(x3);


  gcdquad(x1,x2,x3);
  writeln('The GCD of 7 and 84 is ');
  printquad(x3);writeln;

  x1:= alpha;
  x2:= alphabar;
  gcdquad(x1,x2,x3);
  writeln('The GCD of alpha and alphabar is one, but is calculated as ');
  printquad(x3);
  writeln;

  e1:= tua;
  e2:= tuf;
(*  printmatrix(tua);
  convertmatrix(tua,rtua);
  printratmatrix(rtua);
  ratmatmult(rtua,rtua,rtmp1);
  printratmatrix(rtmp1);
  readln;
  writeln;
  matmult(tua,tua,tmp1);
  printmatrix(tmp1);   *)
  writeln;
  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 compare_quad(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 compare_mult(tmp4,e2) then writeln('Conjugation verified') else
writeln('Conjugation verification failed.');
readln;
tau:=tmp2;
*)
mkscalar(s,zero);
s[1,3]:= unitquad;
s[2,1]:= unitquad;
s[2,3]:= alphabar; multquad(s[2,3],minus1,s[2,3]);
s[3,2]:= unitquad;
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]:= unitquad;
  v1[2]:= alpha; v1[2].r:= 2;
  v1[3]:= alpha ;
  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,zero);
(* 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]:=unitquad;
t[2,1]:= alpha; t[2,1].r:=2;
t[2,2]:= alpha;
t[2,3]:= minus1;
t[3,1]:= alpha;
t[3,2] := minus1;
t[3,3]:= alphabar;

gg:=gorder(t);
writeln('Order of t is ',gg);
writeln('For T = ');
printmatrix(T);
matpower(t,tt,2);
writeln('The square of T = TT = ');
printmatrix(tt);

writeln( 'TSTT = ');
matmult(t,s,tmp1);
matmult(tmp1,tt, tmp2);
writeln;
printmatrix(tmp2);
readln;
if compare_quad(tmp2,ss) then writeln('S conjugated by T is S^2') else
writeln('Wrong conjugate of S by T.');
readln;

v1[1].r := 3;
v1[2].r := 2;
v1[2].a := 2;
v1[3]:= alpha;

v2[1]:= zero;
v2[2]:= unitquad;
v2[3]:= zero;

writeln('The vector v1 should be invariant under T , v1 = ');
printvector(v1);
writeln('Here is Tv1 = ');
vectmat(T,v1,vv);
printvector(vv);
writeln('Then have to find one other vector to generate a rank two cyclic ');
writeln('Take the generator to be v2  = ');
vectmat(T,v2,v3);
printvector(v2);
writeln(' And the other to be T v2 = v3 = ');
printvector(v3);

9998:
readln;

A[1,1]:= v1[1];
A[2,1]:= v1[2];
A[3,1]:= v1[3];

a[1,2]:= v2[1];
A[2,2]:= v2[2];
a[3,2]:= v2[3];

a[1,3]:= v3[1];
a[2,3]:= v3[2];
a[3,3]:= v3[3];

writeln('The transformation matrix is A = ');
printmatrix(A);
writeln;
convertmatrix(A,rA);
if ratinvert(rA,rA_inv)  then begin
   writeln;writeln('It''s inverse is ');
   printratmatrix(rA_inv);
   writeln;
end else writeln('rA did not invert.!');
readln;
writeln('Check on this inverse.   ');
writeln;
ratmatmult(rA_inv,rA,rtmp2);
printratmatrix(rtmp2);
readln;
writeln('A is supposed to put T in companion form. However ');
writeln('A_inv * T * A = ');
writeln; convertmatrix(t,rt);
ratmatconjr(rA,rt,rtnew);

writeln('Here is A conj. T = ');
printratmatrix(rtnew);
ratmatpower(rtnew,rC,3);
writeln;writeln('Here is the conjugate to third power.');
writeln;
printratmatrix(rC);
writeln;
readln;

writeln;
writeln('Here is the A_inv S A = ');
writeln;
convertmatrix(s,rs);
ratmatconjr(rA,rS,rsnew);
printratmatrix(rsnew);
writeln;
readln;

9997:
{ 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 }
writeln;
writeln('Conjugating TAU by A =');

 mkscalar(taunew,zero);
 taunew[1,1]:=minus1;
 taunew[3,2]:=unitquad;
 taunew[2,3]:=unitquad;
 writeln;
 convertmatrix(taunew,rtaunew);
 printratmatrix(rtaunew);
 writeln;
 readln;
(* now check the relations *)

ratmatpower(rsnew,rC,7);
writeln;
writeln('SNEW ^7 =  ID ? ');
printratmatrix(rC);
writeln;
readln;

ratmatconj(rtnew,rSNEW,rC);
ratmatpower(rSnew,rB,2);
writeln;
writeln('Is T S T_inv = SS?');
writeln;
printratmatrix(rB);
writeln;
printratmatrix(rC);
readln;

writeln('TAU ^2 = ? ');
ratmatmult(rtaunew,rtaunew,rB);
printratmatrix(rB);
writeln;readln;

writeln('Is T^3 = ? ');
ratmatmult(rt,rt,rB);
ratmatmult(rB,rt,rC);
printratmatrix(rC);
writeln;
readln;

writeln;
writeln('Is tau T tau = T_inv? ');
ratmatmult(rtaunew,rtnew,rB);
ratmatmult(rB,rtaunew,rC);
printratmatrix(rC);
writeln;
ratmatpower(rtnew,rC,2);

printratmatrix(rC);
writeln;
readln;

writeln;
writeln;
ratmatmult(rtaunew,rsnew,rB);
for ii:= 1 to 9 do begin
  ratmatpower(rB,rc,ii);
  writeln;
  writeln('tau S to ',ii:2,'-th power is ');
  printratmatrix(rC);
  writeln;
  readln;
end;