BeginPackage["TriangleSequences`"];

T::usage = "T[{a,b}] gives the result of the triangle 
	mapping on the point {a,b}.";
Triangle::usage = "Triangle[{a,b}] gives the number
	of the partition triangle where the point {a,b}
	is located.";
TriangleSequence::usage ="TriangleSequence[{a,b},n]
	returns a list containing the first n terms in
	the triangle sequence for the point {a,b}.";
CVectors::usage= "CVectors[{a,b},n] gives you n vectors which
    approach the plane x+ay+cz=0.";
PlotOrbit::usage = "PlotOrbit[{a,b}] plots the orbit of
	the point {a,b} under the triangle iteration.";
NumberPoints::usage = "NumberPoints is an option for
	PlotOrbit.  NumberPoints->True numbers the points
	in the orbit, beginning with 0. Default setting:
	NumberPoints->True.";
DisplayLines::usage = "DisplayLines is an option for
	PlotOrbit.  DisplayLines->True draws lines
	between the points in the orbit. Default setting:
	DisplayLines->False.";
PartitionColor::usage = "PartitionColor is an option
	for PlotOrbit. PartitionColor specifies the color
	used to display the partition lines of the triangle.
	Default setting: PartitionColor->RGBColor[1,.5,0].";
OrbitColor::usage = "OrbitColor is an option for
	PlotOrbit.  OrbitColor specifies the color in which
	lines between orbit points will be displayed.
	Default setting: OrbitColor->RGBColor[0,.25,1].";
OrbitPoints::usage = "OrbitPoints is an option for
	PlotOrbit.  OrbitPoints->n displays the first n
	points in the orbit.  Default setting:
	OrbitPoints->30.";
PartitionLines::usage = "PartitionLines is an option
	for PlotOrbit.  PartitionLines->n displays the 
	first n-1 partition triangles.  Default setting:
	PartitionLines->10.";
FindPeriodicPoint::usage = "FindPeriodicPoint[{a1,a2,...}]
	finds a point with periodic triangle sequence
	{a1,a2,...,a1,a2,...}.  FindPeriodicPoint[{b1,b2,...},
	{a1,a2,...}] finds the point with triangle sequence
	{b1,b2,...,a1,a2,...,a1,a2,...}";
Spread::usage = "Spread[{a,b}] begins with a region
	surrounding the point {a,b}, and plots the effect
	of the iterative triangle mapping on points inside
	this region.  Spread[{a,b},n] plots n iterations.";

Begin["`Private`"];

K[a_,b_] := Floor[N[Simplify[(1-a)/b]]];

T[{a_,b_}] := Simplify[Apart[{b/a,(1-a-K[a,b]*b)/a}]];

Triangle[{a_,b_}]:= K[a,b];

TriangleSequence[{a_,b_},n_:50] := Module[{x,y,seq,c},
	(x = a; y = b; seq = {}; c = 0;
	While[N[y] != 0 && c < n, {seq = Append[seq,K[x,y]];
	{x,y} = T[{x,y}]; c = c+1}]; seq)];

CVectors[{aa_,bb_},n_:50] := Module[{lst,c2,c1,c0,i,ak,new},
  lst={}; ak=TriangleSequence[{aa,bb},n]; 
  c2={1,0,0}; c1={0,1,0}; c0={0,0,1};
  For[i=1, i<=Length[ak], i++,
   {new=c2-c1-ak[[i]]*c0,
    lst=Append[lst,new],
    c2=c1, c1=c0, c0=new}];
  lst];

Options[PlotOrbit]={NumberPoints->True,
	DisplayLines->False,PartitionColor->RGBColor[1,.5,0],
	OrbitColor -> RGBColor[0,.25,1],OrbitPoints->30,
	PartitionLines->10};
	
PlotOrbit[{a_,b_}, opts___]:=
	PlotOrb[{a,b},
	{NumberPoints, DisplayLines,PartitionColor,
	OrbitColor,OrbitPoints,PartitionLines}/.{opts}/.Options[PlotOrbit] ];

PlotOrb[{a_,b_}, {numpts_,lines_,tcol_,ocol_,
	pts_,lns_}] := 
	Module[{len,triangle,points,orbit},
	(len = Length[TriangleSequence[{a,b},pts]];
	triangle = {tcol,Line[{{1,0},{0,0},{1,1}}],
		Table[Line[{{1,0},{1/i,1/i}}],{i,1,lns}]};
	points = NestList[T,{a,b},len];
	orbit = {ocol,Line[points]};
	plot = Table[Point[points[[i]] ], {i,1,len}];
	labels = Table[Text[i-1,points[[i]]+{.018,.008} ],
		{i,1,len}];
	display = {triangle, plot};
	If[lines,display = Append[display,orbit]];
	If[numpts,display = Append[display,labels]];
	Show[Graphics[display,
		Axes->True, PlotRange->{{0,1},{0,1}},
		AspectRatio->1]])];
	
Points[x1_,x2_,y1_,y2_] := 
	Flatten[Table[{x1+(x2-x1)i/50,y1+(y2-y1)j/50},
 	{i,1,10},{j,1,10}],1];
	
Spread[{a_,b_},r_:10]:=
	Module[{p,xch,ych},
	(xch = .002;ych = .002;
	p = Points[a-xch,a+xch,b-ych,b+ych];
	Show[ Graphics[{RGBColor[1,0,.2],
		Table[Point[p[[i]]],
		{i,1,Length[p]}]}],
		Axes->True, PlotRange->{{0,1},{0,1}},
		AspectRatio->1];
	For[n=0,n<r,n=n+1,
		{p = Table[T[p[[i]]],{i,1,Length[p]}];
		Show[ Graphics[{RGBColor[1,0,.2],
		Table[Point[p[[i]]],
	 	{i,1,Length[p]}]}],
		Axes->True, PlotRange->{{0,1},{0,1}},
		AspectRatio->1]}])];

P[k_]:={{0,0,1},{1,0,-1},{0,1,-k}};

InTriangle[{a_,b_}]:=
	If[Im[a] != 0 || Im[b]!=0,False,
		If[ N[a] > 1 || N[b] > N[a] || N[b]<0,
		False,True]];

AlmostInTriangle[{a_,b_}] := InTriangle[N[{a,b}]];

FindPeriodicPoint[per_]:= Module[{m,sol},
	(m = IdentityMatrix[3];
	Do[m = m.P[per[[i]] ], {i,1,Length[per]}];
	sol = Solve[{1,x,y}.m == l*{1,x,y},{x,y},l];
	el = Table[{x,y} /. sol[[i]],{i,1,Length[sol]}];
	ans = Simplify[Flatten[Select[el,InTriangle]]];
	If[ans != {}, ans,
		Flatten[Select[Re[el],AlmostInTriangle]]])];

S[{x_,y_},k_]:= {1/(1+k x+y),x/(1+k x+y)};

FindPeriodicPoint[beg_,per_]:= Module[{pt},
	(pt = FindPeriodicPoint[per];
	Simplify[Fold[S,pt,Table[beg[[-i]],
	{i,1,Length[beg]}]]])];


End[]
EndPackage[]