simulib.c

一个满足各种分布的随机数产生方法,包含了指数分布、几何分布、ELANG分布等等

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <conio.h>
#include <malloc.h>
#include <fstream.h>
#include "SIMULIB.h"

extern struct list_point processing_queue,free_event_list;

void sim_init()
{
	time_t ltime;
	int i;

    time(&ltime);
	srand((unsigned)ltime);
	for(i=0;i<1000;i++) rand();
	
	processing_queue.head=processing_queue.tail=NULL;
	processing_queue.length=0;
    free_event_list.head=free_event_list.tail=NULL;
	free_event_list.length=0;
}

//将当前事件表current_event放回空表队列
void release_event_list(struct event *current_event){
	current_event->link=free_event_list.head;
	free_event_list.head=current_event;
	free_event_list.length++;
	if(free_event_list.tail==NULL) free_event_list.tail=current_event;
}

//得到一个新的事件表
struct event* get_new_event(){
	int i;
	struct event *tem_point;
//	printf("get_new_event\n");
	if(free_event_list.head==NULL){
		tem_point=(struct event *)malloc(1000*sizeof(struct event));
		if( tem_point == NULL ){
            printf( "Insufficient memory available length=%d\n",free_event_list.length);
			exit(1);
		}
		free_event_list.head=tem_point;
		for(i=0;i<1000;i++) {
			tem_point[i].link=&tem_point[i+1]; //+sizeof(struct event);
            //tem_point=tem_point->link;
		}
		tem_point[999].link=NULL;
		free_event_list.tail=&tem_point[999];
		free_event_list.length=1000;
    }
	tem_point=free_event_list.head;
	free_event_list.head=free_event_list.head->link;
	free_event_list.length--;
//	printf("get_new_event complete\n");
    return tem_point; 
}

//将事件current_event按时间序插入事件处理队列processing_queue
void insert_processing_list(struct event *current_event){
	struct event *tem_point_1,*tem_point_2;
	double wait_time;
	processing_queue.length++;
	if(processing_queue.head==NULL){
		processing_queue.head=current_event;
		processing_queue.tail=current_event;
		current_event->link=NULL;
		return;
	}

	tem_point_1=processing_queue.head;
	tem_point_2=processing_queue.head;
	wait_time=0;
	while((wait_time+tem_point_2->timer)<current_event->timer){
		tem_point_1=tem_point_2;
		tem_point_2=tem_point_2->link;
		wait_time+=tem_point_1->timer;
		if(tem_point_2==NULL){    //插到队尾
			current_event->timer-=wait_time;
	        current_event->link=NULL;
	        tem_point_1->link=current_event;
			processing_queue.tail=current_event;
			return;
		}
	}
	current_event->timer-=wait_time;
	tem_point_2->timer-=current_event->timer;
	current_event->link=tem_point_2;
	if(tem_point_1==tem_point_2) processing_queue.head=current_event; //队列长度为一
		else tem_point_1->link=current_event;
	return;
}

//将事件current_event插到队列list尾
void insert_list(struct list_point *list,struct event *current_event){
	if(list->head==NULL) list->head=current_event;
	else (list->tail)->link=current_event;
	list->tail=current_event;
	current_event->link=NULL;
	list->length++;           //队列长度加一
}

//将 list 队列的第一个事件弹出队列， 返回事件地址
struct event *remove_event(struct list_point *list){
	struct event *tem_point;
	tem_point=list->head;
	if(tem_point==NULL){
		printf("Error, remove operation on a NULL list\n");
		exit(1);
	}
	list->head=(list->head)->link;
	list->length--;
	return(tem_point);
}
//清除list 队列
void empty_list(struct list_point *list){
	free_event_list.length+=list->length;
	if(list->head!=NULL){
		(list->tail)->link=free_event_list.head;
		free_event_list.head=list->head;
		if(free_event_list.tail==NULL) free_event_list.tail=list->tail;
		list->head=NULL;
		list->tail=NULL;
		list->average_length=0;
		list->length=0;
	}
	return;
}

void empty_processing_list(){
	free_event_list.length+=processing_queue.length;
	if(processing_queue.head!=NULL){
		(processing_queue.tail)->link=free_event_list.head;
		free_event_list.head=processing_queue.head;
		if(free_event_list.tail==NULL) free_event_list.tail=processing_queue.tail;
		processing_queue.head=NULL;
		processing_queue.tail=NULL;
		processing_queue.average_length=0;
		processing_queue.length=0;
	}
	return;
}

//产生一个指数分布的随机数
double get_exp_random(double lambuda){   
	double x=1.0;
	while(x>=1.0) x=(double)rand()/ (double)RAND_MAX;  //产生一个（0，1）均匀分布的随机数
    return(log(1.0/(1-x))/lambuda);  //返回指数分布的随机数
}

//产生一个几何分布的随机数
int get_geometric_random(double average)
{
	return (int)get_exp_random(log((1.0+average)/average));
}

double ranf()
{
	return (double)rand()/ (double)RAND_MAX;
}
/*--------------------  generate random number in（1，n） --------------------*/
int random_integer (int i, int n)
    {                       
	if (i>n) {i^=n; n^=i; i^=n;}; //exchange i and n
      n-=i;
      n=(int)(((double)n+1.0)*ranf());
      return(i+n);
    }



/*----------------  产生Erlang分布  -----------------*/
double erlang (double x, double s)
    { /* 'erlang' returns a psuedo-random variate from an erlang      */
      /* distribution with mean x and standard deviation s.           */
      int i,k; double z;
      /* if (s>x) then error(0,"erlang Argument Error: s > x"); */
      z=x/s;
      k = (int)(z*z);
      z=1.0; for (i=0; i<k; i++) z *= ranf();
      return(-(x/k)*log(z));
    }

/*-----------  产生超指数分布  -----------*/
double hyperx (double x, double s)
    { /* 'hyperx' returns a psuedo-random variate from Morse's two-   */
      /* stage hyperexponential distribution with mean x and standard */
      /* deviation s, s>x.  */
      double cv,z,p;
      if (s<=x) {printf("s is less than or equal to x");exit(0);}
      cv=s/x; z=cv*cv; p=0.5*(1.0-sqrt((z-1.0)/(z+1.0)));
      z=(ranf()>p)? (x/(1.0-p)):(x/p);
      return(-0.5*z*log(ranf()));
    }

/*-----------  三角随机变量产生器  -----------*/
double triang (double a, double c, double b)
{
 /* triangular distribution with left and right being [a,b] and the
    mode being at point c */
  double sample,point;
  point = (c-a)/(b-a);
  sample = ranf();
  if (sample <= point)
    return(sqrt(sample*(b-a)*(c-a)) + a);
  else
    return(b - sqrt((1.0-sample)*(b-a)*(b-c)));
}

/*-----------------  正态变量产生器  ----------------*/
double normal (double x, double s)
    { /* 'normal' returns a psuedo-random variate from a normal dis-  */
      /* tribution with mean x and standard deviation s.              */
      double v1,v2,w,z1; static double z2=0.0;
      if (z2!=0.0)
		{z1=z2; z2=0.0;}  /* use value from previous call */
        else
          {
            do
              {v1=2.0*ranf()-1.0; v2=2.0*ranf()-1.0; w=v1*v1+v2*v2;}
            while (w>=1.0);
            w=sqrt((-2.0*log(w))/w); z1=v1*w; z2=v2*w;
          }
      return(x+z1*s);
  }
/*-----------------  Pareto分布变量产生器  ----------------*/
//由y=F(x)=1-k^alfa/x^alfa，得x=k*(1-y)^(-1/alfa)
double Pareto(double alfa, double k)
{
	double y=1.0;
	while(y>=1.0) y=(double)rand()/ (double)RAND_MAX;  //产生一个（0，1）均匀分布的随机数
	return k*pow(1-y,-1.0/alfa);
}

void main(void)
{
//Program below are FDM simulation
	int n;//n is the station number;
	double c;//c is the channel rate;
	double u;//u is the lengh of the frame;
	struct event * Arrive_E,*Depart_E,*Current_E; 
    l=3;
    c=600;//kbps
    u=0.6;//us
    double a[500]/*到达间隔*/，t[500]/*到达时刻*/,d[500]/*排队时间*/,s[500]/*服务时间*/,c[500]/*离去时刻*/,tmp;
    int Q,Sev;
	int n;//l离去顾客计数器
  sim_init();
  for (i=1;i<500;i++)
  {
	  if (i==1)
	  {
              a[i]=get_exp_random(c/n);
	       t[i]=a[i];
		d[i]=0;
		s[i]=get_exp_random(u*c/n);
		c[i]=s[i];
		Q=0;Sev=1；
	       Arrive_E->timer=t[1];
		Arrive_E->type_of_event=0;
              insert_processing_list(Arrive_E);
              Depart_E->timer=c[i];
		Depart_E->type_of_event=1;
		insert_processing_list(insert_processing_list(Arrive_E););
		return;
	  }
      else
	  { 
              a[i]=get_exp_random(c/n);
              t[i]=t[i-1]+a[i];
              Arrive_E->type_of_event=0;
		insert_processing_list(Arrive_E);
		////////////////////////////
               此处为查找current时间
                 Current_E=processing_queue.head;
               for(j=1;j<=2*i+1;j++)
                       Current_E=Current_E->link;
		////////////////////////////
		if (Current_E->type_of_event ==0)
		{
			if (Sev==1)
			{	
				Q=Q+1;
			    for (j=1;j<Q+1;j++)
				{
					c[i]=c[i]+get_exp_random(u*c/l);
				}
				c[i]=c[i]+c[n+1]-t[i];
			}
			else
			{
			     Q=0;
				 c[i]=ti[i]+get_exp_random(c/l);
				 Depart_E->timer=c[i];
				 Depart_E->type_of_event=1;
				 insert_processing_list(Depart_E);
			}

			    
		}
		else 
		{
		        	 if(Q==0)
		        	 	{
		        	 	   Sev=0;n=n+1;
				          c[i]=t[i]+get_exp_random(c/l);
		        	 	}
				else
					{
					    Q=Q-1;
					    Sev=1;
					    n=n+1;
					    if (Q==0) c[i]=t[i]+get_exp_random(c/l);
					    for (j=1;j<Q+1;j++)
				               {
					            c[i]=c[i]+get_exp_random(c/l);
				                }
				           c[i]=c[i]+c[n+1]-t[i];//