######################################################################
##KamaShidukhim: Save this file as  KamaShidukhim                    #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read KamaEtzim<Enter>                              #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: April 5, 2011
 
print(`Created:  April 5, 2011`):
print(` This is KamaShidukhim `):
print(`It accompanyies the article `):
print(`Automatic Generation of `):
print(`Generating Functions for Enumerating the Number of Perfect Matchings`):
print(`for Grid Graphs (and much more general creatures) of fixed `):
print(`(but arbitrary!) width `):
print(`by Shalosh B. Ekhad and Doron Zeilberger`):
print(`The article is exclusively published in the Personal Journal of`):
print(`Shalosh B. Ekhad and Doron Zeilberger `):
print(` http://www.math.rutgers.edu/~zeilberg/pj.html `):
print(`and arxiv.org `):
print(`the direcet url being:`):
print(`http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/shidukhim.html`):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package `):
 print(` is  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/tokhniot/KamaShidukhim  .`):

print(`-----------------------------------------------`):
 print(`For a list of the Checking procedures type ezraC();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`-----------------------------------------------`):

print(`-----------------------------------------------`):
 print(`For a list of the Monomer-Dimer procedures type ezraMD();,`):
 print(` for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`-----------------------------------------------`):

 print(`For a list of the supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`-----------------------------------------------`):

 print(`For a list of the Story procedures type ezraS();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`-----------------------------------------------`):




 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`-----------------------------------------------`):

with(combinat):
with(linalg):


ezraMD:=proc()

if args=NULL then
 print(` The Monomer-Dimer procedures are: `):
 print(` GFtMD, MD, GFcylMD, GFGmd, GFkingMD, GFrectMD, `):

else
ezra(args):
fi:

end:

ezraS:=proc()

if args=NULL then
 print(` The story procedures are: `):
 print(`SipurC, SipurCmd, SipurCG, SipurCGmd, SipurK, SipurKmd,`):
  print(`SipurR,  SipurRmd  `):

else
ezra(args):
fi:

end:



ezraC:=proc()

if args=NULL then
 print(` The checking procedures are: `):
 print(` CheckGFt, CheckGFtMD, Cyl, Rect `):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: `):
 print(` CBPgraph, CC, CylGC, Children, CylGraph, `):
 print(` GetEqns, Mn, PathGraph,`):
 print(` PM, PPM, PPMmd, RectGC, SeqRectSlow, SeqSlow, TM , ULgraphs `):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:  GFcyl,`):
print(` GFking, GFG, GFrect `):
 print(` GFt ,  Seq, SeqCyl, SeqRect `):
 print(`  `):



elif nops([args])=1 and op(1,[args])=CBPgraph then
print(`CBPgraph(m): The path graph of width m but with all`):
print(`possible edges between consecutive levels`):
print(`For example, try:`):
print(`CBPgraph(3);`):


elif nops([args])=1 and op(1,[args])=CC then
print(`CC(i,n,G): The connected component of vertex i in the graph G`):

elif nops([args])=1 and op(1,[args])=CheckGFt then
print(`CheckGFt(G,C,m,K): checks, by comparing to the completely`):
print(`brute-force enumeration, the first K terms of GFt(G,C,m,z)`):
print(`For example, type:`):
print(`CheckGFt({{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}},4,4);`):

elif nops([args])=1 and op(1,[args])=CheckGFtMD then
print(`CheckGFtMD(G,C,m,K): checks, by comparing to the completely`):
print(`brute-force enumeration, the first K terms of GFtMD(G,C,m,z)`):
print(`For example, type:`):
print(`CheckGFtMD({{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}},4,4);`):

elif nops([args])=1 and op(1,[args])=Children then
print(`Children(m,G,C,V1): Given a pos. integer m, a graph G`):
print(`on vertices {1,..m}, (given as a set of edges), and`):
print(`a bipartite graph from {1, ...,m} to {m+1, ..., 2m}, and`):
print(`a subset V of {1, ,m} finds all`):
print(`subsets of {1,...,m}`):
print(`at the next level followed. For example, try:`):
print(`Children(4,{{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}},{1,2,3,4});`):

elif nops([args])=1 and op(1,[args])=Cyl then
print(`Cyl(m,n): The m by n cylinder-graph`):
print(`For example, try:`):
print(`Cyl(3,4);`):

elif nops([args])=1 and op(1,[args])=CylGC then
print(`RectGC(m): The m by n cylinder motive`):
print(`For example, try: CylGC(5);`):


elif nops([args])=1 and op(1,[args])=GetEqns then
print(`GetEqns(G,C,m,z,f): The system of equations for GFt(G,C,m,z)`):
print(`For example, try:`):
print(`GetEqns({{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}},4,z,f);`):

elif nops([args])=1 and op(1,[args])=GFcyl then
print(`GFcyl(m,z): The generating function, in z, whose`):
print(`coeff of z^n in its Maclaurin expansion would give`):
print(`you the number of perfect matchings in the m by n cylinder graph.`):
print(`For example, try:`):
print(`GFcyl(3,z);`):


elif nops([args])=1 and op(1,[args])=GFcylMD then
print(`GFcylMD(m,z): The generating function, in z, whose`):
print(`coeff of z^n in its Maclaurin expansion would give`):
print(`you the number of matchings (not necessarily perfect)`):
print(` in the m by n cylinder graph.`):
print(`For example, try:`):
print(`GFcylMD(3,z);`):

elif nops([args])=1 and op(1,[args])=GFG then
print(`GFG(G,m,z): Inputs a graph G on m vertices `):
print(`(given as a set of edges)`):
print(` a variable z and a pos. integer K`):
print(`(large enough for the guessing), and outputs`):
print(`the generating function (a rational function in  z)`):
print(`whose coefficient of z^n in its Maclaurin expansion gives`):
print(`you the number of perfect matchings of GxP_n. `):
print(`For example, try:`):
print(` GFG(choose({1,2,3,4},2) minus {{1,3}},4,z); `):

elif nops([args])=1 and op(1,[args])=GFGmd then
print(`GFGmd(G,m,z): Inputs a graph G on m vertices `):
print(`(given as a set of edges)`):
print(` a variable z and a pos. integer K`):
print(`(large enough for the guessing), and outputs`):
print(`the generating function (a rational function in  z)`):
print(`whose coefficient of z^n in its Maclaurin expansion gives`):
print(`you the total number of matchings (not necessarily perfect)`):
print(`of GxP_n. `):
print(`For example, try:`):
print(` GFGmd(choose({1,2,3,4},2) minus {{1,3}},4,z); `):

elif nops([args])=1 and op(1,[args])=GFking then
print(`GFking(m,z): The generating function, in z, whose`):
print(`coeff of z^n in its Maclaurin expansion would give`):
print(`you the number of perfect matchings in the m by n Chess King graph.`):
print(`For example, try:`):
print(`GFking(3,z);`):

elif nops([args])=1 and op(1,[args])=GFkingMD then
print(`GFkingMD(m,z): The generating function, in z, whose`):
print(`coeff of z^n in its Maclaurin expansion would give`):
print(`you the number of matchings (not necessarily perfect)`):
print(`in the m by n Chess King graph.`):
print(`For example, try:`):
print(`GFkingMD(3,z);`):

elif nops([args])=1 and op(1,[args])=GFrect then
print(`GFrect(m,z): The generating function, in z, whose`):
print(`coeff of z^n in its Maclaurin expansion would give`):
print(`you the number of perfect matchings (alias domino tilings,`):
print(` also called dimer tilings ) `):
print(` in the m by n grid-graph.`):
print(`For example, try:`):
print(`GFrect(3,z);`):


elif nops([args])=1 and op(1,[args])=GFrectMD then
print(`GFrectMD(m,z): The generating function, in z, whose`):
print(`coeff of z^n in its Maclaurin expansion would give`):
print(`you the number of (not necessarily perfect)`):
print(` matchings (alias monomer-domino tilings,`):
print(` in the m by n grid-graph.`):
print(`For example, try:`):
print(`GFrectMD(3,z);`):

elif nops([args])=1 and op(1,[args])=GFt then
print(`GFt(G,C,m,z): The generating function (a rational function in  z)`):
print(`whose coefficient of z^k in its Maclaurin expansion gives`):
print(`you the number of perfect matchings of Mn(G,C,m,k).`):
print(`For example, try:`):
print(`GFt({{1,2},{2,3}},{{1,4},{2,5},{3,6}},3,z);`):

elif nops([args])=1 and op(1,[args])=GFtMD then
print(`GFtMD(G,C,m,z): The generating function (a rational function in  z)`):
print(`whose coefficient of z^k in its Maclaurin expansion gives`):
print(`you the number of monomer-dimer tilings of Mn(G,C,m,k).`):
print(`For example, try:`):
print(`GFtMD({{1,2},{2,3}},{{1,4},{2,5},{3,6}},3,z);`):


elif nops([args])=1 and op(1,[args])=MD then
print(`MD(G): The set of monomer-dimer tiltings of the graph G=[V,E]`):
print(`For example, try MD([{1,2},{{1,2}}]);`):

elif nops([args])=1 and op(1,[args])=Mn then
print(`Mn(G,C,m,n): Given a graph G (a set of edges) on {1, ..., m}`):
print(`and a biparatite graph C from {1, ..., m} to {m+1, ..., 2m}`):
print(`constructs the so-called mixed product.M_n(G,C)`):
print(`on {1, ..., nm}. For example, try:`):
print(`Mn({{1,2}},{{1,3},{2,4}},2,3);`):

elif nops([args])=1 and op(1,[args])=PathGraph then
print(`PathGraph(m): The path-graph of width m`):
print(`For example, try:`):
print(`PathGraph(3);`):

elif nops([args])=1 and op(1,[args])=PM then
print(`PM(G): The set of perfect matchings of the graph G=[V,E]`):
print(`For example, try PM([{1,2},{{1,2}}]);`):

elif nops([args])=1 and op(1,[args])=PPM then
print(`PPM(G,V1): The set of partial perfect matchings of the graph G=[V,E]`):
print(`where V1 is a subset of V, and where all the members of V1`):
print(`are matched. For example, try`):
print(`For example, try PPM([{1,2},{{1,2}}],{1,2});`):

elif nops([args])=1 and op(1,[args])=PPMmd then
print(`PPMmd(G,V1): `):
print(`The set of partial perfect matchings of the graph G=[V,E]`):
print(`where V1 is a subset of V, and where some of the members of V1`):
print(`are matched. For example, try`):
print(`For example, try PPMmd([{1,2},{{1,2}}],{1,2});`):

elif nops([args])=1 and op(1,[args])=Rect then
print(`Rect(m,n): The m by n grid-graph. For example, try:`):
print(`Rect(4,4);`):

elif nops([args])=1 and op(1,[args])=RectGC then
print(`RectGC(m): The m by n rectangle motive`):
print(`For example, try: RectGC(5);`):


elif nops([args])=1 and op(1,[args])=SeqCyl then
print(`SeqCyl(m,N): The first N terms of the`):
print(`enumerating sequence for the number of perfect matchings`):
print(`of the m by n grid-graph. For example, try:`):
print(`SeqCyl(4,30);`):

elif nops([args])=1 and op(1,[args])=SeqRect then
print(`SeqRect(m,N): The first N terms of the`):
print(`enumerating sequence for the number of perfect matchings`):
print(`of the m by n grid-graph. For example, try:`):
print(`SeqRect(4,30);`):

elif nops([args])=1 and op(1,[args])=SeqRectSlow then
print(`SeqRectSlow(m,N): `):
print(`Like SeqRect(m,N) but using the generating function`):
print(`For example, try:`):
print(`SeqRectSlow(4,30);`):

elif nops([args])=1 and op(1,[args])=Seq then
print(`Seq(G,C,m,N): The first N terms of the`):
print(`enumerating sequence for the number of perfect matchings`):
print(`of Mn(G,C,n,m))`):
print(`Seq({{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}},4,30);`):

elif nops([args])=1 and op(1,[args])=SeqSlow then
print(`SeqSlow(G,C,m,N): The first N terms of the`):
print(`enumerating sequence for the number of perfect matchings`):
print(`of Mn(G,C,n,m))`):
print(`SeqSlow({{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}},4,30);`):

elif nops([args])=1 and op(1,[args])=SipurC then
print(`SipurC(M1,M2,z,N): the story of all generating functions for`):
print(`GFcyl(m,z) for m from M1 to M2. displaying the first N terms`):
print(`For example, try:`):
print(`SipurC(2,4,z,30);`):



elif nops([args])=1 and op(1,[args])=SipurCmd then
print(`SipurCmd(M1,M2,z,N): the story of all generating functions for`):
print(`GFcylMD(m,z) for m from M1 to M2. displaying the first N terms`):
print(`For example, try:`):
print(`SipurCmd(2,4,z,30);`):


elif nops([args])=1 and op(1,[args])=SipurCG then
print(`SipurCG(M,z,N): the story of all generating functions for`):
print(`GFG(G,m,z) for all connected graphs with m vertices up to M`):
print(`listing  N terms in the sequencse`):
print(`For example, try:`):
print(`SipurCG(4,z,30);`) :

elif nops([args])=1 and op(1,[args])=SipurCGmd then
print(`SipurCGmd(M,z,N): the story of all generating functions for`):
print(`GFGmd(G,m,z) for all connected graphs with m vertices up to M`):
print(`listing  N terms in the sequencse`):
print(`For example, try:`):
print(`SipurCGmd(4,z,30);`) :

elif nops([args])=1 and op(1,[args])=SipurK then
print(`SipurK(M1,M2,z,N): the story of all generating functions for`):
print(`GFking(m,z) for m from M1 to M2. displaying the first N terms`):
print(`For example, try:`):
print(`SipurK(1,3,z,30);`):

elif nops([args])=1 and op(1,[args])=SipurKmd then
print(`SipurKmd(M1,M2,z,N): the story of all generating functions for`):
print(`GFkingMD(m,z) for m from M1 to M2. displaying the first N terms`):
print(`For example, try:`):
print(`SipurKmd(1,3,z,30);`):

elif nops([args])=1 and op(1,[args])=SipurR then
print(`SipurR(M1,M2,z,N): the story of all generating functions for`):
print(`GFrect(m,z) for m from M1 to M2. displaying the first N terms`):
print(`For example, try:`):
print(`SipurR(2,5,z,30);`):


elif nops([args])=1 and op(1,[args])=SipurRmd then
print(`SipurRmd(M1,M2,z,N): the story of all generating functions for`):
print(`GFrectMD(m,z) for m from M1 to M2. displaying the first N terms`):
print(`For example, try:`):
print(`SipurRmd(2,5,z,30);`):

elif nops([args])=1 and op(1,[args])=TM then
print(`TM(m,G,C): inputs a pos. integer m, a graph G on {1, ..., m}`):
print(`(given as a set of edges) and a bipartite graph from`):
print(`{1, ...,m} to {m+1, ..., 2m}, C outputs the`):
print(`transition matrix, followed by the list of states.`);
print(`For example, try:`):
print(`TM(4,{{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}});`):

elif nops([args])=1 and op(1,[args])=ULgraphs then
print(`ULgraphs(n): all the unlabelled graphs on n vertices`):
print(`For example, try:`):
print(`ULgraphs(3);`):

else
print(`There is no ezra for`,args):
fi:
 
end:





#Mn(G,C,m,n): Given a graph G (a set of edges) on {1, ..., m}
#and a biparatite graph C from {1, ..., m} to {m+1, ..., 2m}
#constructs the so-called mixed product.M_n(G,C)
#on {1, ..., nm}. For example, try:
#Mn({{1,2}},{{1,3},{2,4}},2,3);
Mn:=proc(G,C,m,n) local i,j:
[{seq(i,i=1..m*n)},
{seq(op(subs({seq(i=i+j*m,i=1..m)},G)),j=0..n-1)} union
{seq(op(subs({seq(i=i+j*m,i=1..2*m)},C)),j=0..n-2)}]:
end:


#PathGraph(m): The path-graph of width m
#For example, try:
#PathGraph(3);
PathGraph:=proc(m) local i:
{seq({i,i+1},i=1..m-1)},{seq({i,i+m},i=1..m)}:
end:


#CylGraph(m): The cylinder-graph of width m
#For example, try:
#CylGraph(3);
CylGraph:=proc(m) local i:
{seq({i,i+1},i=1..m-1),{1,m}},{seq({i,i+m},i=1..m)}:
end:


#CBPgraph(m): The path graph of width m but with all
#possible edges between consecutive levels
#For example, try:
#CBPgraph(3);
CBPgraph:=proc(m) local i,j:
{seq({i,i+1},i=1..m-1)},{seq(seq({i,j+m},i=1..m),j=1..m)}:
end:

#GFrect(m,z): The generating function, in z, whose
#coeff of z^n in its Maclaurin expansion would give
#you the number of perfect matchints the m by n grid-graph.
#For example, try:
#GFrect(3,z);
GFrect:=proc(m,z) local i:
if not type(m,integer) or m<=1 then
 RETURN(FAIL):
fi:
GFt({seq({i,i+1},i=1..m-1)},{seq({i,i+m},i=1..m)},m,z);

end:

#GFrectSlow(m,z): The generating function, in z, whose
#coeff of z^n in its Maclaurin expansion would give
#you the number of perfect matchints the m by n grid-graph.
#For example, try:
#GFrectSlow(3,z);
GFrectSlow:=proc(m,z) local i:
if not type(m,integer) or m<=1 then
 RETURN(FAIL):
fi:
GFtSlow({seq({i,i+1},i=1..m-1)},{seq({i,i+m},i=1..m)},m,z);

end:


#GFrectMD(m,z): The generating function, in z, whose
#coeff of z^n in its Maclaurin expansion would give
#you the number of matchints the m by n grid-graph.
#For example, try:
#GFrectMD(3,z);
GFrectMD:=proc(m,z) local i:
if not type(m,integer) or m<=1 then
 RETURN(FAIL):
fi:
GFtMD({seq({i,i+1},i=1..m-1)},{seq({i,i+m},i=1..m)},m,z);

end:




#GFcyl(m,z): The generating function, in z, whose
#coeff of z^n in its Maclaurin expansion would give
#you the number of perfect matchings in the m by n cylinder graph.
#For example, try:
#GFcyl(3,z);
GFcyl:=proc(m,z) local i:
if not type(m,integer) or m<=1 then
 RETURN(FAIL):
fi:

GFt({seq({i,i+1},i=1..m-1),{1,m}},{seq({i,i+m},i=1..m)},m,z);
end:

#GFcylMD(m,z): The generating function, in z, whose
#coeff of z^n in its Maclaurin expansion would give
#you the number of matchings in the m by n cylinder graph.
#For example, try:
#GFcylMD(3,z);
GFcylMD:=proc(m,z) local i:
if not type(m,integer) or m<=1 then
 RETURN(FAIL):
fi:

GFtMD({seq({i,i+1},i=1..m-1),{1,m}},{seq({i,i+m},i=1..m)},m,z);
end:



#GFking(m,z): The generating function, in z, whose
#coeff of z^n in its Maclaurin expansion would give
#you the number of perfect matchings in the m by n king graph.
#For example, try:
#GFking(3,z);
GFking:=proc(m,z) local i:
if not type(m,integer) or m<=1 then
 RETURN(FAIL):
fi:

GFt({seq({i,i+1},i=1..m-1)},
{seq({i,i+m},i=1..m), seq({i,i+m+1},i=1..m-1),seq({i,i+m-1},i=2..m)},m,z);
end:

#GFkingMD(m,z): The generating function, in z, whose
#coeff of z^n in its Maclaurin expansion would give
#you the number of  matchings in the m by n king graph.
#For example, try:
#GFkingMD(3,z);
GFkingMD:=proc(m,z) local i:
if not type(m,integer) or m<=1 then
 RETURN(FAIL):
fi:

GFtMD({seq({i,i+1},i=1..m-1)},
{seq({i,i+m},i=1..m), seq({i,i+m+1},i=1..m-1),seq({i,i+m-1},i=2..m)},m,z);
end:

#SipurR(M1,M2,z,N): the story of all generating functions for
#GFrect(m,z) for m from M1 to M2. displaying the first N terms
#For example, try:
#SipurR(2,3,z,30);
SipurR:=proc(M1,M2,z,N) local A,m,gu,n,f,t0,i:
t0:=time():
print(`The Generating Functions for the Number of Perfect Matchings in Grid Graphs`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`This is the story of the generating functions, in `, z):
print(`of the sequences enumerating the number of perfect matchings`):
print(`(alias domino tilings) `):
print(`in grid-graphs P_n x P_m, for general n, and m from`, M1, ` to `,M2):

for m from M1 to M2 do
print(`-------------------------------------------------`):
 gu:=GFrect(m,z):
if m mod 2=0 then

print(`Theorem: Let`, A(m,n) , `be the number of perfect matchings`):
print(`in the `, m, `by n  grid graph, and let`):

print(f(z)=Sum(A(m,n)*z^n,n=0..infinity)):

print(`Then : `):
print(f(z)=gu):
print(``):
print(`And in Maple-input format, it is:`):
lprint(gu):

gu:=taylor(gu,z=0,N+1):
print(`The first `, N, `terms are:`):
print([seq(coeff(gu,z,i),i=1..N)]):

else


print(`Theorem: Let`, B(m,n) , `be the number of perfect matchings`):
print(`in the `, m, `by 2n  grid graph, and let`):

print(f(z)=Sum(B(m,n)*z^n,n=0..infinity)):
gu:=subs(z=sqrt(z),gu):
print(`Then : `):
print(f(z)=gu):
print(``):
print(`And in Maple-input format, it is:`):
lprint(gu):

gu:=taylor(gu,z=0,N+1):
print(`The first `, N, `terms are:`):
print([seq(coeff(gu,z,i),i=1..N)]):

fi:

od:

print(`The whole thing took`, time()-t0, `seconds.`):
end:



#SipurC(M1,M2,z,N): the story of all generating functions for
#GFcyl(m,z) for m from M1 to M2. displaying the first N terms
#For example, try:
#SipurC(2,4,z,30);
SipurC:=proc(M1,M2,z,N) local B,m,gu,n,f,t0,i:
t0:=time():
print(`The Generating Functions for the Number of Perfect Matchings in Cylinder Graphs`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`This is the story of the generating functions, in `, z):
print(`of the sequences enumerating the number of perfect matchings`):
print(`in cylinder-graphs P_n x C_m, for general n, and m from`, M1,` to `,M2):

for m from M1 to M2 do
 gu:=GFcyl(m,z):
print(`-------------------------------------------------------`):

if m mod 2=0 then

print(`Theorem: Let`, B(m,n) , `be the number of perfect matchings`):
print(`in the `, m, `by n  cylinder graph, and let`):

print(f(z)=Sum(B(m,n)*z^n,n=0..infinity)):

print(`Then : `):
print(f(z)=gu):
print(``):
print(`And in Maple-input format, it is:`):
lprint(gu):

gu:=taylor(gu,z=0,N+1):
print(`The first `, N, `terms are:`):
print([seq(coeff(gu,z,i),i=1..N)]):

else

print(`Theorem: Let`, C(m,n) , `be the number of perfect matchings`):
print(`in the `, m, `by 2n  cylinder graph, and let`):

print(f(z)=Sum(C(m,n)*z^n,n=0..infinity)):
gu:=subs(z=sqrt(z),gu):
print(`Then : `):
print(f(z)=gu):
print(``):
print(`And in Maple-input format, it is:`):
lprint(gu):

gu:=taylor(gu,z=0,N+1):
print(`The first `, N, `terms are:`):
print([seq(coeff(gu,z,i),i=1..N)]):
fi:
od:

print(`The whole thing took`, time()-t0, `seconds.`):
end:



#GFG(G,m,z): Inputs a graph G on m vertices 
#(given as a set of edges)
# a variable z and a pos. integer K
#(large enough for the guessing), and outputs
#the generating function (a rational function in  z)
#whose coefficient of z^n in its Maclaurin expansion gives
#you the number of perfect matchings of GxP_n.
#For example, try:
#GFG({{1,2},{2,3}},3,z);
GFG:=proc(G,m,z) local C,i:
C:={seq({i,i+m},i=1..m)}:
GFt(G,C,m,z):
end:

#GFGmd(G,m,z): Inputs a graph G on m vertices 
#(given as a set of edges)
# a variable z and a pos. integer K
#(large enough for the guessing), and outputs
#the generating function (a rational function in  z)
#whose coefficient of z^n in its Maclaurin expansion gives
#you the number of matchings of GxP_n.
#For example, try:
#GFGmd({{1,2},{2,3}},3,z);
GFGmd:=proc(G,m,z) local C,i:
C:={seq({i,i+m},i=1..m)}:
GFtMD(G,C,m,z):
end:


#PM(G): The set of perfect matchings of the graph G=[V,E]
#For example, try PM([{1,2},{{1,2}}]);
PM:=proc(G) local V,E,gu,lu,e,v,v1,gu1,gu11,e1,E1:
option remember:
V:=G[1]:
E:=G[2]:
if G[1]={} then
  RETURN({{}}):
fi:
if nops(G[1])=1 then
   RETURN({}):
fi:

v:=V[nops(V)]:
lu:={}:
for e in E do
if member(v,e) then
  lu:=lu union {e}:
fi:
od:

gu:={}:

for e in lu do
v1:=(e minus {v})[1]:
E1:=E:
 for e1 in E1 do
   if member(v,e1) or member(v1,e1) then
    E1:=E1 minus {e1}:
   fi:
 od:

gu1:=PM([V minus {v,v1},E1]):

gu:=gu union {seq({op(gu11),e},gu11 in gu1)}:
od:

gu:

 

end:



#Rect(m,n): The m by n grid-graph. For example, try:
#Rect(4,4);
Rect:=proc(m,n) local i,G,C:
G:={seq({i,i+1},i=1..m-1)}:
C:={seq({i,i+m},i=1..m)}:
[{seq(i,i=1..m*n)},Mn(G,C,m,n)]:
end:


#Cyl(m,n): The m by n cylinder-graph. For example, try:
#Cyl(4,4);
Cyl:=proc(m,n) local i,G,C:
G:={seq({i,i+1},i=1..m-1),{m,1}}:
C:={seq({i,i+m},i=1..m)}:
[{seq(i,i=1..m*n)},Mn(G,C,m,n)]:
end:





#PPM(G,V1): The set of partial perfect matchings of the graph G=[V,E]
#where V1 is a subset of V, and where all the members of V1
#are matched. For example, try
#For example, try PPM([{1,2},{{1,2}}],{1,2});
PPM:=proc(G,V1) local V,E,gu,lu,e,v,v1,gu1,gu11,e1,E1:
option remember:
V:=G[1]:
E:=G[2]:
if V1={} then
  RETURN({{}}):
fi:

v:=V1[nops(V1)]:
lu:={}:
for e in E do
if member(v,e) then
  lu:=lu union {e}:
fi:
od:

gu:={}:

for e in lu do
v1:=(e minus {v})[1]:
E1:=E:
 for e1 in E1 do
   if member(v,e1) or member(v1,e1) then
    E1:=E1 minus {e1}:
   fi:
 od:

gu1:=PPM([V minus {v,v1},E1], V1 minus {v,v1}):

gu:=gu union {seq({op(gu11),e},gu11 in gu1)}:
od:

gu:

 

end:


#Children(m,G,C,V1): Given a pos. integer m, a graph G
#on vertices {1,..m}, (given as a set of edges), and
#a bipartite graph from {1, ...,m} to {m+1, ..., 2m}, and
#a subset V of {1, ,m} finds all followers state
#(set of subsets of {m+1, ... 2*m} reduced by m
#For example, try:
#Children(4,{{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}},{1,2,3,4});
Children:=proc(m,G,C,V1) local V1c,i,E1,e1,mu,lu,mu1,V11,mu11,v11,gu:
V1c:={seq(i,i=1..m)} minus V1:
E1:=G union C:
lu:={}:
for e1 in E1 do
 if e1 intersect V1c={} then
  lu:=lu union {e1}:
 fi:
od:

mu:=PPM([V1 union {seq(i,i=m+1..2*m)},lu],V1):
gu:=[]:

for mu1 in mu do
 V11:={seq(op(mu11), mu11 in mu1)} intersect {seq(i,i=m+1..2*m)}:
  V11:= {seq(i,i=m+1..2*m)} minus V11:
  V11:={seq(v11-m,v11 in V11)}:
  gu:=[op(gu),V11]:
od:

gu:
end:


#CheckGFt(G,C,m,K): checks, by comparing to the completely
#brute-force enumeration, the first K terms of GFt(G,C,m,z)
#For example, type:
#CheckGFt({{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}},4,z,4);
CheckGFt:=proc(G,C,m,K) local lu,i,n1,z:
lu:=GFt(G,C,m,z):
lu:=taylor(lu,z=0,K+1):
evalb([seq(nops(PM(Mn(G,C,m,n1))),n1=1..K)]=
[seq(coeff(lu,z,i),i=1..K)]):
end:



#GFt(G,C,m,z): The generating function (a rational function in  z)
#whose coefficient of z^k in its Maclaurin expansion gives
#you the number of perfect mathcings of Mn(G,C,m,k)  
#(see the ezra for Mn for the definition).
#For example, to get the generating function for the number
#of perfect matching (alias domino tilings) for the 4xn grid, type:
#GFt({{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}},4,z);
GFt:=proc(G,C,m,z)   local ToDo,i,V1,AlreadyDone,eq,var,eq1,lu,lu1,f,G1,e1:
ToDo:={{seq(i,i=1..m)}}:
AlreadyDone:={{}}:
var:={f[{}]}:
eq1:=f[{}]-z-z*f[{seq(i,i=1..m)}]:
eq:={eq1}:

while ToDo<>{} do
 V1:=ToDo[1]:
 AlreadyDone:=AlreadyDone union {V1}:
 ToDo:=ToDo minus {V1}:
 lu:=Children(m,G,C,V1):
 ToDo:=ToDo union convert(lu,set) minus AlreadyDone:
 G1:={}:
 for e1 in G do
  if e1 minus V1={} then
   G1:=G1 union {e1}:
  fi:
od:

   
 eq1:=f[V1]-z*nops(PM([V1,G1]))-
 add(z*f[lu1], lu1 in lu):
 eq:=eq union {eq1}:
 var:=var union {f[V1]}:

od:


var:=solve(eq,var):


normal(1+subs(var,f[{seq(i,i=1..m)}])):

end:


#SipurK(M1,M2,z,N): the story of all generating functions for
#GFking(m,z) for m from M1 to M2. displaying the first N terms
#For example, try:
#SipurK(2,4,z,30);
SipurK:=proc(M1,M2,z,N) local m,gu,n,f,t0,i:
t0:=time():
print(`The Generating Functions for the Number of Perfect Mathchings in Chess King Graphs`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`This is the story of the generating functions, in `, z):
print(`of the sequences enumerating the number of perfect matchings`):
print(`in King graphs P_n x P_m, for general n, and m from`, M1 , `to `,M2):

for m from M1 to M2 do
 gu:=GFking(m,z):
print(`-------------------------------------------------------------------`):

if m mod 2 =0 then

print(`Theorem: Let`, K(m,n) , `be the number of perfect mathchings`):
print(`in the `, m, `by n  Chess King graph, and let`):

print(f(z)=Sum(K(m,n)*z^n,n=0..infinity)):

print(`Then : `):
print(f(z)=gu):
print(``):
print(`And in Maple-input format, it is:`):
lprint(gu):

gu:=taylor(gu,z=0,N+1):
print(`The first `, N, `terms are:`):
print([seq(coeff(gu,z,i),i=1..N)]):

else

print(`Theorem: Let`, K(m,n) , `be the number of perfect mathchings`):
print(`in the `, m, `by 2n  Chess King graph, and let`):

print(f(z)=Sum(K(m,n)*z^n,n=0..infinity)):
gu:=subs(z=sqrt(z),gu):
print(`Then : `):
print(f(z)=gu):
print(``):
print(`And in Maple-input format, it is:`):
lprint(gu):

gu:=taylor(gu,z=0,N+1):
print(`The first `, N, `terms are:`):
print([seq(coeff(gu,z,i),i=1..N)]):

fi:

od:

print(`The whole thing took`, time()-t0, `seconds.`):
end:

#MD(G): The set of monomer-dimer tilings of the graph G=[V,E]
#For example, try MD([{1,2},{{1,2}}]);
MD:=proc(G) local V,E,gu,lu,e,v,v1,gu1,gu11,e1,E1:
option remember:
V:=G[1]:
E:=G[2]:
if G[1]={} then
  RETURN({{}}):
fi:
if nops(G[1])=1 then
   RETURN({{{V[1]}}}):
fi:

v:=V[nops(V)]:
lu:={}:
for e in E do
if member(v,e) then
  lu:=lu union {e}:
fi:
od:

gu:=MD([V minus {v},E minus lu]):
gu:={seq(gu1 union {{v}},gu1 in gu)}:


for e in lu do
v1:=(e minus {v})[1]:
E1:=E:
 for e1 in E1 do
   if member(v,e1) or member(v1,e1) then
    E1:=E1 minus {e1}:
   fi:
 od:

gu1:=MD([V minus {v,v1},E1]):

gu:=gu union {seq({op(gu11),e},gu11 in gu1)}:
od:

gu:
end:









#PPMmd(G,V1): The set of partial perfect matchings of the graph G=[V,E]
#where V1 is a subset of V, and where some of the members of V1
#are matched. For example, try
#For example, try PPMmd([{1,2},{{1,2}}],{1,2});
PPMmd:=proc(G,V1) local gu,gu1:
gu:=powerset(V1):
{seq(op(PPM(G,gu1)),gu1 in gu)}:
end:





#ChildrenMD(m,G,C,V1): Given a pos. integer m, a graph G
#on vertices {1,..m}, (given as a set of edges), and
#a bipartite graph from {1, ...,m} to {m+1, ..., 2m}, and
#a subset V of {1, ,m} finds all followers state
#(set of subsets of {m+1, ... 2*m} reduced by m
#For example, try:
#ChildrenMD(4,{{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}},{1,2,3,4});
ChildrenMD:=proc(m,G,C,V1) local V1c,i,E1,e1,mu,lu,mu1,V11,mu11,v11,gu:
V1c:={seq(i,i=1..m)} minus V1:
E1:=G union C:
lu:={}:
for e1 in E1 do
 if e1 intersect V1c={} then
  lu:=lu union {e1}:
 fi:
od:

mu:=PPMmd([V1 union {seq(i,i=m+1..2*m)},lu],V1):
gu:=[]:

for mu1 in mu do
 V11:={seq(op(mu11), mu11 in mu1)} intersect {seq(i,i=m+1..2*m)}:
  V11:= {seq(i,i=m+1..2*m)} minus V11:
  V11:={seq(v11-m,v11 in V11)}:
  gu:=[op(gu),V11]:
od:

gu:
end:

#GFtMD(G,C,m,z): The generating function (a rational function in  z)
#whose coefficient of z^k in its Maclaurin expansion gives
#you the number of monomer-dimer tilints of Mn(G,C,m,k)  
#(see the ezra for Mn for the definition).
#For example, to get that generating function for the number
#for the 4xn grid, type:
#GFtMD({{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}},4,z);
GFtMD:=proc(G,C,m,z)   local ToDo,i,V1,AlreadyDone,eq,var,eq1,lu,lu1,f,G1,e1:
ToDo:={{seq(i,i=1..m)}}:
AlreadyDone:={{}}:
var:={f[{}]}:
eq1:=f[{}]-z-z*f[{seq(i,i=1..m)}]:
eq:={eq1}:

while ToDo<>{} do
 V1:=ToDo[1]:
 AlreadyDone:=AlreadyDone union {V1}:
 ToDo:=ToDo minus {V1}:
 lu:=ChildrenMD(m,G,C,V1):
 ToDo:=ToDo union convert(lu,set) minus AlreadyDone:
 G1:={}:
 for e1 in G do
  if e1 minus V1={} then
   G1:=G1 union {e1}:
  fi:
od:

   
 eq1:=f[V1]-z*nops(MD([V1,G1]))-
 add(z*f[lu1], lu1 in lu):
 eq:=eq union {eq1}:
 var:=var union {f[V1]}:

od:


var:=solve(eq,var):


normal(1+subs(var,f[{seq(i,i=1..m)}])):

end:


#CheckGFtMD(G,C,m,K): checks, by comparing to the completely
#brute-force enumeration, the first K terms of GFtMD(G,C,m,z)
#For example, type:
#CheckGFtMD({{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}},4,z,4);
CheckGFtMD:=proc(G,C,m,K) local lu,i,n1,z:
lu:=GFtMD(G,C,m,z):
lu:=taylor(lu,z=0,K+1):
evalb([seq(nops(MD(Mn(G,C,m,n1))),n1=1..K)]=
[seq(coeff(lu,z,i),i=1..K)]):
end:




#SipurRmd(M1,M2,z,N): the story of all generating functions for
#GFrectMD(m,z) for m from M1 to M2. displaying the first N terms
#For example, try:
#SipurRmd(2,3,z,30);
SipurRmd:=proc(M1,M2,z,N) local A,m,gu,n,f,t0,i:
t0:=time():
print(`The Generating Functions for the Number of Matchings in Grid Graphs`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`This is the story of the generating functions, in `, z):
print(`of the sequences enumerating the number of matchings`):
print(`NOT NECESSARILY PERFECT`):
print(`(alias monomer-dimer tilings) `):
print(`in grid-graphs P_n x P_m, for general n, and m from`, M1, ` to `,M2):

for m from M1 to M2 do
print(`-------------------------------------------------`):
 gu:=GFrectMD(m,z):

print(`Theorem: Let`, A(m,n) , `be the number of matchings`):
print(`(perfect or otherwise) `):
print(`in the `, m, `by n  grid graph, and let`):

print(f(z)=Sum(A(m,n)*z^n,n=0..infinity)):

print(`Then : `):
print(f(z)=gu):
print(``):
print(`And in Maple-input format, it is:`):
lprint(gu):

gu:=taylor(gu,z=0,N+1):
print(`The first `, N, `terms are:`):
print([seq(coeff(gu,z,i),i=1..N)]):

od:

print(`The whole thing took`, time()-t0, `seconds.`):
end:



#SipurCmd(M1,M2,z,N): the story of all generating functions for
#GFcylMD(m,z) for m from M1 to M2. displaying the first N terms
#For example, try:
#SipurCmd(2,3,z,30);
SipurCmd:=proc(M1,M2,z,N) local A,m,gu,n,f,t0,i:
t0:=time():
print(`The Generating Functions for the Number of Matchings in Cylinder Graphs`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`This is the story of the generating functions, in `, z):
print(`of the sequences enumerating the number of matchings`):
print(`NOT NECESSARILY PERFECT`):
print(`(alias monomer-dimer tilings) `):
print(`in the cylinder graph P_n x C_m, `):
print(` for general n, and m from`, M1, ` to `,M2):

for m from M1 to M2 do
print(`-------------------------------------------------`):
 gu:=GFcylMD(m,z):

print(`Theorem: Let`, A(m,n) , `be the number of matchings`):
print(`(perfect or otherwise) `):
print(`in the `, m, `by n  grid graph, and let`):

print(f(z)=Sum(A(m,n)*z^n,n=0..infinity)):

print(`Then : `):
print(f(z)=gu):
print(``):
print(`And in Maple-input format, it is:`):
lprint(gu):

gu:=taylor(gu,z=0,N+1):
print(`The first `, N, `terms are:`):
print([seq(coeff(gu,z,i),i=1..N)]):

od:

print(`The whole thing took`, time()-t0, `seconds.`):
end:


#SipurKmd(M1,M2,z,N): the story of all generating functions for
#GFkingMD(m,z) for m from M1 to M2. displaying the first N terms
#For example, try:
#SipurKmd(2,3,z,30);
SipurKmd:=proc(M1,M2,z,N) local A,m,gu,n,f,t0,i:
t0:=time():
print(`The Generating Functions for the Number of Matchings in `):
print(`the Chess King graph`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`This is the story of the generating functions, in `, z):
print(`of the sequences enumerating the number of matchings`):
print(`NOT NECESSARILY PERFECT`):
print(`(alias monomer-dimer tilings) `):
print(`in the n by m Chess King graph `):
print(` for general n, and m from`, M1, ` to `,M2):

for m from M1 to M2 do
print(`-------------------------------------------------`):
 gu:=GFkingMD(m,z):

print(`Theorem: Let`, A(m,n) , `be the number of matchings`):
print(`(perfect or otherwise) `):
print(`in the `, m, `by n  Chess King graph, and let`):

print(f(z)=Sum(A(m,n)*z^n,n=0..infinity)):

print(`Then : `):
print(f(z)=gu):
print(``):
print(`And in Maple-input format, it is:`):
lprint(gu):

gu:=taylor(gu,z=0,N+1):
print(`The first `, N, `terms are:`):
print([seq(coeff(gu,z,i),i=1..N)]):

od:

print(`The whole thing took`, time()-t0, `seconds.`):
end:



#GetEqns(G,C,m,z,f): The system of equations for GFt(G,C,m,z)
#whose coefficient of z^k in its Maclaurin expansion gives
#you the number of perfect mathcings of Mn(G,C,m,k)  
#(see the ezra for Mn for the definition).
#For example, to get the generating function for the number
#of perfect matching (alias domino tilings) for the 4xn grid, type:
#GetEqns({{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}},4,f);
GetEqns:=proc(G,C,m,z,f)   local ToDo,i,V1,AlreadyDone,eq,var,eq1,lu,lu1,G1,e1:
ToDo:={{seq(i,i=1..m)}}:
AlreadyDone:={{}}:
var:=[f[{}]]:
eq1:=f[{}]-z-z*f[{seq(i,i=1..m)}]:
eq:={eq1}:

while ToDo<>{} do
 V1:=ToDo[1]:
 AlreadyDone:=AlreadyDone union {V1}:
 ToDo:=ToDo minus {V1}:
 lu:=Children(m,G,C,V1):
 ToDo:=ToDo union convert(lu,set) minus AlreadyDone:
 G1:={}:
 for e1 in G do
  if e1 minus V1={} then
   G1:=G1 union {e1}:
  fi:
od:

   
 eq1:=f[V1]-z*nops(PM([V1,G1]))-
 add(z*f[lu1], lu1 in lu):
 eq:=eq union {eq1}:
 var:=[op(var),f[V1]]:

od:


RETURN(eq,var):
end:




#TM(m,G,C): inputs a pos. integer m, a graph G on {1, ..., m}
#(given as a set of edges) and a bipartite graph from
#{1, ...,m} to {m+1, ..., 2m}, C outputs the
#transition matrix, followed by the list of states.
#For example, try:
#TM(4,{{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}});
TM:=proc(m,G,C) local T,ToDo,i,V1,AD,G1,S,e1,j,M,kak,x,t:
ToDo:={{seq(i,i=1..m)}}:

S[{seq(i,i=1..m)}]:=1:

AD:={}:

while ToDo<>{} do
 V1:=ToDo[1]:
 AD:=AD union {V1}:
 ToDo:=ToDo minus {V1}:
 T[V1]:=Children(m,G,C,V1):
 G1:={}:
 for e1 in G do
  if e1 minus V1={} then
   G1:=G1 union {e1}:
  fi:
od:

 S[V1]:=nops(PM([V1,G1])):
 ToDo:=ToDo union convert(T[V1],set) minus AD:
od:

AD:=[{seq(i,i=1..m)},op(convert(AD minus {{seq(i,i=1..m)}},list))]:


for i from 1 to nops(AD) do
 kak:=add(x[t], t in T[AD[i]]):
for j from 1 to nops(AD) do
 M[i,j]:=coeff(kak,x[AD[j]],1):
od:
od:


[seq([seq(M[i,j],j=1..nops(AD))],i=1..nops(AD))],
[seq(S[AD[i]],i=1..nops(AD))]:
end:




#GFtSlow(G,C,m,z): The generating function (a rational function in  z)
#whose coefficient of z^k in its Maclaurin expansion gives
#you the number of perfect mathcings of Mn(G,C,m,k)  
#(see the ezra for Mn for the definition).
#For example, to get the generating function for the number
#of perfect matching (alias domino tilings) for the 4xn grid, type:
#GFtSlow({{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}},4,z);
GFtSlow:=proc(G,C,m,z)   local M,S,gu,M1,i,j:
gu:=TM(m,G,C):
M:=gu[1]:

for i from 1 to nops(M) do
for j from 1 to nops(M[i]) do
  if i<>j then
     M1[i,j]:=-z*M[i][j]:
  else
   M1[i,j]:=1-z*M[i][j]:
  fi:
od:
od:
M1:=[seq([seq(M1[i,j],i=1..nops(M))],j=1..nops(M))]:
M1:=convert(M1,matrix):
M1:=inverse(M1):
normal(evalm(M1&*gu[2])[1]):
end:



#RectGC(m): The m by n rectangle "motive"
RectGC:=proc(m):{seq({i,i+1},i=1..m-1)},{seq({i,i+m},i=1..m)}:end:

#CylGC(m): The m by n cylinder "motive"
CylGC:=proc(m):{seq({i,i+1},i=1..m-1),{m,1}},{seq({i,i+m},i=1..m)}:end:

#SeqSlow(G,C,m,N): The first N terms of the
#enumerating sequence for the number of perfect matchings
#of Mn(G,C,n,m))
#SeqSlow({{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}},4,30);
SeqSlow:=proc(G,C,m,N) local gu,z:
gu:=GFt(G,C,m,z):  
gu:=taylor(gu,z=0,N+1):
[seq(coeff(gu,z,i),i=1..N)]:
end:


#SeqRectSlow(m,N): The first N terms of the
#enumerating sequence for the number of perfect matchings
#of an m by n grid graph. For example, try:
#SeqRectSlow({{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}},4,30);
SeqRectSlow:=proc(m,N) local gu,z,i:
gu:=GFrect(m,z):  
gu:=taylor(gu,z=0,N+1):
[seq(coeff(gu,z,i),i=1..N)]:
end:

#SeqCylSlow(m,N): The first N terms of the
#enumerating sequence for the number of perfect matchings
#of an m by n cylinder graph. For example, try:
#SeqCylSlow({{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}},4,30);
SeqCylSlow:=proc(m,N) local gu,z,i:
gu:=GFcyl(m,z):  
gu:=taylor(gu,z=0,N+1):
[seq(coeff(gu,z,i),i=1..N)]:
end:


#Seq(G,C,m,N): The first N terms of the
#enumerating sequence for the number of perfect matchings
#of Mn(G,C,n,m))
#Seq({{1,2},{2,3},{3,4}},{{1,5},{2,6},{3,7},{4,8}},4,30);
Seq:=proc(G,C,m,N) local mu,gu,M,V,i,j,i1:

mu:=TM(m,G,C):
M:=mu[1]:
V:=mu[2]:
gu:=[V[1]]:

for i1 from 1 to N-1 do
V:=[seq(add(M[i][j]*V[j],j=1..nops(V)),i=1..nops(V))]:
gu:=[op(gu),V[1]]:
od:
gu:
end:



#SeqRect(m,N): The first N terms of the
#enumerating sequence for the number of perfect matchings
#of the m by n  grid-graph . For example try:
#SeqRect(4,30);
SeqRect:=proc(m,N):
Seq(RectGC(m),m,N):
end:


#SeqCyl(m,N): The first N terms of the
#enumerating sequence for the number of perfect matchings
#of the m by n  cylinder graph . For example try:
#SeqCyl(4,30);
SeqCyl:=proc(m,N):
Seq(CylGC(m),m,N):
end:


#CC(i,n,G): The connected component of vertex i in the graph G
CC:=proc(i,n,G) local N,e,NewGuys,i1,vu:
for i1 from 1 to n do
N[i1]:={}:
od:

for e in G do
N[e[1]]:=N[e[1]] union {e[2]}:
N[e[2]]:=N[e[2]] union {e[1]}:
od:

vu:={i}:

NewGuys:=N[i]:


while NewGuys<>{} do
vu:=vu union NewGuys:

NewGuys:={seq(op(N[i1]), i1 in NewGuys)} minus vu:
od:
vu:
end:

#ULgraphs(n): all the unlabelled graphs on n vertices
#For example, try:
#ULgraphs(n)
ULgraphs:=proc(n) local lu,i,ku,kha,mu,pi,gu,mu1:
lu:=choose({seq(i,i=1..n)},2):
lu:=powerset(lu):
mu:={}:

ku:=permute(n):

while lu<>{} do
 kha:=lu[1]:
 mu:= mu union {kha}:
lu:=lu minus {seq( subs({seq(i=pi[i],i=1..n)},kha),pi in ku)}:
od:
mu:

gu:={}:

for mu1 in mu do
if nops(CC(1,n,mu1))=n then
 gu:=gu union {mu1}:
fi:
od:
gu:

end:


#SipurCG(M,z,N): the story of all generating functions for
#GFG(G,m,z) for all connected graphs with m vertices up
#listing  N terms in the sequences
#For example, try:
#SipurCG(4,z,30);
SipurCG:=proc(M,z,N) local A,m,gu,n,f,t0,i,lu,G,sof,vu,Li:
t0:=time():
sof:=[]:
print(`The Generating Functions for the Number of Perfect Matchings`):
print(` in Cartesian products of Connected Graphs with Path Graphs`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`This is the story of the generating functions, in `, z):
print(`of the sequences enumerating the number of PERFECT matchings`):
print(`(alias domino tilings) `):
print(`in Cartesian products P_n x G, for general n, and all connected`):
print(`graphs G with at most`, M, `vertices `):

for m from 2 to M do
print(`-------------------------------------------------------------`):
lu:=ULgraphs(m):

print(`There are`, nops(lu), `different connected graphs`):
print(` with `, m, `vertices (up to isomorphism)`):

 for G in lu do
print(`         --------------------------                       `):
 gu:=GFG(G,m,z):

print(`Theorem: Let`, A(n) , `be the number of perfect matchings`):
print(`in the Cartesian product of the graph`, G):
print(`with the n-path, and let f(z) be its (ordinary) generating  `):
print(`function, i.e. `):

print(f(z)=Sum(A(n)*z^n,n=0..infinity)):

print(`Then : `):
print(f(z)=gu):
print(``):
print(`And in Maple-input format, it is:`):
lprint(gu):

vu:=taylor(gu,z=0,N+1):
print(`The first `, N, `terms are:`):
Li:=[seq(coeff(vu,z,i),i=1..N)]:
print(Li):
sof:=[op(sof),[G,gu,Li]]:
od:
od:

print(`The whole thing took`, time()-t0, `seconds.`):
sof:
end:



#SipurCGmd(M,z,N): the story of all generating functions for
#GFGmd(G,m,z) for all connected graphs with m vertices up
#listing  N terms in the sequences
#For example, try:
#SipurCGmd(4,z,30);
SipurCGmd:=proc(M,z,N) local A,m,gu,n,f,t0,i,lu,G,sof,vu,Li:
t0:=time():
sof:=[]:
print(`The Generating Functions for the Number of Matchings in Cartesian`):
print(`products of Connected Graphs with Path Graphs`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`This is the story of the generating functions, in `, z):
print(`of the sequences enumerating the number of ALL matchings`):
print(`not necessarily perfect, (alias monomer-dimer tilings) `):
print(`in Cartesian products P_n x G, for general n, and all connected`):
print(`graphs G with at most`, M, `vertices `):

for m from 2 to M do
print(`-------------------------------------------------------------`):
lu:=ULgraphs(m):

print(`There are`, nops(lu), `different connected graphs`):
print(` with `, m, `vertices (up to isomorphism)`):

 for G in lu do
print(`         --------------------------                       `):
 gu:=GFGmd(G,m,z):

print(`Theorem: Let`, A(n) , `be the number of all matchings`):
print(`(not necessarily perfect) in the Cartesian product of the graph`, G):
print(`with the n-path, and let f(z) be its (ordinary) generating  `):
print(`function, i.e. `):

print(f(z)=Sum(A(n)*z^n,n=0..infinity)):

print(`Then : `):
print(f(z)=gu):
print(``):
print(`And in Maple-input format, it is:`):
lprint(gu):

vu:=taylor(gu,z=0,N+1):
print(`The first `, N, `terms are:`):
Li:=[seq(coeff(vu,z,i),i=1..N)]:
print(Li):
sof:=[op(sof),[G,gu,Li]]:
od:
od:

print(`The whole thing took`, time()-t0, `seconds.`):
sof:
end: