// "rejfree.c": a C program to illustrate Example 2 of "rejectionfree.tex".
// by Jeffrey S. Rosenthal, May 2019

#define METROPOLIS 0  /* 1 for Metropolis, or 0 for uniform rejection-free. */

#define MAXSTATE 10000
#define INITSTATE 100
#define CYCLE 4
#define PRINTFREQ 50

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

double mainpow = 8.0/9;
double rempow = 2.0;

int main()
{
    int x = INITSTATE, y, i, numits;
    double pi(), L, R, leftprob;
    int seedrnd(), ifloor(), gettimeofday();

/*
for (i=0; i<100; i++)
  printf("pi[%d] = %f\n", i, pi(i));
exit(0);
*/

    seedrnd();
    numits = 0;
    while (1) {
if (METROPOLIS) {
      if (drand48()<0.5)
        y = x-1;
      else
        y = x+1;
      if (drand48()<pi(y)/pi(x))
        x = y;
} else {
      L = pi(x-1)/pi(x);
      R = pi(x+1)/pi(x);
      if ( (R<1) && (L<R) )
        leftprob = L/R/2;
      else if ( (L<1) && (R<L) )
        leftprob = 1 - R/L/2;
      else if (L<1)
        leftprob = L/2;
      else if (R<1)
        leftprob = 1 - R/2;
      else
        leftprob = 0.5;
      if (drand48() < leftprob)
        x = x-1;
      else
        x = x+1;
}
      if (numits == PRINTFREQ * (numits/PRINTFREQ))
        printf("%d\n", x);
      numits++;

      if ((x==0)||(x==MAXSTATE)) {
        printf("REACHED  x=%d  AFTER %d ITERATIONS\n", x, numits);
	exit(0);
      }
    }

    printf("%f\n", pi(x));
    exit(0);
}

double pi(int xx)
{
  int main = xx/CYCLE;
  int rem = xx - CYCLE*main;
  if (xx<0)
    return(0.0);
  else if (xx==0)
    return(100.0);
  else
    return (pow(mainpow,1.0*main) * pow(rempow,rem));
}

int seedrnd()
{
    int seed, gettimeofday(), ifloor();
    struct timeval tmptv;
    gettimeofday (&tmptv, (struct timezone *)NULL);
    seed = (int) (tmptv.tv_usec - 1000000 * ifloor ( ( (double) tmptv.tv_usec )
                                  / ( (double) 1000000.0 ) ) );
    srand48(seed);
}

int ifloor(double xxx)  /* returns floor(xxx) as an integer */
{
    return((int)floor(xxx));
}