#include <stdio.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 FIB ((b=a+b),(a=b-a))
#define KISS ((MWC^CONG)+SHR3)
#define LFIB4 (c++,t[c]=t[c]+t[UC(c+58)]+t[UC(c+119)]+t[UC(c+178)])
#define SWB (c++,bro=(x<y),t[c]=(x=t[UC(c+34)])-(y=t[UC(c+19)]+bro))
#define UNI (KISS*2.328306e-10)
#define VNI ((long)KISS)*4.656613e-10
#define UC (unsigned char)/* a cast operation*/
typedef unsigned long UL;

/*Global static variables:*/
static UL z=362436069,w=521288629,jsr=123456789,jcong=380116160;
static UL a=224466889,b=7584631,t[256];

/*Use random seeds to reset a,w,jsr,jcong,a,b,and the table t[256]*/
static UL x=0,y=0,bro; static unsigned char c=0;

/*Example procedure to set the table,using KISS:*/
void settable(UL i1,UL i2,UL i3,UL i4,UL i5,UL i6)
{int i;a=i1;w=i2;jsr=i3;jcong=i4;a=i5;b=i6;
for(i=0;i<256;i=i+1) t[i]=KISS;}

#define C 12  //The number of channels
#define MIN_CALLING_TIME_LIMIT 1
#define MIU 4
#define CIRCLE_NUM_LAMDA 500
#define CIRCLE_NUM_BLOCKING 400
#include <math.h>
int poisson(double lamda) //generate poisson distribution variable
{int i,x_num;
 double a,b,u_num;
 a=exp(-lamda);
 b=1;
 i=0;
 do
 {i=i+1;
  u_num=UNI;
  b=b*u_num;
 }while(b>=a);
 x_num=i-1;
 return x_num;
}
double exponential(double u)//generate exponential distribution variable
{double x_num,y_num;
 x_num=UNI;
 y_num=-u*log(1-x_num);
 return y_num;
}

void empty_all_channel(double channel_capacity[C])
{int i;
 for(i=0;i<C;i++)
	 channel_capacity[i]=0;
}

int get_empty_channel_num(double channel_capacity[C])
{int i,empty_channel_num;
 empty_channel_num=0;
 for(i=0;i<C;i++)
 {if(channel_capacity[i]==0)
    empty_channel_num=empty_channel_num+1;
 }
 return empty_channel_num;
}

int get_blocking_num(int empty_channel_num,int customer_arrival_num)
{int blocking_num;
 blocking_num=customer_arrival_num-empty_channel_num;
 return blocking_num;
}

void empty_customer_calling_time_record(double customer_calling_time[C])
{int i;
 for(i=0;i<C;i++)
	 customer_calling_time[i]=0;
}

void channel_allocate(int customer_arrival_num,double channel_capacity[C])
{int i,empty_channel_num;
 empty_channel_num=get_empty_channel_num(channel_capacity);
 if(customer_arrival_num<=empty_channel_num)
 {for(i=0;i<customer_arrival_num;i++)
   channel_capacity[i]=1;
 }
 for(i=0;i<C;i++)
 {if(channel_capacity[i]==0)
   channel_capacity[i]=1;
 }
}

int get_allocated_customer_num(int customer_arrival_num,int blocking_num,double channel_capacity[C])
{int empty_channel_num,allocated_customer_num;
 empty_channel_num=get_empty_channel_num(channel_capacity);
 if(blocking_num==0)
	 allocated_customer_num=customer_arrival_num;
 else allocated_customer_num=empty_channel_num;
 return allocated_customer_num;
}

void customer_calling_time(double u,int allocated_customer_num,double customer_calling_time[C])
{int i;
 for(i=0;i<allocated_customer_num;i++)
	 customer_calling_time[i]=exponential(u);
}

void channel_return(double customer_calling_time[C],double min_calling_time_limit,double channel_capacity[C])
{int i;
 double remain_calling_time[C];
 for(i=0;i<C;i++)
 {remain_calling_time[i]=customer_calling_time[i]-min_calling_time_limit;
  if(remain_calling_time[i]<=0)
	  channel_capacity[i]=0;
  else customer_calling_time[i]=remain_calling_time[i];
 }
}

int get_total_blocking_num(double lamda,double u,double channel_capacity[C])
{int customer_arrival_num,empty_channel_num,blocking_num,allocated_customer_num;
 customer_arrival_num=poisson(lamda);
 empty_channel_num=get_empty_channel_num(channel_capacity);
 if(empty_channel_num==0)
	 blocking_num=get_blocking_num(empty_channel_num,customer_arrival_num);
 channel_allocate(customer_arrival_num,channel_capacity);
 if(customer_arrival_num>empty_channel_num)
	 blocking_num=get_blocking_num(empty_channel_num,customer_arrival_num);
 else blocking_num=0;
 allocated_customer_num=get_allocated_customer_num(customer_arrival_num,blocking_num,channel_capacity);
 customer_calling_time(u,allocated_customer_num,channel_capacity);
 channel_return(channel_capacity,MIN_CALLING_TIME_LIMIT,channel_capacity);
 return blocking_num;
}

double get_blocking_probability(double lamda,double u,double channel_capacity[C],int circle_num)
{int customer_total_num,blocking_total_num,i;
 double blocking_probability;
 customer_total_num=0;
 blocking_total_num=0;
 for(i=0;i<circle_num;i++)
 {customer_total_num=customer_total_num+poisson(lamda);
  blocking_total_num=blocking_total_num+get_total_blocking_num(lamda,u,channel_capacity);
 }
 blocking_probability=(double)blocking_total_num/(double)customer_total_num;
 return blocking_probability;
}

void main()
{double lamda,blocking_probability,channel_capacity[C];
 FILE *fp;
 fp=fopen("blocking.txt","w+");
 empty_all_channel(channel_capacity);
 for(lamda=1;lamda<=CIRCLE_NUM_LAMDA;lamda++)
 {blocking_probability=get_blocking_probability(lamda,MIU,channel_capacity,CIRCLE_NUM_BLOCKING);
  fprintf(fp,"%f\n",blocking_probability);
 }
 fclose(fp);
}