/*========================================================================*
 * Program:  Monte Carlo integral estimation
 * Author: Mark Schmelzenbach
 *========================================================================
 * Date: 01/10/98
 *=======================================================================*/

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

/* George Marsaglia's uniform random number generator */
/* (has a very large period, > 2^60, and passes Diehard tests; */
/* uses a multiply-with-carry method) */
#define s1new  ((s1=36969*(s1&65535)+(s1>>16))<<16)
#define s2new  ((s2=18000*(s2&65535)+(s2>>16))&65535)
#define UNI    ((s1new+s2new)*2.32830643708e-10)
static unsigned long s1=362436069, s2=521288629;
#define setseed(seed1,seed2) {s1=seed1;s2=seed2;}
 
int init_random()
{
  long a,b;
  a=time(NULL);
  b=0x01234567^a;
  setseed(a,b);
}

double strat(long i, long n)
{
   return ((double)i+UNI)/(double)n;
}

double uniform(double a, double b)
{
  double y,py;

  y=a+(b-a)*UNI;
  py=1/(b-a);
  return (y*y)/py;
}

double stratified(double a, double b, long i, long n)
{
  double y,py;

  y=a+(b-a)*strat(i,n);
  py=1/(b-a);
  return (y*y)/py;
}

double importance()
{
 double a,y,py;

 y=sqrt(4.0*UNI);
 py=y/2;	
 return (y*y)/py;
}

double imp_strat(long i, long n)
{
 double a,y,py;

 y=sqrt(4.0*strat(i,n));
 py=y/2;	
 return (y*y)/py;
}

int main(int argc, char *argv[])
{
 long n,i,j;
 double total[4];

 n=atol(argv[1]);
 printf("Sampling %d times.\n",n);
 init_random();
 
 for(j=0;j<10;j++) {
   for(i=0;i<4;i++)
     total[i]=0;
   for(i=0;i<n;i++) {
     total[0]+=uniform(0,2);
     total[1]+=stratified(0,2,i,n);
     total[2]+=importance();
     total[3]+=imp_strat(i,n);
   } 
   printf("%g %g %g %g\n",total[0]/(double)n,total[1]/(double)n,
	  total[2]/(double)n,total[3]/(double)n);
 }
 return 0;
}