(*
	(c) 2013 Alex Kontorovich 
*)

(*
	Finds "low-lying" (with A=2) closed geodesics over a number field
*)

(* Take Q(sqrt[discr]) *)
discr=5;

(* Find all solutions to the Pell equation for traces in Q(sqrt[discr]) *)
expr=Reduce[t^2-discr s^2==4&&Element[t|s,Integers]][[4]][[4]][[2]]/.C[1]->c

(* make a list of the first few desired trace values *)
trWants=Table[Simplify[expr],{c,1,6}]

(* The main code *)
Timing[
(* Largest trace to consider *)
TrMax=Max[trWants];

(* generating matrices and their inverses *)
a1={{0,1},{1,1}};
a1i=Inverse[a1];
a2={{0,1},{1,2}};
a2i=Inverse[a2];

(* only want even letters *)
a11=a1.a1;
a12=a1.a2;
a21=a2.a1;
a22=a2.a2;
a11i=a1i.a1i;
a12i=a2i.a1i;
a21i=a1i.a2i;
a22i=a2i.a2i;

(* Iterating function *)
getmore[]:=Module[
{},
(* add a letter to the end *)
m=m.a11;
If[Tr[m]≤TrMax, (* If not too big, enter loop *)
If[MemberQ[trWants,Tr[m]] (* If the trace is one of the desired ones, store the matrix *)
,
AppendTo[mats,m];
];
getmore[]; (* Recursion step *)
];
(* Undo adding the letter *)
m=m.a11i;

(* Now do the same thing for the other letters *)
m=m.a12;
If[Tr[m]≤TrMax,
If[MemberQ[trWants,Tr[m]]
,
AppendTo[mats,m];
];
getmore[];
];
m=m.a12i;

m=m.a21;
If[Tr[m]≤TrMax,
If[MemberQ[trWants,Tr[m]]
,
AppendTo[mats,m];
];
getmore[];
];
m=m.a21i;

m=m.a22;
If[Tr[m]≤TrMax,
If[MemberQ[trWants,Tr[m]]
,
AppendTo[mats,m];
];
getmore[];
];
m=m.a22i;

];

(* Initializing *)
mats={};
m={{1,0},{0,1}};
(* And run the code *)
getmore[];
]
Beep[] (* This may take some time *)
Length[mats]

(* Compute set of fixed points, which are the corresponding closed geodesics *)
fixedPt[m_]:=((m[[1,1]]-m[[2,2]]+(Sqrt[Tr[m]^2-4Det[m]]))/(2m[[2,1]]));
pts=Union[Simplify[fixedPt/@mats]]; (* Take only primitive geodesics *)
Length[pts]

ContinuedFraction/@pts//Sort (* Here are their continued fractions *)

FromContinuedFraction/@%

(* Do it all again, now with Q(sqrt[3])*)

(* Take Q(sqrt[discr]) *)
discr=12;

(* Find all solutions to the Pell equation for traces in Q(sqrt[discr]) *)
expr=Reduce[t^2-discr s^2==4&&Element[t|s,Integers]][[4]][[4]][[2]]/.C[1]->c

(* make a list of the first few desired trace values *)
trWants=Table[Simplify[expr],{c,1,6}]

(* The main code *)
Timing[
(* Largest trace to consider *)
TrMax=Max[trWants];

(* generating matrices and their inverses *)
a1={{0,1},{1,1}};
a1i=Inverse[a1];
a2={{0,1},{1,2}};
a2i=Inverse[a2];

(* only want even letters *)
a11=a1.a1;
a12=a1.a2;
a21=a2.a1;
a22=a2.a2;
a11i=a1i.a1i;
a12i=a2i.a1i;
a21i=a1i.a2i;
a22i=a2i.a2i;

(* Iterating function *)
getmore[]:=Module[
{},
(* add a letter to the end *)
m=m.a11;
If[Tr[m]≤TrMax, (* If not too big, enter loop *)
If[MemberQ[trWants,Tr[m]] (* If the trace is one of the desired ones, store the matrix *)
,
AppendTo[mats,m];
];
getmore[]; (* Recursion step *)
];
(* Undo adding the letter *)
m=m.a11i;

(* Now do the same thing for the other letters *)
m=m.a12;
If[Tr[m]≤TrMax,
If[MemberQ[trWants,Tr[m]]
,
AppendTo[mats,m];
];
getmore[];
];
m=m.a12i;

m=m.a21;
If[Tr[m]≤TrMax,
If[MemberQ[trWants,Tr[m]]
,
AppendTo[mats,m];
];
getmore[];
];
m=m.a21i;

m=m.a22;
If[Tr[m]≤TrMax,
If[MemberQ[trWants,Tr[m]]
,
AppendTo[mats,m];
];
getmore[];
];
m=m.a22i;

];

(* Initializing *)
mats={};
m={{1,0},{0,1}};
(* And run the code *)
getmore[];
]
Beep[] (* This may take some time *)
Length[mats]

(* Compute set of fixed points, which are the corresponding closed geodesics *)
fixedPt[m_]:=((m[[1,1]]-m[[2,2]]+(Sqrt[Tr[m]^2-4Det[m]]))/(2m[[2,1]]));
pts=Union[Simplify[fixedPt/@mats]]; (* Take only primitive geodesics *)
Length[pts]

ContinuedFraction/@pts//Sort (* Here are their continued fractions *)

FromContinuedFraction/@%