//======================================================= file = fail.c =====
//= A simulation of machine failures                                        =
//===========================================================================
//= Problem description:                                                    =
//=   Assume that there are initially N machines in service.  At discrete   =
//=   periodic time intervals each machine has an independent probability   =
//=   p of failing.  When and if there are N_good or less machines left in  =
//=   service, no more failures will occur.  Determine the probability that =
//=   at least one machine remains in service at the completion of M time   =
//=   intervals.  Note that it is possible that there will be no machines   =
//=   in service (i.e., all have failed) at some point.                     =
//===========================================================================
//= Notes:                                                                  =
//=    1) Must set parameters as #define, which are NUM_ITER, N, N_GOOD,    =
//=       M, and P_FAIL                                                     =
//=-------------------------------------------------------------------------=
//= Example execution:                                                      =
//=   --------------------------------------------------- fail.c -----      =
//=   -  Num iterations = 1000000                                           =
//=   -  N              = 2 machines                                        =
//=   -  N_good         = 1 machines                                        =
//=   -  M              = 5 intervals                                       =
//=   -  Pr[fail]       = 0.100000 per machine                              =
//=   ----------------------------------------------------------------      =
//=   -  Pr[at least 1 machines up] = 0.966090                              =
//=   ----------------------------------------------------------------      =
//=-------------------------------------------------------------------------=
//= Build: bcc32 fail.c                                                     =
//=-------------------------------------------------------------------------=
//= Execute: fail                                                           =
//=-------------------------------------------------------------------------=
//= Author: Ken Christensen                                                 =
//=         University of South Florida                                     =
//=         WWW: http://www.csee.usf.edu/~christen                          =
//=         Email: christen@csee.usf.edu                                    =
//=-------------------------------------------------------------------------=
//= History: KJC (12/27/06) - Genesis                                       =
//===========================================================================
//----- Include files -------------------------------------------------------
#include <stdio.h>                 // Needed for printf()

//----- Defines -------------------------------------------------------------
#define GOOD                    1  // Machine is good (operational)
#define FAIL                    0  // Machine is failed
#define NUM_ITER          1000000  // Number of experiments to simulate
#define N                       2  // Number of machines to start
#define N_GOOD                  1  // Number of machines left not to fail
#define M                       5  // Number of time intervals
#define P_FAIL                0.1  // Pr[machine fail per interval]

//----- Function prototypes -------------------------------------------------
double unif(int seed);             // Function to generate unif(0,1)

//===== Main program ========================================================
void main(void)
{
  int      m[N];                   // Array of machines (GOOD or FAIL)
  int      sumGood;                // Sum of good machines after M intervals
  int      sumFinal;               // Sum of final good after NUM_ITER
  double   probFinal;              // Probability of good machines
  int      i, j, k;                // Loop counters

  // Seed the RNG
  unif(1);

  // Main simulation loop
  sumFinal = 0;
  for (i=0; i<NUM_ITER; i++)
  {
    // Set all the machines to be good
    for (j=0; j<N; j++)
      m[j] = GOOD;

    // Run for M intervals
    for (j=0; j<M; j++)
    {
      // Make machines fail
      for (k=0; k<N; k++)
        if (m[k] == GOOD)
          if (unif(0) < P_FAIL) m[k] = FAIL;

      // Sum-up number of good machines
      sumGood = 0;
      for (k=0; k<N; k++)
        if (m[k] == GOOD)
          sumGood = sumGood + 1;

      if (sumGood <= N_GOOD) break;
    }

    // Sum-up case of some machines good at end of M intervals
    if (sumGood > 0)
      sumFinal = sumFinal + 1;
  }

  // Compute probability of there being good machines at end of M intervals
  probFinal = (double) sumFinal / NUM_ITER;

  // Output the results
  printf("--------------------------------------------------- fail.c ----- \n");
  printf("-  Num iterations = %d               \n", NUM_ITER);
  printf("-  N              = %d machines      \n", N);
  printf("-  N_good         = %d machines      \n", N_GOOD);
  printf("-  M              = %d intervals     \n", M);
  printf("-  Pr[fail]       = %f per machine   \n", P_FAIL);
  printf("---------------------------------------------------------------- \n");
  printf("-  Pr[at least %d machines up] = %f  \n", N_GOOD, probFinal);
  printf("---------------------------------------------------------------- \n");
}

//=========================================================================
//= Multiplicative LCG for generating uniform(0.0, 1.0) random numbers    =
//=   - From R. Jain, "The Art of Computer Systems Performance Analysis," =
//=========================================================================
double unif(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
  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 is greater than 0
  if (seed > 0)
  {
    x = seed;
    return(0);
  }

  // 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);
}