```python
from pylab import *
import sympy as sym
from sympy.physics.quantum import Commutator, Dagger, Operator
from sympy.physics.quantum import Bra, Ket, InnerProduct
import numpy as np
init_printing(use_unicode=True)
```


```python
def spinx(s):

    n = int(2.0*s+1)
    sx = np.matrix(zeros((n,n)))
    for a in range(0,n):
        for b in range(0,n):
            if (a==b+1):
                sx[a,b] = sx[a,b] + 0.5*sqrt((s+1)*(a+b+1)-(a+1)*(b+1))
            elif (a==b-1):
                sx[a,b] = sx[a,b] + 0.5*sqrt((s+1)*(a+b+1)-(a+1)*(b+1))
                
    return sx


def spiny(s):
    n = int(2.0*s+1)
    sy = np.matrix(zeros((n,n),dtype='complex'))

    for a in range(0,n):
        for b in range(0,n):
            if (a==b+1):
                sy[a,b] = sy[a,b] + 0.5j*sqrt((s+1)*(a+b+1)-(a+1)*(b+1))
            elif (a==b-1):
                sy[a,b] = sy[a,b] - 0.5j*sqrt((s+1)*(a+b+1)-(a+1)*(b+1))
                
    return sy
    
def spinz(s):
    
    n = int(2.0*s+1)
    sz = np.matrix(zeros((n,n)))

    for a in range(0,n):
        for b in range(0,n):
    
            if (a==b):
                sz[a,b] = (s+1-b-1)
                
    return sz
```


```python
s = 0.5
Sx = spinx(s)
Sy = spiny(s)
Sz = spinz(s)
```


```python
ahat = np.matrix([0,1,0])
bhat = np.matrix([np.cos(np.pi/6),np.sin(-np.pi/6),0])
chat = np.matrix([np.cos(-np.pi/6),np.sin(-np.pi/6),0])
print(ahat,bhat,chat,norm(ahat),norm(bhat),norm(chat))
```

    [[0 1 0]] [[ 0.8660254 -0.5        0.       ]] [[ 0.8660254 -0.5        0.       ]] 1.0 1.0 1.0



```python
Nm = 1000000 #number of measurements
Na = 3 #number of axes

thet = linspace(0,(Na-1.0)*2.0*pi/Na,Na) #array of three values of theta
phi = zeros(3) #set phi = 0 for all three axes

#singlet in z-basis
psi0 = (1.0/sqrt(2.0))*np.matrix([0,1.0,-1.0,0])
```


```python
#single spin operators
Ss = zeros((2,2,Na),dtype='complex')

for i in range(0,Na):
    Ss[:,:,i] = 0.5*matrix([[cos(thet[i]),sin(thet[i])*exp(-1j*phi[i])],[sin(thet[i])*exp(1j*phi[i]),-cos(thet[i])]])
    
#two-spin operators
Sab = zeros((2**2,2**2,Na**2),dtype='complex')
Vab = zeros((2**2,2**2,Na**2),dtype='complex') #eigenvectors
Eab = zeros((2**2,Na**2),dtype='complex')      #eigenvalues

for i in range(0,Na):
    for j in range(0,Na):
        Sab[:,:,i*Na+j] = kron(Ss[:,:,i],Ss[:,:,j])
        Eab[:,i*Na+j],Vab[:,:,i*Na+j] = eigh(Sab[:,:,i*Na+j])
```


```python
#MEASUREMENT FUNCTION

def doublemeasurement(Na,psi0,Vab,Eab):
    #select pair of axes at random
    m = randint(0,Na**2)

    #coefficients of psi in measurement basis
    psi = zeros(2**2,dtype='complex')

    for i in range(0,2**2):
        psi[i] = dot(conj(Vab[:,i,m]),psi0)

    r = rand()
    csum = abs(psi[0])**2
    
    n = 0
    while(r>csum):
        n = n + 1
        csum = csum + abs(psi[n])**2

    #MEASURED STATE is indexed by n;
    #returns sign of the product (+ -> same, - -> different)
    return sign(Eab[n,m])
```


```python
#counters
Nsame = 0
Ndiff = 0

#repeat measurement Nm times
for i in range(0,Nm):
    #if
    if (doublemeasurement(Na,psi0,Vab,Eab)>0):
        Nsame = Nsame + 1
        
    else:
        Ndiff = Ndiff+1
        
print('The probability of same measurement is ' + str(Nsame/Nm))
print('The probability of different measurement is ' + str(Ndiff/Nm))

if Nsame>4.0*Nm/9.0:
    print('Bell inequality VIOLATED!')
else:
    print('Bell inequality satisfied.')
```

    The probability of same measurement is 0.500589
    The probability of different measurement is 0.499411
    Bell inequality VIOLATED!


56c706c827a5205c33682eab163dcc0d74e835e8
ipynb
Jupyter Notebook
physics/quant/BellTest.ipynb
soundofives/hpu-classes
8c690b64e2de1ff7ffadf7aa9aabbfe634cf36a8
MIT
physics/quant/BellTest.ipynb
soundofives/hpu-classes
8c690b64e2de1ff7ffadf7aa9aabbfe634cf36a8
MIT
physics/quant/BellTest.ipynb
soundofives/hpu-classes
8c690b64e2de1ff7ffadf7aa9aabbfe634cf36a8
MIT
Qwen/Qwen-72B
1. YES
2. YES