# For a theoretical background, please see http://arxiv.org/abs/1405.2355
Angles = seq(from = 0, to = 360, by = 7.2) * 2 * pi/360
K = length(Angles) # The total number of angles between 0 and 2pi
corrs = matrix(nrow = K, ncol = K, data = 0) # Container for correlations
Ns = matrix(nrow = K, ncol = K, data = 0) # Container for events A and B
Ls = matrix(nrow = K, ncol = K, data = 0) # Container for initial states
# A canvas of auxiliary pre-ensemble (analogous to pre-reduced phase space):
M = 10^5 # Size of the pre-ensemble. Next one can try 10^6, or even 10^7
r = runif(M, 0, 2*pi) # M uniformly distributed numbers between 0 and 2pi
z = runif(M, -1, +1) # M uniformly distributed numbers between -1 and +1
h = sqrt(1 - z^2)
x = h * cos(r)
y = h * sin(r)
e = rbind(x, y, z) # A 3xM matrix, with M columns of e-vectors representing
# the x, y, z coordinates of uniformly distributed points on an S^2. In what
# follows these M number of pre-states will define an auxiliary pre-ensemble
# Defining the metrical and topological structures on S^3:
# So far we have set the variables in R^3. We now construct a metric on S^3
s = runif(M, 0, pi) # Initial states of the spins are the pairs (e, s) within S^3
f = -1 + (2/sqrt(1 + ((3 * s)/pi))) # For details see the paper arXiv:1405.2355
g = function(u,v,s){ifelse(abs(colSums(u*v)) > f, colSums(u*v), 0)}
# g(u,v,s) is non-vanishing only if |u.v| > f(s) ; g(u,v,s) --> "dot" as s --> pi
# Defines an inner product on S^3, thus changing the space from R^3 to S^3
# colSums(u * v) = u.v = the standard inner product between u and v in R^3
# u and v are orthogonal to each other in S^3 when abs(colSums(u * v)) < f
t = function(u,v,s){abs(sign(g(u,v,s))) > 0} # Ensures sign(g) = +1 or -1
# The pair {g, t} defines together a metric on S^3 constraining its topology
n = function(Q,u1,v1,s1,u2,v2,s2){length(Q[t(u1,v1,s1) & t(u2,v2,s2)])}
# Calculates the number of occurrences of 'Q' while respecting the metric {g, t}
# Computing the "quantum" correlations:
i <- 10
alpha = Angles[i]
a = c(cos(alpha), sin(alpha), 0) # Measurement direction 'a'
j <- 20
beta = Angles[j]
b = c(cos(beta), sin(beta), 0) # Measurement direction 'b'
A = +sign(g(a,e,s)) # Alice's measurement results A(a, e, s) = +/-1
B = -sign(g(b,e,s)) # Bob's measurement results B(b, e, s) = -/+1
Cuu = length((A*B)[A > 0 & B > 0]) # Coincidence count of (+,+) events
Cdd = length((A*B)[A < 0 & B < 0]) # Coincidence count of (-,-) events
Cud = length((A*B)[A > 0 & B < 0]) # Coincidence count of (+,-) events
Cdu = length((A*B)[A < 0 & B > 0]) # Coincidence count of (-,+) events
corrs[i,j] = (Cuu + Cdd - Cud - Cdu) / (Cuu + Cdd + Cud + Cdu)
N = n((A*B),a,e,s,b,e,s) # Total number of simultaneous events observed
o = x[t(a,e,s) & t(b,e,s)]
p = y[t(a,e,s) & t(b,e,s)]
q = z[t(a,e,s) & t(b,e,s)]
w = rbind(o,p,q) # N vectors in R^3 representing the initial state (e, s) in S^3
a.w = colSums(a*w) # Standard inner product in R^3 between the vectors a and w
b.w = colSums(b*w) # Standard inner product in R^3 between the vectors b and w
Ns[i,j] = N # Total number of simultaneous +/- events observed by Alice and Bob
Ls[i,j] = n(s,a,e,s,b,e,s) # The number of initial states w = (e, s) within S^3
(Cou = length((A*B)[A == 0 & B > 0])) # Number of (0,+) events within S^3
## [1] 8764
(Cod = length((A*B)[A == 0 & B < 0])) # Number of (0,-) events within S^3
## [1] 8862
(Cuo = length((A*B)[A > 0 & B == 0])) # Number of (+,0) events within S^3
## [1] 8797
(Cdo = length((A*B)[A < 0 & B == 0])) # Number of (-,0) events within S^3
## [1] 8833
(Coo = length((A*B)[g(a,e,s) & A == 0 & B == 0])) # Number of (0,0) events
## [1] 0
(CoB = length(A[g(a,e,s) & A == 0])) # Number of A = 0 events within S^3
## [1] 0
(CAo = length(B[g(b,e,s) & B == 0])) # Number of B = 0 events within S^3
## [1] 0
(W = length(w)/3) # Total number of w vectors
## [1] 49091
(N) # = length((A*B)[abs(A) > 0 & abs(B) > 0])
## [1] 49091
(J = length((A*B)[abs(A*B) > 0]))
## [1] 49091