/* Compute numerical results using Erlang B loss formula */
/* Usage: cc -o erlangloss erlangloss.c -lm               */
/*        ./erlangloss 1.0 3.0 5                          */

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

#define MAX_CAPACITY 100
/* #define DEBUG 1 */


double factorial(int f)
  {
    if ( f == 1 ) {
      return f * 1.0;
    }
    else return (double) (factorial(f-1) * f); 
  }

double summation(float Y, int M)
  {
    if ( M == 0 )
      return 1;
    else return summation(Y, M-1) + pow(Y, M)/factorial(M);
  }

int main (int argc, char * argv[])
  {
    float lambda, dur, rho;
    int C;
    float loss;
    float prob[MAX_CAPACITY];
    float mean, var;
    int i;

    if (argc != 4)
      {
	printf("Usage: erlangloss lambda dur C\n");
	exit(1);
      }
  
    lambda = atof(argv[1]);
    dur = atof(argv[2]);
    C = atoi(argv[3]);
    rho = lambda * dur;

    if (C < 1)
      {
	printf("Circuit capacity too low!!!\n");
	exit(1);
      }

    if (C > MAX_CAPACITY)
      {
	printf("Circuit capacity too high!!!\n");
	exit(1);
      }

#ifdef DEBUG
    /* Compute and report the state probabilities */
    prob[0] = 1.0;
    for( i = 1; i <= C; i++ )
      {
	prob[i] = (pow(rho, i)/factorial(i))/(summation(rho, C));
	prob[0] -= prob[i];
      }
    printf("State probabilities: \n");
    mean = 0.0;
    var = 0.0;
    for( i = 0; i <= C; i++ )
      {
	printf("  p[%d] = %8.6f\n", i, prob[i]);
	mean += i * prob[i];
	var += i * i * prob[i];
      }
    var -= mean * mean;
    printf("Mean occupancy: %8.6f\n", mean);
    printf("Variance of occupancy: %8.6f\n", var);
    printf("\n");
#endif

    loss = (pow(rho, C)/factorial(C))/(summation(rho, C));

    printf("ErlangBLossFormula(%i, %f) = %f\n", C, rho, loss);
    printf("\n");
    printf("EffectiveUtilization(%i, %f) = %f\n", C, rho, (1 - loss)*rho/C);

    return 0;
  }