随机数及随机数生成器.doc

我所关注的随机数及其生成器化学物理系（003）SA06003030马瑞作为化学物理系从事计算化学和计算物理方面工作的科学工作者，我们在科研中经常用到随机数及其生成器，确切的说应该是“伪随机数”和“伪随机数生成器”！特别是在计算物理学中的蒙特卡罗模拟，还有对其它许多较为复杂的模型的动力学过程模拟。在这些计算所需的程序中几乎都离不开伪随机数生成器，而且许多对计算机硬件、编程、算法不太在意的同行们，往往就把这里的“伪随机数”当作真正的随机数用了。我曾经听说过一个故事：有一位计算物理学家在使用蒙特卡罗模拟时，发现其所得的结果总是无法与实验结果相吻合。在排除了其它一切可能之后，他把注意力集中在自己所用的伪随机数生成器上，经过一系列繁琐的数学变换，他最终发现自己所求解的计算物理模型涉及到高维空间中的一组平行等距平面切分高维球体的问题，而这个伪随机数生成器产生的坐标所确定的点，恰恰不能在这个高维球面上均匀分布！我的亲身经历其实，我自己也亲眼见证过类似的事件：大三的时候，我曾经写过一个包含有随机过程的Turbo C程序，它能一边模拟一个模型的演化过程，一边实时将模型的当前状态绘制在屏幕上。刚开始，我使用一个伪随机数生成器来产生连续的伪随机数序列，以此作为二维平面上点的一对又一对坐标值，结果竟然在屏幕的实时监控动画上发现“随机”确定的一连串点，竟然在平面上近似沿着一条直线散布！这明显会向我的模型演化结果中引入奇异的现象，当然这些现象将是不真实的。后来，我不时用系统时钟初始化伪随机数生成器，从而产生没有明显分布规律的二维随机点，至少在表面上获得了合乎需要的二维均匀分布。但是自此以后，我一直相当关注随机数及其生成器的问题，特别注意如何有效避免“伪随机数数列”的内在规律对实际使用的干扰，当然还留意一些产生真正随机数而不是“伪随机数”的方法；另一方面，我们从事理论化学和统计物理方面研究的人，常常会对涉及到热力学函数“熵”、“焓”、“内能”的问题感兴趣，恰好关于随机数的研究也可以从“信息熵”的角度来分析。下面，我将主要围绕三个方面谈谈自己的认识：伪随机数的产生方法，真随机数的产生方法，从信息熵的角度来考察计算程序中的随机数与伪随机数。谈到伪随机数，就不得不首先讨论随机数与随机性。什么是随机性？如何证明一个数列是随机数序列？答案是：迄今为止，没有人能给“随机性”作一个恰当的定义，我们也无法证明一个数列的随机性；当然，我们可以证明一个数列不是随机数序列，只要找出一个简单的规律能描述该序列，显然它就不会是真正的随机数序列但是，如果你或者全人类都没有找出简单的规律，也不能说它就是真正的随机数序列，因为也许只不过是因为这个规律很难找，所以我们迄今尚未找到规律而已。举个例子：425270376470090434957375449439 ，它是随机数序列吗？也许很难找到规律，但是，请把每一位减去“1”，你将得到：314159265358979323846264338328 ，它只不过是圆周率的前30位有效数字。换句话说，上述序列425270376470090434957375449439是“前30位有效数字逐一加1”这一简单规律构造的序列，它显然不是真正的随机数列。现代自然科学中，有两类问题可以产生随机数。一类是量子力学，由它产生的随机数至少目前来看无法证明为伪随机数，因为我们迄今为止没有办法揭示隐藏在量子不确定性之下的内在规律，甚至在量子不确定性导致的量子随机现象之下有没有更深层的规律也就是“隐变量”是否存在，理论物理学家们对此都还没有统一的认识。另一类是浑沌动力学系统，它带来的是“确定性随机现象”，换句话说就是由明确而简单的规律产生的随机数序列，这显然是伪随机数，但是需要看到，由于非线性系统的初值敏感性，目前的科学技术很难从已知的浑沌系统伪随机数序列中反推出产生它的规律！所以某些时候并不把浑沌动力学系统产生的随机数称为“伪随机数”。除了上述两类情况，其它的随机数产生办法都得归结为伪随机数生成器。下面着重讨论纯软件的伪随机数生成器，也就是各种编程语言中的（伪）随机数库函数。伪随机数与伪随机数生成器计算机是确定性的机器，因此它无法直接生成真正的随机数，而浑沌系统的随机数生成速度又比较慢，在许多情况下不适合作为快速的（伪）随机数库函数算法。快速的伪随机数生成算法中最著名的要数linear-congruential method（线性同余法），也就是：Xn+1 =（a Xn + b ）% c/ % 就是C/C+ 中的MOD（同余）运算符这种方法可以从一个种子X0 = seed开始，连续生成任意长的伪随机数序列Xn 。它的运算过程极其简单，并且如果令c = 2m，其中m为Xn的字长，则连MOD运算都直接省掉了 Xn+1 2m 时高位自动溢出而被截除。用这种办法生成的伪随机数序列，在给定范围和精度内确实满足均匀分布的要求，但是并非连续分布，因为计算机存放数据的精度不是无限！正是由于最小数据间隙的存在，该序列将会以一个相当长的周期循环。更为重要的是，上述方法为代表的伪随机数生成器，其生成的数字序列并非拥有不可辨别的模式，相反，它们的规律都极其容易由少量数据反推得到。这里有两种具体情况，一类是数据加密中要应用到伪随机数，则此种“易于反推出规律的”伪随机数生成器很不合适；另一类是科学计算中用到的伪随机数，这时并不要求一定得掩盖住其内在规律，但是却有另外的要求：伪随机数的内在规律不能与所研究问题的自然规律相似，并且在必要的时候需要使用高精度的方法减小“最小数据间隙”。之所以称其为“线性同余法”，是因为Xn+1 =（a Xn + b ）% c生成的伪随机数实际上是分布在一条经MOD运算截断的直线上，其近似具有“随机性”的核心也就是这个同余运算。但是如果将相继的一对数（x=X2i ，y=X2i+1）作为一对二维坐标确定平面上的一个点，则在伪随机数数列上的一个小区段内显然会连续出现几对数，它们各自确定的点在平面上沿一直线排列，这也就是文章开头我所提到的现象。其中直线的方程设该小区段为Xp，Xp+1，Xp+2，Xp+ m则相应的点对为（Xp，Xp+1），（Xp+2，Xp+3），直线y = Ax + B =（aX + b）MOD c =（aX MOD c）+ b而（aX MOD c）= aX - Floor aX /c c ，根据算法要求，c相比于aX是一个较大的常数，在一个小区段内总有Floor aX /c Floor aXp /c 。用常数P记Floor aXp /c ，则直线方程可写为：y = Ax + B = aX + b - cP ，此方程当且仅当 cP aXp aXp+ m cP+c时成立。这种二维平面上的规律性严重影响到我的模拟计算，因此不得不反复调用randomize()函数来初始化伪随机数生成器即使是使用两个独立的伪随机数生成器也不能避免类似的规律性。但是randomize()函数是利用系统时钟作为种子实现初始化的，这导致“不断初始化”的办法存在严重的问题：其一，读取系统时钟的函数需要访问CPU-RAM总线以外的慢速设备，因此执行时间相当慢，会严重影响高性能计算的速度。当然，外围时钟可以利用现代操作系统的虚拟时钟设备在CPU内部仿真实现，这部分的缓解了此问题。其二，系统时钟的计时精度并不高，但是今日的计算机“CPU-高速缓存”子系统运行速度非常高，以至于CPU在系统时钟计时精度的最小间隙（微秒级）之内就可以完成上千次操作。如果运行速度很快，因而非常频繁地读取系统时钟来初始化的话，可能会导致相继的几次读取获得完全相同的时钟值，换句话说，就是连续多次使用同一个种子来做初始化，那么显然会获得完全相同的伪随机数数值。其三，以线性同余法为代表的简单伪随机数生成算法，最初生成的几个数受初始化种子的影响很大，所以随机性都不好，一般需要舍弃伪随机数序列的前几个值这种不均匀性问题在上述高速运行的情况下特别严重！因为在非常频繁地做初始化时，每次采用的种子即使不完全相同，也是紧密相邻的。如果舍弃数列中的前几个数，却又会增加运算量、降低效率。总之，频繁利用系统时间初始化这一方法在高速运算时的效果相当差。但是从好的方面看，只要研究者对自己处理的科学计算问题思路明确、图像清晰，一般来说总是有办法避免伪随机数序列的内在规律与所研究问题的自然规律相似；同时只要略微施加一些编程上的技巧，就能有效地以存储空间换运行时间，避免反复调用系统时间做初始化的动作。例如对于上述例子，可以开辟一个二维数组T2max（max2n），然后：randomize() ;int i ;for (i=0 ; imax ; i+)T 1 i = Random(Range_We_Need) ;/* 这里的Random函数根据具体情况返回需要的数据类型 */randomize() ;for (i=0 ; i0 ; count - - )RandomOutput( Constant n ) ;/ 连续输出伪随机数Ask( ENDorNOT ) ; until ( ENDorNOT=Yes ) ;这就是利用混沌动力学系统的初值敏感性进行“反反向工程”的伪随机数产生算法，初值上的微小差异将在一小段时间之后导致非常巨大的分歧，所以“由运行一段时间之后才开始发布的伪随机数序列反推初值”的做法几乎不可能实现；至于如何应付暴力穷举法的正向破解只要数组Constantn 的总长度足够长，并且严密保护好其中任何一个初值都不被窃取，那么穷举所需时间的数学期望也是不可思议的漫长。困难的是，如何对这些即时生成的伪随机数作归一化校正，使任意一小段很短的伪随机数序列也能呈现均匀分布。另外，虽然混沌动力学系统所产生的伪随机数其内在规律不容易与一般应用中的规律相混淆，（即使同样是研究混沌动力学系统，只要模型略有差异就不会有问题，）但是万一有所混淆，查找问题根源就变得极为困难。这里我提供一个权宜之计：收集伪随机数子程序的输出，待收集到足够多的值以后对它们进行FFT，可以尝试一维、二维任何一种你觉得必要的傅里叶变换；若发现变换出的结果有显著的规律性，那么立即再次收集一批不同的输出数据，重复刚才的测试过程数次。假如类似的规律始终存在，那么就说明此伪随机数子程序有明显的规律性，不宜使用；另一方面，如果你没有发现任何显著的规律，也不表明这个子程序就很优秀只是你暂时还没发现缺陷罢了。如何使伪随机数的内在规律变得更加隐蔽？如何使它看起来（并且用起来）更像是一个真正的随机数而不是可预测的结果？这些问题，不仅在科学计算领域被经常提出，并且更为迫切的出现在电子商务、机密通讯等涉及到高度保密的信息传递的领域。这是因为现今比较流行的加密算法几乎都基于一些来自于随机数的密钥，而如果采用伪随机数的话，其内禀规律一旦被别有用心的人发现，那么就有可能推算出密钥，进而破译相应的密文！上述的混沌动力系统可以产生基本符合要求的伪随机数序列，但是更符合需要的随机数应该是来自于硬件设备的真随机数硬件随机数生成器伪随机数并非是“坏随机数”，只是某些情况下为稳妥起见我们不得不使用“真随机数”，产生真随机数的办法必须依赖硬件，其原理不外乎于量子力学和混沌动力系统。量子力学的随机现象迄今没有更深层的规律性，因此利用放射源衰变等过程得到的随机数序列肯定是真正随机的，在经过归一化校正后，便可以作为优秀的随机数序列发布了。但是，混沌动力系统不是“伪随机”的确定性系统吗？它怎么能产生“真随机数”呢？我们说纯软件的混沌动力学系统伪随机数生成器产生的是伪随机数，这是因为现有软件运行在完全确定性的计算机上，并且计算机可以存放的数值虽然离散、但却准确无误例如因为有限位二进制存储方式离散的缘故不能准确存放实数X=3.00000，可是却能准确存放例如X=2.9999983 E+00这样一个完全确定的近似值。那么，在各种条件和变化过程都完全确定的计算机中，产生的随机数当然是可以预测的“伪随机数”。现实世界中的硬件设备就不同了，例如有些硬件利用的是电子元件温度的随机波动经电子学放大而产生随机数，理论上它应该满足一个“热传导-流体力学”复合的非线性动力学方程组，但实际上并没有可能利用这样的方程组来预测随机数产生情况：首先，我们无法精确确定该模型的初始条件，许多物理量最多只能测到46位有效数字。这一方面是因为测量工具自身精度的局限性；另一方面是因为客观世界的微观涨落使这些初始值快速小幅变化着；甚至在足够小的微观尺度上还会撞上量子的“不确定性原理”。种种原因都导致很难在客观世界中同时测出一系列初值的精确值，相比之下，我们可以在计算机中的模型上指定任意多个足够准确（高达815位有效数字）的初始值，并强制令其为“同一时刻”的值。其次，计算机模拟过程中很容易限定环境条件（在815位有效数字的精度上）保持稳定，而在现实世界中这是不可能的！以PCI卡形式的随机数生成器为例，机箱内空气的湿度、温度和速度矢量分布上的微小波动，都可能显著影响到它的散热情况，又如果把它们都考虑为扩展模型的变量，那又会大大增加模型的初始条件数目，使“同时测出一系列初值的精确值”变得更加不可能。可见，利用现实世界中测不准初值的混沌动力学系统，因为它的初值敏感性，我们也可以获得几乎完全不能预测的随机数序列。对于基于量子力学的硬件产生的无限随机数序列，预测它是不可能事件；而基于混沌动力系统的硬件产生的无限随机数序列，预测它是零概率事件，亦即说，对于可以预见的未来的任何一种应用，都可以放心的视为“不可能事件”！总之，我们把依赖量子力学或者混沌动力系统的硬件随机数生成器统称为“真正的”随机数生成器。一个长半衰期的放射源（如T1/2=30年的Cs-137）加上一个盖格计数器，辅以有效的防护设备和计算机数据接口卡，还要再加上一个调试好参数的归一化校正程序，就可以组成一套不错的硬件随机数生成器，并且是基于量子力学的完全不可预测的随机数产生装置。但随着放射源的嬗变和仪器设备的慢慢老化，这个归一化校正程序需要每过一段时间就校正一次参数，这倒不是什么大问题；真正的大问题在于，它产生随机数的效率太低了！为保证足够的测量精度，探测器的增益不可以太大，这样该装置的随机数输出带宽受放射源辐射强度的限制，可能只有每秒几个bit甚至几十秒才有一个bit ！（电阻或者二极管中）电子热运动的Johnson噪声是量子效应，根据这个原理产生随机数的硬件装置可以更快地输出随机数，例如ComScire或者Tundra公司的一些设备，可以以20 kb/s的带宽输出质量优异的真随机数。Intel RNG的带宽为75kb/s，并且几乎不需要花费额外的费用，因为这个装置所需的硬件是Intel 82802 Firmware Hub（也就是BIOS设备）中的标准配件。从2003年起，威盛电子（VIA）在其CPU中嵌入了PadLock RNG模块，由于采用CPU内核级别的设计和生产工艺，它的访问速度很快，并且号称能以8001600 kb/s的速率生成真随机数折合成单精度浮点数却也不过2550 kHz 。VIA PadLock RNG的技术原理是联合使用一组四个集体管震荡器，其随机性本源还是来自Johnson噪声。正因为硬件随机数生成器产生的数字完全是随机的，即使是设计和使用它的人也无法找到这些真随机数的生成规律，所以用一个软件来精确的校正其分布变得相当困难，必须使用大量统计数据获得经验校正函数，而且还不能保持非常好的校正精度，这一点对加密应用来说问题不大，但是对于科学计算而言就会大大降低计算结果的精度。硬件随机数生成器最大的缺点在于：它们都太慢了！即使是VIA PadLock RNG的2550 kHz单精度浮点数产生速率，也远远跟不上当前主流x86系统的25GFLOPS / CPU的浮点计算速度：单CPU工作时差不多每十万次浮点操作才能等候到一个单精度随机数，用于计算密集型的科学研究工作显然不够。说到这里，就不得不提一些半硬件半软件的解决方案，它们不需要专门的硬件设备，可以跨平台的应用于几乎任何一种计算机上，比Intel或者VIA的RNG还容易获得；但它们又不是完全的软件伪随机数生成器，而是至少近似不可预测的（半真半伪的）随机数来源，从本质上来说，它们是混沌动力系统硬件随机数生成器的一类：这些生成器中，有些捕捉键盘、鼠标、硬盘驱动器的动作以及外设中断信号（DMA等），并由此作为随机数的来源；有些则把互联网上的不确定行为收集起来作为随机数来源，例如本机的网络信号流量波动情况；还有一些会从系统硬盘上任意选定的声音、图像、动画文件中提取随机信号，或者从声卡、调制解调器之类的AD/DA转换器中提取背景噪音。这些不确定行为的最根本来源是个人或者人类社会的活动，而人或者人类社会是非线性的耗散系统，显然可以认为上述随机数生成方法都是基于混沌动力系统的。被采样的行为虽然看起来将会随机发生，但其实这种随机性却无法证明，如果这些物理动作是静止的（比如说服务器的鼠标和键盘可能根本就没连上）或者是简单重复的，就可能生成一系列低质量的随机数，由此导致随机数硬件源的质量远远低于设计者在最理想条件下的预期值。不过，虽然质量较低且不稳定，这些硬件却比那些直接利用放射源计数或者电子热运动等现象制造的专用硬件廉价、常见得多，并且已经在系统层或者驱动程序层实现了通用的随机数池设备，从使用者的角度来看更接近软件解决方案无论是它的易用性还是实现方式。它们当中最著名的一种，就是Linux 1.3.30版以上内核开始在系统层提供的随机数池虚拟设备“/dev/random”。“/dev/random”设备驱动程序收集来自其它硬件的动作信息：鼠标、键盘、磁盘驱动器、各种接口卡等等，将这些未经校正的随机信息填充到一个通用的“（信息）熵池”中，然后在需要访问“/dev/random”设备时提供合适的统计学函数校正和转换。当然，它收集信息的速率也很慢，甚至比那些基于量子效应的专门硬件还慢！不过一般而言质量确实也能接近基于量子效应的硬件随机数生成器。由于运作速度慢，“/dev/random”的“熵池”在频繁读取时将很快枯竭，因此我们不应该把这个设备当作快速获取大量高质量随机数的来源，而应该从另一个相关的虚拟设备“/dev/urandom”中获取。“/dev/urandom”设备以“/dev/random”中的随机数做种子，利用MD5算法生成纯软件的伪随机数，其效果显然比简单的线性同余法高许多，相应的速度也会低一些。调用这两个虚拟设备的方法很简单，本文最后的附录中给出了相应的C99程序实例。信息熵以及收集它们的硬件前述各种“真随机数”生成器都是硬件的，而纯软件方法生成的随机数无论如何只能叫做“伪随机数”，为什么呢？这得从信息熵的角度来分析。联合考虑信息熵H与热力学熵S，可以理解为：越是混乱无序、无规律的系统，其所蕴含的信息量H越大，或者说描述它所用的数学模型的信息量大；同时，系统的热力学熵S也越大。S最低（=0）的系统最有序，好比全零的信号流序列00000000，显然记录它只要一条简单的规则“全为0的序列”；再如前文提到的衍生自圆周率的那串序列425270376470090434957375449439，也只需要一条简单的规则“前30位有效数字逐一加1”，因此它的信息熵似乎也不大（下面额外讨论一下），同理还有那些利用混沌动力学系统生成的伪随机数序列，也因为有简单的生成规则而只能归类为“伪随机数”。关于“前30位有效数字逐一加1”的额外讨论：这一规则是否为一简单的、低熵的信息？我觉得它其实并不简单！为更加明确，我们讨论无穷长的数列“的全部有效数字”。诚然，无穷长数列“全为0的无穷序列”信息量是很低，但是“的全部有效数字”这一规则中的“”，似乎不能与“0”等同为相同的信息量，我认为“”本身所包含的信息量即为无穷大，因此该无穷长序列的H亦为无穷大。“”只不过是专门用于称呼那一特殊无理数的简称，是因为我们约定了这个简称才使具有“圆周率”的含义，如果认为该无理数只不过是一个希腊字母（的信息量），那就未免太形而上学了。否则，对于给定的一串无穷长真随机数序列，我也可以命名其为例如，难道它就变成低熵序列了？既然可以用较简单的规则来概括“伪随机数序列”看上去大得多的信息量，这就说明该序列所包含的信息存在大量的冗余，或者说这些信息可以被压缩。为什么有的“半硬半软”随机数生成器会选择从磁盘上的图像、动画或者声音文件中读取随机数据，特别是各种压缩格式诸如JPEG、PNG或者MP3、RMVB等等，就是因为这些文件中的信息已经经过较好的压缩算法处理，去除了大部分冗余的信息。不讨论基于量子效应的专用硬件随机数生成器，那些更常见的硬件设备的许多活动都是（近似）不可预测的，因此我们通常记录这些活动以作为随机数的“熵源”，例如最常用的硬件设备活动系统时钟变化。但麻烦的是常见硬件设备的动作大多涉及到一些机械-电子的复合运动，即使是系统时钟也需要石英晶体的压电振动，而机械运动的速度显然远远落后于现代数字式元件的工作速度，因此利用常见硬件设备产生高质量随机数的速度都很慢。而且这里存在着“熵源”质量不稳定的情况，不过从信息熵和冗余信息的关系来看，可以利用压缩算法提高输出的随机数序列的质量，就是说先缓存待输出的序列，然后应用高质量的压缩算法对缓存内容进行压缩，这样就能把因为质量下降而带来的冗余信息尽可能除去，最后再输出处理过的序列。为了高速的产生大量（伪）随机数，我们不得不使用纯软件的伪随机数生成器，而任何纯软件随机数生成器本身是不能自发生成信息熵的，它必须有一个来自硬件的输入也就是一个来自“熵源”的数据。按通常的算法描述，就是要先对伪随机数生成器进行初始化，例如用系统时间作种子来初始化，这一操作在32位的C99标准中会输入32bits的信息熵，此后每一次调用系统时间“重新初始化”理论上都会引入32bits信息熵。这里有两点需要注意：第一，两次初始化之间，该伪随机数生成器输出的序列总共只有32bits的信息熵，而无论到底输出了多少bit的信息！当然，并不妨碍在许多情况下正常的使用这些低熵的数据。第二，如果非常频繁的调用系统时间“重新初始化”，则前后两次获取的系统时间数值之间的相关性太大，就不会重新引入32bits信息熵了。因此有必要适当降低调用系统时间的频率，或者使用其它更快的硬件随机数生成器作为“熵源”。一个优质的硬件随机数生成器“熵源”，适度频繁的为一个伪随机数生成算法“重新初始化”，便能快速输出大量较高质量的伪随机数，其实，这正是Linux设备“/dev/random”和“/dev/urandom”的做法。随机过程课程期中作业19附录：源代码调用Linux虚拟设备“/dev/random”（secure）或“/dev/urandom”（strong）的源代码，C99语言标准。这两个设备中的随机数可以按32bits无符号整数来连续读取，即32位C标准的unsigned int 。/* * Written by John Viega () * Free for any use. */#include stdio.hint secure_inited = 0;int strong_inited = 0;FILE *rand_file;FILE *urand_file;void _init_secure() rand_file = fopen(/dev/random, r); secure_inited = 1;void _init_strong() urand_file = fopen(/dev/urandom, r); strong_inited = 1;/* * Return a secure random number from 0 to 232-1 */inline unsigned int secure_rand_base() unsigned int num = 0; int read_ret = 0; if(!secure_inited) _init_secure(); while(!read_ret) read_ret = fread(num , sizeof(unsigned int), 1, rand_file); return num;/* * Return a cryptographically strong random number from 0 to 232-1 */inline unsigned int strong_rand_base() unsigned int num = 0; int read_ret = 0; if(!strong_inited) _init_strong(); while(!read_ret) read_ret = f