(* Ballas-Thistlethwaite deformation of the figure eight knot group  *)

(* This file claims to show that the traces of elements in this representation are integers if the element V is a unit.

The idea is this (implemented below):

By Burnside, we can find elements g_1,....., g_15 which are a basis for 4 x 4 matrices. Let g_1^*,....., g_15^* be the dual basis
wrt trace, i.e.  tr(g_i.g_j^* ) = \delta_{ij}



Then one can compute the matrices which express how an expression  

g = \sum_i  \alpha_i g_i^*              (*)

is transformed when multiplied by x, y etc. These are the matrices locmatx, locmaty which are below.


The key idea is the following. I've worked out locmatx, etc for TWO different basis choices g_1,....., g_15



* For the first one, one sees locmatx, etc have denominators involving powers of V, (-1+V) and  (1+V)

This means that at worst the denominators of the coefficients \alpha_i in the expressions of type (*) are 
powers of those polynomials.




* For the second one, one sees locmatx, etc have denominators involving powers of V, (1+3 V+(V^2)) and (1+3 V+4 (V^2))

Of course trace of the matrix g doesn't depend on choice of basis: We conclude that the only possible denominators
for traces are powers of V. QED



 *)

(* Details now follow *)

(* Here's the initial set up: *)

Clear[V]
Clear[Q]

y={{ V+1/2, -V+1/2,1/2, (7-4V)/Sqrt[12]},
{1/2,1/2,1/2,(4/V-1)/Sqrt[12]},
{V/2,-V/2, 1, (1-V)/Sqrt[3]},
{Sqrt[3]*V/2,-Sqrt[3]*V/2, 0, 2-V}};

Y=Inverse[y];

x={{3/2,1/2,1/2(V+1)/V, (1/V-1)/Sqrt[12]},
{-1/2,1/2, -1/2(V+1)/V, (1-1/V)/Sqrt[12]},
{1,1,1,0},{0,0,0,1}};
X=Inverse[x];

id4=IdentityMatrix[4];

Factor[ (y.X.Y.x).y- x.(y.X.Y.x)]

tr[Y__]:= Factor[ Sum[  Y[[k,k]] , {k,1,4} ]]

(* FIRST BASIS Computation begins here *)

(* This code checks the 16 elements below are linearly independent *)

id=id4=IdentityMatrix[4];
id16=IdentityMatrix[16];

coeff=Array[ cc, {16} ];

basis={id,x,X,y, Y, x.y, Y.X, y.x, X.Y, x.x.y, Y.X.X, X.X.Y, y.x.x,    x.y.X.y, Y.x.Y.X, x.y.X.Y};

coeff/.Solve[ Sum[ coeff[[kk]]*basis[[kk]] , {kk,1,Length[basis]} ] == (id-id) , coeff ][[1]]

(*   This code generates the dual basis to that basis *)

onevec= Factor[ Sum[ coeff[[kk]]*basis[[kk]] , {kk,1,Length[basis]} ]];
dualbasis={};

Do[  
loccoeff=coeff/.Factor[ Solve[ { tr[ basis[[1]].onevec],tr[ basis[[2]].onevec],tr[ basis[[3]].onevec],tr[ basis[[4]].onevec],tr[ basis[[5]].onevec],tr[ basis[[6]].onevec],tr[ basis[[7]].onevec],tr[ basis[[8]].onevec],tr[ basis[[9]].onevec],tr[ basis[[10]].onevec],tr[ basis[[11]].onevec],tr[ basis[[12]].onevec],tr[ basis[[13]].onevec],tr[ basis[[14]].onevec],tr[ basis[[15]].onevec],tr[ basis[[16]].onevec]}==id16[[AA]] , coeff]][[1]] ;

temp= Factor[ Sum[ loccoeff[[kk]]*basis[[kk]] , {kk,1,Length[basis]} ]    ];

AppendTo[ dualbasis, temp];
, {AA,1,16} ]

(*    Check that it's a dual basis : *)

mymat=Array[ qq, {16,16} ];
Do[   mymat[[aa,bb]] =  tr[   basis[[aa]].dualbasis[[bb]]  ] , {aa,1,16},{bb,1,16} ];
MatrixForm[ mymat]

(* Now compute the change of coefficient matrices locmatx etc

Notice denominators are all V, (V+1), (V-1) and that the matrix has determinant =1 in both cases.

 *)

locmaty=Array[ aa, {16,16} ];
Do[ locmaty[[p,q]]=  Factor[tr[ basis[[p]].y.dualbasis[[q]]] ]; , {p,1,16},{q,1,16} ];
MatrixForm[locmaty]

Factor[Det[locmaty]]

locmatx=Array[ aa, {16,16} ];
Do[ locmatx[[p,q]]=  Factor[tr[ basis[[p]].x.dualbasis[[q]]] ]; , {p,1,16},{q,1,16} ];
MatrixForm[locmatx]

Factor[Det[locmatx]]

(*    Here's a sample computation of how these work. The first two cells below, exhibit the coefficients in the expression
for id4 (= 4x 4 identity matrix) and x.x.Y 

It's then followed by locmatx.Inverse[locmaty].vecid and you'll see they are the same as they should be. *)

vecid=Factor[ coeff/.Solve[id4 ==  Factor[ Sum[ coeff[[kk]]*dualbasis[[kk]] , {kk,1,Length[basis]} ]], coeff][[1]]]

Factor[ coeff/.Solve[x.x.Y ==  Factor[ Sum[ coeff[[kk]]*dualbasis[[kk]] , {kk,1,Length[basis]} ]], coeff][[1]]]

Factor[ locmatx.locmatx.Inverse[locmaty].vecid]

(* FIRST BASIS Computation ends here *)

(* SECOND BASIS Computation starts here *)

(* Do the same, but with different group elements. Notice the denominators are now V, (1+3 V+(V^2)) and (1+3 V+4 (V^2))*)

id=IdentityMatrix[4];
id16=IdentityMatrix[16];

coeff=Array[ cc, {16} ];

basis={id,x,x.x, y, y.y, x.y,y.x, y.x.x, y.x.y,x.y.x, x.y.y, x.y.x.y,y.x.y.x,x.x.y.x, x.x.y.y, y.x.y.y};

coeff/.Solve[ Sum[ coeff[[kk]]*basis[[kk]] , {kk,1,Length[basis]} ] == (id-id) , coeff ][[1]]

onevec= Factor[ Sum[ coeff[[kk]]*basis[[kk]] , {kk,1,Length[basis]} ]];
dualbasis={};

Do[  

loccoeff=coeff/.Factor[ Solve[ { tr[ basis[[1]].onevec],tr[ basis[[2]].onevec],tr[ basis[[3]].onevec],tr[ basis[[4]].onevec],tr[ basis[[5]].onevec],tr[ basis[[6]].onevec],tr[ basis[[7]].onevec],tr[ basis[[8]].onevec],tr[ basis[[9]].onevec],tr[ basis[[10]].onevec],tr[ basis[[11]].onevec],tr[ basis[[12]].onevec],tr[ basis[[13]].onevec],tr[ basis[[14]].onevec],tr[ basis[[15]].onevec],tr[ basis[[16]].onevec]}==id16[[AA]] , coeff]][[1]] ;

temp= Factor[ Sum[ loccoeff[[kk]]*basis[[kk]] , {kk,1,Length[basis]} ]    ];

AppendTo[ dualbasis, temp];
, {AA,1,16} ]

mymat=Array[ qq, {16,16} ];

Do[   mymat[[aa,bb]] =  tr[   basis[[aa]].dualbasis[[bb]]  ] , {aa,1,16},{bb,1,16} ];
MatrixForm[ mymat]

locmatx=Array[ aa, {16,16} ];
Do[ locmatx[[p,q]]=  Factor[tr[ basis[[p]].x.dualbasis[[q]]] ]; , {p,1,16},{q,1,16} ]
MatrixForm[locmatx]

Factor[ Det[locmatx]]

locmaty=Array[ aa, {16,16} ];
Do[ locmaty[[p,q]]=  Factor[tr[ basis[[p]].y.dualbasis[[q]]] ]; , {p,1,16},{q,1,16} ];
MatrixForm[locmaty]

Factor[ Det[locmaty] ]

(* Irreducibility *)

(* Case A: There's no 1 or 3 dimensional invariant subspace*)

com2poly=Factor[ Det[ x.x.y.Inverse[y.x.x]-Q*id4]]

Collect[ Expand[com2poly],Q]

Factor[com2poly/.Q->1]

(* Case B: There's no 2 dimensional invariant subspace*)

el=Factor[ y.X.Y.x.x.Y.X.y]

Factor[ x.el-el.x]

Factor[ Det[ el - Q*id4]]

vec=Simplify[Sqrt[3] (-1+(V^2))NullSpace[ el-1/V^3id4][[1]]]

mymat={vec, y.vec, x.y.vec, y.y.y.vec};
Factor[ Det[mymat]]

Solve[(-1-3 V+(V^2))==0, V]

Simplify[ RowReduce[ Simplify[ mymat/.{V->(1/2) (3+(Sqrt[13]))}]]]

xt=Transpose[X];
Xt=Transpose[x];
yt=Transpose[Y];
Yt=Transpose[y];

elt=Factor[ yt.Xt.Yt.xt.xt.Yt.Xt.yt]

vect=Simplify[(2-2 V)NullSpace[ elt-V^3id4][[1]]]

mymatt={vect, yt.xt.vect, xt.yt.vect, yt.xt.yt.yt.vec};
Factor[ Det[mymatt]]

N[Solve[ (-2+5 V+15 (V^2)-15 (V^3)-6 (V^4)+9 (V^5))==0, V], 20]

mymatt={vect, yt.xt.vect, xt.yt.vect, yt.yt.yt.vec};
Factor[ Det[mymatt]]

N[Solve[ (-8+25 V+12 (V^2)-40 (V^3)+14 (V^4)+6 (V^5))==0, V], 20]