#ExpSP.txt: A Maple package by Pablo Blanco and Doron Zeilberger

print(`This is ExpSP.txt, a Maple package by Pablo Blanco and Doron Zeilberger`):
print(`Accompanying their article`):
print(`Automatic Generation of Generating Functions for Enumerating Spanning Trees`):
print(``):
print(`For a list of the MAIN functions, type Help(); for help with a specific function type: Help(FunctionName);`):
print(``):
print(`For a list of the HELP functions, type Help1(); for help with a specific function type: Help(FunctionName);`):

print(``):
print(`For a list of the STORY functions, type HelpS(); for help with a specific function type: Help(FunctionName);`):

with(combinat):



HelpS:=proc(): 

if nargs=0 then 

print(`The Story functions are:`):

print(` HnrS, GnrS   `):

else
 Help(args):
fi:
end:

Help1:=proc(): 

if nargs=0 then 

print(`The HELP functions are:`):

print(` CtoR, Cycle, Hnr,HypCub, GenerateFunct, Cyli, Gnr, GridGr, GuessRec,  LaplMat, Path, RGF, Rseq ,Torus, HnR, GnR`):

else
 Help(args):
fi:
end:


Help:=proc(): 

if nargs=0 then 

print(`The MAIN functions are`):

print(`  GenerateGraph,   NumSpanTree, ProdGraph,  SeqHnr, SeqGnr, SpTreeSeq `): 

elif nargs=1 and args[1]=Path then
print(`Path(n): given positive integer n, returns n,E, where E is the edge-set of a path graph with vertex set [n].`):
print(`Try: Path(5).`):

elif nargs=1 and args[1]=Cycle then
print(`Cycle(n): given positive integer n, returns n,E, where E is the edge-set of a cycle graph with vertex set [n].`):
print(`Try: Cycle(4).`):


elif nargs=1 and args[1]=Cyli then
print(`Cyli(n,a): given positive integer n, a. Returns n,E, where E is the edge-set of a cylinder graph with vertex set [n].`):
print(`The given cylinder graph is obtained by the product of a path on n-a vertices and a cycle of length a.`):
print(`Try: Cyli(n,a).`):

elif nargs=1 and args[1]=GridGr then
print(`GridGr(n,a): given positive integer n, a. Returns n,E, where E is the edge-set of a grid graph with vertex set [n].`):
print(`The given grid graph has dimensions n-a times a (dimensions are counted by vertices).`):
print(`Try: GridGr(6,3).`):

elif nargs=1 and args[1]=ProdGraph then
print(`ProdGraph(n,E,m,F): given positive integers n,m, and sets of edges E,F on [n],[m] (respectively), returns the product`):
print(`of the graphs defined by E and F on their respective vertex-sets. `):
print(`Try: ProdGraph(3,{{1,2},{2,3},{3,1}},  2, {{1,2}}):`):
print(`Equivalently, try: ProdGraph(Cycle(3), Path(2)):`):



elif nargs=1 and args[1]=Hnr then
print(`Hnr(n,r): given integers n>0 and r>0, generates the graph on {1,..,n} where the neighborhood of a vertex i is {i+1,..,min(i+r-1,n)} `):
print(`Returns n,E, where E is the edge-set of such a graph on {1,..,n}.`):

elif nargs=1 and args[1]=HnrS then
print(`HnrS(R,x,N): inputs a positive integer R>=2 and a positive integer N outputs an article about the  genearting functions`):
print(`for the number of spanning trees of the graphs Hnr(n,r), for n>=r, whose vertices are {1, ...,n} and vertex i is connected to {i+1,..,min(i+r,n)}. Try: `):
print(`It uses rigorous guessing using N data points. If it does not make it all the way to R, it stops sooner`):
print(`HnrS(3,x,100);`):

elif nargs=1 and args[1]=HypCub then
print(`HypCub(s): given positive integer s, returns 2^s,E, where E is the edge-set of the hypercube graph with vertex set [2^s].`):
print(`Try: HypCub(3).`):

elif nargs=1 and args[1]=Gnr then
print(`Gnr(n,r): given integers n>0 and r>0, generates the graph on {1,..,n} where the neighborhood of a vertex i is {i+j mod n} union {i-j mod n}; except mod is taken in {1,..,n}.`): # specific type of circulant graph
print(`Returns n,E, where E is the edge-set of such a graph on {1,..,n}.`):


elif nargs=1 and args[1]=GnrS then
print(`GnrS(R,x,N): inputs a positive integer R>=2 and a positive integer N outputs an article about the  genearting functions`):
print(`for the number of spanning trees of the graphs Gnr(n,r), for n>=r, whose vertices are {1, ...,n} and vertex i is connected to {i+1,..,i+r}, where it circles around, so it is regular of degree 2*r. Try: `):
print(`It uses rigorous guessing using N data points. If it does not make it all the way to R, it stops sooner`):
print(`GnrS(3,x,100);`):

elif nargs=1 and args[1]=GenerateGraph then
print(`GenerateGraph(N,p): inputs N and a probability p (optional argument with default p=1/2). Returns a random graph on N vertices with edge probability p. Output is a pair N,E, where E is a set of edges on {1,..,N}.`):

elif nargs=1 and args[1]=GenerateFunct then
print(`GenerateFunct(N): given an integer N > 0, returns a list of length N with random entries in {0,..,N-1}.`):


elif nargs=1 and args[1]=LaplMat then
print( ` Lapl(N,G): inputs a pos. integer N and a set of edges in the graph labeled {1,...,N} outputs the Laplacian matrix as a list of lists. For example, try:`):
print(`LaplG(5,{{1,2},{2,3},{3,5}};`);

elif nargs=1 and args[1]=NumSpanTree then
print(`NumSpanTree(M): given the laplacian matrix M of a graph, computes the number of spanning trees of the corresponding graph`):

elif nargs=1 and args[1]=RGF then
print(`RGF(L,x): Guesses a rational generating function in x for the list of numbers L starting with x^0*L[1]+..., Try:`):
print(`RGF([seq(fibonacci(i),i=0..20)],x);`):

elif nargs=1 and args[1]=RSeq then
print(`RSeq(m): outputs the terms of the sequence a_1,..,a_m, where a_j is the minimum number such that Gnr(j,a_j) is complete.`):


elif nargs=1 and args[1]=SeqHnr then
print(`SeqHnr(r,N): the first N terms, starting at n=r,  of the sequence enumerating the number of spanning trees of Hnr(n,r). Try:`):
print(`SeqHnr(2,30);`):


elif nargs=1 and args[1]=SeqGnr then
print(`SeGnr(r,N): the first N terms, starting at n=2*r+1,  of the sequence enumerating the number of spanning trees of Gnr(n,r). Try:`):
print(`SeqGnr(2,30);`):

elif nargs=1 and args[1]=SpTreeSeq then
print(`SpTreeSeq(F,argu,a,b): inputs a procedure F (which outputs n,E; number of vertices and edges between them) and F's list of arguments argu.`):
print(`Then, SpTreeSeq outputs a list of the number of spanning trees of the graphs F(a,op(argu)) .. F(b,op(argu)). `):
print(`Try: SpTreeSeq(HnR,[{1,3}],50,60);`):

elif nargs=1 and args[1]=GnR then
print(``):

elif nargs=1 and args[1]=HnR then
print(``):
else
print(`There is no such thing as`, args[1]):

fi:

end:

randomize():






# input graph size (vertex set) and edges between them.
# output laplacian matrix
LaplMatOld:=proc(N,G) local L,u,n,e:
 

 L:=[seq([0$N],u=1..N)]:
 for e in G do:
  u:=e[1]:
  n:=e[2]:
  L[u][u]++:
  L[n][n]++:
  L[u][n]:= -1:
  L[n][u]:= -1:
 od:
 L:
end:


# generates random N-vertex graph. returns N and edge-set. Takes in an optional second argument, p (probability) which is assumed to be rational.
GenerateGraph:=proc(N,p:=1/2) local E,i,j:
E:={}:

for i from 1 to N-1 do:
 for j from i+1 to N do:
   if rand(1..denom(p))()<=numer(p) then
    E:= E union {{i,j}}:
   fi:
 od:
od:


N,E
end:

# generates a function from [N] to {0,..,N-1}, represented as a list of length N.
GenerateFunct:=proc(N) local k,F:
F:=[-1$N]:
for k from 1 to N do:
 F[k]:=rand(0..N-1)():
od:
F:
end:

# 
TestLaplEq:= proc(N,E) local i,j,S1,S2,F,M,L,G:

 F:= GenerateFunct(N):
 L:=LaplMat(G):
 

 for i from 1 to N do:
  for j from 1 to N do:
  if member({i,j},E) then
    M[j]:=1: #indicator function for if i,j is in G.
  else:
   M[j]:=0:
  fi:
  od:

  S1:=add(L[i][j]*F[j],j=1..N):
  S2:=add(M[j]*(F[i]-F[j]),j=1..N):
  if not(S1 mod N=S2 mod N) then
   RETURN(false):
  fi:
  
 od:
 
 true:
end:

# Computer Number of Spanning Trees
with(LinearAlgebra):
NumSpanTree:=proc(n,G) local i,j,M:

M:=LaplMat(n,G):
Determinant(Matrix([seq([seq(M[i][j],j=1..nops(M)-1)]  ,i=1..nops(M)-1)])):
end:

# 
Gnr:=proc(n,r) local i,j,E,R:
 R:=r:
 if r >= n then  # I allow the user to input r >= n by simply returning G(n,n-1): 
  R:=n-1:
 fi:
 if r=0 then
  RETURN(n,{}):
 fi:


 # seq({i,i+1}, i=1..n-1) , {1,n}   <--- this code adds the cycle edges
 E:=[seq({i,i+1}, i=1..n-1) , {1,n}, seq( seq( {i+1, (i+j mod n)+1 }, j=1..R ) ,i=0..n-1)]: # this exploits the symmetry of the graph to correctly index them

 n, {op(E)}:
end:

#generalized Gnr. Inputs a subset R of [n]. Note that these graphs are not necessarily connected anymore! See GnR(8,{2,4}). We exclude loops. 
GnR:=proc(n,R) local i,j,E,Rs:
 if R={} then
  RETURN(n,{}):
 fi:
 Rs:=(R mod n) intersect {seq(i,i=1..n-1)}: # to allow loops, use: {seq(i,i=0..n-1)}.

 E:=[seq( seq( {i+1, (i+j mod n)+1 }, j in Rs) ,i=0..n-1)]: # this exploits the symmetry of the graph to correctly index them

 n, {op(E)}:
end:

# Remark: For some reason, Gnr(7,3) = Gnr(7,6). When does this occur? 

# m terms of min r s.t. G(n,r) is complete. Asume m > 0.
RSeq := proc(m) local G, S,n,r:
 S:=[]:
 for n from 1 to m do:
  r:=1:
  G:= Gnr(n,r):
  while nops(G[2]) < (n*(n-1))/2 do:
  r++:
  G:= Gnr(n,r):
 od:

  S:=[op(S),r]:

 od:
end:

# RSeq suggests that this occurs when r >= floor(n/2). 


# "Conjectures" for number of spanning trees in G(n,r) for r=1,2,3
# Given the earlier statement, this question only makes sense for n >= 2r + 2.
# r = 1. G(n,r) = C_n so number of spanning trees is n.
# r = 2. OEIS says this is n*Fib(n)^2 and there is a recent (2024) paper that relates it to spanning trees. A169630
# r = 3. OEIS also relates to spanning trees. A331905. 


# Question: If G is decomposed into cycles, C_1,..,C_m then how can you build all spanning trees of G by removing edges from these cycles?
# Such a graph G would be Eulerian.
# Answering this question would resolve DST( G(n,r) ) if the counting is straightforward.




# 
Hnr:=proc(n,r) local i,j,E:
 E:=[seq(seq({i,i+j}, i=1..n-j),j=1..r)]:

 n, {op(E)}:
end:

# Generalized. 
HnR:=proc(n,R) local i,j,E,Rs:
 Rs:= R intersect {seq(j,j=1..n)}:

 E:=[seq(seq({i,i+j}, i=1..n-j),j in Rs)]:

 n, {op(E)}:
end:

# Torus(a+b,a)
Torus:= proc(n,a) local b:
b:=n-a:
ProdGraph(Cycle(a),Cycle(b)):
end:

# Cyli(a+b,a)
Cyli:= proc(n,a) local b:
b:=n-a:
ProdGraph(Cycle(a),Path(b)):
end:

# GridGr(a+b,a)
GridGr:= proc(n,a) local b:
b:=n-a:
ProdGraph(Path(a),Path(b)):
end:

# constructs a path on n vertices
Path:=proc(n) local i:
n,{seq({i,i+1}  ,i=1..n-1)}:
end:

# constructs a cycle on n vertices
Cycle:=proc(n) local i:
n,{seq({i,i+1}  ,i=1..n-1),{1,n}}:
end:

# HypCub(2^s)
# 
HypCub:=proc(s) local i,G:
G:=Path(2):
for i from 2 to s do:
G:=ProdGraph(G, Path(2)):
od:

2^s,G[2]:
end:

# inputs two graphs n,E and m,F. Constructs their product. 
ProdGraph:=proc(n,E,m,F) local L,e,i:
L:=[seq(seq({e[1]+i*n,e[2]+i*n} ,e in E), i=0..m-1), seq(seq({(e[1]-1)*n+i,(e[2]-1)*n+i} ,e in F)  ,i=1..n)]:

n*m,{op(L)}
end:

#SeqHn(r,N): the first N terms, starting at n=r  of the sequence enumerating the number of spanning trees of Hnr(n,r). Try:
#SeqHnr(2,30);
SeqHnr:=proc(r,N) local i:
[seq(NumSpanTree(i,Hnr(i,r)[2]),i=r..N+r-1)]:
end:


#SeqGnr(r,N): the first N terms, starting at n=2r+1  of the sequence enumerating the number of spanning trees of Gnr(n,r). Try:
#SeqGnr(2,30);
SeqGnr:=proc(r,N) local i:
[seq(NumSpanTree(i,Gnr(i,r)[2]),i=2*r+1..N+2*r)]:
end:

# SpTreeSeq(F,m,argu,a,b)
# Takes in a procedure F and generates a graph by outputting n,E (number of vertices and the edge-set). The procedure must take n as its first argument.
# SpTreeSeq generates a sequence of the number of spanning trees of graphs generated by F, from n=a to n=b. SpTreeSeq also takes in the fixed arguments of F, argu, as a list.
# Try: SpTreeSeq(HnR,[{1,3}],50,60);
SpTreeSeq:=proc(F,argu,a,b) local i:

[seq(NumSpanTree( F(i, op(argu)) ), i=a..b) ]:
end:




# input graph size (vertex set) and edges between them.
# output laplacian matrix
LaplMat:=proc(N,G) local L,u,n,e,i,j:
 

for i from 1 to N do
for j from 1 to N do
L[i,j]:=0:
od:
od:

 for e in G do:
  u:=e[1]:
  n:=e[2]:
#  L[u][u]++:
 L[u,u]:=L[u,u]+1:

#  L[n][n]++:

 L[n,n]:=L[n,n]+1:

  L[u,n]:= -1:
  L[n,u]:= -1:
 od:

[seq([seq(L[i,j],j=1..N)],i=1..N)]:

end:

###start From Cfinite.txt 

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


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


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

#RGF(L,x): Guesses a rational generating function in x for the list of numbers L starting with x^0*L[1]+..., Try:
#RGF([seq(fibonacci(i),i=0..20)],x);
RGF:=proc(L,x) local gu:
gu:=GuessRec(L):
if gu=FAIL then
 RETURN(FAIL):
else
factor(CtoR(gu,x)):
fi:
end:

###End From Cfinite.txt 





#HnrS(R,x,N): inputs a positive integer R>=2 and a positive integer N outputs an article with genearting functions
#for the number of spanning trees of Hnr(n,r)
HnrS:=proc(R,x,N) local r,c,L,a,n,H,t0,f,i:


t0:=time():
print(`Generating Functions for Enumerating the Number of Spanning Trees in Friendship graphs where n people live in a one-sided STRAIGHT street, and every one is friends with all the neighbors distance at most r`):
print(`for r from 2 to`, R):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):



L:=[]:
c:=0:



for r from 2 to R  do
f:=RGF(SeqHnr(r,N),x):

if f=FAIL then
 RETURN(L):
fi:

c:=c+1:
print(``):
print(`Theorem number`, c):
print(``):
print(`Let a(n) be the number of spanning trees of the graph`, H[n,r], `whose vertices are 1, ...,n, and where every vertex i is connected to i+1,...`, min(i+r,n), ` Then`):
print(``):
print(Sum(a[n+r]*x^n,n=0..infinity)=f):
print(``):
print(`and in Maple notation`):
print(``):
lprint(f):
print(``):
L:=[op(L),f]:
od:

print(``):

print(`--------------------------------`):
print(``):
print(`This took`, time()-t0, `seconds. `):
print(``):

L:
end:


#GnrS(R,x,N): inputs a positive integer R>=2 and a positive integer N outputs an article with genearting functions
#for the number of spanning trees 
GnrS:=proc(R,x,N) local r,c,L,a,n,H,t0,f,i:

t0:=time():


print(`Generating Functions for Enumerating the Number of Spanning Trees in Friendship graphs where n people live in a one-sided CIRCULAR street, and every one is friends with all the neighbors distance at most r`):
print(`for r from 2 to`, R):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):


L:=[]:
c:=0:



for r from 2 to R  do
f:=RGF(SeqGnr(r,N),x):

if f=FAIL then
 RETURN(L):
fi:

c:=c+1:
print(``):
print(`Theorem number`, c):
print(``):
print(`Let a(n) be the number of spanning trees of the graph`, G[n,r], `whose vertices are 1, ...,n, arranged in a circle  and there is an edge between two vertices iff their distane is <=`, r, ` Then`):
print(``):
print(Sum(a[n+2*r+1]*x^n,n=0..infinity)=f):
print(``):
print(`and in Maple notation`):
print(``):
lprint(f):
print(``):
L:=[op(L),f]:
od:

print(``):

print(`--------------------------------`):
print(``):
print(`This took`, time()-t0, `seconds. `):
print(``):

L:
end: