##########################################################################
## JointConductance.txt Save this file as JointConductance.txt          ##
## to use it, stay in the same directory, get into Maple                ##
## (by typing: maple <Enter> ) and then type:                           ##
## read JointConductance.txt <Enter>                                    ##
## Then follow the instructions given there                             ##
##                                                                      ##
## Written by Yukun Yao and Doron Zeilberger, Rutgers University ,      ## 
## yao@math.rutgers.edu and DoronZeil@gmail.com                         ##
##########################################################################

#Version of Oct. 15, 2018
with(combinat):
with(linalg):
with(plots):
with(LinearAlgebra):
with(numtheory):

print(`First Written: Aug 2018`):
print(`Version of Sept. 12, 2018:`):
print():
print(`This is JointConductance.txt, one of the  Maple packages`):
print(`accompanying the  article: `):
print(`Untying The Gordian Knot via Experimental Mathematics`):
print(` by Yukun Yao and Doron Zeilberger `):

print():
print(`The most current version is available on WWW at:`):
print(` https://sites.math.rutgers.edu/~yao/JointConductance.txt .`):
print(`Please report all bugs to: yao at math dot rutgers dot edu .`):
print():
print(`For an overview of all functions type “Help()”`):
print(`For a list of the MAIN functions,`):
print(`type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the supporting functions type: ezraS();`):
print(`For a list of the Checking procedures type ezraC().`):
print():


Help:=proc(): print(`All procedures are: L(a,n), Lg(a,n,G), JC(a,n) , JCpq(a,n,p,q) , JCg(a,n,G), JCgND(n,G)`):
print(`JCgpq(a,n,G,p,q) , JCnum(n,G), JCnumND, JCnumpq(n,G,p,q) , JCgND(a,n,g) , BR(i,n) , Neis(i,j,n),  UC(n), Rn(n), RnF(n) , BRmn(i,m,n) , `):
print(`Neismn(i,j,m,n), UCmn(m,n), Rmn(m,n) , RmnND(m,n), Rmnpq(m,n,p,q) ,  BRkn(i,k,n) , Neiskn(i,j,k,n) , UCkn(k,n), Rkn(k,n) , `):
print(`BRL(i,L) , NeisL(i,j,L) , UCL(L) , RL(L) , ST2n(n) , SF2n(n) , R2n(n) , STkn(k,n) , ReArrow(arrow) , OneRecurrence(L) , OneList(L) , OriginalCase(L) , AllCase(L) , AllCaseL(L) , FinalRecurrence(L) , Coefficient(S,L) , FinalVector(S,L) , FinalMatrix(L) , FinalGF(L,x) , CP(L) ,`):
print(` NeighborV(G,S) , IsConnected(G) , DegreeG(G,n) , AllComponent(G) , NumberComponent(G) , GST(G) , IsAcyclic(G) , GST(G) , gst(G) , LowestDegree(G,S) , GSF(G,S) , OneStepG(G,s,v) , randgraph(N,K) , EdgeToGraph(E)`): 
print(` GuessRat(L,degT,degB) , GuessRec1(L,d) , GuessRec(L) , GuessSymRec1(L,d) , GuessSymRec(L) , CtoR(s,t) , SeqFromRec(s,N) , GridMN(M,N) , Mat11(G,n,a) , Mat1(G,n,v0,a) , IndToEdge(M1,a) , MonoToGraph(M,a) , PolyToSG(p,a) , SpF(G,n,v0) , Kn(n) , Mat11N(G,n) , Mat1N(G,n,v0) , SpFn(G,n,v0) `):
print(` GFGridKN(k,t) , DenomSFKN(k,t) , GridH(m,n,t) , GridV(m,n,t) , GridHV(m,n,x,y) , GridHFactorial(i,k,n) , GridVFactorial(i,k,n) , ListSTPreCal(k) , GFGridPreCal(k,t) , DenomSFPreCal(k,t) , DiagMatrix(a,b,c,p,q,n) , ExpandMatrix(M) , IsDiagonal(a,b,c,p,q,M) , LSUnion(A,B) , MatrixRecurrence(a,b,c,p,q,) `):
end:

ezra:=proc()

if args=NULL then
print(`The main procedures are: L , Lg , JC , JCpq , JCg , JCgND, JCgpq , JCnum ,JCnumND,  JCnumpq , Rn , RnF , Rmn ,RmnND, Rmnpq , Rkn , RL , R2n `):
print(`OneRecurrence , OneList , OriginalCase , AllCase , FinalRecurrence , FinalMatrix , FinalGF , CP , GST , GSF , OneStepG , randgraph`):
print(`GuessRec1 , GuessRec , GuessSymRec1 , GuessSymRec , CtoR , GridMN , SpF , SpFn , GFGridKN , DenomSFKN , GridH , GridV , GridHV , GridHFactorial , GridVFactorial , ListSTPreCal , GFGridPreCal , DenomSFPreCal , GuessRat , DiagMatrix , ExpandMatrix , MatrixRecurrence `):


elif nops([args])=1 and op(1,[args])=GuessRec then
print(`GuessRec(L): inputs a sequence L and tries to guess`):
print(`a recurrence operator with constant coefficients `):
print(`satisfying it. It returns the initial  values and the operator`):
print(` as a list. For example try:`):
print(`GuessRec([1,1,1,1,1,1]);`):

elif nops([args])=1 and op(1,[args])=L then
print(`L(a,n): The Laplacian of the complete graph.`):

elif nops([args])=1 and op(1,[args])=Lg then
 print(`Lg(a,n,G): The Laplacian of the graph G.`):

elif nops([args])=1 and op(1,[args])=JC then
print(`JC(a,n): the joint conductance between nodes 1 and n.`):

elif nops([args])=1 and op(1,[args])=JCpq then
print(`JCpq(a,n,p,q): the joint conductance between nodes p and q.`):

elif nops([args])=1 and op(1,[args])=JCg then
print(`JCg(a,n,G): the joint conductance between nodes 1 and n in the graph G. Try:`):
print(`JCg(a,3,{{1,2},{2,3}});`):

elif nops([args])=1 and op(1,[args])=JCgND then
print(`JCgND(a,n,G): Like JCg(a,n,G) but witk numerator and denominator separate. Try:`):
print(`JCgND(a,3,{{1,2},{2,3}});`):

elif nops([args])=1 and op(1,[args])=JCgpq then
print(`JCgpq(a,n,G,p,q): the joint conductance between nodes p and q in the graph G. Try:`):
print(`JCg(a,3,{{1,2},{2,3}},3,2);`):

elif nops([args])=1 and op(1,[args])=JCnum then
print(`JCnum(n,G): the joint conductance between nodes 1 and n in the graph G where all resistances are with 1 Ohm. Try:`):
print(`JCnum(3,{{1,2},{2,3}});`):

elif nops([args])=1 and op(1,[args])=JCnumND then
print(`JCnumND(n,G): the numerator and denominator of the`):
print(`joint conductance between nodes 1 and n in the graph G where all resistances are with 1 Ohm. Try:`):
print(`JCnumND(3,{{1,2},{2,3}});`):

elif nops([args])=1 and op(1,[args])=JCnumpq then
print(`JCnumpq(n,G,p,q): the joint conductance between nodes p and q in the graph G where all resistances are with 1 Ohm. Try:`):
print(`JCnumpq(3,{{1,2},{2,3}},3,2);`):

elif nops([args])=1 and op(1,[args])=Rn then
print(`Rn(n): the formula for the joint resistance between opposite vertices via the general formula.`):

elif nops([args])=1 and op(1,[args])=RnF then
print(`RnF(n): the joint resistance between opposite vertices of the n-dim unit-cube.`):

elif nops([args])=1 and op(1,[args])=Rmn then
print(`Rmn(m,n): the joint resistance between opposite vertices of a m by n grid via the general formula.`):

elif nops([args])=1 and op(1,[args])=RmnND then
print(`RmnND(m,n): the numerator and denominator for the joint resistance between opposite vertices of a m by n grid via the general formula.`):

elif nops([args])=1 and op(1,[args])=Rmnpq then
print(`Rmnpq(m,n,p,q): the formula for joint resistance between node p and q in a m by n grid via the general formula.`):

elif nops([args])=1 and op(1,[args])=Rkn then
print(`Rkn(k,n): the formula for the joint resistance between opposite vertices of [0,k-1]^n via the general formula.`):

elif nops([args])=1 and op(1,[args])=RL then
print(`RL(L): the formula for the joint resistance between opposite vertices of L-shape hyper-rectangle via the general formula.`):

elif nops([args])=1 and op(1,[args])=R2n then
print(`R2n(n): the resistance of 2 by n grid with 1 Ohm each edge.`):

elif nops([args])=1 and op(1,[args])=OneRecurrence then
print(`OneRecurrence(L): input a list, which represents a data structure of k by n grid (k fixed, e.g. [2,0,2,{}] and outputs a set of lists where each list is a data structure following by a number, e.g. {[[2,0,1,{}],2],[[2,0,1,{[1,2],[2,1]}],1]}. It means the number of spanning trees in L is the sum of the corresponding number times the number of spanning trees in those data structures represented by the list. Here 0 means the length is n, -1 means the length is n-1.`):

elif nops([args])=1 and op(1,[args])=OneList then
print(`OneList(L): List only for the output, similar to the above procedure.  `):

elif nops([args])=1 and op(1,[args])=OriginalCase then
print(`OriginalCase(L): input a list L=[k,n,p,arrow] as a data structure, depending on n, output its similar original case.`):

elif nops([args])=1 and op(1,[args])=AllCase then
print(`AllCase(L): input a list L as a data structure, output a set of all cases when we apply OneRecurrence to L and its output(whose n=0 but not negative).`):

elif nops([args])=1 and op(1,[args])=FinalRecurrence then
print(`FinalRecurrence(L): input a list L as a data structure, output a set of lists where each sublist is [data structure of -1, number].`):

elif nops([args])=1 and op(1,[args])=FinalMatrix then
print(`FinalMatrix(L): input a list L as a data structure, output a transition matrix. The order of the entry depends on the position in the list AllCaseL(L) where L will always be the first one. `):

elif nops([args])=1 and op(1,[args])=FinalGF then
print(`FinalGF(L,x): input a list L as a data structure, output its generating function such that the coefficient of x^n for the generating function means the number of spanning trees/forests in this data structure with size k by n.  `):

elif nops([args])=1 and op(1,[args])=CP then
print(`CP(L): input a list L as a data structure, output its characteristic polynomial for its transition matrix.  `):

elif nops([args])=1 and op(1,[args])=GST then
print(`GST(G): input a graph G=[V,E] output the set of all spanning trees.    `):

elif nops([args])=1 and op(1,[args])=GSF then
print(`GSF(G,S): input a graph G=[V,E] where V is a set of vertices and E is a set of edges, and a set of vertices S, which is a subset of V, output the set of spanning forests with |S| components such that each vertex in S is in different component.  `):

elif nops([args])=1 and op(1,[args])=OneStepG then
print(`OneStepG(G,S,v): input a graph G=[V,E], an arrow vertices set S and a vertex v, then we consider to delete v to get a recurrence relation for the number of spanning forests in G such that there are exactly |S| components and different vertices in S are in different components, output a [G,S]-multiset.  `):

elif nops([args])=1 and op(1,[args])=randgraph then
print(`randgraph(N,K): input two integers N and K, output a random graph with N vertices and K edges.  `):

elif nops([args])=1 and op(1,[args])=EdgeToGraph then
print(`EdgeToGraph(E): input a set of edges, output a graph G=[V,E] where we assume there is no isolated vertex.   `):

elif nops([args])=1 and op(1,[args])=GuessRec1 then
print(` GuessRec1(L,d): inputs a sequence L and tries to guess a recurrence operator with constant coefficients of order d satisfying it. It returns the initial d values and the operator as a list. For example try: `):
print(` GuessRec1([1,1,1,1,1,1],1);   `):

elif nops([args])=1 and op(1,[args])=GuessRec then
print(` GuessRec(L): inputs a sequence L and tries to guess a recurrence operator with constant coefficients  satisfying it. It returns the initial  values and the operator as a list. For example try: `):
print(` GuessRec([1,1,1,1,1,1]);   `):

elif nops([args])=1 and op(1,[args])=GuessSymRec1 then
print(` GuessSymRec1(L,d):input a sequence L and try to guess a recurrence operator with symmetric constant coefficients of order d. It returns the initial d values and the operator as a list. Here symmetric means the last number in the operator is -1, the first one = the (d-1)th, the 2nd = the (d-2)th, etc.For example try: `):
print(` GuessSymRec1([1,1,1,1,1,1],1);   `):

elif nops([args])=1 and op(1,[args])=GuessSymRec then
print(` GuessSymRec(L): input a sequence L and tried to guess a recurrence operator with symmetric constant coefficients. It returns the initial d values and the operator as a list. Here symmetric means the last number in the operator is -1, the first one = the (d-1)th, the 2nd = the (d-2)th, etc. For example try: `):
print(` GuessSymRec([1,1,1,1,1,1]);   `):

elif nops([args])=1 and op(1,[args])=CtoR then
print(` CtoR(S,t), the rational function, in t, whose coefficients are the members of the C-finite sequence S. For example, try:` ):
print(` CtoR([[1,1],[1,1]],t);  `):

elif nops([args])=1 and op(1,[args])=GridMN then
print(` GridMN(M,N): The Grid graph {0, ..., N-1)x{0,..., M-1} followed by the list of vertices. Try:` ):
print(` GridMN(3,5);  `):

elif nops([args])=1 and op(1,[args])=SpF then
print(` SpF(G,n,V0): The set of spanning forests of the graph G with set of vertices {1, ..., n} where the members of V0 belong to different components. Try:`):
print(` SpF({{1,2},{1,3},{2,3}},3,{1});  `):

elif nops([args])=1 and op(1,[args])=SpFn then
print(` SpFn(G,n,V0): The NUMBER of spanning forests of the graph G with set of vertices {1, ..., n} where the members of V0 belong to different components. Try:`):
print(` SpFn({{1,2},{1,3},{2,3}},3,{1});  `):

elif nops([args])=1 and op(1,[args])=GFGridKN then
print(` GFGridKN(k,t): input an integer k and a symbol t, output the generating function f(t) whose coefficient of t^n is the number of spanning trees in k by n grid graph.    `):

elif nops([args])=1 and op(1,[args])=DenomSFKN then
print(` DenomSFKN(k,t): input an integer k and a symbol t, output the denominator of the generating function f(t) whose coefficient of t^n is the number of spanning forests in k by n grid graph with two components such that 1 and k*n are in different components / the square of the denominator of GFGridKN(k,t).    `):

elif nops([args])=1 and op(1,[args])=GridH then
print(` GridH(m,n,t): input two integers m and n, and a symbol t, output the generating function f(t) whose coefficient of t^i is the number of spanning trees in the m by n grid graph which have i horizontal edges.   `):

elif nops([args])=1 and op(1,[args])=GridV then
print(` GridV(m,n,t): input two integers m and n, and a symbol t, output the generating function f(t) whose coefficient of t^i is the number of spanning trees in the m by n grid graph which have i vertical edges.   `):

elif nops([args])=1 and op(1,[args])=GridHV then
print(` GridHV(m,n,x,y): input two integers m and n, and a symbol t, output the generating function f(x,y) whose coefficient of x^i*y^j is the number of spanning trees in the m by n grid graph which have i horizontal edges and j vertical edges.   `):

elif nops([args])=1 and op(1,[args])=GridHFactorial then
print(` GridHFactorial(i,k,n): input three integers i,k and n, output a list of i-th factorial moment of GridH(k,j) for j from 1 to n.   `):

elif nops([args])=1 and op(1,[args])=GridVFactorial then
print(` GridVFactorial(i,k,n): input three integers i,k and n, output a list of i-th factorial moment of GridV(k,j) for j from 1 to n.   `):

elif nops([args])=1 and op(1,[args])=ListSTPreCal then
print(` ListSTPreCal(k): input an integer 2<=k<=7, and output a list of pre-calculated numbers which is long enough to get a recurrence via GuessSymRec for k by n grid graph.   `):

elif nops([args])=1 and op(1,[args])=GFGridPreCal then
print(` GFGridPreCal(k,t): input an integer 1<=k<=7 and output the pre-calculated generating function f(t) whose coefficient of t^n is the number of spanning trees in k by n grid graph.    `):

elif nops([args])=1 and op(1,[args])=DenomSFPreCal then
print(` DenomSFPreCal(k,t): input an integer 1<=k<=4 and a symbol t, output the pre-calculated denominator of the generating function f(t) whose coefficient of t^n is the number of spanning forests in k by n grid graph with two components such that 1 and k*n are in different components/the square of the denominator of GFGridKN(k,t).   `):

elif nops([args])=1 and op(1,[args])=GuessRat then
print(` GuessRat(L,degT,degB): input a list L and the degrees of top and button, try to find a rational function R(x)=P(x)/Q(x) with deg(P)<=degT, deg(Q)<=degB, such that R(i)=L(i).   `):

elif nops([args])=1 and op(1,[args])=DiagMatrix then
print(` DiagMatrix(a,b,c,p,q,n): input symbols a,b,c, two positive integers p and q, and a positive integer n > max(p,q) output a matrix as a list of lists where the diagonal element is a, in the first row after a it’s b[1],…,b[p], in the first column under a it’s c[1],…,c[q]. Try:  `):
print(` DiagMatrix(a,b,c,1,2,5); `):

elif nops([args])=1 and op(1,[args])=ExpandMatrix then
print(` ExpandMatrix(M): input a “general diagonal” matrix, do an expansion along its first row or column. (choose row if the number of nonzero elements in the first row is no more than that in the first column). output a set of [multiplicity, Matrix].   `):

elif nops([args])=1 and op(1,[args])=MatrixRecurrence then
print(` MatrixRecurrence(a,b,c,p,q): input symbols a,b,c and positive integers p,q, output a list L which represents the recurrence relation of the determinant of this matrix. For example, if we use D(n) to denote this kind of matrix of dimension n (large enough), then [a, -b*c] means D(n) = a*D(n-1) - b*c*D(n-2). Try:   `):
print(` MatrixRecurrence(a,b,c,1,2); `):


fi:
end:


ezraS:=proc()

if args=NULL then
 print(` The supporting procedures are: BR , Neis ,  UC , BRmn , Neismn , UCmn , BRkn , Neiskn , UCkn , BRL , NeisL , UCL , ST2n , SF2n , STkn  `):
 print(``):

else
ezra(args):
fi:

end:

###From Cfinite
#GuessRec1(L,d): inputs a sequence L and tries to guess
#a recurrence operator with constant cofficients of order d
#satisfying it. It returns the initial d values and the operator
# as a list. For example try:
#GuessRec1([1,1,1,1,1,1],1);
GuessRec1:=proc(L,d) local eq,var,a,i,n:

if nops(L)<=2*d+2 then
 print(`The list must be of size >=`, 2*d+3 ):
 RETURN(FAIL):
fi:

var:={seq(a[i],i=1..d)}:

eq:={seq(L[n]-add(a[i]*L[n-i],i=1..d),n=d+1..nops(L))}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
else
 RETURN([[op(1..d,L)],[seq(subs(var,a[i]),i=1..d)]]):
fi:

end:

#GuessRec(L): inputs a sequence L and tries to guess
#a recurrence operator with constant cofficients 
#satisfying it. It returns the initial  values and the operator
# as a list. For example try:
#GuessRec([1,1,1,1,1,1]);
GuessRec:=proc(L) local gu,d:

for d from 1 to trunc(nops(L)/2)-2 do
 gu:=GuessRec1(L,d):
 if gu<>FAIL then
   RETURN(gu):
 fi:
od:
FAIL:

end:
###End from Cfinite


#L(a,n): The Laplacian of the complete graph.
L:=proc(a,n)   local j1,i,j,T,gu,i1:

for i from 1 to n do
  for j from 1 to n do
    if i<j then
      T[i,j]:=a[i,j]:
    elif i>j then
      T[i,j]:=a[j,i]:
    else
      T[i,i]:=-add(a[i,j1],j1=i+1..n)-add(a[j1,i],j1=1..i-1):
    fi:
  od:
od:
gu:=[seq([seq(T[i,j],j=1..n)],i=1..n)]:
end:


#Lg(a,n,G): The Laplacian of the graph G.
Lg:=proc(a,n,G)   local i,j,T,gu:

gu:=L(a,n):

for i from 1 to n do
  for j from i+1 to n do
    if not member({i,j},G) then
       gu:=subs(a[i,j]=0,gu):
    fi:
  od:
od:
gu:
end:


#JC(a,n): the joint conductance between nodes 1 and n.
JC:=proc(a,n) local gu,i1,gu1,gu2:
gu:=L(a,n):
gu1:=[seq([op(1..n-1,gu[i1])],i1=1..n-1)]:
gu2:=[seq([op(2..n-1,gu[i1])],i1=2..n-1)]:
normal(-det(gu1)/det(gu2)):
end:


#JCpq(a,n,p,q): the joint conductance between nodes p and q.
JCpq:=proc(a,n,p,q) local i,gu1,gu2,p1,q1:
if p<1 or p>n or q<1 or q>n or p=q then
  return FAIL:
fi:
if p>q then
  p1:=q:
  q1:=p:
else
  p1:=p:
  q1:=q:
fi:

gu1:=L(a,n):
gu1:=[op(1..q1-1, gu1),op(q1+1..n,gu1)]:
gu1:=[seq([op(1..q1-1,gu1[i]), op(q1+1..n,gu1[i])],i=1..n-1)]:
gu2:=gu1:
gu2:=[op(1..p1-1, gu2),op(p1+1..n-1,gu2)]:
gu2:=[seq([op(1..p1-1,gu2[i]), op(p1+1..n-1,gu2[i])],i=1..n-2)]:

normal(-det(gu1)/det(gu2)):
end:


#JCg(a,n,G): the joint conductance between nodes 1 and n in the graph G. Try:
#JCg(a,3,{{1,2},{2,3}});
JCg:=proc(a,n,G) local gu,i1,gu1,gu2:
gu:=Lg(a,n,G):
gu1:=[seq([op(1..n-1,gu[i1])],i1=1..n-1)]:
gu2:=[seq([op(2..n-1,gu1[i1])],i1=2..n-1)]:
normal(-det(gu1)/det(gu2)):
end:


#JCgpq(a,n,G,p,q): the joint conductance between nodes p and q in the graph G. Try:
#JCgpq(a,3,{{1,2},{2,3}},3,2);
JCgpq:=proc(a,n,G,p,q) local i,gu1,gu2,p1,q1:
if p<1 or p>n or q<1 or q>n or p=q then
  return FAIL:
fi:
if p>q then
  p1:=q:
  q1:=p:
else
  p1:=p:
  q1:=q:
fi:

gu1:=Lg(a,n,G):
gu1:=[op(1..q1-1, gu1),op(q1+1..n,gu1)]:
gu1:=[seq([op(1..q1-1,gu1[i]), op(q1+1..n,gu1[i])],i=1..n-1)]:
gu2:=gu1:
gu2:=[op(1..p1-1, gu2),op(p1+1..n-1,gu2)]:
gu2:=[seq([op(1..p1-1,gu2[i]), op(p1+1..n-1,gu2[i])],i=1..n-2)]:

normal(-det(gu1)/det(gu2)):
end:


#JCnum(n,G): the joint conductance between nodes 1 and n in the graph G where all resistances are with 1 Ohm. Try:
#JCnum(3,{{1,2},{2,3}});
JCnum:=proc(n,G) local a,gu,i1,gu1,gu2,i,j:
gu:=Lg(a,n,G):

gu:=subs({seq(seq(a[i,j]=1,j=i+1..n),i=1..n)},gu):

gu1:=[seq([op(1..n-1,gu[i1])],i1=1..n-1)]:
gu2:=[seq([op(2..n-1,gu1[i1])],i1=2..n-1)]:
normal(-det(gu1)/det(gu2)):
end:

#JCnumND(n,G): The numerator and denominator of JCnum(n,G). Try:
#JCnumND(3,{{1,2},{2,3}});
JCnumND:=proc(n,G) local a,gu,i1,gu1,gu2,i,j:
gu:=Lg(a,n,G):

gu:=subs({seq(seq(a[i,j]=1,j=i+1..n),i=1..n)},gu):

gu1:=[seq([op(1..n-1,gu[i1])],i1=1..n-1)]:
gu2:=[seq([op(2..n-1,gu1[i1])],i1=2..n-1)]:
[-det(gu1),det(gu2)]:
end:


#JCnumpq(n,G,p,q): the joint conductance between nodes p and q in the graph G where all resistances are with 1 Ohm. Try:
#JCnumpq(3,{{1,2},{2,3}},3,2);
JCnumpq:=proc(n,G,p,q) local a,gu,i,j,gu1,gu2,p1,q1:
gu:=Lg(a,n,G):
gu:=subs({seq(seq(a[i,j]=1,j=i+1..n),i=1..n)},gu):

if p<1 or p>n or q<1 or q>n or p=q then
  return FAIL:
fi:
if p>q then
  p1:=q:
  q1:=p:
else
  p1:=p:
  q1:=q:
fi:

gu1:=gu:
gu1:=[op(1..q1-1, gu1),op(q1+1..n,gu1)]:
gu1:=[seq([op(1..q1-1,gu1[i]), op(q1+1..n,gu1[i])],i=1..n-1)]:
gu2:=gu1:
gu2:=[op(1..p1-1, gu2),op(p1+1..n-1,gu2)]:
gu2:=[seq([op(1..p1-1,gu2[i]), op(p1+1..n-1,gu2[i])],i=1..n-2)]:

normal(-det(gu1)/det(gu2)):


end:


#BR(i,n): inputes an integer i between 0 to 2^n-1 outputs its binary representation, a 0-1 vector of length n. Try:
#BR(3,3);
BR:=proc(i,n) local gu,i1:
gu:=convert(i,base,2):
[0$(n-nops(gu)),seq(gu[nops(gu)-i1+1],i1=1..nops(gu))]:
end:


#Neis(i,j,n): are vertices i and j of the n-dim unit cube neighbors?
Neis:=proc(i,j,n)  local mu,m:
mu:=BR(j-1,n)-BR(i-1,n):
mu:=[seq(abs(m), m in mu)]:
evalb(convert(mu,`+`)=1):
end:


#UC(n): the graph of the n-dim unit cube.
UC:=proc(n) local i,j,gu:
gu:={}:

for i from 1 to 2^n do
  for j from i+1 to 2^n do
    if Neis(i,j,n) then
      gu:= gu union {{i,j}}:
    fi:
  od:
od:

gu:

end:


#Rn(n): the formula for the joint resistance between opposite vertices via the general formula.
Rn:=proc(n):
1/JCnum(2^n,UC(n)):
end:


#RnF(n): the joint resistance between opposite vertices of the n-dim unit-cube.
#according to the formula in Pippenger's paper. Try:
#RnF(5);
RnF:=proc(n) local k:
1/n*add(1/binomial(n-1,k),k=0..n-1):
end:


#BRmn(i,m,n):inputs an integer i between 0 and m*n-1 outputs its location [x,y] (0<=x<=m-1, 0<=y<=n-1) such that n*x+y=i.
BRmn:=proc(i,m,n) local x,y:
y:=i mod n:
x:=(i-y)/n:
[x,y]:
end:


#Neismn(i,j,m,n): are vertices i and j of m by n grid neighbors?
Neismn:=proc(i,j,m,n) local gu,g:
gu:=BRmn(j-1,m,n)-BRmn(i-1,m,n):
gu:=[seq(abs(g), g in gu)]:
evalb(convert(gu,`+`)=1):
end:


#UCmn(m,n): the graph of a m by n grid (# of nodes are m*n).
UCmn:=proc(m,n) local i,j,gu:
gu:={}:

for i from 1 to m*n do
  for j from i+1 to m*n do
    if Neismn(i,j,m,n) then
      gu:=gu union {{i,j}}:
    fi:
  od:
od:
gu:
end:


#Rmn(m,n): the formula for the joint resistance between opposite vertices of a m by n grid via the general formula.
Rmn:=proc(m,n):
1/JCnum(m*n, UCmn(m,n)):
end:
#seq(Rmn(k,k),k=2..11) = 1, 3/2, 13/7, 47/22, 1171/495, 6385/2494, 982871/360161, 441083/153254, 427854195/142065451, 4769986941/1522703822


#Rmnpq(m,n,p,q): the formula for joint resistance between node p and q in a m by n grid via the general formula.
Rmnpq:=proc(m,n, p,q):
1/JCnumpq(m*n, UCmn(m,n),p,q):
end:


#BRkn(i,k,n): inputes an integer i between 0 to k^n-1 outputs its k-nary representation, a 0,1,…,k-1 vector of length n. Try: BRkn(5,3,3);
BRkn:=proc(i,k,n) local gu,i1:
gu:=convert(i,base,k):
[0$(n-nops(gu)),seq(gu[nops(gu)-i1+1],i1=1..nops(gu))]:
end:


#Neiskn(i,j,k,n): are vertices i and j of [0,k-1]^n neighbors?
Neiskn:=proc(i,j,k,n) local gu,g:
gu:=BRkn(j-1,k,n)-BRkn(i-1,k,n):
gu:=[seq(abs(g), g in gu)]:
evalb(convert(gu,`+`)=1):
end:


#UCkn(k,n): the graph of [0,k-1]^n with nodes labelled from 1 to k^n.
UCkn:=proc(k,n) local i,j,gu:
gu:={}:

for i from 1 to k^n do
  for j from i+1 to k^n do
    if Neiskn(i,j,k,n) then
      gu:=gu union {{i,j}}:
    fi:
  od:
od:

gu:

end:


#Rkn(k,n): the formula for the joint resistance between opposite vertices of [0,k-1]^n via the general formula.
Rkn:=proc(k,n):
1/JCnum(k^n, UCkn(k,n)):
end:
#Of course seq(Rkn(k, 2), k = 2 .. 11) is the same as seq(Rmn(k, k), k = 2 .. 11).

#BRL(i,L): inputs an integer i between 0 to L[1]*…*L[n]-1 (n is the length of the list L) and the list L of positive integers output the n-vector [a_1, a_2, … , a_n], where 0<=a_j<L[j].
BRL:=proc(i,L) local n,j,gu,mu,temp,re:
n:=nops(L):
gu:=[0$n]:
mu:=i:
temp:=1:
for j from n to 1 by -1 do
  re:=mu mod (temp*L[j]):
  mu:=mu-re:
  re:=re/temp:
  temp:=temp*L[j]:
  gu[j]:=re:
od:
gu:
end:


#NeisL(i,j,L): are vertices i and j of L-shape hyper-rectangle neighbors?
NeisL:=proc(i,j,L) local gu,g:
gu:=BRL(j-1,L)-BRL(i-1,L):
gu:=[seq(abs(g), g in gu)]:
evalb(convert(gu,`+`)=1):
end:


#UCL(L): the graph of L-shape hyper-rectangle.
UCL:=proc(L) local i,j,gu,upper:
gu:={}:
upper:=mul(L[i],i=1..nops(L)):
for i from 1 to upper do
  for j from i+1 to upper do
    if NeisL(i,j,L) then
      gu:=gu union {{i,j}}:
    fi:
  od:
od:
gu:
end:


#RL(L): the formula for the joint resistance between opposite vertices of L-shape hyper-rectangle via the general formula.
RL:=proc(L) local i,n:
n:=mul(L[i],i=1..nops(L)):
1/JCnum(n, UCL(L)):
end:
#eq(RL([k, k]), k = 2 .. 11) is the same as seq(Rmn(k, k), k = 2 .. 11).


#ST2n(n): the number of spanning trees in the 2*n grid graph.
ST2n:=proc(n) 
option remember:
if n<1 then
  return FAIL:
elif n=1 then
  return 1:
elif n=2 then
  return 4:
else
  return 4*ST2n(n-1)-ST2n(n-2):
fi:
end:


#SF2n(n): the number of spanning forests with 2 trees such that the two vertices on the same side are in distinct trees.
SF2n:=proc(n)
option remember:
if n<1 then
  return FAIL:
elif n=1 then
  return 1:
elif n=2 then
  return 3:
else
  return 4*SF2n(n-1)-SF2n(n-2):
fi:
end:


#R2n(n): the resistance of 2 by n grid with 1 Ohm each edge.
R2n:=proc(n)
SF2n(n)/ST2n(n):
end:
#the result is the same as seq(Rmnpq(2, n, 1, n+1), n = 2 .. 10).


#STkn(k,n): the number of spanning trees in the k*n grid graph for 2<=k<=5 by using generating function.
STkn:=proc(k,n) local F,x:
if k=2 then
  F:=x/(1-4*x+x^2):
elif k=3 then
  F:=3*x*(1+49*x+1152*x^2)/(1+24*x-24*x^2+x^3):
elif k=4 then
  F:=16*x*(1+12*x+x^2)/(1-204*x+1190*x^2-204*x^3+x^4):
elif k=5 then
  F:=125*x*(1+4656*x+10616686*x^2+23432228161*x^3+51714958501250*x^4)/(1+2255*x-105985*x^2+105985*x^3-2255*x^4+x^5):
fi:
coeff(taylor(F,x=0,n+1),x,n):
end:


#ReArrow(arrow): input a set of arrow vertices, empty it if there is only one element. 
ReArrow:=proc(arrow)
if nops(arrow)=1 then
  return {}:
fi:
return arrow:
end:

#OneRecurrence(L): input a list, which represents a data structure of k by n grid (k fixed, e.g. [2,0,2,{}] and outputs a set of lists where each list is a data structure following by a number, e.g. {[[2,0,1,{}],2],[[2,0,1,{[1,2],[2,1]}],1]}. It means the number of spanning trees in L is the sum of the corresponding number times the number of spanning trees in those data structures represented by the list. Here 0 means the length is n, -1 means the length is n-1.

OneRecurrence:=proc(L) local k,n,arrow,p,S:
k:=L[1]:
n:=L[2]:
p:=L[3]:
arrow:=ReArrow(L[4]):
S:={}:
if (p>1) and (not [k-p+1,1-n] in arrow) then
  if ([k-p+1, 2-n] in arrow) and ([k-p+2, 1-n] in arrow) then
    S:=S union {[[k,n,p-1,ReArrow(arrow)],2]}:
  else
    S:=S union {[[k,n,p-1,ReArrow(arrow)],2]}:
    arrow:=arrow union {[k-p+1, 2-n],[k-p+2, 1-n]}:
    S:=S union {[[k,n,p-1,ReArrow(arrow)],1]}:
  fi:
elif (p>1) and ([k-p+1,1-n] in arrow) then
  if ([k-p+1, 2-n] in arrow) and ([k-p+2, 1-n] in arrow) then
    arrow:=arrow minus {[k-p+1,1-n]}:
    S:=S union {[[k,n,p-1,ReArrow(arrow)],1]}:
  elif ([k-p+1, 2-n] in arrow) and (not [k-p+2, 1-n] in arrow) then
    arrow:=arrow minus {[k-p+1,1-n]}:
    S:=S union {[[k,n,p-1,ReArrow(arrow)],1]}:
    arrow:=arrow union {[k-p+2, 1-n]}:
    S:=S union {[[k,n,p-1,ReArrow(arrow)],1]}:
  elif (not [k-p+1, 2-n] in arrow) and ([k-p+2, 1-n] in arrow) then
    arrow:=arrow minus {[k-p+1,1-n]}:
    S:=S union {[[k,n,p-1,ReArrow(arrow)],1]}:
    arrow:=arrow union {[k-p+1, 2-n]}:
    S:=S union {[[k,n,p-1,ReArrow(arrow)],1]}:
  else
    arrow:=arrow minus {[k-p+1,1-n]}:
    S:=S union {[[k,n,p-1,ReArrow(arrow)],1]}:
    arrow:=arrow union {[k-p+2, 1-n]}:
    S:=S union {[[k,n,p-1,ReArrow(arrow)],1]}:
    arrow:=arrow minus {[k-p+2, 1-n]}:
    arrow:=arrow union {[k-p+1, 2-n]}:
    S:=S union {[[k,n,p-1,ReArrow(arrow)],1]}:
    arrow:=arrow union {[k-p+2, 1-n]}:
    S:=S union {[[k,n,p-1,ReArrow(arrow)],1]}:
  fi:
elif (p=1) and (not [k-p+1,1-n] in arrow) then
  S:= S union {[[k,n-1,k,ReArrow(arrow)],1]}:
elif (p=1) and ([k-p+1,1-n] in arrow) then
  if [k-p+1,2-n] in arrow then
    arrow:=arrow minus {[k-p+1,1-n]}:
    S:=S union {[[k,n-1,k,ReArrow(arrow)],1]}:
  else
    arrow:=arrow minus {[k-p+1,1-n]}:
    S:=S union {[[k,n-1,k,ReArrow(arrow)],1]}:
    arrow:=arrow union {[k-p+1,2-n]}:
    S:=S union {[[k,n-1,k,ReArrow(arrow)],1]}:
  fi:
fi:
return S:
end:
 

#OneList(L): List only for the output, similar to the above procedure.  
OneList:=proc(L) local S,C,s:
S:=OneRecurrence(L):
C:={}:
for s in S do
  C:=C union {s[1]}:
od:
C:
end:
 

#OriginalCase(L): input a list L=[k,n,p,arrow] as a data structure, depending on n, output its similar original case.
OriginalCase:=proc(L) local k,n,p,arrow,i,newarrow:
k:=L[1]:
n:=L[2]:
p:=L[3]:
arrow:=L[4]:
newarrow:={}:
for i in arrow do
  newarrow:=newarrow union {(i+[0,n])}:
od:
return [k,0,p,newarrow]:
end:


#AllCase(L): input a list L as a data structure, output a set of all cases when we apply OneRecurrence to L and its output(whose n=0 but not negative).
AllCase:=proc(L) local S,gu,Candidate,NC,mu:
S:={L}:
Candidate:={seq(OriginalCase(gu), gu in OneList(L))}:
while Candidate minus S <> {} do
  NC:=Candidate minus S:
  S:=S union Candidate:
  Candidate:={}:
  for mu in NC do
    Candidate:=Candidate union {seq(OriginalCase(gu), gu in OneList(mu))}:
  od:
od:
S:
end:


#AllCaseL(L): input a list L as a data structure, output a list of all cases when we apply OneRecurrence to L and its output(whose n=0 but not negative).
AllCaseL:=proc(L) local S,SL,gu,Candidate,NC,mu:
S:=[L]:
SL:={L}:
Candidate:={seq(OriginalCase(gu), gu in OneList(L))}:
while Candidate minus SL <> {} do
  NC:=Candidate minus SL:
  SL:=SL union Candidate:
  S:=[op(S),op(NC)]:
  Candidate:={}:
  for mu in NC do
    Candidate:=Candidate union {seq(OriginalCase(gu), gu in OneList(mu))}:
  od:
od:
S:
end:
  
#Coefficient(S,L): input a set S of lists [L, number], and a list L, output the number corresponding to the list L.
Coefficient:=proc(S,L) local s:
for s in S do
  if s[1]=L then
    return s[2]:
  fi:
od:
return 0:
end:


#FinalRecurrence(L): input a list L as a data structure, output a set of lists where each sublist is [data structure of -1, number].
FinalRecurrence:=proc(L) local S,index,s,Candidate,candidate,n,co,gu,mu:
S:=OneRecurrence(L):
index:=mul(s[1][2], s in S):
while index = 0 do
  for s in S do
    if s[1][2]=0 then
      n:=s[2]:
      Candidate:=OneRecurrence(s[1]):
      S:=S minus {s}:
      for candidate in Candidate do
        gu:=candidate[1]:
        mu:=candidate[2]:
        co:=Coefficient(S,gu):
        S:=S union {[gu,mu*n+co]}:
      od:
    fi:
  od:
  index:=mul(s[1][2], s in S):
od:
S:
end:
          

#FinalVector(S,L): input a list of lists S and a list L as a data structure, output a list LL such that LL[i] is the coefficient of S[i] in FinalRecurrence(L). (regardless of n=0 or -1).
FinalVector:=proc(S,L) local n,gu,F,f,i:
n:=nops(S):
gu:=[0$n]:
F:=FinalRecurrence(L):
for f in F do
  for i from 1 to n do
    if S[i]=OriginalCase(f[1]) then
      gu[i]:=f[2]:
      break:
    fi:
  od:
od:
gu:
end:


#FinalMatrix(L): input a list L as a data structure, output a transition matrix. The order of the entry depends on the position in the list AllCaseL(L) where L will always be the first one. 
FinalMatrix:=proc(L) local S,s,M:
S:=AllCaseL(L):
M:=[]:
for s in S do
  M:=[op(M), FinalVector(S,s)]:
od:
M:
end:


#FinalGF(L,x): input a list L as a data structure, output its generating function such that the coefficient of x^n for the generating function means the number of spanning trees/forests in this data structure with size k by n. 
FinalGF:=proc(L,x) local M,n,var,A,eq,sol,i,j,gu:
M:=FinalMatrix(L):
n:=nops(M):
var:={seq(A[i],i=1..n)}:
eq:={A[1]-1-x*add(M[1][j]*A[j],j=1..n)}:
for i from 2 to n do
  eq:=eq union {A[i]-x*add(M[i][j]*A[j],j=1..n)}:
od:
gu:=solve(eq, var)[1]:
subs(gu, A[1]):
end:




#CP(L): input a list L as a data structure, output its characteristic polynomial for its transition matrix. 
CP:=proc(L) 
CharacteristicPolynomial(FinalMatrix(L),x):
end:
 
 
#NeighborV(G,S): input a graph G=[V,E] and a set S of vertices, output a set of all S’s neighbors in G.
NeighborV:=proc(G,S) local V,E,N,u,v:
V:=G[1]:
E:=G[2]:
N:={}:
for v in S intersect V do
  for u in V do
    if {u,v} in E then
      N:=N union {u}:
    fi:
  od:
od:
N:
end:

   
#IsConnected(G): input a graph G=[V,E] (V is a set of vertex, E is a set of edges {i,j}), if it is connected output true, otherwise false.
IsConnected:=proc(G) local V,E,V1,edge,s,S,T,N,NS:
V:=G[1]:
E:=G[2]:

V1:={}:
for edge in E do
  V1:=V1 union edge:
od:

if V1 minus V <> {} then
  return FAIL:
fi:

if nops(V)=1 or nops(V)=0 then
  return true:
fi:

s:=V[1]:
S:={s}:
T:={}:
for edge in E do
  if s in edge then
    T:=T union edge:
  fi:
od:

if T={} then
  return false:
fi:

while S<>T do
  N:=T minus S:
  S:=T:
  for edge in E do
    if edge[1] in N or edge[2] in N then
      T:=T union edge:
    fi:
  od:
od:
if S=V then
  return true:
else
  return false:
fi:
end:
  

#DegreeG(G,n): input a graph G=[V,E] and a vertex n, output the degree of the vertex. 
DegreeG:=proc(G,n) local co, edge:
co:=0:
for edge in G[2] do
  if n in edge then
    co:=co+1:
  fi:
od:
co:
end:


#AllComponent(G): input a graph G=[V,E] with one or more components, output a set of these connected subgraphs.
AllComponent:=proc(G) local V,E,V1,s,S,T,Edge,edge,N:
V:=G[1]:
E:=G[2]:

V1:={}:
for edge in E do
  V1:=V1 union edge:
od:

if V1 minus V <> {} then
  return FAIL:
fi:

if IsConnected(G) then
  return {G}:
fi:

s:=V[1]:
S:={s}:
T:={}:
Edge:={}:
for edge in E do
  if s in edge then
    T:=T union edge:
    Edge:=Edge union {edge}:
  fi:
od:

if T={} then
  return AllComponent([V minus {s}, E]) union {[{s},{}]}:
fi:

while S<>T do
  N:=T minus S:
  S:=T:
  for edge in E do
    if edge[1] in N or edge[2] in N then
      T:=T union edge:
      Edge:=Edge union {edge}:
    fi:
  od:
od:
return AllComponent([V minus S, E minus Edge]) union {[S,Edge]}:
end:





#NumberComponent(G): input a graph G=[V,E], output the number of components in the graph.
NumberComponent:=proc(G) local V,E,s,S,T,Edge,edge,N,V1:
V:=G[1]:
E:=G[2]:

V1:={}:
for edge in E do
  V1:=V1 union edge:
od:

if V1 minus V <> {} then
  return FAIL:
fi:

if V={} then
  return 0:
elif E={} then
  return nops(V):
fi:

s:=V[1]:
S:={s}:
T:={}:
Edge:={}:
for edge in E do
  if s in edge then
    T:=T union edge:
    Edge:=Edge union {edge}:
  fi:
od:

if T={} then
  return 1+NumberComponent([V minus S, E]):
fi:

while S<>T do
  N:=T minus S:
  S:=T:
  for edge in E do
    if edge[1] in N or edge[2] in N then
      T:=T union edge:
      Edge:=Edge union {edge}:
    fi:
  od:
od:
return 1+NumberComponent([V minus S, E minus Edge]):
end:


#IsAcyclic(G): input a graph G=[V,E], output true if the graph contains no cycle, otherwise output false. 
IsAcyclic:=proc(G) local Comp,C,V,E,S,T,edge,t,s:
Comp:=AllComponent(G):
for C in Comp do
  V:=C[1]:
  E:=C[2]:
  s:=V[1]:
  S:={s}:
  T:=NeighborV(C,S):
  while T<>{} do
    S:=S union T:
    for t in T do
      if nops(NeighborV(C,{t}) intersect S) > 1 then
        return false:
      fi:
    od:
    T:=NeighborV(C,T) minus S:
  od:
od:
return true:
end:
  

#GST(G): input a graph G=[V,E] output the set of all spanning trees. 
GST:=proc(G) local V,E,S,NG,NS,ns,L,L1,L2,l1,l2,edge,flag:
option remember:
if not IsConnected(G) then
  return {}:
fi:
S:={}:
V:=G[1]:
E:=G[2]:

if nops(V)=1 then
  return {G}:
fi:

flag:=0:

for edge in E do
  if DegreeG(G,edge[1])=1 then
    NG:=[V minus {edge[1]}, E minus {edge}]:
    NS:=GST(NG):
    S:=S union {seq([ns[1] union {edge[1]},ns[2] union {edge}], ns in NS)}:
    flag:=1:
    break:
  elif DegreeG(G,edge[2])=1 then
    NG:=[V minus {edge[2]}, E minus {edge}]:
    NS:=GST(NG):
    S:=S union {seq([ns[1] union {edge[2]},ns[2] union {edge}], ns in NS)}:
    flag:=1:
    break:
  elif DegreeG(G,edge[1])>1 and DegreeG(G,edge[2])>1 then
    NG:=[V, E minus {edge}]:
    if not IsConnected(NG) then
      L:=AllComponent(NG):
      L1:=L[1]:
      L2:=L[2]:
      S:=S union {seq(seq([l1[1] union l2[1], l1[2] union l2[2] union {edge}], l1 in GST(L1)),l2 in GST(L2))}:
      flag:=1:
      break:
    fi:
  fi:
od:

if flag=0 then
  for edge in E do
    NG:=[V, E minus {edge}]:
    S:=S union GST(NG):
  od:
fi:
S:
end:


#gst(G): input a graph G=[V,E] output the number of all spanning trees.
gst:=proc(G) local V,E,S,NG,NS,ns,L,L1,L2,l1,l2,edge,flag,SS:
option remember:
if not IsConnected(G) then
  return 0:
fi:
S:=0:
V:=G[1]:
E:=G[2]:

if nops(V)=1 then
  return 1:
fi:

flag:=0:

for edge in E do
  if DegreeG(G,edge[1])=1 then
    NG:=[V minus {edge[1]}, E minus {edge}]:
    NS:=gst(NG):
    S:=S + NS:
    flag:=1:
    break:
  elif DegreeG(G,edge[2])=1 then
    NG:=[V minus {edge[2]}, E minus {edge}]:
    NS:=gst(NG):
    S:=S + NS:
    flag:=1:
    break:
  elif DegreeG(G,edge[1])>1 and DegreeG(G,edge[2])>1 then
    NG:=[V, E minus {edge}]:
    if not IsConnected(NG) then
      L:=AllComponent(NG):
      L1:=L[1]:
      L2:=L[2]:
      S:=S + gst(L1)*gst(L2):
      flag:=1:
      break:
    fi:
  fi:
od:    

SS:={}:
if flag=0 then
  for edge in E do
    NG:=[V, E minus {edge}]:
    SS:=SS union GST(NG):
  od:
  S:=nops(SS):
fi:
S:
end: 


#LowestDegree(G,S): input a graph G=[V,E] where V is a set of vertices and E is a set of edges, and S a subset of vertices, output the vertex in S which has the lowest degree. 
LowestDegree:=proc(G,S) local V,E,s,low:
V:=G[1]:
E:=G[2]:
if S minus V <> {} or S={} then
  return FAIL:
fi:
low:=S[1]:
for s in S do
  if DegreeG(G,s) < DegreeG(G,low) then
    low:=s:
  fi:
od:
low:
end:


#GSF(G,S): input a graph G=[V,E] where V is a set of vertices and E is a set of edges, and a set of vertices S, which is a subset of V, output the set of spanning forests with |S| components such that each vertex in S is in different component. 
GSF:=proc(G,S) local V,E,V1,edge,v,low,lowneib,P,p,SF,gu,neighbors,neighboredge,existedge,candidate,mu:
option remember:
V:=G[1]:
E:=G[2]:
V1:={}:
SF:={}:

for edge in E do
  V1:=V1 union edge:
od:

if V1 minus V <> {} then
  return FAIL:
elif V=S then
  return {[S,{}]}:
fi:

for low in V minus S do
  neighbors:={}:
  neighboredge:={}:
  if DegreeG(G,low)=0 then
    return {}:
  fi:

  for v in V do
    if {v, low} in E then
      neighbors:=neighbors union {v}:
      neighboredge:=neighboredge union {{v,low}}:
    fi:
  od:

  for P in powerset(neighbors) do
    if nops(P)>=1 and nops(P intersect S)<=1 then
      existedge:={seq({p,low},p in P)}:
      gu:=GSF([V minus {low}, E minus neighboredge],S):
      for mu in gu do
        candidate:=[mu[1] union {low}, mu[2] union existedge]:
        if NumberComponent(candidate)=nops(S) and IsAcyclic(candidate) then
            SF:=SF union {candidate}:
        fi:
      od:
    fi:
  od:
od:

return SF:
end:


#OneStepG(G,S,v): input a graph G=[V,E], an arrow vertices set S and a vertex v, then we consider to delete v to get a recurrence relation for the number of spanning forests in G such that there are exactly |S| components and different vertices in S are in different components, output a [G,S]-multiset. 
OneStepG:=proc(G,S,v) local V,E,L,N,NV,nv,NE,n,gu,mu:
V:=G[1]:
E:=G[2]:
if not v in V or S minus V <> {} then 
  return FAIL:
fi:
L:=[]:
NV:=NeighborV(G,{v}):
NE:={seq({v,nv},nv in NV)}:
for N in powerset(NV) do
  if v in S and N intersect S = {} then
    L:=[op(L), [[V minus {v}, E minus NE], S union N minus {v}]]:
  elif not v in S and  nops(N intersect S)=1 then     #>1 is impossible
    L:=[op(L), [[V minus {v}, E minus NE], S union N]]:
  elif not v in S and N intersect S = {} then
    gu:=[seq(S union N minus {n}, n in N)]:
    L:=[op(L), seq([[V minus {v}, E minus NE], mu], mu in gu)]:
  fi:
od:
convert(L, multiset):
end:
    
      
#randgraph(N,K): input two integers N and K, output a random graph with N vertices and K edges. 
randgraph:=proc(N,K) local V,E,i,ra,a,b:
V:={seq(i,i=1..N)}:
E:={}:
ra:=rand(1..N):
while nops(E) < K do
  a:=ra():
  b:=ra():
  if a<>b then
    E:=E union {{a,b}}:
  fi:
od:
[V,E]:
end: 
    

#EdgeToGraph(E): input a set of edges, output a graph G=[V,E] where we assume there is no isolated vertex. 
EdgeToGraph:=proc(E) local V,edge:
V:={}:
for edge in E do
  V:=V union edge:
od:
[V,E]:
end:


#GuessSymRec1(L,d):input a sequence L and try to guess a recurrence operator with symmetric constant coefficients of order d. It returns the initial d values and the operator as a list. Here symmetric means the last number in the operator is -1, the first one = the (d-1)th, the 2nd = the (d-2)th, etc.For example try:
#GuessSymRec1([1,1,1,1,1,1],1);

GuessSymRec1:=proc(L,d) local eq, var, a, i, n,de:
de:=floor(d/2):
if nops(L)<=d+de then
  print(`The list must be of size >=floor(3*d/2)+1`):
  return(FAIL):
fi:

var:={seq(a[i],i=1..de)}:
eq:={seq(L[n]-add(a[i]*L[n-i],i=1..de)-add(a[d-i]*L[n-i],i=de+1..d-1)+L[n-d], n=d+1..nops(L))}:

var:=solve(eq, var):

if var=NULL then
  return(FAIL):
else
  return [[op(1..d,L)],[seq(subs(var,a[i]),i=1..de), seq(subs(var, a[d-i]), i=de+1..d-1), -1]]:
fi:

end:


#GuessSymRec(L): input a sequence L and tried to guess a recurrence operator with symmetric constant coefficients. It returns the initial d values and the operator as a list. Here symmetric means the last number in the operator is -1, the first one = the (d-1)th, the 2nd = the (d-2)th, etc. For example try:
#GuessSymRec([1,1,1,1,1,1]);

GuessSymRec:=proc(L) local gu,d:

for d from 1 to nops(L)*2/3 do
  gu:=GuessSymRec1(L,d):
  if gu<>FAIL then
    return gu:
  fi:
od:
FAIL:
end:


#CtoR(S,t), the rational function, in t, whose coefficients
#are the members of the C-finite sequence S. For example, try:
#CtoR([[1,1],[1,1]],t);

CtoR:=proc(S,t) local D1,i,N1,L1,f,f1,L:
if not (type(S,list) and  nops(S)=2 and type(S[1],list) and type(S[2],list)
        and nops(S[1])=nops(S[2]) and type(t, symbol) ) then
   print(`Bad input`):
   RETURN(FAIL):
fi:

D1:=1-add(S[2][i]*t^i,i=1..nops(S[2])):
N1:=add(S[1][i]*t^(i-1),i=1..nops(S[1])):
L1:=expand(D1*N1):
L1:=add(coeff(L1,t,i)*t^i,i=0..nops(S[1])-1):
f:=L1/D1:
L:=degree(D1,t)+10:
f1:=taylor(f,t=0,L+1):
if expand([seq(coeff(f1,t,i),i=0..L)])<>expand(SeqFromRec(S,L+1)) then
print([seq(coeff(f1,t,i),i=0..L)],SeqFromRec(S,L+1)):
 RETURN(FAIL):
else
 RETURN(f):
fi:
end:


#SeqFromRec(S,N): Inputs S=[INI,ope]
#where INI is the list of initial conditions, a ope a list of
#size L, say, and a recurrence operator ope, codes a list of
#size L, finds the first N0 terms of the sequence satisfying
#the recurrence f(n)=ope[1]*f(n-1)+...ope[L]*f(n-L).
#For example, for the first 20 Fibonacci numbers,
#try: SeqFromRec([[1,1],[1,1]],20);
SeqFromRec:=proc(S,N) local gu,L,n,i,INI,ope:
INI:=S[1]:ope:=S[2]:
if not type(INI,list) or not type(ope,list) then
 print(`The first two arguments must be lists `):
 RETURN(FAIL):
fi:
L:=nops(INI):
if nops(ope)<>L then
 print(`The first two arguments must be lists of the same size`):
RETURN(FAIL):
fi:

if not type(N,integer) then
 print(`The third argument must be an integer`, L):
 RETURN(FAIL):
fi:

if N<L then
 RETURN([op(1..N,INI)]):
fi:

gu:=INI:

for n from 1 to N-L do
 gu:=[op(gu),expand(add(gu[nops(gu)-i+1]*ope[i],i=1..L))]:
od:

gu:
end:


#GridMN(M,N): The Grid graph {0, ..., N-1)x{0,..., M-1} followed by the list of vertices. Try:
#GridMN(3,5);
GridMN:=proc(M,N) local G,gu,i,j,T:
G:={
#horizontal edges
seq(seq({[i,j],[i+1,j]},i=0..N-2),j=0..M-1),
#vertical edges
seq(seq({[i,j],[i,j+1]},i=0..N-1),j=0..M-2)
}:

gu:=[seq(seq([i,j],j=0..M-1),i=0..N-1)]:

for i from 1 to nops(gu) do
 T[gu[i]]:=i:
od:

G:=subs({seq(gu[i]=T[gu[i]],i=1..nops(gu))},G):

[G,gu]:

end:



#Mat11(G,n,a): Given a graph G on {1,2, ...,n} described as a set of edges, finds 
#the matrix whose determinant would give the Laplacian matrix, whose determinant of any n-1 by n-1 principal minor would
#give the sum of the weights of all the spanning trees. Try
#Mat11({{1,2},{2,3},{3,4},{4,1}},4,a);
Mat11:=proc(G,n,a) local M,i,j,lu:

if G minus {seq(seq({i,j},i=1..j-1),j=1..n)}<>{} then
 print(`The set of edges`, G, `is not good `):
 RETURN(FAIL):
fi:



for i from 1 to n do
lu:=0:
for j from 1 to n do

 if i<>j and member({i,j},G) then
 M[i,j]:=-a[i,j]:
  lu:=lu+a[i,j]:
  else
  M[i,j]:=0:
 fi:
od:
M[i,i]:=lu:
od:

[seq([seq(M[i,j],j=1..n)],i=1..n)]:
end:


#Mat1(G,n,V0,a): Given a graph G on {1,2, ...,n} described as a set of edges, and a subset of {1, .., n},finds the principal minor of Mat11(G,n,a) (q.v.)
#with the members of V0 removed. Try
#Mat1({{1,2},{2,3},{3,4},{4,1}},{4},a);
Mat1:=proc(G,n,V0,a) local i,j,gu,lu,M:

M:=Mat11(G,n,a):

gu:=[]:

for i from 1 to n do
if not member(i,V0) then
  lu:=[]:
 for j from 1 to n do
    if not member(j,V0) then
     lu:=[op(lu),M[i][j]]:
    fi:
 od:
gu:=[op(gu),lu]:
fi:
od:
gu:
end:

#IndToEdge(M1,a): inputs an indexed expression of the form a[i,j] outputs the edge {i,j}. Try:
#IndToEdge(a[4,7],a);

IndToEdge:=proc(M1,a):

if not type(M1, indexed) then
  RETURN(FAIL):
fi:
if op(0,M1)<>a then
RETURN(FAIL):
fi:
if nops(M1)<>2 then
 RETURN(FAIL):
fi:

if op(1,M1)=op(2,M1) then
  RETURN(FAIL):
fi:
{op(1,M1),op(2,M1)}:
end:

#MonoToGraph(M,a): inputs a monomial in a[i,j] and outputs the set of {i,j}-s that show up. Try:
#MonoToGraph(a[1,2]*a[2,4],a);
MonoToGraph:=proc(M,a) local i,M1,lu,gu:

if type(M,indexed) then
 RETURN({IndToEdge(M,a)}):
fi:
if not type(M,`*`) then
  RETURN(FAIL):
fi:

gu:={}:

for i from 1 to nops(M) do
M1:=op(i,M):

lu:=IndToEdge(M1,a):
 if lu=FAIL then
   RETURN(FAIL):
 else
 gu:=gu union {lu}:
 fi:

od:

gu:

end:



 
#PolyToSG(P,a): inputs a polynnomial in a[i,j] and outputs the set of graphs corresponding to each monomial.
#PolyToSG(a[1,2]*a[2,4]+a[1,3]*a[2,5],a);
PolyToSG:=proc(P,a) local M,gu,i,lu:

if not type(P,`+`) then
 RETURN({MonoToGraph(P,a)}):
fi:
gu:={}:

for i from 1 to nops(P) do
 M:=op(i,P):

 lu:=MonoToGraph(M,a):

 if lu=FAIL then
  RETURN(FAIL):
 else
 gu:=gu union {lu}:
 fi:

od:

gu:

end:


#SpF(G,n,V0): The set of spanning forests of the graph G with set of vertices {1, ..., n} where the members of V0 belong to different components.
#Try:
#SpF({{1,2},{1,3},{2,3}},3,{1});
SpF:=proc(G,n,V0) local a:
PolyToSG(det(Mat1(G,n,V0,a)),a):
end:

#Kn(n): the complete graph on {1, ..., n}. Try:
#Kn(5);
Kn:=proc(n) local i,j:
{seq(seq({i,j},j=i+1..n),i=1..n)}:
end:


#Mat11N(G,n): Given a graph G on {1,2, ...,n} described as a set of edges, finds 
#the matrix whose determinant would give the Laplacian matrix, whose determinant of any n-1 by n-1 principal minor would
#give the number of the spanning trees. Try
#Mat11N({{1,2},{2,3},{3,4},{4,1}},4);
Mat11N:=proc(G,n) local M,i,j,lu:

if G minus {seq(seq({i,j},i=1..j-1),j=1..n)}<>{} then
 print(`The set of edges`, G, `is not good `):
 RETURN(FAIL):
fi:



for i from 1 to n do
lu:=0:
for j from 1 to n do

 if i<>j and member({i,j},G) then
 M[i,j]:=-1:
  lu:=lu+1:
  else
  M[i,j]:=0:
 fi:
od:
M[i,i]:=lu:
od:

[seq([seq(M[i,j],j=1..n)],i=1..n)]:
end:


#Mat1N(G,n,V0): Given a graph G on {1,2, ...,n} described as a set of edges, and a subset of {1, .., n},finds the principal minor of Mat11N(G,n) (q.v.)
#with the members of V0 removed. Try
#Mat1N({{1,2},{2,3},{3,4},{4,1}},{4});
Mat1N:=proc(G,n,V0) local i,j,gu,lu,M:

M:=Mat11N(G,n):

gu:=[]:

for i from 1 to n do
if not member(i,V0) then
  lu:=[]:
 for j from 1 to n do
    if not member(j,V0) then
     lu:=[op(lu),M[i][j]]:
    fi:
 od:
gu:=[op(gu),lu]:
fi:
od:
gu:
end:

#SpFn(G,n,V0): The NUMBER of spanning forests of the graph G with set of vertices {1, ..., n} where the members of V0 belong to different components.
#Try:
#SpFn({{1,2},{1,3},{2,3}},3,{1});
SpFn:=proc(G,n,V0)
det(Mat1N(G,n,V0)):
end:


#GFGridKN(k,t): input an integer k and a symbol t, output the generating function f(t) whose coefficient of t^n is the number of spanning trees in k by n grid graph. 
GFGridKN:=proc(k,t) local d,L,n:
d:=3/2*2^(k-1):
L:=[0,1,seq(SpFn(GridMN(k,n)[1],k*n,{1}),n=2..d)]:
CtoR(GuessSymRec(L),t):
end:


#DenomSFKN(k,t): input an integer k and a symbol t, output the denominator of the generating function f(t) whose coefficient of t^n is the number of spanning forests in k by n grid graph with two components such that 1 and k*n are in different components / the square of the denominator of GFGridKN(k,t).

DenomSFKN:=proc(k,t) local L,f,g,gu:
L:=[0,k-1,seq(SpFn(GridMN(k,n)[1],k*n,{1,k*n}),n=2..3^k+2*k)]:
g:=add(t^(n-1)*L[n],n=1..nops(L)):
f:=denom(GFGridPreCal(k,t)):
g:=expand(g*f^2):
L:=[seq(coeff(g,t,i),i=2^k+k..3^k+k)]:
gu:=simplify(CtoR(GuessSymRec(L),t)):
denom(gu):
end:


#GridH(m,n,t): input two integers m and n, and a symbol t, output the generating function f(t) whose coefficient of t^i is the number of spanning trees in the m by n grid graph which have i horizontal edges. 
GridH:=proc(m,n,t) local G,S,f,s,co,edge:
G:=EdgeToGraph(GridMN(m,n)[1]):
S:=GST(G):
f:=0:
for s in S do
  co:=0:
  for edge in s[2] do
    if edge[2]-edge[1]<>1 then
      co:=co+1:
    fi:
  od:
  f:=f+t^co:
od:
f:
end:

#GridV(m,n,t): input two integers m and n, and a symbol t, output the generating function f(t) whose coefficient of t^i is the number of spanning trees in the m by n grid graph which have i vertical edges. 
GridV:=proc(m,n,t) local G,S,f,s,co,edge:
G:=EdgeToGraph(GridMN(m,n)[1]):
S:=GST(G):
f:=0:
for s in S do
  co:=0:
  for edge in s[2] do
    if edge[2]-edge[1]=1 then
      co:=co+1:
    fi:
  od:
  f:=f+t^co:
od:
f:
end:


#GridHV(m,n,x,y): input two integers m and n, and a symbol t, output the generating function f(x,y) whose coefficient of x^i*y^j is the number of spanning trees in the m by n grid graph which have i horizontal edges and j vertical edges. 
GridHV:=proc(m,n,x,y) local G,S,f,s,co,edge:
G:=EdgeToGraph(GridMN(m,n)[1]):
S:=GST(G):
f:=0:
for s in S do
  co:=0:
  for edge in s[2] do
    if edge[2]-edge[1]=1 then
      co:=co+1:
    fi:
  od:
  f:=f+x^(m*n-1-co)*y^co:
od:
f:
end:


#GridHFactorial(i,k,n): input three integers i,k and n, output a list of i-th factorial moment of GridH(k,j) for j from 1 to n.
GridHFactorial:=proc(i,k,n) local L,j,f,t,num:
L:=[]:
for j from 1 to n do
  f:=GridH(k,j,t):
  num:=subs(t=1,f):
  L:=[op(L), subs(t=1,diff(f, t$i))/num]:
od:
L:
end:
#for example GridHFactorial(1,2,9) = [0, 3/2, 44/15, 61/14, 1208/209, 2809/390, 25108/2911, 13645/1358, 155024/13515].


#GridVFactorial(i,k,n): input three integers i,k and n, output a list of i-th factorial moment of GridV(k,j) for j from 1 to n.
GridVFactorial:=proc(i,k,n) local L,j,f,t,num:
L:=[]:
for j from 1 to n do
  f:=GridV(k,j,t):
  num:=subs(t=1,f):
  L:=[op(L), subs(t=1,diff(f, t$i))/num]:
od:
L:
end:
#for example GridVFactorial(1,2,9) = [1, 3/2, 31/15, 37/14, 673/209, 1481/390, 12735/2911, 6725/1358, 74731/13515]




#ListSTPreCal(k): input an integer 2<=k<=7, and output a list of pre-calculated numbers which is long enough to get a recurrence via GuessSymRec for k by n grid graph. 
ListSTPreCal:=proc(k)
if k<2 then
  return FAIL:
fi:
if k=2 then
  return [0, 1, 4, 15]:
elif k=3 then
  return [0, 1, 15, 192, 2415, 30305, 380160]:
elif k=4 then
  return [0, 1, 56, 2415, 100352, 4140081, 170537640, 7022359583, 289143013376, 11905151192865, 490179860527896, 20182531537581071, 830989874753525760]:
elif k=5 then
  return [0, 1, 209, 30305, 4140081, 557568000, 74795194705, 10021992194369, 1342421467113969, 179796299139278305, 24080189412483072000, 3225041354570508955681, 431926215138756947267505, 57847355494807961811035009, 7747424602888405489208931601, 1037602902862756514154816000000, 138964858389586339640223412108401, 18611389483734199394023624777573409, 2492600085599977923424220468405177105, 333830807688353225138019865387722924481, 44709541971379003103897461691112357888000, 5987892960038182781131697625354150226327105, 801951004627869433685025226859351146717402769, 107404293649401297954327034703922488508540561569, 14384522530358739351890623742584897464468359377905]:
elif k=6 then
  return [0, 1, 780, 380160, 170537640, 74795194705, 32565539635200, 14143261515284447, 6136973985625588560, 2662079368040434932480, 1154617875754582889149500, 500769437567956298239402223, 217185579535490113365186969600, 94193702839904633186530210863025, 40851869157273984726590135085017940, 17717469746416280095776019395706656000, 7684070867169415429692559499446691755680, 3332583081296808509759455619848802299528513, 1445341907485491645328460310146924377335398400, 626845049313054375044367343971643549398400207439, 271862811296852944176805652529210910158678393501000, 117906950273496738441417982275113800569185953508401920, 51136265564260810934348514677358578182539093981985579140, 22177807578786871771975685968088359261519721558786693144351, 9618519137859358557484551263002772100493327681623268157030400, 4171553480859669166764609074130946419690927672433431476417301025, 1809203495268378561427549408448430750759030036616841097919332985500, 784651881439272378819594430425368142664597935508249012776181518286080, 340303662167183966505537475986213779843792970000024727053163782358717080, 147589759514663135524738411416132378729266135389120483181714254960506345457, 64009705258157188975303879517278636748220260668937846227901842986863001600000, 27761020688084252861438737654487965787656058577609755240697147365225357658880511, 12039959667616423898273302782896010357082899393078794364724702223338929390272343360, 5221732674261867556341716244029432951359105841004959005308354662058227559497202420480, 2264666400402646369686180160329283928881123430191123975339161311624342834915309308782540, 982186225348592455606198909858639570930247856454715089145947291692834447312970875836923535, 425974342663890239552204750871916025730482119948962370321564289059227721374521817656878694400, 184745149061251184750862045211510331356516624797499136987274537982693388581135194499647643913041, 80124004390082193434085642426318991468753999560516408431413166233551810878934366969111416133330660, 34749795121134381681821839317664554206982305871476178101095508198643466384328317097198320510995744000, 15070992396759017420870386902511562713661550891173802019350772669184640978192637384627777779333515514000, 6536292114281491680639421330783959347585757245049468192958544444845571990082859048577271547340334575596897, 2834791066075112686863202316024234550237660565354148599781112233550332088730359206236952468623286542879948800, 1229449395436428336597016386379198373474847531314379432973439673094182365634137324929779722563621274783913620207, 533212424022415807911061451425420454123254290970263295391190628314462857450089007654039033611294842931513917846840, 231254324242385438301103992134638746932293869484870852254306060133513286690122058148739547644691024996747106183648000, 100295042034793516035937040154695457921537232804819506882483180218880229602947954326808953644968025427360497478207467540, 43497977777132166820598001749673834018899288729543261941443597022506674629218089699381370066419505795784783072595426620223, 18865080788774198326107739046655419206613363795216314368754269542747051403827127437331851904725251589964297145427834594918400]:
elif k=7 then
  return [0, 1, 2911, 4768673, 7022359583, 10021992194369, 14143261515284447, 19872369301840986112, 27873182693625548898079, 39067130344394503972142977, 54740416599810921320592441119, 76692291658239649098972455530913, 107441842254735898225957962027174559, 150517199699838971875005120330439121217, 210860422397100784567572728149075575177216, 295394542455170446994290978914210208936856257, 413817852397080403148772215737537261617562185247, 579716568874689646501757437321208369959290374622113, 812123531133033168142883595132576030926156951442592543, 1137701760634799958175744040984009442816388920028037144833, 1593803394109789292797576569479470546332292105269153758133407, 2232754928999722159130672264241999812738025912226830847970115584, 3127860409030804832077623110864882699311021626916094040318841846367, 4381811268884140282919155969069868192446292013819070322586801748900801, 6138467661942008213109992731398215374098731636135137023549118053901775903, 8599362890097379978412302513147331911895856972299296233430895989256768574881, 12046824414846834643642577662840027453702639483076290726194886922508852622617567, 16876364019516071743543634751182092459347416721241149154591948423807265185066276481, 23642053100541343365851738096759129095486700260436370416654508151840816551118573142016, 33120088791277080288248658818564781485527366292764043666969724060788126481246419278707073, 46397843574339855146335797833427214880938469066803965077370775337455327420387155128809566943, 64998614644724866766508458788115362403573819875257030517535696736143493773026573423827861342625, 91056384958058708845315575858125224492307112362287733149646139876719145156313268462028201438370847, 127560645514456012592885670138342569757279573008931567530944258478826887839655731130009409050484948673, 178699366239463426503803563145487841371388433855926407134017619501096060963675707863065876050636779025759, 250339462971473686581456694581019284734049621732814282114622752443365734598445247048637763780896611307880448, 350699882375990465592853080807886127871182090166852581826909424266851361056701025879490751403298734479454286239, 491294524797102255984197885637849839375660034194337934805829163959236629983910453715119873414903447529762781107713, 688253182351605652390729628944931238096615656187114423584459860234129114665168108703544634427919285877185907045241631, 964172037562885144989408438839753915806603893844745195745521345487653144871487603546571604871821228848892593621408362401, 1350706021934888594174070772657494242101733732831966475192205527373104234798197890648677498633778642054745858119516894048031, 1892200444126839502369865739940932382339626691118801453153015194101882010591016530300302593029141117148800063187299837457119169, 2650778528124755789288820034779453562299580104941557155236912513786973117675464236754671005700456961701004430622574907625008267264, 3713468531823382110418816226747511016372584877763257842835567157221168974380677811397538755437254853176286283497796188834896913681473, 5202188108335809492924453363090508717157315770374277809447047329436440329876773545322504214096434634119277832968586607719607319265171871, 7287731370978439032401734992835593513361351429437810106705504956477184381447830117741363678632798413462015901811663978584105297262237941153, 10209363335101237594613046494853287238832751340799193436587195994508821007274498623302212340523619618014295294366752192464411425363482742017823, 14302269719104035231972666468785479862335151812644263288523486091664094468203988177674933827737670014565642118740905313811339462851122569047898753, 20036011297071916885178889933689101704898531192182570745608908110011900890075467858574307343581261575115769727509230850972819925883600797414170978847, 28068394498264417411011910034206677958264952968747763691723760486541879605143392067116198233471629002325891845417650571822474563683775015796544603160576, 39320938585482590171895141483187054341524215989813588752188968459992042985282101558904683119171708366581179941147438610348114619194604988811087613788382431, 55084597422873526569853327801631720594970017549558838720123499800909166154017336054994700433472503773577751100913083201901976116726999225748256191498712813121, 77167864816948247618828046724147228638846326171153085138630315180735963767024319811065277614710774937872227084673176451314696528391077057906594305036106339146783, 108104254891660035176093278405208765331387354887090890393001310940009233174714478469471007160801999302986551052697312345607382655092206569819958411816824087112630177, 151442960789480206585596889009332184133851296196189867980689337832217755447144059865408654449914394804221439710364821082925661388990146084968837926385192789109315884127, 212156037666315881596502384102285213886165976087926711964481142099538207272734951566151124500164678767226402152725226355621591607554307082662257971340384056388919027973377, 297208824257203894652197995394073346166223388431070935961315576717760841596734303363042881269289203924333340711415737006869988501710489199805624496431781064196876956118024192, 416359044918070684486185068744862405064062422323116139736390112097411774079553065755764802550323116078126356692096279492978673251665771820870571144794130804290496630733499532033, 583276269533192169606911766304778963001274713950533556392449259616900554711806957929645935605662071561078121662755727795393649090361127535809450560883991168614381106760622756061791, 817110161897653301238127799538180679610507556638874723925732819226660731388243385563142396252875709524044706872968412761729724266599590522867160568132663188801954993439909436956092321, 1144687434670980608651359694101074570897333028012093056977604188233679129693124911129804955370701950907586217991245330866453367015533564911110932609776584599255065457951676161396774579231, 1603589557680415429969370346604251701380041918528599174160062629007738816636495118894127971695577456692920044282323605625438469003369558155068258167074643760247769226980083831644617563903041, 2246464311229904192729569202502413603917585207680755357785759244120173678746767482731586672339085462510468843289679091146415037023766513208301547797027108313426591457358282415499319106684017375, 3147065829568966006503653886053659407965504392044320065633527832229201415782236052875629798299179201327934515777587896191132540703396670501471815429052688817799359109473254096975255911666909184000, 4408716081591479101611260735594046173972481638015450640631745790050084239814504589543342769673314128103083788810413904095443691213253333286093050665712778615742790158801057183776345245557781098746911, 6176158536456626868044666080895207162780019545465534613224129326438379294103447771976362524616833360753365733842772502917795376236292634299658238951689021213722939350799818762303422872399786601286820993, 8652163931970943584201968809926551832936906486804238225017148766987851512941032726213726080520826366848222915363354508742797215323929649050156311822556895843809523891320771195827782951236939019947125755679, 12120793250985327203405298739760434898648763299743343868180360930577678435753267593407964090879198530068652466726819263765738497041951644247082088030093660403184275629909690453492317608178580059197785255757217, 16979986762648504414484329219329986897832697235195099671160573401587460480433432280407307852470211871601812528999607058959254920850342427960386751939570643081843811986109585160965669040256146863656995854275959711, 23787217922909479911946580018071893278490341729325259344668951668858818040176741826208337721812734853560841170929108629768825137721572799444639191200870108669708206234534270543263194212720595108996794946613073618497, 33323449801307659435567558291694985963816765350726998476204987069543540195214946147722116648041129925123274727929252631391693622069203616709498756392994072188617125746971278751248227564977663690314309252707468751929344, 46682731467760018640915929320638420350195172223836062219052339478021068155199344850038406251846682075691593475383745568698231882487996139722123425732579846586005275707735971662123315439648583270395879465282110516264607169, 65397713330552994769344768446296375409623298987774200790112346322286299999731399119552866330797147930796101723491888733204712570197435427114202906344637534073734999221897280713126951967187404028851547025072415337495302720799, 91615481236758185046950464052454891019893176170990292349655697198242931911047934983891053028116737947448086655754200462443083596725824891877422586817296392771681135478026965765135354998896112895360647485302350911688423883229089, 128343881991994681997579434340210453796610264864654374299136112218198426176932067964996783465932493962103529851665838611503763212688640647249165454160482776979262940822804436392345868911856577373166985666532412712222467957699556127, 179796600120527011147444756803260253945812342412193402988094554145668937074888897286410336271998892179779983149440419210209221419458965031567027219784668445970667476211352823891019126976021891460494554135658470058138107256033463176193, 251876574973141654851086836626832851595730870269923379472165051335584032299636128049430926521697658547134208384155408107913149016489167711728141384291719043084449692285649620267966085967564888317794252363330463829865317774521646041892767, 352853218457258397732545569059862631289314631996679693254843631948801068076593586144744872894769271014738117711143571927222846636238868046488968133637253261720870540720188022549768697233284021440542902969011289650243426554718380510506909696, 494311127539042077566691082037509296506139449067611336571649226688969453215703501679945108982244072940993884282995292486710279361236078320980635184836899738145626403727189966375877519639785864358069880569186859819247040697086037696430821064543, 692479133043579679572773056676716122384544562800368136097584749713654462767547568576066122572788495010401552347748886435057941388957300163601827086957775528835014816197372714246565463848087415032905470289816017525253742572141559023161537456216257, 970092160555121861133119191758709805251265641252978813687179797740633080971523076598456013215017798566921642578001668116072605987355329411610480210377559628367407738531529228518562710866575082548297440417944652703853269628637353430257228706637130783, 1358999506359536133403098840650186967708948312218961172018932220289467684104625433529348604949524618948151365701590305651319010490641404226122818751316489529369048477344342473391533376986925059593418770698813549829563112050152712106363518652566468002209, 1903818764217836247480269086604743470570701380428075528146733846863835259195429864449489101600621939703473600431093945405033226430679030877286298163471117821998992890547234772861461669713979976366461718598426100001867649269225987440363127739960396540724447, 2667054601584988852735617891828986131450748612227626415183582330447714568118433487589511141750409576124828379636597971489692121454372200745022371733866565652136362461449978568506145199115423874608990681420445206680868296674252636963458210028314876312365937537, 3736269639488524689184151881203318401672828861575282369153675329688157863940525147419498555112942207932616095228946349198473591199060284821844231965479808766981311849525064248959433869077909740613666947016126108314924483062044291000392695452577083321925544968192, 5234130118921327053890053720778922794915451510534676653134022644939892458798501996429636025626795539692787498127061800224321523619775580876198648542646454088579340094719922187946436015371416287982905027143299324108487693789077837647445822557001932641163154553986177, 7332478847953213351390940767321445442685323439993603499656421880564697405307510915273522491858668453446737080893044833479352758278143735322420905104741918440589182070314217555247530584273474799946088479341434794238007967621025691294483945040261678652916404552334383583, 10272049955602071518545004659169976809928189302064136439787413738134270510735867772210848927638059205445826914989859637931646904010319555350906636468915999828831521547080717255732929701363518931215320522620028042576101544696686946707286666616228004571913543453866187857313, 14390087237665603070330036916126385858121282258931613638495581347841159295262975160276616182971337305780452069896025995580367876099895706606686758496319574083924181821880386187509754917594344067417527505296248614670132841473311360813413952678102259671959861245773594346759199, 20159034623336709237790883008434733142100984612986977220681613091041338187609136106920885309566166661501013977445714683489673102859660625802494893866136153702460239140956667968767866652836857820155621793482251023133557206797994357328606345457534983353112215489462761752109605313, 28240737546134105960669963159055315765230316301090748532288812637819345665161865319953266784313199114335301645448468915391209821499627659414517415325552075442355501985922105196320364989176933284635789628616621809371138276321458004218965887305042621692745293798934223556531781742687, 39562373499094693626252860509013904751983254336472822509058487249374075164539226500904303123435213864193962294193360371507723919844538708670101565499247749849383763355725681119256586000203952225887514527779904360880563805043716589029303240671531120926445690855642595284952001644593152, 55422822945930067364338533912275237604302577278185868984336704906925604275628301165869704801617925286065651041274876488868323950631304368866169067529194433106108130518286636855414955422449343869701580477224027337148079022270464688221338707540912556639054711292188980411718798981802878623, 77641684045225494425449069632950938851417197256483701220345981068330113640451275110607673197884016807291567684953155795111756975879462337543149641700796617981632805377865114089416184054857781880736719533506356127056562369700615289000536444306245609375937557952162257842380215036850501626625, 108768026977978061267420184525956880297279465188409317104419877773363025438859934959160224305977296095113144962794550096380340558429327628814728012403489475899284519279878461049820437562736665305237191865890802550335308018598320956677633401703383607221660944017038697139496677953048941178794783, 152372837325257221908356767120523766414532702793676154020284762407267727142880423059489778405165996817154433785709250793440687629584637734521349380004181103295746613644675531272180213634653591418428219065216618363086645346487596944030200362122356270756287441404211431340321171365706583687305311137, 213458699211764453680856689453367270810052327548199294455849932293932589552821737613449290383513542640020041034223287136838096287713817304278262411192783303040928093975675966868274528604202017886734717741204886878805306709161932122116514067032987113750188797165208223820575936807405241852994996441119]:
fi:
end:




#GFGridPreCal(k,t): input an integer 1<=k<=7 and output the pre-calculated generating function f(t) whose coefficient of t^n is the number of spanning trees in k by n grid graph. 
GFGridPreCal:=proc(k,t) 
if k<1 then
  return FAIL:
fi:
 
if k=1 then
  return t/(1-t):
elif k=2 then
  return t/(t^2-4*t+1):
elif k=3 then
  return (-t^3+t)/(t^4-15*t^3+32*t^2-15*t+1):
elif k=4 then
  return (t^7-49*t^5+112*t^4-49*t^3+t)/(t^8-56*t^7+672*t^6-2632*t^5+4094*t^4-2632*t^3+672*t^2-56*t+1):
elif k=5 then
  return (-t^15+1440*t^13-26752*t^12+185889*t^11-574750*t^10+708928*t^9-708928*t^7+574750*t^6-185889*t^5+26752*t^4-1440*t^3+t)/(t^16-209*t^15+11936*t^14-274208*t^13+3112032*t^12-19456019*t^11+70651107*t^10-152325888*t^9+196664896*t^8-152325888*t^7+70651107*t^6-19456019*t^5+3112032*t^4-274208*t^3+11936*t^2-209*t+1):
elif k=6 then
  return (t^31-33359*t^29+3642600*t^28-173371343*t^27+4540320720*t^26-70164186331*t^25+634164906960*t^24-2844883304348*t^23-1842793012320*t^22+104844096982372*t^21-678752492380560*t^20+2471590551535210*t^19-5926092273213840*t^18+9869538714631398*t^17-11674018886109840*t^16+9869538714631398*t^15-5926092273213840*t^14+2471590551535210*t^13-678752492380560*t^12+104844096982372*t^11-1842793012320*t^10-2844883304348*t^9+634164906960*t^8-70164186331*t^7+4540320720*t^6-173371343*t^5+3642600*t^4-33359*t^3+t)/(t^32-780*t^31+194881*t^30-22377420*t^29+1419219792*t^28-55284715980*t^27+1410775106597*t^26-24574215822780*t^25+300429297446885*t^24-2629946465331120*t^23+16741727755133760*t^22-78475174345180080*t^21+273689714665707178*t^20-716370537293731320*t^19+1417056251105102122*t^18-2129255507292156360*t^17+2437932520099475424*t^16-2129255507292156360*t^15+1417056251105102122*t^14-716370537293731320*t^13+273689714665707178*t^12-78475174345180080*t^11+16741727755133760*t^10-2629946465331120*t^9+300429297446885*t^8-24574215822780*t^7+1410775106597*t^6-55284715980*t^5+1419219792*t^4-22377420*t^3+194881*t^2-780*t+1):
elif k=7 then
  return (-t^47-142*t^46+661245*t^45-279917500*t^44+53184503243*t^43-5570891154842*t^42+341638600598298*t^41-11886702497030032*t^40+164458937576610742*t^39+4371158470492451828*t^38-288737344956855301342*t^37+7736513993329973661368*t^36-131582338768322853956994*t^35+1573202877300834187134466*t^34-13805721749199518460916737*t^33+90975567796174070740787232*t^32-455915282590547643587452175*t^31+1747901867578637315747826286*t^30-5126323837327170557921412877*t^29+11416779122947828869806142972*t^28-18924703166237080216745900796*t^27+22194247945745188489023284104*t^26-15563815847174688069871470516*t^25+15563815847174688069871470516*t^23-22194247945745188489023284104*t^22+18924703166237080216745900796*t^21-11416779122947828869806142972*t^20+5126323837327170557921412877*t^19-1747901867578637315747826286*t^18+455915282590547643587452175*t^17-90975567796174070740787232*t^16+13805721749199518460916737*t^15-1573202877300834187134466*t^14+131582338768322853956994*t^13-7736513993329973661368*t^12+288737344956855301342*t^11-4371158470492451828*t^10-164458937576610742*t^9+11886702497030032*t^8-341638600598298*t^7+5570891154842*t^6-53184503243*t^5+279917500*t^4-661245*t^3+142*t^2+t)/(t^48-2769*t^47+2630641*t^46-1195782497*t^45+305993127089*t^44-48551559344145*t^43+5083730101530753*t^42-366971376492201338*t^41+18871718211768417242*t^40-709234610141846974874*t^39+19874722637854592209338*t^38-422023241997789381263002*t^37+6880098547452856483997402*t^36-87057778313447181201990522*t^35+862879164715733847737203343*t^34-6750900711491569851736413311*t^33+41958615314622858303912597215*t^32-208258356862493902206466194607*t^31+828959040281722890327985220255*t^30-2654944041424536277948746010303*t^29+6859440538554030239641036025103*t^28-14324708604336971207868317957868*t^27+24214587194571650834572683444012*t^26-33166490975387358866518005011884*t^25+36830850383375837481096026357868*t^24-33166490975387358866518005011884*t^23+24214587194571650834572683444012*t^22-14324708604336971207868317957868*t^21+6859440538554030239641036025103*t^20-2654944041424536277948746010303*t^19+828959040281722890327985220255*t^18-208258356862493902206466194607*t^17+41958615314622858303912597215*t^16-6750900711491569851736413311*t^15+862879164715733847737203343*t^14-87057778313447181201990522*t^13+6880098547452856483997402*t^12-422023241997789381263002*t^11+19874722637854592209338*t^10-709234610141846974874*t^9+18871718211768417242*t^8-366971376492201338*t^7+5083730101530753*t^6-48551559344145*t^5+305993127089*t^4-1195782497*t^3+2630641*t^2-2769*t+1):
fi:

end:


#DenomSFPreCal(k,t): input an integer 1<=k<=4 and a symbol t, output the pre-calculated denominator of the generating function f(t) whose coefficient of t^n is the number of spanning forests in k by n grid graph with two components such that 1 and k*n are in different components/the square of the denominator of GFGridKN(k,t).
DenomSFPreCal:=proc(k,t)
if k<2 then
  return FAIL:
fi:
if k=2 then
  return t-1:
elif k=3 then
  return t^4-8*t^3+17*t^2-8*t+1:
elif k=4 then
  return t^12-46*t^11+770*t^10-6062*t^9+24579*t^8-55388*t^7+72324*t^6-55388*t^5+24579*t^4-6062*t^3+770*t^2-46*t+1:
fi:
end:

#GuessRat(L,degT,degB): input a list L and the degrees of top and button, try to find a rational function R(x)=P(x)/Q(x) with deg(P)<=degT, deg(Q)<=degB, such that R(i)=L(i).
GuessRat:=proc(L, degT, degB) local n,a,b,var,eq,i,j,k,sol,F,x:
n:=nops(L):
if n<degT+degB+3 then
  return FAIL:
fi:
var:={}:
eq:={}:
for i from 0 to degT do
  var:=var union {a[i]}:
od:
for j from 0 to degB do
  var:=var union {b[j]}:
od:
for k from 1 to n do
  eq:=eq union {add(a[i]*k^i,i=0..degT)/add(b[j]*k^j, j=0..degB)-L[k]}:
od:
sol:=solve(eq, var):
F:=add(a[i]*x^i,i=0..degT)/add(b[j]*x^j,j=0..degB):
if {sol} <> {} then
  return simplify(subs(sol,F)):
else
  return FAIL:
fi:
end:



#JCgND(a,n,G): Like JCg(a,3,{{1,2},{2,3}}) but with numerator and denominator separate
JCgND:=proc(a,n,G) local gu,i1,gu1,gu2:
gu:=Lg(a,n,G):
gu1:=[seq([op(1..n-1,gu[i1])],i1=1..n-1)]:
gu2:=[seq([op(2..n-1,gu1[i1])],i1=2..n-1)]:
[normal(-det(gu1)), normal(det(gu2)) ]:
end:

#RmnND(m,n): the numerator and denominator formula for the joint resistance between opposite vertices of a m by n grid via the general formula.
RmnND:=proc(m,n) local gu:
gu:=JCnumND(m*n, UCmn(m,n)):
[gu[2],gu[1]]:
end:
#seq(RmnND(k,k),k=2..11) = 1, 3/2, 13/7, 47/22, 1171/495, 6385/2494, 982871/360161, 441083/153254, 427854195/142065451, 4769986941/1522703822


#DiagMatrix(a,b,c,p,q,n): input symbols a,b,c, two positive integers p and q, and a positive integer n > max(p,q) output a matrix as a list of lists where the diagonal element is a, in the first row after a it’s b[1],…,b[p], in the first column under a it’s c[1],…,c[q].
DiagMatrix:=proc(a,b,c,p,q,n) local M,L,i:
if n<= max(p,q) then
  return FAIL:
fi:
L:=[0$(n-1-q), seq(c[q-i+1],i=1..q), a, seq(b[i], i=1..p), 0$(n-1-p)]:
M:=[0$n]:
for i from 1 to n do
  M[i]:=[op(n+1-i..2*n-i,L)]:
od:
M:
end:

#ExpandMatrix(M): input a “general diagonal” matrix, do an expansion along its first row or column. (choose row if the number of nonzero elements in the first row is no more than that in the first column). output a set of [multiplicity, Matrix]. 
ExpandMatrix:=proc(M) local n,r1, r2,i,S,T,k,j:
n:=nops(M):
S:={}:
r1:=n+1:
r2:=n+1:
for i from 1 to n do
  if M[1][i] = 0 then
    r1:=i:
    break:
  fi:
od:
for i from 1 to n do
  if M[i][1] = 0 then
    r2:=i:
    break:
  fi:
od:
if r1<=r2 then
  for i from 1 to r1-1 do
    T:=[op(2..n,M)]:
    T:=[seq([seq(T[j][k],k=1..i-1), seq(T[j][k],k=i+1..n)],j=1..n-1)]:
    S:=S union {[M[1][i]*(-1)^(1+i), T]}:
  od:
else
  for i from 1 to r2-1 do
    T:=[op(1..i-1,M), op(i+1..n,M)]:
    T:=[seq([op(2..n,T[j])],j=1..n-1)]:
    S:=S union {[M[i][1]*(-1)^(1+i), T]}:
  od:
fi:
S:
end:


#ExpandMatrixR(M): input a “general diagonal” matrix, do an expansion along its first row. output a set of [multiplicity, Matrix]. 
ExpandMatrixR:=proc(M) local n,i,S,T,k,j:
n:=nops(M):
S:={}:
for i from 1 to n do
  if M[1][i]=0 then
    break:
  fi:
  T:=[op(2..n,M)]:
  T:=[seq([seq(T[j][k],k=1..i-1), seq(T[j][k],k=i+1..n)],j=1..n-1)]:
  S:=S union {[M[1][i]*(-1)^(1+i), T]}:
od:
S:
end:



#IsDiagonal(a,b,c,p,q,M): determine whether a matrix M is a (a,b,c,p,q)-type diagonal matrix by only examine its first row and first column.
IsDiagonal:=proc(a,b,c,p,q,M) local i,n,M1,M2,L1,L2:
n:=nops(M):
if n<=max(p,q) then
  return false:
fi:
M1:=[op(1..p+1,M[1])]:
M2:=[seq(M[i][1], i=1..q+1)]:
L1:=[a, seq(b[i],i=1..p)]:
L2:=[a, seq(c[i],i=1..q)]:
if M1=L1 and M2=L2 then
  return true:
else
  return false:
fi:
end:


#LSUnion(A, B): input two sets A and B which contain lists like [multiplicity, a data structure], return a set containing the same kind of list for which the number is the sum of the multiplicity in A and B.
LSUnion:=proc(A,B) local A1, B1,C,C1,a,b,c,co:
A1:={seq(a[2], a in A)}:
B1:={seq(b[2], b in B)}:
C1:=A1 union B1:
C:={}:
for c in C1 do
  co:=0:
  for a in A do
    if a[2]=c then
      co:=co+a[1]:
    fi:
  od:
  for b in B do
    if b[2]=c then
      co:=co+b[1]:
    fi:
  od:
  C:=C union {[co, c]}:
od:
C:
end:




#MatrixRecurrence(a,b,c,p,q): input symbols a,b,c and positive integers p,q, output a list L which represents the recurrence relation of the determinant of this matrix. For example, if we use D(n) to denote this kind of matrix of dimension n (large enough), then [a, -b*c] means D(n) = a*D(n-1) - b*c*D(n-2). 
MatrixRecurrence:=proc(a,b,c,p,q) local M, n,S,L,s,dim,C,cc,SC,i,gu,SD:
n:=4*p*q+5:
L:=[0$n]:
M:=DiagMatrix(a,b,c,p,q,n):
S:=ExpandMatrixR(M):
SC:=S:
for i from 1 to n do
  if S={} then
    break:
  else
    for s in S do
      if IsDiagonal(a,b,c,p,q,s[2]) then
        dim:=nops(s[2]):
        L[n-dim]:=L[n-dim]+s[1]:
        SC:=SC minus {s}:
      else
        C:=ExpandMatrixR(s[2]):
        C:={seq([cc[1]*s[1],cc[2]],cc in C)}:
        SC:=SC minus {s}:
        SC:=LSUnion(SC,C):
      fi:
    od:
    S:=SC:
  fi:
od:
gu:=0:
for i from n to 1 by -1 do
  if L[i]<>0 then
    gu:=i:
    break:
  fi:
od:
[op(1..i,L)]:    
end: