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

(* This file claim to show that the image is Zariski dense with finitely many exceptions in V 

MEMO: MatrixLog doesn't work well with a variable in it.*)

(* 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];

el=Factor[ y.X.Y.x.x.Y.X.y];
Factor[ el.x-x.el]
Factor[ Det[ el-Q*id4]]

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

(* Here's the vector which exponentiates to x *)

myvec= {{(-1+V)/(4 V),(-1+V)/(4 V),(1+V)/(2 V),-((-1+V)/(2 Sqrt[3] V))},{1/4 (-1+1/V),1/4 (-1+1/V),-((1+V)/(2 V)),(-1+V)/(2 Sqrt[3] V)},{1,1,0,0},{0,0,0,0}};

Factor[ myvec.myvec]
Factor[ myvec.myvec.myvec]

Factor[ id4 + myvec + 1/2myvec.myvec - x]

(* Now compute the representation of the figure eight knot group on the Lie algebra of SL(4,R)*)

basis={    
{{0,1,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}},  
{{0,0,1,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}}, 
{{0,0,0,1},{0,0,0,0},{0,0,0,0},{0,0,0,0}},
 
{{0,0,0,0},{1,0,0,0},{0,0,0,0},{0,0,0,0}},
{{0,0,0,0},{0,0,1,0},{0,0,0,0},{0,0,0,0}}, 
{{0,0,0,0},{0,0,0,1},{0,0,0,0},{0,0,0,0}},

 {{0,0,0,0},{0,0,0,0},{1,0,0,0},{0,0,0,0}}, 
{{0,0,0,0},{0,0,0,0},{0,1,0,0},{0,0,0,0}},
 {{0,0,0,0},{0,0,0,0},{0,0,0,1},{0,0,0,0}}, 

{{0,0,0,0},{0,0,0,0},{0,0,0,0},{1,0,0,0}}, 
{{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,1,0,0}}, 
{{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,1,0}},

 {{1,0,0,0},{0,-1,0,0},{0,0,0,0},{0,0,0,0}},
 {{0,0,0,0},{0,1,0,0},{0,0,-1,0},{0,0,0,0}},
 {{0,0,0,0},{0,0,0,0},{0,0,1,0},{0,0,0,-1}}};
Length[basis]

coeff=Array[ cc, {15}];

aa={};
Do[ AppendTo[ aa , coeff/.Solve[  Together[x.basis[[pp]].X] == Sum[ coeff[[pp]]*basis[[pp]] , {pp,1,15} ], coeff][[1]]    ] , {pp,1,15} ]

MatrixForm[Transpose[aa]]

Det[aa]

bb={};
Do[ AppendTo[ bb , Simplify[coeff/.Solve[  Together[y.basis[[pp]].Y] == Sum[ coeff[[pp]]*basis[[pp]] , {pp,1,15} ], coeff][[1]]]    ] , {pp,1,15} ]

MatrixForm[Transpose[bb]]

coeff=Array[cc,{15}]

AA=Inverse[aa];
BB=Inverse[bb];
uno=Together[ (bb.AA.BB.aa)];
Factor[ (uno).bb- aa.(uno)]

(* Compute myvec in this basis*)

myvec1=coeff/. Solve[ myvec== Sum[ coeff[[pp]]*basis[[pp]] , {pp,1,15} ] , coeff][[1]]

Clear[V]
V=N[12,200];

MatrixForm[  N[ Factor[ RowReduce[{myvec1, aa.myvec1,  aa.aa.myvec1, aa.aa.aa.myvec1,bb.myvec1,bb.bb.myvec1,aa.bb.myvec1,bb.aa.myvec1,aa.bb.aa.myvec1, bb.aa.bb.myvec1,aa.bb.aa.aa.myvec1,bb.bb.aa.myvec1, bb.aa.bb.aa.myvec1, aa.aa.aa.bb.myvec1, bb.bb.bb.aa.myvec1}]]   ,10]]

(* Show that with finitely many exceptions, there's 15 independent elements in the orbit of myvec1 under figure eight*)

Clear[V];

Factor[ Det[{myvec1, aa.myvec1,  aa.aa.myvec1, aa.aa.aa.myvec1,bb.myvec1,bb.bb.myvec1,aa.bb.myvec1,bb.aa.myvec1,aa.bb.aa.myvec1, bb.aa.bb.myvec1,aa.bb.aa.aa.myvec1,bb.bb.aa.myvec1, bb.aa.bb.aa.myvec1, aa.aa.aa.bb.myvec1, bb.bb.bb.aa.myvec1}]]