#include <stdio.h>
#include <math.h>
#include <stdlib.h>
double drand48();

#define NN 1000

int v[(2*NN+1)*(2*NN+1)];
int le_x[(2*NN+1)*(2*NN+1)];
int le_y[(2*NN+1)*(2*NN+1)];
int lerw[(2*NN+1)*(2*NN+1)];
int ct[4];

void print_preamble(FILE *file, int time_to_exit, int loop_erase_time, int seed);
int run_walk(int seed);

int main(int argc, char *argv[])
{
	FILE *file;
	int i,j;
	char file_name[256];
	int seed;
	int time;
	
	for(seed=100;seed<200;seed++)
	{
		time = run_walk(seed);
		if(ct[0]+ct[1]+ct[2]+ct[3]>5000000)
		{
			printf("printing file %d (num_steps = %d)\n",seed,ct[0]+ct[1]+ct[2]+ct[3]);
			sprintf(file_name,"./results/lerw-%d.eps",seed);
			file = fopen(file_name,"w");
			print_preamble(file,ct[0]+ct[1]+ct[2]+ct[3],time,seed);
			//fprintf(file,"0 0 m\n",50.);
			
			for(i=0;i<2*NN+1;i++)
			{
				for(j=0;j<2*NN+1;j++)
				{
					if(v[(2*NN+1)*i+j]==1)
						fprintf(file,"%d %d b\n",i,j);
				}
			}
			
			fprintf(file,"S\n");
			fprintf(file,"%d %d m\n",NN,NN);
			for(i=0;i<time;i++)
			{
				fprintf(file,"%d %d l\n",le_x[i]+NN,le_y[i]+NN);
				//printf("%d %d\n",le_x[i],le_y[i]);
			}
			
			fprintf(file,"2 SLW\n0 SG\nS\nshowpage\n%%EOF");
			fclose(file);
		}
	}
	return 1;
}

void print_preamble(FILE *file, int time_to_exit, int loop_erase_time, int seed)
{
	fprintf(file, "%%!PS-Adobe EPSF-3.0\n");
	fprintf(file, "%%%%BoundingBox: %d %d %d %d\n",-5,-5,2*NN+5,2*NN+5);
	fprintf(file, "%%%%Parameters: seed = %d, box_side = %d, time_to_exit = %d, LERW_length = %d\n\n",seed, 2*NN,time_to_exit,loop_erase_time);
	fprintf(file, "/m {moveto} def\n/l {lineto} def\n/rm {rmoveto} def\n/rl {rlineto} def\n/SG {setgray} def\n");
	fprintf(file, "/SLW {setlinewidth} def\n/S {stroke} def\n/F {fill} def\n/NP {newpath} def\n/CP {closepath} def\n");
	//fprintf(file, "\n/b {m\ncurrentpoint %g 0 360 arc\nF\n}def\n",0.5);	
	fprintf(file, "\n/b {m\n-0.5 -0.5 rm\n1 0 rl\n0 1 rl\n-1 0 rl\nCP\n%g SG\nF }def\n",0.6);
	//fprintf(file, "\n/bb {m\n-0.5 -0.5 rm\n1 0 rl\n0 1 rl\n-1 0 rl\nCP\n0 SG\nF }def\n",0.0);
	fprintf(file, "\n0 0 m\n%d 0 l\n%d %d l\n0 %d l\nCP\n2 SLW\nS\n\n",2*NN,2*NN,2*NN,2*NN);
}

int run_walk(int seed)
{
	int x_sum,y_sum;
	int i,j,time;
	
	srand48(seed);
	x_sum=0; 
	y_sum=0;
	time=0;
	
	// setting all arrays to zero
	for(i=0;i<4;i++)
		ct[i]=0;
	for(i=0;i<2*NN+1;i++)
	{
		for(j=0;j<2*NN+1;j++)
		{
			// storage for positions ever visited
			v[(2*NN+1)*i+j]=0;
			// storage for loop erasure
			lerw[(2*NN+1)*i+j]=0;
			// auxiliary arrays
			le_x[(2*NN+1)*i+j]=0;
			le_y[(2*NN+1)*i+j]=0;
		}
	}
	// setting origin as visited
	v[(2*NN+1)*NN+NN]=1;
	lerw[(2*NN+1)*NN+NN]=1;
	le_x[0]=0;
	le_y[0]=0;
	
	// running walk & loop-erasing
	while(x_sum>-NN && x_sum < NN && y_sum>-NN && y_sum<NN)
	{
		double s = drand48();
		int i,f;
		// generate a step
		f = (int)4*s;
		if(f==0)
			++x_sum;
		if(f==1)
			--x_sum;
		if(f==2)
			++y_sum;
		if(f==3)
			--y_sum;
		// check if visited and loop-erase otherwise
		if(lerw[(2*NN+1)*(x_sum+NN)+y_sum+NN]==1)
		{
			// if true then run backwards
			while(time >=0 && ((le_x[time] != x_sum) || (le_y[time] != y_sum)))
			{
				//printf("erasing %d %d (time = %d)\n",le_x[time],le_y[time],time);
				// erase
				lerw[(2*NN+1)*(le_x[time]+NN)+le_y[time]+NN] = 0;
				// and decrease time
				--time;
			}
		}
		// if not visited then
		else
		{
			// increment time
			++time;
			// and add point to loop-erasure and aux-arrays
			lerw[(2*NN+1)*(x_sum+NN)+y_sum+NN]=1;
			le_x[time] = x_sum;
			le_y[time] = y_sum;
		}
		// mark point as visited by SRW
		v[(2*NN+1)*(x_sum+NN)+y_sum+NN]=1;
		// record the direction of the step
		++ct[f];
		//printf("s=%.2g\tf = %d,\t x_sum = %d,\t y_sum = %d\n",s,f,x_sum,y_sum);
	}
	//printf("time = %d\n",time);
	return time;
}