//=================================================== file = blockError.c =====
//= A Monte Carlo simulation of independent bit errors in a block             =
//=   - For Assignment #3, problem #3                                         =
//=============================================================================
//= Notes: See problem #3 in assignment #3                                    =
//=---------------------------------------------------------------------------=
//= Build: bcc32 blockError.c                                                 =
//=---------------------------------------------------------------------------=
//= Execute: blockError                                                       =
//=---------------------------------------------------------------------------=
//= History: KJC (05/28/13) - Genesis                                         =
//=============================================================================
//----- Include files ---------------------------------------------------------
#include <stdio.h>             // Needed for printf()

//----- Constants -------------------------------------------------------------
#define    P_ERROR     0.01    // Probability of a bit error
#define  BLOCK_LEN       40    // Number of bits in a block
#define NUM_TRIALS 10000000    // Number of trials to run

//----- Function prototypes ---------------------------------------------------
double rand_val(int seed);     // RNG

//===== Main program ==========================================================
int main(void)
{
  int    error_count;          // Error count for a message
  int    error[BLOCK_LEN + 1]; // Error count vector
  int    i, j;                 // Loop counters

  // Initialization error vector to zero
  for (i=0; i<BLOCK_LEN; i++)
    error[i] = 0;

  // Seed the RNG
  rand_val(1);

  // Main simulation loop
  printf("Simulation is running for %d trials...  \n", NUM_TRIALS);
  for (i=0; i<NUM_TRIALS; i++)
  {
    error_count = 0;
    for (j=0; j<BLOCK_LEN; j++)
      if (rand_val(0) < P_ERROR) error_count++;

    // Increment appropriate element in error vector
    error[error_count]++;
  }

  // Determine and output probabilities of interest
  printf("  Pr[0 errors] = %f \n", (double) error[0] / NUM_TRIALS);
  printf("  Pr[1 errors] = %f \n", (double) error[1] / NUM_TRIALS);
  printf("  Pr[2 errors] = %f \n", (double) error[2] / NUM_TRIALS);
  printf("  Pr[3 errors] = %f \n", (double) error[3] / NUM_TRIALS);

  return(0);
}

//=========================================================================
//= Multiplicative LCG for generating uniform(0.0, 1.0) random numbers    =
//=   - x_n = 7^5*x_(n-1)mod(2^31 - 1)                                    =
//=   - With x seeded to 1 the 10000th x value should be 1043618065       =
//=   - From R. Jain, "The Art of Computer Systems Performance Analysis," =
//=     John Wiley & Sons, 1991. (Page 443, Figure 26.2)                  =
//=   - Use seed == 0 to generate random value, use seed > 0 for seeding  =
//=========================================================================
double rand_val(int seed)
{
  const long  a =      16807;  // Multiplier
  const long  m = 2147483647;  // Modulus
  const long  q =     127773;  // m div a
  const long  r =       2836;  // m mod a
  static long x;               // Random int value (this is the seed)
  long        x_div_q;         // x divided by q
  long        x_mod_q;         // x modulo q
  long        x_new;           // New x value

  // Seed the RNG
  if (seed > 0) x = seed;

  // RNG using integer arithmetic
  x_div_q = x / q;
  x_mod_q = x % q;
  x_new = (a * x_mod_q) - (r * x_div_q);
  if (x_new > 0)
    x = x_new;
  else
    x = x_new + m;

  // Return a random value between 0.0 and 1.0
  return((double) x / m);
}