#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>

// k-nearest neighbor Monte Carlo high confidence bounds ---------------------

#define	K	4	// number of nearest neighbors
#define	X	'B'	// type of component U/O/B
#define BND	'L'	// upper/lower bound U/L
#define R	50.	// radius of central disk
#define	S	0.	// buffer zone around central disk
#define ITER	100	// number of iterations

#define	M	100000	// maximum number of points in rectangle
#define SEED	1	// random number generator seed

// For simplicity, mostly use global variables -------------------------------

double vx[M],vy[M];	// coordinates of points
int st[M];		// status 0-normal 1-left 2-right 3-done
int pc[M];		// vertex is in cluster or comes from boundary
int qq[M],qs,qe;	// queue of vertices
int nb[M][6*K+1];	// neighbor list
int nbb[M][K];		// neighbor backup list
double ds[K+2];		// list of distances^2 from given vertex
int nt[K+2];		// temp neighbor list
int nv;			// total number of vertices
int cl,cr;		// left/right vertex count
int sl,sr;		// left/right success count
int erc;		// error count (too many vertices in rectangle)

// Random number generator ---------------------------------------------------

unsigned char rs[256],ri,rj;		// random number generator status

// Random number, uniform in [0,1] (never =0)
double rand_uniform(){
 unsigned int s,t,i;
 double w;
 for(w=.5,i=0;i<7;i++){			// generates 56 bits of random data
  ri=(ri+1)&0xFF;rj=(rj+rs[ri])&0xFF;	// &0xFF to fix compiler bug
  s=rs[ri];rs[ri]=t=rs[rj];rs[rj]=s;
  w=(rs[(s+t)&0xFF]+w)/256;		// rs[(s+t)&0xFF]=random byte
 }
 return w;				// w=0.{56 random bits}1 in binary
}

// Initialize random number generator
// Call before first use of rand_uniform()
void rand_setup(){
 unsigned int s,t;
 int i;
 ri=1; rj=SEED;				// initialize with seed
 for(i=0;i<256;i++)rs[i]=i;		// setup permutation array
 for(i=0;i<64;i++)rand_uniform();	// run algorithm to scramble rs[]
}

// Generate instance ---------------------------------------------------------

// generate graph and modify edges according to X = U/O/I/B
//
int generate(){
 double x,y,d;
 int i,j,k,t,w;
// first generate the points
 nv=0; x=-2*(R+S);
 cl=cr=0;
 for(;x<2*(R+S) && nv<M;nv++){
  x-=log(rand_uniform())/(2*(R+S)); vx[nv]=x;
  vy[nv]=y=(2*rand_uniform()-1)*(R+S);
  st[nv]=pc[nv]=0;
  if((x+(R+S))*(x+(R+S))+y*y<R*R){cl++;st[nv]=1;} // in C_1
  if((x-(R+S))*(x-(R+S))+y*y<R*R){cr++;st[nv]=2;} // in C_2
 }
 if(nv==M){erc++;return 0;}		// too many points
// now find nearest neighbors
 for(i=0;i<nv;i++){
  x=2*(R+S)-fabs(vx[i]);y=(R+S)-fabs(vy[i]); // fake vertices on boundary
  if(x>y)x=y; d=x*x;
  for(j=0;j<K;j++){nb[i][j]=-1;ds[j]=d;}
  nb[i][K]=-1;
  for(j=i-1;j>=0;j--){
   x=vx[i]-vx[j]; if(x*x>=ds[K-1])break;
   y=vy[i]-vy[j]; d=x*x+y*y;
   if(d<ds[K-1]){
    for(k=K-1;k>0 && ds[k-1]>d;k--){
     nb[i][k]=nb[i][k-1]; ds[k]=ds[k-1];
    }
    nb[i][k]=j; ds[k]=d;
   }
  }
  for(j=i+1;j<nv;j++){
   x=vx[i]-vx[j]; if(x*x>=ds[K-1])break;
   y=vy[i]-vy[j]; d=x*x+y*y;
   if(d<ds[K-1]){
    for(k=K-1;k>0 && ds[k-1]>d;k--){
     nb[i][k]=nb[i][k-1]; ds[k]=ds[k-1];
    }
    nb[i][k]=j; ds[k]=d;
   }
  }
 }
 if(X=='U'||X=='I'){		// undirected / in-component case
  for(i=0;i<nv;i++){
   for(k=0;k<K;k++)nbb[i][k]=nb[i][k]; // copy lists
   if(nb[i][K-1]<0)pc[i]=1;	// comes from boundary
  }
  if(X=='I'){
   for(i=0;i<nv;i++)nb[i][0]=-1; // clear out-going edges for I
  }
  for(i=0;i<nv;i++){
   for(k=0;k<K && nbb[i][k]>=0;k++){
    j=nbb[i][k];
    for(t=0;nb[j][t]>=0 && nb[j][t]!=i;t++);
    if(nb[j][t]<0){
     if(t>=6*K){printf("\nFatal error - too many neighbors\n");exit(1);}
     nb[j][t]=i;
     nb[j][t+1]=-1;
    }
   }
  }
 }
 if(X=='B'){			// bidirectional case
  for(i=0;i<nv;i++){
   if(nb[i][K-1]<0)pc[i]=2;	// possibly from boundary
   for(k=w=0;nb[i][k]>=0;k++){
    j=nb[i][k];
    for(t=0;nb[j][t]>=0 && nb[j][t]!=i;t++);
    if(nb[j][t]==i){nb[i][w]=j;w++;}
   }
   nb[i][w]=-1;
  }
 }
// check for incoming edges from boundary
// (only lower bound U/O/B cases)
 if(BND!='U' && X!='I'){
  if(nv<K+2){erc++;return 0;}	// need at least K+2 vertices
  x=-2*(R+S); y=-(R+S);		// start at corner
  while(1){			// creep round boundary

//printf("(%f,%f)",x,y);

   for(k=0;k<K+2;k++)ds[k]=1e99;
   for(i=0;i<nv;i++){		// find nearest K+2 vertices
    d=(x-vx[i])*(x-vx[i])+(y-vy[i])*(y-vy[i]);
    if(d<ds[K+1]){
     for(k=K+1;k>0 && ds[k-1]>d;k--){
      nt[k]=nt[k-1]; ds[k]=ds[k-1];
     }
     nt[k]=i; ds[k]=d;
    }
   }
// we mark the nearest K+1 points
// and move (d_{k+2}-d_K)/2 around boundary
// this guarantees we get K nearest points of any
// boundary point.
   for(k=0;k<K+1;k++)pc[nt[k]]|=1;
   d=(sqrt(ds[K+1])-sqrt(ds[K-1]))/2;
   if(y==-(R+S)){x+=d;d=0;}			// bottom
   if(x>=2*(R+S)){y+=d+x-2*(R+S);x=2*(R+S);d=0;} // right
   if(y>=(R+S)){x-=d+y-(R+S);y=R+S;d=0;}	// top
   if(x<=-2*(R+S)){y-=d-x-2*(R+S);x=-2*(R+S);d=0;} // left
   if(y<-(R+S))break;
  }
  if(X=='B'){
   for(i=0;i<nv;i++)pc[i]=(pc[i]==3); // need both
  }
 }

//for(j=i=0;i<nv;i++)j+=pc[i];
//printf("%d/%d from boundary\n",j,i);

 return 1;
}



// Upper bound trial ---------------------------------------------------------

// Test whether most points in central disk (radius R)
// inside left square (side 2R+2S) -> most points
// in central disk on right.
//
int ubtrial(){
 int i,j,k,t,w,b;
// generate graph
 if(!generate())return 0;
// Now percolate from each vertex
 for(b=0;b<2;b++){			// both directions
  for(i=0;i<nv;i++)pc[i]=0;
  sl=0; w=0;
  for(i=0;i<nv;i++)if(st[i]==1){	// pick vertex in C_1
   sr=0; w++; pc[i]=w;			// vertex identifier
   qs=0; qe=1; qq[0]=i;			// setup queue
   while(qe>qs){j=qq[qs++];
    for(k=0;nb[j][k]>=0;k++){
     if(pc[nb[j][k]]<w){
      if(qe>=nv){printf("Fatal queue error\n"); exit(1);}
      qq[qe++]=nb[j][k];
      pc[nb[j][k]]=w;
     }
    }
    if(st[j]==2)sr++;
   }
   if(X=='U' || X=='B'){
    t=0;
    for(j=0;j<nv;j++){
     if(st[j]==1 && pc[j]==w){t++; st[j]=3;}
    }
   }
   else {st[i]=3;t=1;}
   if(sr*2>cr)sl+=t;
  }
  if(2*sl<=cl)return 0;
  for(i=0;i<nv;i++){				// swap L<->R
   if(st[i])st[i]--;
  }
  i=cr;cr=cl;cl=i;
 }
 return 1;
}

// Lower bound trial ---------------------------------------------------------

// Test whether there exists path from boundary cutting
// line between centers of squares
//
int lbtrial(){
 int i,j,k;
// generate graph
 if(!generate())return 0;
// percolate from boundary
 qs=0;qe=0;
 for(i=0;i<nv;i++)if(pc[i])qq[qe++]=i;
 while(qe>qs){i=qq[qs++];
  for(k=0;nb[i][k]>=0;k++){
   j=nb[i][k];
   if(vy[i]*vy[j]<0){		// test if crosses center line
    if(fabs((vx[i]*vy[j]-vx[j]*vy[i])/(vy[j]-vy[i]))<=R+S)return 0;
   }
   if(!pc[j]){
    if(qe>=nv){printf("Fatal queue error\n"); exit(1);}
    qq[qe++]=j; pc[j]=1;
   }
  }
 }
 return 1;
}

// Compute statistics --------------------------------------------------------

main(){
 time_t start,end;
 double f,s;
 int c,t,i;

 time(&start);
 rand_setup();
 erc=c=0;
 for(t=0;t<ITER;t++){
  if(BND=='U')c+=ubtrial(); else c+=lbtrial();
  if(!(31&~t)){printf(".");fflush(NULL);}
 }
 time(&end);

 printf("\n%c k=%d r=%f s=%f: %d/%d successes\n",X,K,R,S,c,ITER);

 f=1;s=0;t=0;
 for(i=0;i<=ITER-c;i++){
  s=s*0.8639+f;
  f=f*(ITER-i)/(i+1)*0.1361;
  if(f>1){f/=1000;s/=1000;t+=3;}
 }
 for(;i<=ITER;i++){
  s*=0.8639;
  if(s<1){s*=1000;t-=3;}
 }
 while(s>=10-5e-7){s/=10;t++;}
 while(s<1 && s>0){s*=10;t--;}
 printf("Confidence=%fE%d\n",s,t);
 if(erc)printf("%d errors due to too many points\n",erc);
 printf("%d seconds\n",(int)(end-start));
}