This is a little test of Joy Christian’s one-page paper: http://arxiv.org/pdf/1103.1879v1.pdf

I use the R package “onion” to implement the quaternions, through which I am able to define the various bivector algebras used by Joy Christian in his paper.

First of all we load the library “onion” and we define a function for sampling points on S2, uniformly at random. We are going to test Christian’s model by picking two measurement directions a, b, completely at random

library(onion)


S2 <- function(N) {    
    z <- runif(N, -1, 1)
    t <- runif(N, 0, 2*pi)
    r <- sqrt(1 - z^2)
    x <- r * cos(t)
    y <- r * sin(t)
    rbind(x, y, z)
}

set.seed(1234)
a <- S2(1)
b <- S2(1)

a
##         [,1]
## x -0.4564959
## y -0.4412610
## z -0.7725932
b
##         [,1]
## x -0.6970031
## y -0.6829515
## z  0.2185495

Next we will create the bivector basis elements and test that they have the desired properties; Christian’s formula (3).

Note that in “onion”, H1, Hi, Hj, Hk stand for the quaternion (real) unit and the three quaternionic roots of minus one. To get the bivector algebra, we simply have to multiply Hi, Hk and Hk each by minus 1.

## test 

beta1 <- -Hi
beta2 <- -Hj
beta3 <- -Hk

(beta1 * beta2 == -beta3)
## [1] TRUE
(beta2 * beta3 == -beta1)
## [1] TRUE
(beta3 * beta1 == -beta2)
## [1] TRUE
(beta2 * beta1 == beta3)
## [1] TRUE
(beta3 * beta2 == beta1)
## [1] TRUE
(beta1 * beta3 == beta2)
## [1] TRUE

Now we will also test the two derived algebras generated through choice of lambda = +1 or -1, in particular, testing Christian’s formulas (1) and (2).

My “abeta” is Christian’s \(\sum a_j\beta_j\) (note: he uses Einstein summation convention in his paper; here I write out the summartion).

My “abetaL” is Christian’s \(\sum a_k\beta_k(\lambda)\).

My “bbeta” is Christian’s \(\sum b_k\beta_k\).

My “bbetaL” is Christian’s \(\sum b_j\beta_j(\lambda)\).

I do not explicitly test (4) but the reader is welcome to do that also. Christian’s new bivector basis elements \(\beta_j(\lambda)\) are equal to my “betaL1”, “betaL2”, “betaL3.”

lambda <- 1
betaL1 <- beta1 * lambda
betaL2 <- beta2 * lambda
betaL3 <- beta3 * lambda
abeta <- a[1]*beta1 +a[2]*beta2 +a[3]*beta3
abetaL <- a[1]*betaL1 +a[2]*betaL2 +a[3]*betaL3
bbeta <- b[1]*beta1 +b[2]*beta2 +b[3]*beta3
bbetaL <- b[1]*betaL1 +b[2]*betaL2 +b[3]*betaL3

all.equal((-abeta) * abetaL, lambda*H1)
## [1] TRUE
all.equal((bbetaL) * bbeta, -lambda*H1)
## [1] TRUE
lambda <- -1
betaL1 <- beta1 * lambda
betaL2 <- beta2 * lambda
betaL3 <- beta3 * lambda
abeta <- a[1]*beta1 +a[2]*beta2 +a[3]*beta3
abetaL <- a[1]*betaL1 +a[2]*betaL2 +a[3]*betaL3
bbeta <- b[1]*beta1 +b[2]*beta2 +b[3]*beta3
bbetaL <- b[1]*betaL1 +b[2]*betaL2 +b[3]*betaL3

all.equal((-abeta) * abetaL, lambda*H1)
## [1] TRUE
all.equal((bbetaL) * bbeta, -lambda*H1)
## [1] TRUE

Now we will test formulas (5) through (7). I only do 1000 samples, but the reader is welcome to test with a much larger number. For that purpose it would be useful to “vectorise” the code; i.e., do not do an explict “for” loop, but instead create a long vector of values of lambda and use R’s built-in vector and matrix algebra.

Here, I want to the code to be as close as possible to the formulas in the paper.

Nrep <- 1000
S <- 0 * H1
for (i in 1:1000){
  lambda <- sample(c(-1, 1), 1)
  betaL1 <- beta1 * lambda
  betaL2 <- beta2 * lambda
  betaL3 <- beta3 * lambda
  abeta <- a[1]*beta1 +a[2]*beta2 +a[3]*beta3
  abetaL <- a[1]*betaL1 +a[2]*betaL2 +a[3]*betaL3
  bbeta <- b[1]*beta1 +b[2]*beta2 +b[3]*beta3
  bbetaL <- b[1]*betaL1 +b[2]*betaL2 +b[3]*betaL3
  
  A <- (-abeta) * abetaL
  B <- bbetaL * bbeta
  S <- S + A * B
}

(abeta)^{-1} * (S/Nrep) * bbeta   ## E(a, b)
##           [1]
## Re -0.4506891
## i  -0.6240811
## j   0.6382668
## k   0.0042043

We did not get a real number, and certainly not the real number a.b

What did we get? I calculate the bivectorial product a b and the real number a.b, so we can see what has happened.

abeta * bbeta
##           [1]
## Re -0.4506891
## i  -0.6240811
## j   0.6382668
## k   0.0042043
sum(a * b)
## [1] 0.4506891