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

/*constants for ran11(long *)*/
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

long seed;

double ran11(long *idum);
double gasdev(long *idum);

main()
{
	int i,N;
	int k,M;
	double dt,sdt,dw;
	double t,w;
	char filename[25];

	double *db;
	double *sum,*sum2;

	FILE *outfile;

	// user io
	printf("Please input random seed [-xxx]:\n");
	scanf("%ld",&seed);
	printf("Length of data [1000]:\n");
	scanf("%d",&N);

	dt=(double)(1.0/N);
	sdt=sqrt(dt);

	// for use in calculating mean and variance of Brownian process
	M=10000;
	sum=dvector(1,N);
	sum2=dvector(1,N);

	// vector storage for Brownian increments
	db=dvector(1,N);

	for(k=1;k<=M;k++)//this outer loop will generate M different Brownian Paths
	{

	// Generate Brownian Path
	w=0.0;
	for(i=1;i<=N;i++)
		{
		dw=sdt*(double)gasdev(&seed);
		w+=dw;

		sum[i]+=w;
		sum2[i]+=w*w;

		db[i]=dw;
		}
	
	if(k<=5)
		{
		// Brownian Bridge - only plot the first 5 of the M different Paths

			sprintf(filename,"bbridge%d.dat",k);
			outfile=fopen(filename,"w");
			fprintf(outfile,"tt\tww\n");

			dw=w;
			w=0.0;
			t=0.0;
			for(i=1;i<=N;i++)
			{
				t+=dt;
				w+=db[i];
				fprintf(outfile,"%lg\t%lg\n",t,w-t*dw);
			}

			fclose(outfile);
		}


	}

// Mean and Variance - calcuated using the M different Brownian Paths
// variance is calculated by <x^2>-<x>^2

	sprintf(filename,"mean.dat");
	outfile=fopen(filename,"w");
	fprintf(outfile,"tt\tmean\tvar\n");
	fprintf(outfile,"%lg\t%lg\t%lg\n",0.0,0.0,0.0);

	t=0.0;
	for(i=1;i<=N;i++)
	{
		t+=dt;
		sum[i]/=M;
		sum2[i]/=M;

		fprintf(outfile,"%lg\t%lg\t%lg\n",t,sum[i],sum2[i]-sum[i]*sum[i]);
	}
	
	fclose(outfile);

	return 0;
}

//
// random number generators
//

double ran11(long *idum)
{
	int j;
	long k;
	static long iy=0;
	static long iv[NTAB];
	double temp;
	
	if(*idum <= 0 || !iy)
		{
		if(-(*idum) < 1)
			*idum = 1;
		else
			*idum = -(*idum);
		for(j=NTAB+7;j>=0;j--)
			{
			k=(*idum)/IQ;
			*idum=IA*(*idum-k*IQ)-IR*k;
			if(*idum<0)
				*idum += IM;
			if(j<NTAB)
				iv[j]= *idum;
			}
		iy=iv[0];
		}
	k=(*idum)/IQ;
	*idum=IA*(*idum-k*IQ)-IR*k;
	if(*idum < 0)
		*idum += IM;
	j=iy/NDIV;
	iy=iv[j];
	iv[j]= *idum;
	if((temp=AM*iy) >RNMX)
		return RNMX;
	else
		return temp;
}

double gasdev(long *idum)
{
	static int iset=0;
	static double gset;
	double fac,r,v1,v2;
    
	if  (iset == 0) {
		do {
			v1=2.0*ran11(idum)-1.0;
			v2=2.0*ran11(idum)-1.0;
			r=v1*v1+v2*v2;
		} while (r >= 1.0);
		fac=sqrt(-2.0*log(r)/r);
		gset=v1*fac;
		iset=1;
		return v2*fac;
	} else {
		iset=0;
		return gset;
	}
}