/* 
   monty_hall.c
   by Martin Bazant, 2/13/02
   Intro to Modeling and Simulation

To compile type:
  
gcc monty_hall.c -o monty_hall -lm

or

g++ -o monty_hall.c monty_hall

To run for 1000 trials of 3 doors with Monty opening 1 type:

monty_hall 3 1 1000

(Without any arguments, you just get instructions.)
*/

/* INCLUDE C LIBRARIES */

#include <stdio.h>    /* library for input/output */
#include <stdlib.h>   /* library with rand() and srand() */
#include <math.h>     /* library with floor()*/
#include <time.h>     /* library to get the current time */

/* DEFINITIONS */

#define MAX_DOORS 100  /* maximum number of doors n */

/* FUNCTIONS */

int my_rand(int p,int q)
     /* my_rand(p,q) returns an integer between p and q inclusive. */
{
  /* NOTE: rand() returns a random integer between 0 and a very large
     number RAND_MAX.  rand()/(RAND_MAX+1) is between 0 and 1 (not
     inclusive), so this times |p-q|+1 is between 0 and |p-q|+1 (not
     inclusive).  Rounding down to the nearest integer gives us a
     random integer between 0 and |p-q| (inclusive).  We use this
     method over something similar like rand() % |p-q| + 1, as the
     latter uses the low order bits of the random number.  These are,
     in general, somewhat less random than the high order bits.  It is
     not usually good to treat _bits_ of a pseudo-random number as
     random.  For many purposes, the difference is academic, but for
     some applications it makes a difference.  See _Numerical Recipes
     in C_ for some discussion and more references. */

  int delta;
  int shift_from_zero;

  if(p > q) {
    delta = p - q;
    shift_from_zero = q;
  } else {
    delta = q - p;
    shift_from_zero = p;
  }

  return shift_from_zero+(int)floor((float)rand()
				    /((float)RAND_MAX+1.0)*(float)(delta+1));
}

/* MAIN PROGRAM */

int main(int argc, char *argv[])
{

  /* VARIABLE DECLARATIONS */

  int m,n,trials=100;
  int wins=0;
  int num_open,test_door,found_door,i,j;
  int door[MAX_DOORS];   /* array entries for doors = 0 if unopened
			    and no prize, = 1 if prize, =2 if opened */
  float p,p_stay;  

  /* PROCESS COMMAND-LINE INPUT */

  if(argc < 3) {
    printf("  USAGE:  monty_hall <n> <m> [<trials>] \n");
    printf("  n = number of doors available\n");
    printf("  m = number of doors opened by Monty\n");
    printf("  trials = number of computer experiments\n");
    exit(1);
  }
  n = atoi(argv[1]);    /* convert command-line characters to integers */
  if(n < 3 || n > MAX_DOORS) {
    printf("Error:  need 3 < n < MAX_DOORS.\n");
    exit(1);    
  }
  m = atoi(argv[2]);
  if(m > n-2 || m < 0) {
    printf("Error: need 0 <= m <= n-2.\n");
    exit(1);    
  }
  if(argc > 3) {
    trials = atoi(argv[3]);
  }
  
  /* SEED RANDOM NUMBER GENERATOR */

  srand(time(NULL));  /* "Seed" the generator with the current time
			 to get a fresh random number. (If not, you 
			 always get the same sequence starting 
			 from seed=1 when you run the program again,
			 which means you are not getting new data 
			 if you run the program more than once.) */


  /****** LOOP OVER COMPUTER EXPERIMENTS ******/
  
  for(i=0;i<trials;i++) {

    /* PUT THE PRIZE BEHIND A RANDOM DOOR 0,...,n-1 */

    for(j=0;j<n;j++) door[j] = 0;   /* initialize door array */
    door[my_rand(0,n-1)]++;     /* label a door with the prize */

    /* ASSUME YOU PICKED DOOR 0. LET MONTY OPEN m OTHER DOORS. 
       Note: this "blind" algorithm is not efficient if m is large!  */ 
    
    num_open = 0;
    while(num_open < m){
      test_door = my_rand(1,n-1);      /* pick a door */
      if(door[test_door] < 1) {        /* if not open and no prize */
	door[test_door] = 2;           /* label the door as open */
	num_open++;
      }
    }
   
    /* SWITCH TO ANOTHER RANDOM UNOPENED DOOR. */
    
    found_door = 0;
    while(!found_door) {    
      test_door = my_rand(1,n-1);       /* pick a door */
      if(door[test_door]<2) found_door++;
    }

    /* SEE IF YOU GET THE PRIZE! */

    wins += door[test_door];     
    
  }
 
  /******** END OF LOOP OVER EXPERIMENTS **********/
  
  /* ESTIMATE OF PROBABILITY OF WINNING IF YOU SWITCH DOORS */

  p = (float)wins/(float)trials;
  

  /* OUTPUT RESULTS */

  printf("\n Percentage of wins in %d trials = %f.\n",trials,100.0*p);

  p_stay = 1.0/(float)n;
  
  if(p > p_stay) {
    printf("\n Switching doors appears to be helpful,\n");
    printf(" since you expect to win only %f \n",100*p_stay);
    printf(" of the time if you stick with your original choice.\n");
  }
    
}