//本程序用于M/M/C/C系統的阻塞率仿真計算
#include "stdio.h"
#include "math.h"   
#define znew (z=36969*(z&65535)+(z>>16))
#define wnew (w=18000*(w&65535)+(w>>16))
#define MWC    ((znew<<16)+wnew )
#define SHR3  (jsr^=(jsr<<17), jsr^=(jsr>>13), jsr^=(jsr<<5))
#define CONG  (jcong=69069*jcong+1234567)
#define KISS  ((MWC^CONG)+SHR3)    //用KISS算法生成U(0,1)分布隨機數
#define NUM 50000            //每次仿真到達的總用戶數
typedef unsigned long UL;

static UL z=362436069, w=521288629, jsr=123456789, jcong=380116160;//全局靜態變量

void settable(UL i1,UL i2,UL i3,UL i4)  //數表初始化函數
{ z=i1;
  w=i2,
  jsr=i3;
  jcong=i4; 
}

double expo(float lamda)  //生成參數為lamda的指數分布隨機數
{  UL k;
   double U01;
   k=KISS;
   U01=(double)(k-1372460312U)/4294967296;
   return(-lamda*log(1-U01));
} 

double theory(float arrival_rate,float server_rate,int capacity) //應用Erlang公式計算阻塞率的理論值
{long double a,factor=0,block,c=1,temp=1;
  int i,j;
  a=arrival_rate/server_rate;
  for(i=capacity;i>0;i--)
	c=c*i;
  for(i=1;i<=capacity;i++)
  {
   temp=1;
   for(j=i;j>0;j--)
     temp=temp*j;
   factor=factor+pow(a,i)/temp;
  }
  block=pow(a,capacity)/(c*factor);
  return(block);
}

  void main()                 //主函數
  {
   int i,j,resource,capacity; //resource剩余資源數,capacity系統容量
   float arrival_rate,server_rate,lamda,mu,block_probability,block;
   double interarrival,holdtime,a[500];
   FILE *out,*out2,*out3;    //輸出文件指針
   settable(12345,65435,34221,12345);  //對KISS數表進行初始化
   arrival_rate=3600;                        //到達率(單位時間為小時)
   server_rate=20;                         //服務率
   lamda=1/arrival_rate;                    //到達時間間隔
   mu=1/server_rate;                        //每個用戶的保持時間
   out=fopen("block.txt","w");
   out2=fopen("capacity.txt","w");
   out3=fopen("theory.txt","w");
   for(capacity=1;capacity<=231;capacity+=10) //統計各種不同容量的情況
   {
    resource=capacity;
    for(i=0;i<capacity;i++)
	  a[i]=0;
	block=0;
     for(j=0;j<NUM;j++)
	 {
      interarrival=expo(lamda);      //生成參數為lamda的指數分布隨機數,代表用戶到達時間間隔
      holdtime=expo(mu);           //生成參數為mu的指數分布隨機數，代表該用戶的通話時間
	  for(i=0;i<capacity;i++)        //判斷是否有用戶離開，有則釋放資源
		if(a[i]!=0)
		{
		 a[i]=a[i]-interarrival;
		 if(a[i]<=0)                 //每次有一個用戶離開，則資源數加1
		 {
		  a[i]=0;
		  resource++;
		 }
		}
	 if(resource==0)                 //沒有資源時，阻塞次數加1，等待下一個用戶到達
	 {
	  block++;
	  continue;
	 }
	 for(i=0;i<capacity;i++)         //有資源時，把該用戶的通話時間記錄到數組c[i]
	  if(a[i]==0)
	  {
	   a[i]=holdtime;
	   resource--;                   //資源數減1
	   break;
      }
	}
   fprintf(out2,"%d\n",capacity);    //分別輸出容量與它對應的理論阻塞率和仿真阻塞率
   block_probability=block/NUM;
   fprintf(out,"%f\n",block_probability);
   fprintf(out3,"%f\n",theory(arrival_rate,server_rate,capacity));
   printf("capacity %d success %f %f\n",capacity,block_probability,theory(arrival_rate,server_rate,capacity));
   }
   fclose(out);                    //關閉文件
   fclose(out2);
   fclose(out3);
 }