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

#define NN 100000
#define alpha 1.99

#define SRW 1

void print_preamble(FILE *file, double max, double min);
double rw_step(double u);

int main(int argc, char *argv[])
{
	int i;
	double s[NN];
	double max = 0, min = 0;
	double sum = 0;
	double scale;
	FILE *file;
	
	srand48(10000);
	
	for(i=0;i<NN;i++)
	{
		s[i] = drand48();
		//printf("%g\n",rw_step(s[i]));
		sum += rw_step(s[i]);
		if(sum>max)
			max = sum;
		else if(sum<min)
			min = sum;
		//printf("%d\n", sum);
	}
	scale = 1.*(max>-min ? max : -min);

	printf("max = %g, min = %g, scale = %g\n",max,min,scale);
	
	file = fopen("./srw.eps","w");
	print_preamble(file,max,min);
	fprintf(file,"0 %g m\n",50.);
	sum = 0;
	for(i=0;i<NN;i++)
	{
		sum += rw_step(s[i]);
		fprintf(file,"%g %g l\n",(i+1)*100./NN,50.+45.*sum/scale);
		
	}
	fprintf(file,"%g SLW\nS\nshowpage\n%%EOF",0.2);
	fclose(file);
	return 1;
}

void print_preamble(FILE *file, double max, double min)
{
	fprintf(file, "%%!PS-Adobe EPSF-3.0\n");
	fprintf(file, "%%%%BoundingBox: %g %g %g %g\n",-5.,0.,105.,100.);
	fprintf(file, "%%%%Parameters: num_steps = %d, max_val = %g, min_val = %g, alpha = %g\n\n",NN,max,min,1.0*alpha);
	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/LW 2 def\n/CW 1.5 def\n/BW 1 def\n");	
	
	fprintf(file, "\n-3 50 m\n105 0 rl\n0 0 m\n0 100 rl\n0.4 SLW\nS\n\n");
}

double rw_step(double u)
{
	if(SRW)
		// SRW
		return 2*(u>0.5)-1;
	else
	{ 
		
		// alpha-stable
		if(u<0.5)
			return -pow(2*u,-1/alpha);
		if(u>0.5)
			return pow(2*(1-u),-1/alpha);
	}
}