随机数是许多应用中的关键因素, 在数值模拟计算、机器学习、密码学等方面起到非常重要的作用. 任何确定性算法都无法产生真正的随机数, 因此伪随机数发生器(pseudo random number generator, PRNG) 产生的随机数序列, 严格意义上来说不是真正的随机数序列, 但在实际应用中希望这些序列在某种意义上类似于真正的随机数序列, 也就是具有真正的随机数序列的一些统计性质, 比如独立性、均匀性和不可预测性等. 有两类测试方法用于评估PRNG的性能——理论测试和统计测试. 理论测试检查给定发生器的内在结构, 不一定需要生成随机数序列, 例如晶格测试^[1]和光谱测试^[2], 这类测试对每一类发生器都是不同的. 而统计测试是在PRNG产生的序列上进行的, 不需要知道它是如何产生的. 这类测试的主要目的是检查由PRNG产生的伪随机数序列是否具有类似于真正随机数序列的统计性质, 比如独立性、均匀性等.

统计测试试图找到反对零假设($ H_{0} $: 生成的伪随机数序列服从独立且均匀分布)的统计证据. 因为任何关于有限个服从独立且均匀分布的随机变量的函数变量, 在零假设下的分布或者近似分布是已知的, 可以用作统计测试. 统计测试旨在检测生成序列中的各种偏差, 以及揭示发生器的缺陷. 一些测试方法中包括常见的统计工具, 如K-S检验、卡方检验等, 用于比较PRNG产生的伪随机数序列通过适当计算所得统计量的理论和经验分布. 常见的统计测试方法有参数检验、游程检验、扑克检验、卡方检验以及序列检验等, 这些测试方法利用了独立且均匀分布的随机数序列的特性. 然而单一的统计测试只关注随机数序列的均匀性或者独立性等某个特定的属性. 因此, 许多研究者出于实用目的创建了各种测试套件, 测试套件由一系列独立的统计测试方法组成. 如果PRNG产生的伪随机数序列通过了给定测试套件中的所有测试, 则认为它是好的. 但是这类测试无法证明发生器生成的伪随机数序列具有良好的随机性, 只是增加了对模拟计算结果的信心. 正如文献[3]指出, 它们实际上是对非随机性的测试. 此类测试套件包括DIEHARD^[4], TestU01^[5]以及NIST测试套件^[6]等, 其中NIST测试套件是由美国国家标准与技术研究所设计, 目前被认为是最先进的测试套件之一.

PRNG的随机性通常以能否通过检验均匀性、独立性或者不可预测性的统计测试来衡量, 但所有PRNG相比较, 仍然无法确定哪种发生器最好, 因为它们没有定量分析, 现有的随机性测试仍缺乏PRNG之间的比较手段^[7]. 因此, 分析每个PRNG的随机性性能, 或者比较PRNG的随机性仍是一项重要任务. 单个统计测试的结果通常以$ p $值的形式给出, 若$ p\geqslant \alpha $, 则接受$ H_{0} $; 若$ p<\alpha $, 则拒绝$ H_{0} $. 其中$ \alpha $是显著性水平($ \alpha $通常取0.01、0.05等), 这种测试方法通常被称为一级测试^[8], 一级测试一般只关注随机数序列的均匀性或独立性等某一个属性, 只有非常差的发生器无法通过这些简单一级测试方法. 另一种测试方法是二级测试, 考虑给定发生器生成的多个不相交序列的一级测试结果, 通过拟合优度检验比较这些结果$ p $值的理论分布和经验分布^[9]. 二级测试方法在一级测试结果$ p $值的基础上进行, 在测试过程中产生的误差较小, 测试结果可靠, 常用的测试套件包含了大量的一级测试方法. 虽然单个测试的结果有明确的统计解释, 但如何解释一个测试套件的结果却没有明确的说明, 已有研究表明NIST测试套件具有一定的局限性, 存在第二类错误, 以及具有推荐参数的NIST测试套件无法检测到一些PRNG的已知缺陷^[10].

近年来, 一些新的测试方法被引入用来对PRNG进行严格的测试. 这类方法的核心思想是: 选择一些适用于真正随机数序列的概率定律, 并将相应统计量的理论分布与给定PRNG生成的伪随机数序列的经验分布进行比较, 比如可以根据在零假设下经过适当计算所得统计量的$ p $值来完成比较. 2015年, Wang等人根据重对数定律LIL提出了基于统计距离的伪随机数序列测试方法, 可用于识别几种常用PRNG实现中的缺陷, 这些缺陷尚未被NIST SP800-22测试套件检测出来^[10]; 2020年, Lorek等人提出了一种基于反正弦定律的随机游动的二级统计测试方法, 该方法可以检测出许多常见的PRNG的弱点, 还能够揭示生成序列中其他工具(如NIST SP800-22测试套件)无法识别的一些缺陷和规律^[8]. 然而这几种测试方法都是测试二进制随机数序列, 属于二级测试方法. 本文尝试提出一种检验$ (0, 1) $上的实数序列的一级测试方法.

1998年, Rudnick和Sarnak提出了$ [0, 1) $上实数序列$ \left\{ x_{n}\right\}_{n\geq 1} $的泊松对相关(Poissonian pair correlations, PPC) 的概念^[11], 认为服从独立且均匀分布的$ [0, 1) $上的实数序列都满足泊松对相关, 文献[12] 证明了满足泊松对相关的$ [0, 1) $上的实数序列一定服从均匀分布, 并且认为随机数序列产生的PPC折线图的收敛速度可以用来衡量序列的均匀性与独立性. 本文在此基础上提出了一个能够同时检验$ (0, 1) $上伪随机数序列的均匀性和独立性的一级统计测试方法, 并且给出利用幂函数拟合选取收敛判别标准的具体方法. 测试结果显示LCG、RANDU、MT19937-32、SPSS.RV.UNIFORM函数、多重递归发生器、BBS发生器、Matlab.rand函数、Python.random库中random函数和Excel.RAND函数这些发生器的性能差异不显著, 而这几种发生器的性能明显优于基于无理数$ \pi $重叠产生的PRNG. 此外, 测试结果表明该一级统计测试方法能有效地同时检验伪随机数序列的均匀性和独立性.

论文结构组织如下: 第二节介绍了典型的伪随机数发生器、泊松对相关的概念以及相关定理; 第三节引入基于泊松对相关的PRNG的测试方法; 第四节介绍分析PPC测试对多个PRNG的测试结果; 最后总结全文.

伪随机数发生器是一种能够产生具有某些随机性的数字序列的确定性算法, 利用数学公式迭代产生随机数. 下面按照Asmussen和Glynn的描述给出一个严格的定义^[13].

定义2.1  一个伪随机数发生器定义在一个5元组$ \left \langle E, V, s_{0}, f, g \right \rangle $上, 其中$ E $是一个有限状态空间, $ V $是一组值, $ s_{0}\in E $是一个种子, 即序列$ \left\{ s_{i}\right\} _{i=0}^{\infty } $中的初始状态, 函数$ f:E\rightarrow E $描述了连续状态$ s_{n+1}=f\left ( s_{n} \right ) $和$ g:E\rightarrow V $之间的转换, 将发生器的状态映射到输出.

对于$ M\in N $, 通常$ V=\left ( 0, 1 \right ) $或$ V=\left \{ 0, 1, 2, \cdots , M-1 \right \} $, 本文主要研究$ (0, 1) $上的随机数序列, 所以对于每一个发生器都取$ V=\left ( 0, 1 \right ) $. 比如线性同余发生器是一个根据公式$ s_{n+1}=\left ( as_{n}+c \right ) $mod$ M $更新其状态的发生器. 理想的真随机数发生器产生的二进制随机位序列, 输出序列的每个位都是独立的, 所以1和0出现的概率都为$ \frac{1}{2} $; 同样对于$ (0, 1) $上的真随机数序列, 每个数都是独立生成的, 出现的概率相等, 并且服从$ (0, 1) $上的均匀分布, 即如果将$ (0, 1) $等分为$ n $个子区间, 那么序列中所有的数落在每个子区间的个数都为$ \frac{N}{n} $, 其中$ N $为该随机数序列的长度, 并且任意一个数落在某子区间的概率与其他数无关.

本文提出的PPC测试方法被用于测试不同的伪随机数发生器, 包括线性同余发生器、RANDU、MT19937-32^[14]、Matlab.rand函数、BBS发生器、Python.random库中random函数、Excel.RAND函数、SPSS.RV.UNIFORM函数、多重递归发生器以及基于无理数$ \pi $重叠产生的PRNG, 此外, 对由$ \sqrt{n} $的小数部分生成的数字序列进行测试, 反向检验PPC测试方法是否能够检测出该"发生器"的缺陷.

线性同余发生器(Linear Congruential Generator, LCG)是利用求余运算来产生随机数, 其递推公式为

(2.1) $ \begin{equation} \left\{\begin{array}{lll} s_{n+1}=\left ( as_{n}+c \right ){\rm mod} M, \\ { } r_{n}=\frac{s_{n}}{M}, \end{array}\right. \end{equation} $

其中$ n=0, 1, 2, \cdots $, $ r_{n} $是$ (0, 1) $上的伪随机数序列, $ M $是模数, $ a $是乘数, $ c $是增量, 并且$ a, c< M $, $ s_{0} $为初始种子. 线性同余发生器中参数$ c=0 $时, 称之为乘线性同余发生器, 或积式发生器. LCG已被应用在各种编程语言中, 本文使用的参数为$ a=25214903917 $, $ c=11 $, $ M=2^{48}-1 $, 这是Java库java.util.Random中的默认值. RANDU是一个乘线性同余发生器, 参数为$ a=65539, M=2^{31} $.

多重递归发生器(Multiple Recursive Generator, MRG) 是线性同余发生器的一个自然推广, 其递推公式为

(2.2) $ \begin{equation} \left\{\begin{array}{lll} s_{n+1}=\left ( a_{1}s_{n}+\cdots + a_{k}s_{n-k+1}\right ){\rm mod} M\\ { } r_{n}=\frac{s_{n}}{M} \end{array}\right. n\geqslant k-1, \end{equation} $

其初始值$ \left ( s_{0}, \cdots, s_{k-1} \right ) $为非负向量, 当$ k=1 $时, MRG便转换成了LCG. 与LCG相比, MRG的最大周期更长, 但是运行时间也相对较长. 本文所测试的多重递归发生器是MRGk5-93^[15], 递推公式为$ s_{n}=\left ( 107374182s_{n-1}+104480s_{n-5} \right ) $mod$ \left ( 2^{31}-1 \right ) $.

Blum Blum Shub(BBS)生成器是由Lenore Blum, Manuel Blum和Michael Shub于1986年共同提出的一种非线性同余发生器, 其递推公式为

(2.3) $ \begin{equation} \left\{\begin{array}{lll} s_{n+1}=s_{n}^{2}{\rm mod} M, \\ { } r_{n}=\frac{s_{n}}{M}. \end{array}\right. \end{equation} $

其中$ M=pq $, 而$ p, q $为两个大素数, 且$ p, q $皆为$ 4z+3 $的形式, 本文使用的参数为$ p=11234773052181932039 $, $ q=15755662711309472467 $. BBS发生器的安全性很好, 并且通过了几乎所有的理论检验, 但由于其运行速度很慢, 很少用于模拟领域.

基于无理数$ \pi $重叠产生的PRNG利用$ \pi $的十进制扩展的小数部分重叠产生伪随机数序列, 每次处理8位数字的数字串, 并将所得的8位数除以$ 10^{8} $, 这样得到的每个随机数都在$ (0, 1) $范围内. 由于存储空间的限制, 本文选择重叠法产生伪随机数序列, 也就是每次处理完一个8位的数字串后, 指针仅向后移动1位, 得到一个新的8位数, 依次向后产生数字序列. 本文从网站http://stuff.mit.edu/afs/sipb/下载了所需要的$ \pi $的小数部分.

1998年Rudnick和Sarnak提出了在$ [0, 1) $上实数序列$ \left\{ x_{n}\right\}_{n\geq 1} $的泊松对相关(Poissonian pair correlations)的概念^[11], 引起了学者们的广泛关注. 如果对于任意$ s>0 $, $ [0, 1) $上的实数序列$ \left\{ x_{n}\right\}_{n\geq 1} $满足

(2.4) $ \begin{equation} \lim\limits_{N\rightarrow \infty }\frac{1}{N}\sharp\left \{ 1\leq l\neq m\leq N:\left \| x_{l}-x_{m} \right \| < \frac{s}{N}\right \}=2s, \end{equation} $

则称该序列满足泊松对相关, 其中$ \left \| \cdot \right \| $表示到最近整数的距离, $ \sharp A $表示集合$ A $中元素的个数. 为方便后文表述, 记

(2.5) $ \begin{equation} f_{N}\left ( s \right )=\frac{1}{N}\sharp\left \{ 1\leq l\neq m\leq N:\left \| x_{l}-x_{m} \right \| < \frac{s}{N}\right \}. \end{equation} $

泊松对相关是对$ [0, 1) $上实数序列中各元素间隔行为的描述, 本文关注$ [0, 1) $上实数序列的泊松对相关性与服从独立且均匀分布的性质之间的关系. 对于在$ (0, 1) $上服从独立且均匀分布的随机实数序列$ \left\{ x_{n}\right\}_{n\geq 1} $有如下结论.

定理2.1^[16]  若实数序列$ \left\{ x_{n}\right\}_{n\geq 1} $在$ [0, 1) $上相互独立且服从均匀分布, 那么对于任意的$ s>0 $, 该实数序列满足

(2.6) $ \begin{equation} \lim\limits_{N\rightarrow \infty }\frac{1}{N}\sharp\left \{ 1\leq l\neq m\leq N:\left \| x_{l}-x_{m} \right \| < \frac{s}{N}\right \}=2s{\quad}{\rm a.s., } \end{equation} $

即当$ N\rightarrow \infty $时, $ f_{N}\left ( s \right ) $几乎必然收敛到$ 2s $.

证   令$ X=\left| x_{l}-x_{m}\right| $, 记$ X $的分布函数为$ F\left ( z \right ) $, 则

$ F\left ( z \right )=P(X\leq z)=P(\left| x_{l}-x_{m}\right|\leq z)=P(-z\leq x_{l}-x_{m}\leq z ). \nonumber $

记$ (x_{l}, x_{m}) $的联合密度函数为$ f(x_{l}, x_{m}) $, 由于序列$ \left\{ x_{n}\right\}_{n\geq 1} $相互独立, 所以有

$ f(x_{l}, x_{m})=f(x_{l})f(x_{m})=\left\{\begin{array}{lll} 1, & 0\leq x_{l}, x_{m}\leq 1, \\ 0, & {\rm else}. \end{array}\right. $

因此

$ \begin{eqnarray*} F(z) &=&P(-z\leq x_{l}-x_{m}\leq z)\\ &=& {{\int\!\!\!\int}}_{-z\leq x_{l}-x_{m}\leq z } {\rm d}F(x_{l}, x_{m})\\ &= & {{\int\!\!\!\int}}_{-z\leq x_{l}-x_{m}\leq z } f(x_{l}, x_{m}){\rm d}x_{l}{\rm d}x_{m}\\ &=&\left\{\begin{array}{lll} 1-(1-z)^{2}, & z\in (0, 1), \\ 0, & {\rm else.} \end{array}\right. \end{eqnarray*} $

于是

$ f(z)=\left\{\begin{array}{lll} 2(1-z), & z\in (0, 1), \\ 0, &{\rm else}. \end{array}\right. $

令随机变量$ Y=\left\| x_{l}-x_{m}\right\| $, 易得

$ Y=\left\{\begin{array}{lll} X, &{ } X\in [0, \frac{1}{2}], \\ [3mm] 1-X, {\quad} & { } X\in (\frac{1}{2}, 1]. \end{array}\right. $

那么

$ P(Y\leq a)=P(X\leq a)+P(1-X\leq a)=\int_{0}^{a}2(1-x){\rm d}x + \int_{1-a}^{1}2(1-x){\rm d}x=2a.\nonumber $

所以随机变量$ Y=\left\| x_{l}-x_{m}\right\| $在$ [0, \frac{1}{2}] $上均匀分布, 即

$ P(Y\leq \frac{s}{N})=\frac{2s}{N}, \nonumber $

则有

$ \lim\limits_{ N\to \infty}f_{N}(s)= \lim\limits_{ N\to \infty}\frac{1}{N}\cdot (N^{2}-N)\cdot \frac{2s}{N}=2s. $

即对任意$ \varepsilon>0 $, 存在$ k>0 $, 当$ N>k $时, $ \left| f_{N}(s)-2s\right|<\varepsilon , $可得

$ \lim\limits_{k \to \infty}P\left\{ \bigcup\limits_{N=k}^{\infty}(\left| f_{N}(s)-2s\right|)\geq \varepsilon \right\}=0. $

等价于

$ P\left\{ \lim\limits_{N\to \infty} f_{N}(s)=2s \right\}=1.\nonumber $

即当$ N\rightarrow \infty $时, $ f_{N}\left ( s \right ) $几乎必然收敛到$ 2s $. 证毕.

泊松对相关与均匀分布之间的关系同样也备受关注, Larcher和Grepstad以及Aistleitner等人都证明了满足泊松对相关的$ [0, 1) $上的实数序列$ \left\{ x_{n}\right\}_{n\geq 1} $一定是服从均匀分布的^{[12, 17]}. 服从均匀分布是一个整体的统计性质(测试间隔始终保持不变), 而泊松对相关是一个高度局部化的统计性质(测试间隔的大小与$ N $成比例缩小)^[18]. 但是这两个属性不是独立的, 满足泊松对相关的序列一定服从均匀分布, 但是服从均匀分布不一定意味着泊松对相关^[19], 例如Kronecker序列$ \left\{ (n\alpha) \right\} _{n\geq 1} $, 其中$ \left ( \cdot\right ) $表示取小数部分, 对于每个无理数$ \alpha $都服从均匀分布, 但对于任何$ \alpha $值都不满足泊松对相关^{[20, 21]}. 然而$ [0, 1) $上服从独立且均匀分布的实数序列满足泊松对相关, 因此$ [0, 1) $上的实数序列满足泊松对相关这一行为除了可以说明该序列是均匀的之外, 也体现了该序列一定程度上的独立性^[22].

序列服从均匀分布可以看作是一种伪随机性, 而$ (0, 1) $上独立且均匀分布的实数序列随着$ N\rightarrow \infty $, $ f_{N}\left ( s \right ) $几乎必然收敛到$ 2s $, 所以同均匀分布一样, 泊松对相关也可以看作是一种伪随机性. 显然序列$ \left\{ x_{n}\right\}_{n\geq 1} $服从独立且均匀分布是比满足泊松对相关更强的性质. 所以本文类似地做出零假设($ H_{0} $: 生成的伪随机数序列服从独立且均匀分布), 则在零假设的情况下, 该序列一定满足泊松对相关; 若该生成序列不满足泊松对相关, 则意味着该序列在$ (0, 1) $上不服从独立且均匀分布, 不具有真随机数序列的特性. 此外如果生成的伪随机数序列不是均匀的, 那么该序列一定不满足泊松对相关; 如果伪随机数序列只是在$ (0, 1) $上均匀分布但不独立, 由Kronecker序列可知, 该序列也不一定会满足泊松对相关. 所以基于泊松对相关提出的随机性统计测试方法在一定程度上能够兼顾检验生成序列的均匀性和独立性, 可用于揭示伪随机数序列中的一些非随机性.

本节利用$ [0, 1) $上实数序列的泊松对相关性来设计一个测试PRNG随机性的实用方法, 称之为PPC测试. 该方法仍然采用假设检验的方法, 首先做出零假设($ H_{0} $: 生成的伪随机数序列服从独立且均匀分布). 测试的总体思路如下: 由给定的PRNG生成若干个长为$ M $的随机实数序列, 选取一个合适的间隔$ k $, 则每一个实数序列可以被分为$ m $个重叠子序列$ \left\{ x_{n}\right\} _{1\leq n\leq k} $, $ \left\{ x_{n}\right\}_{1\leq n\leq 2k} $, $ \cdots $, $ \left\{ x_{n}\right\}_{1\leq n\leq mk} $, 多余的数可以忽略, 相对应每个子序列的长度分别为$ k, 2k, \cdots, mk, $ ($ mk\leqslant M $, $ m $应该尽可能大一点). 在零假设下实数序列随着$ N\rightarrow \infty $, 对于任意的$ s>0 $, $ f_{N}\left ( s \right ) $几乎必然收敛到$ 2s $, 则$ f_{N}\left ( s \right )-2s $几乎必然收敛到0, 所以对于在$ (0, 1) $上的随机数序列$ \left \{ x_{n} \right \}_{1\leq n\leq N} $, 任意给定$ s $, 则当$ N\rightarrow \infty $时, 有

(3.1) $ \begin{equation} R_{N}:=f_{N}\left ( s \right )-2s\mathop{\rightarrow }\limits^{\rm a.s.} 0, \end{equation} $

则PRNG的每一组生成序列都可以产生一个PPC序列$ \left\{R_{N} \right\}_{N=k, 2k, \cdots , mk} $, 得到一条PPC折线, 因此由给定发生器生成的若干随机数序列可以得到类似于图 1的PPC折线图.

此时判别伪随机数序列是否满足泊松对相关可转化成判断当$ N\rightarrow \infty $时, PPC序列$ \left\{R_{N} \right\}_{N=k, 2k, \cdots , mk} $是否几乎必然收敛到0. 数列几乎必然收敛于0, 相当于对于任意的$ \varepsilon>0 $, 当$ n $足够大时, 数列至多有有限个点落在0的$ \varepsilon $邻域之外. 受该定义的启发, 本文尝试通过控制PPC序列$ \left\{ R_{N}\right\} $中的每个元素来建立PPC序列几乎必然收敛于0的判别标准^[23]. 选定极限值0的一个邻域, 如果PPC序列落入该区域内的点占总数的百分比$ p $超过给定的阈值$ \beta $, 则认为该PPC序列几乎必然收敛于0^[24-25]. 因此本文需要选取一个收敛判别区域, 将PPC序列落入该区域内的点占总数的百分比$ p $作为测试统计量, 若$ p\geq \beta $, 则认为PPC序列$ \left\{ R_{N}\right\} $收敛于0, 测试的伪随机数序列具有泊松对相关性, 接受零假设, 认为该测试序列在$ (0, 1) $上服从独立且均匀分布; 否则, 认为测试的伪随机数序列不具有泊松对相关性, 拒绝零假设, 认为测试序列是非随机的. 该方法主要是根据不同判别标准下$ p $值的大小来衡量发生器的性能.

由图 1可以看出, PPC折线图的波动随着$ N $增大而逐渐变缓, 图中第一象限的收敛趋势与幂函数$ y=x^{b}, (x>0, b<0) $图像的走势相近, 并且只要选择合适的参数, 幂函数图像可以包含PPC折线图中的大部分波动. 因此本文用参数$ \alpha(\alpha>0) $控制收敛速度, 用参数$ c $来控制收敛域的宽度, 即以$ \left ( -\frac{c}{N^{\alpha }}, \frac{c}{N^{\alpha }} \right ) $区域作为判别区域. 总的来说, 收敛判别标准依赖于收敛速度$ \alpha $, 收敛域的宽度$ c $以及阈值$ \beta $.

PPC测试方法的重点在于参数$ s $和收敛判别标准的选取, 即参数$ (\alpha, c, \beta) $的确定. 首先考虑参数s的选取, 当$ s $取不同的值时, 对PRNG生成的伪随机数序列产生的PPC序列进行分析. 由不同的发生器生成若干个长为$ 2^{15} $的伪随机数序列, 选取间隔$ k=2^{7} $, 任取$ s>0(0<s<k/2) $, 本文$ s $的取值为$ 10, 20, 30, 40, 50 $, 同一个伪随机数序列在$ s $不同时, 会得到不同的PPC序列. 图 2给出分别由LCG, MT生成的两个伪随机数序列在$ s $取不同的值时产生的不同的PPC折线. 此外在$ s $不同时, 我们研究了由其他十余种常见发生器生成的随机数序列产生的PPC折线, 结果显示$ s $的不同取值并不影响折线图的最终走向, 也就是PPC序列是否收敛的趋势. 但是通常当$ s $取值变大时, 产生的PPC序列$ \left\{ R_{N}\right\}_{N=k, 2k, \cdots , mk} $上的每一个值在一定程度上也有变大的趋势, 所以$ s $的不同取值会影响收敛域的宽度.

为了更直观地感受$ s $的不同取值对收敛判别区域的影响, 图 3给出了当$ s=20, 40 $时, MT生成的100组伪随机数序列产生的PPC折线图, 以$ \left ( -\frac{7}{N^{2/5}}, \frac{7}{N^{2/5}}\right ) $区域做对照, 显然$ s=40 $时, 产生的PPC折线图有更多的点落在区域外, 该折线图实际收敛区域较宽.

在任意取定$ s $之后(本文取$ s=20 $), 接下来选择合适的收敛判别区域$ \left ( -\frac{c}{N^{\alpha }}, \frac{c}{N^{\alpha }} \right ) $, 即参数$ \alpha, c $. 不同发生器产生的PPC折线图的收敛速度不同, 而且真正的随机数序列产生的PPC折线图的收敛速度不能被确切地计算出来, 同时真随机数序列难以获得, 从而导致在零假设情况下的理论收敛速度难以得到. 所以本文利用一些广泛应用的PRNG生成类似于真随机数序列的伪随机数序列, 然后产生PPC折线图, 再通过聚类分析剔除PPC序列中的异常值, 尽可能地消除偶然出现的异常值对下一步的拟合产生的消极影响, 最后利用幂函数拟合的方式确定合适的收敛判别区域. 本文综合多种发生器的实际情况选取了几个最优的判别区域. 具体操作方法如下.

(1) 由PRNG产生PPC折线图.

由PRNG生成若干个相同长度的伪随机数序列, 则可以产生若干个$ \left\{ R_{N}\right\}_{N=k, 2k, \cdots , mk} $序列. 本文中每个发生器都产生100个长为$ 2^{15} $的伪随机数序列, 在选取判别标准时, 取$ k=2^{7} $, 则$ m=256 $, 可以得到100个长度为256的$ \left\{ R_{N}\right\}_{N=k, 2k, \cdots , 128k} $序列, 构造矩阵

(3.2) $ \begin{eqnarray} R_{100\times m}&=&\left[\begin{array}{ccccc} \left\{ R_{N}^{1}\right\}_{N=k, 2k, \cdots, mk}\\ \left\{ R_{N}^{2}\right\}_{N=k, 2k, \cdots, mk}\\ \vdots\\ \left\{ R_{N}^{100}\right\}_{N=k, 2k, \cdots, mk} \end{array}\right] =\left[\begin{array}{ccccc} R_{k}^{1} &R_{2k}^{1} &\cdots &R_{mk}^{1}\\ R_{k}^{2} &R_{2k}^{2} &\cdots &R_{mk}^{2}\\ \vdots &\vdots & \ddots &\vdots \\ R_{k}^{100} &R_{2k}^{100} &\cdots &R_{mk}^{100} \end{array}\right] _{100\times m}{}\\ &=&\left[\begin{array}{ccccc} \left\{ R_{k}^{i}\right\} &\left\{ R_{2k}^{i}\right\} &\cdots &\left\{ R_{mk}^{i}\right\} \end{array}\right] _{i=1, 2, \cdots, 100}, \end{eqnarray} $

其中$ \left\{ R_{N}^{i}\right\}_{N=k, 2k, \cdots, mk} $表示给定发生器生成的第$ i $个伪随机数序列产生的PPC序列, $ R_{jk}^{i} $表示$ \left\{ R_{N}^{i}\right\}_{N=k, 2k, \cdots, mk} $序列中的第$ j $个元素, $ \big\{ R_{jk}^{i}\big\}_{i=1, 2, \cdots, 100} $表示当$ N=jk $时, 100个实数序列产生的$ R_{jk}^{i} $值组成的序列.

(2) 处理$ R_{100\times m} $中的数据, 剔除异常值.

图 1显示由生成的伪随机数序列产生的折线图振动幅度初始时较大后来逐渐变小, 且PPC折线图在$ N $接近10000时基本收敛到0, 此时的振动幅度非常小, 因此可知PPC折线图的前一部分对收敛判别区域的选取影响较大. 所以主要对所有PPC序列的前一部分, 即子序列$ \left\{ R_{N}\right\}_{N=k, 2k, \cdots , hk} $进行分析, 其中$ h $为满足$ 1< h\leqslant m $, $ hk>10000 $的一个非负整数. 本文在选取判别标准时取$ h=80 $, 分析100个长度为80的$ \left\{ R_{N}\right\}_{N=k, 2k, \cdots , 80k} $序列, 即分析$ \left [ \left\{ R_{k}^{i}\right\}, \left\{ R_{2k}^{i}\right\}, \cdots, \left\{ R_{80k}^{i}\right\} \right ]_{i=1, 2, \cdots, 100} $中的数据. 将矩阵$ R_{100\times m} $的每一列数据看作一组, 则当$ N=jk(j=1, 2, \cdots , 80) $时, 对应着$ \big\{ R_{jk}^{i}\big\}_{i=1, 2, \cdots , 100} $, 因此可得相应的80组数据, 利用SPSS的系统聚类法分别对每一组数据进行分析, 剔除每组中少量的特殊点(尤其针对异常的较大值和较小值), 将较为异常的数值变为0.

(3) 利用幂函数进行拟合.

接下来取修改后每组数据$ \big\{ R_{jk}^{i}\big\}_{i=1, 2, \cdots , 100} $的最大值和最小值, 得到最大值序列$ \left\{ M_{jk}\right\}_{j=1, 2, \cdots, 80} $和最小值序列$ \left\{ m_{jk}\right\}_{j=1, 2, \cdots, 80} $, 然后以$ N(=k, 2k, \cdots , 80k) $为自变量进行幂函数形式的曲线拟合. 每个发生器都会通过拟合得到不同的幂函数表达式

$ y=\pm c/n^{\alpha }, $

即不同的$ (\alpha, c) $.

(4) 综合所有发生器的拟合结果, 选取合适的参数$ \alpha, c $.

对于收敛速度$ \alpha $, 本文选取了三个大多数发生器的拟合结果都接近的有理数$ \frac{3}{7}, \frac{2}{5}, \frac{4}{9} $, 相对应收敛判别区域的宽度$ c $分别为10, 7, 9. 图 4给出了当$ N=k, 2k, \cdots , 80k $时, 由LCG生成的100个伪随机数序列产生的PPC折线图和选取的三个不同的收敛判别区域.

图 4中两条黑色线表示收敛区域$ \left ( -10/N^{3/7}, 10/N^{3/7} \right ) $, 两条蓝色线表示收敛区域$ \left ( -7/N^{2/5}, 7/N^{2/5} \right ) $, 两条绿色线表示收敛区域$ \left ( -9/N^{4/9}, 9/N^{4/9} \right ) $. 这三个收敛区域在靠近0的地方差异较大, 显然两条黑色线代表的收敛区域较宽, 绿色线代表的收敛区域收敛速度最快, 到后面三个收敛域基本上重合, 差异较小. 并且可以看到PPC折线图上只有有限个点落在三个收敛域之外, 这表明选取的三个收敛域都是可行的.

关于阈值$ \beta $的选取, 通常取0.95或0.9. 收敛判别标准的参数$ \alpha, c, \beta $以及参数$ s $之间的关系如下.

(1) 固定$ \alpha $时, 参数$ s $越大, 收敛域宽度$ c $越大.

(2) 固定$ \alpha, c $时, 阈值$ \beta $越大, 判别标准越严, PPC测试也就越严格; 反之越松.

(3) 固定$ \alpha, \beta $时, 收敛域宽度$ c $越大, 判别标准越松, PPC测试也就越松; 反之越严.

本文认为可以通过PPC统计量$ p $值的大小来衡量伪随机数序列以及伪随机数发生器的随机性. 下面利用PPC测试方法对PRNG的性能进行测试和比较, 此时测试结果更侧重于发生器性能的比较, 因此可以忽略阈值$ \beta $的选取.

本节展示一些广泛使用的PRNG的测试结果, 这些PRNG是用Matlab编程语言实现的. 我们将PPC测试应用于不同的PRNG, 以及由$ \sqrt{n} $的小数部分生成的数字序列, 每个发生器都是用来自{http://www.random.org}的随机种子初始化的, 生成100组长为$ 2^{15} $的伪随机数序列, $ \sqrt{n} $生成一个长为$ 100\times 2^{15} $的实数序列, 分成100个长为$ 2^{15} $子序列, 每一位生成的伪随机数都保留8位小数. 对PRNG生成的伪随机数序列进行PPC测试, 计算每一条PPC折线分别落在三个不同的收敛判别区域$ \left ( -10/N^{3/7}, 10/N^{3/7} \right ) $、$ \left ( -7/N^{2/5}, 7/N^{2/5} \right ) $以及$ \left ( -9/N^{4/9}, 9/N^{4/9} \right ) $内的点的占比$ p_{1} $、$ p_{2} $和$ p_{3} $, 然后对得到100组测试结果$ \left ( p_{1}, p_{2}, p_{3} \right ) $求均值, 相当于计算由给定发生器生成的总PPC折线图分别落在三个不同的收敛域内的点占总数的百分比$ P_{1} $、$ P_{2} $以及$ P_{3} $.

具体的PPC测试PRNG的实验步骤如下.

(1) 由给定的PRNG生成100组长为$ 2^{15} $的伪随机数序列;

(2) 确定参数$ k=2^{7}, s=20 $, 根据(2.5)式和(3.1)式产生100个长为256的$ \left\{ R_{N}\right\}_{N=k, 2k, \cdots , mk} $序列, 得到相应的PPC折线图;

(3) 选择合适的收敛判别标准, 计算PPC折线图落在判别区域内的点的占比$ P $.

表 1给出了每个发生器生成的100个伪随机数序列产生的PPC序列落在区域$ \left ( -\frac{10}{N^{3/7}}, \frac{10}{N^{3/7}} \right ) $内的占比$ p_{1} $的均值, 即$ P_{1} $, 以及标准差和变异系数.

从表 1中可以看出, 前9种发生器生成的100组伪随机数序列的PPC测试结果的标准差和变异系数都小于0.1. 变异系数较小说明同一个发生器生成的伪随机数序列的PPC测试结果$ p $值离散程度较小, 用多组伪随机数序列的PPC测试结果的均值$ P $来衡量发生器的随机性, 结果可靠性较高.

基于无理数$ \pi $重叠产生的PRNG生成序列的PPC测试结果变异系数较大, 相比较其他发生器, 生成的序列随机性差异较大. 而且基于无理数$ \pi $重叠产生的PRNG的最终测试结果$ P_{1}=0.9453 $, 明显小于其他9种发生器的测试结果, 更加体现了该发生器性能较差. 由$ \sqrt{n} $的小数部分生成的数字序列的测试结果. $ P_1=0.1478 $, 与预期结果相符, 远远小于其他PRNG的测试结果. 图 5给出由$ \sqrt{n} $的小数部分生成的数字序列产生的PPC折线图, 显然该图与PRNG生成的伪随机数序列产生的PPC折线图有较大的差异.

本文选用基本的一级测试方法卡方检验和序列检验来测试生成序列的均匀性, 游程检验和自相关检验来测试伪随机数序列的独立性, 表 2显示了每个发生器生成的的100个实数序列通过基本测试的比例以及在三个不同判别标准下的PPC测试结果.

由表 2可知前面几种发生器的PPC测试结果相差不大, 为了判断这些发生器的PPC测试结果是否有明显的差异, 对前9种发生器的PPC测试结果进行单因素方差分析. 结果显示在0.01的显著性水平下, 所研究的PRNG性能之间不存在显著差异. 但是这9个发生器的性能明显优于基于无理数$ \pi $重叠产生的PRNG. 表 2中的结果显示在显著性水平取0.01的情况下, 除了基于无理数$ \pi $重叠产生的PRNG和由$ \sqrt{n} $的小数部分生成的数字序列, 其余的发生器生成的伪随机数序列几乎全部通过了这四种检验均匀性和独立性的一级测试方法, 而基于无理数$ \pi $重叠产生的PRNG所生成的序列均未通过序列检验和游程检验, 这说明它在均匀性和独立性方面均有缺陷, 由$ \sqrt{n} $的小数部分生成的数字序列也只通过了卡方检验, 这和PPC测试结果相吻合. 这表明了基于泊松对相关的伪随机数发生器的PPC测试能够同时检验生成序列的均匀性和独立性, 可以用于揭示随机数序列中的一些缺陷.

总之, 当生成的伪随机数序列的长度较短时(本文只检验了长度为$ 2^{15} $的实数序列), 无论是本文的测试方法还是现有的随机性测试方法, 测试结果都显示LCG、RANDU、MT19937-32、SPSS.RV.UNIFORM函数、多重递归发生器、BBS发生器、Matlab.rand函数、Python.random库中random函数和Excel.RAND函数这些发生器性能接近, 这个结果和文献[7]对LCG、RANDU、MT的测试结果一致, 文献[7]认为无论是它所建议的可视化方法还是现有的统计测试方法都显示LCG、RANDU、MT等数值PRNG方法的性能接近. 表 2中的PPC测试结果显示表中LCG、RANDU以及MT19937-32等前9种发生器的性能相较于基于无理数$ \pi $重叠产生的PRNG有一定的优势.

通过对比三个不同的收敛判别区域下各发生器的PPC测试结果, 可以看出阈值$ \beta $一定时, 收敛判别准则中当$ \alpha=4/9, c=9 $时, PPC测试相较于其他两组参数更加严格. 测试结果充分说明了本文提出的基于泊松对相关的PPC测试能够同时检验生成的伪随机数序列的均匀性和独立性, 检验能力与卡方检验、序列检验、游程检验以及自相关检验四种检验方法的综合能力相当, 甚至在均匀性和独立性方面, 还可以检测出这四种检验方法不能检测出来的缺点. PPC测试有较强的理论基础, 可操作性强, 对任意长度的伪随机数序列都可以测试, 其中参数$ k, s $可以灵活选取, 而且还可以根据实际需要选取合适的或严格或宽松的收敛判别准则, 且测试结果可靠.

PRNG在数值模拟计算^[26]、密码学^[27]等许多方面都起到非常重要的作用. 本文提出了一种基于泊松对相关的PRNG的测试方法, 该方法灵活稳健, 可操作性强, 并且具有强有效的理论支撑, 作为一个测试$ (0, 1) $上实数序列的一级测试方法, 主要是通过比较落在判别区域内的点的占比来衡量发生器的性能. 该测试方法的重点在于收敛判别标准的选取, 即参数$ (\alpha, c) $的确定, 本文利用SPSS系统聚类法去掉每组数据中的异常点, 用幂函数拟合的方式选取了三组不同的参数. 测试结果表明, PPC测试方法可以用来测试许多常见的PRNG, 然而具有良好随机性的经典PRNG仍然难以区分. 将PPC测试结果与卡方检验、序列检验、游程检验以及自相关检验的检验结果相比较, 结果显示该测试方法可以同时检验生成序列中独立性和均匀性方面的缺陷, 所以PPC测试方法作为一级测试是一个非常重要的工具. 但是像其他统计测试一样, PPC测试并不能捕捉所有已知的缺陷. 因此, 该测试方法应该与其他测试方法一起使用, 以便更严格地评估伪随机发生器.

此外不同的PRNG收敛速度都会有一定的差异, 文献[12] 认为实数序列产生PPC折线图的收敛速度可以用来衡量序列的均匀性和独立性, 所以在后续的工作中将继续研究参数s, 收敛速度以及收敛判别区域之间的关系, 对判别标准做出调整, 使得该测试方法更加严格, 更容易揭示生成序列以及发生器的一些缺陷, 有助于设计更稳健可靠的测试套件, 进而评估由PRNG生成的伪随机数序列的质量.