# OK to post homework
# Aurora Hiveley, 3/12/25, Assignment 14

Help := proc(): print(`RandomBell1(K), RandomBell(K,N)`): end:

with(combinat):
with(linalg):


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


## given:
ra:=proc(K): rand(-K..K)():end:
A0:=C1QG()[4]: A1:=C1QG()[2]: B0:= -expand((C1QG()[2]+C1QG()[4])/sqrt(2)): B1:= expand((C1QG()[2]-C1QG()[4])/sqrt(2)):


RandomBell1 := proc(K) local i,V:
V := [ra(K)+I*ra(K),ra(K)+ra(K)*I,ra(K)+I*ra(K),ra(K)+I*ra(K)]:
V := expand(V/ sqrt(add( V[i]*conjugate(V[i]), i=1..nops(V))) ):

Bell(A0,A1,B0,B1,V):
end:


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

RandomBell := proc(K,N) local tot,gr2,a,i:
tot := 0:
gr2 := 0:

for i from 1 to N do 
    a := evalf(RandomBell1(K)): 
    tot := tot + a:

    if a > 2 then 
        gr2 ++:
    fi:
od:

[tot/N, gr2/N]:
end:


# What did you get for RandomBell(10,1000)?
# output: [-0.019, 1/500]



### copied from C14.txt 


#C14.txt, March 10, 2025, Bell's Theorem
Help14:=proc(): print(`Bell(A0,A1,B0,B1), `):
print(`Bell(C1QG()[4],C1QG()[2], -(C1QG()[2]+C1QG()[4])/sqrt(2),(C1QG()[2]-C1QG()[4])/sqrt(2),ES()[1]); `): 
end:


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

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

#From C8.txt
#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:

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