(注:这里生成的随机数所处的分布为 0-1 区间上的均匀分布。不是 0-1 区间怎么办? 除以 (high-low), 再加上 low 就可以完成任务。我们需要的随机数序列应具有非退化性,周期长,相关系数小等优点。)
这里在迭代取中法中介绍平方取中法 , 其迭代式如下 :
Xn+1=(Xn^2/10^s)(mod 10^2s)
Rn+1=Xn+1/10^2s
其中,Xn+1是迭代算子,而 Rn+1 则是每次需要产生的随机数。
第一个式子表示的是将 Xn 平方后右移 s 位,并截右端的 2s 位。
而第二个式子则是将截尾后的数字再压缩 2s 倍,显然 :0=<Rn+1<=1。
迭代取中法有一个显著的不良特性就是它比较容易退化成 0。
- #include <stdio.h>
- #include <math.h>
- #define S 2
- float Xn=12345; //Seed & Iter
- float Rn; //Return Val
- void InitSeed(float inX0)
- {
- Xn=inX0;
- }
- float MyRnd()
- {
- //implementation: Xn+1=(Xn^2/10^S)(mod 10^2S)
- Xn=(int)fmod((Xn*Xn/pow(10,S)),pow(10,2*S));//here can's use %
- //implementation: Rn+1=Xn+1/10^2S
- Rn=Xn/pow(10,2*S);
- return Rn;
- }
#include <stdio.h>
#include <math.h>
#define S 2
float Xn=12345; //Seed & Iter
float Rn; //Return Val
void InitSeed(float inX0)
{
Xn=inX0;
}
float MyRnd()
{
//implementation: Xn+1=(Xn^2/10^S)(mod 10^2S)
Xn=(int)fmod((Xn*Xn/pow(10,S)),pow(10,2*S));//here can's use %
//implementation: Rn+1=Xn+1/10^2S
Rn=Xn/pow(10,2*S);
return Rn;
}
- int main()
- {
- int i;
- FILE * debugFile;
- if((debugFile=fopen("outputData.txt","w"))==NULL)
- {
- fprintf(stderr,"open file error!");
- return -1;
- }
- printf("\n");
- for(i=0;i<100;i++)
- {
- tempRnd=MyRnd();
- fprintf(stdout,"%f ",tempRnd);
- fprintf(debugFile,"%f ",tempRnd);
- }
- getchar();
- return 0;
- }
int main()
{
int i;
FILE * debugFile;
if((debugFile=fopen("outputData.txt","w"))==NULL)
{
fprintf(stderr,"open file error!");
return -1;
}
printf("\n");
for(i=0;i<100;i++)
{
tempRnd=MyRnd();
fprintf(stdout,"%f ",tempRnd);
fprintf(debugFile,"%f ",tempRnd);
}
getchar();
return 0;
}
前一百个测试生成的随机数序列(明显可见其容易退化为0):
0.399000 0.920100 0.658400 0.349000 0.180100 0.243600 0.934000 0.235600 0.550700 0.327000 0.692900 0.011000 0.012100 0.014600 0.021300 0.045300 0.205200 0.210700 0.439400 0.307200 0.437100 0.105600 0.115100 0.324800 0.549500 0.195000 0.802500 0.400600 0.048000 0.230400 0.308400 0.511000 0.112100 0.256600 0.584300 0.140600 0.976800 0.413800 0.123000 0.512900 0.306600 0.400300 0.024000 0.057600 0.331700 0.002400 0.000500 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
2、乘同余法:
乘同余法的迭代式如下:
Xn+1=Lamda*Xn(mod M)
Rn+1=Xn/M
当然,这里的参数选取是有一定理论基础的,否则所产生的随机数的周期将较小,相关性会较大。经过前人检验的两组性能较好的素数取模乘同余法迭代式的系数为:
Lamda=5^5,M=2^35-31
Lamda=7^5,M=2^31-1
实现代码(C语言)关键部分:
- double long M;//请注意,这里一定要用到double long,否则计算2^32会溢出
- float MyRnd()
- {
- Xn=fmod(Lamda*Xn,M);//here can's use %
- Rn=Xn/M;
- return Rn;
- }
- 另外初始化段应有:
- Lamda=pow(5,5);
- M=pow(2,35)-31;
double long M;//请注意,这里一定要用到double long,否则计算2^32会溢出
float MyRnd()
{
Xn=fmod(Lamda*Xn,M);//here can's use %
Rn=Xn/M;
return Rn;
}
另外初始化段应有:
Lamda=pow(5,5);
M=pow(2,35)-31;
前三百个测试生成的随机数序列(可见该随机数生成方法所生成的随机序列比较符合0-1上的均匀分布,不过在某些数据段还有些起伏):
图1: 乘同余法生成的300随机数的产生序列图
![](http://blogimg.chinaunix.net/blog/upfile2/071201174936.jpg)
- float MyRnd()
- {
- Xn=fmod(Lamda*Xn+Miu,M);
- Rn=Xn/M;
- return Rn;
- }
- 另外初始化段应有:
- M=pow(2,32);
- Lamda=pow(2,16)+1;
- Miu=(0.5+sqrt(3)/6)/M;
float MyRnd()
{
Xn=fmod(Lamda*Xn+Miu,M);
Rn=Xn/M;
return Rn;
}
另外初始化段应有:
M=pow(2,32);
Lamda=pow(2,16)+1;
Miu=(0.5+sqrt(3)/6)/M;
![](http://blogimg.chinaunix.net/blog/upfile2/071201175547.jpg)
![](http://blogimg.chinaunix.net/blog/upfile2/071201175635.jpg)