#include<iostream>
#include<cstdlib>
#include<ctime>
#include<vector>
#include<math.h>
#include</usr/local/include/root/TCanvas.h>
#include</usr/local/include/root/TGraph.h>
//#include<TCanvas.h>
//#include<TGraph.h>

using namespace std;

int main()
{
  int NSAMPLES=10000, x=0, y=0;
  double X2=0, X4=0, Y2=0, Y4=0;
  vector<double> TIME, R2X, R2Y, SR2X, SR2Y;
  vector<double> NSTEP;
  NSTEP.push_back(6);
  NSTEP.push_back(10);
  NSTEP.push_back(25);
  NSTEP.push_back(50);
  NSTEP.push_back(75);
  NSTEP.push_back(100);
  NSTEP.push_back(1000);

  srand(time(0));
  //for(int i=0;i<NSTEP;i++) cout << rand()%10 << endl;
  
  //cout << "extern" << endl;
  //ensemble averages vs step length
  for(unsigned int k=0;k<NSTEP.size();k++)
    {
      //one ensemble average with fixed step length
      for(int i=0;i<NSAMPLES;i++)
	{
	  x=0; y=0; X2=0; X4=0; Y2=0; Y4=0;
	  //single system evolution
	  for(int j=0;j<NSTEP[k];j++)
	    {
	      int RND=0;
	      RND = rand()%4;
	      //cout << RND << endl;
	      if(RND==0)      {x+=1;} //cout << x << endl;} 
	      else if(RND==1) y+=1; //cout << x << endl;
	      else if(RND==2) x-=1; //cout << y << endl;
	      else if(RND==3) y-=1; //cout << y << endl;
	    }
	  //cout << x << endl;
	  //cout << pow(x,2) << endl;
	  X2 += pow(x,2);
	  cout << X2 << endl;
	  X4 += pow(x,4);
	  Y2 += pow(y,2);
	  Y4 += pow(y,4);
	}
      //cout << "X2 " << (X2/NSAMPLES) << endl << "pow " << pow((X2/NSAMPLES),2) << endl;
      R2X.push_back(X2/=NSAMPLES);
      X4/=NSAMPLES;
      //cout << "X4 " << X4 << endl;
      R2Y.push_back(Y2/=NSAMPLES);
      Y4/=NSAMPLES;
      SR2X.push_back((X4-pow((X2/NSAMPLES),2))/sqrt(NSAMPLES));
      SR2Y.push_back((Y4-pow((Y2/NSAMPLES),2))/sqrt(NSAMPLES));
    }
  
  TCanvas* c1 = new TCanvas("c1","Graph",200,10,700,500);
  Int_t n = NSTEP.size();
  Double_t ord[n], absc[n];
  for(Int_t i=0;i<n;i++)
    {
      ord[i]=R2X[i];
      absc[i]=NSTEP[i];
    }
  TGraph* gr = new TGraph(n,ord,absc);
  gr->Draw("ACP");
  c1->Update();
  
  for(unsigned int i=0;i<R2X.size();i++)
    cout << NSTEP[i] <<" "<< R2X[i] <<" "/*<< SR2X[i] << " "*/ << R2Y[i] <</* " " << SR2Y[i] <<*/ endl;

  return 0;
}