//===================================================== file = chain1.c =====
//=  Iterative solver for M/M/1/K Markov chain                              =
//=   - Assumes integer valued constant arrival and service rates           =
//===========================================================================
//=  Notes:                                                                 =
//=    1) Solves for M/M/1/K steady state probailities using a novel (?)    =
//=       iterative method                                                  =
//=    2) Must manually set K, LAMBDA, and MU values in #defines            =
//=    3) Assumes that LAMBDA and MU are integer so that proper fractions   =
//=       can be output                                                     =
//=    4) Notation is pi[i][j] where i is size of chain (i+1 states) and j  =
//=       is index for state (0,... i).                                     =
//=    5) Output is to stdout                                               =
//=-------------------------------------------------------------------------=
//= Example execution: (for K = 3, LAMBDA = 3, and MU = 4)                  =
//=                                                                         =
//=   -------------------------------------------------------------         =
//=   K = 3   Lambda = 3   Mu = 4                                           =
//=   -------------------------------------------------------------         =
//=   pi[ 0][ 0] =    1 /    1 	 = 1.000000                                 =
//=   -------------------------------------------------------------         =
//=   pi[ 1][ 0] =    4 /    7 	 = 0.571429                                 =
//=   pi[ 1][ 1] =    3 /    7 	 = 0.428571                                 =
//=   -------------------------------------------------------------         =
//=   pi[ 2][ 0] =   16 /   37 	 = 0.432432                                 =
//=   pi[ 2][ 1] =   12 /   37 	 = 0.324324                                 =
//=   pi[ 2][ 2] =    9 /   37 	 = 0.243243                                 =
//=   -------------------------------------------------------------         =
//=   pi[ 3][ 0] =   64 /  175 	 = 0.365714                                 =
//=   pi[ 3][ 1] =   48 /  175 	 = 0.274286                                 =
//=   pi[ 3][ 2] =   36 /  175 	 = 0.205714                                 =
//=   pi[ 3][ 3] =   27 /  175 	 = 0.154286                                 =
//=   -------------------------------------------------------------         =
//=-------------------------------------------------------------------------=
//=  Build: bcc32 chain1.c, gcc chain1.c                                    =
//=-------------------------------------------------------------------------=
//=  Execute: chain1                                                        =
//=-------------------------------------------------------------------------=
//=  Author: Kenneth J. Christensen                                         =
//=          University of South Florida                                    =
//=          WWW: http://www.csee.usf.edu/~christen                         =
//=          Email: christen@csee.usf.edu                                   =
//=-------------------------------------------------------------------------=
//=  History: KJC (01/01/02) - Genesis                                      =
//===========================================================================

//----- Include files -------------------------------------------------------
#include <stdio.h>                 // Needed for printf()

//----- Defines -------------------------------------------------------------
#define         K     3            // Size of chain (states = 0, 1,... K)
#define    LAMBDA     3            // Arrival rate (assumed to be an int)
#define        MU     4            // Service rate (assumed to be an int)

//----- Function prototypes -------------------------------------------------
int pow_int(int x, int n);         // Returns integer x raised to n

//===========================================================================
//=  Main program                                                           =
//===========================================================================
void main(void)
{
  double       pi[K+1][K+1];     // Steady state probabilities
  int          num[K+1][K+1];    // Numerator of pi
  int          den[K+1];         // Demoninator of pi
  int          i, j;             // Loop counters

  // Output initialized values
  printf("------------------------------------------------------------- \n");
  printf("K = %d   Lambda = %d   Mu = %d \n", K, LAMBDA, MU);

  // Initialize numerator and denominator for single state case
  //  - Note that (obviously) pi[0][0] = 1.0
  num[0][0] = 1;
  den[0] = 1;

  // Compute demoninators and numerators for states 1 to K
  for (i=1; i<=K; i++)
  {
    // Compute the demoninator
    den[i] = (MU * den[i-1]) + pow_int(LAMBDA, i);

    // Compute the numerators
    for (j=0; j<i; j++)
      num[i][j] = MU * num[i-1][j];

    // Compute the numerator for the new ith state
    num[i][i] = LAMBDA * num[i-1][i-1];
  }

  // Compute all pi values
  for (i=0; i<=K; i++)
    for (j=0; j<=i; j++)
      pi[i][j] = (double) num[i][j] / den[i];

  // Output all results
  for (i=0; i<=K; i++)
  {
    printf("------------------------------------------------------------- \n");
    for (j=0; j<=i; j++)
      printf("pi[%2d][%2d] = %4d / %4d \t = %f \n",
        i, j, num[i][j], den[i], pi[i][j]);
  }
  printf("------------------------------------------------------------- \n");
}

//===========================================================================
//=  Function to return integer x raised to n                               =
//===========================================================================
int pow_int(int x, int n)
{
  int          pow_val;      // Returned power value
  int          i;            // Loop counter

  // Loop to generate pow_val
  pow_val = 1;
  for (i=0; i<n; i++)
    pow_val = pow_val * x;

  return(pow_val);
}