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

double tvalue(double pred, double obs, double var, int N)
{
	/*t-value*/
	return (pred - obs)*sqrt((double)N)/sqrt(var);;
}

double pvalue(double t)
{
	double x;
	if (t < .1)
		x = 0.5;
	else if (t < .2)
		x = 0.46;
	else if (t < .3)
		x = 0.421;
	else if (t < .4)
		x = .382;
	else if (t < .5)
		x = .345;
	else if (t < .6)
		x = .309;
	else if (t < .7)
		x = .274;
	else if (t < .8)
		x = .242;
	else if (t < .9)
		x = .212;
	else if (t < 1.0)
		x = .184;
	else if (t < 1.1)
		x = .159;
	else if (t < 1.2)
		x = .136;
	else if (t < 1.3)
		x = .115;
	else if (t < 1.4)
		x = .097;
	else if (t < 1.5)
		x = .081;
	else if (t < 1.6)
		x = .067;
	else if (t < 1.7)
		x = .055;
	else if (t < 1.8)
		x = .045; 
	else if (t < 1.9)
		x = .036;
	else if (t < 2.0)
		x = .029;
	else if (t < 2.1)
		x = .023;
	else if (t < 2.2)
		x = .018;
	else if (t < 2.3)
		x = .014;
	else if (t < 2.4)
		x = .011;
	else if (t < 2.5)
		x = .008;
	else if (t < 2.6)
		x = .006;
	else if (t < 2.7)
		x = .005;
	else if (t < 2.9)
		x = .003;
	else if (t < 3.0)
		x = .002;
	else if (t < 3.3)
		x = .001;
	else
		x = 0.0;

	/* 2-tailed */
	return x*2;
}


int main()
{


  // int names[] = {110017, 120041, 130021, 140053, 150001, 160009, 170101, 180001, 190093, 200041};
  // int trials[] = {18240, 24000, 15680, 21200, 20000, 25984, 19440, 24000, 25920, 26656};

  // int names[] = {99923, 99961, 99971, 99989, 99991, 100003, 100019, 100049, 100069, 100103};
  // int trials[] = {48852, 10752, 36864, 21420, 24000, 28560, 48804, 22464, 16080, 50050};

	int names[] = {110459, 120167, 130127, 140123, 150083, 160079, 170099, 180023, 190523, 200087};
	int trials[] = {55228, 60082, 65062, 70060, 75040, 80038, 85048, 90010, 95260, 100042};

  //	int names[] = {100043, 100057, 106261};
  //	int trials[] = {50020, 15120, 10560};
	int size = 10;

	int i = 0;
	for (i=0; i < size; i++)
	{
		printf("%d\n", i);

		char fileStr[20];
		char fileStr2[20];
		sprintf(fileStr, "results_%d.dat", names[i]);
		sprintf(fileStr2, "processed_%d.dat", names[i]);

		FILE * in;
		in = fopen(fileStr, "r");

		FILE * out;
		out = fopen(fileStr2, "w");

		/* GAMMA(k + 1) / GAMMA(k + 0.5) */
		double r1over5;
		/* GAMMA(k + 1.5) / GAMMA(k) */
		double r15over0;
		/* 2k = n */
		double k = (names[i] - 1)/2;

		
		/* You have to input the values above at the top of the .dat to get this to work.
		   You'll also have to fill in the components variance and components mean (no way around it) */
		fscanf(in, "%lf %lf\n", &r1over5, &r15over0);

		double mean = 0.0;
		double var = 0.0;
		double var2 = 0.0;
		double t = 0.0;
		double p = 0.0;
		double vart = 0.0;
		double varp = 0.0;
		double predMean;
		double predVar;

		fprintf(out, "\t\t\tObserved\tPredicted\tt-value\t\tp-value\n");

		/* ======== COMPONENTS ============*/
		fscanf(in, "components: %lf %lf \n", &mean, &predMean);
		fscanf(in, "components variance: %lf %lf \n", &var, &predVar);
		fscanf(in, "components variance (var): %lf\n", &var2);

		t = tvalue(predMean, mean, var, trials[i]);
		p = pvalue(fabs(t));

		vart = tvalue(predVar, var, var2, trials[i]);
		varp = pvalue(fabs(vart));


		fprintf(out, "components: \t\t%lf\t%lf\t%lf\t%lf \n\n", mean, predMean, t, p );
		fprintf(out, "components variance: \t%lf\t%lf\t%lf\t%lf \n\n", var, predVar, vart, varp);
		fprintf(out, "components var var: \t%lf\n\n", var2);



		/* ======== CYCLIC ============*/
		fscanf(in, "cyclic nodes: %lf \n", &mean);
		fscanf(in, "cyclic nodes variance: %lf\n", &var);
		fscanf(in, "cyclic nodes variance (var): %lf\n", &var2);

		/* -(1.*(GAMMA(k+.5)-1.772453851*GAMMA(k+1.)))/GAMMA(k+.5) */
		predMean = -1*(1 - 1.772453851*r1over5);
		
		/* -(1.*(-3.*GAMMA(k+.5)-4.*GAMMA(k+.5)*k+5.317361553*GAMMA(k+1.)))/GAMMA(k+.5)-
		1.*(-1.*GAMMA(k+.5)+1.772453851*GAMMA(k+1.))^2/GAMMA(k+.5)^2 */
		predVar = -1*(-3 - 4*k + 5.317361553*r1over5) - (-1 + 1.772453851*r1over5)*(-1 + 1.772453851*r1over5);

		t = tvalue(predMean, mean, var, trials[i]);	
		p = pvalue(fabs(t));

		vart = tvalue(predVar, var, var2, trials[i]);
		varp = pvalue(fabs(vart));

		fprintf(out, "cyclic nodes: \t\t%lf\t%lf\t%lf\t%lf\n\n", mean, predMean, t, p);
		fprintf(out, "cyclic nodes variance: \t%lf\t%lf\t%lf\t%lf\n\n", var, predVar, vart, varp);
		fprintf(out, "cyclic nodes var var: \t%lf\n\n", var2);


		/* ======== AVG CYCLE ============*/
		fscanf(in, "avg cycle: %lf\n", &mean);
		fscanf(in, "avg cycle variance: %lf\n", &var);
		fscanf(in, "avg cycle variance (var): %lf\n", &var2);


		/* .8862269255*GAMMA(k+1.)/GAMMA(k+.5) */
		predMean =  .8862269255*r1over5;
				
		/*-(.1666666667*(-8.*GAMMA(k+1.5)+5.317361553*k*GAMMA(k)))*GAMMA(k+1.)/(GAMMA(k)*k*GAMMA(k+.5))
		  -.7853981635*GAMMA(k+1.)^2/GAMMA(k+.5)^2 */
		predVar = -r1over5*(-8*r15over0 + 5.317361553*k)/(k*6) - .7853981635*(r1over5*r1over5);
		
		t = tvalue(predMean, mean, var, trials[i]);
		p = pvalue(fabs(t));

		vart = tvalue(predVar, var, var2, trials[i]);
		varp = pvalue(fabs(vart));

		fprintf(out, "avg cycle: \t\t%lf\t%lf\t%lf\t%lf\n\n", mean, predMean, t, p);
		fprintf(out, "avg cycle variance: \t%lf\t%lf\t%lf\t%lf\n\n", var, predVar, vart, varp);
		fprintf(out, "avg cycle var var: \t%lf\n\n", var2);


		/* ======== AVG TAIL ============*/
		fscanf(in, "avg tail: %lf\n", &mean);
		fscanf(in, "avg tail variance: %lf\n", &var);
		fscanf(in, "avg tail variance (var): %lf\n", &var2);

		/*(.5*(1.772453851*GAMMA(k+1.)+1.772453851*k*GAMMA(k+1.)-2.*GAMMA(k+1.5)))/(GAMMA(k+.5)*k)*/
		predMean = .5*(1.772453851*r1over5 + 1.772453851*k*r1over5 - 2 * (k + .5))/k;
		
		/*((1/6)*(18.*GAMMA(k+1.5)+8.*GAMMA(k+1.5)*k-15.95208466*k*GAMMA(k+1.)-15.95208466*GAMMA(k+1.)))/(k*GAMMA(k+.5))
		-.25*(1.772453851*GAMMA(k+1.)+1.772453851*k*GAMMA(k+1.)-2.*GAMMA(k+1.5))^2/(k^2*GAMMA(k+.5)^2)*/
		predVar = (18*(k + .5) + 8*k*(k + .5) - 15.95208466*k*r1over5 - 15.95208466*r1over5)/(k*6) 
			- .25*(1.772453851*r1over5 + 1.772453851*k*r1over5 - 2*(k+0.5))* (1.772453851*r1over5 + 1.772453851*k*r1over5 - 2*(k+0.5)) / (k * k);

		t = tvalue(predMean, mean, var, trials[i]);
		p = pvalue(fabs(t));

		vart = tvalue(predVar, var, var2, trials[i]);
		varp = pvalue(fabs(vart));

		fprintf(out, "avg tail variance: \t%lf\t%lf\t%lf\t%lf\n\n", mean, predMean, t, p);	
		fprintf(out, "avg tail: \t\t%lf\t%lf\t%lf\t%lf\n\n", var, predVar, vart, varp);
		fprintf(out, "avg tail var var: \t%lf\n\n", var2);

		/* ======== IMAGE NODES ============*/
		fscanf(in, "image nodes: %lf\n", &mean);
		fprintf(out, "image nodes: \t\t%lf\t\%lf\n\n", mean, k);


		/* ======== MAX CYCLE ============*/
		fscanf(in, "max cycle: %lf\n", &mean);
		fscanf(in, "max cycle variance: %lf\n", &var);
		fscanf(in, "max cycle variance (var): %lf\n", &var2);
		
		/* 1.106596092*GAMMA(k+1.)/GAMMA(k+.5) */
		predMean = 1.106596092*r1over5;

		t = tvalue(predMean, mean, var, trials[i]);
		p = pvalue(fabs(t));

		fprintf(out, "max cycle: \t\t%lf\t%lf\t%lf\t%lf\n\n", mean, predMean, t, p);
		fprintf(out, "max cycle variance: \t%lf\n\n", var);
		fprintf(out, "max cycle var var: \t%lf\n\n", var2);

		
		/* ======== MAX Tail ============*/
		fscanf(in, "max tail: %lf\n", &mean);
		fscanf(in, "max tail variance: %lf\n", &var);
		fscanf(in, "max tail variance (var): %lf\n", &var2);
		
		/* NOTE : ASSYMPTOTIC */
		/* was: 1.737462322*sqrt(n)-1.613705639 */
		/* predMean = 1.737462322*sqrt(names[i] - 1)-1.613705639; */

		/* new: 1.737462321*sqrt(n)-2.248887061-.5000000000*ln(n) */
		predMean = 1.737462321*sqrt(names[i] - 1)-2.248887061-.5000000000*log(names[i] - 1);

		t = tvalue(predMean, mean, var, trials[i]);
		p = pvalue(fabs(t));

		fprintf(out, "max tail: \t\t%lf\t%lf\t%lf\t%lf\n\n", mean, predMean, t, p);
		fprintf(out, "max tail variance: \t%lf\n\n", var);
		fprintf(out, "max tail var var: \t%lf\n\n", var2);

		fclose(in);
		fclose(out);

	}

}