
(*-----------------------Examples - SL_2(Z) -------------------------*)
<<cox.m

SL0:=DiSL[1]
SL1:=DiSL[x-3]
SL2:=DiSL[2+x]

LP0:=DiLP[1]
LP1:=DiLP[x-3]
LP2:=DiLP[2+x]

(* ------------------SL(2,3,inf) Graphics routines -------------------*)
Color = Red

(* Graphics for geodesic in Poincare' model from Minkowski vector of positive norm *)
PoincareCircle[{t_,x_,y_}]:=If[t != 0, Circle[{x/t,y/t},Sqrt[x^2+y^2-t^2]/Abs[t]],
					Line[{{-y,x}/Sqrt[x^2+y^2],{y,-x}/Sqrt[x^2+y^2]}]]

GraphRoots[rs_]:=Graphics[Join[{Circle[{0,0},1]},Map[PoincareCircle[QtoM[#]]&,rs]]]

(* Convert from QSL[2,1] to Minkowski space *)
QtoM[v_]:={{0, 1/2, 0}, {1, -1/2, 0}, {0, -1, 1}}.v

(* Graph for SL_2(R) with Coxeter element *)
GraphSL[rs_]:=Block[{geod,g1,g2,g,q},
		q = QSL[2,1];
		geod = Eigenvectors[Cox[q]];
		geod = Select[geod,(And @@ Map[IntegerQ,#])&];
		g1=GraphRoots[rs];
		g2=GraphRoots[geod];
		g = Graphics[Join[{Color},g1[[1]],{Black,AbsoluteThickness[2]},g2[[1]]]]; 
		Return[g]
]

(*-----------------------SL[2,R] specific routines -------------------------*)
DiSL[poly_]:= Block[{f,g,q,rs,u},
		q = QSL[2,1];
		f = Cox[q];
		u = PofM[poly,f];
		rs = GenRoots[q,20];
		rs = Select[rs,RangeQ[#,u]&];
		g = Show[GraphSL[rs],PlotRange->{{-1.1,1.1},{-1.1,1.1}}];
		Export["f2.ps",g,"eps"];
]

(*-----------------------Demonstrate linear programming test ---------------*)
(* (a+2)(a-3) = 0 *)
Astar[v_]:={v.q.p1.v,v.q.p2.v}

DiLP[poly_]:=Block[{a,basis,e1,e2,e3,f,g,i,m,m1,m2,ns,p1,p2,q,r,u,vs},
		q = QSL[2,1];
		f = Cox[q];
		a = f + Inverse[f];
		u = PofM[poly,f];
		p1 = (a+2 Id[a])/5;
		p2 = (3Id[a]-a)/5;
		(* Linear programming minimum *)
		e1={1,0,0}; e2={0,1,0}; e3={0,0,1};
		basis = {e1,e2,e3,e1+e2,e1+e3,e2+e3};
		basis = Map[u.#&,basis];
		{m1,m2} = LatticeReduce[Map[Astar,basis]];
		m = Flatten[Table[i*m1 + j*m2,{i,-10,10},{j,-10,10}],1];
		r = Min[Select[Map[Norm,m],#>0&]];
		g = ListPlot[m,PlotRange->{{-.6r,3.1r},{-.3r,3.1r}},AspectRatio->Automatic,PlotMarkers->{Automatic,Large},BaseStyle->{FontSize->24}];
		Export["f2.ps",g,"eps"];
		Print["Lattice defined by ",poly," has index: ",Abs[Det[u]]];
		Print["A* Lattice vectors: ",{m1,m2}];
		Print["Minimum by LP:      ",SalemOb[QSub[q,poly],f]];

	(* Search cube for minimal solution *)
	vs = IntegerBox[{{0,5},{-5,5},{-5,5}}];
	vs = Map[u.#&,vs];
	vs = Select[vs,Astar[#][[1]] > 0 && Astar[#][[2]] >= 0&];
	ns = Map[#.q.#&,vs];
	i = Position[ns,Min[ns]][[1,1]];
	v = vs[[i]];
		Print["Lattice min: ",v.q.v," at v = ",v," = u( ",Inverse[u].v,") with psi = ",Astar[v]];

		(* Obstructing roots *)
		o1 = Inverse[u].e1;
		o2 = Inverse[u].e2;
		o3 = Inverse[u].e3;
		Print["If integral, obstructing root remains:  ",{o1,o2,o3}];
]
ForPaper:=(Color=Black; SL0; Run["mv -f f2.ps ../pics/sl0.ps"];
		        SL1; Run["mv -f f2.ps ../pics/sl1.ps"];
		        SL2; Run["mv -f f2.ps ../pics/sl2.ps"];
		        LP2; Run["mv -f f2.ps ../pics/lp2.ps"];)

q = QSL[2,1];
f = Cox[q];

