# Ok to post homework
# Lucy Martinez, 03-12-2025, Assignment 14

with(combinat):

######################
#Problem 1: Write a procedure RandomBell1(K)
# (1) After defining ra:=proc(K): rand(-K..K)():end:
#     Keeping the same A0,A1,B0,B1, as in the wikipedia article
#     A0:=C1QG()[4],A1:=C1QG()[2], B0:= -(C1QG()[2]+C1QG()[4])/sqrt(2),
#     B1:=(C1QG()[2]-C1QG()[4])/sqrt(2)
# (2) Picks a random state (in the 2-qubit space generated by 00,01,10,11,
#     sets V=[ra(K)+I*ra(K),ra(K)+ra(K)*I,ra(K)+I*ra(K),ra(K)+I*ra(K)]:
# (3) NORMALIZE it! , i.e. divide V by its NORM. rename the new state V.
# Computes Bell(A0,A1,B0,B1,V)
RandomBell1:=proc(K) local A0,A1,B0,B1,V,v,i,k:
A0:=C1QG()[4]:
A1:=C1QG()[2]:
B0:=-expand((C1QG()[2]+C1QG()[4])/sqrt(2)):
B1:=expand((C1QG()[2]-C1QG()[4])/sqrt(2)):
V:=[ra(K)+I*ra(K),ra(K)+ra(K)*I,ra(K)+I*ra(K),ra(K)+I*ra(K)]:
k:=sqrt(add(abs(V[i])^2,i=1..nops(V))):
V:=[seq(V[i]/k,i=1..nops(V))]:
Bell(A0,A1,B0,B1,V):
end:

#ra(K): random number between -K and K
ra:=proc(K):
rand(-K..K)():
end:



#Problem 2: Write a procedure RandomBell(K,N)
# that runs RandomBell1(K) N times and returns
# (1) The average value over all N trials.
# (2) The fraction of times it exceeds 2,
#     thereby refuting Albert (and Boris and Nathan)

# What did you get for RandomBell(10,1000)?
# Answer: I ran it a few times, here are the results:
# [0.03855372122, 0.005000000000]
# [0.006640886106, 0.008000000000]
# [0.01375946829, 0.007000000000]

RandomBell:=proc(K,N) local R,i,co,co1:
co:=0:
co1:=0:
for i from 1 to N do
 R:=RandomBell1(K):
 co:=co+R:
 if evalf(R)>2 then
  co1:=co1+1:
 fi:
od:
[evalf(co/N),evalf(co1/N)]:
end:



####################################previous classes:
#C14.txt, March 10, 2025
#In honor of Albert Einstein's bday, born on 3/14
#Bell's theorem
with(combinat):
with(linalg):


Help14:=proc(): 
print(`Bell(A0,A1,B0,B1,V)`):
#Observe that V=ES()[1] for Bell
print(`Bell(C1QG()[4],C1QG()[2],-expand((C1QG()[2]+C1QG()[4])/sqrt(2)),expand((C1QG()[2]-C1QG()[4])/sqrt(2)),ES()[1]) `):
end:


#We have 4 observables: A0,A1,B0,B1
#Bell(A0,A1,B0,B1,V): the expression due to Bell, that is the quantum analog
# of a0*b0+a0*b1+a1*b0-a1*b1 (the expectation of its absolute value is <=2)
Bell:=proc(A0,A1,B0,B1,V):
ExpMV(TP(A0,B0),V)+ExpMV(TP(A0,B1),V)+ExpMV(TP(A1,B0),V)-ExpMV(TP(A1,B1),V):
end:



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


####################################################
#C8.txt, Feb. 17, 2025
Help8:=proc():
print(`ES(),TP(A,B),ExpMV(M,V)`):
end:


#ES(): the four fully entangled states in the Alice-Bob combined system
# based on pg. 166 of 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 Bob
# (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
# the resulting matrix (after applying tensor prod) is dim a1*b1 by a2*b2
TP:=proc(A,B) local a1,a2,b1,b2,i,j,i1,j1,AB1,AB:
a1:=nops(A):
a2:=nops(A[1]):
b1:=nops(B):
b2:=nops(B[1]):

AB:=[]:
#the rows of TP(A,B) are called [i,i1] where 1<=i<=a1 and 1<=i1<=b1
#the cols of TP(A,B) are called [j,j1] where 1<=j<=a2 and 1<=j1<=b2
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:


#Added after class:
#ExpMV(M,V): inputs  a Hermitian matrix M and a state vector V finds the expected value of the observable M in state V
#<conjugate(V),M,V>: 
ExpMV:=proc(M,V) local i,n,j:
add(add(conjugate(V[i])*M[i][j]*V[j],j=1..nops(V)),i=1..nops(V)):
end:

#################################################
#C6.txt, Feb. 10, 2025
Help6:=proc():
print(`Mul(A,B), Com(A,B),RSV(n,K), CheckUP(A,B,V)`):
end:

#Mul(A,B): the product of matrix A and B (assumming that it exists)
Mul:=proc(A,B) local i,j,k,n:
n:=nops(A):
[seq([seq( add(A[i][k]*B[k][j],k=1..n), j=1..n)],i=1..n)]:
end:


#Com(A,B): the commutator of A and B
Com:=proc(A,B):
Mul(A,B)-Mul(B,A):
end:


#RSV(n,K): A (normalized) random state vector with n complex components
# whose real and imaginary parts are from -K to K
RSV:=proc(n,K) local ra,v,i,a:
ra:=rand(-K..K):
v:=[seq(ra()+I*ra(),i=1..n)]:
a:=sqrt(add(v[i]*conjugate(v[i]),i=1..n)):
[seq(v[i]/a,i=1..n)]:
end:


#SOME NOTES:
# mu:=Sum(a*Pr(X=a),a in A):
# av=average
# sigma^2=Sum((a-av)^2*Pr(X=a),a in A) = Sum(a^2*Pr(), a in A)
# The standard deviation of a random variable is the square root of its variance

#NOTE: Com(PM()[1],PM()[2])=2*I*PM()[3]

#Del(L,A): the "average error", i.e. standard deviation of the prob. distribution
#  of measuring the observable L when the state is A
Del:=proc(L,A) local pd,ev,av,i,n:
n:=nops(L):
pd:=ProbDist(L,A):
ev:=pd[1]:
#the average (third element in ProbDist)
av:=pd[3]:
pd:=pd[2]:
#the first term in the following is the second moment
sqrt(add(ev[i]^2*pd[i],i=1..n)-av^2):
end:


#Here, we will change notation by making 
#CheckUP(A,B,V): the ratio of Del(A,V)*Del(B,V)/Expectation(Com(A,B),V)
# it should be >=1/2
CheckUP:=proc(A,B,V) local C:
C:=Com(A,B):

if abs(ProbDist(C,V)[3])=0 then
 RETURN(FAIL):
fi:

evalf(Del(A,V)*Del(B,V)/abs(ProbDist(C,V)[3])):
end:


#################################################
#C5.txt, Feb. 6, 2025
Help5:=proc(): 
print(` PM(), RHM(n,K) , EV(L), EVC1(L,g), EVC(L), IP(A,B), ProbDist(L,A) `):
end:


#PM(): the three famous Pauli matrices sigma_x,sigma_y, sigma_z
PM:=proc():
[
[[0,1],[1,0]],
[[0,-I],[I,0]],
[[1,0],[0,-1]],
[[1,0],[0,1]]
]:
end:



#RHM(n,K): a random n by n Hermitian matrix with entriries from -K to K
RHM:=proc(n,K) local ra,i,j,L:
ra:=rand(-K..K):
for i from 1 to n do
 for j from 1 to i-1 do
  L[i,j]:=ra()+I*ra():
  L[j,i]:=conjugate(L[i,j]):
 od:
L[i,i]:=ra():
od:
[seq([seq(L[i,j],j=1..n)],i=1..n)]:
end:



#EV(L): inputs a Hermitian matrix (say n by n) and outputs the LIST of eigenvalues
EV:=proc(L) local  x,L1,i,n:
n:=nops(L):
L1:=[seq(L[i]-[0$(i-1),x,0$(n-i)],i=1..n)]:
[solve(det(L1),x)]:
end:



#EVC1(L,g): inputs a matrix L and an eigenvalue g finds the NORMALIZED eigenvector
EVC1:=proc(L,g) local v,i,j,var,eq,V,n,L1,V1,va,k:
n:=nops(L):
L1:=[seq(L[i]-[0$(i-1),g,0$(n-i)],i=1..n)]:
if simplify(det(L1))<>0 then
  RETURN(FAIL):
fi:

V:=[seq(v[i],i=1..n)]:
var:={seq(v[i],i=1..n)}:
eq:={seq(add(L1[i,j]*V[j],j=1..n),i=1..n)}:
var:=solve(eq,var):
#V1 the set of unknowns that is arbitrary
V1:={}:
for va in var do
  if op(1,va)=op(2,va) then
  V1:=V1 union {op(1,va)}:
  fi:
od:
if nops(V1)<>1 then
  RETURN(FAIL):
fi:
V:=subs(var,V):
V:=subs(V1[1]=1,V):
k:=sqrt(add(abs(V[i])^2,i=1..n)):
[seq(V[i]/k,i=1..n)]:
end:



#EVC(L): inputs a Hermitian matrix and outputs the pair ListOfEigenvalues, ListOfEigenvectors
EVC:=proc(L) local n,ev,evc,i:
n:=nops(L):
ev:=EV(L):
evc:=[seq(EVC1(L,ev[i]),i=1..n)]:
[ev,evc]:
end:



#IP(A,B): the inner product of state-vectors A and B
IP:=proc(A,B) local n,i:
n:=nops(A):
if nops(B)<>n then
  RETURN(FAIL):
fi:
simplify(add(A[i]*conjugate(B[i]),i=1..n)):
end:



#ProbDist(L,A): Given an obervalble represented by the Hermitian matrix L and a state vector A
#outputs first all the eigenvalues, the prob. of winding in the respective pure states
#followed by the expected value of A
#(i.e. if you make the observation L and the state is A the average in the long run)
ProbDist:=proc(L,A) local evc,i,A1,k,PD,n:
n:=nops(L):
k:=sqrt(add(abs(A[i])^2,i=1..n)):
A1:=[seq(A[i]/k,i=1..n)]:
#evc outputs the pair ListOfEigenvalues and ListOfEigenvectors
evc:=EVC(L):
PD:=[seq(abs(IP(A1,evc[2][i]))^2,i=1..n)]:
#evc[1] is the list of eigenvalues
[evc[1],PD,add(PD[i]*evc[1][i],i=1..n)]:
end: