######################################################################
## AL.txt Save this file as  AL.txt  to use it,                      #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read `AL.txt` <Enter>                             #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Joseph Koutsoutis as project in Dr. Z.'s class         #
# (Doron Zeilberger), Rutgers University , Spring 2025               #
######################################################################


print(`Version  : May 7, 2025  `):
print():
print(`Written by Joseph Koutsoutis as a project in the Experimental Math class, taught by Dr. Z. (Doron Zeilberger), Rutgers University , Spring 2025`):


print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/EM25/AL.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 "Help();". For specific help type "Help(procedure_name);" `):
print(`For a list of the supporting functions type: Help1();`):
print():
 
Digits:=20:


Help1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are:`):
print(`Mul, Pow, LC, Vk`):
print(`RGM, RGMs, RSM, RSMs, RASM, RASMs,`):
print(`RTDM, RTDMs, RUTM, RUTMs, RToeplitzM, RToeplitzMs,`):
print(`RHankelM, RHankelMs, RCM, RCMs, RACM, RACMs,`):
print(`RCGM, RCGMs, RTDToeplitzMs, RTDToeplitzMs,`):
print(`RTDHankelM, RTDHankelMs`):
print(`TestAC, TestCG0,`):

else
Help(args):

fi:


end:

Help:=proc()
if args=NULL then

print(` A Maple package for discovering Amistur-Levitsky type identities `):

print(`The MAIN procedures are:`):
print(`ALnm,`):
print(`ALnmG, ALnmS, ALnmAS, ALnmTD, ALnmUt,`):
print(`ALnmToeplitz, ALnmHankel, ALnmC, ALnmAC, ALnmCG,`):
print(`ALnmTDToeplitz, ALnmTDHankel,`):


elif nargs=1 and args[1]=ALnm then
print(`ALnm(X,n,m,K,matrix_gen,k,S,pow): Tries to find a relation between m n by n matrices generated using matrix_gen.`):
print(`matrix_gen takes as input n,K,m where K is usually used so that the elements of the matrices are randomly chosen in [-K,K].`):
print(`matrix_gen is run ceil(m!/n^2) + k times, where k is an optional parameter with default value 0`):
print(`S is also an optional parameter that zeros out the entries of the generated matrices in the positions specified by S. Its default value is {}`):
print(`pow is a third optional parameter where we raise the matrices in the sum to the powers specified. Its default value is [1$m]`):
print(`Try: ALnm(X,1,2,10,RGMs); ALnm(X,2,4,10,RGMs);`):
print(`This will take a long time:  ALnm(X,3,6,5,RGMs);`):

elif nargs=1 and args[1]=RGM then
print(`RGM(n,K): a random n by n matrix with entries in [1,K]. Try:`):
print(`RGM(5,100);`):

elif nargs=1 and args[1]=RGMs then
print(`RGMs(n,K,m): A list of m random n by n matrices with entries from [1,K]. Try: `):
print(`RGMs(3,100,3);`):

elif nargs=1 and args[1]=ALnmG then
print(`ALnmG(X,n,m,K,k,S,pow): Tries to find a relation between m n by n matrices.`):
print(`k,S,pow are optional parameters (see Help(ALnm) for their description). Try:`):
print(`ALnmG(X,3,5,10,30,{[1,1],[1,2],[1,3]});`):

elif nargs=1 and args[1]=RSM then
print(`RSM(n,K): a random symmetric n by n matrix with entries in [1,K]. Try:`):
print(`RSM(5,100);`):

elif nargs=1 and args[1]=RSMs then
print(`RSMs(n,K,m): A list of m random symmetric n by n matrices with entries from [1,K]. Try: `):
print(`RSMs(3,100,3);`):

elif nargs=1 and args[1]=ALnmS then
print(`ALnmS(X,n,m,K,k,S,pow): Tries to find a relation between m n by n symmetric matrices.`):
print(`k,S,pow are optional parameters (see Help(ALnm) for their description). Try:`):
print(`ALnmS(X,3,3,10,30,{[1,1],[2,2],[3,3]});`):

elif nargs=1 and args[1]=RASM then
print(`RASM(n,K): a random antisymmetric n by n matrix with entries in [-K,K]. Try:`):
print(`RASM(5,100);`):

elif nargs=1 and args[1]=RASMs then
print(`RASMs(n,K,m): A list of m random antisymmetric n by n matrices with entries from [-K,K]. Try: `):
print(`RASMs(3,100,3);`):

elif nargs=1 and args[1]=ALnmAS then
print(`ALnmAS(X,n,m,K,k,S,pow): Tries to find a relation between m n by n antisymmetric matrices.`):
print(`k,S,pow are optional parameters (see Help(ALnm) for their description). Try:`):
print(`ALnmAS(X,3,4,10,30,{[1,1],[2,2],[3,3]}); ALnmAS(X,4,5,10,30,{[1,1],[2,2],[3,3],[4,4]});`):

elif nargs=1 and args[1]=RTDM then
print(`RTDM(n,K): a random tridiagonal n by n matrix with entries in [0,K]. Try:`):
print(`RTDM(5,100);`):

elif nargs=1 and args[1]=RTDMs then
print(`RTDMs(n,K,m): A list of m random tridiagonal n by n matrices with entries from [0,K]. Try: `):
print(`RTDMs(3,100,3);`):

elif nargs=1 and args[1]=ALnmTD then
print(`ALnmTD(X,n,m,K,k,S,pow): Tries to find a relation between m n by n tridiagonal matrices.`):
print(`k,S,pow are optional parameters (see Help(ALnm) for their description). Try:`):
print(`ALnmTD(X,3,5,10,30,{[1,1],[2,2],[3,3]});`):

elif nargs=1 and args[1]=RUTM then
print(`RUTM(n,K): a random uppertriangular n by n matrix with entries in [0,K]. Try:`):
print(`RUTM(5,100);`):

elif nargs=1 and args[1]=RUTMs then
print(`RUTMs(n,K,m): A list of m random uppertriangular n by n matrices with entries from [0,K]. Try: `):
print(`RUTMs(3,100,3);`):

elif nargs=1 and args[1]=ALnmUT then
print(`ALnmUT(X,n,m,K,k,S,pow): Tries to find a relation between m n by n uppertriangular matrices.`):
print(`k,S,pow are optional parameters (see Help(ALnm) for their description). Try:`):
print(`ALnmUT(X,3,3,10,30,{[1,1],[2,2],[3,3]}); ALnmUT(X,4,4,10,30,{[1,1],[2,2],[3,3],[4,4]}); ALnmUT(X,5,5,10,30,{[1,1],[2,2],[3,3],[4,4],[5,5]});`):

elif nargs=1 and args[1]=RToeplitzM then
print(`RToeplitzM(n,K): a random Toeplitz n by n matrix with entries in [1,K]. Try:`):
print(`RToeplitzM(5,100);`):

elif nargs=1 and args[1]=RToeplitzMs then
print(`RToeplitzMs(n,K,m): A list of m random Toeplitz n by n matrices with entries from [1,K]. Try: `):
print(`RToeplitzMs(3,100,3);`):

elif nargs=1 and args[1]=ALnmToeplitz then
print(`ALnmToeplitz(X,n,m,K,k,S,pow): Tries to find a relation between m n by n Toeplitz matrices.`):
print(`k,S,pow are optional parameters (see Help(ALnm) for their description). Try:`):
print(`ALnmToeplitz(X, 3, 2, 10, 30, {[2,1],[3,1],[3,2]}); ALnmToeplitz(X, 4, 2, 10, 30, {[2,1],[3,1],[3,2],[4,1],[4,2],[4,3]});`):

elif nargs=1 and args[1]=RHankelM then
print(`RHankelM(n,K): a random Hankel n by n matrix with entries in [1,K]. Try:`):
print(`RHankelM(5,100);`):

elif nargs=1 and args[1]=RHankelMs then
print(`RHankelMs(n,K,m): A list of m random Hankel n by n matrices with entries from [1,K]. Try: `):
print(`RHankelMs(3,100,3);`):

elif nargs=1 and args[1]=ALnmHankel then
print(`AALnmHankel(X,n,m,K,k,S,pow): Tries to find a relation between m n by n Hankel matrices.`):
print(`k,S,pow are optional parameters (see Help(ALnm) for their description). Try:`):
print(`ALnmHankel(X, 3, 4, 10, 30, {[3,2],[3,3],[2,3]}); ALnmHankel(X, 4, 5, 10, 30, {[4,4],[3,4],[2,4],[4,3],[3,3],[4,2]}); ALnmHankel(X, 5, 6, 10, 30, {[5,2],[5,3],[5,4],[5,5],[4,3],[4,4],[4,5],[3,4],[3,5],[2,5]});`):

elif nargs=1 and args[1]=RCM then
print(`RCM(n,K): a random circulant n by n matrix with entries in [1,K]. Try:`):
print(`RCM(5,100);`):

elif nargs=1 and args[1]=RCMs then
print(`RCMs(n,K,m): A list of m random circulant n by n matrices with entries from [1,K]. Try: `):
print(`RCMs(3,100,3);`):

elif nargs=1 and args[1]=ALnmC then
print(`ALnmC(X,n,m,K,k,S,pow): Tries to find a relation between m n by n circulant matrices.`):
print(`k,S,pow are optional parameters (see Help(ALnm) for their description). Try:`):
print(`ALnmC(X,3,2,10,30); ALnmC(X,4,2,10,30);`):

elif nargs=1 and args[1]=RACM then
print(`RACM(n,K): a random anticirculant n by n matrix with entries in [1,K]. Try:`):
print(`RACM(5,100);`):

elif nargs=1 and args[1]=RACMs then
print(`RACMs(n,K,m): A list of m random anticirculant n by n matrices with entries from [1,K]. Try: `):
print(`RACMs(3,100,3);`):

elif nargs=1 and args[1]=ALnmAC then
print(`ALnmAC(X,n,m,K,k,S,pow): Tries to find a relation between m n by n anticirculant matrices.`):
print(`k,S,pow are optional parameters (see Help(ALnm) for their description). Try:`):
print(`ALnmAC(X,4,3,10,30); ALnmC(X,5,3,10,30);`):

elif nargs=1 and args[1]=RCGM then
print(`RCGM(n,K,s): a random circulant n by n matrix with shift s with entries in [1,K]. Try:`):
print(`RCGM(5,100,2);`):

elif nargs=1 and args[1]=RCGMs then
print(`RCGMs(n,K,m,s): A list of m random circulant n by n matrices with shift s with entries from [1,K]. Try: `):
print(`RCGMs(3,100,3,2);`):

elif nargs=1 and args[1]=ALnmCG then
print(`ALnmCG(X,n,m,K,s,k,S,pow): Tries to find a relation between m n by n circulant matrices with shift s.`):
print(`k,S,pow are optional parameters (see Help(ALnm) for their description). Try:`):
print(`ALnmCG(X,4,3,10,0,30); ALnmCG(X,5,3,10,0,30);`):

elif nargs=1 and args[1]=RTDToeplitzM then
print(`RTDToeplitzM(n,K): a random tridiagonal Toeplitz n by n matrix with entries in [0,K]. Try:`):
print(`RTDToeplitzM(5,100);`):

elif nargs=1 and args[1]=RTDToeplitzMs then
print(`RTDToeplitzMs(n,K,m): A list of m random tridiagonal Toeplitz n by n matrices with entries from [0,K]. Try: `):
print(`RTDToeplitzMs(3,100,3);`):

elif nargs=1 and args[1]=ALnmTDToeplitz then
print(`ALnmTDToeplitz(X,n,m,K,k,S,pow): Tries to find a relation between m n by n tridiagonal Toeplitz matrices.`):
print(`k,S,pow are optional parameters (see Help(ALnm) for their description). Try:`):
print(`ALnmTDToeplitz(X, 4, 4, 10, 30); ALnmTDToeplitz(X, 5, 4, 10, 30);`):

elif nargs=1 and args[1]=RTDHankelM then
print(`RTDHankelM(n,K): a random tri(skew)diagonal Hankel n by n matrix with entries in [0,K]. Try:`):
print(`RTDHankelM(5,100);`):

elif nargs=1 and args[1]=RTDHankelMs then
print(`RTDHankelMs(n,K,m): A list of m random tri(skew)diagonal Hankel n by n matrices with entries from [0,K]. Try: `):
print(`RTDHankelMs(3,100,3);`):

elif nargs=1 and args[1]=ALnmTDHankel then
print(`ALnmTDHankel(X,n,m,K,k,S,pow): Tries to find a relation between m n by n tri(skew)diagonal Hankel matrices.`):
print(`k,S,pow are optional parameters (see Help(ALnm) for their description). Try:`):
print(`ALnmTDHankel(X, 4, 4, 10, 30); ALnmTDHankel(X, 5, 4, 10, 30);`):

elif nargs=1 and args[1]=Mul then
print(`Mul(A,B): the product of matrix A and B (assuming that it exists). Try:`):
print(`Mul([[a1,a2]],[[b1],[b2]]);`):

elif nargs=1 and args[1]=Pow then
print(`Pow(M,k): the matrix M^k. Try:`):
print(`Pow([[a1,a2],[a3,a4]], 2);`):

elif nargs=1 and args[1]=LC then
print(`LC(p): returns true with probability p and false otherwise. Try:`):
print(`add(LC(1/4), i=1..100);`):

elif nargs=1 and args[1]=Vk then
print(`Vk(k, K): returns a set of all length k lists with entries in [K]. Try:`):
print(`Vk(5, 2)`):

elif nargs=1 and args[1]=TestAC then
print(`TestAC(n): tests if 3 n by n anticirculant matrices satisfy the identity suggested by ALnmAC(X,4,3,10,30).`):
print(`TestAC(50)`):

elif nargs=1 and args[1]=TestCG0 then
print(`TestCG0(n): tests if 3 n by n circulant matrices with shift 0 satisfy the identity suggested by ALnmCG(X,4,3,10,0,30). Try:`):
print(`TestCG0(50)`):

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

fi:


end:

with(combinat):

#Mul(A,B): the product of matrix A and B (assuming that it exists)
Mul:=proc(A,B) local i,j,k:
[seq([seq(add(A[i][k]*B[k][j],k=1..nops(A[i])),j=1..nops(B[1]))],i=1..nops(A))]:
end:


#Pow(M,k): the matrix M^k
Pow:=proc(M,k) local i,M1:
  M1:=M:
  for i from 2 to k do
    M1:=Mul(M1,M):
  od:
  M1:
end:


#LC(p): returns true with probability p and false otherwise
LC:=proc(p) local a,b,ra:
  a:=numer(p):
  b:=denom(p):
  ra:=rand(1..b)():
  if ra<=a then
    true:
  else
    false:
  fi:
end:


#ALnm(X,n,m,K,matrix_gen,k,S,pow): Tries to find a relation between m n by n matrices generated using matrix_gen
#matrix_gen takes as input n,K,m where K is usually used so that the elements of the matrices are randomly chosen in [-K,K]
#matrix_gen is run ceil(m!/n^2) + k times, where k is an optional parameter with default value 0
#S is also an optional parameter that zeros out the entries of the generated matrices in the positions specified by S. Its default value is {}
#pow is a third optional parameter where we raise the matrices in the sum to the powers specified. Its default value is [1$m]
ALnm:=proc(X,n,m,K,matrix_gen,k:=0,S:={},pow:=[]) local A, pi, eqs, M, var, i, var1, v, JK,M1,c,i1, pow1, s:
  if pow = [] then:
    pow1 := [1$m]:
  else:
    pow1 := pow:
  fi:
  var := {seq(c[pi], pi in permute(m))}:
  JK:=add(c[pi]*X[pi],pi in permute(m)):

  eqs := {}:

  for i from 1 to ceil(m!/n^2) + k do
    A := matrix_gen(n,K,m): 
    for s in S do: 
      for i1 from 1 to m do:
        A[i1][s[1]][s[2]] := 0:
      od:
    od:
    M:=[[0$n]$n]: 
    for pi in permute(m) do
      M1:=Pow(A[pi[1]], pow1[1]):
      for i1 from 2 to m do
        M1:=Mul(M1, Pow(A[pi[i1]], pow1[i1])):
      od:

      M:=expand(M+c[pi]*M1):
    od:

    eqs := eqs union {seq(op(M[i1]),i1=1..nops(M))}:
  od:

  var := solve(eqs, var):

  if var=NULL then
    RETURN(FAIL):
  fi:

  if subs(var,JK)=0 then
    RETURN(0):
  fi: 

  var1:={}:
  for v in var do
    if op(1,v)=op(2,v) then
      var1:=var1 union {op(1,v)}:
    fi:
  od:

  subs({seq(var1[i]=1,i=1..nops(var1))},subs(var,JK)):
end:


#RGM(n,K): a random n by n matrix with entries in [1,K]
RGM:=proc(n,K) local ra,i,j:
  ra := rand(1..K):
  [seq([seq(ra(),i=1..n)],j=1..n)]:
end:


#RGMs(n,K,m): A list of m random n by n matrices with entries from [1,K]
RGMs:=proc(n,K,m) local i:
  [seq(RGM(n,K), i=1..m)]:
end:


#ALnmG(X,n,m,K,k,S): Tries to find a relation between m n by n matrices with entries
ALnmG:=proc(X,n,m,K,k:=20,S:={},pow:=[]):
  ALnm(X,n,m,K,RGMs,k,S,pow):
end:


#RSM(n,K): a random symmetric n by n matrix with entries in [1,K]
RSM:=proc(n,K) local ra,i,j,M,c:
  ra := rand(1..K):
  for i from 1 to n do:
    for j from i to n do:
      c := ra():
      M[i,j] := c:
      M[j,i] := c:
    od:
  od:
  [seq([seq(M[i,j],i=1..n)],j=1..n)]:
end:


#RSMs(n,K,m): A list of m s random n by n symmetric matrices with entries from [1,K]
RSMs:=proc(n,K,m) local i:
  [seq(RSM(n,K), i=1..m)]:
end:


#ALnmS(X,n,m,K): Tries to find a relation between m n by n symmetric matrices
ALnmS:=proc(X,n,m,K,k:=0,S:={},pow:=[]):
  ALnm(X,n,m,K,RSMs,k,S,pow):
end:


#RASM(n,K): a random antisymmetric n by n matrix with entries in [-K,K]
RASM:=proc(n,K) local ra,i,j,M,c:
  ra := rand(1..K):
  for i from 1 to n do:
    for j from i to n do:
      c := ra():
      if LC(1/2) then:
        M[i,j] := c:
        M[j,i] := -1*c:
      else:
        M[i,j] := -1*c:
        M[j,i] := c:
      fi:
    od:
  od:
  [seq([seq(M[i,j],i=1..n)],j=1..n)]:
end:


#RASMs(n,K,m): A list of m s random n by n antisymmetric matrices with entries from [-K,K]
RASMs:=proc(n,K,m) local i:
  [seq(RASM(n,K), i=1..m)]:
end:


#ALnmAS(X,n,m,K): Tries to find a relation between m n by n antisymmetric matrices
ALnmAS:=proc(X,n,m,K,k:=0,S:={},pow:=[]):
  ALnm(X,n,m,K,RASMs,k,S,pow):
end:


#RTDM(n,K): a random tridiagonal n by n matrix with entries in [0,K]
RTDM:=proc(n,K) local ra,i,j,M:
  ra := rand(1..K):
  for i from 1 to n do:
    for j from 1 to n do:
      if j >= i-1 and j <= i+1 then:
        M[i,j] := ra():
      else:
        M[i,j] := 0:
      fi:
    od:
  od:
  [seq([seq(M[i,j],i=1..n)],j=1..n)]:
end:


#RTDMs(n,K,m): A list of m random n by n tridiagonal matrices with entries in [0,K]
RTDMs:=proc(n,K,m) local i:
  [seq(RTDM(n,K), i=1..m)]:
end:


#ALnmTD(X,n,m,K): Tries to find a relation between m n by n tridiagonal matrices
ALnmTD:=proc(X,n,m,K,k:=0,S:={},pow:=[]):
  ALnm(X,n,m,K,RTDMs,k,S,pow):
end:


#RUTM(n,K): a random uppertriangular n by n matrix with entries in [0,K]
RUTM:=proc(n,K) local ra,i,j,M:
  ra := rand(1..K):
  for i from 1 to n do:
    for j from 1 to n do:
      if j <= i then:
        M[i,j] := ra():
      else:
        M[i,j] := 0:
      fi:
    od:
  od:
  [seq([seq(M[i,j],i=1..n)],j=1..n)]:
end:


#RUTMs(n,K,m): A list of m s random n by n uppertriangular matrices with entries in [0,K]
RUTMs:=proc(n,K,m) local i:
  [seq(RUTM(n,K), i=1..m)]:
end:


#ALnmUT(X,n,m,K): Tries to find a relation between m n by n uppertriangular matrices
ALnmUT:=proc(X,n,m,K,k:=0,S:={},pow:=[]):
  ALnm(X,n,m,K,RUTMs,k,S,pow):
end:


#RToeplitzM(n,K): a random Toeplitz n by n matrix with entries in [1,K]
RToeplitzM:=proc(n,K) local ra,i,j,M,c:
  ra := rand(1..K):
  for i from 1 to n do:
    c := ra():
    for j from 1 to n+1-i do:
      M[i+j-1,j] := c:
    od:
  od:
  for j from 2 to n do:
    c := ra():
    for i from 1 to n+1-j do:
      M[i,j+i-1] := c:
    od:
  od:
  [seq([seq(M[i,j],i=1..n)],j=1..n)]:
end:


#RToeplitzMs(n,K,m): A list of m s random n by n Toeplitz matrices with entries in [1,K]
RToeplitzMs:=proc(n,K,m) local i:
  [seq(RToeplitzM(n,K), i=1..m)]:
end:


#ALnmToeplitz(X,n,m,K): Tries to find a relation between m n by n Toeplitz matrices
ALnmToeplitz:=proc(X,n,m,K,k:=0,S:={},pow:=[]):
  ALnm(X,n,m,K,RToeplitzMs,k,S,pow):
end:


#RHankelM(n,K): a random Hankel n by n matrix with entries in [1,K]
RHankelM:=proc(n,K) local ra,i,j,M,c:
  ra := rand(1..K):
  for i from 1 to n do:
    c := ra():
    for j from 1 to i do:
      M[i+1-j,j] := c:
    od:
  od:
  for j from 2 to n do:
    c := ra():
    for i from j to n do:
      M[i, j+n-i] := c:
    od:
  od:
  [seq([seq(M[i,j],i=1..n)],j=1..n)]:
end:


#RHankelMs(n,K,m): A list of m s random n by n Hankel matrices with entries in [1,K]
RHankelMs:=proc(n,K,m) local i:
  [seq(RHankelM(n,K), i=1..m)]:
end:


#ALnmHankel(X,n,m,K): Tries to find a relation between m n by n Hankel matrices
ALnmHankel:=proc(X,n,m,K,k:=0,S:={},pow:=[]):
  ALnm(X,n,m,K,RHankelMs,k,S,pow):
end:


#RCM(n,K): a random circulant n by n matrix with entries in [1,K]
RCM:=proc(n,K) local ra,i,j,M,L:
  ra := rand(1..K):
  for i from 0 to n-1 do:
    L[i] := ra():
  od:
  for i from 0 to n-1 do:
    for j from 0 to n-1 do:
      M[i,j] := L[(j-i) mod n]:
    od:
  od:
  [seq([seq(M[i,j],j=0..n-1)],i=0..n-1)]:
end:


#RCMs(n,K,m): A list of m s random n by n circulant matrices with entries in [1,K]
RCMs:=proc(n,K,m) local i:
  [seq(RCM(n,K), i=1..m)]:
end:


#ALnmC(X,n,m,K): Tries to find a relation between m n by n circulant matrices
ALnmC:=proc(X,n,m,K,k:=0,S:={},pow:=[]):
  ALnm(X,n,m,K,RCMs,k,S,pow):
end:


#RACM(n,K): a random anticirculant n by n matrix with entries in [1,K]
RACM:=proc(n,K) local ra,i,j,M,L:
  ra := rand(1..K):
  for i from 0 to n-1 do:
    L[i] := ra():
  od:
  for i from 0 to n-1 do:
    for j from 0 to n-1 do:
      M[i,j] := L[(i+j) mod n]:
    od:
  od:
  [seq([seq(M[i,j],j=0..n-1)],i=0..n-1)]:
end:


#RACMs(n,K,m): A list of m s random n by n anticirculant matrices with entries in [1,K]
RACMs:=proc(n,K,m) local i:
  [seq(RACM(n,K), i=1..m)]:
end:


#ALnmAC(X,n,m,K): Tries to find a relation between m n by n anticirculant matrices
ALnmAC:=proc(X,n,m,K,k:=0,S:={},pow:=[]):
  ALnm(X,n,m,K,RACMs,k,S,pow):
end:


#TestAC(n): tests if 3 n by n anticirculant matrices satisfy the identity suggested by ALnmAC(X,4,3,10,30)
TestAC := proc(n) local pi, T, i, j, k, row, M, Ms:
  T[[1,2,3]] := 1:
  T[[1,3,2]] := 1:
  T[[2,1,3]] := 1:
  T[[2,3,1]] := -1:
  T[[3,1,2]] := -1:
  T[[3,2,1]] := -1:
  

  Ms := []:
  for i from 1 to 3 do:
    for j from 0 to n-1 do:
      row[j] := x[i][j]:
    od:

    for j from 0 to n-1 do:
      for k from 0 to n-1 do:
        M[j,k] := row[(j+k) mod n]:
      od:
    od:
    Ms := [op(Ms), [seq([seq(M[i,j],j=0..n-1)],i=0..n-1)]]:
  od:

  evalb(expand(add(T[pi]*Mul(Mul(Ms[pi[1]],Ms[pi[2]]), Ms[pi[3]]), pi in permute(3))) = [[0$n]$n]):
end:


#RCGM(n,K,s): a random circulant n by n matrix with shift s with entries in [1,K]
RCGM:=proc(n,K,s) local ra,i,j,M,L:
  ra := rand(1..K):
  for i from 0 to n-1 do:
    L[i] := ra():
  od:
  for i from 0 to n-1 do:
    for j from 0 to n-1 do:
      M[i,j] := L[(j - i*s) mod n]:
    od:
  od:
  [seq([seq(M[i,j],j=0..n-1)],i=0..n-1)]:
end:


#RCGMs(n,K,m,s): A list of m s random n by n circulant matrices with shift s with entries in [0,K]
RCGMs:=proc(n,K,m,s) local i:
  [seq(RCGM(n,K,s), i=1..m)]:
end:


#ALnmCG(X,n,m,K): Tries to find a relation between m n by n circulant with shift s matrices
ALnmCG:=proc(X,n,m,K,s,k:=0,S:={},pow:=[]):
  ALnm(X,n,m,K,(n1,K1,m1) -> RCGMs(n1, K1, m1, s),k,S,pow):
end:


#TestCG0(n): tests if 3 n by n circulant matrices with shift 0 satisfy the identity suggested by ALnmCG(X,4,3,10,0,30)
TestCG0 := proc(n) local pi, T, i, j, k, row, M, Ms:
  T[[1,2,3]] := -1:
  T[[1,3,2]] := -1:
  T[[2,1,3]] := 1:
  T[[2,3,1]] := -1:
  T[[3,1,2]] := 1:
  T[[3,2,1]] := 1:
  

  Ms := []:
  for i from 1 to 3 do:
    for j from 0 to n-1 do:
      row[j] := x[i][j]:
    od:

    for j from 0 to n-1 do:
      for k from 0 to n-1 do:
        M[j,k] := row[k]:
      od:
    od:
    Ms := [op(Ms), [seq([seq(M[i,j],j=0..n-1)],i=0..n-1)]]:
  od:
  evalb(expand(add(T[pi]*Mul(Mul(Ms[pi[1]],Ms[pi[2]]), Ms[pi[3]]), pi in permute(3))) = [[0$n]$n]):
end:


#RTDToeplitzM(n,K): a random tridiagonal Toeplitz n by n matrix with entries in [0,K]
RTDToeplitzM:=proc(n,K) local ra,i,j,M,c:
  ra := rand(1..K):
  for i from 1 to n do:
    c := ra():
    for j from 1 to n+1-i do:
      if i <> 1 and i <> 2 then:
        M[i+j-1, j] := 0:
      else:
        M[i+j-1,j] := c:
      fi:
    od:
  od:
  for j from 2 to n do:
    c := ra():
    for i from 1 to n+1-j do:
      if j <> 2 then:
        M[i, j+i-1] := 0:
      else:
        M[i,j+i-1] := c:
      fi:
    od:
  od:
  [seq([seq(M[i,j],i=1..n)],j=1..n)]:
end:


#RTDToeplitzMs(n,K,m): A list of m s random n by n tridiagonal Toeplitz matrices with entries in [0,K]
RTDToeplitzMs:=proc(n,K,m) local i:
  [seq(RTDToeplitzM(n,K), i=1..m)]:
end:


#ALnmTDToeplitz(X,n,m,K): Tries to find a relation between m n by n Toeplitz matrices
ALnmTDToeplitz:=proc(X,n,m,K,k:=0,S:={},pow:=[]):
  ALnm(X,n,m,K,RTDToeplitzMs,k,S,pow):
end:


#RTDHankelM(n,K): a random tridiagonal Hankel n by n matrix with entries in [1,K]
RTDHankelM:=proc(n,K) local ra,i,j,M,c:
  ra := rand(1..K):
  for i from 1 to n do:
    c := ra():
    if i <> n-1 and i <> n then:
      c := 0:
    fi:
    for j from 1 to i do:
      M[i+1-j,j] := c:
    od:
  od:
  for j from 2 to n do:
    c := ra():
    if j <> 2 then:
      c := 0:
    fi:
    for i from j to n do:
      M[i, j+n-i] := c:
    od:
  od:
  [seq([seq(M[i,j],i=1..n)],j=1..n)]:
end:


#RTDHankelMs(n,K,m): A list of m s random n by n tridiagonal Hankel matrices with entries in [1,K]
RTDHankelMs:=proc(n,K,m) local i:
  [seq(RTDHankelM(n,K), i=1..m)]:
end:


#ALnmTDHankel(X,n,m,K): Tries to find a relation between m n by n tri(skew)diagonal Hankel matrices
ALnmTDHankel:=proc(X,n,m,K,k:=0,S:={},pow:=[]):
  ALnm(X,n,m,K,RTDHankelMs,k,S,pow):
end:


Vk := proc(k, K:=3) local S, s, i:
  if k = 1 then:
    return {seq([i], i=1..K)}:
  fi:

  S := Vk(k-1, K):
  {seq(seq([op(s), i], i=1..K), s in S)}:
end:


#FindRelations(N): an example of how to search for identities with matrices of size n with 2 <= n <= N.
#Do not run this without modifying (it will take too much time as written).
FindRelations := proc(N) local n, i, j, k, U, pow, S, P, f, m, result, found:
  for n from 2 to N do:
    for f in [ALnmG, ALnmS, ALnmAS, ALnmTD, ALnmToeplitz, ALnmHankel, ALnmC, ALnmAC, ALnmTDToeplitz] do:
      U := {seq(seq([i,j], i=1..n), j=1..n)}:
      for S in powerset(U) do:
        found := false:
        for m from 1 to 2*n-1 while not found do:
          P := Vk(m):
          for pow in P do:
            result := f(X, n, m, 10, 30, S, pow):
            if result <> 0 then:
              found := true:
              fprintf(`out.txt`, `f=%s, n=%d, m=%d, S=%a, pow=%a, relation=%a\n`,  
                       f, n, m, S, pow, result):
              fflush(`out.txt`):
            fi:
          od:
        od:
      od:
    od:
  od:
end: