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

const int N = 20;
const int L = 2*N;

struct Edge
{
  int i;
  int j;
};


void initE(Edge E[L]);
void prE(Edge E[L], char* fname);
int newp(int p[N+L], int n);
int set(int m[N], int e[L], int pa, Edge E[L], int R[N][N]);
void init(int p[N+L], int m[N], int e[L], int R[N][N]);
void prR(int m[N], Edge E[L], int e[L], char* fname);
int ck(int m[N], int R[N][N]);
void explore(int i, int m[N], int R[N][N]);
int recv(int p[N+L], int m[N], int e[L], int R[N][N], Edge E[L]);
double rate(Edge E[L], int suc[N+L], int T);
void torus3x3(Edge E[L]);
void K_ring_6de(Edge E[L]);
void K_ring_4de(Edge E[L]);

main()
{
  int ti=time(0);
  printf("srand=%d\r\n",ti);
  srand(ti);

  Edge E[L];


  int suc[N+L+1];

  int T=50000;

  char fname[1000];
  sprintf(fname,"out-N=%d-L=%d.csv",N,L);
  FILE *f=fopen(fname,"w");
  fprintf(f,",P\r\n");

  double maxS=0;
  double minS=1e10;
  for(int t=0;t<10000;t++)
  {
    if(t==0)
      K_ring_4de(E);
    else
      initE(E);

    double S=rate(E,suc,T);
    printf("%2d, %f\r\n",t,S);

    if(S<minS || S>maxS)
    {
      for(int i=N;i<=N+L;i++)
        fprintf(f,"%2d, %f\r\n",i,(double)suc[i]/T);
      fprintf(f,",\r\n");
    }
    if(S<minS)
    {
      minS=S;
      printf("saving the min graph\r\n");
      prE(E,"out-minS.eps");
    }
    if(S>maxS && t>0)
    {
      maxS=S;
      printf("saving the max graph\r\n");
      prE(E,"out-maxS.eps");
    }

  }

  fclose(f);
}

void K_ring_4de(Edge E[L])
{
  if(L!=2*N)
  {
    initE(E);
    return;
  }

  int G[N][N];

  for(int i=0;i<N;i++)
    for(int j=0;j<N;j++)
      G[i][j]=0;

  for(int i=0;i<N;i++)
  {
    G[i][(N+i+1)%N]=1;
    G[i][(N+i-1)%N]=1;

    G[i][(N+i+2)%N]=1;
    G[i][(N+i-2)%N]=1;
  }

  int l=0;
  for(int i=0;i<N-1;i++)
    for(int j=i+1;j<N;j++)
      if(G[i][j])
      {
        E[l].i=i;
        E[l].j=j;
        l++;
      }

}

void K_ring_6de(Edge E[L])
{
  if(!(N==15 && L==45))
  {
    initE(E);
    return;
  }

  int G[N][N];

  for(int i=0;i<N;i++)
    for(int j=0;j<N;j++)
      G[i][j]=0;

  for(int i=0;i<N;i++)
  {
    G[i][(N+i+1)%N]=1;
    G[i][(N+i-1)%N]=1;

    G[i][(N+i+2)%N]=1;
    G[i][(N+i-2)%N]=1;

    G[i][(N+i+3)%N]=1;
    G[i][(N+i-3)%N]=1;
  }

  int l=0;
  for(int i=0;i<N-1;i++)
    for(int j=i+1;j<N;j++)
      if(G[i][j])
      {
        E[l].i=i;
        E[l].j=j;
        l++;
      }

}

void torus3x3(Edge E[L])
{
  if(!(N==9 && L==18))
  {
    initE(E);
    return;
  }

  int G[9][9];

  for(int i=0;i<9;i++)
    for(int j=0;j<9;j++)
      G[i][j]=0;

  for(int i=0;i<3;i++)
    for(int j=0;j<3;j++)
    {
      G[i*3+j][3*( i       ) + (j+1)%3]=1;
      G[i*3+j][3*( (i+1)%3 ) + j      ]=1;
      G[i*3+j][3*( i       ) + (j+2)%3]=1;
      G[i*3+j][3*( (i+2)%3 ) + j      ]=1;
    }

  int l=0;
  for(int i=0;i<8;i++)
    for(int j=i+1;j<9;j++)
      if(G[i][j])
      {
        E[l].i=i;
        E[l].j=j;
        l++;
      }
}

double rate(Edge E[L], int suc[N+L+1], int T)
{
  int p[N+L];
  int m[N];
  int e[L];
  int R[N][N];

  for(int i=0;i<=N+L;i++) suc[i]=0;

  for(int t=0;t<T;t++)
  {
    int n=recv(p,m,e,R,E);
    for(int i=n;i<=N+L;i++)
      suc[i]++;
  }

  double S=0;
    for(int j=0;j<=N+L;j++)
      S+=(double)suc[j]/T;

  return S;
}

int recv(int p[N+L], int m[N], int e[L], int R[N][N], Edge E[L])
{
  init(p,m,e,R);

  int n;
  for(n=1;n<=N-1;n++)
  {
    set(m,e,newp(p,n),E,R);
  }

  for(n=N;n<=N+L;)
  {
     while(!set(m,e,newp(p,n++),E,R))
       ;

     int dec=ck(m,R);

     if(dec)
       break;
  }
  return n-1;
}

void explore(int i, int m[N], int R[N][N])
{
  for(int j=0;j<N;j++)
    if(R[i][j])
      if(!m[j])
      {
        m[j]=1;
        explore(j,m,R);
      }
}

int ck(int m[N], int R[N][N])
{
  for(int i=0;i<N;i++)
    if(m[i])
      explore(i,m,R);

  for(int i=0;i<N;i++)
    if(!m[i])
      return 0;
  return 1;
}

void prR(int m[N], Edge E[L], int e[L], char* fname)
{
  const int R = 70;
  const int r = 7;
  const int margin = 1;

  double x[N];
  double y[N];

  for(int i=0;i<N;i++)
  {
    double a=2*M_PI*i/N;
    x[i]=R*cos(a);
    y[i]=R*sin(a);
  }

  FILE *f;
  f=fopen(fname,"w");
  fprintf(f,"%%!PS-Adobe-3.0 EPSF-3.0\r\n");
  fprintf(f,"%%%%BoundingBox: 0 0 %d %d\r\n",2*(R+r+margin),2*(R+r+margin));
  fprintf(f,"%d %d translate\r\n",R+r+margin,R+r+margin);

  for(int l=0;l<L;l++)
  {
    int i=E[l].i;
    int j=E[l].j;
    if(e[l])
      fprintf(f,"3 setlinewidth 0 setgray\r\n");
    else
      fprintf(f,"1 setlinewidth 0.5 setgray\r\n");
    fprintf(f,"%f %f moveto %f %f lineto stroke\r\n",x[i],y[i],x[j],y[j]);
  }

  for(int i=0;i<N;i++)
  {
    fprintf(f,"1 setgray\r\n");
    fprintf(f,"newpath %f %f %d 0 360 arc fill\r\n",x[i],y[i],r);   
    if(m[i])
      fprintf(f,"3 setlinewidth 0 setgray\r\n");
    else
      fprintf(f,"1 setlinewidth 0.5 setgray\r\n");
    fprintf(f,"newpath %f %f %d 0 360 arc stroke\r\n",x[i],y[i],r);   
  }
  fprintf(f,"showpage\r\n");   
  fclose(f);
}

void init(int p[N+L], int m[N], int e[L], int R[N][N])
{
  for(int i=0;i<N;i++)
    for(int j=0;j<N;j++)
      R[i][j]=0;
  for(int i=0;i<N+L;i++)
    p[i]=0;
  for(int i=0;i<N;i++)
    m[i] = 0;
  for(int i=0;i<L;i++)
    e[i] = 0;
}

int set(int m[N], int e[L], int pa, Edge E[L], int R[N][N])
{
  int res=1;
  if(pa<N)
  {
    if(m[pa])
      res=0;
    else
      m[pa]=1;
  }
  else
  {
    e[pa-N]=1;
    int i=E[pa-N].i;
    int j=E[pa-N].j;
    R[i][j]=1;
    R[j][i]=1;
  }
  return res;
}

int newp(int p[N+L], int n)
{
  int i = rand()%(N+L-n+1);
  int j = -1;
  for(;i>=0;i--)
  {
    j++;
    for(;p[j];j++)
      ;
  }
  p[j]=1;
  return j;
}

void initE(Edge E[L])
{
  const int K = N*(N-1)/2;
  int g[K];

  for(int i=0;i<K;i++)
    g[i] = 0;

  for(int n=1; n<=L; n++) //adding 18 random edges
  {
    int i = rand()%(K-n+1);
    int j = -1;
    for(;i>=0;i--)
    {
      j++;
      for(;g[j];j++)
        ;
    }
    g[j]=1;
  }

  int l=0;
  for(int k=0;k<K;k++)
    if(g[k])
    {
      int n=k;
      for(int i=0;i<N-1 && n>-1;i++)
        for(int j=i+1;j<N && n>-1;j++,n--)
        {
          if(n==0)
          {
            E[l].i=i;
            E[l].j=j;
            l++;
          }
        }
    }
}

void prE(Edge E[L], char* fname)
{
  const double R = 70;
  const double r = 7;
  const double margin = 1;

  double x[N];
  double y[N];

  for(int i=0;i<N;i++)
  {
    double a=2*M_PI*i/N;
    x[i]=R*cos(a);
    y[i]=R*sin(a);
  }

  FILE *f;
  f=fopen(fname,"w");
  fprintf(f,"%%!PS-Adobe-3.0 EPSF-3.0\r\n");
  fprintf(f,"%%%%BoundingBox: 0 0 %f %f\r\n",2*(R+r+margin),2*(R+r+margin));
  fprintf(f,"%f %f translate\r\n",R+r+margin,R+r+margin);

  for(int l=0;l<L;l++)
  {
    int i=E[l].i;
    int j=E[l].j;
    fprintf(f,"%f %f moveto %f %f lineto stroke\r\n",x[i],y[i],x[j],y[j]);
  }

  fprintf(f,"1 setgray\r\n");
  for(int i=0;i<N;i++)
    fprintf(f,"newpath %f %f %f 0 360 arc fill\r\n",x[i],y[i],r);   

  fprintf(f,"0 setgray\r\n");
  for(int i=0;i<N;i++)
    fprintf(f,"newpath %f %f %f 0 360 arc stroke\r\n",x[i],y[i],r);   

  fprintf(f,"showpage\r\n");   
  fclose(f);
}