(*  LamanGraphs -- construction of Laman graphs and their embedding numbers
 *  Copyright (C) 2016  Christoph Koutschan
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version, see http://www.gnu.org/licenses/
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *)

If[$Notebooks,
  CellPrint[Cell[TextData[{#1, StyleBox[#2, Small]}], "Text", FontColor -> Black, CellFrame -> 0.5,
    Background -> RGBColor[0.796887, 0.789075, 0.871107]]],
  Print[#1 <> #2]]&[
  "LamanGraphs   \[LongDash]   Copyright \[Copyright] 2016   Christoph Koutschan\n",
  "This program comes with absolutely no warranty; it is free software,\n" <>
  "and you are welcome to redistribute it and/or modify it under the\n" <>
  "terms of the GNU General Public License (http://www.gnu.org/licenses/)."];


(*****************************************************************************
 *           Procedures for generating and counting Laman graphs             *
 *****************************************************************************)

(* We represent a graph by the integer obtained by flattening the upper triangle
 * of its adjacency matrix and interpreting this binary sequence as an integer.
 * Conversely, G2Mat retrieves the adjacency matrix from this integer representation.
 *)
Mat2G[mat_] := FromDigits[Flatten[MapIndexed[Drop[#1, #2[[1]]]&, mat]], 2];
G2Mat[g_Integer, n_Integer] := (# + Transpose[#])&[PadLeft[
  Table[Take[#, {(2n-i)*(i-1)/2+1, (2n-i-1)*i/2}], {i, n}]&[PadLeft[IntegerDigits[g, 2], n*(n-1)/2]], {n, n}]];
G2Mat[g_Integer] := G2Mat[g, Floor[(3 + Sqrt[8*Floor[Log[2, g]] + 1]) / 2]];

(* Transform a graph to a list of edges and vice versa. *)
G2Edges[g_Integer] := Select[Position[G2Mat[g], 1], Less @@ # &];
Mat2Edges[mat_] := Select[Position[mat, 1], Less @@ # &];
Edges2Mat[ed_List] :=
Module[{n = Max[ed], mat},
  mat = Table[0, {n}, {n}];
  (mat[[##]] = 1) & @@@ ed;
  Return[mat + Transpose[mat]];
];
Edges2G[ed_List] := Mat2G[Edges2Mat[ed]];


(* Input g is an undirected graph represented by its symmetric n*n adjacency matrix.
 * Among all graphs that differ from g only by relabelling the vertices, compute a
 * unique representative; we chose this representative to be the lexicographically
 * largest among all equivalent graphs with decreasing vertex neighborhoods.
 * "Lexicographically" in the sense that we interpret the adjacency matrix as a
 * {0,1}-word of length n^2 (not by flattening but by taking blocks of increasing size).
 *)
GraphNormalForm[g_Integer] := Mat2G[GraphNormalForm[G2Mat[g]]];
GraphNormalForm[gr1_List] :=
Module[{gr = gr1, rn = Range[Length[gr1]], nbs, perm, b, p1},

  (* Apply some permutation to the graph s.t. the neighborhoods of its vertices are
     decreasing (by length and lexicographically). By a "neighborhood" of a vertex
     we mean the sorted list of valencies of the neighboring vertices. *)
  nbs = Total /@ gr; nbs = Sort[Pick[nbs, #, 1]]& /@ gr;
  {nbs, perm} = Transpose[Reverse[Sort[Transpose[{nbs, rn}]]]];
  gr = gr[[perm,perm]];

  (* Among all graphs whose vertices are ordered as above, pick the one with
     "lexicographically largest" matrix. *)
  b = Map[Last, SplitBy[Transpose[{nbs, rn}], First], {2}];
  perm = {{}};
  Do[
    perm = Join @@ Outer[Join, perm, Permutations[b[[k]]], 1];
    (* There seems to be a memory problem with Outer. The below two lines *)
    (* achieve the same result with less memory consumption, but are slightly slower. *)
    (*p1 = Permutations[b[[k]]];
    perm = Join @@ (Function[p, Join[p, #]& /@ p1] /@ perm);*)
    p1 = FromDigits[Flatten[gr[[#,#]]], 2]& /@ perm;
    perm = Pick[perm, p1, Max[p1]];
  , {k, Length[b]}];
  Return[gr[[#,#]]& @ perm[[1]]];
];


(* Given g, an encoded graph with n-1 vertices, compute all graphs (with n vertices)
 * that can be obtained from g by a single Henneberg move (type I or type II).
 *)
Henneberg[g_Integer] := Henneberg[Floor[(3 + Sqrt[8*Floor[Log[2, g]] + 1]) / 2], g];
Henneberg[n_Integer, g_Integer] :=
Module[{gr, ed},
  gr = PadRight[G2Mat[g, n - 1], {n, n}];
  ed = G2Edges[g];
  Union[Flatten[{

    (* Perform all possible type I Henneberg moves. *)
    Table[Mat2G[GraphNormalForm[ReplacePart[gr, {{v1, n} -> 1, {n, v1} -> 1, {v2, n} -> 1, {n, v2} -> 1}]]],
      {v1, n - 1}, {v2, v1 + 1, n - 1}],

    (* Perform all possible type II Henneberg moves. *)
    Table[
      {v1, v2} = ed[[k]];
      Mat2G[GraphNormalForm[ReplacePart[gr, {{v1, v2} -> 0, {v2, v1} -> 0, {v1, n} -> 1, {n, v1} -> 1,
        {v2, n} -> 1, {n, v2} -> 1, {#, n} -> 1, {n, #} -> 1}]]]& /@ Complement[Range[n - 1], {v1, v2}],
      {k, Length[ed]}]
  }]]];


(* Similar as before, but only with Henneberg type I moves. *)
Henneberg1[g_Integer] := Henneberg1[Floor[(3 + Sqrt[8*Floor[Log[2, g]] + 1]) / 2], g];
Henneberg1[n_Integer, g_] :=
Module[{gr},
  gr = PadRight[G2Mat[g, n - 1], {n, n}];
  Union[Flatten[
    (* Perform all possible type I Henneberg moves. *)
    Table[Mat2G[GraphNormalForm[ReplacePart[gr, {{v1, n} -> 1, {n, v1} -> 1, {v2, n} -> 1, {n, v2} -> 1}]]],
      {v1, n - 1}, {v2, v1 + 1, n - 1}]
  ]]];


(* Compute a list of Laman graphs with n vertices. *)
LamanGraphs[3] = {7};
LamanGraphs[n_ /; n > 3] := LamanGraphs[n] = Union @@ (Henneberg[n, #]& /@ LamanGraphs[n - 1]);

(* Compute a list of H1-Laman graphs (i.e., obtained by only H1 moves) with n vertices. *)
H1LamanGraphs[3] = {7};
H1LamanGraphs[n_ /; n > 3] := H1LamanGraphs[n] = Union @@ (Henneberg1[n, #]& /@ H1LamanGraphs[n - 1]);


(* Test whether g is a Laman graph by performing inverse Henneberg moves. *)
LamanGraphQ[g_Integer] := LamanGraphQ[G2Mat[g]];
LamanGraphQ[gr1_List] :=
Module[{gr = gr1, p, s, vs, lq},
  If[gr === {{0, 1}, {1, 0}}, Return[True]];
  p = Position[Total /@ gr, 2];
  If[p =!= {}, (* Henneberg type I move *) p = p[[1, 1]]; Return[LamanGraphQ[Delete[Transpose[Delete[gr, p]], p]]]];
  p = Position[Total /@ gr, 3];
  If[p =!= {}, (* Henneberg type II move *)
    p = p[[1, 1]];
    vs = Flatten[Position[gr[[p]], 1]];
    s = Select[Subsets[vs, {2}], gr[[#[[1]], #[[2]]]] === 0 &];
    If[s === {}, Return[False]];
    Do[
      gr = gr1;
      (gr[[##]] = 1)& @@@ Permutations[s[[i]]];
      If[lq = LamanGraphQ[Delete[Transpose[Delete[gr, p]], p]], Break[]];
    , {i, Length[s]}];
    Return[lq];
  ];
  Return[False];
];


(*****************************************************************************
 *            Procedures for computing the number of embeddings              *
 *****************************************************************************)

MyRI := 10^6 + RandomInteger[10^6];

(* Number of embeddings of a Laman graph *)
NumberOfEmbeddings[g_Integer] :=
Module[{n, gr, p, s, t, perm, opt, ed, ed1, eqs = {}, vars, fix, emb = 1, gb, eq1},
  gr = G2Mat[g];
  n = Length[gr];

  (* Remove vertices of valency 2. *)
  While[n > 3 && (p = Position[Total /@ gr, 2]) =!= {},
    gr = Delete[Transpose[Delete[gr, p]], p];
    p = Length[p];
    emb *= 2^p;
    n -= p;
  ];

  (* Test whether the graph is separable. *)
  (* The smallest graph which (after the above step) is separable has 7 vertices. *)
  Do[
    s = t = {First[Complement[{1, 2, 3}, {i, j}]]};
    While[(t = Complement[Flatten[Position[Total[gr[[t]]], a_?Positive]], s, {i, j}]) =!= {}, s = Join[s, t]];
    If[Length[s] < n - 2,
      s = Join[{i, j}, #]& /@ {s, Complement[Range[n], s, {i, j}]};
      (p = gr[[#,#]];
       t = 2 * Length[#] - Total[Flatten[p]] / 2;
       Which[
         t > 4 || t < 3, Throw["Problem with separation of graphs."],
         t === 4, p = ReplacePart[p, {{1, 2} -> 1, {2, 1} -> 1}];
       ];
       emb *= NumberOfEmbeddings[Mat2G[p]])& /@ s;
      gr = {}; Break[];
    ];
  , {i, If[n < 7, 0, n-1]}, {j, i+1, n}];
  If[gr === {}, Return[emb]];

  (* Identify an edge that is part of the largest triangulation subgraph. *)
  (* Relabel the vertices s.t. this edge becomes (1,2). *)
  ed = Mat2Edges[gr];
  opt = -1;
  Do[
    fix = ed[[i]];
    While[(p = Complement[Flatten[Position[Plus @@ gr[[fix]], 2]], fix]) =!= {}, fix = Join[fix, p]];
    If[Length[fix] > opt, opt = Length[fix]; ed1 = ed[[i]]];
  , {i, Length[ed]}];
  perm = Range[n];
  MapThread[(perm[[{#1, #2}]] = perm[[{#2, #1}]])&, {{1, 2}, ed1}];
  gr = gr[[perm,perm]];

  (* Generate the equations *)
  vars = {};
  Do[
    If[gr[[i, j]] === 1, AppendTo[vars, z[i, j]];
(*      AppendTo[eqs, (x[i] - x[j])^2 + (y[i] - y[j])^2 - z[i, j]];*)
      (* in the following x[i], y[i] stand for v[i], w[i], resp. *)
      AppendTo[eqs, (x[i] - x[j]) * (y[i] - y[j]) - z[i, j]]];
  , {i, 1, n - 1}, {j, i + 1, n}];
  eqs = eqs /. {x[1] -> 0, y[1] -> 0, y[2] -> MyRI};

  Return[emb * NumberOfSolutions[Rest[eqs] /. x[2] -> MyRI]];
];


(* Compute the number of solutions (complex, counted with multiplicity) of the polynomial system eqs. *)
NumberOfSolutions[eqs_List] :=
Module[{eq1, p, sol, vars, gb},
  If[eqs === {}, Return[1]];

  (* Identify vertices that are part of a triangle whose other two vertices have already been fixed. *)  
  eq1 = Select[eqs, Length[Variables[#]] === 3 &];
  p = Cases[Tally[Cases[eq1, x[a_] :> a, Infinity]], {a_, 2} :> a];

  If[p =!= {},
    eq1 = eqs /. Thread[{x[p[[1]]], y[p[[1]]]} -> Table[MyRI, {2}]];
    Return[2 * NumberOfSolutions[Select[eq1, Not[FreeQ[#, x | y]] &]]];
  ,
    (* If no such triangle exists, we determine the solutions by a GB computation. *)
    eq1 = eqs /. z -> (MyRI &);
    vars = First /@ SortBy[Tally[Cases[eq1, x[_], Infinity]], Last];
    vars = Riffle[vars, y @@@ vars];
    (* Lexicographic *)
    (*gb = GroebnerBasis[eq1, vars, Modulus -> NextPrime[2^31, -2]];
    If[Length[gb] =!= Length[vars] || Not[And @@ MapIndexed[FreeQ[#1, Alternatives @@ Drop[vars, -#2[[1]]]]&, gb]],
      Throw["Groebner basis is not of required form."]];
    Return[Times @@ MapThread[Exponent, {gb, Reverse[vars]}]];*)
    (* DegreeReverseLexicographic *)
    gb = GroebnerBasis[eq1, vars, MonomialOrder -> DegreeReverseLexicographic, Modulus -> NextPrime[2^31, -2]];
    sol = UnderTheStaircase[Exponent[LMDRL[#, vars], vars] & /@ gb];
    If[Not[MatchQ[sol, {(_List) ..}]], Throw["Problem with gb."]];
    Return[Length[sol]];
  ];
];


(* Given the leading monomials of some Groebner basis (represented as exponent vectors),
 * compute the set of irreducible monomials ("monomials under the staircase").
 * (This procedure is copied from the HolonomicFunctions package.)
 *)
UnderTheStaircase[leadingMonomials_?(MatrixQ[#, IntegerQ]&)] :=
Module[{i, dim, max, vec, pos, vsbasis},

  dim = Length[First[leadingMonomials]];

  (* Check whether the ideal is zero-dimensional,
     i.e., for each variable there is a leading monomial pure in this variable *)
  If[Or @@ Table[Cases[leadingMonomials, Table[If[j === i, _, 0], {j, dim}]] === {}, {i, dim}],
    Return[Infinity]];

  max = Max[Cases[leadingMonomials, Prepend[Table[0, {dim-1}], i_] :> i]];
  vec = Table[0, {dim}];
  pos = dim;

  vsbasis = {};
  While[First[vec] < max && pos > 0,
    If[And @@ ((Or @@ Thread[# > vec]) & /@ leadingMonomials),
      AppendTo[vsbasis, vec];
      pos = dim;
      vec[[pos]]++;
    ,
      vec[[pos]] = 0;
      pos--;
      If[pos > 0, vec[[pos]]++];
    ];
  ];
  Return[vsbasis];
];


(* Leading monomial w.r.t. DegreeReverseLexicographic ordering. *)
LMDRL[poly_, vars_List] :=
Module[{terms, t},
  terms = If[Head[#] === Plus, List @@ #, {#}]&[Expand[poly]];
  t = Total[Exponent[#, vars]]& /@ terms;
  terms = Pick[terms, t, Max[t]];
  (t = Exponent[terms, #]; terms = Pick[terms, t, Min[t]])& /@ Reverse[vars];
  Return[First[terms / (terms /. Thread[vars -> 1])]];
];


(*** Subroutines for Josef's combinatorial algorithm ***)

(* Components of a graph, given as lists of vertices. *)
Components[e_List] := Components[Union[Flatten[e]], e];
Components[v_List, e_List] :=
Module[{comps = {}, c},
  While[(c = Complement[v, Union[Flatten[comps]]]) =!= {},
    AppendTo[comps, FixedPoint[Union[#, Cases[e, ({a_, Alternatives @@ #} | {Alternatives @@ #, a_}) :> a]]&, Take[c, 1]]];
  ];
  Return[comps];
];

(* Quotient a graph (given by a list of edges) by a subset of the edges (given by their positions in the list). *)
QuotientEdges[e_List, p_List] :=
  With[{e1 = Delete[e, Transpose[{p}]] /. Flatten[MapIndexed[(#1 -> #2[[1]])&, Components[Union[Flatten[e]], e[[p]]], {2}]]},
    If[Cases[e1, {a_, a_}] =!= {}, {}, Sort /@ e1]];

(* List of triangles in a graph (given by its edge list).
 * The triangles are given by the positions of their edges in the list. *)
Triangles[e1_List] :=
Module[{e = e1, k = Length[e1], tr = {}, vs, t},
  While[k > 2,
    vs = e[[-1]];
    e = Most[e];
    t = Tuples[Position[e, #, {2}]& /@ vs];
    tr = Join[tr, Append[First /@ #, k]& /@ Pick[t, SameQ @@@ Apply[e[[#1, 3 - #2]]&, t, {2}]]];
    k--;
  ];
  Return[tr];
];

(* List of quadrilaterals in a graph (given by its edge list).
 * The quadrilaterals are given by the positions of their edges in the list. *)
Quadrilaterals[e1_List] :=
Module[{e = e1, k = Length[e1], qu = {}, v1, v2, p, t},
  While[k > 3,
    {v1, v2} = e[[-1]];
    e = Most[e];
    qu = Join[qu, Join @@ (Function[{i, j},
      t = Tuples[Position[e, #, {2}]& /@ {e[[i, 3 - j]], v1}];
      Join[{k, i}, First /@ #]& /@ Pick[t, SameQ @@@ Apply[e[[#1, 3 - #2]]&, t, {2}]]
    ] @@@ Position[e, v2, {2}])];
    k--;
  ];
  Return[Select[qu, Length[Union[#]] === 4 &]];
];


(* Laman number (= number of embeddings) of a graph, using Josef's combinatorial algorithm. *)
(* LamanNumber[g_Integer] := LamanNumber[G2Edges[g]]; *)
SubsetsEx[t_, i_] := Total /@ (2^(Join[Subsets[t, i - 1], Subsets[t, {i + 1, Length[t]}]] - 1));
LamanNumber[e_List] := LamanNumber[e, e];
LamanNumber[e1a_List, e2a_List] :=
Module[{e1 = e1a, e2 = e2a, k, e, sum, l, v1, v2, t1, t2, q1, q2, pick, subsets},
  If[(k = Length[e1]) === 1, Return[1]];
  {v1, v2} = Union[Flatten[#]]& /@ {e1, e2};
  {t1, t2} = Length /@ {v1, v2};
  If[k <= 4,
    (* If there are few edges, normalize the input so that we have fewer cases to remember. *)
    If[v1 =!= Range[t1] || v2 =!= Range[t2] || e1 =!= Sort[e1],
      q1 = e1 /. Thread[v1 -> Range[t1]];
      q2 = e2 /. Thread[v2 -> Range[t2]];
      {q1, q2} = Transpose[SortBy[Transpose[{q1, q2}], First]];
      Return[LamanNumber[q1, q2] = LamanNumber[q1, q2]];
    ];
  ];
  If[t1 - Length[Components[e1]] + t2 - Length[Components[e2]] =!= k + 1, Return[0]];
  If[k >= 5,
    (* Choose a "good" edge e. *)
    {t1, t2} = Triangles /@ {e1, e2};
    e = SortBy[Tally[Flatten[{t1, t2, Range[k]}]], Last][[1,1]];
    (* Relabel s.t. the chosen edge becomes the last (i.e., k-th) one. *)
    With[{p = ReplacePart[Range[k], {e -> k, k -> e}]}, e1 = e1[[p]]; e2 = e2[[p]]];
    {t1, t2} = {t1, t2} /. {e -> k, k -> e};
    (* Subsets encoded as integers. *)
    subsets = Range[1, 2^(k-1) - 2];
    pick = Function[{t, p}, subsets = Pick[subsets, BitAnd[subsets, Total[2^(t - 1)]], Alternatives @@ p];];
    (* Exclude subsets for which M resp. N involves some but not all members of a multiple edge. *)
    {q1, q2} = Function[e3, Flatten[Position[e3, #]]& /@ Cases[Tally[e3], {a_, b_ /; b > 1} :> a]] /@ {e1, e2};
    pick[#, {0, Total[2^(# - 1)]}]& /@ q1;
    pick[#, If[MemberQ[#, k], {Total[2^(DeleteCases[#, k] - 1)]}, {0, Total[2^(# - 1)]}]]& /@ q2;
    (* Exclude subsets for which M resp. N involves precisely 2 edges of a triangle in e2 resp. e1. *)
    pick[#, SubsetsEx[#, 1]]& /@ t1;
    pick[#, SubsetsEx[#, If[MemberQ[#, k], 1, 2]]]& /@ t2;
  ];
  If[k >= 9 && Length[subsets] >= 30,
    (* Exclude subsets for which M resp. N involves precisely 3 edges of a quadrilateral in e2 resp. e1. *)
    {q1, q2} = Quadrilaterals /@ {e1, e2};
    pick[#, SubsetsEx[#, 1]]& /@ q1;
    pick[#, SubsetsEx[#, If[MemberQ[#, k], 2, 3]]]& /@ q2;
  ];
  If[k < 5,
    subsets = Subsets[Range[k - 1], {1, k - 2}];
  ,
    subsets = Flatten[Position[Reverse[IntegerDigits[#, 2]], 1]]& /@ subsets;
  ];
  sum = With[{M = Append[#, k], N = Complement[Range[k], #]},
      If[(Length[M] < k / 2 && (q1 = QuotientEdges[e1, N]) =!= {} && (q2 = QuotientEdges[e2, M]) =!= {}) ||
         (Length[M] >=  k / 2 && (q2 = QuotientEdges[e2, M]) =!= {} && (q1 = QuotientEdges[e1, N]) =!= {}),
        If[(l = LamanNumber[q1, Delete[e2, Transpose[{N}]]]) === 0, 0,
          l * LamanNumber[Delete[e1, Transpose[{M}]], q2]], 0]]& /@ subsets;
  sum = Total[sum];
  l = If[(q2 = QuotientEdges[e2, {k}]) =!= {}, LamanNumber[Most[e1], q2], 0];
  sum += If[e1 === e2, 2 * l, l + If[(q1 = QuotientEdges[e1, {k}]) =!= {}, LamanNumber[q1, Most[e2]], 0]];
  Return[sum];
];


(* Laman number (= number of embeddings) of a graph, using Josef's combinatorial algorithm
 * together with some preprocessing (code copied from NumberOfEmbeddings, see above). *)
LamanNumber[g_Integer] :=
Module[{n, gr, p, s, t, emb = 1},
  gr = G2Mat[g];
  n = Length[gr];

  (* Remove vertices of valency 2. *)
  While[n > 3 && (p = Position[Total /@ gr, 2]) =!= {},
    gr = Delete[Transpose[Delete[gr, p]], p];
    p = Length[p];
    emb *= 2^p;
    n -= p;
  ];

  (* Test whether the graph is separable. *)
  (* The smallest graph which (after the above step) is separable has 7 vertices. *)
  Do[
    s = t = {First[Complement[{1, 2, 3}, {i, j}]]};
    While[(t = Complement[Flatten[Position[Total[gr[[t]]], a_?Positive]], s, {i, j}]) =!= {}, s = Join[s, t]];
    If[Length[s] < n - 2,
      s = Join[{i, j}, #]& /@ {s, Complement[Range[n], s, {i, j}]};
      (p = gr[[#,#]];
       t = 2 * Length[#] - Total[Flatten[p]] / 2;
       Which[
         t > 4 || t < 3, Throw["Problem with separation of graphs."],
         t === 4, p = ReplacePart[p, {{1, 2} -> 1, {2, 1} -> 1}];
       ];
       emb *= LamanNumber[Mat2G[p]])& /@ s;
      gr = {}; Break[];
    ];
  , {i, If[n < 7, 0, n-1]}, {j, i+1, n}];
  If[gr === {}, Return[emb]];

  Return[emb * LamanNumber[Mat2Edges[gr]]];
];


(* This procedure can be used to draw the computation tree of the combinatorial algorithm. *)
Clear[QE, LN];
QE[e_List, p_List] :=
  Sort /@ (Delete[e, Transpose[{p}]] /. Flatten[MapIndexed[(#1 -> #2[[1]])&, Components[Union[Flatten[e]], e[[p]]], {2}]]);
LN[e_List] := LN[e, e];
LN[e1_List, e2_List] :=
Module[{k, e = 1, sum, l, v1, v2, d1, d2, t1, t2, q1, q2, subsets, ret = Return},
  If[Cases[{e1, e2}, {a_, a_}, {2}] =!= {}, ret[L[e1, e2][0]]; Return[0]];
  If[(k = Length[e1]) === 1, ret[L[e1, e2][1]]; Return[1]];
  {v1, v2} = Union[Flatten[#]]& /@ {e1, e2};
  {t1, t2} = Length /@ {v1, v2};
  If[t1 - Length[Components[e1]] + t2 - Length[Components[e2]] =!= k + 1, ret[L[e1, e2][0]]; Return[0]];
  subsets = Subsets[Rest[Range[k]], {1, k - 2}];
  sum = With[{M = Append[#, e], N = Complement[Range[k], #]},
      q1 = QE[e1, N];
      q2 = QE[e2, M];
      d1 = Delete[e1, Transpose[{M}]];
      d2 = Delete[e2, Transpose[{N}]];
    times[LN[q1, d2], LN[d1, q2]]]& /@ subsets;
  q1 = QE[e1, {e}];
  q2 = QE[e2, {e}];
  d1 = Delete[e1, e];
  d2 = Delete[e2, e];
  sum = plus @@ Join[{LN[d1, q2], LN[q1, d2]}, sum];
  ret[L[e1, e2][sum]];
  Return[sum /. L[__] -> Identity /. plus -> Plus /. times -> Times];
];
