// AutTest is modified from Dr. S. Allen Broughton's 
// ThetaQuest as part of tileclass.mgm

AutTest := function(G,a,b,ta,tb)
   	aP := {@ <a^i,ta^i> : i in [1..(Order(a)-1)] @};
   	bP := {@ <b^i,tb^i> : i in [1..(Order(b)-1)] @};
	
	W := {@<Id(G),Id(G)>@};
	W1 := W;
	W2 := W;
	P1 := aP;
	P2 := bP;
	while #W lt #G do
	     	W1 := {@ <g[1]*x[1],g[2]*x[2]> : g in P1, x in W1 @};
	     	W2 := {@ <g[1]*x[1],g[2]*x[2]> : g in P2, x in W2 @};
	     	W := W join (W1 join W2);
	    	P3 := P1;
		P1 := P2;
		P2 := P3;
	end while;

   	X   := {@ z[1]:z in W @};
   	thX := {@ z[2]:z in W @};
   	LaX := {@ a*x :x in X @};
   	LbX := {@ b*x :x in X @};
   	LtaX := {@ ta*x :x in X @};
   	LtbX := {@ tb*x :x in X @};

	if not (#X eq #G) then 
		return false; 
	end if;
   	if not (#thX eq #G) then 
		return false; 
	end if;

   	iY  := [Position(X,z):z in thX];
   	iLa := [Position(X,z):z in LaX];
   	iLb := [Position(X,z):z in LbX];
   	iLta := [Position(X,z):z in LtaX];
   	iLtb := [Position(X,z):z in LtbX];

   	S := SymmetricGroup(#G);
   	th := S!iY;
   	La := S!iLa;
   	Lb := S!iLb;
   	Lta := S!iLta;
   	Ltb := S!iLtb;


	if (La^th eq Lta) and (Lb^th eq Ltb) then 
		return true;
	else 
		return false;
      	end if;
end function;

// IsSplitS.mgm is a modification of the script IsSplitS.mgm  
// originally developed by J.M. Belk - RHIT REU 1997
// the script determines whether a tiling splits at a reflection R
// the full tiling group Gs and its generators p,q, and r must be supplied
// R is the reflection we are testing for separability, it is uusally one of p,q, or r.

IsSplitS := function(Gs,p,q,r,R)
  	Xold := {};
	X := {Id(Gs)};
   	repeat
        	bX  := X diff Xold;
        	Xold := X;
        	Xnew := {{x*p,x*q,x*r} diff {R*x} : x in bX}; // in triples and doubles 
        	Xnew := {xn : xn in yn,yn in Xnew};           // as single elements 
        	X := Xold join Xnew;
        until (#Xold eq #X) or (R in Xnew);
   	if (#X eq #Gs) or (R in Xnew) then
		return false;
      	else 
		return true;
      	end if;
   end function;

GroupS := function(G,l,m,n)

	aleph := [];
	beis := [];

  	K := {@ g : g in G @};
	for p in K do
		if Order(p) eq l then
			for q in K do
				if Order(q) eq m then
					if Order(q^-1*p^-1) eq n then
						if #G eq #sub<G | p, q> then
							aleph := Append(aleph,p);
							beis := Append(beis,q);
						end if;
					end if;
				end if;
			end for;
		end if;
	end for;

	Alpha := [aleph[1]];
	Beta := [beis[1]];

	for w in [2..#aleph] do
		automorphic := false;
		for k in [1..#Alpha] do
			if  automorphic then
				break k;
			end if;
			automorphic := AutTest(G,Alpha[k],Beta[k],aleph[w],beis[w]);
		end for;
		if not automorphic then
			Alpha := Append(Alpha,aleph[w]);
			Beta := Append(Beta,beis[w]);
		end if;
	end for;


	for w in [1..#Alpha] do
	
		a := Alpha[w];
		b := Beta[w];

		ia := a^-1;
		ib := b^-1;

		aP := {@ <a^k, ia^k> : k in [1..(l-1)] @};
		bP := {@ <b^k, ib^k> : k in [1..(m-1)] @};

		Q := {@<Id(G),Id(G)>@};
		Q1 := Q;
		Q2 := Q;
		W1 := aP;
		W2 := bP;
		while #Q lt #G do
			Q1 := {@ <g[1]*x[1],g[2]*x[2]> : g in W1, x in Q1 @};
			Q2 := {@ <g[1]*x[1],g[2]*x[2]> : g in W2, x in Q2 @};
		        Q := Q join (Q1 join Q2);
	        	W3 := W1;
			W1 := W2;
			W2 := W3;
		end while;
	
		X := {@ q[1] : q in Q @};
		thX := {@ q[2] : q in Q @};

		Z := CyclicGroup(2);
		
		D := {@ <g,z> : g in X, z in Z @};
	
		P := SymmetricGroup(#D);
	
		S := [];
		for d in D do
			if d[2] eq Id(Z) then
				Y := [ Position(D, <d[1]*y[1], d[2]*y[2]>) : y in D];
			else
				Y := [ Position(D, <d[1]*thX[Position(X,y[1])], d[2]*y[2]>) : y in D];
			end if;
			S := Append(S, P!Y);
		end for;		
		
		Gs := PermutationGroup< #D | S >;

		for g in Gs do
			if Order(g) eq 2 then
				for h in Gs do
					if Order(h) eq 2 then
						for i in Gs do
							if Order(i) eq 2 then
								Dp := D[Position(S,g)];
								Dq := D[Position(S,h)];
								Dr := D[Position(S,i)];
								if not (Dp[2] eq Id(Z)) and not (Dq[2] eq Id(Z)) then
									if  Dp[1]*thX[Position(X,Dq[1])] eq a and Dq[1]*thX[Position(X,Dr[1])] eq b then
										p := g;
										q := h;
										r := i;
										break g;
									end if;
								end if;
							end if;
						end for;
					end if;
				end for;
			end if;
		end for;
		
		if IsSplitS(Gs,p,q,r,p) then
			//print "The tiling seperates along the p reflection";
		else
			//print "The tiling does not seperate along the p reflection";
		end if; 

		if IsSplitS(Gs,p,q,r,q) then
			//print "The tiling seperates along the q reflection";
		else
			//print "The tiling does not seperate along the q reflection";
		end if; 	
	
		if IsSplitS(Gs,p,q,r,r) then
			//print "The tiling seperates along the r reflection";
		else
			//print "The tiling does not seperate along the r reflection";
		end if; 

	end for;

	return <#G, #Alpha, "S">;
end function;



StarGroupS := function(GstarList)
   for GL in GstarList do
	Gs := GL[1];
	p := GL[2];
	q := GL[3];
	r := GL[4];

	if IsSplitS(Gs,p,q,r,p) then
		print "The tiling seperates along the p reflection";
	else	
		print "The tiling does not seperate along the p reflection";
	end if; 
	if IsSplitS(Gs,p,q,r,q) then
		print "The tiling seperates along the q reflection";
	else
		print "The tiling does not seperate along the q reflection";
	end if; 	
	if IsSplitS(Gs,p,q,r,r) then
		print "The tiling seperates along the r reflection";
	else
		print "The tiling does not seperate along the r reflection";
	end if; 
   end for;
return "IsSplitS Method -- Star";
end function;