# Ok to post homework
# Lucy Martinez, 03-07-2025, Assignment 13

with(combinat):

######################
#Problem 1: On p. 130 of Thomas G. Wong's wonderful book
# Introduction to Classical and Quantum Computing it is proved that HXH=Z
# This can be rephrased, in our program as
# Mul(Mul(C1QG()[7],C1QG()[2]),C1QG()[7])=C1QG()[4]

# Find all triples 2 ≤ i,j,k ≤ 7 (there are 6^3=216 such triples) such that
# Mul(Mul(C1QG()[i],C1QG()[j]),C1QG()[k])
# is equal to a member of C1QG(), say C1QG()[r].
# Output the set of all such quadruples [i,j,k,r]
# (in particular [7,2,7,4] should be a member of that set).
# In other words find all such identities where a product of three
# (not necessarily distinct) mmembers of C1QG() equals to one of them.
# How big is that set? for some 1 ≤ r ≤ 7

#The following is the output with 54 members.
# Observe that [7, 2, 7, 4] is a member.
Ta:=[[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]]:

FindTriples:=proc() local i,j,k,r,S,M:
S:=[]:
for i from 2 to 7 do
 for j from 2 to 7 do
  for k from 2 to 7 do
   M:=C1QG():
   if member(Mul(Mul(C1QG()[i],C1QG()[j]),C1QG()[k]),M,'r') then
    S:=[op(S),[i,j,k,r]]:
   fi:
  od:
 od:
od:
S: 
end:




#Problem 2: Adapt the above problem to find, more generally
# all quadruples 2 ≤ i,j,k ≤ 7 and c constant (obviously of absolute value 1)
# such that Mul(Mul(C1QG()[i],C1QG()[j]),C1QG()[k])=c*C1QG()[r]
# Output the set of all such [i,j,k,r,c]

#Answer: The following is the output with 80 members.
# Observe that [7, 2, 7, 4, 1] is a member of Q
Q:= [[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, 6, 6, 4, 1], [5, 7, 7, 5, 1],
	 [6, 2, 2, 6, 1], [6, 2, 6, 2, (1/2 + I/2)*sqrt(2)],
	 [6, 3, 3, 6, 1], [6, 3, 6, 3, (1/2 + I/2)*sqrt(2)],
	 [6, 4, 4, 6, 1], [6, 5, 6, 4, 1], [6, 6, 5, 4, 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]]

FindQuad:=proc() local S,i,j,k,r,c,M1,M2,i1,i2:
S:=[]:
for i from 2 to 7 do
 for j from 2 to 7 do
  for k from 2 to 7 do
  M1:=Mul(Mul(C1QG()[i],C1QG()[j]),C1QG()[k]):
   for r from 1 to 7 do
    M2:=C1QG()[r]:
	c:={}:
	for i1 from 1 to 2 do
	 for i2 from 1 to 2 do
	  if M2[i1][i2]<>0 then
	   c:=c union {simplify(M1[i1][i2]/M2[i1][i2])}:
	  fi:
	 od:
	od:
	if nops(c)=1 and abs(c[1])=1 then
	 S:=[op(S),[i,j,k,r,c[1]]]:
	fi:
   od:
  od:
 od:
od:
S:
end:



#Problem 3: Write a procedure ApplyT(T,X,n,u)
# where T=[T1,V], T1 is a table and V is nBasis(n,X),
# and u is an arbitrary linear combination of members of nBasis(n,X),
# that uses the linearity of the table, and outputs
# the result of applying the linear transformation T
# (given as a table defined on the canonical basis),
# as a linear combination of the basis elements.

# If T=UMtoT(TP(H,Z),2,X), what is 
# ApplyT(T,X,2,X[[0,0]]+I*X[[0,1]]+2*I*X[[1,0]]+5*X[[1,1]]); ?

# Answer:
# (((5 - I)*X[[1, 1]] - (5 + I)*X[[0, 1]] + (1 - 2*I)*X[[1, 0]] + (1 + 2*I)*X[[0, 0]])*sqrt(2))/2

ApplyT:=proc(T,X,n,u) local T1,V,A,B,i:
T1:=T[1]:
V:=T[2]:
A:=[coeffs(u,{op(V)},'terms')]:
B:=[terms]:
if nops(A)<>nops(B) then
 return(FAIL):
fi:
simplify(add(A[i]*T1[B[i]] , i=1..nops(A))):
end:



####################################previous classes:
#C13.txt, March 06, 2025
#Computing Quantum gates
with(combinat):
with(linalg):

Help13:=proc(): 
print(`C1QG(), C1QGnames(),BV(n),nBasis(n,X),UMtoT(M,n,X)`):
end:


C1QGnames:=proc():
[` Identity `, ` PauliX `, `PauliY `, `PauliZ `, ` PhaseS `, `TPi/8` , ` Hadamard `]:
end:


C1QG:=proc():
[
[[1,0],[0,1]],
[[0,1],[1,0]],
[[0,-I],[I,0]],
[[1,0],[0,-1]],
[[1,0],[0,I]],
[[1,0],[0,sqrt(2)/2+sqrt(2)/2*I]],
expand([[1,1],[1,-1]]/sqrt(2))
]:
end:


#BV(n):inputs a non-neg. integer and outputs the list of length 2^n-i
# of all the 0-1 vectors in lex. order
BV:=proc(n) local V,i:
option remember:
if n=0 then
 return([[]]):
fi:
V:=BV(n-1):
[seq([0, op(V[i])], i=1..nops(V)),seq([1, op(V[i])], i=1..nops(V))]:
end:


#nBasis(n,X): the basis of n-quantum gates circuits X[e1,e2,...,en]
# corresponding to |e1>|e2>...|en>
nBasis:=proc(n,X) local V,i:
V:=BV(n):
[seq(X[V[i]],i=1..nops(V))]:
end:


#UMtoT(M,n,X): inputs a unitary matrix M of size 2^n by 2^n 
# (corresponding to a quantum circuit with n qbits)
# and outputs the corresponding linear transformation as it acts
# on the canonical basis nBasis(n,X)
UMtoT:=proc(M,n,X) local T,V,i,j:
if not IsUM(M) then
 print(M, `is not unitary`):
 return(FAIL):
fi:

if nops(M)<>2^n then
 return(FAIL):
fi:

V:=nBasis(n,X):
for j from 1 to nops(V) do
 T[V[j]]:=add(V[i]*M[i][j],i=1..nops(M)):
od:
[op(T),V]:
end:


#TtoUM(T,n): goes from the table to unitary matrix 
TtoUM:=proc(T,n) local T1,B,i,j:
B:=T[2]:
T1:=T[1]:
if nops(B)<>2^n then
 return(FAIL):
fi:
[seq([seq( coeff(T1[B[j]],B[i],1) ,j=1..2^n)],i=1..2^n)]:
end:


####from previous classes
HelpOld:=proc(): print(`CT(A), IsUM(A), Mul(A,B), ES(), TP(A,B)`):end:

#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:

#CT(A): the conjugate transpose of A
CT:=proc(A): local i,j,n:
n:=nops(A):
[seq([seq(conjugate(A[j][i]),j=1..n)],i=1..n)]:
end:


Con:=proc(A): local i,j,n:
n:=nops(A):
[seq([seq(A[j][i],j=1..n)],i=1..n)]:
end:


#IsUM(A): Is the matrix A unitary?
IsUM:=proc(A) local n,i,j, P:
n:=nops(A):
P:=simplify(Mul(A,CT(A))):
evalb(P=[seq([0$(i-1),1,0$(n-i)],i=1..n)]):
end:


#ES(): the four fullly entangled states in the Alice-Bob combined system
#based on p. 166of Susskind-Friedman's book
#[singlet, T1,T2,T3]
#[uu.ud,du,dd]
ES:=proc()
[
[0,1/sqrt(2),-1/sqrt(2),0],
[0,1/sqrt(2),1/sqrt(2),0],
[1/sqrt(2),0,0, 1/sqrt(2)],
[1/sqrt(2),0,0,-1/sqrt(2)]
]
end:

#TP(A,B):the tensor product of matrix A and matrix B (using our data structure
#of lists-of-lists
#Let A be an a1 by a2 matrix and 
#Let B be a b1 by b2 matrix
#TP(A,B): is a1*b1 by a2*b2 matrix
TP:=proc(A,B) local a1,a2,b1,b2,AB,i,i1,j,j1,AB1:
a1:=nops(A): a2:=nops(A[1]):
b1:=nops(B): b2:=nops(B[1]):
#The rows of TP(A,B) are called [i,i1] where 1<=i<=a1 and 1<=i1<=b1
#The columns of TP(A,B) are called [j,j1] where 1<=j<=a2 and 1<=j1<=b2
AB:=[]:
for i from 1 to a1 do
 for i1 from 1 to b1 do
 #AB1 is the [i,i1]-row of TP(A,B)
  AB1:=[]: 
  for j from 1 to a2 do
   for j1 from 1 to b2 do
    AB1:=[op(AB1),A[i][j]*B[i1][j1]]:
   od:
  od:
  AB:=[op(AB),AB1]:
 od:
od:
AB:
end:

##End Old stuff