
(* started by MK 2018-02-06 *)

(* counting how many ways there are to place zero or more self-avoiding polygons
   on a mxn torus, without that they share edges (but they may share vertices)

   such an arrangement is the same as a subgraph of the torodial lattice in which 
   every vertex has exactly two edges attached to it

   examples: 

   In[1]:= SAPRatfun[1]
   Out[1]= (1 - e*v^2 - e^2*v^2)/((-1 + e*v)*(-1 + v + e*v))

   In[2]:= Series[%, {v,0,3}]
   Out[2]= 1 + (1 + 2*e)*v + (1 + 2*e + 2*e^2)*v^2 + (1 + 3*e + 3*e^2 + 2*e^3)*v^3 + O[v]^4

   In[3]:= SAPSeries[1, 3]
   Out[3]= 1 + (1 + 2*e)*v + (1 + 2*e + 2*e^2)*v^2 + (1 + 3*e + 3*e^2 + 2*e^3)*v^3

 *)

MatrixTrace[m_] := Sum[m[[i,i]], {i, 1, Length[m]}]

SAPAlphabet[m_] :=
  Module[ {omega},
     (* a bit-vector of length 2m encodes a building block consisting of m horizontal 
        and m vertical edge positions. the k-th bit corresponds to an edge conneting the
        k-th vertex to its horizontal neighbor (in the past). the (m+k)th bit corresponds
        to an edge connecting the k-th vertex and the (k-1)-th vertex (mod m). *)
     omega = Table[Int2Vec[k, 2m] /. -1 -> 0, {k, 0, 2^(2m) - 1}];
     (* Select[omega, Max[Table[#[[k]] + #[[m + k]] + #[[m + 1 + Mod[k,m]]], {k, 1, m}]] <= 2&] *)
     omega
  ];

PrettyPrint[block1_] :=
  Module[ {m},
    m = Length[block1]/2;
    Print /@ StringJoin@@@Transpose[
     {Take[block1, m] /. {1 -> "_", 0 -> " "}, Drop[block1, m] /. {1 -> "|", 0 -> " "}}];
  ];

PrettyPrint[block1_, block2_] :=
  Module[ {m},
    m = Length[block1]/2;
    Print /@ StringJoin@@@Transpose[
     {Take[block1, m] /. {1 -> "_", 0 -> " "}, Drop[block1, m] /. {1 -> "|", 0 -> " "},
      Take[block2, m] /. {1 -> "_", 0 -> " "}, Drop[block2, m] /. {1 -> "|", 0 -> " "}}];
  ];

FitQ[block1_, block2_] := (* torus *)
  Module[ {m, out},
    (* determines whether two blocks are compatible, i.e., can be attached to each other
       so that all vertices at the cut line have even degree  *)
    m = Length[block1];
    If[ m != Length[block2], Return[False] ];
    m = m/2; out = True;
    Do[
      If[ OddQ[block1[[k]] + block2[[k]] + block1[[m + k]] + block1[[m + 1 + Mod[k, m]]]], out = False]
    , {k, 1, m}];
    Return[out]
  ];

(* transfer matrix *)
SAPMatrix[m_] := Table[If[FitQ[b1, b2], e^Count[b2, 1], 0], {b1, #}, {b2, #}]&[SAPAlphabet[m]]

(* the first n terms of the expansion for tori of height m *)
(* exponents of v encode the number of vertices in the underlying torus;
   exponents of e encode the number of edges of the polygon arrangement *)
SAPSeries[m_Integer, n_] := Module[ {M},
  M = SAPMatrix[m];
  1 + Sum[Expand[MatrixTrace[MatrixPower[M, k]]]v^(m k), {k, 1, n}]
]

(* generating function for tori of height m *)
(* exponents of v encode the number of vertices in the underlying torus;
   exponents of e encode the number of edges of the polygon arrangement *)
Clear[SAPRatfun]   
SAPRatfun[m_Integer] := SAPRatfun[m] = Module[ {M, rat},
  M = SAPMatrix[m];
  M = IdentityMatrix[Length[M]] - v^m M;
  rat = Together[Sum[Det[Delete[#, i]& /@ Delete[M, i]], {i, 1, Length[M]}]/Det[M]];
  Together[1 + rat - (rat /. {v -> 0, e -> 0})]
]

Rat2Poly[rat_, x_, n_] := Module[ {},
  (* given a rational function of the form poly(x)/(1-x)^d, return the polynomial
     that describes the coefficients of large index *)
  Expand[PolynomialQuotient[(Expand[rat /. x -> 1-x] /. (x^i_/;i>=1) -> 0 /. x -> 1/x) x, x^2, x] /. x^i_. -> Binomial[n + i, i] // FunctionExpand]
]

(* number of placements with k edges on a torus of height m with N vertices in total (viz width=N/m),
   as a polynomial in N *)
KZPoly[m_, k_, N_] := Module[ {n},
  Rat2Poly[Coefficient[Series[SAPRatfun[m], {e, 0, k}], e, k] /. v^n_ -> v^(n/m), v, N] /. N -> N/m
]