常見偽隨機數常生

#include "stdlib.h"
#include "stdio.h"
#include "math.h"

double uniform(double a,double b,long int* seed);
double gauss(double mean,double sigma,long int *seed);
double exponent(double beta,long int *seed);
double laplace(double beta,long int* seed);
double rayleigh(double sigma,long int *seed);
double weibull(double a,double b,long int*seed);
int bn(double p,long int*seed);
int bin(int n,double p,long int*seed);
int poisson(double lambda,long int *seed);

void main()
{
double a,b,x,mean;
int i,j;
long int s;

a=4;
b=0.7;
s=13579;
mean=0;
for(i=0;i<10;i++)
{
for(j=0;j<5;j++)
{
x=poisson(a,&s);
mean+=x;
printf("%-13.7f",x);
}
printf("/n");
}
mean/=50;
printf("平均值為:%-13.7f/n",mean);
}

/*******************************************************************
* 求[a,b]上的均勻分布
* 輸入: a--雙精確度實型變數，給出區間的下限
*       b--雙精確度實型變數，給出區間的上限
*    seed--長整型指標變數，*seed為隨機數的種子  
********************************************************************/
double uniform(double a,double b,long int*seed)
{
double t;
*seed=2045*(*seed)+1;
*seed=*seed-(*seed/1048576)*1048576;
t=(*seed)/1048576.0;
t=a+(b-a)*t;

return(t);
}

/*******************************************************************
* 常態分佈
* 輸入: mean--雙精確度實型變數，常態分佈的均值
*      sigma--雙精確度實型變數，常態分佈的均方差
*       seed--長整型指標變數，*seed為隨機數的種子  
********************************************************************/
double gauss(double mean,double sigma,long int*seed)
{
int i;
double x,y;
for(x=0,i=0;i<12;i++)
x+=uniform(0.0,1.0,seed);
x=x-6.0;
y=mean+x*sigma;

return(y);
}

/*******************************************************************
* 指數分布
* 輸入: beta--指數分布均值
*       seed--種子
*******************************************************************/
double exponent(double beta,long int *seed)
{
double u,x;
u=uniform(0.0,1.0,seed);
x=-beta*log(u);

return(x);
}

/*******************************************************************
* 拉普拉斯隨機分布
* beta--拉普拉斯分布的參數
* *seed--隨機數種子
*******************************************************************/
double laplace(double beta,long int* seed)
{
double u1,u2,x;

u1=uniform(0.,1.,seed);
u2=uniform(0.,1.,seed);
if(u1<=0.5)
x=-beta*log(1.-u2);
else
x=beta*log(u2);

return(x);
}

/********************************************************************
* 瑞利分布
*
********************************************************************/
double rayleigh(double sigma,long int *seed)
{
double u,x;
u=uniform(0.,1.,seed);
x=-2.0*log(u);
x=sigma*sqrt(x);
return(x);
}

/************************************************************************/
/* 韋伯分布                                                                     */
/************************************************************************/
double weibull(double a,double b,long int*seed)
{
double u,x;

u=uniform(0.0,1.0,seed);
u=-log(u);
x=b*pow(u,1.0/a);

return(x);
}

/************************************************************************/
/* 貝努利分布                                                           */
/************************************************************************/
int bn(double p,long int*seed)
{
int x;
double u;
u=uniform(0.0,1.0,seed);
x=(u<=p)?1:0;
return(x);
}

/************************************************************************/
/* 二項式分布                                                           */
/************************************************************************/
int bin(int n,double p,long int*seed)
{
int i,x;
for(x=0,i=0;i<n;i++)
x+=bn(p,seed);
return(x);
}

/************************************************************************/
/* 泊松分布                                                             */
/************************************************************************/
int poisson(double lambda,long int *seed)
{
int i,x;
double a,b,u;

a=exp(-lambda);
i=0;
b=1.0;
do {
u=uniform(0.0,1.0,seed);
b*=u;
i++;
} while(b>=a);
x=i-1;
return(x);
}

 

/************************************************************************
* 離散傅立葉變換與反變換
* 輸入: x--要變換的資料的實部
* y--要變換的資料的虛部
*       a--變換結果的實部
*       b--變換結果的虛部
*       n--資料長度
*    sign--sign=1時，計算離散傅立葉正變換;sign=-1時;計算離散傅立葉反變換
************************************************************************/
void dft(double x[],double y[],double a[],double b[],int n,int sign)
{
int i,k;
double c,d,q,w,s;

q=6.28318530718/n;
for(k=0;k<n;k++)
{
w=k*q;
a[k]=b[k]=0.0;
for(i=0;i<n;i++)
{
d=i*w;
c=cos(d);
s=sin(d)*sign;
a[k]+=c*x+s*y;
b[k]+=c*y-s*x;
}
}
if(sign==-1)
{
c=1.0/n;
for(k=0;k<n;k++)
{
a[k]=c*a[k];
b[k]=c*b[k];
}
}
}