(*
	(c) 2011 Alex Kontorovich

	This file complements the expository note on an arithmetic surface
*)

(*
	We work throughout with the following indefinite, anisotropic quadratic form (represented by a matrix)
*)
Q={{1,0,0},{0,1,0},{0,0,-3}};
Q//MatrixForm

(* 
	m_x are the symmetric matrices of Det = - Q(x) 
*)
mx={{z Sqrt[3]-y,x},{x,z Sqrt[3]+y}} ;
mx//MatrixForm

Det[mx]

(* 
	One can read off the elements x,y,z from m_x, e.g. : 
*)
Tr[mx]/(2Sqrt[3])

(*
	Let g be an element of SL(2,R)
*)
g={{a,b},{c,d}};
g//MatrixForm

(* 
	Define the action of g on mx 
*)
mxp=g.mx.Transpose[g]//Expand; 
(* I use p for the prime symbol, p = ' *)
mxp//MatrixForm

(*
	Output to TeX
*)
mxp//TeXForm

(* Pick off x' element *)
xp=mxp[[2,1]]//Expand

(* Pick off z' element *)
zp=Tr[mxp]/(2Sqrt[3])//Expand

(* Pick off y' element *)
yp=mxp[[2,2]]-zp Sqrt[3]//Expand

(*
	Action is linear in x,y,z, so can be represented as a matrix; this gives iota
*)
iotaG={
{D[xp,x],D[yp,x],D[zp,x]},
{D[xp,y],D[yp,y],D[zp,y]},
{D[xp,z],D[yp,z],D[zp,z]}
}//Simplify;
iotaG//MatrixForm

(* 
	Check that iotaG is in SO_Q
*)
iotaG.Q.Transpose[iotaG]/(a d-b c)^2//Simplify//MatrixForm

(* 
	Defining quaternions as matrices
*)
e={{1,0},{0,1}};
i=3^(1/2){{1,0},{0,-1}};
j=3^(1/2){{0,-1},{-1,0}};
k=1/3i.j;

i.i//MatrixForm

j.j//MatrixForm

k.k//MatrixForm

(* 
	The general form of a quaternion
*)
u=(a e + b i + c j+ d k);
u//MatrixForm

Det[u]

(*
	Where does iota send the quaternion u?
*)
iotaU=iotaG/.a->aa+Sqrt[3]bb/.d->aa-Sqrt[3]bb/.b->-dd-Sqrt[3]cc/.c->dd-Sqrt[3]cc/.aa->a/.bb->b/.cc->c/.dd->d//Simplify;
iotaU//MatrixForm

(* 
	Now try to find some generators of Gamma
*)
bound=9;
Solns=FindInstance[(a^2)-3 (b^2)-3 (c^2)+(d^2)==1
&&
Abs[a]<bound
&&
Abs[b]<bound
&&
Abs[c]<bound
&&
Abs[d]<bound
,{a,b,c,d},Integers,400];
Length[Solns]

(* 
	This turns solutions into matrices 
*)
listToMat[list_]:=
Simplify[u/.list];

(*
	Now we have 2x2 matrices in SL(2,R)
*)
Gens=listToMat/@Solns;
(* we want a symmetric set of generators *)
GensI=Inverse/@Gens;
Gammas=Union[Gens,GensI];
Length[Gammas]

(*
	We now construct a Dirichlet domain corresponding to the action of the above matrices on some base point
*)

(* Pick a point *)
basept=2I;
(* Action of SL(2,R) on H *)
zTogz[z_,g_]:=Simplify[Re[(g[[1,1]]z+g[[1,2]])/(g[[2,1]]z+g[[2,2]])]+I  Im[(g[[1,1]]z+g[[1,2]])/(g[[2,1]]z+g[[2,2]])]];
(* Going between C and R^2 *)
CtoVect[z_]:={Re[z],Im[z]};
VectToC[p_]:=p[[1]]+I p[[2]];

(* Action of matrix on basept *)
MatTimesBase[mat_]:=zTogz[basept,mat];
(* This is the orbit of Gammas found above on basept *)
OrbitAll=MatTimesBase/@Gammas;

(* Don't waste time with elements too far down *)
ImInf=.11;
ImSup=10;
ReInf=-4;
ReSup=4;
Orbit=Select[Union[OrbitAll],ReInf<Re[#]<ReSup && ImInf<Im[#]<ImSup &];
Length[Orbit]

(* Lets see the orbit *)
Cpts=Simplify[CtoVect/@Orbit];
ptsPic=ListPlot[Cpts,PlotRange->All]

(* Given this orbit, we'll try to draw a Dirichlet domain *)

(* Hyperbolic distance functions *)
dist[u_,v_]:=Simplify[Norm[u-v]^2/(4 Im[u]Im[v])];
(* Compute distance of basept to others *)
setofdistances[z_]:=Module[{},
Table[dist[Orbit[[j]],z],{j,1,Length[Orbit]}]
];
(* Only minimal distance is in the Dirichlet domain *)
IsMin[z_]:=(dist[basept,z]<=Min[setofdistances[z]]);

(* Plot Dirichlet domain *)
regPic=RegionPlot[IsMin[VectToC[{x,y}]],{x,-4,4},{y,0,4},AspectRatio->Automatic]
Beep[]

pic=Show[regPic,ptsPic,PlotRange->{{-4,4},{0,4}}]

(*
	Formula defining a geodesic of all points equidistant to two given points, x0 + I y0 and x1 + I y1
*)
(Reduce[Assuming[
x∈Reals &&y∈Reals&&
x0∈Reals &&y0>0&&
x1∈Reals &&y1>0
,dist[x+I y,x0+I y0]==dist[x+I y,x1+I y1]]&&y>0,{x,y},Reals]//FullSimplify)[[2]][[1]]

(*
	So the "y" value of a geodesic between two points z0 and z1 in the upper half plane is given by:
*)
geodEqui[x_,z0_,z1_]:=(Sqrt[((-(((x-x1))^2) y0+((((x-x0))^2)+(y0^2)) y1-y0 (y1^2))/(y0-y1))])/.x0->Re[z0]/.y0->Im[z0]/.x1->Re[z1]/.y1->Im[z1];

(*
	Make a function that plots the above geodesic for z0 and z1
*)
plotGeo[z0_,z1_]:=Plot[geodEqui[x,z0,z1],{x,-10,10}];
(*
	And a function with z0=basept
*)
plotGeoBase[z1_]:=plotGeo[basept,z1];

(*
	We want the orbit with no base point
*)
OrbitNoBase=Drop[Orbit,Position[Orbit,basept][[1]]];

(*
	Plot all geodesic pictures
*)
geoPics=Show[plotGeoBase/@OrbitNoBase, PlotRange->{{-4,4},{0,4}},AspectRatio->Automatic]

(*
	Show them atop the Dirichlet domain and orbit
*)
Show[pic, geoPics]

(*
	Which Matrices did we end up using?
*)
GammasUsed={};
For[j=1,j≤Length[Orbit],j++,
GammasUsed=Append[GammasUsed,Gammas[[
(Flatten[Position[OrbitAll,Orbit[[j]]]][[1]])
]]
];
Print[TeXForm[GammasUsed[[j]]
]];
];

(*
	From a matrix back to (a,b,c,d)
*)
matToQuat[mat_]:=FullSimplify[{
(mat[[1,1]]+mat[[2,2]])/(2),
(mat[[1,1]]-mat[[2,2]])/(2Sqrt[3]),
(mat[[2,1]]+mat[[1,2]])/(-2Sqrt[3]),
(mat[[2,1]]-mat[[1,2]])/(2)
}];

(*
	And the corresponding elements of the quaternion algebra
*)
matToQuat/@GammasUsed//TeXForm