# OK to post homework
# Aurora Hiveley, 2/17/25, Assignment 8

Help := proc(): print(`EntangelmentAlice(V), RandomProductState4(K), EntanglementBob(V)`): end:

with(combinat):
with(linalg):

# Read and understand procedure ExpMV(M,V) that is supposed to give you the same 
# output as ProbDist(M,V)[3]. Take three take a random 3 by 3 Hermitian matrices 
# (using RHM(3,10)) and a random 3-state vetor (using RSV(3,10)) and check that 
# if you take evalf(ProbDist(M,V)[3]) you get the same thing as evalf(ExpMV(M,V)). 
# Do it three times.

# M := RHM(3,10):
# V := RSV(3,10):
# print([evalf(ProbDist(M,V)[3]), evalf(ExpMV(M,V))]);

# sure enough, they're the same! (although the prob dist has a basically irrelevant imaginary component
# on the order of e-9 or e-10)


# By taking V:=RSV(2,100); verify that add(ExpMV(PM()[i],V)^2,i=1..3)=1. Do it three times. 
# In other words the sum of the squares of the expectation of the three Pauli matrices 
# adds up to 1 for every state vector.

# V:=RSV(2,100):
# add(ExpMV(PM()[i],V)^2,i=1..3);     # expect = 1

# sure enough, we get a sum of 1, yay!



# The Alice-versions of the Pauli matrices are obtained by taking the tensor product 
# with the identity, getting a 4 by 4 Hermitian matrix, in other words TP(PM()[i],PM()[4])
# Let's call them Ax,Ay,Az.

# EntangelmentAlice(V): For a random state vector (obtained by RSV(4,10), say), 
# outputs the sum of the squares of the expected values of Alice's measurement in that state.

EntanglementAlice := proc(V) local i:
add(ExpMV( TP(PM()[i],PM()[4]), V)^2,i=1..3):
end:



# RandomProductState4(K): outputs a random product state vector of length 4 
# obtained by taking the tensor product of two RSV(2,K)

RandomProductState4 := proc(K):
op(TP([RSV(2,K)],[RSV(2,K)])):
end:



# what are EntangelmentAlice(V) for ES()[i], i=1..4?
# [seq(EntanglementAlice(ES()[i]), i=1..4)];
# output vector is [0,0,0,0] as expected 


# What is it for a random product state?
# EntanglementAlice(RandomProductState4(3));
# output is 1, cool!


# By using RSV(4,10) twenty times, find the vector with maximal entanglemt 
# (i.e. the smallest value of EntanglemtAlice(V)) and the least (closest to 1)

# L:=[]:

# for i from 1 to 20 do 
#     V := RSV(4,10):
#     L := [op(L), [EntanglementAlice(V), V]]:
# od:

# M := min(seq(L[i][1], i=1..20));

# for l in L do 
#     if l[1] = M then 
#         print(l);
#     fi:
# od:

# output: [4904/38809, [(-7/394 - 2*I/197)*sqrt(394), (5/197 + 4*I/197)*sqrt(394), (4/197 + 3*I/197)*sqrt(394), (1/394 + 4*I/197)*sqrt(394)]]



# EntangelmentBob(V): For a random state vector (obtained by RSV(4,10), say), 
# outputs the sum of the squares of the expected values of Bob's measurement in that state.

EntanglementBob := proc(V) local i:
add(ExpMV( TP(PM()[i],[[0,1],[1,0]]), V)^2,i=1..3):  # not sure how to calculate Bx,By,Bz, tried to adapt basis vectors??
end:








### copied from C8.txt 
#C8.txt, Feb. 17, 2025
#Feb. 10, 2025 C6.txt
Help8:=proc(): print(`ES(), TP(A,B), ExpMV(M,V) `):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:



#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
#<conjuate(V),M,V>: 
ExpMV:=proc(M,V) local i,j,n:
add(add(conjugate(V[i])*M[i][j]*V[j],j=1..nops(V)),i=1..nops(V)):
end:


#old stuff
Help6:=proc(): print(`Mul(A,B), Com(A,B), RSV(n,K), Del(L,A), CheckUP(A,B,V) `): end:

#Mul(A,B): the product of matrix A and B (assuming 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):
RSV:=proc(n,K) local ra,i,v,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:

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

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

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

with(linalg):

#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 eigenvectors (pure states) and 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:=EVC(L):
PD:=[seq(abs(IP(A1,evc[2][i]))^2,i=1..n)]:
[evc[1],PD,add(PD[i]*evc[1][i],i=1..n)]:
end: