/* The following C program finds the value of A (specifically for our frog
   example) used in eigenvalue comparison in Diaconis and Saloff-Coste
   (1993), "Comparison theorems for reversible Markov chains."  Ann. Appl.
   Prob. 3, 696-730 */

#define NODES 200 /* the size of the state space */
#define M1 10     /* P1 has edges between nodes within M1 of each other */
#define M2 9      /* P2 has edges between nodes within M2 of each other */
#define BETA 10   /* the target distribution is proportional to
                     (dist(x,0)+1)^-BETA where dist(x,0) is the distance
                     from node x to node 0 */

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

int findLength(int, int, int, int);
int distance (int);

int main()
{
  int min, i, j, z, w, x, y, length;
  int searchstart, searchstart2, searchend, searchend2;
  double alpha, max=0, sum, denom=0;
  double P1[NODES][NODES], P2[NODES][NODES], pi1[NODES], pi2[NODES];

  for (i=0; i<NODES; i++)
  {
    min = distance(i);
    denom = denom + pow(min+1.0, -BETA);
  }
  for (i=0; i<NODES; i++)
  {
    min = distance(i);
    pi1[i] = pow(min+1.0, -BETA)/denom;
    pi2[i] = pi1[i]; /* the stationary distributions of P2 and P1 */
  }

  for (i=0; i<NODES; i++)
  {
    P1[i][i] = 1.0/(2*M1+1);
    P2[i][i] = 1.0/(2*M2+1);
    for (j=0; j<NODES; j++)
    {
      if (pi1[j] >= pi1[i])
        alpha = 1;  /* alpha is the probability that the proposed state, j,
                       is accepted given that we are in state i; this is
                       the same for P1 and P2 since they have the same
                       stationary distribution.  However, the proposals
                       are different. */
      else
        alpha = pi1[j]/pi1[i];

      if ((0 < abs(i-j)) && (abs(i-j) < M1+1))
      {
        P1[i][j] = alpha/(2*M1+1);
        P1[i][i] = P1[i][i] + (1-alpha)/(2*M1+1);
      }
      else if (abs(i-j) >= NODES-M1) /* the highest nodes are connected to
                                        the lowest nodes */
      {
        P1[i][j] = alpha/(2*M1+1);
        P1[i][i] = P1[i][i] + (1-alpha)/(2*M1+1);
      }
      else if ((M1+1 <= abs(i-j)) && (abs(i-j) < NODES-M1))
        P1[i][i] = 0;

      if ((0 < abs(i-j)) && (abs(i-j) < M2+1))
      {
        P2[i][j] = alpha/(2*M2+1);
        P2[i][i] = P2[i][i] + (1-alpha)/(2*M2+1);
      }
      else if (abs(i-j) >= NODES-M2)
      {
        P2[i][j] = alpha/(2*M2+1);
        P2[i][i] = P2[i][i] + (1-alpha)/(2*M2+1);
      }
      else if ((M2+1 <= abs(i-j)) && (abs(i-j) < NODES-M2))
        P2[i][j] = 0;
    }
  }

  for (z=0; z<NODES; z++)
  {
    /* an edge of P cannot be outside the following boundaries */
    searchstart = (z+NODES-M1) % NODES;
    searchend = (z+M1) % NODES;
    if (searchstart > searchend)
      searchend += NODES; /* allows wrapping around from NODES-1 to 0 */
    for (i=searchstart; i<=searchend; i++)
    {
      w = i % NODES;
      sum = 0;
      for (x=0; x<NODES; x++)
      {
        searchstart2 = (x+NODES-M2) % NODES;
        searchend2 = (x+M2) % NODES;
        if (searchstart2 > searchend2)
          searchend2 += NODES;
        for (j=searchstart2; j<=searchend2; j++)
        {
          y = j % NODES;
          length = findLength(z, w, x, y);
          sum = sum + length*pi2[x]*P2[x][y];
        }
      }
      sum = sum/(pi1[z]*P1[z][w]);
      if (sum > max)
        max = sum;
    }
  }
  printf ("A is %g\n", max);
  exit(0);
}

/* returns the path length from x to y if (z,w) was one of the steps on
   the path, otherwise returns 0 */
int findLength(int z, int w, int x, int y)
{
  int left, right, indicator, rem, min, len;

  indicator = 0; /* changes to 1 when (z,w) is found to be a step on the
                    path from x to y */
  /* the edge comes by moving clockwise, e.g. 1->2*/
  if (((x<y) && (y-x < NODES/2)) || ((y<x) && (x-y >= NODES/2)))
  {
     left = x;
     right = (left + M1) % NODES;
     while ((abs(left-y) > M1) && (abs(left-y)<NODES-M1))
     {
        if ((left==z) && (right==w))
        {
           indicator = 1;
           left = y; /* ending the while loop */
        }
        right = (right+M1) % NODES; /* updating the values for the next
                                                       iteration */
        left = (left+M1) % NODES;
     }
     if ((left==z) && (y==w))
        indicator = 1;
  }

  /* the edge comes by moving counterclockwise, e.g. 2->1 */
  else if (((x<y) && (y-x >= NODES/2)) || ((y<x) && (x-y < NODES/2)))
  {
     right = x;
     left = (x+NODES-M1) % NODES;
     while ((abs(right-y)>M1) && (abs(right-y)<NODES-M1))
     {
        if ((right==z) && (left==w))
        {
           indicator = 1;
           right = y;
        }
        right = (right+NODES-M1) % NODES;
        left = (left+NODES-M1) % NODES;
     }
     if ((right==z) && (y==w))
        indicator = 1;
  }

  min = distance(abs(x-y)); /* the shortest distance from x to y */
  rem = min % M1;
  if (rem==0)
    len = indicator*min/M1; /* len is the quickest way to get from x to y
                               jumping at most M1 nodes at time, but only
                               if indicator=1 */
  else
    len = indicator*(min/M1+1);
  return len;
}

/* finds the shortest distance to node 0 from node i*/
int distance(int i)
{
  if (NODES-i < i)
    return NODES-i;
  return i;
}