# Ok to post homework
# Lucy Martinez, 04-18-2025, Assignment 23


######################
#Problem 1: A rooted ternary tree is either a single vertex,
# or an ordered tree where every vertex is either a leaf (has no children)
# or has exactly three children. 
# Write ternary analogs for BT(n), T(N), Images(T), UBT(n), and UT(N),
# calling them: TT(n), T3(N), Images3(T), UBT3(n), UT3(N), respectively.

# What is UT3(40)? Is it in the OEIS? What is its A number.
# [Hint, the functional equation, according to Polya is
# r(x)=x+1/6*(r(x)^3+ 3*r(x)*r(x^2)+ 2*r(x^3))
# Note that the number of leaves is always odd so only extract
# the odd coefficients]

# Answer: The OEIS entry is A000598


#TT(n): the set of ternary trees with n leaves 
TT:=proc(n) local T1,T2,T3,k,j,i,T,t1,t2,t3:
if n=1 then
 return({[]}):
fi:


T:={}:
for k from 1 to n-2 do
 for j from 1 to n-k-1 do
  #i:=n-k-j:
  #if i<1 then
  # next:
  #fi:
  T1:=TT(k):
  T2:=TT(j):
  T3:=TT(n-k-j):
  T:=T union {seq(seq(seq([t1,t2,t3],t1 in T1),t2 in T2),t3 in T3)}
 od:
od:
T:
end:


# OEIS entry A001764
#T3(N): the first N terms of the sequence enumerating the number
# of ternary trees with n leaves
T3:=proc(N) local x,f,i:
f:=x:
for i from 1 to N do
 f:=expand(x+f^3):
 f:=add(coeff(f,x,i)*x^i,i=1..N):
od:
[seq(coeff(f,x,i),i=1..N,2)]:
end:


#Images3(T):=[Images(T1),Images(T2),Images(T3)]
Images3:=proc(T) local T1,T2,T3,t1,t2,t3,S:
if T=[] then
 return({[]}):
fi:
T1:=Images3(T[1]):
T2:=Images3(T[2]):
T3:=Images3(T[3]):
S:={}:
for t1 in T1 do
 for t2 in T2 do
  for t3 in T3 do
   S:=S union {[t1,t2,t3],[t1,t3,t2],[t2,t1,t3],
        [t2,t3,t1],[t3,t1,t2],[t3,t2,t1]}:
  od:
 od:
od:
S:
end:


#UBT3(n): the set of all unlabeled ternary trees
# (one representative of each)
UBT3:=proc(n) local S1,S2,S3:
S1:=TT(n):
S2:={}:
while S1<>{} do
 S3:=Images3(S1[1]):
 S2:=S2 union {S1[1]}:
 S1:=S1 minus S3:
od:
S2:
end:


# OEIS entry A000598
#UT(n): the first N terms of the set of all 
# unlabeled rooted ternary trees
UT3:=proc(N) local x,r,i:
r:=x:
for i from 1 to N do
 #r(x)=x+1/6*(r(x)^3+ 3*r(x)*r(x^2)+ 2*r(x^3))
 r:=expand(x+1/6*(r^3+ 3*r*subs(x=x^2,r)+ 2*subs(x=x^3,r))):
 r:=add(coeff(r,x,i)*x^i,i=1..N):
od:
[seq(coeff(r,x,i),i=1..N,2)]:
end:




#Problem 2: If you type ALnmG(X,3,4,10,{[1,3],[2,1],[2,3],[3,1],[3,2]});
# In the Maple package AL.txt you would get six identities
# that any four 3x3 matrices A1,A2,A3,A4 where
# the only diagonal entries as well as the [1,2] entry are 0,
# i.e. matrices of the form
# [[a11,a12,0],[0,a22,0],[0,0,a33]]

# One of them is:
# X[[1, 2, 3, 4]]-X[[1, 2, 4, 3]]-X[[2, 1, 3, 4]]+X[[2, 1, 4, 3]]

# Convince yourself that this is stating that for any four such matrices
# (A1A2-A2A1)(A3A4-A4A3)= 0 matrix
# Or in our language
# Mul(Mul(A1,A2)-Mul(A2,A1), Mul(A3,A4)-Mul(A4,A3))=
# [[0,0,0],[0,0,0],[0,0,0]] .

# Prove it rigorously.




####################################from previous classes:
#C23.txt: April 17, 2025
Help23:=proc():
print(`BT(n), T(N), Images(T), UBT(n), UT(N)`):
end:

#Def. A full binary tree is either a single vertex called the root,
# or it has a root that has exactly two children
# T=*  or
#       *
#    T1   T2
#T=[] or T=[T1,T2]
#BT(1)={[]}, BT(2)={[[],[]]}, BT(3)={[[],[[],[]]],...}

#BT(n): the SET of full binary trees with n leaves
BT:=proc(n) local T1,T2,k,T,t1,t2:
if n=1 then
 return({[]}):
fi:

#k is the number of leaves in the left children
T:={}:
for k from 1 to n-1 do
 T1:=BT(k):
 T2:=BT(n-k):
 T:=T union {seq(seq([t1,t2], t1 in T1),t2 in T2)}:
od:
T:
end:

#Wt(T):=x^NumberOfLeaves
#WT([T1,T2])=WT(T1)*WT(T2)
#T=[] or T=[T1,T2]
#f(x)=Sum(weight(T), all trees)
#f(x)=x+f(x)*f(x)=x+f(x)^2
# t_n:=coeff of x^n in -(1-4*x)^(1/2)/2
# -binomial(1/2,n)*4^n/2

#T(N): the first N terms of the sequence enumerating the number
# of full binary trees with n leaves
T:=proc(N) local x,f,i:
f:=x:
for i from 1 to N do
 f:=expand(x+f^2):
 f:=add(coeff(f,x,i)*x^i,i=1..N):
od:
[seq(coeff(f,x,i),i=1..N)]:
end:



#Images(T):=[ Images(T1), Images(T2) ]
Images:=proc(T) local T1,T2,t1,t2,S:
if T=[] then
 return({[]}):
fi:
T1:=Images(T[1]):
T2:=Images(T[2]):
S:={}:
for t1 in T1 do
 for t2 in T2 do
  S:=S union {[t1,t2],[t2,t1]}:
 od:
od:
S:
end:


#A001190: Wedderburn-Etherington numbers
#UBT(n): the set of all unlabeled full binary trees
# (one representative of each)
UBT:=proc(n) local S1,S2,S3:
S1:=BT(n):
S2:={}:
while S1<>{} do
 S3:=Images(S1[1]):
 S2:=S2 union {S1[1]}:
 S1:=S1 minus S3:
od:
S2:
end:


#UT(N): the first N terms of the set of all unlabeled full binary trees
UT:=proc(N) local x,f,i:
f:=x:
for i from 1 to N do
 f:=expand(x+(f^2+subs(x=x^2,f))/2):
 f:=add(coeff(f,x,i)*x^i,i=1..N):
od:
[seq(coeff(f,x,i),i=1..N)]:
end:



######################################################################
## 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(` RMs, RMgs`):

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


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]=ALnmG then
print(`ALnmG(X,n,m,K,S): Tries to find a relation between m  n by n matrices with entries in S that are zero. Try:`):
print(`ALnmG(X,2,3,10,{[1,2]}); `):

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]=RMg then
print(`RMg(n,K,S): inputs a positive integer n and a pos. integer K, and a subset of [1,n],[1,n], outputs random n by n matrix with entries in [1,K] that is 0 at S.`):
print(`For example to get a random lower-triangular 2 by 2 matrix Try:`):
print(`RMg(2,10,{[1,2]});`):


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);`):

elif nargs=1 and args[1]=RMgs then
print(`RMgs(n,K,S,m): A list of m random nxn matrices with entries from [1,K] that vanish at S. Try:`):
print(`RMgs(2,10,{[1,2]},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 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(var1[1]=1,subs(var,JK)):
end:






#RMg(n,K,S): inputs a positive integer n and a pos. integer K, and a subset of [1,n],[1,n], outputs random n by n matrix with entries in [1,K] that is 0 at S.
#For example to get a random lower-triangular 2 by 2 matrix Try:
#RMg(2,10,{[1,2]});
RMg:=proc(n,K,S) local ra,i,j,T:
if S minus {seq(seq([i,j],j=1..n),i=1..n)}<>{} then
 RETURN(FAIL):
fi:

 ra := rand(1..K):

for i from 1 to n do
 for j from 1 to n do
  if member([i,j],S) then
     T[i,j]:=0:
  else
   T[i,j]:=ra():
  fi:
 od:
od: 

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


#RMgs(n,K,S,m): A list of m random nxn matrices with entries from [1,K] that vanish at S. Try:
#RMgs(2,10,{[1,2]},3);
RMgs:=proc(n,K,S,m) local i:
  [seq(RMg(n,K,S), i=1..m)]:
end:


#ALnmG(X,n,m,K,S): Tries to find a relation between m  n by n matrices with entries in S that are zero. Try:
#ALnmG(X,2,3,10,{[1,2]}); 
ALnmG:=proc(X,n,m,K,S) local A, pi, eqs, M, var, i, var1, v, JK,M1,c,i1,v1:
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)+30 do
    A:=RMgs(n,K,S,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 var=NULL then
  RETURN(FAIL):
fi:

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

JK:=subs(var,JK):

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


{seq(coeff(JK,v1), v1 in var1)}:

end: