//this program integrates f(x) using the Monte Carlo with importance sampling


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

//constants for ran1 and 2(int *)
#define M1 259200
#define IA1 7141
#define IC1 54773
#define RM1 (1.0/M1)
#define M2 124456
#define IA2 8121
#define IC2 28411
#define RM2 (1.0/M2)
#define M3 243000
#define IA3 4561
#define IC3 51349

//accuracy for integrator
#define EPS 1.0e-4
#define MaxN 1000000000

double pi;
int seed1,seed2;//seeds for random generator with uniform distribution

double integrateMC(double xstart, double xend, double tol, long *count);

double func(double xx);
float rann1(int *idum);
float rann2(int *idum);

main()
	{
	int i,N;
	double a,b;//range of x to be plotted
	double x,mcval;
	double xstep;
	long counter;

	FILE *outfile;

	pi=4.0*atan(1.0);

	seed1=-1;//initialize random generator1
	seed2=-3;//initialize random generator2

	a=0.0;
	b=2.0;
	N=100;
	xstep=(b-a)/(double)N;

	outfile=fopen("importanceMC.dat","w");
	fprintf(outfile,"xx\tmcval\n");
	fprintf(outfile,"%lg\t%lg\n",0,0);//the zero point
	x=xstep;//start out near zero
	for(i=1;i<=N;i++)
	{
		mcval=integrateMC(0,x,EPS,&counter);
		fprintf(outfile,"%lg\t%lg\n",x,mcval);
		printf("%lg\t%ld\n",x,counter);

		x+=xstep;
	}

	return 0;
	}

//
//
//Integrating Routinues
//
//


// Monte Carlo Integrator
double integrateMC(double xstart, double xend, double tol, long *count)
	{
	long Num,i;
	double u,v;
	double val1,val2;
	int inflag;
	

	Num=10000;
	val1=val2=0.0;

	//calculate val1 with Num
	for(i=1;i<=Num;i++)
		{
		//find a random points with [0,x]X[0,1]
		u=rann1(&seed1)*xend;
		v=rann2(&seed2);

		//rejection step
		if(v<=func(u))
			val1+=1.0;
		}
	val1=val1*xend/(double)Num;
	
	//calculate val2 with a different set of Num
	inflag=1;
	do
		{
		for(i=1;i<=Num;i++)
			{
			//find a random points with [-1,1]X[-1,1]
			u=rann1(&seed1)*xend;
			v=rann2(&seed2);

			//rejection step
			if(v<=func(u))
				val2+=1;
			}
		val2=val2*xend/(double)Num;

		//check for convergence
		if(fabs((val1-val2)/val1)<tol)
			{
			inflag=0;
			*count=2*Num;
			}
		else
			{
			val1=(val1+val2)/2.0;
			
			val2=0.0;
			Num+=Num;
			}
		}
	while (inflag && Num<MaxN);

	
	if(Num<MaxN)
		return (val1+val2)/2.0;
	else
		{
		*count=Num;
		return val1;
		}

	}


//Evaluation function of the integrant
double func(double xx)
	{
	double temp;

	temp=sin(1.0/xx);
	return temp*temp;
	}


//
//
//Random Number generators
//
//
//to generate uniform distribution in 2D, we used two identical
//	INDEPENDENT generators ran1 and ran2
float rann1(int *idum)
{
	static long ix1,ix2,ix3;
	static float r[98];
	float temp;
	static int iff=0;
	int j;
	
	extern void nrerror(char []);

	if (*idum < 0 || iff == 0)
		{
		iff=1;
		ix1=(IC1-(*idum)) % M1;
		ix1=(IA1*ix1+IC1) % M1;
		ix2=ix1 % M2;
		ix1=(IA1*ix1+IC1) % M1;
		ix3=ix1 % M3;
		for (j=1;j<=97;j++)
			{
			ix1=(IA1*ix1+IC1) % M1;
			ix2=(IA2*ix2+IC2) % M2;
			r[j]=(ix1+ix2*RM2)*RM1;
			}
		*idum=1;
		}
	ix1=(IA1*ix1+IC1) % M1;
	ix2=(IA2*ix2+IC2) % M2;
	ix3=(IA3*ix3+IC3) % M3;
	j=1 + ((97*ix3)/M3);
	if (j > 97 || j < 1)
		nrerror("ran1: This cannot happen.");
	temp=r[j];
	r[j]=(ix1+ix2*RM2)*RM1;
	return temp;
}


float rann2(int *idum)
{
	static long ix1,ix2,ix3;
	static float r[98];
	float temp;
	static int iff=0;
	int j;
	
	extern void nrerror(char []);

	if (*idum < 0 || iff == 0)
		{
		iff=1;
		ix1=(IC1-(*idum)) % M1;
		ix1=(IA1*ix1+IC1) % M1;
		ix2=ix1 % M2;
		ix1=(IA1*ix1+IC1) % M1;
		ix3=ix1 % M3;
		for (j=1;j<=97;j++)
			{
			ix1=(IA1*ix1+IC1) % M1;
			ix2=(IA2*ix2+IC2) % M2;
			r[j]=(ix1+ix2*RM2)*RM1;
			}
		*idum=1;
		}
	ix1=(IA1*ix1+IC1) % M1;
	ix2=(IA2*ix2+IC2) % M2;
	ix3=(IA3*ix3+IC3) % M3;
	j=1 + ((97*ix3)/M3);
	if (j > 97 || j < 1)
		nrerror("ran1: This cannot happen.");
	temp=r[j];
	r[j]=(ix1+ix2*RM2)*RM1;
	return temp;
}