/*************************************************************************/
/* 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 nsize 10000
#define nskip 15
#define ntotal nsize*nskip
#define neq 10000
#define g(x)  (x*x/(exp(x*x)-1)) 
#define weight(x)  (exp(x*x)-1)
int seed, iaccept;
double ranf();
double x, h=0.4, w;

// An example of Monte Carlo simulation with the
// Metropolis scheme with integrand f(x) = x*x.

main() {
  int i;
  int t1, t2, t3, t4, t5, t6;
  time_t ts, tp;
  struct tm t;
  double z = 0.46265167, f, s0, ds, accept;
  void metropolis();
  extern int seed, iaccept;
  extern double x, h, w;

// 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;

  x = ranf();
  w = weight(x);
  iaccept = 0;
  for (i=0; i<neq; ++i) metropolis();

  s0 = 0;
  ds = 0;
  iaccept = 0;
  for (i=0; i<ntotal; ++i) {
    metropolis();
    if (i%nskip == 0) {
      f = g(x);
      s0 += f;
      ds += f*f;
    }
  }
  s0 /= nsize;
  ds /= nsize;
  ds = sqrt(fabs(ds-s0*s0)/nsize);
  s0 *= z;
  ds *= z;
  accept = 100*iaccept;
  accept /= ntotal;
  printf("S = %16.8lf +- %16.8lf\n", s0, ds);
  printf("Accept rate = %16.8lf\n", accept);
}

void metropolis() {
  extern int seed, iaccept;
  extern double x, h, w;
  double wnew, xold = x;
  xold = x;
  x += 2*h*(ranf()-0.5);
  if ((x<0) || (x>1)) x = xold;
  else {
    wnew = weight(x);
    if (wnew > (w*ranf())) {
      w = wnew;
      ++iaccept;
    }
    else x = xold;
  }
}

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;
}