#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>

/* Do Monte Carlo integral of a function over distance (0,1)
 */

#include "mersenne.h"

char usage[]  = " usage: mc_integrate N N2\n\
  Integrates function using N mc-hits, repeats N2 times\n\
  Prints all results (with errors) to standard output\n"; 

/* Integration range from XMIN to XMAX
 */

#define XMIN 0.0
#define XMAX 1.0

/* the function to be integrated defined here
*/
double f(double x)
{
  return(1.5 * (1.0 - x*x));
  //return(1.0/sqrt(x));
  //return ( pow( x, (-1.0/3.0) ) + 0.1*x );
  //return( x*x );
}

/*
 * This routine does one monte carlo integration.
 * Returns the mc-estimate of the value integral, and err
 * contains the estimate of the error.
 */

double mcint(int n, double *err)
{
  register int i;
  register double s,s2,t;

  s = s2 = 0;

  for (i=0; i<n; i++) {
    t = f( XMIN + (XMAX-XMIN)* mersenne() );
    s += t;
    s2 += t*t;
  }

  s2 /= n;
  s  /= n;
  *err = sqrt( (s2 - s*s) / (n-1) );
  return( s );
}
 
/* The main routine
 */

int main(int argc,char *argv[])
{
  struct timeval tv;
  struct timezone tz;
  double res,err;
  int loops,n,i;

  /* require correct number of arguments 
   */
  if (argc != 3) {
    fprintf(stderr,usage);
    exit(0);
  }

  n     = atoi(argv[1]);    /* first int */
  loops = atoi(argv[2]);    /* 2nd int */
  
  /* seed the random number using system time
   * tv.tv_sec contains time in seconds, tv.tv_usec in microseconds
   */

  gettimeofday( &tv, &tz );
  seed_mersenne( tv.tv_sec + tv.tv_usec );
  srand48( tv.tv_sec + tv.tv_usec );


  for (i=0; i<100000; i++) mersenne();

  for (i=0; i<loops; i++) {
    res = mcint(n,&err);
    printf("%d  %g  %g\n",i,res,err);
  }
 
  return(1);
}