with(Bits):
with(combinat):
with(LinearAlgebra):
with(GraphTheory):
with(RandomGraphs):


###############
# version August 21, 2026
###############

##################
##### HELP ! #####
##################
Help:=proc():
 if args=NULL then:
	print(`The main procedures are:`):
	print(`IsAlmostRegular, IsoComplex, IsAlmostRegular, LocalMaxNST, RandNST, MaxNST, MinNST, MaxNSTGraphs, MaxNSTSeq, MaxTheta, NST, SetNST, Theta`):
 elif nargs=1 and args[1]=MaxNST then:
	print(`MaxNST(n,m): returns the max number of spanning trees over all graphs with`):
	print(`n vertices and m edges`):
 elif nargs=1 and args[1]=MaxNSTGraphs then:
	print(`MaxNSTGraphs(n,m): returns the (connected) graphs satisfying max NST and`):
 	print(`which have n vertices and m edges`):
 elif nargs=1 and args[1]=MaxNSTSeq then:
	print(`MaxNSTSeq(N,c): for c fixed, computes the sequence MaxNST(n,n+c) up to n=N.`):
 elif nargs=1 and args[1]=NST then:
	print(`NST(F): Computes number of spanning trees of F. F may be a Graph object or an`):
	print(`adjacency Matrix`):
 elif nargs=1 and args[1]=Theta then:
	print(`Theta(l1,l2,l3): Given path lengths l1,l2,l3 > 1, returns the theta graph with those path`):
	print(`lengths as an expr seq (n,E), where n is number of vertices and E is the edge-set.`).
	print(`Try: Theta(2,2,2):`):
 fi:

end:

# Find the lengths l1,l2,l3 that maximize NST for a theta graph where l1+l2+l3 = N and l1,l2,l3>=2. N >= 6
MaxTheta:=proc(N) local p,tooshort,champ,num,num1,len,Q:
 p:=firstpart(N,N-3):
 tooshort := false:
 Q:=DEQueue():
 num:=0:

 while p <> FAIL do:
	len:=nops(p):	

	if len = 3 and p[1] >= 2 then:
		
		num1:=NumberOfSpanningTrees(Graph(Theta(p[1],p[2],p[3]))):
		if num1 > num then:
			#print(num1):
			Q:=DEQueue([p]):
			num:=num1:
			champ:=[p[1],p[2],p[3]]:
		elif num1 = num then:
			Q:-push_back([p]):
		fi:
	elif len < 3 then:
		return convert(Q,set):
	fi:
		
	p:=nextpart(p):

 od:

 return convert(Q,set):
end:


# This version uses SubDivEdg more cleverly to obtain a nicer labeling
Theta:=proc(l1::posint,l2::posint,l3::posint) local n,E,N,F,resu:
 option remember:
 if min(l1,l2,l3) < 2 then:
	return(FAIL):
 fi:

 N:=(l1-1)+(l2-1)+(l3-1)+2:
 
 E:=DEQueue():
 resu:=SubdivEdg(2,[{1,2}],table([{1,2} = l1-1])): # create first path
 F:=resu[2]:
 n:=resu[1]: # number of vtcs so far
 E:-push_back(op(F)):
 
 resu:=SubdivEdg(n,[{1,2}],table([{1,2} = l2-1])): # create second path
 F:=resu[2]:
 n:=resu[1]:
 E:-push_back(op(F)):

 resu:=SubdivEdg(n,[{1,2}],table([{1,2} = l3-1])): # create third path
 F:=resu[2]:
 n:=resu[1]:
 E:-push_back(op(F)):

 if n <> N then:
	error "code is bugged":
 fi:

 #print(N):
 return(n,convert(E,set)):
end:


# Generates a theta graph formed from taking two cycles of lengths m and n, which share a path on k vertices
Theta1:=proc(m,n,k) local N,E,F,T,resu:
 option remember:
 N:=m+n-k:
 if N < 5 then:
	print(`Need more vertices (m+n-k).`):
	return(FAIL):
 fi:
 E:=[{1,3},{1,4},{1,5},{2,3},{2,4},{2,5}]: # a theta graph is just a subdivided K_{2,3}
 T:=table(sparse=0,[]): 
 T[{2,3}]:=m-k-1:
 if m-k-1 < 0 then:
	return FAIL:
 fi:

 T[{2,4}]:=k-3:
 if k < 3 then:
	return FAIL:
 fi:

 T[{2,5}]:=n-k-1:
 if n-k-1 < 0 then:
	return FAIL:
 fi:

 resu:=SubdivEdg(5,E,op(T)):

 resu[1], convert(resu[2],set):
end:

# Given list L of edges (for vertices {1..n}) and a table T that maps those edges
# to the number of subdivisions, return a list F
# of the subdivided edges and N, the number of total vertices
SubdivEdg:=proc(n::posint,L::list, T::table) local newvtx,e,num,Q,i:
 Q:=DEQueue(): # this will be the output.
 newvtx:=n+1: # the next vertex that will be added
 for e in L do:
	if nops(e) <> 2 then:
		error "Error in SubdivEdg: the input list L is not a list of edges!":	
	fi:
	
	num:=T[e]:	

	# if no subdivisions, then keep edge as it is
	if num = 0 then:
		Q:-push_back(e):
	else:	
		# otherwise, start and finish off the path
		Q:-push_back({e[1], newvtx}):
		Q:-push_back({newvtx+num-1, e[2]}):
		newvtx++:
	fi:

	# form a path from e[1] to e[2]:
	for i from 2 to num do:
		Q:-push_back({newvtx-1,newvtx}):
		newvtx++:
	od:
 od:


 return(newvtx-1,convert(Q,list)):
end:


TwoCycles:=proc(n1) local n,E,i:
 n:=n1-2:
 E:={seq({i,i+1},i=1..n-1), {1,n}, {1,n+1},{n+1,n+2},{n+2, 1+floor(n/2)}}:
 return n1,E:	
end:


MinNST:=proc(n,m:=-1) local prc,out1,numT,champ:
 option remember:
 if n >= 9 then:
	print(`Warning: input is large... computation may take a while`);
 fi:


 champ:=n**(n-2):
 if m = -1 then:
 	prc:=NonIsomorphicGraphs(n,0, output = iterator, outputform = graph, restrictto = connected):	# the stable set on n vertices has the least spanning trees
 else: 
	prc:=NonIsomorphicGraphs(n,m, output = iterator, outputform = graph, restrictto = connected):
 fi:
 out1:=prc():

 while out1 <> FAIL do:
	numT:=NumberOfSpanningTrees(out1):
	if numT < champ then:
 		champ:=numT:
 	fi:
	out1:=prc():
 od: 

 champ:
end:


#####################
# 
###################
LocalMaxNST:=proc(n,m) local G,A,L:
	G:=RandomGraph(n,m,connected):
	A:=AdjacencyMatrix(G):
	while SwapEdge(A) <> FAIL do:
		L:=SwapEdge(A):
		A:=Matrix(L):
		#print(L):
	od:

	return A:
end:



#############
# Given the adjacency matrix of a graph, 
# swaps an edge to improve its number of spanning trees 
#############
SwapEdge:=proc(A) local n,m,F,E,i,j,f,e,oriNST,G,champ,B,T,L,nst,CHAMP,G2:
 E:=DEQueue():
 F:=DEQueue():
 n:=Dimensions(A)[1]:
 G:=Graph(A):
 oriNST:=true: 				# tracks whether there is a neighbor with strictly more spanning trees
 champ:=NumberOfSpanningTrees(G):
 B:=Copy(A):
 
 for i from 1 to n do:
	for j from i+1 to n do:
		if B[i,j]=1 then:
			E:-push_back([i,j]):	# edges
		elif B[i,j]=0 then:
			F:-push_back([i,j]):	# non-edges
		else:
			error "the given matrix is not an adjancecy matrix":
		fi:

	od:
 od:

 F:=convert(F,list):
 while not(empty(E)) do: 
	e:= E:-pop_back(): 	# swap edge e
	B[op(e)]:=0:
	B[e[2],e[1]]:=0:
	for f in F do:		
		B[op(f)]:=1:		# swap with a new edge
		B[f[2],f[1]]:=1:
		L:=convert(B,list,nested):	
		nst:=EdgeNST(L):
		if nst > champ then:
			champ:=nst:
			oriNST:=false:
			CHAMP:=L:
		fi:
		
		B[op(f)]:=0:
		B[f[2],f[1]]:=0:
	od:
	
	B[op(e)]:=1:		# put the edge back
	B[e[2],e[1]]:=1:
 od:

 if oriNST then:
	return FAIL:
 fi:

 G2:=Graph(Matrix(CHAMP)):
 if IsIsomorphic(G,G2) then:	# check if the best graph is actually just the same graph
	return FAIL:
 fi:
 
 return CHAMP:	# return the best neighbor
end:

####################
# remember this adjacency matrix for NST computation for later
####################
EdgeNST:=proc(L::list) local G,M:
 option remember:
  M:=Matrix(L):
  G:=Graph(M):
  NumberOfSpanningTrees(G):
end:


#########
# for c fixed, computes the sequence MaxNST(n,n+c)
########
MinNSTSeq:=proc(N,c) local n:
 [seq(MinNST(n,n+c),n=ceil(sqrt(2*c))+3..N)]:
end:


##############
# Returns the graphs with min NST for connected graphs on n vertices and m edges
##############
MinNSTGraphs:=proc(n,m:=-1) local prc,out1,numT,Q,champ:
 option remember:
 if n >= 9 then:
	print(`Warning: input is large... computation may take a while`);
 fi:


 champ:=n**(n-2):
 Q:=DEQueue();
 if m = -1 then:
 	prc:=NonIsomorphicGraphs(n,0, output = iterator, outputform = graph, restrictto = connected):	# the stable set on n vertices has the least spanning 
 else: 
	prc:=NonIsomorphicGraphs(n,m, output = iterator, outputform = graph, restrictto = connected):
 fi:
 out1:=prc():

 while out1 <> FAIL do:
	numT:=NumberOfSpanningTrees(out1):

	if numT < champ then:
 		champ:=numT:
		Q:=DEQueue():  	# reset the set of maximizers
		Q:-push_back(out1):
	elif numT = champ then:
		Q:-push_back(out1):
 	fi:

	out1:=prc():
 od: 

 convert(Q,set):
end:

#########
# for c fixed, computes the sequence MaxNST(n,n+c)
########
MaxNSTSeq:=proc(N,c) local n:
 [seq(MaxNST(n,n+c),n=ceil(sqrt(2*c))+3..N)]:
end:



###########
# Inputs a Graph object (from GraphTheory package) G and outputs true if G is Almost Regular
# and false if it is not.
##########
IsAlmostRegular:=proc(G) local S:
 S:={op(DegreeSequence(G))}:
 if S = {} then:
	error "Graph is empty.":
 fi:

 # the set of "degrees" is at most 2
 if nops(S) > 2 then:
	return(false):
 fi:

 # degrees do not differ by more than one
 if S[-1] - S[1] > 1 then:
	return(false):
 fi:

 return(true):
end:


######
UpToN:=proc(N) local k,S:
 S:={1}:
 for k from 2 to N do:
 	S:=S union kNST(k):
 od:

 S:
end:


##################
# sample labeled graphs randomly and construct the set of the number of their spanning trees
##################
RandNST:=proc(k,K) local r,A,i,Q:
 Q:=DEQueue():
 r:=rand(0..2^(binomial(k,2))-1):
 for i from 1 to K do:
 	A:=NumToGraph(k, r()): # generate random labeled graph
	Q:-push_back(NST(A)): # record how many spanning trees it has
 od:

 convert(Q,set):
end:

###############
# The set of number of spanning trees for all simple graphs on k vertices
###############
kNST:=proc(k) local Q,t,A:
 option remember:
 Q:=DEQueue():
 for t from 0 to 2^(binomial(k,2))-1 do:
 	A:=NumToGraph(k,t):			# generate a labeled graph
	Q:-push_back(NST(A)):			# compute num of span trees and add it to the list
 od:

 convert(Q,set):
end:

#################
# NumToGraph. Given positive integers n,t, convert the number t to one of the 2^(n choose 2) labeled simple graphs (represented as an adj matrix.
#################
NumToGraph:=proc(n,t) local A,L,i,j,idx,diff:
 
 if t > 2^(binomial(n,2)) then:
	print(t,2^binomial(n,2));
	return(FAIL):
 fi: 

 A:=Matrix(1..n,1..n,shape=symmetric):
 idx:=1:
 if t<>0 then:
  L:=[GetBits(t,-1..0)]:
 else:
  L:=[]:
 fi:
 diff:=binomial(n,2)-nops(L):
 L:=[0$diff,op(L)]:
 #print(L):

 # convert the binary string into a 0-1 matrix
 for i from 1 to n do:  # row coordinate
 	for j from i+1 to n do: # column coordinate
		A[i,j]:=L[idx]:
		idx++:
 	od:
 od:

 A: 
end:

#############
# Returns the number of Non-Isomorphic Connected graphs on n vertices
#############
NumGraphs:=proc(n):
 NonIsomorphicGraphs(n,restrictto=connected):
end:


##############
# Returns the graphs with max NST for connected graphs on n vertices and m edges
##############
MaxNSTGraphs:=proc(n,m:=-1) local prc,out1,numT,Q,champ:
 option remember:
 if n >= 9 then:
	print(`Warning: input is large... computation may take a while`);
 fi:


 champ:=1:
 Q:=DEQueue();
 if m = -1 then:
 	prc:=NonIsomorphicGraphs(n,(n^2-n)/2, output = iterator, outputform = graph, restrictto = connected):
 else: 
	prc:=NonIsomorphicGraphs(n,m, output = iterator, outputform = graph, restrictto = connected):
 fi:
 out1:=prc():

 while out1 <> FAIL do:
	numT:=NumberOfSpanningTrees(out1):

	if numT > champ then:
 		champ:=numT:
		Q:=DEQueue():  	# reset the set of maximizers
		Q:-push_back(out1):
	elif numT = champ then:
		Q:-push_back(out1):
 	fi:

	out1:=prc():
 od: 

 convert(Q,set):
end:

##############
# Returns the max NST for connected graphs on n vertices and m edges
##############
MaxNST:=proc(n,m:=-1) local prc,out1,numT,champ:
 option remember:
 if n >= 9 then:
	print(`Warning: input is large... computation may take a while`);
 fi:


 champ:=1:
 if m = -1 then:
 	prc:=NonIsomorphicGraphs(n,(n^2-n)/2, output = iterator, outputform = graph, restrictto = connected):	# the complete graph has the most spanning trees
 else: 
	prc:=NonIsomorphicGraphs(n,m, output = iterator, outputform = graph, restrictto = connected):
 fi:
 out1:=prc():

 while out1 <> FAIL do:
	numT:=NumberOfSpanningTrees(out1):
	if numT > champ then:
 		champ:=numT:
 	fi:
	out1:=prc():
 od: 

 champ:
end:


##############
# Returns the set of NST for connected graphs on n vertices. Optionally, can specify number of edges.
##############
SetNST:=proc(n,m:=-1) local prc,Q,out1,numT:
 option remember:
 if n >= 9 then:
	print(`Warning: input is large... computation may take a while`);
 fi:

 Q:=DEQueue():
 if m = -1 then:
 	prc:=NonIsomorphicGraphs(n, output = iterator, outputform = graph, restrictto = connected):	# use all graphs
 else: 
	prc:=NonIsomorphicGraphs(n,m, output = iterator, outputform = graph, restrictto = connected):
 fi:
 out1:=prc():

 while out1 <> FAIL do:
	numT:=NumberOfSpanningTrees(out1):
 	Q:-push_back(numT):
	out1:=prc():
 od: 

 convert(Q,set):
end:

###########
# outputs a table T where T[t] is the set of connected graphs with exactly t spanning trees. Note that all graphs are non-isomorphic
###########
IsoComplex:=proc(n) local prc,Q,out1,numT,T:
 option remember:
 T:=table(sparse={},[]);
 if n >= 9 then:
	print(`Warning: input is large... computation may take a while`);
 fi:

 prc:=NonIsomorphicGraphs(n, output = iterator, outputform = graph, restrictto = connected):
 out1:=prc():

 while out1 <> FAIL do:
	numT:=NumberOfSpanningTrees(out1):
 	T[numT]:=T[numT] union {out1}:
	out1:=prc():
 od: 

 op(T):
end:

###########
# outputs a table T where T[t] is the set of connected graphs with exactly t spanning trees. Note that all graphs are non-isomorphic
###########
RegIsoComplex:=proc(n) local prc,Q,out1,numT,T:
 option remember:
 T:=table(sparse={},[]);
 if n >= 9 then:
	print(`Warning: input is large... computation may take a while`);
 fi:

 prc:=NonIsomorphicGraphs(n, output = iterator, outputform = graph, restrictto = connected,restrictto=regular):
 out1:=prc():

 while out1 <> FAIL do:
	numT:=NumberOfSpanningTrees(out1):
 	T[numT]:=T[numT] union {out1}:
	out1:=prc():
 od: 

 op(T):
end:

######################################
### Included from MarkedGraphs.txt ###
######################################

##############
# Given an Adjacency matrix, returns the Laplacian Matrix
############
AdjToLapl:=proc(M) local M2,i,j,n:
	n:=Dimensions(M)[1]:
	M2:=Matrix(n,n):	# compute the degree matrix
	for i from 1 to n do:
		M2[i,i]:=add(M[j,i],j=1..n):
	od:
	
	M2-M:	# return the laplacian
end:


############
# Returns the number of spanning trees given an adjacency matrix
############
NST:=proc(A) local L,d: 
	if type(A,Matrix) then:
		L:=AdjToLapl(A);
		return Minor(L, 1, 1,output='determinant');
	fi:

	if type(A,function) then:
		return NumberOfSpanningTrees(A):
	fi:


	return whattype(A):
end: