迭代取中法、乘同余法及混合同余法产生随机数方法

    在用计算机编制程序时，经常需要用到随机数，尤其在仿真等领域，更对随机数的产生提出了较高的要求，仅仅使用 C 语言类库中的随机函数已难以胜任相应的工作。现实中，用投色子计数的方法产生真正的随机数，但电脑若也这样做，将会占用大量内存；虽然用噪声发生器或放射性物质也可产生真正的随机数，但操作不可重复。而用数学方法产生随机数则最适合计算机，这就是周期有限，易重复的“伪随机数”。
（注：这里生成的随机数所处的分布为 0-1 区间上的均匀分布。不是 0-1 区间怎么办？ 除以 (high-low), 再加上 low 就可以完成任务。我们需要的随机数序列应具有非退化性，周期长，相关系数小等优点。）

 

1、迭代取中法：
这里在迭代取中法中介绍平方取中法 , 其迭代式如下 :
                  Xn+1=(Xn^2/10^s)(mod 10^2s) 
                  Rn+1=Xn+1/10^2s
其中，Xn+1是迭代算子，而 Rn+1 则是每次需要产生的随机数。
第一个式子表示的是将 Xn 平方后右移 s 位，并截右端的 2s 位。
而第二个式子则是将截尾后的数字再压缩 2s 倍，显然 :0=迭代取中法有一个显著的不良特性就是它比较容易退化成 0。

 

实现代码（C语言）：

#include

#include

#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;

}

 

前一百个测试生成的随机数序列（明显可见其容易退化为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;

 

 

前三百个测试生成的随机数序列（可见该随机数生成方法所生成的随机序列比较符合0-1上的均匀分布,不过在某些数据段还有些起伏）：

图1：  乘同余法生成的300随机数的产生序列图

 

 

 图2: 乘同余法生成的300随机数的分布情况

 

3、混合同余法：

混合同余法是加同余法和乘同余法的混合形式,其迭代式如下:

                Xn+1=(Lamda*Xn+Miu)%M

                Rn+1=Xn/M

经前人研究表明，在M=2^q的条件下，参数lamda,miu,X0按如下选取，周期较大，概率统计特性好:

                Lamda=2^c+1,c取q/2附近的数

                Miu=(1/2+sqrt(3))/M

                X0为任意非负整数

      

实现代码（C语言）关键部分：

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;

 

前三百个测试生成的随机数序列（可见该种随机数生成方法已相当接近0-1上的均匀分布。但在图3中可以看出它的一个致命的弱点，那就是随机数的生成在某一周期内成线性增长的趋势，显然，在大多数场合，这种极富“规律”型的随机数是不应当使用的。）：

图3: 乘同余法生成的300随机数的产生序列图

图4: 乘同余法生成的300随机数的分布情况