#Ok to post

#Salman Manzoor, 03-13-25, Assignment 14



RandomBell1:=proc(K) local A0, A1, B0, B1, V, k, i:

A0:=C1QG()[4]: A1:=C1QG()[2]: B0:= -(C1QG()[2]+C1QG()[4])/sqrt(2): B1:=(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:=add(abs(V[i]),i=1..4):

Bell(A0,A1,B0,B1,V/k):

end:



RandomBell:=proc(K,N) local sum, count, instance:

sum:=0:

count:=0:

for i from 1 to N do

	instance:=RandomBell1(K):

	sum:=sum+instance:

	if evalb(instance>2) then count++: end:

od:

evalf([sum,count]/N):

end:



ra:=proc(K): rand(-K..K)(): 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:



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: