##########################################################################
## GFMatrix.txt Save this file as GFMatrix.txt                          ##
## to use it, stay in the same directory, get into Maple                ##
## (by typing: maple <Enter> ) and then type:                           ##
## read GFMatrix.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. 27, 2018
with(combinat):
with(linalg):
with(plots):
with(LinearAlgebra):
with(numtheory):

print(`First Written: Oct 2018`):
print(`Version of Oct. 27, 2018:`):
print():
print(`This is GFMatrix.txt, one of the  Maple packages`):
print(`accompanying Yukun Yao's and Doron Zeilberger's  article: `):
print(`Untying The Gordian Knot via Experimental Mathematics`):


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 DiagMatrix(a,b,c,p,q,n) , DiagMatrixL(L) , ExpandMatrixL(L,L1) , ExpandMatrixPer(L,L1) , ChildrenMatrixL(L) , ChildrenMatrixPer(L) , CountRank(L,l) ,  EQMatrix(LMT) , GFMatrixL(L,t) , GFMatrixPer(L,t) , GuessRec1(L,d) , GuessRec(L) , GFfamliyDet(A,t,m,n) , GFfamilyPer(A,t,m,n) , MatrixGPn(G,n) , MatrixMinor(M) , GFGPn(G,t,m,n) , Ver11N(G,n,v) , VerFn(G,n,v) , VerST(G,v,N) , VerGF(G,N,v,t) , VerGFPreCal(n,v,t) , CtoR(S,t) , SeqFromRec(S,N)`):
end:

ezra:=proc()

if args=NULL then
print(`The main procedures are: DiagMatrix, DiagMatrixL , ExpandMatrixL , ExpandMatrixPer , ChildrenMatrixL , ChildrenMatrixPer , GFMatrixL, GFMatrixPer , GuessRec , GFfamilyDet , GFfamilyPer , MatrixGPn , GFGPn , VerGF`):

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,3,2,6);`):

elif nops([args])=1 and op(1,[args])=DiagMatrixL then
print(`DiagMatrixL(L): input a data structure L=[dim, first_row=[], first_col=[]], where first_row and first_col are two lists of entries in first row and first column from the first one to the last non-zero one. Try: `):
print(`DiagMatrixL([10,[1,2,3],[1,4]]);`):

elif nops([args])=1 and op(1,[args])=ExpandMatrixL then
print(`ExpandMatrixL(L,L1): input a data structure L=[dim, first_row=[], first_col=[]] as the matrix we start and the other data structure L1 as the current submatrix we have, expand L1 along its first row. Output a list of [multiplicity, data structure]. (for simplicity, if the elements in L[2] and L[3] are numbers, they are unique. Otherwise they are symbolic). Try: `):
print(`ExpandMatrixL([20,[1,2,3],[1,4]],[20,[1,2,3],[1,4]]);`):

elif nops([args])=1 and op(1,[args])=ExpandMatrixPer then
print(`ExpandMatrixPer(L,L1): input a data structure L=[dim, first_row=[], first_col=[]] as the matrix we start and the other data structure L1 as the current submatrix we have, expand L1 along its first row for permanent. Output a list of [multiplicity, data structure]. (for simplicity, if the elements in L[2] and L[3] are numbers, they are unique. Otherwise they are symbolic). Try: `):
print(`ExpandMatrixPer([20,[1,2,3],[1,4]],[20,[1,2,3],[1,4]]);`):

elif nops([args])=1 and op(1,[args])=ChildrenMatrixL then
print(`ChildrenMatrixL(L): input a data structure L=[dim, first_row=[], first_col=[]], for large enough dim, output the set of all children despite the dim. Try:  `):
print(`ChildrenMatrixL([20,[1,2,3],[1,4]]);`):

elif nops([args])=1 and op(1,[args])=ChildrenMatrixPer then
print(`ChildrenMatrixPer(L): input a data structure L=[dim, first_row=[], first_col=[]], for large enough dim, output the set of all children despite the dim for the permanent instead of determinant. Try: `):
print(`ChildrenMatrixPer([20,[1,2,3],[1,4]]);`):

elif nops([args])=1 and op(1,[args])=GFMatrixL then
print(`GFMatrixL(L,t): input a data structure L=[dim, first_row=[], first_col=[]], and a symbol t, for large enough dim output the generating function for the determinant. Try: `):
print(`GFMatrixL([20,[1,2,3],[1,4]],t);`):

elif nops([args])=1 and op(1,[args])=GFMatrixPer then
print(`GFMatrixPer(L,t): input a data structure L=[dim, first_row=[], first_col=[]], and a symbol t, for large enough dim output the generating function for the permanent. Try: `):
print(`GFMatrixPer([20,[1,2,3],[1,4]],t);`):

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])=GFfamilyDet then
print(`GFfamilyDet(A,t,m,n): inputs (i) A: a name of a Maple procedure that inputs an integer n and outputs an n by n matrix according to some rule (ii) a variable name t, (iii) two integers m and n which are the lower and upper bounds of the sequence of determinants we consider. output: A (conjectured) rational function in t, let's call it R(t), that is the generating function of the sequence. Try: `):
print(`GFfamilyDet(SampleA, t, 5, 30);`):

elif nops([args])=1 and op(1,[args])=GFfamilyPer then
print(`inputs (i) A: a name of a Maple procedure that inputs an integer n and outputs an n by n matrix according to some rule (ii) a variable name t, (iii) two integers m and n which are the lower and upper bounds of the sequence of determinants we consider.  output: A (conjectured) rational function in t, let's call it R(t), that is the generating function of the permanent sequence. Try: `):
print(`GFfamilyPer(SampleA, t, 5, 30);`):

elif nops([args])=1 and op(1,[args])=MatrixGPn then
print(`MatrixGPn(G, n): input a graph G=[V,E] where V is a set of vertices (starts from 1 to nops(v)) and E is a set of edges, output its Laplacian matrix (a_ij = -1 if {i,j} in E and the diagonal entry is (-1)*the sum of the rest entries on that row). Try: `):
print(`MatrixGPn([{1,2,3,4}, {{1,2},{2,3},{3,4}}], 3);`):

elif nops([args])=1 and op(1,[args])=GFGPn then
print(`GFGPn(G,t,m,n): input a graph G=[V,E], a symbol t and two integers m and n which are the lower and upper bounds of the sequence of determinants we consider, output a (conjectured) rational function in t, let's call it R(t), that is the generating function of the sequence. Try: `):
print(`GFGPn([{1,2,3}, {{1,2},{2,3}}],t,3,20);`):

elif nops([args])=1 and op(1,[args])=VerGF then
print(`VerGF(G,N,v,t): Input a graph G=[V,E], a number N for the length of lists we get for VerST to guess the recurrence, two symbols v and t, output the generating function add(VerFn(G,n,v)*t^n, n=1..infinity), i.e., the coefficient of t^n, as a polynomial of v, is a weighted counting of spanning trees in G*Pn. The weight for each spanning tree is v^(#of vertical edges in this spanning tree). Try:`):
print(`VerGF([{1,2,3,4},{{1,2},{2,3},{3,4}}],20,v,t); `):


fi:
end:


#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:


#DiagMatrixL(L): input a data structure L=[dim, first_row=[], first_col=[]], where first_row and first_col are two lists of entries in first row and first column from the first one to the last non-zero one. 
DiagMatrixL:=proc(L) local n, r1, c1,p,q,S,M,i:
n:=L[1]:
r1:=L[2]:
c1:=L[3]:
p:=nops(r1)-1:
q:=nops(c1)-1:
if r1[1] <> c1[1] then
  return FAIL:
fi:
S:=[0$(n-1-q), seq(c1[q-i+1],i=0..q-1), op(r1), 0$(n-1-p)]:
M:=[0$n]:
for i from 1 to n do
  M[i]:=[op(max(0,n-1-q)+q+2-i..max(0,n-1-q)+q+1+n-i,S)]:
od:
return M:
end:


#ExpandMatrixL(L,L1): input a data structure L=[dim, first_row=[], first_col=[]] as the matrix we start and the other data structure L1 as the current submatrix we have, expand L1 along its first row. Output a list of [multiplicity, data structure]. (for simplicity, if the elements in L[2] and L[3] are numbers, they are unique. Otherwise they are symbolic). 
ExpandMatrixL:=proc(L,L1) local n,R,C,dim,R1,C1,i,r,S,candidate,newrow,newcol,gu,mu,temp,p,q,j,i1:
n:=L[1]:
R:=L[2]:
C:=L[3]:
p:=nops(R)-1:
q:=nops(C)-1:
dim:=L1[1]:
R1:=L1[2]:
C1:=L1[3]:

if R1=[] or C1=[] then
  return {}:
elif R[1]<>C[1] or R1[1]<>C1[1] or dim>n then 
  return FAIL:
else

  S:={}:

  gu:=[0$(n-1-q), seq(C[q-i+1],i=0..q-1), op(R), 0$(n-1-p)]:
  candidate:=[0$nops(R1),R1[-1]]:
  for i from 1 to nops(R1) do
    mu:=R1[i]:
    for j from n-q to nops(gu) do
      if gu[j]=mu then
        candidate[i]:=gu[j-1]:
      fi:
    od:
  od:
      
  for i from n-q to nops(gu) do
    if gu[i] = R1[2] then
      temp:=i:
      break:
    fi:
  od:

  for i from 1 to nops(R1) do
    if i = 1 then
      mu:=[R1[i]*(-1)^(i+1), [dim-1,[op(i+1..nops(candidate), candidate)], [seq(gu[temp-i1],i1=1..temp-n+q)]]]:
      S:=S union {mu}:
    else
      mu:=[R1[i]*(-1)^(i+1), [dim-1, [op(1..i-1, candidate), op(i+1..nops(candidate), candidate)], [op(2..nops(C1), C1)]]]:
      S:=S union {mu}:
    fi:
  od:

  return S:

fi:

end:


#ExpandMatrixPer(L,L1): input a data structure L=[dim, first_row=[], first_col=[]] as the matrix we start and the other data structure L1 as the current submatrix we have, expand L1 along its first row for permanent. Output a list of [multiplicity, data structure]. (for simplicity, if the elements in L[2] and L[3] are numbers, they are unique. Otherwise they are symbolic). 
ExpandMatrixPer:=proc(L,L1) local n,R,C,dim,R1,C1,i,r,S,candidate,newrow,newcol,gu,mu,temp,p,q,j,i1:
n:=L[1]:
R:=L[2]:
C:=L[3]:
p:=nops(R)-1:
q:=nops(C)-1:
dim:=L1[1]:
R1:=L1[2]:
C1:=L1[3]:

if R1=[] or C1=[] then
  return {}:
elif R[1]<>C[1] or R1[1]<>C1[1] or dim>n then 
  return FAIL:
else

  S:={}:

  gu:=[0$(n-1-q), seq(C[q-i1+1],i1=0..q-1), op(R), 0$(n-1-p)]:
  candidate:=[0$nops(R1),R1[-1]]:
  for i from 1 to nops(R1) do
    mu:=R1[i]:
    for j from n-q to nops(gu) do
      if gu[j]=mu then
        candidate[i]:=gu[j-1]:
      fi:
    od:
  od:
      
  for i from n-q to nops(gu) do
    if gu[i] = R1[2] then
      temp:=i:
      break:
    fi:
  od:

  for i from 1 to nops(R1) do
    if i = 1 then
      mu:=[R1[i], [dim-1,[op(i+1..nops(candidate), candidate)], [seq(gu[temp-i1],i1=1..temp-n+q)]]]:
      S:=S union {mu}:
    else
      mu:=[R1[i], [dim-1, [op(1..i-1, candidate), op(i+1..nops(candidate), candidate)], [op(2..nops(C1), C1)]]]:
      S:=S union {mu}:
    fi:
  od:

  return S:

fi:

end:




#ChildrenMatrixL(L): input a data structure L=[dim, first_row=[], first_col=[]], for large enough dim, output the set of all children despite the dim. 
ChildrenMatrixL:=proc(L) local S,t,T,dim,U,u,s:
dim:=L[1]:
S:={[op(2..3,L)]}:
T:={seq([op(2..3,t[2])],t in ExpandMatrixL(L,L))}:
while T minus S <> {} do
  U:=T minus S:
  S:=S union T:
  T:={}:
  for u in U do
    T:=T union {seq([op(2..3,t[2])],t in ExpandMatrixL(L,[dim,op(u)]))}:
  od:
od:
for s in S do
  if s[1]=[] or s[2]=[] then
    S:=S minus {s}:
  fi:
od:
S:
end:


#ChildrenMatrixPer(L): input a data structure L=[dim, first_row=[], first_col=[]], for large enough dim, output the set of all children despite the dim for the permanent instead of determinant. 
ChildrenMatrixPer:=proc(L) local S,t,T,dim,U,u,s:
dim:=L[1]:
S:={[op(2..3,L)]}:
T:={seq([op(2..3,t[2])],t in ExpandMatrixPer(L,L))}:
while T minus S <> {} do
  U:=T minus S:
  S:=S union T:
  T:={}:
  for u in U do
    T:=T union {seq([op(2..3,t[2])],t in ExpandMatrixPer(L,[dim,op(u)]))}:
  od:
od:
for s in S do
  if s[1]=[] or s[2]=[] then
    S:=S minus {s}:
  fi:
od:
S:
end:


#CountRank(L, l): input a list of unique elements L and an element l, output the index of the element in the list.
CountRank:=proc(L,l) local i:
for i from 1 to nops(L) do
  if L[i]=l then
    return i:
  fi:
od:
return FAIL:
end:


#EQMatrixL(L,t): input a data structure L=[dim, first_row=[], first_col=[]], and a symbol t, for large enough dim output the equations for its determinant.
EQMatrixL:=proc(L,t) local S,dim,var,eq,n,A,i,result,gu,mu:
dim:=L[1]:
S:=ChildrenMatrixL(L):
S:=[[op(2..3,L)], op(S minus {[op(2..3,L)]})]:
n:=nops(S):
var:={seq(A[i],i=1..n)}:
eq:={}:
for i from 1 to 1 do
  result:=ExpandMatrixL(L,[dim,op(S[i])]):
  for gu in result do
  if gu[2][2]=[] or gu[2][3]=[] then
    result:=result minus {gu}:
  fi:
  od:
  eq:=eq union {A[i] - 1 - add(gu[1]*t*A[CountRank(S, [op(2..3, gu[2])])], gu in result)}:
od:
for i from 2 to n do
  result:=ExpandMatrixL(L,[dim,op(S[i])]):
  for gu in result do
  if gu[2][2]=[] or gu[2][3]=[] then
    result:=result minus {gu}:
  fi:
  od:
  eq:=eq union {A[i] - add(gu[1]*t*A[CountRank(S, [op(2..3, gu[2])])], gu in result)}:
od:
eq:
end:


#GFMatrixL(L,t): input a data structure L=[dim, first_row=[], first_col=[]], and a symbol t, for large enough dim output the generating function for the determinant.
GFMatrixL:=proc(L,t) local S,dim,var,eq,n,A,i,result,gu,mu:
dim:=L[1]:
S:=ChildrenMatrixL(L):
S:=[[op(2..3,L)], op(S minus {[op(2..3,L)]})]:
n:=nops(S):
var:={seq(A[i],i=1..n)}:
eq:={}:
for i from 1 to 1 do
  result:=ExpandMatrixL(L,[dim,op(S[i])]):
  for gu in result do
  if gu[2][2]=[] or gu[2][3]=[] then
    result:=result minus {gu}:
  fi:
  od:
  eq:=eq union {A[i] - 1 - add(gu[1]*t*A[CountRank(S, [op(2..3, gu[2])])], gu in result)}:
od:
for i from 2 to n do
  result:=ExpandMatrixL(L,[dim,op(S[i])]):
  for gu in result do
  if gu[2][2]=[] or gu[2][3]=[] then
    result:=result minus {gu}:
  fi:
  od:
  eq:=eq union {A[i] - add(gu[1]*t*A[CountRank(S, [op(2..3, gu[2])])], gu in result)}:
od:
gu:=solve(eq, var)[1]:
subs(gu, A[1]):

end:


#GFMatrixPer(L,t): input a data structure L=[dim, first_row=[], first_col=[]], and a symbol t, for large enough dim output the generating function for the permanent.
GFMatrixPer:=proc(L,t) local S,dim,var,eq,n,A,i,result,gu,mu:
dim:=L[1]:
S:=ChildrenMatrixPer(L):
S:=[[op(2..3,L)], op(S minus {[op(2..3,L)]})]:
n:=nops(S):
var:={seq(A[i],i=1..n)}:
eq:={}:
for i from 1 to 1 do
  result:=ExpandMatrixPer(L,[dim,op(S[i])]):
  for gu in result do
  if gu[2][2]=[] or gu[2][3]=[] then
    result:=result minus {gu}:
  fi:
  od:
  eq:=eq union {A[i] - 1 - add(gu[1]*t*A[CountRank(S, [op(2..3, gu[2])])], gu in result)}:
od:
for i from 2 to n do
  result:=ExpandMatrixPer(L,[dim,op(S[i])]):
  for gu in result do
  if gu[2][2]=[] or gu[2][3]=[] then
    result:=result minus {gu}:
  fi:
  od:
  eq:=eq union {A[i] - add(gu[1]*t*A[CountRank(S, [op(2..3, gu[2])])], gu in result)}:
od:
gu:=solve(eq, var)[1]:
subs(gu, A[1]):

end:



#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:


#GFfamilyDet(A,t,m,n): inputs (i) A: a name of a Maple procedure that inputs an integer n and outputs an n by n matrix according to some rule (ii) a variable name t, (iii) two integers m and n which are the lower and upper bounds of the sequence of determinants we consider. output: A (conjectured) rational function in t, let's call it R(t), that is the generating function of the sequence.
GFfamilyDet:=proc(A,t,m,n) local i,rec,GF,B,gu,Denom,L,Numer:
L:=[seq(det(A(i)),i=1..n)]:
rec:=GuessRec([op(m..n,L)])[2]:
gu:=solve(B-1-add(t^i*rec[i]*B,i=1..nops(rec)), {B}):
Denom:=denom(subs(gu,B)):
Numer:=Denom*(1+add(L[i]*t^i, i=1..n)):
Numer:=add(coeff(Numer,t,i)*t^i, i=0..degree(Denom,t)):
Numer/Denom:
end: 


#GFfamilyPer(A,t,m,n): inputs (i) A: a name of a Maple procedure that inputs an integer n and outputs an n by n matrix according to some rule (ii) a variable name t, (iii) two integers m and n which are the lower and upper bounds of the sequence of determinants we consider.  output: A (conjectured) rational function in t, let's call it R(t), that is the generating function of the sequence.
GFfamilyPer:=proc(A,t,m,n) local i,rec,GF,B,gu,Denom,L,Numer:
L:=[seq(permanent(A(i)),i=1..n)]:
rec:=GuessRec([op(m..n,L)])[2]:
gu:=solve(B-1-add(t^i*rec[i]*B,i=1..nops(rec)), {B}):
Denom:=denom(subs(gu,B)):
Numer:=Denom*(1+add(L[i]*t^i, i=1..n)):
Numer:=add(coeff(Numer,t,i)*t^i, i=0..degree(Denom,t)):
Numer/Denom:
end: 



#SampleA(n): sample input for GFfamilyDet and GFfamilyPer.
SampleA:=proc(n) local L,M:
L:=[n, [1,2], [1,3]]:
M:=DiagMatrixL(L):
end:


#SampleB(n): sample input for GFfamilyDet and GFfamilyPer.
SampleB:=proc(n) local L,M:
L:=[n, [2,3], [2,4,5]]:
M:=DiagMatrixL(L):
end:


#SampleC(n): sample input for GFfamilyDet and GFfamilyPer.
SampleC:=proc(n) local L,M:
L:=[n, [1,2,3,4,5], [1,6,7]]:
M:=DiagMatrixL(L):
end:


#MatrixGPn(G, n): input a graph G=[V,E] where V is a set of vertices (starts from 1 to nops(v)) and E is a set of edges, output its Laplacian matrix (a_ij = -1 if {i,j} in E and the diagonal entry is (-1)*the sum of the rest entries on that row). 
MatrixGPn:=proc(G,n) local nv,dim,i,j,M,edge,v1,v2:
if n<1 then
  return FAIL:
fi:
nv:=nops(G[1]):
dim:=nv*n:
M:=[[0$dim]$dim]:
for edge in G[2] do
  v1:=edge[1]:
  v2:=edge[2]:
  for i from 0 to n-1 do
    M[i*nv+v1][i*nv+v2]:=-1:
    M[i*nv+v2][i*nv+v1]:=-1:
  od:
od:

for i from 1 to n-1 do
  for j from 1 to nv do
    v1:=j+i*nv:
    v2:=j+(i-1)*nv:
    M[v1][v2]:=-1:
    M[v2][v1]:=-1:
  od:
od:

for i from 1 to dim do
  M[i][i]:=(-1)*(add(M[i][j], j=1..i-1) + add(M[i][j], j=i+1..dim)):
od:

M:
end:


#MatrixMinor(M): input a square matrix, output the minor of the matrix by deleting the first row and the first coloum.
MatrixMinor:=proc(M) local n,i:
n:=nops(M):
[seq([op(2..n, M[i])], i=2..n)]:
end:


#GFGPn(G,t,m,n): input a graph G=[V,E], a symbol t and two integers m and n which are the lower and upper bounds of the sequence of determinants we consider, output a (conjectured) rational function in t, let's call it R(t), that is the generating function of the sequence.
GFGPn:=proc(G,t,m,n) local i,L,rec,gu,B,Denom,Numer:
L:=[seq(det(MatrixMinor(MatrixGPn(G,i))), i=1..n)]:
rec:=GuessRec([op(m..n,L)])[2]:
gu:=solve(B-1-add(t^i*rec[i]*B,i=1..nops(rec)), {B}):
Denom:=denom(subs(gu,B)):
Numer:=Denom*(1+add(L[i]*t^i, i=1..n)):
Numer:=add(coeff(Numer,t,i)*t^i, i=0..degree(Denom,t)):
(Numer-Denom)/Denom:
end: 


#Ver11N(G,n,v): Given a connected graph G=[V,E] and a positive integer n , finds the matrix for the graph G*Pn whose determinant would give the Laplacian matrix, especially for vertical edges let the entry be -v. Try:
#Ver11N([{1,2,3,4},{{1,2},{2,3},{3,4}}],3,v);
Ver11N:=proc(G,n,v) local nv,dim,i,j,M,edge,v1,v2:
if n<1 then:
  return FAIL:
fi:
nv:=nops(G[1]):
dim:=nv*n:
M:=[[0$dim]$dim]:
for edge in G[2] do
  v1:=edge[1]:
  v2:=edge[2]:
  for i from 0 to n-1 do
    M[i*nv+v1][i*nv+v2]:=-v:
    M[i*nv+v2][i*nv+v1]:=-v:
  od:
od:

for i from 1 to n-1 do
  for j from 1 to nv do
    v1:=j+i*nv:
    v2:=j+(i-1)*nv:
    M[v1][v2]:=-1:
    M[v2][v1]:=-1:
  od:
od:

for i from 1 to dim do
  M[i][i]:=(-1)*(add(M[i][j], j=1..i-1) + add(M[i][j], j=i+1..dim)):
od:

M:
end:


#VerFn(G,n,v): Given a connected graph G=[V,E] and a positive integer n , finds the determinant of the minor of the matrix for the graph G*Pn whose determinant would give the Laplacian matrix, especially for vertical edges let the entry be -v. Try:
#VerFn([{1,2,3,4},{{1,2},{2,3},{3,4}}],3,v);
VerFn:=proc(G,n,v)
det(MatrixMinor(Ver11N(G,n,v))):
end:


#VerST(G,v,N): inputs a graph  G=[V,E], a symbol v and a positive integer N, outputs
#the first N terms of the sequence enumerating the spanning trees of Gx{1,...,N1} (closed version) for N1 from 1 to N as a weight counting for the number of vertical edges in the spanning trees.
#Try:
#VerST([{1,2,3,4},{{1,2},{2,3},{3,4}}],v,10);
VerST:=proc(G,v,N) local n:
[seq(VerFn(G,n,v), n=1..N)]:
end:


#VerGF(G,N,v,t): Input a graph G=[V,E], a number N for the length of lists we get for VerST to guess the recurrence, two symbols v and t, output the generating function add(VerFn(G,n,v)*t^n, n=1..infinity), i.e., the coefficient of t^n, as a polynomial of v, is a weighted counting of spanning trees in G*Pn. The weight for each spanning tree is v^(#of vertical edges in this spanning tree). Try:
#VerGF([{1,2,3,4},{{1,2},{2,3},{3,4}}],20,v,t);
VerGF:=proc(G,N,v,t) local L:
L:=VerST(G,v,N):
CtoR(GuessRec(L), t)*t:
end:


#VerGFPreCal(n,v,t): input a positive integer 2<=n<=4, output a pre-calculated generating function such that the coefficients of t^m (as a polynomial of v) is weighted counting with weight v^(#of vertical edges) among the set of spanning trees of the n by m grid. 
VerGFPreCal:=proc(n,v,t) 
if n=2 then
  return v*t/(1-(2*v+2)*t+t^2):
elif n=3 then
  return (-t^2*v^2+v^2)*t/(1-(3*v^2+8*v+4)*t-(-10*v^2-16*v-6)*t^2-(3*v^2+8*v+4)*t^3+t^4):
elif n=4 then
  return t*(v^3+(-16*v^5-24*v^4-9*v^3)*t^2+(8*v^6+40*v^5+48*v^4+16*v^3)*t^3+(-16*v^5-24*v^4-9*v^3)*t^4+v^3*t^6)/(1-(4*v^3+20*v^2+24*v+8)*t-(-52*v^4-192*v^3-256*v^2-144*v-28)*t^2-(64*v^5+416*v^4+892*v^3+844*v^2+360*v+56)*t^3-(-16*v^6-160*v^5-744*v^4-1408*v^3-1216*v^2-480*v-70)*t^4-(64*v^5+416*v^4+892*v^3+844*v^2+360*v+56)*t^5-(-52*v^4-192*v^3-256*v^2-144*v-28)*t^6-(4*v^3+20*v^2+24*v+8)*t^7+t^8):
fi:
end:



  

############ From CFinite.txt ############

#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:


  
############ End of procedures from Finite.txt ############