#C20.txt,April 8, 2019
Help:=proc(): print(`BT(n) ,NuLeaves(T) , RandMat(m1,m2,K) `): 
print(`RandMatList(P,K) , EvalBT(L,T), PriceM(P,T) `):
end:

with(linalg):

#Lemma: The cost (in multiplication operations) of multiplying an mxn matrix by an nxp matrix
#is mnp
#Cor. The cost of A=[m1xm2], B=[m2xm3] C=[m3xm4] 
#the cost of (AB)C=m1m2m3+ m1m3m4
#The cost of A(BC)=m2m3m4+ m1m2m4
#A1, A2, ..., An, [m1,m2,m3, ..., m[n+1]] 

#input a positive integer n
#returns the set of all binary trees with exactly n leaves 
BT:=proc(n) local S,a,b,i:

	option remember:

	if n<=0 then
		return {}:
	fi:

	if n=1 then 
		return {[]}:
	fi:
	
	S:={}:
	
	#to get a tree on n leaves, you can join any two trees whose numbers of leaves sum to n 
	for i from 1 to n-1 do 
	
		for a in BT(i) do
		
			for b in BT(n-i) do
			
				S:=S union {[a,b]}:
			
			od:
		
		od:
	
		
	
	od:
	
	return S:

end:

#NuLeaves(T): the number of leaves in the binary tree T
NuLeaves:=proc(T)
option remember:
if T=[] then
 RETURN(1):
else
 NuLeaves(T[1])+NuLeaves(T[2]):
fi:

end:



#RandMat(m1,m2,K): a random m1 by m2 matrix with entries from -K to K
#given as a list of lists 
RandMat:=proc(m1,m2,K) local ra,i,j:
ra:=rand(-K..K):
[seq( [seq( ra(), j=1..m2)],i=1..m1)]:
end:

#RandMatList(P,K): given a list of positive integers P of length n+1
#and a pos. integer K outputs a list of n compatible matrices ready to
#be multiplied
RandMatList:=proc(P,K) local n,i:
n:=nops(P)-1:
[seq(RandMat(P[i],P[i+1],K),i=1..n)]:
end:


#EvalBT(inputs a list L of matrices and a binary tree with the number of leaves equal to
#the length of L, outputs the matrix product L[1]...L[n] suggested by the tree
EvalBT:=proc(L,T) local n,P,i,P1,k,A,B:
n:=nops(L):
if NuLeaves(T)<>n then
 RETURN(FAIL):
fi:

P:=[seq([nops(L[i]),nops(L[i][1])],i=1..n)]:

for i from 1 to n-1 do
 if P[i][2]<>P[i+1][1] then
   ERROR(L[i], `is not compatible with`, L[i+1] ):
 fi:
od:

#P1:=[seq(P[i][1],i=1..n-1),P[n][1],P[n][2]]:

if n=1 then
  RETURN(L[1]):
fi:

#k:=number of leaves of the left tree of the root
k:=NuLeaves(T[1]):
A:=EvalBT( [op(1..k,L)] ,  T[1]):
B:=EvalBT([op(k+1...n ,L)] ,T[2]):
convert(multiply(A,B),list):
end:

#PriceM(P,T): inputs a profile of lenght n+1 and a binary tree with n leaves
#outputs the number of multiplications
PriceM:=proc(P,T) local n,k:
n:=NuLeaves(T):
if n=1 then
 RETURN(0):
fi:

if n=2 then
 RETURN(P[1]*P[2]*P[3]):
fi:

k:=NuLeaves(T[1]):


PriceM([op(1..k+1,P)],T[1])+PriceM([op(k+1..n+1,P)],T[2])+ P[1]*P[k+1]*P[n+1]:
   
end: