无周期伪随机数.doc

无周期伪随机数生成方法素数的应用 郭 凤 鸣(中国地质大学.武汉，430074)摘要 混合同余法生成的伪随机数序列最大周期为M，改进的混合同余法生成的伪随机数序列最大周期扩大到(M-1)M。 无周期伪随机数生成方法解决了伪随机数生成的周期问题。关键词 随机数，混合同余法，改进的混合同余法，素数分类号: O29 A Generating method of cycle-free pseudorandom numbersapplication of prime numberGuo Fengming(China University of Geosciences,Wuhan 430074)Abstract The maximum cycle of pseudorandom number sequence generated by mixed congruence method is M，The maximum cycle of pseudorandom number sequence generated by improved mixed congruence method is (M-1)M，the generating method of cycle-free pseudorandom numbers resolves the cycle problem of pseudo random number generator。Key words random numbers, Mixed congruence method, Improved mixed congruence method, Prime numbers.一 引言Lehmer于1951年提出的混合同余法，具有速度快、内存省、周期长、统计特性好等优点。在各个领域得到广泛应用，几乎所有的计算机语言生成随机数的库函数都选用这一方法。过去的几十里，对于混合同余法计算公式 x i+1= (axi+b) modM (1)ri+1=xi+1/M （2）人们进行了很多研究，旨在寻找较好的参数x0、a、b和M，以生成较长周期的均匀伪随机数序列。到目前为止，由Hull和Dobell给出的证明是最好的。当且仅当x0、a、b和M满足下列条件时，可以得到周期为模数M的伪随机数序列。1. b和M互质；2. 设q为某一质数，取M分别能被q和4整除。总之，对于给定的M，无论如何改变x0、a、b的值，生成的伪随机数序列最大周期不会超过模数M。在实际应用中，为了得到数量较多的伪随机数不得不选用较大的M。例如在字长为32位的计算机系统中取M为-1。我们有个想法：以前的研究，在生成伪随机数序列时，x0、a、b和M都是基于固定不变的常数。假若把x0、a、b和M看作参变量，在伪随机数生成过程中，让它们可以变化，能否生成较好的伪随机数序列呢？按照这一想法，对x0、a、b和M的不同取值进行了大量计算、比较和统计分析，1992年我们提出了改进的混合同余法。其基本思想是：对于(1)式取M为素数，b为变量，b=i=1,2,3,，经过适当选择a，使得生成伪随机数序列的周期可以达到(M-1)M。下面给出一个实例，取M=7，a=3，x0=5，b=1,2,3,50，计算50个伪随机数如表一。表一 50个伪随机数表 5 2 1 6 1 1 2 6 5 35 5 6 3 2 0 2 2 3 06 4 6 6 0 4 3 1 3 34 1 0 5 0 0 1 5 4 24 4 5 2 1 6 1 1 2 6从表一中可以看出，其中前42个数无周期性，第43个数之后，开始周期重复。用这一方法得到的伪随机数，经检验满足均匀性和独立性要求。新的计算数据，引进素数作为增量b,又有了新的进展，使得生成无周期伪随机数成为可能，介绍于此。二 无周期伪随机数生成方法1用素数作为增量将(1)式中b的取值变为b= x1，x2，x3 =2,3,5,7,11,13 ，其中xi为素数。这时生成的伪随机数序列周期会更长。表二是取M=7，a=3，x0=3，b=3,5,7,11,13,17,19,时，计算得到的350个伪随机数。表中伪随机数的数量350已经大于M，没有看到重复周期。我们计算到3000MM个数据，也没有看到重复周期，以此推测周期达到无限大。进一步计算证明这种推测是正确的。表二 350个伪随机数表(M=7,a=3)3 5 6 4 2 5 4 3 4 5 0 2 5 2 4 2 3 0 4 6 0 2 4 6 3 5 6 3 6 5 2 4 2 5 3 6 0 2 5 6 1 2 1 0 1 6 5 0 3 0 2 4 3 1 0 3 5 6 1 4 1 3 5 4 3 4 0 1 0 6 0 5 3 4 6 2 3 0 2 3 1 4 2 5 6 3 4 0 6 5 6 5 4 6 3 1 2 3 0 2 0 4 1 5 6 0 6 2 3 1 2 3 4 1 4 2 5 4 0 1 6 5 6 1 2 0 2 4 5 6 1 4 0 1 2 5 4 1 2 3 1 5 6 5 4 5 0 3 0 2 1 2 0 5 6 5 3 0 6 0 2 0 1 2 3 5 6 0 1 2 3 1 5 6 0 6 5 6 3 0 2 5 2 4 2 3 6 0 2 3 0 1 5 6 1 6 1 5 0 5 6 3 0 3 2 6 2 1 4 1 5 0 1 2 3 0 4 1 5 2 1 2 1 5 0 1 2 3 0 2 3 5 2 3 1 6 2 1 4 5 6 2 3 4 5 6 3 5 6 5 2 4 2 5 4 0 4 1 2 4 2 0 3 1 6 2 0 6 5 6 2 3 5 3 1 2 1 0 2 4 1 0 6 3 0 2 3 4 1 2 3 5 3 1 4 1 2 5 6 30 2 3 4 0 1 0 6 3 4 3 6 3 6 5 4 5 6 3 41 5 0 1 0 5 3 0 2 0 2 3 1 6 1 2 0 5 0 50 5 0 1 5 2 3 4 5 6三 无周期伪随机数生成方法2用奇合数作为增量将(1)式中b的取值变为b= x1，x2，x3，其中xi为奇合数，这时生成的伪随机数序列周期也是无限的。表三是取M=11，a=3，x0=3，b=9,15,21,25,27,33,35,，计算得到的350个伪随机数。四 无周期伪随机数生成方法3用合数作为增量将(1)式中b的取值变为b= x1，x2，x3，其中xi为合数，也可以生成无周期伪随机数序列。例如取M=11，a=3，x0=3，b=4,6,8,9,10,12,14,15,。计算结果不再展示。三种无周期伪随机数生成方法都是素数的应用，因为要得到合数必须从整数中剔除素数。除了这三种方法以外，还可以构造其他数列，用以替换 (1)式中的增量b，都可以生成无周期伪随机数序列。根据以上思想，生成伪随机数的方法可以有很多，也很容易。在做随机试验和模拟时，可以自己构造一些随机数序列，选择更多的试验方案，也节约时间，实现起来并不复杂。例如，需要数字在100以内的随机数序，取M=101，可以轻易生成10000个，甚至更多伪随机数都不困难。表三 350个伪随机数表(M=11,a=3)3 7 3 8 5 9 5 6 2 7 4 8 2 8 10 7 2 4 1 78 1 6 1 10 8 9 6 1 10 6 7 1 7 8 3 1 2 6 910 7 0 5 0 1 7 5 3 10 4 10 8 4 9 8 7 6 5 45 10 5 3 10 2 6 0 6 8 5 9 1 3 0 6 4 2 9 39 10 7 0 3 4 10 8 4 5 10 5 6 3 9 7 3 4 9 45 9 3 9 7 3 8 5 9 3 9 7 5 1 4 5 0 1 5 105 3 1 10 6 7 3 4 9 4 6 3 7 10 3 6 8 7 6 71 7 8 5 0 9 5 6 0 8 1 4 6 5 4 3 4 0 1 89 3 9 7 8 1 4 5 9 1 2 6 0 6 7 0 3 4 10 86 2 7 2 3 7 10 0 6 4 2 9 3 9 10 3 8 7 6 54 3 4 9 4 2 3 9 7 5 1 6 1 10 6 7 3 6 7 20 7 10 0 4 7 9 6 10 4 10 8 4 7 8 1 6 1 2 60 8 1 4 6 3 9 10 3 6 8 7 6 5 4 3 2 3 8 31 10 6 7 1 7 9 6 10 2 6 2 3 8 3 4 8 0 1 58 10 0 7 8 2 8 6 2 5 9 3 9 7 3 4 9 4 2 91 3 0 4 9 4 2 9 1 2 6 0 8 1 6 1 2 6 9 103 6 8 7 6 7 1 9 2 5 7 4 10 8 4 5 10 7 2 01 5 8 9 6 10 4 10 8 4五 百万电子随机数表百万电子随机数表是根据无周期伪随机数生成方法1编写的一个Foxpro程序。程序短小，使用方便。运行此程序，生成随机数的数量可以达到百万个，不会出现周期性。随机数的位数可以是1、2、3、4位，任你选择。程序简略框图1计算09之间的随机数请选择随机数类型1 092 0992 09994 099995 结束2计算099之间的随机数初始化计算0999之间的随机数34计算09999之间的随机数5结 束 程序清单* 增量为素数的混合同余法计算伪随机数素数应用set defa to c:program filesmicrosoft visual studiovfp98同余法select 1use ssbww.dbf & 文件中存放1亿以内素数s=1do while s5 clear 8,20 say 请选择您需要随机数的变化范围 10,30 prompt 1 09 12,30 prompt 2 099 14,30 prompt 3 0999 16,30 prompt 4 09999 18,30 prompt 5 结束 menu to s do case case s=1 m=11 & 选择模数M mm=10 & 给定随机数的最大边界 do js & 调用计算过程 do dy & 调用打印过程 case s=2 m=101 mm=100 do js do dy case s=3 m=1009 mm=1000 do js do dy case s=4 m=10007 mm=10000 do js do dy case s=5 clear wait 计算完成，谢谢使用！ window at 20,40 timeout 10 endcaseenddo procedure js & 计算过程clear 20,20 say 请输入需要随机数的数量 get kk default 100reada=3x=2select 2 create table ss1(sjshu n ) &准备存放计算产生的随机数select 1go 1 k=1for i=1 to 3*kk/2if xmm insert into ss1 (sjshu) values (x) if k=kk exit endif k=k+1endify=xx=(a*x+ssbqw.shuju)%mskipdo case case x=y and xm-3 x=int(x/2)endcase endforselect 2copy to sjshu.dbf & 将计算结果保存 clear? 已计算随机数数量 k=,kreturnprocedure dy * 打印过程wait 计算完成，打印吗(y/n)？ to yn if upper(yn)=Y select 2 go bottomm=recno()goto top ?for i=1 to m? sjshu skip if i%10=0 ? endif if i%100=0 ? endif endfor wait 随机数打印完毕，按任意键继续！elseclearwait 不打印，按任意键继续！endifreturn几点说明1 运行本程序产生随机数的数量可以达到100万，最大数量可以达到500万以上。如果需要更大数量的随机数，需要用更大的素数文件替换ssbww.dbf。2 如果需要与前次运行产生随机数的位数相同，而分布不同的随机数，你只要修改程序中一个参数m的值就行了。mmm，选择一个接近mm的素数。当然，你也可以修改其它参数来实现。3 如果需要更多位数的随机数，如6位、7位等。只需对程序菜单部分略加修改，就能实现。4 运行程序产生的随机数被复制到