/*      Metropolis Monte Carlo method (with p(x)=3/4(1-x^2)) for the
	integration of the function f(x)=(1-x^2)*exp(-x^2) on
	the [-1,1] interval  */

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

#define  Nt 1000
#define  Nmax 1000000
#define x1  -1.0
#define x2   1.0

double randf()
{return( (double)(rand())/(double)(RAND_MAX+1.0));
}


main()
{long int i;
 double I3, x, xn, xi, xf, w, r;
 I3=0;

 x=randf()*(x2-x1)+x1;
 if (x<0) xi=-x;
 else xi=x;
 for(i=1;i<=Nt; i++)
	{xn=randf()*(x2-x1)+x1;
	 if (xn<0) xf=-xn;
	 else xf=xn;
	 if (xf>xi) w=(1.0-xf*xf)/(1.0-xi*xi);
	 else w=1;
	 r=randf();
	 if (r<=w) {x=xn; xi=xf;}
	}

  for(i=1;i<=Nmax; i++)
	{xn=randf()*(x2-x1)+x1;
	 if (xn<0) xf=-xn;
	 else xf=xn;
	 if (xf>xi) w=(1.0-xf*xf)/(1.0-xi*xi);
	 else w=1;
	 r=randf();
	 if (r<=w) {x=xn; xi=xf;}
	 I3+=4.0/3.0*exp(-x*x);
	}

  I3=I3/(double)(Nmax);

 printf("%lf\n", I3);
 getchar();
}