#include <iostream>
#include <cmath>
#include <fstream>
#include <ctime>
#include "mersenne.h"

using namespace std;

#define ARRAYELEMCOUNT(x) sizeof(x) / sizeof(x[0])
#define TAB "\t"

inline double Pow2(const double a) {return a*a;}

#define Rand01() mersenne()
#define SeedRand(x) seed_mersenne(x)

struct Point
//==========
{
	double x,y;
};

inline double H(const Point& point)
//---------------------------------
{
	return point.x*point.x + point.y*point.y + 5*(point.x-point.y)*(point.x-point.y);
}

inline Point GenerateTrialPoint(const Point& point, const double C)
//--------------------------------------------
{
	Point a = {point.x + C*(Rand01()-0.5), point.y + C*(Rand01()-0.5)};
	return a;
}

typedef unsigned int uint;

int main()
//--------
{
	const double Cs[] = {0.01, 0.1, 1, 10, 100};
	//Trial constant values to be tested.

	const uint Ns[] = {0, 1000, 10000, 50000, 100000, 1000000};
	//Values at which print values of <a^2>, <b^2> expect for first element, which should be zero.
	
	cout.precision(12);
	const int seed = time(0);
	SeedRand(seed);
	ofstream datafile("data.txt", ios::trunc);
	ostream& datastrm = datafile;

	datastrm << "#Seed: " << seed << endl;

	datastrm << "#N\t<(x - y)^2>\tError\tn1se\t<(x + y)^2>\tError\tn1se" << endl;
	//n1se <-> Naive 1 Sigma Error

	for(uint c = 0; c<ARRAYELEMCOUNT(Cs); c++)
	{
		Point point = {0,0};	//Set initial point.
		double Hprev = H(point);
		uint accepted = 0;
		uint counter = 0;
		double asum = 0, a2sum = 0, bsum = 0, b2sum = 0;
		datastrm << "#C   : " << Cs[c] << endl;
		for(uint j = 1; j<ARRAYELEMCOUNT(Ns); j++)
		{
			for(uint i = Ns[j-1]; i<Ns[j]; i++)
			{
				counter++;
				Point trialpoint = GenerateTrialPoint(point, Cs[c]);
				const double deltaH = H(trialpoint) - Hprev;
				if(deltaH <= 0 || exp(-deltaH) > Rand01())
				{
					accepted++;
					Hprev += deltaH;
					point = trialpoint;
				}

				double temp = Pow2(point.x - point.y);
				asum += temp;
				a2sum += temp*temp;
				temp = Pow2(point.x + point.y);
				bsum += temp;
				b2sum += temp*temp;
			}
			const double a = asum / counter;
			const double b = bsum / counter;
			datastrm << counter << TAB << a << TAB << fabs((a - 1.0/11.0)*11.00) << TAB << sqrt((a2sum/counter - a*a)/(counter-1)) << TAB
				 << b << TAB << fabs(b - 1.0) << TAB << sqrt((b2sum/counter - b*b)/(counter-1)) << endl;
		}
		datastrm << "#Accept rate: " << double(accepted) / counter << endl << endl;
	}
	datafile.close();
	return 0;
}