# hw13 - Pablo Blanco
# OK to post.

read(`C13.txt`):

####
# kMul(P): Inputs a list of k matrices (lists of lists) P=[P1,..,Pk] and multiplies them in order.
####
kMul:=proc(P) local k, M, i:
 k:=nops(P):
 if k < 2 then
  error "need at least two matrices in this list of matrices":
 fi:

 M:=P[1]:
 for i from 1 to k-1 do:
  M:=Mul(M,P[i+1]):
 od: 

 M:
end:

#######################
# As described in the homework page.
# Output: 
# {[2, 2, 2, 2], [2, 2, 3, 3], [2, 2, 4, 4], [2, 2, 5, 5], [2, 2, 6, 6], [2, 2, 7, 7], [2, 3, 3, 2], [2, 4, 4, 2], [2, 7, 4, 7], [2, 7, 7, 2], [3, 2, 2, 3], [3, 3, 2, 2], [3, 3, 3, 3], [3, 3, 4, 4], [3, 3, 5, 5], [3, 3, 6, 6], [3, 3, 7, 7], [3, 4, 4, 3], [3, 5, 2, 5], [3, 7, 7, 3], [4, 2, 2, 4], [4, 3, 3, 4], [4, 4, 2, 2], [4, 4, 3, 3], [4, 4, 4, 4], [4, 4, 5, 5], [4, 4, 6, 6], [4, 4, 7, 7], [4, 5, 4, 5], [4, 5, 5, 1], [4, 6, 4, 6], [4, 7, 2, 7], [4, 7, 7, 4], [5, 2, 2, 5], [5, 3, 3, 5], [5, 4, 4, 5], [5, 4, 5, 1], [5, 5, 4, 1], [5, 7, 7, 5], [6, 2, 2, 6], [6, 3, 3, 6], [6, 4, 4, 6], [6, 7, 7, 6], [7, 2, 2, 7], [7, 2, 7, 4], [7, 3, 3, 7], [7, 4, 4, 7], [7, 4, 7, 2], [7, 7, 2, 2], [7, 7, 3, 3], [7, 7, 4, 4], [7, 7, 5, 5], [7, 7, 6, 6], [7, 7, 7, 7]}
# Size of output: 54
#######################
FindAllQuads:=proc() local S, P,i,j,k,memIndex,Out:
 S:=C1QG():
 Out:={}:
 for i from 2 to 7 do:
  for j from 2 to 7 do:
   for k from 2 to 7 do:
     if member(kMul([S[i],S[j],S[k]]),S,'memIndex') then:
      Out:= Out union {[i,j,k,memIndex]}:
     fi:
   od:
  od:
 od:
end:


#######################
# membG(M,S,'r','c'): Given a collection (list) of matrices S and a matrix M. All matrices should be square and of the same dimension.
# Returns TRUE if there is M=cA for some member A of S and a non-zero constant c. Stores the result in the constant c and index r in the collection.
# ASSUMPTION: S should not have members that are multiples of each other.
# ASSUMPTION: All matrices in S should correspond to bijective linear maps.
# Returns FALSE otherwise.
#######################
membG:=proc(M,S,`r`,`c`) local A,n,i:
 n:=nops(M): # assume M is a list of lists.
 if n < 2 then
  error "matrix in first argument should be at least dim 2 (and square).":
 fi: 
 r:=0:

 for A in S do:
  r++:  # index of A
  i:=1:
  while A[1][i]=0 do: # no row is all zeros
   i++:
  od:
   
  if M[1][i] = 0 then
   continue:
  else: 
   c:=(M[1][i])/(A[1][i]): # this is a CANDIDATE c value

   if evalb(M = expand((M[1][i])/(A[1][i])*A)) then: 
    return(true): 
   fi:
  fi:
  
  
 od: 

 false:
end:


####################
# Output: {[2, 2, 2, 2, 1], [2, 2, 3, 3, 1], [2, 2, 4, 4, 1], [2, 2, 5, 5, 1], [2, 2, 6, 6, 1], [2, 2, 7, 7, 1], [2, 3, 2, 3, -1], [2, 3, 3, 2, 1], [2, 3, 4, 1, I], [2, 4, 2, 4, -1], [2, 4, 3, 1, -I], [2, 4, 4, 2, 1], [2, 5, 3, 5, -1], [2, 5, 5, 3, -I], [2, 7, 4, 7, 1], [2, 7, 7, 2, 1], [3, 2, 2, 3, 1], [3, 2, 3, 2, -1], [3, 2, 4, 1, -I], [3, 3, 2, 2, 1], [3, 3, 3, 3, 1], [3, 3, 4, 4, 1], [3, 3, 5, 5, 1], [3, 3, 6, 6, 1], [3, 3, 7, 7, 1], [3, 4, 2, 1, I], [3, 4, 3, 4, -1], [3, 4, 4, 3, 1], [3, 5, 2, 5, 1], [3, 5, 5, 2, I], [3, 7, 3, 7, -1], [3, 7, 7, 3, 1], [4, 2, 2, 4, 1], [4, 2, 3, 1, I], [4, 2, 4, 2, -1], [4, 3, 2, 1, -I], [4, 3, 3, 4, 1], [4, 3, 4, 3, -1], [4, 4, 2, 2, 1], [4, 4, 3, 3, 1], [4, 4, 4, 4, 1], [4, 4, 5, 5, 1], [4, 4, 6, 6, 1], [4, 4, 7, 7, 1], [4, 5, 4, 5, 1], [4, 5, 5, 1, 1], [4, 6, 4, 6, 1], [4, 7, 2, 7, 1], [4, 7, 7, 4, 1], [5, 2, 2, 5, 1], [5, 2, 5, 2, I], [5, 3, 3, 5, 1], [5, 3, 5, 3, I], [5, 4, 4, 5, 1], [5, 4, 5, 1, 1], [5, 5, 2, 3, I], [5, 5, 3, 2, -I], [5, 5, 4, 1, 1], [5, 7, 7, 5, 1], [6, 2, 2, 6, 1], [6, 2, 6, 2, sqrt(2)/2 + sqrt(2)*I/2], [6, 3, 3, 6, 1], [6, 4, 4, 6, 1], [6, 7, 7, 6, 1], [7, 2, 2, 7, 1], [7, 2, 7, 4, 1], [7, 3, 3, 7, 1], [7, 3, 7, 3, -1], [7, 4, 4, 7, 1], [7, 4, 7, 2, 1], [7, 7, 2, 2, 1], [7, 7, 3, 3, 1], [7, 7, 4, 4, 1], [7, 7, 5, 5, 1], [7, 7, 6, 6, 1], [7, 7, 7, 7, 1]}
# Output Size: 76
#
####################
FindAllQuints:=proc() local S,Out,i,j,k,r,c:
 S:=C1QG():
 Out:={}:
 for i from 2 to 7 do:
  for j from 2 to 7 do:
   for k from 2 to 7 do:
     if membG(kMul([S[i],S[j],S[k]]),S,'r','c') then:
      Out:= Out union {[i,j,k,r,c]}:
     fi:
   od:
  od:
 od:
end:


####################
# NOTE: the arguments X,n were redundant for the implementation of this procedure.
# Note that the symbol X is stored in T anyway. 
# They can be removed without any functional loss.
#####################
ApplyT:=proc(T,X,n,u) local T1,V,L,m,i:
 T1:=T[1]:
 V:=T[2]:
 L:=[op(u)]:
 m:=nops(L):
 for i from 1 to m do:
  if nops(L[i]) = 1 then: 
   L[i]:=1+L[i]:                 # a little trick. This is to handle entries like X[[0,0]]. This way op(L[i]) will always extract coefficients correctly.
  fi:
 od:
 collect(add(op(1,L[i])*T1[op(2,L[i])]  ,i=1..m),V):
end:


################
# Output test: 
# ApplyT(T,X[[0,0]]+I*X[[0,1]]+2*I*X[[1,0]]+5*X[[1,1]]):
# (sqrt(2)*I/2 + (5*sqrt(2))/2)*X[[0, 0]] + (sqrt(2)/2 + sqrt(2)*I)*X[[0, 1]] + (sqrt(2)*I/2 - (5*sqrt(2))/2)*X[[1, 0]] + (sqrt(2)/2 - sqrt(2)*I)*X[[1, 1]]
###############