######################################################################
## 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  : April 2, 2025  `):
print():
print(`Written by the 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(` `):

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, Mul, RM, RMs  `):


elif nargs=1 and args[1]=ALnm then
print(`ALnm(X,n,m,K): Tries to find a relation between ANY m  n by n matrices by picking random n by n matrices with entries in [1,K]. Try:`):
print(`ALnm(X,1,2,10); Alnm(X,2,4,10);`):
print(`This will take a long time:  ALnm(X,3,6,5);`):

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]=RM then
print(`RM(n,K): a random n by n matrix with entries in [1,K]. Try:`):
print(`RM(5,100);`):

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

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). Try:
#Mul([[a1,a2]],[[b1],[b2]]);
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:



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

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

#ALnm(X,n,m,K): Tries to find a relation between m  n by n matrices. Try:
#ALnm(X,1,2,5); Alnm(X,2,4); Alnm(X,3,6);
ALnm:=proc(X,n,m,K) local A, pi, eqs, M, var, i, var1, v, JK,M1,c,i1:
  var := {seq(c[pi], pi in permute(m))}:
JK:=add(c[pi]*X[pi],pi in permute(m)):

  eqs := {}:

  for i from 1 to trunc(m!/n^2)+5 do
    A := RMs(n,K,m):
   M:=[[0$n]$n]:
    for pi in permute(m) do
        M1:=A[pi[1]]:
    for i1 from 2 to m do
     M1:=Mul(M1,A[pi[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 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(var1[1]=1,subs(var,JK)):
end: