//===================================================== file = chain2.c =====
//=  Iterative solver for M/M/1/K Markov chain                              =
//=   - Assumes floating point 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) Notation is pi[i][j] where i is size of chain (i+1 states) and j  =
//=       is index for state (0,... i).                                     =
//=    4) Output is to stdout                                               =
//=-------------------------------------------------------------------------=
//= Example execution: (for K = 3, LAMBDA = 1.1 and MU = 2.25               =
//=                                                                         =
//=   -------------------------------------------------------------         =
//=   K = 3   Lambda = 1.200000   Mu = 2.250000                             =
//=   -------------------------------------------------------------         =
//=   pi[ 0][ 0] = 1.000000                                                 =
//=   -------------------------------------------------------------         =
//=   pi[ 1][ 0] = 0.652174                                                 =
//=   pi[ 1][ 1] = 0.347826                                                 =
//=   -------------------------------------------------------------         =
//=   pi[ 2][ 0] = 0.550122                                                 =
//=   pi[ 2][ 1] = 0.293399                                                 =
//=   pi[ 2][ 2] = 0.156479                                                 =
//=   -------------------------------------------------------------         =
//=   pi[ 3][ 0] = 0.507748                                                 =
//=   pi[ 3][ 1] = 0.270799                                                 =
//=   pi[ 3][ 2] = 0.144426                                                 =
//=   pi[ 3][ 3] = 0.077027                                                 =
//=   -------------------------------------------------------------         =
//=-------------------------------------------------------------------------=
//=  Build: bcc32 chain2.c, gcc chain2.c -lm                                =
//=-------------------------------------------------------------------------=
//=  Execute: chain2                                                        =
//=-------------------------------------------------------------------------=
//=  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 (from chain1.c)                      =
//===========================================================================

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

//----- Defines -------------------------------------------------------------
#define         K     3            // Size of chain (states = 0, 1,... K)
#define    LAMBDA  1.20            // Arrival rate
#define        MU  2.25            // Service rate

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

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

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

  // 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(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] = num[i][j] / den[i];

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