######################################################################
## Larcombe.txt Save this file as  Larcombe.txt  to use it,          #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read `Lacrombe.txt` <Enter>                       #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
##  DoronZeil at gmail dot com                                       # 
######################################################################

with(linalg):

print(`First Written: July 2020: tested for Maple 2018 `):
print(`Version  : Sept. 5, 2021 [adding procedures RELn and PaperN to deal with numeric (or partially symbolic) matrices`):
print():
print(`This is Larcombe.txt, A Maple package`):
print(`accompanying Shalsoh B. Ekhad and Doron Zeilberger's  article: `):
print(` On Invariance Properties of Entries of Matrix Powers `):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/Larcombe.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():
print(`For general help, and 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: ezra1();`):
print():
 



ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are`):
print(` CheckT, CheckTi, CheckWt, CheckWtTD, DG, GenAD, GenMat, GenREL, GenTD, PeteG,  T12, T21, Tikva, , WDSji, Wt `):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` : A Maple package for  discovering invariance properties inspired by the Larcombe-Fennessey amazing invariance property`):

print(`The MAIN procedures are`):
print(`  NuLarcombeRel, Paper, PaperN, PeteTD, Ui,Ti, WDS, WDSij, REL,  RELn, WDStd, WDSijTD `):

elif nargs=1 and args[1]=CheckT then
print(`CheckT(m): Checks that T12 and T21 are inverses of each other or words in {1,2} of length m. Try:`):
print(`CheckT(6);`):

elif nargs=1 and args[1]=CheckTi then
print(`CheckTi(n,m,i): Checks that Ti(i,w) is indeed a bijection between WDSijTD(n,m,i,i+1) and WDSijTD(n,m,i+1,i). Try:`):
print(`CheckTi(4,6,1);`):

elif nargs=1 and args[1]=CheckWt then
print(`CheckWt(n,m,i,j): Check that GenMat(n,a)^m[i,j] is the weight-enumerator of WDS(n,m+1,i,j). Try:`):
print(`CheckWt(4,6,1,3);`):


elif nargs=1 and args[1]=CheckWtTD then
print(`CheckWtTD(n,m,i,j): Check that GenTD(n,a)^m[i,j] is the weight-enumerator of WDStd(n,m+1,i,j). Try:`):
print(`CheckWtTD(4,6,1,3);`):

elif nargs=1 and args[1]=DG then
print(`DG(n): all directed loopless graphs on n vertices with n edges up to isomorphism. Try:`):
print(`DG(3);`):

elif nargs=1 and args[1]=GenAD then
print(`GenAD(n,a,k): A generic n by n matrix with 2*k+1 diagonals. try:`):
print(`GenAD(5,a,2);`):

elif nargs=1 and args[1]=GenMat then
print(`GenMat(n,a): A generic n by n matrix whose entries are a[i,j]. Try:`):
print(`GenMat(3,a);`):


elif nargs=1 and args[1]=GenREL then
print(`GenREL(X,n,a,S,M): Finds all linear relations among the entries X[s] , s in S, where S is a subset of {seq([i,j],i=1..n,j=1..n) of the m-th power of a generic a[i,j] n by n matrix.`):
print(`Then it investigates whether it is still true for the powers of an n-diagonal  M by M matrix where M>m. It fails! Try:`):
print(`GenREL(X,2,a,{[1,2],[2,1]},4);`):

elif nargs=1 and args[1]=GenTD then
print(`GenTD(n,a): A generic n by n triadia-gonal matrix. Try`):
print(`GenTD(4,a);`):

elif nargs=1 and args[1]=NuLarcombeRel then
print(`NuLarcombeRel(N): The first N terms of the sequence, starting at n=1 of the sequence "number of inequivalent Larcombe relations for n by n matrices" (that equals the number of unlabeld loopless graphs`):
print(`on n vertices with exactly n edges. Try:`):
print(`NuLarcombeRel(20);`):

elif nargs=1 and args[1]=Paper then
print(`Paper(n): inputs a postive integer n>=2 and outputs a paper about all  Larcombe-Fennessey relations for a general n by n matrix. Try:`):
print(`Paper(3):`):

elif nargs=1 and args[1]=PeteG then
print(`PeteG(n,a,m,i,j): if A is a generic n by n matrix,  m a positive integer a[j,i]*A^m[i,j]-a[i,j]*A^m[j,i]. Try:`):
print(`PeteG(2,a,5,1,2);`):

elif nargs=1 and args[1]=PaperN then
print(`PaperN(a,G): inputs a matrix a and a subset of cardinality n of SUBDIAGONAL entries outputs a paper about the Larcombe relation for the matrix a`):
print(`PaperN([[1,4],[2,5]],{[1,2],[2,1]}):`):
print(`ra:=rand(1..10): a:=[seq([seq(ra(),i=1..5)],j=1..5)];PaperN(a,{[1,2],[2,3],[3,4],[4,5],[5,1]});`):

elif nargs=1 and args[1]=PeteTD then
print(`PeteTD(n,a,m,i,j): if A is a generic Tri-diagonal n by n matrix,  m a positive integer a[j,i]*A^m[i,j]-a[i,j]*A^m[j,i]. Try:`):
print(`PeteTD(2,a,5,1,2);`):

elif nargs=1 and args[1]=REL then
print(`REL(X,n,a,S): Finds all linear relations among the entries X[s] , s in S, where S is a subset of {seq([i,j],i=1..n,j=1..n) of the m-th power  (m arbitrary) of a generic a[i,j] n by n matrix. Try:`):
print(`REL(X,2,a,{[1,2],[2,1]});`):

elif nargs=1 and args[1]=RELn then
print(`RELn(X,a,S): Finds all linear relations among the entries X[s] , s in S, where S is a subset of {seq([i,j],i=1..n,j=1..n) of the n-th power of a NUMERIC n by n matrix, a. Try:`):
print(`RELn(X,[[1,3],[5,6]],{[1,2],[2,1]});`):
print(`ra:=rand(1..1000): a:=[seq([seq(ra(),i=1..5)],j=1..5)];RELn(X,a,{[1,2],[2,3],[3,4],[4,5],[5,1]});`):


elif nargs=1 and args[1]=Ui then
print(`Ui(i,w):  applies the mapping T_i to the word w . try:`):
print(`Ui(1,[2,1,3,4,3,2,1,2,3,3,2,2]);`):

elif nargs=1 and args[1]=T12 then
print(`T12(w): given a word in WDSij(2,m,1,2) for some m outputs the word in WDSji(2,m,1,2)`):
print(` 1w21^a -> 1^aw21 (w empty or starts with a 2)`):
print(`T12([1,2,1,2,1,2,2,1,1,1]);`);

elif nargs=1 and args[1]=T21 then
print(`T21(w): given a word in WDSji(2,m,1,2) for some m outputs the word in WDSij(2,m,1,2)`):
print(`1^aw21 -> 1w21^a  (w empty or starts with a 2)`):
print(`T21([1,1,1,2,1,2,2,1]);`);

elif nargs=1 and args[1]=Ti then
print(`Ti(i,w):  applies the mapping T_i to the word w . try:`):
print(`Ti(1,[1,2,3,4,3,2,1,2,3,3,2,1]);`):

elif nargs=1 and args[1]=Tikva then
print(`Tikva(n,S): Finds whether there is likely to be a relation with the directed graph S on n vertices`):
print(`Tikva(3,{[1,2],[2,1]});`):

elif nargs=1 and args[1]=WDS then
print(`WDS(n,m,i,j): The set of m-letter words in the alphabet {1,...,n} that start with a letter i and ends with a letter j. Try:`):
print(`WDS(2,4,1,2);`):


elif nargs=1 and args[1]=WDStd then
print(`WDStd(n,m,i,j): The set of m-letter words in the alphabet {1,...,n} that start with a letter i and ends with a letter j and after a letter k cacn only follow k-1,k,or k+1`):
print(`This corresponds to a tri-diagonal matrix.. Try:`):
print(`WDStd(3,4,1,2);`):


elif nargs=1 and args[1]=WDSijTD then
print(`WDSijTD(n,m,i,j): The set of m-letter words in the alphabet {1,...,n} that start and end with the letter  and after a letter k cacn only follow k-1,k,or k+1`):
print(`WDSijTD(3,4,1,2);`):


elif nargs=1 and args[1]=WDSij then
print(`WDSij(n,m,i,j): The set of m-letter words in the alphabet {1,...,n} that start and and end with the letter i AND the second letter is j. Try:`):
print(`WDSij(2,4,1,2);`):

elif nargs=1 and args[1]=WDSji then
print(`WDSji(n,m,i,j): The set of m-letter words in the alphabet {1,...,n} that start and and end with the letter i AND the penultimate letter is j. Try:`):
print(`WDSji(2,4,1,2);`):

elif nargs=1 and args[1]=Wt then
print(`Wt(w,a): The weight of the word w with respect to the symbol a. Try:`):
print(`Wt([1,2,3,1,2],a);`):

else
 print(`There is no such thing as`, args):

fi:


end:



ez:=proc(): print(` PeteTD(n,a,m,i,j) `): end:

with(combinat):



#RELn(X,a,S): Finds all linear relations among the entries X[s] , s in S, where S is a subset of {seq([i,j],i=1..n,j=1..n) of the n-th power of a NUMERIC n by n matrix, a. Try:
#RELn(X,[[1,3],[5,6]],{[1,2],[2,1]});
RELn:=proc(X,A,S) local n,i,j,b,eq,var,B,eq1,s,gu,var1,var11:
n:=nops(A):

gu:=add(b[op(s)]*X[op(s)],s in S):

var:={seq(b[op(s)],s in S)}:

eq:={}:

for i from 1 to nops(S)+2 do
B:=evalm(A&^i):
eq1:=expand(add(b[op(s)]*B[op(s)], s in S)):
eq:=eq union {eq1}:
od:

var:=solve(eq,var):



var1:={}:

for var11 in var do

if op(1,var11)=op(2,var11) then
 var1:=var1 union {op(1,var11)}:
fi:

od:

if nops(var1)<>1 then
 RETURN(FAIL):
fi:

gu:=subs(var,gu):

gu:=numer(subs(var1[1]=1,gu)):

gu:=add(factor(coeff(gu,X[op(s)],1))*X[op(s)], s in S):


for i from 1 to nops(S)+5 do
B:=evalm(A&^i):
if expand(subs({seq(X[op(s)]=B[op(s)], s in S)},gu))<>0 then
print(gu, `did not work out `):
RETURN(FAIL):
fi:
od:
gu:


end:


#Tikva1(n,S): Finds whether there is likely to be a relation with the directed graph S on n vertices
#Tikva1(3,{[1,2],[2,1]});
Tikva1:=proc(n,S) local X,a,i,j,A,b,eq,var,B,eq1,s,gu,ra:


ra:=rand(20..100):

gu:=add(b[op(s)]*X[op(s)],s in S):

var:={seq(b[op(s)],s in S)}:





A:=matrix([seq([seq(ra(),j=1..n)],i=1..n)]):

eq:={}:

for i from 1 to nops(S)+2 do
B:=evalm(A&^i):
eq1:=expand(add(b[op(s)]*B[op(s)], s in S)):
eq:=eq union {eq1}:
od:

var:=solve(eq,var):


gu:=subs(var,gu):

if gu=0 then
false:
else
true:
fi:
end:


#Tikva(n,S): Finds whether there is likely to be a relation with the directed graph S on n vertices
#Tikva(3,{[1,2],[2,1]});
Tikva:=proc(n,S) local i:
for i from 1 to 10 do
 if not Tikva1(n,S) then
  RETURN(false):
 fi:
od:
true:
end:



#DG(n): all directed loopless graphs on n vertices with n edges up to isomorphism
DG:=proc(n) local lu,U,gu,i,j,pi:

lu:=permute(n):

U:=choose({seq(seq([i,j],j=1..n),i=1..n)} minus {seq([i,i],i=1..n)},n):


gu:={}:

while U<>{} do
gu:=gu union {U[1]}:

U:=U minus {seq(subs({seq(i=pi[i],i=1..n)},U[1]),pi in lu)}:
od:

gu:

end:


#GenMat(n,a): A generic n by n matrix
GenMat:=proc(n,a)  local i,j:
matrix([seq([seq(a[i,j],j=1..n)],i=1..n)]):
end:


#PeteG(n,a,m,i,j): if A is a generic n by n matrix,  m a positive integer a[j,i]*A^m[i,j]-a[i,j]*A^m[j,i]. Try:
#PeteG(2,a,5,1,2);
PeteG:=proc(n,a,m,i,j) local A:
A:=GenMat(n,a) :
A:=evalm(A&^m):
expand(A[i,j]*a[j,i]-A[j,i]*a[i,j]):
end:


#GenTD(n,a): A generic n by n triadia-gonal matrix
GenTD:=proc(n,a)  local i:
matrix([ [a[1,1],a[1,2],0$(n-2)], seq([0$(i-2),a[i,i-1],a[i,i],a[i,i+1],0$(n-i-1)],i=2..n-1),  [0$(n-2),a[n,n-1],a[n,n]]    ]):
end:



#GenAD(n,a,k): A generic n by n matrix with 2*k+1 diagonals. try:
#GenAD(5,a,2);
GenAD:=proc(n,a,k)  local i,j,M:
M:=[seq([seq(a[i,j],j=1..n)],i=1..n)]:

for i from 1 to n do
 for j from 1 to n do
  if abs(i-j)>k then
   M:=subs(a[i,j]=0,M):
  fi:
 od:
od:
matrix(M):
end:


#PeteTD(n,a,m,i,j): if A is a generic Tri-diagonal n by n matrix,  m a positive integer a[j,i]*A^m[i,j]-a[i,j]*A^m[j,i]. Try:
#PeteTD(2,a,5,1,2);
PeteTD:=proc(n,a,m,i,j) local A,A1:


A1:=GenTD(n,a) :
A:=evalm(A1&^m):
expand(A[i,j]*A1[j,i]-A[j,i]*A1[i,j]):
end:


#WDS(n,m,i,j): The set of m-letter words in the alphabet {1,...,n} that start with a letter i and ends with a letter j. Try:
#WDS(2,4,1,2);
WDS:=proc(n,m,i,j) local gu,lu,k,lu1:
option remember:
if m<2 then
 RETURN({}):
fi:
if m=2 then
 RETURN({[i,j]}):
fi:

gu:={}:

for k from 1 to n do
 lu:=WDS(n,m-1,k,j):
 gu:=gu union {seq([i,op(lu1)],lu1 in lu)}:
od:

gu:

end:


#WDSji(n,m,i,j): The set of m-letter words in the alphabet {1,...,n} that start and and end with the letter i AND the penultimate letter is j. Try:
#WDS(2,4,1,2);
WDSji:=proc(n,m,i,j) local gu,gu1:
gu:=WDS(n,m-1,i,j):
{seq([op(gu1),i],gu1 in gu)}:
end:



#WDSij(n,m,i,j): The set of m-letter words in the alphabet {1,...,n} that start and and end with the letter i AND the second letter is j. Try:
#WDSij(2,4,1,2);
WDSij:=proc(n,m,i,j) local gu,gu1:
gu:=WDS(n,m-1,j,i):
{seq([i,op(gu1)],gu1 in gu)}:
end:


#T21(w): given a word in WDSji(2,m,1,2) for some m outputs the word in WDSij(2,m,1,2)
#1^aw21 -> 1w21^a  (w empty or starts with a 2)
T21:=proc(w) local a,m:
m:=nops(w):

if not member(w, WDSji(2,m,1,2)) then
 print(`bad input`):
 RETURN(FAIL):
fi:

for a from 1 while w[a]=1 do od:
a:=a-1:
[op(a..m-1,w),1$a]:
end:



#T12(w): given a word in WDSij(2,m,1,2) for some m outputs the word in WDSji(2,m,1,2)
# 1w21^a ->1^aw21   (w empty or starts with a 2)
T12:=proc(w) local a,m:
m:=nops(w):

if not member(w, WDSij(2,m,1,2)) then
 print(`bad input`):
 RETURN(FAIL):
fi:

for a from 1  while w[m-a+1]=1 do od:
a:=a-1:
[1$(a-1),op(1..m-a+1,w)]:
end:




#CheckT(m): Checks that T12 and T21 are inverses of each other or words in {1,2} of length m. Try:
#CheckT(6);
CheckT:=proc(m) local lu21,lu12,w:

lu21:=WDSji(2,m,1,2):
lu12:=WDSij(2,m,1,2):

if {seq(T21(w),w in lu21)}<>lu12 then
RETURN(false):
fi:


if {seq(T12(w),w in lu12)}<>lu21 then
RETURN(false):
fi:

true:

if {seq(evalb(T21(T12(w))=w), w in lu12)}<>{true} then
 RETURN(false):
fi:

if {seq(evalb(T12(T21(w))=w), w in lu21)}<>{true} then
 RETURN(false):
fi:


true:

end:

#Paper(n): inputs a postive integer n>=2 and outputs a paper about all  Larcombe-Fennessey relations for a general n by n matrix. Try:
#Paper(3):
Paper:=proc(n) local X, a, gu,gu1,gu11,M,t0,i,mu,co:

co:=0:
t0:=time():

M:=GenMat(n,a):

gu:=DG(n):

print(`There are`, nops(gu), `directed loop-less graphs with`, n , `vertices and `, n, ` edges up to isomprphism here they are`):
print(``):
lprint(gu):
print(``):
print(`For each of these there is a Larcombe-Fennessey type relations for an n by n generic matrix`):
print(``):
print(`Let M be the generic`, n, `by `, n, ` matrix `):
print(``):
print(M):

print(``):

for i from 1 to nops(gu) do
 gu1:=gu[i]:
co:=co+1:
print(``):
print(`------------------------------------`):
print(``):
print(`Theorem number`, co,`: let`, {seq(X[op(gu11)],gu11 in gu1)}, `be the entries of an ARBITRARY power of the general matrix M, we always have the linear relation `):
print(``):
mu:=REL(X,n,a,gu1):
print(mu=0):
print(``):
print(`in Maple notation `):
print(``):
lprint(mu=0):
od:

print(``):
print(`---------------------------`):
print(``):
print(`This ends this paper that took`, time()-t0, `seconds to produce.`):
end:

#Wt(w,a): The weight of the word w with respect to the symbol a. Try:
#Wt([1,2,3,1,2],a);
Wt:=proc(w,a) local j: mul(a[w[j],w[j+1]],j=1..nops(w)-1):end:



#WDStd(n,m,i,j): The set of m-letter words in the alphabet {1,...,n} that start with a letter i and ends with a letter j and after i can only come i-1 or i+1 
#correpsondint to a tri-diagonal matrix). Try:
#WDStd(2,4,1,2);
WDStd:=proc(n,m,i,j) local gu,lu,k,lu1:
option remember:
if m<2 then
 RETURN({}):
fi:

if m=2 then

if abs(i-j)<=1 then
 RETURN({[i,j]}):
else
 RETURN({}):
fi:

fi:

gu:={}:

for k from 1 to n do
 if abs(k-i)<=1 then
 lu:=WDStd(n,m-1,k,j):
 gu:=gu union {seq([i,op(lu1)],lu1 in lu)}:
 fi:
od:

gu:

end:


#CheckWt(n,m,i,j): Check that GenMat(n,a)^m[i,j] is the weight-enumerator of WDS(n,m+1,i,j). Try:
#CheckWt(4,6,1,3);
CheckWt:=proc(n,m,i,j)  local a,  M,lu,lu1:
M:=GenMat(n,a):
lu:=WDS(n,m+1,i,j):
evalb(expand(evalm(M&^m)[i,j]-add(Wt(lu1,a),lu1 in lu))=0):
end:


#CheckWtTD(n,m,i,j): Check that GenTD(n,a)^m[i,j] is the weight-enumerator of WDStd(n,m+1,i,j). Try:
#CheckWtTD(4,6,1,3);
CheckWtTD:=proc(n,m,i,j)  local a,  M,lu,lu1:
M:=GenTD(n,a):
lu:=WDStd(n,m+1,i,j):
evalb(expand(evalm(M&^m)[i,j]-add(Wt(lu1,a),lu1 in lu))=0):
end:



#WDSijTD(n,m,i,j): The subset of WDStd(n,m,i,i) whose second letter is j. Try:
#WDSijTD(2,4,1,2);
WDSijTD:=proc(n,m,i,j) local gu,gu1:
gu:=WDStd(n,m-1,j,i):
{seq([i,op(gu1)],gu1 in gu)}:
end:



#Ti(i,w):  applies the mapping T_i to the word w . try:
#Ti(1,[1,2,3,4,3,2,1,2,3,3,2,1]);
Ti:=proc(i,w) local a,w1,w2:

if nops(w)<2 then
 RETURN(FAIL):
fi:

if [w[1],w[2]]<>[i,i+1] then
 print(w, `should have started with`, [i,i+1]):
 RETURN(FAIL):
fi:

if w=[i,i+1,i] then
 RETURN([i+1,i,i+1]):
fi:

for a from 2 to nops(w)-1 while [w[a],w[a+1]]<>[i+1,i] do od:

if a=2 then
w2:=[op(4..nops(w)-1,w)]:
RETURN([i+1,i,op(w2),i,i+1]):

elif a=nops(w)-1 then
 w1:=[op(3..nops(w)-2,w)]:
 RETURN([i+1,i,i+1,op(w1),i+1]):

else
w1:=[op(3..a-1,w)]:
w2:=[op(a+2..nops(w)-1,w)]:

RETURN([i+1,i,op(w2),i,i+1,op(w1),i+1]):
fi:

end:





#Ui(i,w):  applies the mapping S_i to the word w . try:
#Ui(1,[2,1,3,4,3,2,1,2,3,3,2,2]);
Ui:=proc(i,w) local a,w1,w2:

if nops(w)<2 then
 RETURN(FAIL):
fi:

if [w[1],w[2]]<>[i+1,i] then
 print(w, `should have started with`, [i+1,i]):
 RETURN(FAIL):
fi:

if w=[i+1,i,i+1] then
 RETURN([i,i+1,i]):
fi:

for a from nops(w)-1 by -1 to 2  while [w[a],w[a+1]]<>[i,i+1] do od:

if a=nops(w)-1 then
w2:=[op(3..nops(w)-2,w)]:
RETURN([i,i+1,i,op(w2),i]):

elif a=2 then
 w1:=[op(4..nops(w)-1,w)]:
 RETURN([i,i+1,op(w1),i+1,i]):

else
w1:=[op(a+2..nops(w)-1,w)]:
w2:=[op(3..a-1,w)]:

RETURN([i,i+1,op(w1),i+1,i,op(w2),i]):
fi:

end:


#CheckTi(n,m,i): Checks that Ti(i,w) is indeed a bijection between WDSijTD(n,m,i,i+1) and WDSijTD(n,m,i+1,i). Try:
#CheckTi(4,6,1);
CheckTi:=proc(n,m,i) local gu,mu,gu1,mu1,w:

if not (type(n, integer) and type(m,integer) and n>=2 and m>=3 and i>=1 and i<n) then
 print(`bad input`):
 RETURN(FAIL):
fi:


gu:=WDSijTD(n,m,i,i+1):

mu:=WDSijTD(n,m,i+1,i):

evalb({seq(evalb(Ui(i,Ti(i,w))=w), w in gu)}={true}) and evalb({seq(evalb(Ti(i,Ui(i,w))=w), w in mu)}={true}) and
evalb({seq(Ti(i,gu1), gu1 in gu)}=mu) and evalb({seq(Ui(i,mu1), mu1 in mu)}=gu):


end:



#GenREL(X,n,a,S,M): Finds all linear relations among the entries X[s] , s in S, where S is a subset of {seq([i,j],i=1..n,j=1..n) of the m-th power of a generic a[i,j] n by n matrix.
#Then it investigates whether it is still true for the powers of an n-diagonal  M by M matrix where M>m. Try:
#GenREL(X,2,a,{[1,2],[2,1]},4);

GenREL:=proc(X,n,a,S,M) local gu,A,k,Ak,gu1,s:

if M<=n then
 print(M, `should be larger than n`):
 RETURN(FAIL):
fi:

gu:=REL(X,n,a,S):

A:=GenAD(M,a,n-1):


for k from 1 to M+2 do

print(`k is`, k):

 Ak:=evalm(A&^k):

print(expand(subs({seq(X[op(s)]=Ak[op(s)], s in S)},gu))):

od:


gu:

end:









###Begin graph counting
ez:=proc(): print(` PeterS(n,x), Peter(n,x), PeterA(n,x), IndPer(pi), Conve(P), IndPerG(pi), PtoF(p,n):  `):end:

with(combinat):

##start general procedures
#IndPer(pi): The induced permutation of length n*(n-1) on directed loop-less graphs
IndPer:=proc(pi) local n,T,i,j,S,U:

n:=nops(pi):

S:=[seq(op([seq([i,j], j=1..i-1),seq([i,j], j=i+1..n)]),i=1..n)]:

for i from 1 to nops(S) do
U[S[i]]:=i:
T[S[i]]:=[pi[S[i][1]],pi[S[i][2]]]:
od:

[seq(U[T[S[i]]],i=1..nops(S))]:

end:

#GrabCycle(P,i): inputs a permutation P of {1, ...,n} and outputs the cycle belonging to i. 
#We use the convention that it starts with the largest entry
#Try
#GrabCycle([3,1,4,2],2);
GrabCycle:=proc(P,i) local C,j:

#We view the cycle belnging to i as a list, starting at i
C:=[i]:

#We find where i points to, call it j
j:=P[i]:

while j<>i do
#sooner or later we will wind back at i (since the permutation is finite) so we append the new
#arrivals to the list C and rename j
C:=[op(C),j]:
j:=P[j]:
od:

#We look at the place where it it is max
j:=max[index](C):

#we rewrite the cycle with the largest entry first

[op(j..nops(C),C),op(1..j-1,C)]:

end:

#PtoC(P): Inputs a permutation P outputs its list of cycles. Try:
#PtoC([2,1,3,4]);
PtoC:=proc(P) local n,StillToDo,L,i,C,L1,T:

n:=nops(P):

StillToDo:={seq(i,i=1..n)}:

#L is a dynamically construcete list of cycles in the permutatib P, we start with the empty list
#until no one is left
L:=[]:

while StillToDo<>{} do
#we pick the smallest survivor
 i:=min(op(StillToDo)):
 #we grab its cycle
 C:=GrabCycle(P,i):

# We append C to L
 L:=[op(L),C]:

#We update StillToDo by kicking out the members of C

StillToDo:= StillToDo minus convert(C,set):
od:

#So far L is a list of cycles but arranged in a random order
#Now we arrange them according to the convention that the first (largest) entries are increasing
#L1 is the list of first (largest) entries of each cycle:
L1:=[seq(L[i][1],i=1..nops(L))]:

#We now make a table where T[a] is the cycle that starts with  a
for i from 1 to nops(L) do
 T[L[i][1]]:=L[i]:
od:

#Now we sort L1
L1:=sort(L1):

#Now we rewrite L with the above convention
[seq(T[L1[i]],i=1..nops(L1))]:

end:

#CtoP(C): Inputs a permutation in cycle structre , outputs it in 1-line notation
#Try([[1],[3,2]]);
CtoP:=proc(C) local i,L,n,T,C1,j:
L:=[seq(op(C[i]),i=1..nops(C))]:
n:=nops(L):

if convert(L,set)<>{seq(i,i=1..n)} then
  RETURN(FAIL):
fi:

for i from 1 to nops(C) do
 C1:=C[i]:
 for j from 1 to nops(C1)-1 do
   T[C1[j]]:=C1[j+1]:
 od:
 T[C1[nops(C1)]]:=C1[1]:
od:

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


#PtoF(p,n):converts the partition p of n to frequency notation. Try
#PtoF([3,2],5);
PtoF:=proc(p,n) local x,gu,i:
gu:=mul(x[p[i]],i=1..nops(p)):
[seq(degree(gu,x[i]),i=1..n)]:
end:


##end general procedures

####start loopless directed graphs

#Conve(P): converts the cycle-structre of a permutation pi on the vertices of a loopless directed graph to those of the edges. Try:
#Conve([3]);
Conve:=proc(P) local i,co,C,i1:

co:=1:

C:=[]:

for i from 1 to nops(P) do
C:=[op(C),[seq(co+i1,i1=1..P[i]-1),co]]:
co:=co+P[i]:

od:

C:=PtoC(IndPer(CtoP(C))):

sort([seq(nops(C[i]),i=1..nops(C))],`>`):

end:



#Peter(n,x): The weight-enumerator, according to the number of edges, of  directed loopless graphs graphs on n vertices. Try:
#Peter(5,x);
Peter:=proc(n,x) local mu,i,gu,k1,j1,lu,ja:

mu:=partition(n):
gu:=0:

for i from 1 to nops(mu) do
ja:=PtoF(mu[i],n):
lu:=mul((1/k1)^ja[k1]/ja[k1]!,k1=1..nops(ja)):

j1:=Conve(mu[i]):

gu:=expand(gu+lu*mul(1+x^j1[k1],k1=1..nops(j1))):

od:

gu:

end:


####end loopless directed graphs

####start simple graphs graphs

#IndPerS(pi): The induced permutation of length n*(n-1) on directed loop-less graphs
IndPerS:=proc(pi) local n,T,i,j,S,U:

n:=nops(pi):

S:=[seq(op([seq({i,j}, j=i+1..n)]),i=1..n)]:

for i from 1 to nops(S) do
U[S[i]]:=i:
T[S[i]]:={pi[S[i][1]],pi[S[i][2]]}:
od:

[seq(U[T[S[i]]],i=1..nops(S))]:

end:




#ConveS(P): converts the cycle-structre of a permutation pi on the vertices of a simple graph to those of the edges. Try:
#ConveS([3]);
ConveS:=proc(P) local i,co,C,i1:

co:=1:

C:=[]:

for i from 1 to nops(P) do
C:=[op(C),[seq(co+i1,i1=1..P[i]-1),co]]:
co:=co+P[i]:

od:

C:=PtoC(IndPerS(CtoP(C))):

sort([seq(nops(C[i]),i=1..nops(C))],`>`):

end:




#PeterS(n,x): The weight-enumerator, according to the number of edges, of  simple graphs on n vertices. Try:
#PeterS(5,x);
PeterS:=proc(n,x) local mu,i,gu,k1,j1,lu,ja:

mu:=partition(n):
gu:=0:

for i from 1 to nops(mu) do
ja:=PtoF(mu[i],n):
lu:=mul((1/k1)^ja[k1]/ja[k1]!,k1=1..nops(ja)):

j1:=ConveS(mu[i]):

gu:=expand(gu+lu*mul(1+x^j1[k1],k1=1..nops(j1))):

od:

gu:

end:

####end simple graphs 


####start all directed graphs including loops






#IndPerA(pi): The induced permutation of length n*(n-1) on directed loop-less graphs
IndPerA:=proc(pi) local n,T,i,j,S,U:

n:=nops(pi):

S:=[seq(op([seq([i,j], j=1..n)]),i=1..n)]:

for i from 1 to nops(S) do
U[S[i]]:=i:
T[S[i]]:=[pi[S[i][1]],pi[S[i][2]]]:
od:

[seq(U[T[S[i]]],i=1..nops(S))]:

end:

#ConveA(P): converts the cycle-structre of a permutation pi on the vertices of a loopless directed graph to those of the edges. Try:
#ConveA([3]);
ConveA:=proc(P) local i,co,C,i1:

co:=1:

C:=[]:

for i from 1 to nops(P) do
C:=[op(C),[seq(co+i1,i1=1..P[i]-1),co]]:
co:=co+P[i]:

od:

C:=PtoC(IndPerA(CtoP(C))):

sort([seq(nops(C[i]),i=1..nops(C))],`>`):

end:



#PeterA(n,x): The weight-enumerator, according to the number of edges, of  directed  graphs (with loops) on n vertices. Try:
#PeterA(5,x);
PeterA:=proc(n,x) local mu,i,gu,k1,j1,lu,ja:

mu:=partition(n):
gu:=0:

for i from 1 to nops(mu) do
ja:=PtoF(mu[i],n):
lu:=mul((1/k1)^ja[k1]/ja[k1]!,k1=1..nops(ja)):

j1:=ConveA(mu[i]):

gu:=expand(gu+lu*mul(1+x^j1[k1],k1=1..nops(j1))):

od:

gu:

end:


#NuLarcombeRel(N): The first N terms of the sequence, starting at n=1 of the sequence "number of inequivalent Larcombe relations for n by n matrices" (that equals the number of unlabeld loopless graphs
#on n vertices with exactly n edges. Try:
#NuLarcombeRel(20);
NuLarcombeRel:=proc(N) local x,n:
[seq(coeff(Peter(n,x),x,n),n=1..N)]:
end:

####end all  directed graphs, including loops

###End graph counting



#PaperN(a,G): inputs a matrix a and a subset of cardinality n of SUBDIAGONAL entries outputs a paper about the Larcombe relation for the matrix a
#PaperN([[1,4],[2,5]],{[1,2],[2,1]}):
PaperN:=proc(M,G) local X, n,mu,A,m,i,j:
n:=nops(M):
if not (type(M,list) and type(G,set) and nops(G)=nops(M) and G minus { seq(op([seq([i,j],j=1..i-1),seq([i,j],j=i+1..n)]),i=1..n)}={} ) then
 print(`bad input`):
RETURN(FAIL):
fi:


print(``):
print(`------------------------------------`):
print(``):
print(`Fact: For the following set of entries`, G ):
print(``):
print(`and the matrix A:=`):
print(``):
print(matrix(M)):
print(``):

mu:=RELn(X,M,G):

print(`Then for every integer m we have the relation`):
print(subs(X=A^m,mu)=0):
print(``):
print(`in Maple notation `):
print(``):
lprint(subs(X=A^m,mu)=0):
print(``):

print(``):
print(`---------------------------`):
print(``):

end: