Help15:=proc(): print(` UM(alpha,beta,gamma,delta), ApplyT(T,u), RM(a,b,K), RQC(n,r,k),MergeV(v1,v2,S1,S2), Embed1(n,S,M) `):end:

with(combinat):

#MergeV(v1,v2,S1,S2): Given vector v1 supported on a list S1, and a vector v2 supported on a list S2, let n:=nops(v1)+nops(v2)  such that {op(S1) union op(S2)}={1,.,n}
#outputs the vector v such that v[]=v1[i] and v[S2[i]]:=v2[i]. Try:
#MergeV([0,1],[1,0],[1,3],[2,4]);
MergeV:=proc(v1,v2,S1,S2) local n,v,i:
n:=nops(v1)+nops(v2):
if {op(S1),op(S2)}<>{seq(i,i=1..n)} then
 RETURN(FAIL):
fi:

if {op(S1)} intersect {op(S2)}<>{} then
 RETURN(FAIL):
fi:

for i from 1 to nops(S1) do
v[S1[i]]:=v1[i]:
od:

for i from 1 to nops(S2) do
v[S2[i]]:=v2[i]:
od:

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

#Embed1(n,S,M): inputs a pos. integer n and a subset S given as list of length k, and 2 2^k by 2^k matrix acting on the qubits in S\
#outputs the  correponding 2^n by 2^n matrices where the other qubits are the same. Try:
#Embed1(4,{1},[[a11,a12],[a21,a22]]);
Embed1:=proc(n,S,M) local k,B1,T,v1,S1:
k:=nops(S):
S1:={seq(i,i=1..n)} minus S:
if not (type(M,list) and nops(M)=2^k and nops(M[1])=2^k) then
 RETURN(FAIL):
fi:

Bk:=BV(k):

Bk1:=BV(n-k):

for i from 1 to nops(Bk) do
for j from 1 to nops(Bk) do
 T1[Bk[i],Bk[j]]:=M[i][j]:
od:
od:

Bn:=BV(n):

for i from 1 to nops(Bn) do
for j from 1 to nops(Bn) do
 T[Bn[i],Bn[j]]:=0:
od:
od:

for i from 1 to nops(BK1) do
for i1 from 1 to BK do
for j1 from 1 to BK do
T[MergeV(BK[i1],BK1[i],S,Sk),MergeV(BK[j1],BK1[i],S,Sk)]:=T1[BK[i1],BK[j1]]:
od:
od:
od:
op(T):
end:


#From hw4, Adapted from Auorora Hively's hw4  (See https://sites.math.rutgers.edu/~zeilberg/EM25/hw4posted/hw4AuroraHiveley.txt )
# UM(alpha,beta,gamma,delta):  implements box 1.1. on p. 20 of Nielsen-Chuang 
# note: alpha, beta, delta, gamma are all real valued

UM := proc(alpha, beta, gamma, delta) local A,X,Y,Z:
A := [[exp(I*alpha),0],[0,exp(I*alpha)]]:
X := [[exp(-I*beta/2),0],[0,exp(I*beta/2)]]:
Y := [[cos(gamma/2), -sin(gamma/2)], [sin(gamma/2), cos(gamma/2)]]:
Z :=[[exp(-I*delta/2),0],[0,exp(I*delta/2)]]:

simplify(evalc(Mul(Mul(A,X),Mul(Y,Z)))):

end:

#Adapted From Joseph Koutsoutis's hw13 (see https://sites.math.rutgers.edu/~zeilberg/EM25/hw13posted/hw13JosephKoutsoutis.txt)
ApplyT := proc(T,u) local T1, V, v, val:
  T1, V := op(T):
  val := 0:
  for v in V do:
    val += coeff(u, v, 1) * T1[v]:
  od:
  expand(simplify(val)):
end:



#RQC(n,r,k): generates a random quantum circuit with n qubits of length r, with gates involving at most k qubits at a time
#obtained by randomly tensor-producting C1QG()[i]. The outputs is a list of length r where S_i is the support of M_i
#[n,[seq([S_i,M_i],i=1..r)]]:
RQC:=proc(n,r,k) local ra,M,k1,S1,M1,i1,i:
ra:=rand(2..7):
M:=[]:
for i from 1 to r do
 k1:=rand(1..k)():
 S1:=randcomb(n,k1):
 M1:=C1QG()[ra()]:
 for i1 from 2 to k1 do
  M1:=TP(M1,C1QG()[ra()]):
 od:
M:=[op(M),[S1,M1]]:
od:
[n,M]:
end:

#Old stuff
#March 6, 2025 Computing Quantum gates
Help13:=proc(): print(`C1QG(), C1QGnames(), BV(n), nBasis(n,X), UMtoT(M,n,X) , TtoUM(T,n) `): 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 n and outputs the list of length 2^n of all 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] corresponds 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 (correpsonding to a quantum circuit with n
#qbits 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:=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