上帝到底掷不掷骰子？

"random"（随机）这个词看似简单易懂，但若深究，问题就会上升到量子力学和哲学的高度。

爱因斯坦说出的那句名言： "God does NOT play dice with the Universe!" 一直被人们视为他否定量子力学的证据（虽然事实好像并非如此），而量子力学把随机性看做物理世界的内在性质。

很多人持有这样的观点，即宇宙中所发生的一切从宇宙诞生的那一刻（The Bing Bang）就已经确定下来，所有发生的事件都是确定性的，而人类之所以无法预测事件的发生是因为人类自身的认知水平的限制。所以问题的焦点就在于一件事情的发生到底是 indeterminate（不可确定的）还是indeterminable （不能确定的）。举个简单的例子，对于轮盘游戏来说，当你得到所有变量(比如出手的力量，速度，轮盘的摩擦系数，空气的阻力，等等）的时候，轮盘转过的总角度就变成一个可计算的确定性结果，这时候就已经没有了随机性可言。这就能解释为什么“经典物理中不存在随机性”的观点能够得到普遍接受。

我们既不是理论物理学家，也并非哲学家。很多时候我们只需要有方法能够产生足够“随机”的数字供我们使用就可以。彩票销售，抽奖活动，密码学加密，模型的验证等等，我们都需要随机数来帮助我们实现想要的功能。

人类对随机数生成器的研究可以追溯到很久远的历史，发展到现在大致可以分为两大类研究方向：pseudorandom number generator（伪随机数生成器），True random number generator（“真”随机数生成器）

一、伪随机数

顾名思义就是假的随机数，所有的伪随机数都是用确定性的算法得到的，言外之意就是，可以通过某种算法计算出接下来要生成什么数字。所以一般的伪随机数算法都不是加密安全的。

伪随机数的生成算法有很多，但是很多算法到现在很多都不实用了，比如线性同余等许多同余类的算法，所以再去研究没有太大意义，目前使用最广泛的是Mersenne Twister（马特赛特旋转演算法），这完全得益于它的长周期以及弥补了经典算法的很多不足之处。

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html

1.1Mersenne Twister算法的基本过程

Mersenne Twister算法是利用线性反馈移位寄存器（LFSR）产生随机数的，n级的反馈移位寄存器理论上来说最大周期为2^n-1，也就是说理论上可以产生[0,2^n-1]这个范围的随机数。Mersenne Twister算法的n为19937，其周期为梅森素数2^19937-1，这也是该算法名称的由来。

当然计算机的位数最大仅为64位，所以代码实现的方法就是构造一个长度为19937/位数的数组。

需要注意的是，该算法实际上生成的是[0,2^32-1]范围内的随机数，长周期更好地保证了随机性。

关于线性反馈移位寄存器的简单介绍：https://blog.csdn/ACdreamers/article/details/44656743

该算法可以分成三大步（更详细介绍参看：https://blog.csdn/tick_tock97/article/details/78657851）

第一阶段：初始化，获得基础的梅森旋转链； 
第二阶段：对于旋转链进行旋转算法； 
第三阶段：对于旋转算法所得的结果进行处理；

算法中用到的变量如下所示：
    ·w：长度（以bit为单位）
    ·n：递归长度
    ·m：周期参数，用作第三阶段的偏移量
    ·r：低位掩码/低位要提取的位数
    ·a：旋转矩阵的参数
    ·f：初始化梅森旋转链所需参数
    ·b,c：TGFSR的掩码
    ·s,t：TGFSR的位移量
    ·u,d,l：额外梅森旋转所需的掩码和位移量


MT19937-32的参数列表如下：
    ·(w, n, m, r) = (32, 624, 397, 31)
    ·a = 9908B0DF（16）
    ·f = 1812433253
    ·(u, d) = (11, FFFFFFFF16)
    ·(s, b) = (7, 9D2C568016)
    ·(t, c) = (15, EFC6000016)
    ·l = 18

初始化：（以32位为例）

首先将传入的seed赋给MT[0]作为初值，然后根据递推式：MT[i] = f × (MT[i-1] ⊕ (MT[i-1] >> (32-2))) + i递推求出梅森旋转链。

MT19937规定了算法1中的f值为1812433253 ，算法将输入的随机种子赋值给第0个状态，剩余的623个状态则用前一状态值的一系列操作进行赋值更新。当执行完这一算法后，全部的624个状态已经填充完毕，之后梅森算法便可以不断地生成伪随机数

旋转生成随机数：

遍历旋转链，对每个MT[i]，根据递推式：MT[i] = MT[i+m]⊕((upper_mask(MT[i]) || lower_mask(MT[i+1]))A）进行旋转链处理。

其中，“||”代表连接的意思，即组合MT[i]的高 w-r （32-31=1）位和MT[i+1]的低 r （31）位，设组合后的数字为x，则xA的运算规则为（x0是最低位）： 

注：从图上可以看到integer0的后31位没有用到，所以一次运算用到的位总数为624*（1）+623*31=19937。

第三阶段：对于旋转算法所得的结果进行处理

b = 0x9d2c5680 
c = 0xefc60000 
y = mt[i] 
y = y XOR (y >> 11) 
y = y XOR ((y << 7) AND b) 
y = y XOR ((y << 15) AND c) 
y = y XOR (y >> 18)

 

1.2Mersenne Twister算法的相关定义简介

伪随机数发生器质量的度量——k-维 v-比特准确度（原文：https://liam.page/2018/01/12/Mersenne-twister/）

这里介绍评价 PRNG 最严格指标[Tootill et al. 1973][Fushimi and Tezuka 1983][Couture et al. 1993][Tezuka 1995][Tezuka 1994][Tezuka and L’Ecuyer 1991][L’Ecuyer 1996]：k-维 v-比特准确度（k-distributed to v-bit accuracy）

以梅森旋转算法为例简单讲解一下：该算法产生的是周期为P=2^19927-1，w=32比特的随机序列，现在假设用该算法产生了P个随机数组成的序列{Xi}，i的取值范围为[0,P).   已经知道该算法是623维32位比特准确的。也就是说对任何k<623,v<32的选取，都是k-维v-比特准确的。   假设k=3，v=2，即取每相邻的三个数（如X0,X1,X2）的最高两个比特位组成一个新的二进制数（如：010011），则这个二进制数的可能的取值有2^(6)=64中可能性，按此方法遍历{Xi}，生成P+1个这样的6比特二进制数，则每种可能的取值在所有P+1个数中所占的比例为1/64.

1.3Mersenne Twister算法代码（官网最新修改版本）

/*
A C-program for MT19937, with initialization improved 2002/1/26.
Coded by Takuji Nishimura and Makoto Matsumoto.

Before using, initialize the state by using init_genrand(seed)
or init_by_array(init_key, key_length).

Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:

1. Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.

3. The names of its contributors may not be used to endorse or promote
products derived from this software without specific prior written
permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


Any feedback is very welcome.
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
*/

#include <stdio.h>

/* Period parameters */
#define N 624
#define M 397
#define MATRIX_A 0x9908b0dfUL   /* constant vector a */
#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
#define LOWER_MASK 0x7fffffffUL /* least significant r bits */

static unsigned long mt[N]; /* the array for the state vector  */
static int mti = N + 1; /* mti==N+1 means mt[N] is not initialized */

						/* initializes mt[N] with a seed */
void init_genrand(unsigned long s)
{
	mt[0] = s & 0xffffffffUL;
	for (mti = 1; mti<N; mti++) {
		mt[mti] =
			(1812433253UL * (mt[mti - 1] ^ (mt[mti - 1] >> 30)) + mti);
		/* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
		/* In the previous versions, MSBs of the seed affect   */
		/* only MSBs of the array mt[].                        */
		/* 2002/01/09 modified by Makoto Matsumoto             */
		mt[mti] &= 0xffffffffUL;
		/* for >32 bit machines */
	}
}

/* initialize by an array with array-length */
/* init_key is the array for initializing keys */
/* key_length is its length */
/* slight change for C++, 2004/2/26 */
void init_by_array(unsigned long init_key[], int key_length)
{
	int i, j, k;
	init_genrand(19650218UL);
	i = 1; j = 0;
	k = (N>key_length ? N : key_length);
	for (; k; k--) {
		mt[i] = (mt[i] ^ ((mt[i - 1] ^ (mt[i - 1] >> 30)) * 1664525UL))
			+ init_key[j] + j; /* non linear */
		mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
		i++; j++;
		if (i >= N) { mt[0] = mt[N - 1]; i = 1; }
		if (j >= key_length) j = 0;
	}
	for (k = N - 1; k; k--) {
		mt[i] = (mt[i] ^ ((mt[i - 1] ^ (mt[i - 1] >> 30)) * 1566083941UL))
			- i; /* non linear */
		mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
		i++;
		if (i >= N) { mt[0] = mt[N - 1]; i = 1; }
	}

	mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
}

/* generates a random number on [0,0xffffffff]-interval */
unsigned long genrand_int32(void)
{
	unsigned long y;
	static unsigned long mag01[2] = { 0x0UL, MATRIX_A };
	/* mag01[x] = x * MATRIX_A  for x=0,1 */

	if (mti >= N) { /* generate N words at one time */
		int kk;

		if (mti == N + 1)   /* if init_genrand() has not been called, */
			init_genrand(5489UL); /* a default initial seed is used */

		for (kk = 0; kk<N - M; kk++) {
			y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
			mt[kk] = mt[kk + M] ^ (y >> 1) ^ mag01[y & 0x1UL];
		}
		for (; kk<N - 1; kk++) {
			y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
			mt[kk] = mt[kk + (M - N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
		}
		y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
		mt[N - 1] = mt[M - 1] ^ (y >> 1) ^ mag01[y & 0x1UL];

		mti = 0;
	}

	y = mt[mti++];

	/* Tempering */
	y ^= (y >> 11);
	y ^= (y << 7) & 0x9d2c5680UL;
	y ^= (y << 15) & 0xefc60000UL;
	y ^= (y >> 18);

	return y;
}

/* generates a random number on [0,0x7fffffff]-interval */
long genrand_int31(void)
{
	return (long)(genrand_int32() >> 1);
}

/* generates a random number on [0,1]-real-interval */
double genrand_real1(void)
{
	return genrand_int32()*(1.0 / 4294967295.0);
	/* divided by 2^32-1 */
}

/* generates a random number on [0,1)-real-interval */
double genrand_real2(void)
{
	return genrand_int32()*(1.0 / 4294967296.0);
	/* divided by 2^32 */
}

/* generates a random number on (0,1)-real-interval */
double genrand_real3(void)
{
	return (((double)genrand_int32()) + 0.5)*(1.0 / 4294967296.0);
	/* divided by 2^32 */
}

/* generates a random number on [0,1) with 53-bit resolution*/
double genrand_res53(void)
{
	unsigned long a = genrand_int32() >> 5, b = genrand_int32() >> 6;
	return(a*67108864.0 + b)*(1.0 / 9007199254740992.0);
}
/* These real versions are due to Isaku Wada, 2002/01/09 added */

int main(void)
{
	int i;
	unsigned long init[4] = { 0x123, 0x234, 0x345, 0x456 }, length = 4;
	init_by_array(init, length);
	printf("1000 outputs of genrand_int32()\n");
	for (i = 0; i<1000; i++) {
		printf("%10lu ", genrand_int32());
		if (i % 5 == 4) printf("\n");
	}
	printf("\n1000 outputs of genrand_real2()\n");
	for (i = 0; i<1000; i++) {
		printf("%10.8f ", genrand_real2());
		if (i % 5 == 4) printf("\n");
	}
	return 0;
}

 

MATLAB中有一页主要介绍其使用的各种伪随机数生成算法

MATLAB随机数生成算法介绍

对随机数颇有价值的研究方向是对衡量随机性和检测随机数的方法的研究，这方面相关论文有很多。

需要注意的是，尽管Mersenne Twister算法以及各种伪随机数算法都有广泛的应用，但其方法本质是用有限的状态来近似无限的取值范围，从本质上来说都不是真正的随机数，但是在一般的应用过程中不需要那么高的精度，所以这种伪随机数算法在实际应用中是可行的

二、“真”随机数

单凭计算机的算法是无法得到意义上的“随机数”的，必须借助外部的物理现象来得到。通过采集外部噪声（比如大气噪声，量子波动等等）的方法制作真随机数生成器

有关硬件随机数的介绍请参看wiki词条硬件随机数生成器

硬件随机数生成器种类很多，比如量子随机数生成器，具体的产品很多并且广泛的用于信息安全领域来实现加密需求,当你在google中搜索quantum random number generator 会出现各种各样的产品以及对应的论文。

1.https://www.random

2.http://qrng.anu.edu.au

3.http://www.fourmilab.ch/hotbits/

以上三个都是通过物理实验的方法得到随机数的网站，可以通过提供的API等获取实验室时时产生的随机数。

2.1真随机数的生成过程

在概率论中，对于任何一个连续的随机变量X，设其密度函数为p(x),分布函数为F(x)，令随机变量Y=F(x)，则根据随机变量函数的相关原理和计算公式可知随机变量Y一定服从U(0,1),即在区间（0,1）上的均匀分布。

根据上面的定理，从理论上我们就可以根据均匀分布的随机数 y 来生成服从任意分布的随机数 x ，其计算公式为,其中为随机变量X的分布函数的反函数

但事实上，这种方法在实际应用中有诸多困难，首先，即使随机变量x的密度函数p(x)为初等函数，其分布函数F(x)也不一定为初等函数，通常无法写出其具体的表达式。其次，即使写出其表达式，对于复杂的函数来说求其反函数是一件很困难的事情。所以在一般的应用中，你所见到的大多是正态分布和均与分布以及一些其他的简单的离散分布。

暂且抛开上面不谈，我们来看一下真随机数生成器是如何产生随机数的。

上面三个网址分别介绍了三个不同的随机数生成器，虽然他们的原理各有一些不同（大气噪声，量子物理），但他们产生的随机数都是0,1二元序列，也就是我们所说的比特流（例如“00010101010010100101”）

而从上面的概率论的原理分析来看，U(0,1)也就是区间（0,1）上的均匀分布是核心，得到了（0,1）均匀分布的随机变量，我们就可以通过逆累积函数，或者Box–Muller transform来得到标准正态分布随机数，进而得到服从其他分布的随机数。

借鉴伪随机数生成中所用到的方法，可以通过规定周期M来生成在某种精度下的U(0,1)的随机数。

比如我们取定每个序列的长度为10,则一个随机序列块所能表示的最大随机数为 1111111111(2)=2^(10)-1=1023.将随机数生成器生成的长串01序列按每10个为一块进行分割。则每一块就代表了一个随机数字ai，然后用ai除以1023，就可以得到开区间（0,1）上的小数ai/1023.且每个数字的概率都为1/1023(约等于0.001),也就是说可以近似的达到千分之一的精度。

或者为了提高精度，选定更长的块长度，但这样带来的问题就是01序列的生成速度是有限的，所以真随机数生成器的效率往往没有伪随机生成器的效率高。可以采取每次只移动一位来方法来生成下一个随机数，例如序列X=“010101010010100101001”长度为21，若以20为块长度，则可以用后20位为第一个随机数，以前20位为第二个随机数。

 

2.2数据处理

利用像量子随机数生成器这样的设备生成随机数也存再一定的缺点，这些通过外部设备得到的随机数数据不可避免的会带有一定的的噪声，所以原始的数据必须要经过一定的处理削减其自相关性和偏差的存在。

三、随机性测试

The problems with randomness is you can never be sure

对于随机数的研究，最重要并不是如何产生随机数，而是如何测试生成器所生成的数字是否真的随机。这并不像听上去那么简单，证明一个随机数序列真的随机是一件非常困难的事情，因为即使是真的随机数序列也会存在测试不通过的情况，而测试通过也并不意味着所产生的序列一定是随机的而并非精巧的算法构造出来的。尽管如此，这些测试对于衡量序列的随机性依然必不可少。

关于随机性的测试有很多，目前主要用到的比如NIST（美国国家标准技术研究所）的标准：

A Statistical Test Suite for Random and Pseudorandom Number Generators for Cryptographic Applications

以及我国的国家密码管理局公布的标准规范：

国密规范   随机性检测规范

当然测试方法远不止上面所列的规范中的方法，其他方法可以在google学术、百度学术等论文查询网站进行查找

目前普遍的测试软件主要有diehard，TestU01，nist sp800-22

https://en.wikipedia/wiki/Diehard_tests

http://simul.iro.umontreal.ca/testu01/tu01.html

https://csrc.nist.gov/Projects/Random-Bit-Generation/Documentation-and-Software

需要注意的是，规范中的很多方法都是用于检测0,1二进制序列（如：000101010010这样的随机序列），这和生成随机数的方法和数学原理有关。对于其他分布的随机性检测还有很多的方法，可以参看相关文献查找。

3.1 TestU01的安装与使用（windows）

官网的安装教程网址：http://simul.iro.umontreal.ca/testu01/install.html

注：你可以按照官网的教程，安装msys，然后在Bourne shell中运行，这样做确实简单，但这样每次运行的操作流程很麻烦，所以我打算用其他的方法来实现

第一步：下载mingw：http://www.mingw/，在安装选项里面可以选择安装msys，这样会稍微简单一点。如果不安装的话通过其他方法也可以。TestU01需要GNU C编译器来编译，所以只要电脑里面装了gcc或其他类似的编译器就可以了。如果你电脑里面安装有CodeBlocks等带有gcc编译器的IDE，那你完全可以跳过这一步

第二步：下载：Binaries for MinGW under MS Windows，里面有一个usr文件夹，内容如下：

lib里面是一些库文件，include里面是各种算法的头文件，share里面有参数文件夹param（各种算法的参数设置），和各种例子examples

下面开始创建工程，运行一个例子，我的电脑里面用的是CodeBlocks，所以以这个为例

第三步：打开CodeBlocks，创建一个新工程，将include，lib，examples，param文件夹拷贝到工程文件夹下，然后在gcc的lib文件夹里面找到libwsock32.a，拷贝到该工程的lib文件夹里面

第四步：打开一个examples中的文件，比如birth1

第五步：设置工程的build option

然后在linker setting下的link libraries里面add工程文件夹下的lib文件夹中的.a文件，（顺序一定不能错，test，prodist，mylib，wsock的顺序不能变）

然后search directories的compiler选项 add  include文件夹，linker add lib文件夹，点击ok，搞定

第六步：编译，运行