#### Catalan Polygons

###################### functions #####################################



## packages and basic functions

# packages
 with(plots):
 with(numtheory):
 with(combinat):
 
# catalan numbers
 catalan := proc(K)
   RETURN(binomial(2*K,K)/(K+1));
   end:
 
 Ocatalan := proc(K)
   local e,n,fp,s;
   e := K mod 6;
   n := K+2;
   if (e = 3) or (e = 5) 
      then fp := catalan(K)
                 +n*catalan((K-1)/2);
      fi; 
   if (e = 1)  
      then fp := catalan(K)+2*(n/3)*catalan((K-1)/3)
                 +n*catalan((K-1)/2);
      fi; 
   if (e = 0) or (e = 2)  
      then fp := catalan(K)+(n/2)*catalan(K/2)
               + (n/2)*catalan(K/2)
               + (n/2)*(add(catalan(K/2-1-s)*catalan(s),s=0..K/2-1));
      fi; 
   if (e = 4)  
      then fp := catalan(K)+(n/2)*catalan(K/2)+2*(n/3)*catalan((K-1)/3)
                 +(n/2)*catalan(K/2)
                 +(n/2)*(add(catalan(K/2-1-s)*catalan(s),s=0..K/2-1));
      fi; 
     RETURN(fp/(2*n));
   end:
 
 

## basic polygon functions

# polygon functions
 NumTriangles := proc(HP)
   RETURN(CP[1]-2);
   end:
 
 NumVertices := proc(CP)
   RETURN(CP[1]);
   end:
 
 NumDiags := proc(CP)
   RETURN(nops(CP[2]));
   end:
   
 BoundaryEdges := proc(CP)
   local n,k;
   n := NumVertices(CP); 
   RETURN({seq({k,1+(k mod n)}, k=1..n)});
   end:
 
 Diags := proc(CP)
   RETURN(CP[2]);
   end:
 
 AllEdges := proc(CP)
   local i;
   RETURN({op(BoundaryEdges(CP)),op(Diags(CP))});
   end:
 
# geometric type 
 mu := proc(L)
   local x;
   RETURN(nops(L)-2-add(1/x,x=L));
   end:  
 
 hyperbolicpoly := proc(L)
   local m,vars,x,subsseq,mi;
   m := mu(L);
   vars := indets(m);
   subsseq := seq(1/x=0,x=vars); 
   mi := subs(x=x,subsseq,m);
   RETURN(evalb(mi > 0));
   end:  
 
## polygon drawing 
 DrawPoly := proc()
   local CP,labelopt,pic,labels,a,n,e,pnts,edges,diags,k;
   CP := args(1);
   n := NumVertices(CP);
   pnts := evalf([seq([cos(2*k*Pi/n),sin(2*k*Pi/n)],k=1..n)]);
   labelopt := false;
   for a in args[2..nargs] do
      if type(a,list) then pnts := a; fi;
      if type(a,relation) then labelopt := rhs(a); fi;
      od;
   edges := {seq([pnts[e[1]],pnts[e[2]]], e=BoundaryEdges(CP))};
   diags := {seq([pnts[e[1]],pnts[e[2]]], e=Diags(CP))};
   pic := plot({op(edges),op(diags)},
                view =[-1.2..1.2,-1.2..1.2],
                scaling=constrained,colour = black,axes=none);
   labels := textplot({seq([1.1*pnts[i][1],1.1*pnts[i][2],``.i],
                       i=1..n)},
                       view =[-1.2..1.2,-1.2..1.2],
                       scaling=constrained,colour = black,axes=none);
   if labelopt = true 
      then RETURN(display({pic,labels}));
      else RETURN(pic);
      fi;
   end: 
 
 DrawQuad := proc()
   local CP,QA,n,nd,pnts,k,j,QC,QS,e1,e2,s;
   CP := args[1];
   QA := args[2];
 
   n := NumVertices(CP);
   QC := [];
   QS := array(1..4,[0,0,0,0]);
   k := 0;
   for j from 1 to nops(QA) do 
     if QA[j] <> 1 
        then QC := [op(QC),j];
             k := k+1;
        else QS[(k+3 mod 4)+1] := QS[(k+3 mod 4)+1]+1;
        fi; 
     od; 
   QS[4] := n-4-add(QS[k],k=1..3);
   pnts:= array(1..n);
   pnts[QC[1]] := [1,1];
   pnts[QC[2]] := [-1,1]; 
   pnts[QC[3]] := [-1,-1];
   pnts[QC[4]] := [1,-1];
   for j from 1 to 4 do
     s := QS[j];
     e1 := pnts[QC[j]];
     e2 := pnts[QC[(j mod 4)+1]];
     for k from 1 to s do
         pnts[((QC[j]+k+n-1) mod n)+1] := 
                [(e1[1]*(s+1-k)+e2[1]*k)/(s+1),
                 (e1[2]*(s+1-k)+e2[2]*k)/(s+1)]; 
        od;
     od;   
  DrawPoly(CP,convert(pnts,list),args[3..nargs]);
  end: 
## polygon construction and transformation 
# relabelling
 TranslateDiags := proc(diags,tab)
   local ldiags,e,newP;
   ldiags := [];
   for e in diags do;
       ldiags :={op(ldiags),{tab[e[1]],tab[e[2]]}};
       od;
   RETURN(ldiags);
   end:
 
# construction
 AttachTriangle := proc(CP,k)
   local n,ldiags,tab,j;
   n := NumVertices(CP);
   tab := [seq(j,j=1..k),seq(j+1,j=k+1..n)];
   ldiags := {op(TranslateDiags(Diags(CP),tab)), 
                 {k,((k+1) mod (n+1))+1}};
   RETURN([n+1,ldiags]);   
   end:
 
# dihedral action
 DihedralTable:= proc(n)
   local k,i,rots,refls;
   rots := [seq([seq(((i+k-1) mod n)+ 1, i=1..n)],k=0..n-1)];
   refls := [seq([seq(((-i+k-1) mod n)+ 1, i=1..n)],k=0..n-1)];
   RETURN([op(rots),op(refls)]);
   end:
 
 PermutePoly := proc(CP,tab)
   RETURN([NumVertices(CP),TranslateDiags(Diags(CP),tab)]);
   end:
  
 DihedralOrbit:= proc(CP,DHtab)
   RETURN({seq(PermutePoly(CP,tab),tab=DHtab)});
    end:
 
 DHOrbitSpace:= proc(CPL)
   local OS,LCPL,DHtab,n,CP,DHO;
   LCPL := {op(CPL)};
   n :=  CPL[1][1];
   OS := [];
   DHtab := DihedralTable(n);
   while nops(LCPL) > 0 do
     CP := LCPL[1];
     OS := [op(OS),CP];
     DHO := DihedralOrbit(CP,DHtab);
     LCPL := LCPL minus DHO;
     od;     
   RETURN(OS);
   end:
 
## polygon analysis 
# vertex types
 ConnectedVertices := proc(edgeset,v)
   local VL,e;
   VL := {};
   for e in edgeset do
       if member(v,e) then VL := VL union e; fi;
       od;
   RETURN(VL minus {v});
   end:
 
 ThirdVertexSet := proc(edgeset,v1,v2)
   if not member({v1,v2},edgeset) then RETURN({}); fi;
   RETURN(ConnectedVertices(edgeset,v1) intersect
          ConnectedVertices(edgeset,v2));
   end:
 
 VertexTypes := proc(CP)
   local i,j,k,n,verts,vtypes,t,u,v,w,TVS,es;
   n := NumVertices(CP);
   es := AllEdges(CP);
   vtypes := array(1..n);
   vtypes[1]:= 1; vtypes[2] := 2;
   verts := {1,2};
   while nops(verts) < n do; 
         for u in verts do
         for v in verts do
             TVS := ThirdVertexSet(es,u,v);
             t := 6-vtypes[u]-vtypes[v];
             for w in TVS do 
                 vtypes[w] := t;
                 verts := {op(verts),w}; 
                od;
             od;od;
         od;               
    RETURN([seq(vtypes[k],k=1..n)]);
   end:
 
 VertexTypeMults := proc(CP)
   local AML,VTL,s,i,j,k;
   AML := AnglesinCornerList(CP);
   VTL := VertexTypes(CP);
   s := array(1..3,[[],[],[]]); 
   for i from 1 to NumVertices(CP) do
      j := VTL[i];
      s[j] := [op(s[j]),AML[i]]; 
      od;
   RETURN([seq(s[k],k=1..3)]);
   end:
 
# angle functions
 AnglesinCorner := proc(CP,n)
   local diag,s;
   s := 1;
   for diag in Diags(CP) do
       if member(n,diag) then s := s+1; fi;
       od;
   RETURN(s);
   end:
 
 AnglesinCornerList := proc(CP)
   local k;
   RETURN([seq(AnglesinCorner(CP,k),k=1..NumVertices(CP))]);
   end:
 
 VertexAngles := proc(CP,L)
   local VTL,ACL,i;
   VTL := VertexTypes(CP);
   ACL := AnglesinCornerList(CP); 
   RETURN([seq(ACL[i]/L[VTL[i]],i=1..NumVertices(CP))]);
   end:
 
# angle multiplicities and corner constraints
 MinMaxCorners := proc(ML)
   local L,G,m,t;
   L := lcm(op(ML));
   if L=1 then RETURN([nops(ML)]); fi;
   if L=G then RETURN([0,nops(ML)]); fi;
   m := 0;
   for t in ML do if t<L then m := m+1; fi; od;
   if m = nops(ML) 
      then RETURN([nops(ML)]);
      else RETURN([m,nops(ML)]); 
      fi;
    end:
 
 MinPolyCorners := proc(CP)
   local VTL,ML,L,i;
   VTL := VertexTypeMults(CP);
   L := [seq(MinMaxCorners(ML),ML=VTL)];
   RETURN(add(L[i][1],i=1..nops(L)));
   end:
 
# compatibilty and validity functions
 ValidPolys := proc(CP,parms)
   local lcmL,VTL,MML,VA,ML,L,a,i,j,k,stuv,polypairs,lmn,QA;
   VTL := VertexTypeMults(CP);
   MML := [seq(MinMaxCorners(ML),ML=VTL)];
   lcmL := [seq(lcm(op(ML)),ML=VTL)];
   lmn := array(1..3);
   polypairs := [];
   for i from 1 to nops(MML[1]) do
   for j from 1 to nops(MML[2]) do
   for k from 1 to nops(MML[3]) do
     if (MML[1][i]+MML[2][j]+MML[3][k] = 4)
        then
             if i = nops(MML[1])
                then lmn[1] := lcmL[1]*parms[1]; 
                else lmn[1] := lcmL[1] 
                fi; 
             if j = nops(MML[2]) 
                then lmn[2] := lcmL[2]*parms[2]; 
                else lmn[2] := lcmL[2] 
                fi;
             if k = nops(MML[3]) 
                then lmn[3] := lcmL[3]*parms[3]; 
                else lmn[3] := lcmL[3] 
                fi; 
             VA := VertexAngles(CP,lmn);
             stuv := [];
             for a in VA do if a<>1 then stuv :=[op(stuv),1/a]; fi od;
             if (nops(stuv) = 4) and hyperbolicpoly(stuv)
                then QA := simplify(subs(parms[1]=lmn[1],parms[2]=lmn[2],
                                parms[3]=lmn[3],VA));
                     polypairs := [op(polypairs),[[lmn[1],lmn[2],lmn[3]],
                                  stuv,QA]];
                     fi;
             fi;
      od;od;od;
   RETURN(polypairs);
   end:


###################### driver ################

restart;
   
# main loop - polygon creation
 top := top;
 [seq([K,catalan(K),Ocatalan(K)],K=1..top)];
 polys.1 := {[3,{}]}:
 orbspa.1 := {[3,{}]}:
 for K from 2 to top do 
   print(`------------------------------------------------`);
   print(`# triangs ` = K);
   st := time():
   polys.K := [];
   orbspa.K := [];
   for CP in orbspa.(K-1) do
     for k from 1 to K+1 do 
       polys.K := {op(polys.K),AttachTriangle(CP,k)}   
       od;
     od;
    print(`# polys found ` = nops(polys.K),
          `time taken `= time()-st, `catalan # ` =catalan(K));
   st := time():
   orbspa.K := DHOrbitSpace(polys.K);
   save orbspa.K,polyorbs.K.`.mpl`; 
   print([`# orbits found `= nops(orbspa.K),
           `time taken `= time()-st,`Ocatalan # ` = Ocatalan(K)]);
   od: 
 
# main loop - quadrilateral search
 [seq([k,Ocatalan(k)],k=1..12)];
print("finding quadrilaterals");
 otop := top;
 unassign('l','m','n');
 for K from 2 to otop do
   read polyorbs.K.`.mpl`;
   OS := orbspa.K:    
 print(`------------------------------------------------`);
   print(`# triangs` = K);
   VPL :=[];
   for CP in OS do 
        VP := ValidPolys(CP,[l,m,n]);
        if nops(VP) > 0 then
           VTL := VertexTypes(CP);
           print(VTL);
           VPL := [op(VPL),seq([CP,op(x)],x=VP)];
        for P in VP do  
           print(P);
           od;
        fi;
      od;
   save VPL, hquads.K.`.mpl`;
   od:


quit;
#