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

/* Play around with those constants */
#define U0 1.0
#define R0 1.0
#define Q0 1.0
#define QMAX 100.0
#define NBINS 1000

/* Do not change those important constants */
#define STRING_BUFFER 1024
#define TRUE 1
#define FALSE 0
#define EPSILON 1e-7
#define LOWEST_ORDER 0
#define HIGHEST_ORDER 2
#define PI 3.14159265358979323846

/* Inline helper functions */
#define MAX(x, y) ((x) < (y) ? (y) : (x))
#define MIN(x, y) ((x) < (y) ? (x) : (y))
#define ABS(x) (SIGN(x) * x)
#define SIGN(x) ((x) > 0.0 ? 1.0 : ((x) < 0.0 ? -1.0 : 0.0))

/* Variables */
double answer;
double Dsigma = 2.0;
int order;
int Zsigma;
double *H;
double s;
double qmax;
double dq;
double q;
double qp;
double chi;
int rejects;

/* Generates a random number */
double Random()
{
	double r = (rand() / (double)RAND_MAX);
	return r;
}

/* Value of U(q) */
double U(double q)
{
	//Use limit of analytic expression for small q s
	if(q < EPSILON)
	{
		return 2.0 * U0 * (R0 * R0 * R0) / 3.0;
	}

	return 2.0 * U0 * (sin(q * R0) - q * R0 * cos(q * R0)) / (q * q * q);
}

/* Value of f(q) */
double f(double q)
{
	/* If q is out of bounds or not enough statistics - just return the first term... f(q) = -u(q) */
	if(q > qmax || Zsigma < 100000)
	{
		return -U(q);
	}

	int index = (int)(q / dq);
	return H[index] * Dsigma / (double)Zsigma / dq;
}

/* Value of q' distribution function at q' */
double WQp(double q)
{
	double Q = (Q0 + q) * (Q0 + q);
	return Q0 / Q;
}

/* Value of q distribution function at q */
double WQ(double q)
{
	double R = (Q0 + qmax) / qmax;
	return R * WQp(q);
}

/* Value of chi distribution function at chi */
double WChi(double chi)
{
	/* This is just a constant - I = 2 => I^(-1) = 1 /2 */
	return 0.5;
}

/* Generate variables at order */
void Variables(int order)
{
	double r = Random();

	switch(order)
	{
		case 1:
			q   = Q0 * r / (Q0 / qmax + 1.0 - r);
			break;
		case 2:
			qp  = Q0 * r / (1.0 - r);
			r   = Random();
			chi = 2.0 * r - 1.0;
			break;
	}
}

/* Compute the value of the diagram at order */
double Level(int order)
{
	/* Value of the diagram at order */
	double value = 0.0;
	/* Normlization factor (due to detailled balance) at given order */
	double norm  = 1.0;

	switch(order)
	{
		case 0:
			value = Dsigma;
			break;
		case 1:
			norm  = 0.5;
			value = -U(q) / WQ(q);
			break;
		case 2:
			value = sqrt(q * q + qp * qp - 2.0 * q * qp * chi);
			value = -U(value) * f(qp) / (PI * WQp(qp) * WChi(chi) * WQ(q));
			break;
	}

	return norm * value;
}

/* Try a Metropolis update from srcOrder (source) to destOrder (destination) */
void Update(int srcOrder, int destOrder)
{
	double news    = 0.0;
	double r       = Random();
	double ratio   = 1.0;
	int high_order = MAX(srcOrder, destOrder);
	int low_order  = MIN(srcOrder, destOrder);

	/* If we increase the order we have to generate new variables / values */
	if(destOrder > srcOrder)
	{
		Variables(destOrder);
	}

	/* Compute the diagram values */
	double higher  = Level(high_order);
	double lower   = Level(low_order);

	/* Compute the ratio and new sign value */
	if(destOrder < srcOrder)
	{
		ratio = ABS(lower / higher);
		news  = SIGN(lower);
	}
	else
	{
		ratio = ABS(higher / lower);
		news  = SIGN(higher);
	}

	/* Reject the step ?! */
	if(ratio < 1.0 && r > ratio)
	{
		rejects++;
		return;
	}

	s     = news;
	order = destOrder;
}

/* Perform Metropolis update step */
void Step()
{
	if(LOWEST_ORDER == order)
	{
		Update(order, order + 1);
	}
	else if(HIGHEST_ORDER == order)
	{
		Update(order, order - 1);
	}
	else
	{
		double r    = Random();
		int dest    = r > 0.5 ? order + 1 : order - 1;
		Update(order, dest);
	}
}

/* Measure histogram */
void Measure()
{
	/* Non-physical diagram won't be measured */
	if(order == 0)
	{
		Zsigma++;
	}
	else if(q < qmax)
	{
		/* If we are inbound we add the value to the histogram */
		int index = (int)(q / dq);
		H[index]  = H[index] + s;
	}
}

/* Print statistics */
void Print()
{
	printf("Rejected diagram changes  =  %d\n", rejects);
	printf("Normalization statistics  =  %d\n", Zsigma);
	printf("MC -f(0) value            =  %lf\n", -f(0.0));
	printf("Exact answer              =  %lf\n", answer);
	printf("-------------------------------------\n");
}

/* Extract statistics into a file */
void Extract(int nbins)
{	
	/* Write everything to the statistics.dat */
	FILE* file = fopen("statistic.dat", "w");
	int i = 0;

	for(i = 0; i < nbins; i++)
	{
		double x = dq * ((double)i + 0.5);
		double y = f(x);
		char* line = malloc(sizeof(char) * STRING_BUFFER);
		sprintf(line, "%lf\t%lf\n", x, y);
		fputs(line, file);
	}

	fclose(file);
}

/* Metropolis iteration */
void Metropolis(int nbins, double _qmax)
{
	/* Initialize some values */
	int i;
	int step  = 0;
	s         = 1.0;
	rejects   = 0;
	qmax      = _qmax;
	dq        = _qmax / (double)nbins;
	H         = malloc(sizeof(double) * nbins);
	Zsigma    = 0;
	order     = 0;

	/* Start up histogram with 0-values */
	for(i = 0; i < nbins; i++)
	{
		H[i] = 0.0;
	}

	/* Infinite loop of Metropolis steps */
	while(TRUE)
	{
		step++;
		Measure();
		Step();

		if(step % 10000000 == 0)
		{
			Print();
			Extract(nbins);
		}
	}
}

/* Main entry point */
int main()
{
	/* Number of bins */
	int nbins   = NBINS;
	/* Maximum of q */
	double qmax = QMAX;
	/* Initialize random generator */
	uint iseed  = (uint)time(NULL);
	srand (iseed);
	/* Calculate analytic answer (for comparison) */
	answer = 1.0 - tanh(sqrt(2.0 * U0)) / sqrt(2.0 * U0);
	/* Start Metropolis algorithm / iteration */
	Metropolis(nbins, qmax);
	return 0;
}