For a version of Chantal's algorithm written in R, see
http://rpubs.com/gill1109/13494
The plot shows the sum of the true outcomes (+/-1) over the runs, which I took equal to 100 thousand per setting pair. Curiously its maximum is at 25 thousand and its minimum is at - 5 thousand.
(Actually I reverse one sign: my curve is upside down relative to Chantal's)
Doing this exercise I noticed some strange features in the Java code at
http://libertesphilosophica.info/eprsim/eprsim.txt
(Javascript version of Chantal's code by Daniel Sabsay).
Some of them duplicate similar odd features of Chantal Roth's Java code at
https://github.com/chenopodium/JCS2
others seem to be original.
Common to both simulations:
In each run, two particle pairs are generated, not one. But the pair is only counted once when reporting the number of runs.
Alongside these two, another pair is generated with equal and opposite angle between them.
The generation of uniform random points on the sphere is done in an extremely primitive way by Sabsay. It is computationally ineffient and it is inaccurate. He uses the multivariate normal method and his normal variables are sums of 6 uniforms.
As far as I can see Chantal computes the final correlation properly, dividing the sum of -1 and +1 outcomes by the number of runs. However Sabsay simply normalizes the raw sum of outcomes by scaling it to the interval [-1, 1]