(*Algorithm for finding the generating function F(A, d) described in Theorem \
6.1 of the paper Matjaz Konvalinka : A generalization of Foata's fundamental \
transformation and its applications to the right - quantum algebra and \
enumeration*)

(*March 4, 2007*)

BeginPackage["genmacmahon`"];

F::usage = "F[A_,d_,x_] gives the generating function G(A,p,r)x^p, where A is a square matrix, d is a vector of the same size with entries summing up to 0, p and r non-zero integer vectors with p = r + d and G(A,p,r) the coefficient [x^r](A.x)^p.";

Begin["Private`"];

Del[A_, i_, j_,x_] := Det[(IdentityMatrix[
    Length[A]] - 
          DiagonalMatrix[Table[x[[k]], {k, Length[A]}]] . A)[[
            Complement[Range[Length[A]], i], 
                Complement[Range[Length[
                  A]], j]]]]; F[A_,d_,x_]:=Which[ !MatrixQ[A] || Length[A] != Length[First[A]], 
  Print[A, " is not a square matrix."],  !VectorQ[d, IntegerQ] || Plus @@ d != 0 || 
   Length[d] != Length[A], Print[d, " is not an integer vector with the sum of \
entries equal to 0, or is not of the same size as A."],!VectorQ[x] || Length[x] != Length[A],Print[x," is
not a vector of the same size as A and d."], True, 
  Module[{q = Flatten[(Table[#1, {-d[[#1]]}] & ) /@ Select[Range[Length[A]], 
        d[[#1]] < 0 & ]], s = Flatten[(Table[#1, {d[[#1]]}] & ) /@ 
       Range[Length[A]]]}, 
   Plus @@ (((-1)^(Plus @@ s + Plus @@ q)*Product[
         (-1)^Length[Intersection[Union[Range[s[[i]] + 1, #1[[i]] - 1], 
              Range[#1[[i]] + 1, s[[i]] - 1]], Union[Take[#1, i - 1]]]]*
          (Del[A, Union[Take[#1, i]], Complement[Union[Take[#1, i - 1], {s[[i]]}], 
             {#1[[i]]}],x]/Del[A, Union[Take[#1, i]], Union[Take[#1, i]],x]), 
         {i, Length[#1]}])/Del[A, {}, {},x] & ) /@ Permutations[q]]];
         
End[];

EndPackage[];