//===================================================== file = chain3.c =====
//=  Iterative solver for M/M/1/K Markov chain                              =
//=   - Assumes arrival and service rates are in Lambda[] and Mu[] arrays   =
//===========================================================================
//=  Notes:                                                                 =
//=    1) Solves for M/M/1/K steady state probailities using a novel (?)    =
//=       iterative method.  Note the need for normalizing pi.              =
//=    2) Must manually set K value in #define                              =
//=    3) Must manually load K values into Lambda[] and mu[] in globals.    =
//=       Notation is Lambda[i] and Mu[i] are between states i and i + 1.   =
//=    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:                                                      =
//=                                                                         =
//=    -------------------------------------------------------------        =
//=    K = 4                                                                =
//=    Lambda = 2.000000 2.000000 5.000000 2.000000                         =
//=    Mu     = 3.000000 3.000000 4.000000 3.000000                         =
//=    -------------------------------------------------------------        =
//=    pi[ 0][ 0] = 1.000000                                                =
//=    -------------------------------------------------------------        =
//=    pi[ 1][ 0] = 0.600000                                                =
//=    pi[ 1][ 1] = 0.400000                                                =
//=    -------------------------------------------------------------        =
//=    pi[ 2][ 0] = 0.473684                                                =
//=    pi[ 2][ 1] = 0.315789                                                =
//=    pi[ 2][ 2] = 0.210526                                                =
//=    -------------------------------------------------------------        =
//=    pi[ 3][ 0] = 0.375000                                                =
//=    pi[ 3][ 1] = 0.250000                                                =
//=    pi[ 3][ 2] = 0.166667                                                =
//=    pi[ 3][ 3] = 0.208333                                                =
//=    -------------------------------------------------------------        =
//=    pi[ 4][ 0] = 0.329268                                                =
//=    pi[ 4][ 1] = 0.219512                                                =
//=    pi[ 4][ 2] = 0.146341                                                =
//=    pi[ 4][ 3] = 0.182927                                                =
//=    pi[ 4][ 4] = 0.121951                                                =
//=    -------------------------------------------------------------        =
//=-------------------------------------------------------------------------=
//=  Build: bcc32 chain3.c, gcc chain3.c -lm                                =
//=-------------------------------------------------------------------------=
//=  Execute: chain3                                                        =
//=-------------------------------------------------------------------------=
//=  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 chain2.c                       =
//===========================================================================

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

//----- Defines -------------------------------------------------------------
#define         K     4            // Size of chain (states = 0, 1,... K)

//----- Globals -------------------------------------------------------------
double       Lambda[K] =           // Lambda values (K values)
{2.0, 2.0, 5.0, 2.0};
double       Mu[K] =               // Mu values (K values)
{3.0, 3.0, 4.0, 3.0};

//===========================================================================
//=  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
  double       sum;              // Sum variable for normalizing pi[][]
  int          i, j;             // Loop counters

  // Output initialized values
  printf("------------------------------------------------------------- \n");
  printf("K = %d \n", K);
  printf("Lambda = ");
  for (i=0; i<K; i++)
    printf("%f ", Lambda[i]);
  printf("\n");
  printf("Mu     = ");
  for (i=0; i<K; i++)
    printf("%f ", Mu[i]);
  printf("\n");

  // 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[i-1] * den[i-1]) + pow(Lambda[i-1], i);

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

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

  // Compute all pi values (un-normalized and then normalize)
  for (i=0; i<=K; i++)
  {
    sum = 0.0;
    for (j=0; j<=i; j++)
    {
      pi[i][j] = num[i][j] / den[i];
      sum = sum + pi[i][j];
    }

    for (j=0; j<=i; j++)
      pi[i][j] = pi[i][j] / sum;
  }

  // Output all results (normalized 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");
}