/*************************************************************************/
/* Program file name: monte.c                                            */
/*                                                                       */
/* © Tao Pang 2006                                                       */
/*                                                                       */
/* Last modified: July 3, 2006                                           */
/*                                                                       */
/* (1) This C program is created for the book, "An Introduction to       */
/*     Computational Physics, 2nd Edition," written by Tao Pang and      */
/*     published by Cambridge University Press on January 19, 2006.      */
/*                                                                       */
/* (2) No warranties, express or implied, are made for this program.     */
/*                                                                       */
/*************************************************************************/

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

#define n 1000000
int seed;

// An example of integration with direct Monte
// Carlo scheme with integrand f(x) = x*x.

main() {
  int i, m;
  int t1, t2, t3, t4, t5, t6;
  time_t ts, tp;
  struct tm t;
  double x, f, s0, ds;
  double ranf();
  extern int seed;

// Initiate the seed from the current time
  ts = time(&tp);
  t  = *gmtime(&tp);
  t1 = t.tm_sec;
  t2 = t.tm_min;
  t3 = t.tm_hour;
  t4 = t.tm_mday;
  t5 = t.tm_mon;
  t6 = t.tm_year;
  seed = t6+70*(t5+12*(t4+31*(t3+23*(t2+59*t1))));
  if (seed%2 == 0) seed = seed-1;


  s0 = 0;
  ds = 0;
  for(i=0; i<n; ++i) {
    x = ranf();
    f = x*x;
    s0 += f;
    ds += f*f;
  }
  s0 /= n;
  ds /= n;
  ds = sqrt(fabs(ds-s0*s0)/n);
  printf("S = %16.8lf +- %16.8lf\n", s0, ds);
}

double ranf() {
// Function to generate a uniform random number in [0,1]
// following x(i+1)=a*x(i) mod c with a=pow(7,5) and
// c=pow(2,31)-1.  Here the seed is a global variable.

  const int a = 16807,  c = 2147483647, q = 127773, r = 2836;
  int l, h, t;
  double cd = c;
  extern int seed;

  h = seed/q;
  l = seed%q;
  t = a*l - r*h;
  if (t > 0) seed = t;
  else seed = c + t;
  return seed/cd;
}