//======================================================= file = mm1_qdc.c =====
//=  A Queue Departure Computation M/M/1 queue simulation                      =
//==============================================================================
//=  Notes:                                                                    =
//=   1) Uses the Queue Departure Computation method found used queuecomputer  =
//=      in R                                                                  =
//=   2) The values of SIM_TIME, ARR_TIME, and SERV_TIME need to be set.       =
//=----------------------------------------------------------------------------=
//=  Build: Standard C build                                                   =
//=----------------------------------------------------------------------------=
//=  Execute: mm1_qdc                                                          =
//=----------------------------------------------------------------------------=
//=  History: KJC (12/08/19) - Genesis (from mm1.c)                            =
//==============================================================================
//----- Include files ----------------------------------------------------------
#include <stdio.h>                 // Needed for printf()
#include <stdlib.h>                // Needed for max(), exit(),and rand()
#include <math.h>                  // Needed for log()

//----- Constants --------------------------------------------------------------
#define SIM_TIME   1.0e7           // Simulation time
#define ARR_TIME   1.11111111      // Mean time between arrivals
#define SERV_TIME  1.00            // Mean service time

//----- Function prototypes ----------------------------------------------------
double rand_val(int seed);         // RNG for unif(0,1)
double exponential(double x);      // Generate exponential RV with mean x

//===== Main program ===========================================================
int main(void)
{
  double   end_time = SIM_TIME;    // Total time to simulate
  double   Ta = ARR_TIME;          // Mean time between arrivals
  double   Ts = SERV_TIME;         // Mean service time
  double   time;                   // Simulation time
  double   interarrival_time;      // Current interarrival time
  double   service_time;           // Current service time
  double   arrival_time_now;       // Current arrival time
  double   arrival_time_last;      // Last arrival time
  double   completion_time_now;    // Current completion time
  double   completion_time_last;   // Last completion time
  double   enter_service_time;     // Enter servive time
  double   residence_time;         // Residence time (wait + service)
  double   sum_residence_time;     // Sum of residence time
  double   sum_service_time;       // Sum of service time
  unsigned int c = 0;              // Number of service completions
  double   x;                      // Throughput
  double   u;                      // Utilization
  double   w;                      // Mean residence time
  double   l;                      // Mean number in the system

  // Seed the RNG
  rand_val(1);

  // Initialize key variables for QDC
  arrival_time_last = completion_time_last = 0.0;
  sum_residence_time = sum_service_time = 0.0;

  // Main simulation loop
  while (completion_time_last < end_time)
  {
    // Generate interarrival time and service time
    interarrival_time = exponential(Ta);
    service_time = exponential(Ts);

    // QDC steps
    arrival_time_now = arrival_time_last + interarrival_time;
    arrival_time_last = arrival_time_now;
    enter_service_time = max(completion_time_last, arrival_time_now);
    completion_time_now = enter_service_time + service_time;
    completion_time_last = completion_time_now;
    residence_time = completion_time_now - arrival_time_now;

    // Sum the service_time and residence_time
    sum_service_time = sum_service_time + service_time;
    sum_residence_time = sum_residence_time + residence_time;

    // increment completions
    c++;
  }

  // Set time to completion_time_last
  time = completion_time_last;

  // Compute outputs
  x = c / time;                    // Compute throughput rate
  u = (sum_service_time / c) * x;  // Compute utilization
  w = sum_residence_time / c;      // Compute mean residence of system time
  l = w / Ta;                      // Little's Law for mean number in system

  // Output results
  printf("=============================================================== \n");
  printf("=          *** Results from QDC M/M/1 simulation ***          = \n");
  printf("=============================================================== \n");
  printf("=  Total simulated time         = %f sec         \n", time);
  printf("=============================================================== \n");
  printf("=  INPUTS:                                       \n");
  printf("=    Mean time between arrivals = %f sec         \n", Ta);
  printf("=    Mean service time          = %f sec         \n", Ts);
  printf("=============================================================== \n");
  printf("=  OUTPUTS:                                      \n");
  printf("=    Number of completions      = %d cust        \n", c);
  printf("=    Throughput rate            = %f cust/sec    \n", x);
  printf("=    Server utilization         = %f %%          \n", 100.0 * u);
  printf("=    Mean number in system      = %f cust        \n", l);
  printf("=    Mean residence time        = %f sec         \n", w);
  printf("=============================================================== \n");

  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)                  =
//=   - Seed the RNG if seed > 0, return a unif(0,1) if seed == 0         =
//=========================================================================
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 (seed is set to 1)
  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);
}

//==============================================================================
//=  Function to generate exponentially distributed RVs using inverse method   =
//=    - Input:  x (mean value of distribution)                                =
//=    - Output: Returns with exponential RV                                   =
//==============================================================================
double exponential(double x)
{
  double z;                    // Uniform random number from 0 to 1

  // Pull a uniform RV (0 < z < 1)
  do
  {
    z = rand_val(0);
  }
  while ((z == 0) || (z == 1));

  return(-x * log(z));
}