平方取中法是冯·纽曼(John van Neumann)在40年代中期提出的。这个方法首先从某个初始的种子数开始,求出这个数的平方。取这个平方数的中间几位作为随机数序列中的第2个数;再求出第2个数的平方,又取这个平方数的中间几位作为随机数序列中的第3个数;不断按这个方式继续此算法,即可得到相应的伪随机序列。
例1:利用平方取中法产生4位数的随机数序列,序列的种子数取为x0= 3187,通过计算得到
x0= 3187
(3187)2=10156969; x1=1569
(1569)2=02461761; x2=4617
(4617)2=21316689; x3=3166
(3166)2=10023556; x4=0235
(0235)2=00055225; x5=0552
(0552)2=00304704; x6=3047
(3047)2=09284209; x7=2842
将此过程继续下去,还可以得到0769,5913,9635,8332,4222,8542,…,
这个方法有几个缺点。首先,利用这个方法产生的伪随机数序列的重复周期通常较短。第二,对于较长的伪随机数序列,利用这种方法可能无法通过随机性的统计检验。第三,当在任何时候生成之后,其后产生的数都将为0。如果这种现象在一个较复杂的仿真研究过程中出现,它将会使仿真分析人员误入歧途。这种现象如下例所示。
例2:利用平方取中法产生两位数的随机数序列,种子数取为x0=44。通过计算得到:
x0=44
(44)2=1936; x1=93
(93)2=8649; x2=64
(64)2=4096; x3=09
(09)2=0081; x4=08
(08)2=0064; x5=06
(06)2=0036; x6=03
(03)2=0009; x7=00
(00)2=0000; x8=00
这样,随机数就无法继续产生了。
利用平方取中法的另一个问题是这个方法可能产生退化,即总是得到相同的xi值。
例3:设在产生四位数的随机数过程中,得到了一个xi值为4500,即
xi=4500
从而, (4500)2=20250000; xi+1=2500
(2500)2=06250000; xi+2=2500
相继产生的数值总为2500。
由于这些原因,平方取中法已经被许多新的能提供更可靠的随机数序列的算法所取代。
线性同余法在1951年由菜默尔(Lehmer)首先提出。目前大多数随机数发生器都采用这种方法。在这个算法中,随机数序列中的数由如下的递推关系产生
n≥0 (1)
初始值x0称为种子,常数a称为乘子,常数c称为增量,而常数m称为模数。
对于(1)式,当c=0时,该算法称为乘同余法;当c≠0时,该算法称为混合同余法。从例4中可以看出利用这种方法产生的序列的重复性,一般来讲任何由此方法产生的序列都存在重复性。在大多数情况下,合理地选择常数a,c,x0和m,可以使重复周期充分的长。
在(1)式中,显然可以得出
0≤xn ≤m-1 (2)
为了得到[0,1]区间上分布的随机数,可以令
(3)
Rn为满足要求的随机数。
常数m,a,c的值对所产生的随机数序列的周期长度有很大的影响。
例4:设a=5,c=3,m=15,取x0=7,利用线性同余法产生随机数序列。
x0=7
x1=(5*7+3)mod16=6 R1=0.375
x2=(5*6+3)mod16=1 R2=0.063
x3=(5*1)+3mod16=8 R3=0.500
x4=(5*8+3)mod16=11 R4=0.688
x5=(5*11+3)mod16=10 R5=0.625
x6=(5*10+3)mod16=5 R6=0.313
x7=(5*5+3)mod16=12 R7=0.750
x8=(5*12+3)mod16=15 R8=0.938
x9=(5*15+3)mod16=14 R9=0.875
x10=(5*14+3)mod16=9 R10=0.563
x11=(5*9+3)mod16=0 R11=0.000
x12=(5*0+3)mod16=3 R12=0.188
x13=(5*3+3)mod16=2 R13=0.125
x14=(5*2+3)mod16=13 R14=0.813
x15=(5*13+3)mod16=4 R15=0.250
x16=(5*4+3)mod16=7 R16=0.438
x17=(5*7+3)mod16=6 R17=0.375
x18=(5*6+3)mod16=1 R18=0.063
x19=(5*1+3)mod16=8 R19=0.500
在n=16时出现循环。
对于常数的选择讨论如下:
① m的选择
由于重复周期的长度总是小于m,因此需要将m取大的数值,更进一步,所选用的m的值应能简化同余关系的解,在计算机中数字都是用二进制表达的,因此已经证明m取值为(2k)是很好的。
② a和c的选取
当且仅当下列条件满足时,一个由线性同余法产生的随机数序列的最大可能重复周期为m。
l c与m互质,即同时能被c和m整除的正整数只有1。
l 如果m能被4整除,则(a-1)也能被4整除,即a=1+4k。
特别地,当选择a=216+5=65541或a=216+3=65539时可以得到满意的结果。至子c的选择,只要满足c与m互为质数的条件即可。
③ x0的选取
如果随机数序列的周期为m,因为能产生完全的序列,即在一个周期内可以取到0至(m-1)的所有值,因此x0的选取是不重要的。但仍然要小心,例如取x0=0时会产生退化的序列。
对于乘同余法,由于c=O,无论怎样选择m,都无法满足c与m互质的条件,因而不可能得到满周期。若选择m=2k,则所产生的随机数序列的周期p≤2k-2,即在0至m-1之间的整数至多只有四分之一可能成为xn的值,而且这四分之一的整数在0至(m-1)之间是如何分布的尚难确定。这与种子数x0的选取有关。若取乘子为a=8L十3或a=8L+5形式的整数,种子x0取为奇数.则可以达到最长的周期p=2k-2。
例5:使用乘同余法,对a=13,m=2^6=64,且x0=1,2,3,4,求产生器的周期。可以看出当种子为1或3时,序列的周期都为16。而当种子为2或4时,序列的周期分别为8和4。
i |
X0=1 |
X0=2 |
X0=3 |
X0=4 |
0 |
1 |
2 |
3 |
4 |
1 |
13 |
26 |
39 |
52 |
2 |
41 |
18 |
59 |
36 |
3 |
21 |
42 |
63 |
20 |
4 |
17 |
34 |
51 |
4 |
5 |
29 |
58 |
23 |
|
6 |
57 |
50 |
43 |
|
7 |
37 |
10 |
47 |
|
8 |
33 |
2 |
35 |
|
9 |
45 |
|
7 |
|
10 |
9 |
|
27 |
|
11 |
53 |
|
31 |
|
12 |
49 |
|
19 |
|
13 |
61 |
|
55 |
|
14 |
25 |
|
11 |
|
15 |
5 |
|
15 |
|
16 |
1 |
|
3 |
|
从(1)式中得到的序列χ0,χ1…实质上完全不是随机的,因为若设
(4)
则有
=
┉
=
即一旦m,a,c,x0确定,xi就完全被唯一地确定下来了。而且在[0,1]区间上得到的数值Rn最多只能取到m个值。
加同余法需要n个数的序列作为它的种子,这n个数的序列可以应用其它的方法产生,应用加同余法可以使这个序列不断扩大,加同余法的算法是
(5)
这种方法的主要好处是速度快,因为它不需要作乘法运算。它可以得到大于m的周期但这种方法产生随机数的过程不象混合同余法那样清楚,因此由此方法产生的随机数序列需要经过仔细的确认。
二次平方同余法适用于m为2的幂次的情况,这种方法的递推关系式为
n>0 (6)
种子数x0必须满足条件