/* The result of this program is to compute two vectors, each number of
 * which corresponds to a square on the Monopoly board.  The numbers
 * represent the long term probability of ending up on that square.
 * There are two vectors to model the two strategies of paying to
 * immediately get out of jail and staying in as long as possible.  There
 * are times in the game where one or the other is preferable.  All of
 * the rules of the game, including Chance and Community Chest cards and
 * being sent to Jail on the third double, are taken into account for
 * these results.
 * 
 * Jail is split into two spaces.  One for just landing on and one for
 * landing in.  The space that is just landed on is at the normal place
 * in the array.  The one representing being in jail is at the end of the
 * array and is treated separately.
 * 
 * In the process of calculating these two vectors, we have to compute two
 * matrices, one for prefered short jail stay and one for prefered long
 * jail stay.  These matrices hold, for each square, an array holding the
 * probability of reaching each other square on the board in one roll.
 * 
 * These matrices are used to generate the final probability vectors by
 * taking a vector initialized with a probability of 1.0 in the Go square
 * and repeatedly multiplying by the appropriate matrix until the values
 * in the vector converges.
 *
 * By Truman Collins
 * April 9, 1997
 *
 * Copyright 1997 By Truman Collins
 */

#include <stdio.h>


/* The index of the in jail square. */

const  int  jail_index = 40;

/* Indexes for the financial data for properties. */

const int prop_price   = 0;
const int house_price  = 1;
const int hotel_price  = 2;
const int rent_prop    = 3;
const int rent_1_house = 4;
const int rent_2_house = 5;
const int rent_3_house = 6;
const int rent_4_house = 7;
const int rent_hotel   = 8;

/* Indexes for financial data on the railroads. */

const int rr_cost      = 0;
const int rr_rent_1    = 1;
const int rr_rent_2    = 2;
const int rr_rent_3    = 3;
const int rr_rent_4    = 4;

/* Indexes for financial data on the utilities. */

const int util_cost    = 0;
const int util_rent_1  = 1;
const int util_rent_2  = 2;

/* Static data   */

static int      jail_long;
static double   eigenvector[41];
static double   jail_long_eigenvector[41];
static double   jail_short_eigenvector[41];
static double   simple_roll_probs[13];
static double   prob_roll_is_double[13];
static double   markov_matrix[41][41];
static double   jail_long_markov_matrix[41][41];
static double   jail_short_markov_matrix[41][41];
static int      is_a_chance_square[41];
static int      is_a_comm_chest_square[41];
static double   prob_last_two_rolls_doubles[41];

/* Computed by another program through simulation. */

static double   jail_short_prob_last_two_rolls_doubles[41] = {
    0.019165, 0.029872, 0.017706, 0.028649, 0.020851, 0.026115, 0.020261,
    0.026939, 0.019900, 0.024021, 0.022502, 0.023941, 0.022877, 0.021994,
    0.022415, 0.020188, 0.021668, 0.017492, 0.021605, 0.017042, 0.023329,
    0.019216, 0.025720, 0.020805, 0.024754, 0.021709, 0.024125, 0.023024,
    0.024075, 0.024473, 0.000000, 0.024763, 0.018811, 0.026052, 0.017904,
    0.027403, 0.020718, 0.030218, 0.018410, 0.028543, 0.000000  
};
static double   jail_long_prob_last_two_rolls_doubles[41] = {
    0.019554, 0.029881, 0.017728, 0.028759, 0.020882, 0.026610, 0.020306,
    0.027174, 0.019967, 0.024284, 0.022582, 0.024579, 0.021549, 0.022856,
    0.023235, 0.021294, 0.024598, 0.018795, 0.025794, 0.018126, 0.028700,
    0.019693, 0.031988, 0.020806, 0.033137, 0.022303, 0.031569, 0.022806,
    0.029914, 0.024055, 0.000000, 0.024347, 0.021044, 0.025596, 0.018391,
    0.027055, 0.020542, 0.029982, 0.018184, 0.028726, 0.000000  
};


int mod_index(int ind)
   {
      return(ind % 40);
   }


void initialize_variables(void)
   {
      int     i,j;
      double  one_thirtysixth = (1.0 / 36.0);
      double  prob_value;

      /* Zero simple probability array. */

      for(i = 0; i < 13; i++) {
         simple_roll_probs[i] = 0.0;
         prob_roll_is_double[i] = 0.0;
      }

      /* Zero the markov matrix and the Chance and Community Chest flags. */

      for(i = 0; i < 41; i++) {
         for(j = 0; j < 41; j++) {
            markov_matrix[i][j] = 0.0;
         }
         is_a_chance_square[i] = 0;
         is_a_comm_chest_square[i] = 0;
      }

      /* Depending on if we're computing for the player wanting to stay    */
      /* in jail or not, copy in the appropriate array for the probability */
      /* that the last two rolls were doubles for each square.  The        */
      /* probabilities are based on empirical probability computed with    */
      /* another program. */

      if(jail_long) {
         for(i = 0; i < 41; i++) {
            prob_last_two_rolls_doubles[i] = jail_long_prob_last_two_rolls_doubles[i];
         }
      } else {
         for(i = 0; i < 41; i++) {
            prob_last_two_rolls_doubles[i] = jail_short_prob_last_two_rolls_doubles[i];
         }
      }

      /* Now put in the appropriate probabilities in the simple roll probs array. */
      /* These represent the probabilities of the various rolls of two dice.  For */
      /* example, simple_roll_probs[4] is the probability of rolling a 4 with two */
      /* dice. */

      for(i = 1; i < 6; i++) {
         prob_value = i * one_thirtysixth;
         simple_roll_probs[i + 1]  = prob_value;
         simple_roll_probs[13 - i] = prob_value;
      }
      simple_roll_probs[7] = 6 * one_thirtysixth;

      /* Now fill in the array of probabilities that a roll summing to a particular */
      /* value is a double. */

      prob_roll_is_double[2]  = 1.0;
      prob_roll_is_double[12] = 1.0;
      prob_roll_is_double[4]  = 1.0 / 3.0;
      prob_roll_is_double[10] = 1.0 / 3.0;
      prob_roll_is_double[6]  = 0.2;
      prob_roll_is_double[8]  = 0.2;

      /* Note where the Chance and Community Chest squares are. */

      is_a_comm_chest_square[2] = 1;
      is_a_chance_square[7] = 1;
      is_a_comm_chest_square[17] = 1;
      is_a_chance_square[22] = 1;
      is_a_comm_chest_square[33] = 1;
      is_a_chance_square[36] = 1;
   }


void distribute_chance_probabilities(
      int     orig_square,
      int     new_square
   )
   /* There are 16 Chance cards.  10 of them send the player to */
   /* another square.  Distribute these probabilities.          */
   {
      double  one_sixteenth;

      one_sixteenth = markov_matrix[orig_square][new_square] / 16.0;
      markov_matrix[orig_square][new_square] = 6.0 * one_sixteenth;
      markov_matrix[orig_square][39] += one_sixteenth;
      markov_matrix[orig_square][5]  += one_sixteenth;
      markov_matrix[orig_square][24] += one_sixteenth;
      markov_matrix[orig_square][11] += one_sixteenth;
      markov_matrix[orig_square][0]  += one_sixteenth;
      markov_matrix[orig_square][40] += one_sixteenth;

      /* Some cards are dependent on the Chance square landed on. */

      switch(new_square) {
         case  7 :

            /* Next Railroad, 3 squares back, and next utility.   */

            markov_matrix[orig_square][15] += 2.0 * one_sixteenth;
            markov_matrix[orig_square][4]  += one_sixteenth;
            markov_matrix[orig_square][12] += one_sixteenth;
            break;

         case 22 :

            /* Next Railroad, 3 squares back, and next utility.   */

            markov_matrix[orig_square][25] += 2.0 * one_sixteenth;
            markov_matrix[orig_square][19] += one_sixteenth;
            markov_matrix[orig_square][28] += one_sixteenth;
            break;

         case 36 :

            /* Next Railroad, 3 squares back, and next utility.   */
            /* Note that for 3 squares back, we land on a         */
            /* Community Chest square, so we have to distribute   */
            /* this probability over the likelyhood that the      */
            /* the card drawn there will send him to Go or Jail.  */

            markov_matrix[orig_square][5]  += 2.0 * one_sixteenth;

            markov_matrix[orig_square][33] += one_sixteenth * 14.0 / 16.0;
            markov_matrix[orig_square][0]  += one_sixteenth / 16.0;
            markov_matrix[orig_square][40] += one_sixteenth / 16.0;

            markov_matrix[orig_square][12] += one_sixteenth;
            break;

         default :

           /* This should never happen. */

           fprintf(stderr, "Bad Chance square.  It was %d.\n", new_square);
      }
   }


void distribute_comm_chest_probabilities(
      int     orig_square,
      int     new_square
   )
   /* There are 16 Community Chest cards.  2 of them send the player to */
   /* another square.  Distribute these probabilities.                  */
   {
      double  one_sixteenth;

      one_sixteenth = markov_matrix[orig_square][new_square] / 16.0;
      markov_matrix[orig_square][new_square] = 14.0 * one_sixteenth;
      markov_matrix[orig_square][0]  += one_sixteenth;
      markov_matrix[orig_square][40] += one_sixteenth;
   }


void compute_markov_array(void)
   {
      int    i, j;
      int    new_index;
      double prob_doubled_to_jail;
      double prob_first_or_second;
      double prob_third;
      double sum_of_other_jail_probs = 0;
      double temp_prob;


      /* Put in the simple probabilities from each square.             */
      /* Correct for probability of rolling a double with two previous */
      /* doubles and being sent to the in jail square.                 */

      for(i = 0; i < 40; i++) {
         for(j = 2; j < 13; j++) {
            new_index = mod_index(j + i);
            prob_doubled_to_jail = prob_roll_is_double[j] * prob_last_two_rolls_doubles[i];
            markov_matrix[i][40] += prob_doubled_to_jail * simple_roll_probs[j];
            markov_matrix[i][new_index] +=
                    (1.0 - prob_doubled_to_jail) * simple_roll_probs[j];

            /* Now take care of the Go to Jail square. */

            if(new_index == 30) {
               markov_matrix[i][jail_index] += markov_matrix[i][30];
               markov_matrix[i][30] = 0.0;
            }

            /* Now, correct for Chance and Community Chest squares. */

            if(is_a_chance_square[new_index]) {
               distribute_chance_probabilities(i, new_index);
            }
            if(is_a_comm_chest_square[new_index]) {
               distribute_comm_chest_probabilities(i, new_index);
            }
         }
      }

      /* Now compute the probabilities for leaving the in jail square.         */

      /* We have two possibilities here.  Either we are making the computation */
      /* assuming that the player wants to spend time in jail or that he       */
      /* wants out right away. */

      if(jail_long) {

         /* Here the player waits to get out of jail until he rolls doubles */
         /* or has to get out on the third roll.  Likely when all or most   */
         /* properties have been purchased.                                 */

         /* Note that we have to separate the possibilites of having just   */
         /* landed in jail, having been there for one roll, and for two.    */
         /* The chances of being there on the first roll in jail are more   */
         /* than one third and the chances of being there on the third are  */
         /* less than one third.  Once sent to jail, there is a 36/36 chance*/
         /* that you will see your first roll, a 30/36 chance you will see  */
         /* your second, and a 25/36 chance you will see your third.  So,   */
         /* the probability that it is your third roll is 25/(36+30+25).    */

         prob_third = 25.0 / 91.0;
         prob_first_or_second = 1.0 - prob_third;
   
         /* Put in the probabilities for the third roll where we leave */
         /* the jail square reguardless of the roll.                   */
   
         for(j = 2; j < 13; j++) {
            new_index = j + 10;
            temp_prob = simple_roll_probs[j] * prob_third;
            markov_matrix[jail_index][new_index] += temp_prob;
         }
         sum_of_other_jail_probs = prob_third;
   
         /* Now put in the probabilities for the first or second roll */
         /* where we only get out on a double.  We just loop through  */
         /* even rolls, and there is a 1 in 36 chance of getting one  */
         /* of these with a double. */
   
         for(j = 2; j <= 12; j = j + 2) {
            new_index = j + 10;
            temp_prob = (1.0 / 36.0) * prob_first_or_second;
            markov_matrix[jail_index][new_index] += temp_prob;
            sum_of_other_jail_probs += temp_prob;
         }

         /* Now distribute any Chance or Community Chest squares     */
         /* that might have been landed on.                          */

         distribute_comm_chest_probabilities(jail_index, 17);
         distribute_chance_probabilities(jail_index, 22);

         /* Now the probability that we stay in jail.                */
         /* Note that it is possible to get immediately back to jail */
         /* after rolling out if you land on Community Chest and get */
         /* the Go To Jail card.  Because of that, we have to add    */
         /* to the in jail square on this next line.                 */
   
         markov_matrix[jail_index][jail_index] += 1.0 - sum_of_other_jail_probs;

      } else {

         /* Put in the probabilities for leaving jail.  Assume the player */
         /* pays to get out immediately.  Likely early in the game when   */
         /* there are still some properties to buy.                       */
   
         for(j = 2; j < 13; j++) {
            new_index = j + 10;
            markov_matrix[jail_index][new_index] += simple_roll_probs[j];
   
            /* Correct for Chance and Community Chest squares. */
   
            if(is_a_chance_square[new_index]) {
               distribute_chance_probabilities(jail_index, new_index);
            }
            if(is_a_comm_chest_square[new_index]) {
               distribute_comm_chest_probabilities(jail_index, new_index);
            }
         }
      }
   }


void compute_eigenvector(void)
   {
      int     i, j, k;
      double  last_prob_vector[41];
      double  max_diff;
      double  max_neg_diff;
      double  max_pos_diff;
      double  next_prob_vector[41];
      double  prob_diffs[41];


      /* First zero out the vectors we'll use to calculate with. */

      for(i = 0; i < 41; i++) {
         last_prob_vector[i] = 0;
         next_prob_vector[i] = 0;
         prob_diffs[i] = 0;
      }

      /* Now put a probability of one into the Go square and start       */
      /* iterating with the Markov matrix.  Each iteration will          */
      /* advance the probabilities through the vector for a single roll. */
      /* for example, after five passes, the probabilities for each      */
      /* square are the probabilities after five rolls of the dice       */
      /* with the initial position starting at Go.  Eventually, the      */
      /* resulting probability vector will settle out revealing the      */
      /* long term probabilites of landing on each square.  We compute   */
      /* an array of differences, which holds the difference for each    */
      /* element of the vector from one iteration to the next.  We use   */
      /* this to determine when to stop the iterations.                  */

      last_prob_vector[0] = 1.0;
      do {

         /* Zero new vector since we're adding into it. */

         for(i = 0; i < 41; i++) {
            next_prob_vector[i] = 0.0;
         }

         /* Compute new vector with matrix multiplication. */

         for(i = 0; i < 41; i++) {
            for(j = 0; j < 41; j++) {
               next_prob_vector[j] += markov_matrix[i][j] * last_prob_vector[i];
            }
         }

         /* Find differences. */

         max_neg_diff = 0.0;
         max_pos_diff = 0.0;
         for(i = 0; i < 41; i++) {
            prob_diffs[i] = last_prob_vector[i] - next_prob_vector[i];
            if(prob_diffs[i] < max_neg_diff) {
               max_neg_diff = prob_diffs[i];
            }
            if(prob_diffs[i] > max_pos_diff) {
               max_pos_diff = prob_diffs[i];
            }
         }
         max_diff = max_pos_diff;
         if(-max_neg_diff > max_diff) {
            max_diff = -max_neg_diff;
         }

         /* Copy to last array in preparation for next iteration. */

         for(i = 0; i < 41; i++) {
            last_prob_vector[i] = next_prob_vector[i];
         }

      } while(max_diff > 1E-15);

      /* Copy to eigenvector and markov matrix. */

      for(i = 0; i < 41; i++) {
         eigenvector[i] = next_prob_vector[i];
         if(jail_long) {
            jail_long_eigenvector[i] = eigenvector[i];
         } else {
            jail_short_eigenvector[i] = eigenvector[i];
         }
         for(j = 0; j < 41; j++) {
            if(jail_long) {
               jail_long_markov_matrix[i][j] = markov_matrix[i][j];
            } else {
               jail_short_markov_matrix[i][j] = markov_matrix[i][j];
            }
         }
      }
   }


void print_result(
      char  *title
   )
   {
      int     i;
      double  summ = 0.0;


      summ = 0.0;
      printf("\n%s\n", title);
      for(i = 0; i < 41; i++) {
         printf("%7.5f  ", eigenvector[i]);
         summ += eigenvector[i];
         if(i == 9 || i == 19 || i == 29 || i == 39) {
            printf("\n");
         }
      }
      printf("\n");
   }


void main()
   {
      /* Do first run assuming the player buys his way out of */
      /* jail immediatly.                                     */

      jail_long = 0;

      initialize_variables();

      compute_markov_array();

      compute_eigenvector();

      print_result("Short jail stay eigenvector is:");

      /* Do second run assuming the player stays in jail as long */
      /* as he can.                                              */

      jail_long = 1;

      initialize_variables();

      compute_markov_array();

      compute_eigenvector();

      print_result("Long jail stay eigenvector is:");
   }