//Magma package of functions for classifying tilings 
// by S. Allen Broughton


//************** general functions

OrbitReps := function(G,X);
   OX := Orbits(G,X);
   return {@Rep(O): O in OX @};
   end function;

OrbitSizes := function(G,X);
   OX := Orbits(G,X);
   return [#O: O in OX ];
   end function;


SumList := function(L);
   S := 0;
   for x in L do
      S := S+x;
      end for;
   return S;
   end function;

HomImage := function(f);
   return sub<Codomain(f)| f(Generators(Domain(f)))>;
   end function;


//*********** String functions

Strip :=  function(s,t);
   u := "";
   for i := 1 to #s do
       if Position(t,s[i]) eq 0
          then u := u*s[i];
          end if;
       end for;
   return u;
   end function;


ParseBData :=  function(s);
   t := Strip(s," ");
   t := Split(s,"[]");
   if #t lt 4 then t[4] := ""; end if;
   u := Split(t[2],", ");
   sig := StringToInteger(u[1]);
   g := StringToInteger(u[2]);
   ords := Split(t[3],", ");
   ords := [StringToInteger(c): c in ords];
   note := Strip(t[4],",\" ");
   return sig, g, ords, note;
   end function;

//************** compatibility test
RHtest := function(sig,g,l,m,n);
  if 2*sig-2 eq g*(1-1/l-1/m-1/n)
     then return true;
     else return false;
     end if;
  end function;

//******** group constructions
// construct product of permutation reprsentations of a group
// and the "block diagonal" normalizer and centralizer
ProductPermRep := function(G,flist);
   Glist := [Domain(f): f in flist];
   ngens := Ngens(G);

   Ggens := [G.i : i in [1..ngens]];
   degs := [Degree(Codomain(f)): f in flist];
   td := SumList(degs);
   BSG := SymmetricGroup(td);
   T := [];
   rs := 0;
   for d in degs do
      L := [((j+rs-1) mod td)+1: j in [1..td]];
      Append(~T,BSG!L);
      rs := rs+d;
      end for;
   PGgens := [];
   homlist := [];
   for i in [1..ngens] do
       fg := Id(BSG);
       for j in [1..#flist] do;
            f := flist[j];
            g := Glist[j].i;
            fg :=  fg*(BSG!f(g))^T[j];
           end for;
       Append(~PGgens,fg);
       Append(~homlist,fg);
       end for;
   f := hom<G -> BSG | homlist>;
   P := sub<BSG|PGgens>;
   if #P eq #G
      then fi := hom<P -> G | Ggens>;
      else fi := "fail";
      end if;
   NPgens := [];
   ZPgens := [];
   for j in [1..#flist] do;
       L := HomImage(flist[j]);
       S := SymmetricGroup(degs[j]);
       N := Normalizer(S,S!!L);
       Z := Centralizer(S,S!!L);
       for g in Generators(N) do Include(~NPgens,(BSG!g)^T[j]); end for;
       for g in Generators(Z) do Include(~ZPgens,(BSG!g)^T[j]); end for;
       end for;
   BNP := sub<BSG|NPgens>;
   NP := Normalizer(BNP,P);
   BZP := sub<BSG|ZPgens>;
   ZP := Centralizer(BZP,P);
   return f,P,fi,NP,ZP;
   end function;

SylowPermRep := function(G,p);
   H := Sylow(G,p);
   f,P,K := CosetAction(G,H);
   n :=  Degree(P);
   S := SymmetricGroup(n);
   NP := Normalizer(S,P);
   return f,P,K,NP;
   end function;

HallPermRepAtp := function(G,p);
   H := HallSubgroup(G,-p);
   f,P,K := CosetAction(G,H);
   n :=  Degree(P);
   S := SymmetricGroup(n);
   NP := Normalizer(S,P);
   return f,P,K,NP;
   end function;

HallPermRep := function(G);
   prL := PrimeDivisors(#G);
   flist := [* *];
   for p in prL do
     f,P := HallPermRepAtp(G,p);
     Append(~flist,f);
     end for;
   f,P,fi,NP,ZP := ProductPermRep(G,flist);
   return f,P,fi,NP,ZP;
   end function;

BlockNormalCentralizer := function(G);
//it is assumed that the natural GSet for G is 1..n for some n.
//the "block diagonal" normalizer and centralizer are returned (does not switch orbits)
   td := Degree(G);
   BSG := SymmetricGroup(td);
   NPgens := [];
   ZPgens := [];
   OL := Orbits(G);
   for O in OL do 
      f,M := Action(G,O);
      S := Sym(O);
      N := Normalizer(S,M);
      Z := Centralizer(S,M);
      NPgens := NPgens cat [BSG!g : g in Generators(N)];
      ZPgens := ZPgens cat [BSG!g : g in Generators(Z)];
      end for;
   BNP := sub<BSG|NPgens>;
   NG := Normalizer(BNP,G);
   BZP := sub<BSG|ZPgens>;
   ZG := Centralizer(BZP,G);
   return NG,ZG;
   end function;



//************* finding elements of a given order
OrderClassReps := function(G,n)
   KG := Classes(G);
   R := {@ @};
   for K in KG do
      if K[1] eq n
         then Include(~R,K[3]);
         end if;
      end for;
   return R;
   end function;

OrderClass := function(G,n)
   R := OrderClassReps(G,n);
   OC := {@ @};
   for g in R do
      OC := OC join {@ h : h in g^G @};
      end for;
   return OC;
   end function;

// NPG is a group containing G as a normal subgroup
OrderClassRepsbyNPG := function(G,NPG,n)
   R := OrderClassReps(G,n);
   RAG := {@ @};
   for g in R do
      found := false;
      for h in RAG do
          if g in h^NPG
             then found := true;
             break;
             end if;
          end for;
          if not found
             then Include(~RAG,g);
             end if;
       end for;
   return RAG;
   end function;


// H is a subgroup of G
OrderClassRepsbyH := function(G,H,n)
   Kn := OrderClass(G,n);
   OC :={@ @};
   while #Kn gt 0 do
     Include(~OC,Kn[1]);
     Kn := Kn diff {@ x: x in Kn[1]^H @};
     end while;
   return OC;
   end function;


// ************* finding automorphisms
ThetaQuest := function(G,a,b,ta,tb)
   trivth := Id(SymmetricGroup(1));
   trivX := {@ @};
// basic test
   if not (Order(a) eq Order(ta)) then return trivth,trivX; end if;
   if not (Order(b) eq Order(tb)) then return trivth,trivX; end if;
// powers of generators
   aP := {@ <a^i,ta^i> : i in [1..(Order(a)-1)] @};
   bP := {@ <b^i,tb^i> : i in [1..(Order(b)-1)] @};
//create the graph of theta
   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;
// create image sequences for th, La, Lb, Lta and Ltb
   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 @};
// test for 1-1 property
   if not (#X eq #G) then return trivth,trivX; end if;
   if not (#thX eq #G) then return trivth,trivX; end if;
// transform to standard form
   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];
// compute cycles in standard form
   S := SymmetricGroup(#G);
   th := S!iY;
   La := S!iLa;
   Lb := S!iLb;
   Lta := S!iLta;
   Ltb := S!iLtb;
// test and return
   if (La^th eq Lta) and (Lb^th eq Ltb)
      then return th,X;
      else return trivth,trivX;
      end if;
  end function;

IsKaleidoscopic := function(G,tri);
   a := tri[1];
   b := tri[2];
   th,X := ThetaQuest(G,a,b,a^(-1),b^(-1));
   return (#X gt 0);
   end function;

//*********** finding triples
lmnTriplesbyNPG := function(G,NPG,l,m,n);
   OKl := OrderClassRepsbyNPG(G,NPG,l);
   triset :={@ @};
   for a in OKl do
      H := Centralizer(NPG,a);
      OKm := OrderClassRepsbyH(G,H,m);
      for b in OKm do
          c := (a*b)^(-1);
          if Order(c) eq n
            then if #sub<G|a,b> eq #G
                    then Include(~triset,<a,b,c>);
                 end if;
            end if;
          end for;
      end for;
    return triset;
    end function;

lmnTriplesbyAutG:= function(G,triset);
  triset2 := {@ @};
  for t in triset do
      found := false;
      for u in triset2 do;
         th,X := ThetaQuest(G,t[1],t[2],u[1],u[2]);
         if #X gt 0
           then found := true;
           break;
           end if;
         end for;
      if not found
         then Include(~triset2, t);
         end if;
      end for;
  return triset2;
  end function;


NumlmPairs := function(G,NPG,l,m);
   OKl := OrderClassRepsbyNPG(G,NPG,l);
   S := 0;
   for x in OKl do
     H := Centralizer(NPG,x);
     S := S + #OrderClassRepsbyH(G,H,m);
     end for;
   return S;
   end function;

//************* equivalence of triples
ConformalEquivalence := function(G,T1,T2);
   a1 := T1[1]; b1 := T1[2]; c1 := T1[3];
   a2 := T2[1]; b2 := T2[2]; c2 := T2[3];
   l := Order(a1);
   m := Order(b1);
   n := Order(c1);
   if (l eq m) and (m eq n)
      then a3 := b1; b3 := c1; c3 :=a1;
      th,X := ThetaQuest(G,a2,b2,a3,b3);
      if #X gt 0 then
         return [2,3,1],th,X;
         end if;
      end if;
   if (l eq m) and (m eq n)
      then a3 := c1; b3 := a1; c3 :=b1;
      th,X := ThetaQuest(G,a2,b2,a3,b3);
      if #X gt 0 then
         return [3,1,2],th,X;
         end if;
      end if;
   if l eq m
      then a3 := b1; b3 := a1; c3 :=c1^a1;
      th,X := ThetaQuest(G,a2,b2,a3,b3);
      if #X gt 0 then
         return [2,1,3],th,X;
         end if;
      end if;
   if l eq n
      then a3 := c1; b3 := b1^c1; c3 :=a1;
      th,X := ThetaQuest(G,a2,b2,a3,b3);
      if #X gt 0 then
         return [3,2,1],th,X;
         end if;
      end if;
   if m eq n
      then a3 := a1^b1; b3 := c1; c3 :=b1;
      th,X := ThetaQuest(G,a2,b2,a3,b3);
      if #X gt 0 then
         return [1,3,2],th,X;
         end if;
      end if;

      a3 := a1; b3 := b1; c3 :=c1;
      th,X := ThetaQuest(G,a2,b2,a3,b3);
      if #X gt 0 then
         return [1,2,3],th,X;
         end if;

    return [],th,X;
    end function;

AntiConformalEquivalence := function(G,T1,T2);
   a1 := T1[1]; b1 := T1[2]; c1 := T1[3];
   a2 := T2[1]; b2 := T2[2]; c2 := T2[3];
   l := Order(a1);
   m := Order(b1);
   n := Order(c1);
   if (l eq m) and (m eq n)
      then a3 := b1^(-1); b3 := c1^(-1); c3 :=c1*a1^(-1)*c1^(-1);
      th,X := ThetaQuest(G,a2,b2,a3,b3);
      if #X gt 0 then
         return [2,3,1],th,X;
         end if;
      end if;
   if (l eq m) and (m eq n)
      then a3 := c1^(-1); b3 := a1^(-1); c3 :=a1*b1^(-1)*a1^(-1);
      th,X := ThetaQuest(G,a2,b2,a3,b3);
      if #X gt 0 then
         return [3,1,2],th,X;
         end if;
      end if;
   if l eq m
      then a3 := b1^(-1); b3 := a1^(-1); c3 :=c1^(-1);
      th,X := ThetaQuest(G,a2,b2,a3,b3);
      if #X gt 0 then
         return [2,1,3],th,X;
         end if;
      end if;
   if l eq n
      then a3 := c1^(-1); b3 := b1^(-1); c3 :=a1^(-1);
      th,X := ThetaQuest(G,a2,b2,a3,b3);
      if #X gt 0 then
         return [3,2,1],th,X;
         end if;
      end if;
   if m eq n
      then a3 := a1^(-1); b3 := c1^(-1); c3 :=b1^(-1);
      th,X := ThetaQuest(G,a2,b2,a3,b3);
      if #X gt 0 then
         return [1,3,2],th,X;
         end if;
      end if;

      a3 := a1^(-1); b3 := b1^(-1); c3 :=b1*c1*b1^(-1);
      th,X := ThetaQuest(G,a2,b2,a3,b3);
      if #X gt 0 then
         return [1,2,3],th,X;
         end if;

    return [],th,X;
    end function;

IsometricClassesoflmnTriples := function(G,NPG,l,m,n);
   triples := lmnTriplesbyNPG(G,NPG,l,m,n);
   print "found", #triples, "NPG class(es) of triples";
   triples := lmnTriplesbyAutG(G,triples);
   print "found", #triples, "Aut(G) class(es) of triples";
   triples2 := {@ @};
   for t in triples do
      found := false;
      for u in triples2 do;
         pi,th,X := ConformalEquivalence(G,t,u);
         if #pi gt 0
           then found := true;
           break;
           end if;
         pi,th,X := AntiConformalEquivalence(G,t,u);
         if #pi gt 0
           then found := true;
           break;
           end if;
         end for;
      if not found
         then Include(~triples2, t);
         end if;
      end for;
  print "found", #triples2,  "isometric class(es) of triples";
  return triples2;
  end function;


//********** Finding triples for specific groups

//**Cyclic groups
CyclicGroupTripleNormalForm := function(g,tri);
   G := SymmetricGroup(g);
   a := G!tri[1];
   b := G!tri[2];
   c := G!tri[3];
   if Order(c) eq g
     then d := c;
     else e := Gcd(Order(b),Order(c));
          d := (b^e)*c;
     end if;
   oa := 0;
   repeat oa := oa+1; until (a eq d^oa) or (oa eq g);
   ob := 0;
   repeat ob := ob+1; until (b eq d^ob) or (ob eq g);
   oc := 0;
   repeat oc := oc+1; until (c eq d^oc) or (oc eq g);
   return <oa,ob,oc>;
   end function;

CyclicGroupTriples := function(sig,g,ords);
   l := ords[1]; m := ords[2]; n := ords[3];
   G := CyclicGroup(g); PG :=SymmetricGroup(g);
   NPG := Normalizer(PG,G);
   print "bdata: ", <g,sig,ords>, "|NPG| = ", #NPG; 
   triset := IsometricClassesoflmnTriples(G,NPG,l,m,n);
   normforms := [<CyclicGroupTripleNormalForm(g,tri),true>: tri in triset];
   return triset,normforms;
   end function;

//**abelian groups
AbelianPairs := function(n);
   D := Divisors(n);
   P := {@ @};
   for d in D do
      q := n div d;
      s := Gcd(q,d);
      Include(~P,<s,n div s>);
      end for;
   return P;
   end function;

ProperAbelianPairs := function(n);
   AP := AbelianPairs(n);
   PAP := {@ x : x in AP | x[1] gt 1 @};
   return PAP;
   end function;

AbelianGroups2Gen := function(n);
   AL := {@ @};
   PAP := ProperAbelianPairs(n);
   for p in PAP do
     Include(~AL, AbelianGroup([p[1],p[2]]));
     end for;
   return AL;
   end function;

AbelianGroup2GenTripleNormalForm := function(A,tri);
   return <ElementToSequence(tri[1]),ElementToSequence(tri[2]),ElementToSequence(tri[3])>;
   end function;

AbelianGroup2GenTriples := function(A,sig,g,ords);
   l := ords[1]; m := ords[2]; n := ords[3];
   inv := Invariants(A);
   P1 := CyclicGroup(inv[1]);
   P2 := CyclicGroup(inv[2]);
   f1 := hom<A -> P1 | [P1.1,Id(P1)]>;
   f2 := hom<A -> P2 | [Id(P2),P2.1]>;
   f,G,fi,NPG := ProductPermRep(A,[f1,f2]);
   print "bdata: ", <g,sig,ords>, "|NPG| = ", #NPG; 
   triset := IsometricClassesoflmnTriples(G,NPG,l,m,n);
   normforms :=[];
   for tri in triset do
      fitri := <fi(tri[1]),fi(tri[2]), fi(tri[3])>;
      Append(~normforms,<AbelianGroup2GenTripleNormalForm(A,fitri),true>);
      end for;
   return triset,normforms;
   end function;


//**polycyclic groups and pgroups
PolycyclicTripleNormalForm := function(tri);
   return <ElementToSequence(tri[1]),ElementToSequence(tri[2]),ElementToSequence(tri[3])>;
   end function;

// the following script is experimental and does not work well
PolycyclicGroupTriples := function(G,sig,g,ords);
   l := ords[1]; m := ords[2]; n := ords[3];
   f,P,fi,NPG,ZPG := HallPermRep(G);
   print "bdata: ", <g,sig,ords>, "|NPG| = ", #NPG; 
   triset := IsometricClassesoflmnTriples(P,NPG,l,m,n);
   normforms :=[];
   for tri in triset do
      fitri := <fi(tri[1]),fi(tri[2]), fi(tri[3])>;
      Append(~normforms,<PolycyclicTripleNormalForm(fitri),IsKaleidoscopic(P,tri)>);
      end for;
   return triset,normforms;
   end function;

PolycyclicGroupTriples2 := function(G,sig,g,ords);
   l := ords[1]; m := ords[2]; n := ords[3];
   NPG := G;
   print "bdata: ", <g,sig,ords>, "|NPG| = ", #NPG; 
   triset := IsometricClassesoflmnTriples(G,NPG,l,m,n);
   normforms :=[];
   for tri in triset do
      Append(~normforms,<PolycyclicTripleNormalForm(tri),IsKaleidoscopic(G,tri)>);
      end for;
   return triset,normforms;
   end function;

Pgroup3TripleNormalForm := function(tri);
   return <ElementToSequence(tri[1]),ElementToSequence(tri[2]),ElementToSequence(tri[3])>;
   end function;

Pgroup3Triples := function(G,sig,g,ords);
   l := ords[1]; m := ords[2]; n := ords[3];
   NPG := G;   
   print "bdata: ", <g,sig,ords>," |NPG| = ", #NPG; 
   triset := IsometricClassesoflmnTriples(G,NPG,l,m,n);
   normforms :=[];
   for tri in triset do
      Append(~normforms,<PolycyclicTripleNormalForm(tri),IsKaleidoscopic(G,tri)>);
      end for;
   return triset,normforms;
   end function;

// the following script is experimental and does not work well
//Permutation Groups - Non-solvable groups
PermGroupTriples := function(G,sig,g,ords);
   l := ords[1]; m := ords[2]; n := ords[3];
   d := Degree(G);
   NPG := BlockNormalCentralizer(G);  
   print "bdata: ", <g,sig,ords>, "|NPG| = ", #NPG; 
   triset := IsometricClassesoflmnTriples(G,NPG,l,m,n);
   normforms :=[];
   for tri in triset do
      Append(~normforms,<tri,IsKaleidoscopic(G,tri)>);
      end for;
   return triset, normforms;
   end function;

PermGroupTriples2 := function(G,sig,g,ords);
   l := ords[1]; m := ords[2]; n := ords[3];
   d := Degree(G);
   NPG := G;  
   print "bdata: ", <g,sig,ords>, "|NPG| = ", #NPG; 
   triset := IsometricClassesoflmnTriples(G,NPG,l,m,n);
   normforms :=[];
   for tri in triset do
      Append(~normforms,<tri,IsKaleidoscopic(G,tri)>);
      end for;
   return triset, normforms;
   end function;


//experimental only not used
MetaCyclic := function(p,q,k);
   if not ((k^p mod q) eq 1)
      then return {@ @};
      end if;
   H := SymmetricGroup(q);
   b := H![1+(i mod q): i in [1..q]] ;
   t,a := IsConjugate(H,b,b^k);
   G := SymmetricGroup(q+p);
   a := G!a;
   b := G!b;
   c := G!([i: i in [1..q]] cat [1+q+(i mod p): i in [1..p]]);
   return(sub<G|b,a*c>);
   end function;