#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>
#include <iostream>
using namespace std;

struct Edge
{
  int i;
  int j;
};

struct Side
{
  int i;
  int j;
  int k;
};

class Graph
{
protected:

  int Nn;
  int En;
  int Sn;

  Edge *E;
  Side *S;

  void randomE();
  void randomS();

public:

  Graph(int Nn_, int En_, int Sn_);
  ~Graph();
  void prG(char* fname);
  void saveG();
  double degree_variance();
  void kring_half(int s1, int s2);
  void kring_sides(int s1, int s2, int s3);
  void torus_sides();
  void prGtorus();
};

class Receive: public Graph
{
  int *nodes;
  int *edges;
  int *sides;
  int *packets;
  int pn;
  int newp();
  FILE * ps;
  FILE * eps;
  char fout[256];
  static const double di_R = 200;
  static const double di_r = 5;
  static const double di_margin = 1;
  int exploreE();
  int exploreS();
  void fprR(FILE *f);
protected:
  int Pn;
public:
  Receive(int Nn_, int En_, int Sn_);
  ~Receive();
  void init();
  int rx();
  void propen(char *fname);
  void prR();
  void prclose();
  int decode();
};

class Rate: public Receive
{
  int excess;
  int *succ;
  int passes;
  double perf;
  void compperf();
public:
  Rate(int Nn_, int En_, int Sn_);
  ~Rate();
  void rate_and_draw(int T, char* fname);
  void rate(int T);
  void prsucc(char *fname);
  double getperf();
  void saveChart();
};

void randseed();
void subset(int *S, int N, int n);

main()
{
  randseed();
  Rate* R;

  const int Nn=18;
  const int Sn=2*Nn;
  const int En=0;

  R=new Rate(Nn,En,Sn);
  R->torus_sides();
  R->rate(10000);
  cout << "perf = " << R->getperf() << endl;
  cout << "degree variance = " << R->degree_variance() << endl;
  R->saveG();
  R->prGtorus();
  R->saveChart();
  delete R;

}

void Graph::prGtorus()
{
  if( !(En==0 && Sn==2*Nn) )
    return;
  int n=(int)sqrt((double)Nn/2);
  if( 2*n*n != Nn )
    return;

  const double R = 20;
  const double r = 4;
  const double margin = 1;

  double x[Nn];
  double y[Nn];

  double x0=0;
  double y0=n*R;
  double x1=R*sin(M_PI/3);
  double y1=R/2+n*R;

  double dy=R;
  double dx=R*2*sin(M_PI/3);

  double he = y1;
  double wi = dx*n;

  for(int i=0;i<n;i++)
    for(int j=0;j<n;j++)
    {
       int me=i*n+j;
       x[2*me+0]=x0+j*dx;
       y[2*me+0]=y0-i*dy;

       x[2*me+1]=x1+j*dx;
       y[2*me+1]=y1-i*dy;
    }

  char fname[256];
  sprintf(fname,"torus-sides-%dx%d.eps",n,n);

  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",wi+2*r+2*margin,he+2*r+2*margin);
  fprintf(f,"%f %f translate\r\n",margin+r,margin+r);

  fprintf(f,"0.9 setgray\r\n");
  for(int l=0;l<Sn;l++)
  {
    int i=S[l].i;
    int j=S[l].j;
    int k=S[l].k;
    fprintf(f,"%f %f moveto %f %f lineto %f %f lineto closepath fill\r\n",
      x[i],y[i],x[j],y[j],x[k],y[k]);
  }

  fprintf(f,"0.8 setgray\r\n");
  for(int l=0;l<Sn;l++)
  {
    int i=S[l].i;
    int j=S[l].j;
    int k=S[l].k;
    fprintf(f,"%f %f moveto %f %f lineto %f %f lineto closepath stroke\r\n",
      x[i],y[i],x[j],y[j],x[k],y[k]);
  }

  fprintf(f,"0 setgray\r\n");
  for(int l=0;l<En;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<Nn;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<Nn;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);
}

void Graph::torus_sides()
{
  if( !(En==0 && Sn==2*Nn) )
    return;
  int n=(int)sqrt((double)Nn/2);
  if( 2*n*n != Nn )
    return;
  for(int i=0;i<n;i++)
    for(int j=0;j<n;j++)
    {
       int me=i*n+j;
       int rig=i*n+(j+1)%n;
       int dow=((i+1)%n)*n+j;
       int cor=((i+1)%n)*n+(j+1)%n;

       S[4*me+0].i=2*dow+1;
       S[4*me+0].j=2*dow+0;
       S[4*me+0].k=2*me+0;

       S[4*me+1].i=2*dow+1;
       S[4*me+1].j=2*me+0;
       S[4*me+1].k=2*me+1;

       S[4*me+2].i=2*dow+1;
       S[4*me+2].j=2*me+1;
       S[4*me+2].k=2*rig+0;

       S[4*me+3].i=2*dow+1;
       S[4*me+3].j=2*rig+0;
       S[4*me+3].k=2*cor+0;

    }
}

void Graph::kring_sides(int s1, int s2, int s3)
{
  if( !(En==0 && Sn==2*Nn) )
    return;

  for(int i=0;i<Nn;i++)
  {
    S[i].i=i;
    S[i].j=(i+1)%Nn;
    S[i].k=(i+s1)%Nn;

    S[Nn+i].i=i;
    S[Nn+i].j=(i+s2)%Nn;
    S[Nn+i].k=(i+s3)%Nn;
  }

}

void Graph::kring_half(int s1, int s2)
{
  if(! (En==Nn && Sn==Nn) )
    return;
  for(int i=0;i<Nn;i++)
  {
    E[i].i=i;
    E[i].j=(i+1)%Nn;

    S[i].i=i;
    S[i].j=(i+s1)%Nn;
    S[i].k=(i+s2)%Nn;
  }
}

double Graph::degree_variance()
{
  int edge_degree[Nn];
  int side_degree[Nn];
  for(int i=0;i<Nn;i++)
  {
    edge_degree[i]=0;
    side_degree[i]=0;
  }
  for(int l=0;l<En;l++)
  {
    edge_degree[E[l].i]++;
    edge_degree[E[l].j]++;
  }
  for(int l=0;l<Sn;l++)
  {
    side_degree[S[l].i]++;
    side_degree[S[l].j]++;
    side_degree[S[l].k]++;
  }

  double vari=0;
  for(int i=0;i<Nn;i++)
    vari += abs(edge_degree[i]/2.0 + side_degree[i]/3.0 - (double)(En+Sn)/Nn);

  return vari;
}

void Rate::saveChart()
{
  char fname[256];
  sprintf(fname,"out-chart-N%02d-E%02d-S%02d.csv",Nn,En,Sn);
  FILE *f=fopen(fname,"w");

  fprintf(f,",N=%d E=%d S=%d\r\n",Nn,En,Sn);

  for(int i=0;i<=excess;i++)
    fprintf(f,"%3d, %f\r\n",i,(double)succ[i]/passes);

  fclose(f);
}

double Rate::getperf()
{
  return perf;
}

void Graph::saveG()
{
  char fname[256];
  sprintf(fname,"out-graph-N%02d-E%02d-S%02d.eps",Nn,En,Sn);
  prG(fname);
}

void Rate::compperf()
{
  perf=0;
  for(int i=0;i<=excess;i++)
    perf+=(double)succ[i];
  perf/=(double)passes;
}

void Rate::prsucc(char *fname)
{
  FILE *f;
  f=fopen(fname,"a");

  fprintf(f,"\r\n");

  for(int i=0;i<=excess;i++)
    fprintf(f,"%3d, %f\r\n",i,(double)succ[i]/passes);

  fclose(f);
}

void Rate::rate(int T)
{
  for(int i=0;i<=excess;i++)
    succ[i]=0;

  for(int t=0;t<T;t++)
  {
    init();
    int dec=0;
    int n=0;
    do
    {
      int notbad=rx();
      n++;
      if(n<Nn)
        continue;
      if(!notbad)
        continue;
      dec=decode();
    }
    while(!dec);

    for(int i=n-Nn;i<=excess;i++)
      succ[i]++;

  }
  passes=T;

  compperf();
}

void Rate::rate_and_draw(int T, char* fname)
{
  for(int i=0;i<=excess;i++)
    succ[i]=0;

  propen(fname);

  for(int t=0;t<T;t++)
  {
    init();
    prR();
    prR();
    int dec;
    int n=0;
    do
    {
      int notbad=rx();
      n++;
      if(!notbad)
        continue;
      prR();
      dec=decode();
      prR();
    }
    while(!dec);

    for(int i=n-Nn;i<=excess;i++)
      succ[i]++;

  }
  passes=T;

  prclose();
}

Rate::Rate(int Nn_, int En_, int Sn_): Receive(Nn_,En_,Sn_)
{
  //cout << "constructing a Rate object" << endl;
  excess=Pn-Nn;
  succ=new int[excess+1];
}

Rate::~Rate()
{
  //cout << "destructing the Rate object" << endl;
  delete[] succ;
}

int Receive::decode()
{
  while(exploreE()+exploreS())
    ;

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

int Receive::exploreS()
{
  int n=0;
  for(int l=0;l<Sn;l++)
    if(sides[l])
    {
      int i=S[l].i;
      int j=S[l].j;
      int k=S[l].k;
      if(nodes[i]+nodes[j]+nodes[k]==2)
      {
        nodes[i]=nodes[j]=nodes[k]=1;
        n++;
      }
    }
  return n;
}

int Receive::exploreE()
{
  int n=0;
  for(int l=0;l<En;l++)
    if(edges[l])
    {
      int i=E[l].i;
      int j=E[l].j;
      if(nodes[i]^nodes[j])
      {
        nodes[i]=nodes[j]=1;
        n++;
      }
    }
  return n;
}

void Receive::propen(char* fname)
{
  char fname_ps[256];
  char fname_eps[256];
  sprintf(fname_ps,"%s.ps",fname);
  sprintf(fname_eps,"%s.eps",fname);

  strcpy(fout,fname);

  ps=fopen(fname_ps,"w");
  eps=fopen(fname_eps,"w");

  fprintf(ps,"%%!PS\r\n");
  fprintf(ps,"[/View [/XYZ null null 1]\r\n");
  fprintf(ps,"/Page 1\r\n");
  fprintf(ps,"/PageMode /FullScreen\r\n");
  fprintf(ps,"/DOCVIEW pdfmark\r\n");

  fprintf(eps,"%%!PS-Adobe-3.0 EPSF-3.0\r\n");
  fprintf(eps,"%%%%BoundingBox: 0 0 %f %f\r\n",
    2*(di_R+di_r+di_margin),2*(di_R+di_r+di_margin));
}

void Receive::prclose()
{
  fclose(ps);
  fclose(eps);

  char mkgif[1024];
  char mkpdf[1024];

  sprintf(mkgif,"convert -delay 100 %s.eps %s.gif",fout,fout);
  sprintf(mkpdf,
    "gswin32c -dNOPAUSE -dBATCH"
    " -dDEVICEWIDTHPOINTS=%.0f -dDEVICEHEIGHTPOINTS=%.0f"
    " -sDEVICE=pdfwrite -sOutputFile=%s.pdf -f %s.ps > mkpdf.log",
    2*(di_R+di_r+di_margin),2*(di_R+di_r+di_margin),fout,fout);

  system(mkgif);
  system(mkpdf);
}

void Receive::prR()
{
  fprR(ps);
  fprR(eps);
}

void Receive::fprR(FILE *f)
{
  const float width1=0.6;
  const float graye1=0.4;
  const float grayl1=0.9;
  const float grayf1=0.98;

  const float width2=1.8;
  const float graye2=0;
  const float grayl2=0.8;
  const float grayf2=0.92;

  double x[Nn];
  double y[Nn];

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

  fprintf(f,"\r\n%f %f translate\r\n",di_R+di_r+di_margin,di_R+di_r+di_margin);

  //drawing the sides (the redundant packets for triplets of nodes)

  fprintf(f,"%f setlinewidth\r\n",width1);
  fprintf(f,"%f setgray\r\n",grayf1);
  for(int l=0;l<Sn;l++)
    if(!sides[l])
    {
      int i=S[l].i;
      int j=S[l].j;
      int k=S[l].k;
      fprintf(f,"%f %f moveto %f %f lineto %f %f lineto closepath fill\r\n",
        x[i],y[i],x[j],y[j],x[k],y[k]);
    }

  fprintf(f,"%f setlinewidth\r\n",width2);
  fprintf(f,"%f setgray\r\n",grayf2);
  for(int l=0;l<Sn;l++)
    if(sides[l])
    {
      int i=S[l].i;
      int j=S[l].j;
      int k=S[l].k;
      fprintf(f,"%f %f moveto %f %f lineto %f %f lineto closepath fill\r\n",
        x[i],y[i],x[j],y[j],x[k],y[k]);
    }

  fprintf(f,"%f setlinewidth\r\n",width1);
  fprintf(f,"%f setgray\r\n",grayl1);
  for(int l=0;l<Sn;l++)
    if(!sides[l])
    {
      int i=S[l].i;
      int j=S[l].j;
      int k=S[l].k;
      fprintf(f,"%f %f moveto %f %f lineto %f %f lineto closepath stroke\r\n",
        x[i],y[i],x[j],y[j],x[k],y[k]);
    }

  fprintf(f,"%f setlinewidth\r\n",width2);
  fprintf(f,"%f setgray\r\n",grayl2);
  for(int l=0;l<Sn;l++)
    if(sides[l])
    {
      int i=S[l].i;
      int j=S[l].j;
      int k=S[l].k;
      fprintf(f,"%f %f moveto %f %f lineto %f %f lineto closepath stroke\r\n",
        x[i],y[i],x[j],y[j],x[k],y[k]);
    }

  //drawing the edges (the redundant packets for pairs of nodes)

  fprintf(f,"%f setlinewidth\r\n",width1);
  fprintf(f,"%f setgray\r\n",graye1);
  for(int l=0;l<En;l++)
    if(!edges[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,"%f setlinewidth\r\n",width2);
  fprintf(f,"%f setgray\r\n",graye2);
  for(int l=0;l<En;l++)
    if(edges[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]);
    }

  //drawing the nodes (the source packets)

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

  fprintf(f,"%f setlinewidth\r\n",width1);
  fprintf(f,"%f setgray\r\n",graye1);
  for(int i=0;i<Nn;i++)
    if(!nodes[i])
      fprintf(f,"newpath %f %f %f 0 360 arc stroke\r\n",x[i],y[i],di_r);   

  fprintf(f,"%f setlinewidth\r\n",width2);
  fprintf(f,"%f setgray\r\n",graye2);
  for(int i=0;i<Nn;i++)
    if(nodes[i])
      fprintf(f,"newpath %f %f %f 0 360 arc stroke\r\n",x[i],y[i],di_r);   

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

int Receive::rx()
{
  int p=newp();
  if(p==-1)
    return 0;
  if(p<Nn)
  {
    if(nodes[p])
      return 0;
    nodes[p]=1;
  }
  else if(p<Nn+En)
  {
    edges[p-Nn]=1;
  }
  else if(p<Nn+En+Sn)
  {
    sides[p-Nn-En]=1;
  }
  return 1;
}

int Receive::newp()
{
  if(pn>=Pn)
    return -1;

  int i = rand()%(Pn-pn);
  int j = -1;
  for(;i>=0;i--)
  {
    j++;
    for(;packets[j];j++)
      ;
  }
  packets[j]=1;
  pn++;

  return j;
}

void Receive::init()
{
  for(int i=0;i<Nn;i++)
    nodes[i]=0;
  for(int i=0;i<En;i++)
    edges[i]=0;
  for(int i=0;i<Sn;i++)
    sides[i]=0;
  for(int i=0;i<Pn;i++)
    packets[i]=0;
  pn=0;
}

Receive::~Receive()
{
  //cout << "destructing the Receive object" << endl;
  delete[] nodes;
  delete[] edges;
  delete[] sides;
  delete[] packets;
}

Receive::Receive(int Nn_, int En_, int Sn_): Graph(Nn_,En_,Sn_)
{
  //cout << "constructing a Receive object" << endl;
  nodes=new int[Nn];
  edges=new int[En];
  sides=new int[Sn];
  Pn=Nn+En+Sn;
  packets=new int[Pn];
  init();
}

void Graph::randomS()
{
  const int K = Nn*(Nn-1)*(Nn-2)/6;
  int Set[K];
  subset(Set,K,Sn);

  int l=0;
  for(int m=0;m<K;m++)
    if(Set[m])
    {
      int n=m;
      for(int i=0;i<Nn-2 && n>-1;i++)
        for(int j=i+1;j<Nn-1 && n>-1;j++)
          for(int k=j+1;k<Nn && n>-1;k++,n--)
            if(n==0)
            {
              S[l].i=i;
              S[l].j=j;
              S[l].k=k;
              l++;
            }
    }
}

void Graph::randomE()
{
  const int K = Nn*(Nn-1)/2;
  int Set[K];
  subset(Set,K,En);

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

Graph::Graph(int Nn_, int En_, int Sn_)
{
  Nn=Nn_;
  En=En_;
  Sn=Sn_;

  /*
  cout << "constructing a Graph with " <<
     Nn << " nodes, " <<
     En << " edges and " <<
     Sn << " sides" << endl;
  */

  E=new Edge[En];
  S=new Side[Sn];
  randomE();
  randomS();
}

Graph::~Graph()
{
  //cout << "destructing the Graph" << endl;
  delete[] E;
  delete[] S;
}

void subset(int *Set, int N, int n)
{
  for(int i=0;i<N;i++)
    Set[i] = 0;

  for(int m=0; m<n; m++)
  {
    int i = rand()%(N-m);
    int j = -1;
    for(;i>=0;i--)
    {
      j++;
      for(;Set[j];j++)
        ;
    }
    Set[j]=1;
  }
}

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

void Graph::prG(char* fname)
{
  const double R = 70;
  const double r = 4;
  const double margin = 1;

  double x[Nn];
  double y[Nn];

  for(int i=0;i<Nn;i++)
  {
    double a=2*M_PI*i/Nn;
    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);

  fprintf(f,"0.9 setgray\r\n");
  for(int l=0;l<Sn;l++)
  {
    int i=S[l].i;
    int j=S[l].j;
    int k=S[l].k;
    fprintf(f,"%f %f moveto %f %f lineto %f %f lineto closepath fill\r\n",
      x[i],y[i],x[j],y[j],x[k],y[k]);
  }

  fprintf(f,"0.8 setgray\r\n");
  for(int l=0;l<Sn;l++)
  {
    int i=S[l].i;
    int j=S[l].j;
    int k=S[l].k;
    fprintf(f,"%f %f moveto %f %f lineto %f %f lineto closepath stroke\r\n",
      x[i],y[i],x[j],y[j],x[k],y[k]);
  }

  fprintf(f,"0 setgray\r\n");
  for(int l=0;l<En;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<Nn;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<Nn;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);
}