Needs["Combinatorica`"]


Needs["GraphUtilities`"]



Clear[EdgesWithVertices];
EdgesWithVertices[v_, e_] := 
 Module[{i, j, n = Length[e], m = Length[v], edges = {}},
  For[i = 1, i <= n, i++,
   For[j = 1, j <= m, j++,
     If[MemberQ[e[[i]], v[[j]]], AppendTo[edges, e[[i]]]; 
       Break[]];];]; edges]


Clear[SubGraph];
SubGraph[graph_, v_] := 
 Combinatorica`Graph[EdgesWithVertices[v, Edges[graph]], v]


Clear[BivariatePolynomialAGM1];

BivariatePolynomialAGM1[graph_, x_, y_] := 
 x^Length[Vertices[graph]] /; Length[Edges[graph]] == 0

BivariatePolynomialAGM1[graph_, x_, y_] := 
 Apply[Times, 
   Map[BivariatePolynomialAGM1[SubGraph[graph, #], x, y] &, 
    Combinatorica`ConnectedComponents[graph]]] /; 
  Not[Combinatorica`ConnectedGraphQ[graph]]

BivariatePolynomialAGM1[graph_, x_, y_] := Module[{e, e1, e2, g1, g2, g3},
   e = Edges[graph][[1]];
   g1 = DeleteEdge[graph, e];
   g2 = MakeSimple[Contract[graph, e]];
   e1 = Min[e];
   e2 = Max[e];
   g3 = DeleteVertices[graph, {e1, e2}];
   Return[
    Expand[BivariatePolynomialAGM1[g1, x, y] - 
      BivariatePolynomialAGM1[g2, x, 
       y] + (x - y) BivariatePolynomialAGM1[g3, x, y]]];] /; 
  Length[Edges[graph]] > 0

$RecursionLimit = Infinity


Clear[BivariatePolynomialAGM2];

BivariatePolynomialAGM2[graph_, x_, y_] :=
 Sum[Binomial[Length[Vertices[graph]], k]*(x - y)^k*
    FunctionExpand[
     FactorialPower[y, Length[Vertices[graph]] - k]], {k, 0, 
    Length[Vertices[graph]]}] /; 
  Length[Edges[graph]] == 
   Length[Vertices[graph]]*(Length[Vertices[graph]] - 1)/2

BivariatePolynomialAGM2[graph_, x_, y_] := 
 Apply[Times, 
   Map[BivariatePolynomialAGM2[SubGraph[graph, #], x, y] &, 
    Combinatorica`ConnectedComponents[graph]]] /; 
  Not[Combinatorica`ConnectedGraphQ[graph]]

BivariatePolynomialAGM2[graph_, x_, y_] := Module[{e, e1, e2, g1, g2, g3},
   e = Complement[Subsets[VertexList[graph], {2}], EdgeList[graph]][[
     1]];
   g1 = AddEdge[graph, e];
   g2 = MakeSimple[Contract[AddEdge[graph, e], e]];
   e1 = Min[e];
   e2 = Max[e];
   g3 = DeleteVertices[AddEdge[graph, e], {e1, e2}];
   Return[
    Expand[BivariatePolynomialAGM2[g1, x, y] + 
      BivariatePolynomialAGM2[g2, x, y] - (x - y)*
       BivariatePolynomialAGM2[g3, x, y]]];] /; 
  Length[Edges[graph]] < 
   Length[Vertices[graph]]*(Length[Vertices[graph]] - 1)/2


Clear[BivariatePolynomialGerling];

BivariatePolynomialGerling[r_,1,x_,y_]:=x^r
BivariatePolynomialGerling[r_,k_,x_,y_]:=Sum[Binomial[r,i]*(x-y)^i*Sum[StirlingS2[r-i,j]*FunctionExpand[FactorialPower[y,j]]*BivariatePolynomialGerling[r,k-1,x-j,y-j],{j,0,r-i}],{i,0,r}]