with(linalg):
with(combinat):

Help:=proc():print(`CAL(d,k,L),Pro2(pi,Lim,d),Id(c,d),Pro(pi,A)`):end:

Id:=proc(c,d) 
matrix(d,d,[seq( op([0$(i-1),c,0$(d-i)])
, i=1..d)]):
end:


#Pro2(pi,Lim): Given a permutation pi
#and a list of matrices of the same size
#outputs the matrix LiM[pi[1]] LiM[pi[2]]...
Pro2:=proc(pi,LiM,d)
local i, T:

T:=matrix(d,d,[seq( op([0$(i-1),1,0$(d-i)])
, i=1..d)]):


 for i from 1 to nops(pi) do
  T:=multiply(T,LiM[pi[i]]):
  od:
T:

end:


#Pro(pi,A): Given a permutation pi
#and a list of matrices of the same size
#outputs the matrix LiM[pi[1]] LiM[pi[2]]...
Pro:=proc(pi,A)
local i, T:

T:=A[pi[1]]:

 for i from 2 to nops(pi) do
  T:=T.A[pi[i]]:
 od:
T:

end:

#CAL(d,k,L): Trying to conj. a linear dependence
#between all the k! permutations of a product
#of k matrices
CAL:=proc(d,k,L)
local A,i,Per,c,T,T1,LiM,eq,var,ra,j,ANS	:

Per:=permute(k):

ANS:=add(c[i]*Pro(Per[i],L),i=1..nops(Per)):



var:={seq(c[i],i=1..nops(Per))}:

ra:=rand(1..10):


T:=matrix(d,d,[0$(d^2)]):
eq:={}:

for j from 1 to 2*nops(var) do

for i from 1 to k do
  A[i]:=matrix(d,d,ra):
od:

LiM:=[seq(A[i],i=1..k)]:



for i from 1 to nops(Per) do


 T1:=Id(c[i],d):
  
 T1:=multiply(T1, Pro2(Per[i],LiM,d)):
  T:=evalm(T+T1):

od:

eq:=eq union{seq(seq(T[i,j],i=1..d),j=1..d)}:
od:

var:=solve(eq,var):

factor(subs(var,ANS)):


end: