(*
	(c) 2007 Alex Kontorovich
*)

(* First show all Pythagorean triples, sorted according to number of prime factors of the product xyz *)

(* Parameters *)

ptSize=Large;
range=100;


(* The Code *)

(* make a cone *)
coneBot=Plot3D[-(Sqrt[(x^2)+(y^2)]),{x,-(range^2),(range^2)},{y,-(range^2),(range^2)},PlotPoints->2 range,Mesh->False,Axes->False,Boxed->False,FaceGrids->None,AspectRatio->1,PlotRange->{{-(range^2),(range^2)},{-(range^2),(range^2)},{-(range^2),0}},ViewPoint->{3,-1,1},SphericalRegion->False,ClippingStyle->None];

gcd[v_]:=GCD[v[[1]],v[[2]],v[[3]]];
primTriples={};

(* make all triples *)
For[i=- range,i≤ range,i=i+2,
For[j=- range-1,j≤  range,j=j+2,
If[(Abs[i]≠Abs[j])&&(i j ≠0)&& (GCD[i,j]==1),
AppendTo[primTriples,{i^2-j^2,2i j,-i^2-j^2}];
AppendTo[primTriples,{2i j,i^2-j^2,-i^2-j^2}];
];
];
];

(* Sort triples according to number of prime factors of xyz/60 *)
CTriplesDown2 = Select[primTriples, Total[Last/@ FactorInteger[ Abs[#[[1]]  #[[2]] #1[[3]]/60]]]≤ 4&];
CTriplesDown3 = Select[primTriples, Total[Last/@ FactorInteger[ Abs[#[[1]]  #[[2]]#1[[3]] /60]]] ==5&];CTriplesDown4=Select[primTriples,Total[Last/@FactorInteger[Abs[(1/60) #1[[1]] #1[[2]]#1[[3]]]]]≥6&];

C4tripsDown=Graphics3D[{RGBColor[0.6`,0.6`,1],PointSize[ptSize],Point[CTriplesDown4]}];C2tripsDown=Graphics3D[{Red,PointSize[ptSize],Point[CTriplesDown2]}];
C3tripsDown=Graphics3D[{Black,PointSize[ptSize],Point[CTriplesDown3]}];

pic=Show[{coneBot,C4tripsDown,C2tripsDown,C3tripsDown},ViewPoint->{3,1,1}]

(* 
	The fact that there are infinitely many "black" points above is a consequence of Green-Tao's theorem. 
The infinitude of "red" points above would follow from Schinzel's hypothesis.
*)

(* Now for a thin orbit *)

(* The Parameters: *)

(* Group Gamma is thin, generated by these matrices *)
a={{2,1},{3,2}};
b={{1,2},{2,5}};
c={{2,3},{1,2}};
d={{5,2},{2,1}};
ai=Inverse[a];
bi=Inverse[b];
ci=Inverse[c];
di=Inverse[d];
Smat={a,ai,b,bi,c,ci,d,di};

inv[j_]:=If[Mod[j,2]==0,j-1,j+1]; (* Location of inverse matrix *)

levelmax=5; (* How far in the wordlength metric to go in Gamma *)

(* Base point for the orbit *)
u={2,1};


(* The Code: *)

numMats=Length[Smat];

m={{1,0},{0,1}};
matrices={m};
vect[m_]:=u.m;
triple[u_]:={u[[1]]^2-u[[2]]^2,2u[[1]] u[[2]], -u[[1]]^2-u[[2]]^2};
norm[m_]:=m[[1,1]]^2+m[[1,2]]^2+m[[2,1]]^2+m[[2,2]]^2;
tr[mat_]:=Abs[mat[[1,1]]+mat[[2,2]]];

vec3[v_]:=v[[3]];
gcd[v_]:=GCD[v[[1]],v[[2]],v[[3]]];

moremat[from_,level_]:=If[level<levelmax,
Module[{j},
For[j=1,j≤numMats,j++,
If[from≠j,
m=Smat[[j]].m;
AppendTo[matrices,m];
moremat[inv[j],level+1];
m=Smat[[inv[j]]].m;
]
]
]
];

(* Run the Code *)

moremat[0,1];
vectors=triple /@ (vect/@ matrices);

(* Sort triples according to number of prime factors of xyz/60 *)

CTriplesDown2O = Select[vectors, Total[Last/@ FactorInteger[ Abs[#[[1]]  #[[2]] #1[[3]]/60]]]≤ 4&];
CTriplesDown3O = Select[vectors, Total[Last/@ FactorInteger[ Abs[#[[1]]  #[[2]]#1[[3]] /60]]] ==5&];CTriplesDown4O=Select[vectors,Total[Last/@FactorInteger[Abs[(1/60) #1[[1]] #1[[2]]#1[[3]]]]]≥6&];
C4tripsDownO=Graphics3D[{RGBColor[0.6`,0.6`,1],PointSize[ptSize],Point[CTriplesDown4O]}];C2tripsDownO=Graphics3D[{Red,PointSize[ptSize],Point[CTriplesDown2O]}];
C3tripsDownO=Graphics3D[{Black,PointSize[ptSize],Point[CTriplesDown3O]}];

pic=Show[{coneBot,C4tripsDownO,C2tripsDownO,C3tripsDownO} (*,ViewPoint->{3,1,1}*)]

(*
We conjecture that the number of "black" points is infinite (5 is the saturation number), while the number of reds is finite!
*)