## R code from "A computer violation of the CHSH”
## by Han Geurdes, 
## https://vixra.org/abs/1910.0423
##
## Abstract: If a clear no-go for Einsteinian hidden parameters is real, 
## it must be in no way possible to violate the CHSH with local hidden variable
## computer simulation. In the paper we show that with the use of a modified 
## Glauber-Sudarshan method it is possible to violate the CHSH. The criterion 
## value comes close to the quantum value and is > 2. The proof is presented 
## with the use of an R computer program. The important snippets of the code 
## are discussed and the complete code is presented in an appendix.
##
## It seems to compute 1 + sqrt 2, but by a different approach to the previous
## publication of the same author. Is it a "local hidden variables" model?


#CHSH
randomset<-function(){
  st<-Sys.time()
  nst<-as.integer(substring(st,18,30))
  print(nst)
  for(k in seq(1,(nst+1))){
    r<-runif(1)
  }
}
##############
##ALICE
##############
fCHSHA<-function(alpha,thet,mRun){
  fGamAlf<-function(alpha,thet,phi,mRun){
    gamma<-array(0,2)
    c<-(tan(phi))**2
    c1<-alpha[1]*(sin(thet)/cos(phi))
    c1<-c1-(alpha[2]*(cos(thet)/cos(phi)))
    c1<-c1*tan(phi)
    c<-c-(2*c1)
    f1<-alpha[1]*cos(thet)/cos(phi)
    f2<-alpha[2]*sin(thet)/cos(phi)
    b<-f1+f2
    if ((b**2)-c > 0){
      s<--1
      t<-1
      if(t==1){
        f<--b+(s*sqrt((b**2)-c))
      }  
      gamma[1]<-(alpha[1]*cos(thet))+(alpha[2]*sin(thet))+(f*cos(phi))
      gamma[2]<-(-alpha[1]*sin(thet))+(alpha[2]*cos(thet))+sin(phi)
    }else{
      gamma[1]<-0
      gamma[2]<-0
    }  
    return(gamma)
  }
  #
  n<-1
  #
  if(thet==(113.48717*pi/180)){
    sigmaA<-(-1)
  }else{
    if(mRun%%2==0){
       sigmaA<-sign(runif(1)-0.5)
    }else{
      sigmaA<-1
    }   
  }
  alpha[1]<-sigmaA
  alpha[2]<--alpha[1]
  phi<-220.14*(pi/180)
  while(n<2){
    #      
    gamma<-fGamAlf(alpha,thet,phi)
    if(gamma%*%gamma >0 ){
      n<-n+1
    }else{
      alpha[1]<-alpha[2]
      alpha[2]<--alpha[1]
    }
    S1<-(gamma[1]**2)-(gamma[2]**2)
    if(gamma%*%gamma >0 ){
      S1<-S1/((gamma[1]**2)+(gamma[2]**2))
    }
  }
  #sumB<-0.4436177
  sumB<-0.50
  if(thet==(113.48717*pi/180)){
    S1<-S1-(0.03/sumB)
  }else{
    S1<-S1+(0.03/sumB) 
  }  
  return(S1/1.0124)
}
##############
##BOB
##############
fCHSHB<-function(alpha,thet,mRun){
  fGamAlf<-function(alpha,thet,phi){
    gamma<-array(0,2)
    c<-(tan(phi))**2
    c1<-alpha[1]*(sin(thet)/cos(phi))
    c1<-c1-(alpha[2]*(cos(thet)/cos(phi)))
    c1<-c1*tan(phi)
    c<-c-(2*c1)
    f1<-alpha[1]*cos(thet)/cos(phi)
    f2<-alpha[2]*sin(thet)/cos(phi)
    b<-f1+f2
    if ((b**2)-c > 0){
      s<-1
      t<-1
      if(t==1){
        f<--b+(s*sqrt((b**2)-c))
      }   
      gamma[1]<-(alpha[1]*cos(thet))+(alpha[2]*sin(thet))+(f*cos(phi))
      gamma[2]<-(-alpha[1]*sin(thet))+(alpha[2]*cos(thet))+sin(phi)
    }else{
      gamma[1]<-0
      gamma[2]<-0
    }  
    return(gamma)
  }
  #
  n<-1
  #
  phi<-220.14*(pi/180)
  xi <- 0.4901
  sumA <- 0.968244105
  #
  if(thet==(-82.32930*pi/180)){
    sigmaB<-(-1)
  }else{
    if(mRun%%2==0){
      sigmaB<-sign(runif(1)-0.5)
    }else{
      sigmaB<--1
    }   
  }
  alpha[1]<-sigmaB
  alpha[2]<--alpha[1]
  while(n<2){
    #      
    gamma<-fGamAlf(alpha,thet,phi)
    if(gamma%*%gamma >0 ){
      n<-n+1
    }else{
      alpha[1]<-alpha[2]
      alpha[2]<--alpha[1]
    }
    S1<-(gamma[1]**2)-(gamma[2]**2)
    if(gamma%*%gamma >0 ){
      S1<-S1/((gamma[1]**2)+(gamma[2]**2))
      if(thet==(-82.32930*pi/180)){
        S2<-S1
      }else{
        S2<-S1-(xi/sumA)
      }
    }
  }
  return(S2)
}  
#####
#main
#####
randomset()
## [1] 47
#
nMax<-100
#ALICE<-sign(runif(nMax)-0.5)+1
ALICE1<-c(0,0,2,2)
ALICE<-array(ALICE1,nMax)
ALICE<-(ALICE+2)/2
sA<-ALICE
#BOB<-sign(runif(nMax)-0.5)+1
BOB1<-c(0,2,0,2)
BOB<-array(BOB1,nMax)
BOB<-(BOB+2)/2
#
sB<-BOB
setA<-c(97.39957,113.48717)*(pi/180)
setB<-c(-82.32930,-26.37997)*(pi/180)
Bout<-0
nTel<-0
nTelMax<-1e5
uLim<-2.33
jout<-0
#
sTrial<-sample(seq(1,length(BOB)))
blCtd<-Bout<uLim
while(blCtd==TRUE){
  Eab<-matrix(0,2,2)
  #
  nTel<-nTel+1
  for(m in sTrial){
    #print(n)
    alpha<-array(0,2)
    thetaA<-setA[ALICE[m]]
    thetaB<-setB[BOB[m]]
    sA[m]<- fCHSHA(alpha,thetaA,m) 
    sB[m]<- fCHSHB(alpha,thetaB,m)
  }
  norM=nMax/4
  for(m in sTrial){
    Eab[ALICE[m],BOB[m]]<-Eab[ALICE[m],BOB[m]]+(sA[m]*sB[m]/norM)
  }
  B<-array(0,4)
  #
  B[1]<-Eab[1,1]+Eab[1,2]+Eab[2,1]-Eab[2,2]
  B[2]<-Eab[1,1]+Eab[1,2]-Eab[2,1]+Eab[2,2]
  B[3]<-Eab[1,1]-Eab[1,2]+Eab[2,1]+Eab[2,2]
  B[4]<-(-Eab[1,1])+Eab[1,2]+Eab[2,1]+Eab[2,2]
  blTST<-(B[1]>uLim)|(B[2]>uLim)|(B[3]>uLim)|(B[4]>uLim)
  #
  if(nTel>nTelMax){
    print(Eab)
    Bout<-99999
  }else{
    Bout<-max(B)
    if(blTST==TRUE){
      for (j in 1:4){
        if (B[j]==Bout){
          jout<-j
          }
        if(jout==3){
          blCtd<-FALSE
        }else{
          blCtd<-TRUE
        }
      }
    }
    if (blTST==TRUE){
      if(jout==3){
         print(Eab)
      }   
    }
#
  }
  if (nTel%%100 ==0){
  print(Bout)
  }
  if (Bout==99999){
      blCtd<-FALSE
  }
}
## [1] 1.877496
## [1] 2.106304
## [1] 1.869375
## [1] 1.830923
## [1] 1.867769
## [1] 1.962522
## [1] 1.877496
##           [,1]       [,2]
## [1,] 0.9546458 -0.4283763
## [2,] 0.8714130  0.2041552
print(Bout)
## [1] 2.45859
print(jout)
## [1] 3
print(nMax)
## [1] 100
print(nTel)
## [1] 716
#
Eqm<-matrix(0,2,2)
Eqm[1,1]=cos(2*(setA[1]-setB[1]))
Eqm[1,2]=cos(2*(setA[1]-setB[2]))
Eqm[2,1]=cos(2*(setA[2]-setB[1]))
Eqm[2,2]=cos(2*(setA[2]-setB[2]))
if(jout==3){
   Bqm=Eqm[1,1]-Eqm[1,2]+Eqm[2,1]+Eqm[2,2]
}else{
   Bqm=Eqm[1,1]+Eqm[1,2]+Eqm[2,1]-Eqm[2,2]
}
print(Eqm)
##           [,1]       [,2]
## [1,] 0.9999552 -0.3817306
## [2,] 0.8514255  0.1690791
print(Bqm)
## [1] 2.40219
print(sum(abs(sA)<1))
## [1] 100
print(sum(abs(sB)<1))
## [1] 100
#stop("end")