本文篇幅较长，如果想直接看js实现可定位ECMAScriptJavaScript。

（github：https://github.com/MichealWayne，个人博客地址：https://blog.michealwayne.cn/）

无论开发哪种编程语言，都可能满足随机数的需求。当然，我们会在开发过程中调用它API实现方法，如js的Math.random()、Python的random()，但实现原则也是一个值得深入研究的问题。

首先，回顾大学概率论中与随机数相关的概念：

• 频率：n次实验是在相同的条件下定义的。在这个n次实验中，事件A发生的次数nA称为事件A的频率，比值nA/n被称为事件A的频率。
• 概率：设E是随机试验，S是它的样本空间，给E的每个事件A一个实数，记录为P(A)，称为事件概率。
• 随机变量：随机试验样本空间为S={e}，若X=X(e)将样本空间S上的实值单值函数定义为X=X(e)随机变量。
• 随机事件：可能发生或不发生在条件S下的事件事件称为随机事件。

确定分布的独立随机数序列：每个数字的出现只是偶然的，与该序列中的其他数字无关。每个数字都有确定的概率进入任何给定的范围。

随机事件是否发生在试验中是随机的，但随机性是规律的。随机事件A的概率是频率的稳定值，频率是概率的近似值。

像js中的Math.random()在概率分布理论中，它属于连续和离散的均匀分布。均匀分布的特点是：所有基本事件的可能性相等。

可以参考Federal Office for Information Security (德国联邦信息安全办公室)给出了随机数发生器质量评估的四个标准-BSI评估标准：

• K1 - 随机数序列不同的概率应该很高
• K2 - 根据一些统计测试，无法区分产生的序列和真实的随机序列。这些测试包括： monobit测试(序列中0和1的数量相等)，poker测试（ chi-squared特殊情况)、runs来自测试(不同长度运行的频率)BSI 和 NIST 的longruns测试、autocorrelation测试。
• K3 - 给定任何子序列，任何攻击者都无法计算后续序列或生成器的内部状态
• K4 - 任何攻击者都内部状态，任何攻击者都无法计算生成器前的序列或状态 只满足密码学的应用K3和K标准生成器是可以接受的。

确定性方法不能真正创造随机性，正如《计算机程序设计艺术（第二卷）》所说：今天使用的大多数随机数生成器都不够好，开发人员倾向于使用它们，而不了解具体的生成策略。因此，我们经常发现，一些有缺陷的旧随机数生成器将盲目地用于一个又一个程序，但没有人关心它们的局限性。

• 在早期阶段，那些在科学工作中需要随机数量的人通过在一个完全旋转的罐子中抓球或骰子和分牌来获得随机数量。

• 1927年，L.H.C.Tippett从人口统计调查报告中随机获得发布了4万多个随机数字的表格。从那以后，一些特殊的设备被制造出来机械地生成随机数。

• 1939年，M.G.Kendall和B.Babington-Smith使用第一台这样的机器创建了一个有1万个随机数字的表

• 1949年，D.H.Lehmer提出线性同余法（LCG）算法方案。

• 1951年首次安装Ferranti Mark Ⅰ计算机有一个内置指令，它使用电阻噪声生成器将20个随机位置放入累加器中A.M.Turing的推荐。

• 1955年，RAND该公司还发布了一张广泛使用的100万个随机数字表，该表是由另一个特殊设备获得的。一个叫做ERNIE(厄尼)著名的随机数机器已经使用多年，用于生成英国政府有奖债券的中奖号码。

• 计算机问世后不久，人们开始探索在计算机程序中随机数量的有效方法，但由于内存空间和输入时间的要求，该方法的实用性有限。随着20世纪90年代技术的发展，由于CDROM可存储)1GB测试的随机数。

• 1995年，George Marsaglia 通过扰频输出和确定噪声二极管的电流rap音乐叠加在一起生成650MB随机值(称为黑白噪声)并将其制成演示光盘，从而促进了随机数表的复兴。

• 2014年，计算机PRNG算法-PCG诞生时，它以更小、更快的代码和更小的状态实现了优异的统计性能。

早期机械生成随机数的缺点引起了人们对使用计算机的普通算术操作产生随机数的兴趣。1946年左右，Jogn von Neumann(冯诺伊曼)第一个建议是用这种方法。他的方法是取前面随机数的平方，取中间数。

例如，如果正在生成10位数字，且之前的值为577215649，则将其平方得到33317792380594909201，因此下一个值为7923805949。

对于这种技术，有十分明显的异议：既然每个数完全由它先前的数所决定，那么以这样的方法产生的序列怎么会是随机的呢？当然，这个序列不是随机的，只是随机的。在典型的应用中，一个数与跟在它后面的那个数之间的真实关系并无客观的意义，因此，这种非随机性的特征并不真是不可取的。直观地说，平方取中似乎完全扰乱了前面的数字。

已经证明，冯诺伊曼的方法不是寻求随机数量的好方法，危险在于该序列容易出现重复元素的短循环。例如，如果0作为该序列的数量出现，它将继续重现自己。

20世纪50年代初，有人用平方取中法进行了实验。G.E.Forsythe用4位数代替10位数，测试了16个不同的初始值，发现其中12个导致循环6100、2100、4100、8100、6100，…为结局序列，其中两个退化为0。N.Metropolis广泛测试了平方取中法，其中大部分是二进数系统。他证明，当使用20位数时，序列可能会退化为13个不同的循环，其中最长的循环周期为142。

1965年，根据I.J.Matrix博士说，数学家说π十进制作为随机序列，对于现代数值逻辑学家来说，它有着非常丰富的模式值得注意。

例如，Matrix医生指出，在π第一次重复的两位数是26，它的第二次出现在一个奇妙的重复模式中间

在列出许多数字或其他性质后，人们会发现，如果正确解释，π能反映人类经历的整个历史。

（其中X是一个10位的十进制数）

考虑到算法K的复杂设计，这一算法似乎能产生无穷尽的令人难以置信的随机数，但事实上，当把这个算法头一次放到计算机上时，它几乎立即收敛到10位数值6065038420——非常巧合，本算法把这个数转换成它自己。若以另外一个数开始，则这个序列在7401个值之后，开始以长度为3178的周期进行循环

这个事告诉我们，随机数并不是通过随机地选择的方法就可生成的，应当使用某些理论。

后文算法便是通过理论算法实现的，统称为PRNG（pseudo random number generator）伪随机数生成器，又被称为确定性随机比特生成器（deterministic random bit generator，DRBG）。其中适合于加密应用程序的PRNG称为加密安全的PRNG（CSPRNG）。CSPRNG的一个必要条件是不知道初始种子的敌人在分辨生成器的输出序列和真随机序列时只有可忽略的优势。换句话说，若PRNG只需要通过特定统计测试时，则CSPRNG必须通过种子规模的多项式复杂度内的所有统计测试

PRNG是一个生成数字序列的算法，其特性近似于随机数序列的特性。PRNG生成的序列并不是真随机，因此它完全由一个初始值决定，这个初始值被称为PRNG的随机种子（seed，但这个种子可能包含真随机数）。尽管接近于真随机的序列可以通过硬件随机数生成器生成，但伪随机数生成器因为其生成速度和可再现的优势，在实践中也很重要。

PRNG的周期定义为：所有初始值的最大长度的无重复前缀序列。周期受状态数的限制，通常用比特位数表示。然而，每增加一个比特位，周期长度就可能增加一倍，所以构建周期足够长的PRNG对于许多实际应用程序来说是很容易的。

当今使用的流行的随机数生成程序是D.H.Lehmer 1949年介绍过的下列方案的特殊情况。我们选择4个魔术整数：

m, 模*数； 0<m
a, 乘数；  0<=a<m
c, 增量；  0<=c<m
X0,开始值; 0<=X0<m

然后通过置公式（1）

Xn+1 = (aX0 + c) mod m,  n>=0

而得到所求的随机数序列<Xn>。这个序列称做线性同余序列。对m求余有点像确定转动的轮盘上球的落点。

例如，当m = 10和X0 = a = c = 7时，得到的序列是：

7,6,9,0,7,6,9,0...

如此例所示，对于m, a, c和X0的所有选择，这个序列并不总是“随机”的。这个例子有一个长度为4的周期，也说明了同余序列总是进入一个循环的事实。

改进：

定义b = a - 1，a>=2, b>=1，公式（2）：

Xn+k = (a^k * Xn + (a^k - 1) * c / b) mod m,  k>=0, n>=0

它借助于第n项来直接表达第(n + k)项，由此得出，由<Xn>的每第k项组成的子序列是另一个线性同余序列，该子序列具有乘数a^k mod m和增量((a^k - 1) * c / b) mod m。

这个公式的一个重要推论是：由m、a、c和X0定义的一般序列，能非常简单地借助c=1和X0=0的特殊情况来表示，设

Y0 = 0, Yn+1 = (a * Yn + 1) mod m

根据公式（2），我们得到Yk等于(a^k - 1) / b(modulo m)，因此公式（1）中定义的一般序列满足

Xn = (A * Yn + X0) mod m，其中A = (X0 * b + c) mod m

并且可以通过分析（详见《计算机程序设计艺术（第二卷）》3.2.1.1）得出m可以取2^31 - 1（2^32），乘数a选择在0.01m和0.99m之间，而且它的二进制表示或十进制表示数字不应有一个简单的正规模式，推荐值是（7^5，即16807）。

但是对于现代的计算机而言，这样的周期还是太短了。生成器的内部结构如栅栏结构和超平面点的分布也很重要。对于不同的生成器有特定的检测方法。结构检测用到最多的就是谱检验，谱检验就是基于相邻平行超平面之间最大距离的检验，该距离越大，生成器越差。

如图：

1991年，MWC（multiply-with-carry）由 George Marsaglia 和 Zaman 发明，非常类似于 LCG。MWC以其最简单的形式使用与线性同余生成器类似的公式，但是c（“加数”）在每次迭代中都不同。

其优势在于调用简单的计算机整数算法，使得非常快速地生成具有巨大周期的随机数序列；生成的循环周期更长（2^60 ~ 2^2000000），接近于 CPU 的循环周期。

它被广泛运用于游戏类的开发中，被非正式地称为“所有PRNG的母亲”，这个名字最初由Marsaglia自己创造。

以下是使用128位乘法并具有64位输出的小型MWC生成器的c语言实现：

// C99 + __uint128_t MWC, 256 bits of state, period approx. 2^255

/* The state must be neither all zero, nor x = y = z = 2^64 - 1, c = MWC_A3 - 1. The condition 0 < c < MWC_A3 - 1 is thus sufficient. */
   
uint64_t x, y, z, c = 1;

#define MWC_A3 0xff377e26f82da74a

uint64_t inline next() { 
        
    const __uint128_t t = MWC_A3 * (__uint128_t)x + c;
    x = y;
    y = z;
    c = t >> 64;
    return z = t;
}

一个高速、简便的MWC实现，其中js的V8引擎早期就是用了它（见后文）。

#define ulong unsigned long

static ulong mwc1616_x = 1;
static ulong mwc1616_y = 2;

void seed_rand_mwc1616(ulong seed) { 
        
    mwc1616_x = seed | 1;   /* Can't seed with 0 */
    mwc1616_y = seed | 2;
}

ulong rand_mwc1616(void) { 
        
    mwc1616_x = 18000 * (mwc1616_x & 0xffff) + (mwc1616_x >> 16);
    mwc1616_y = 30903 * (mwc1616_y & 0xffff) + (mwc1616_y >> 16);
    return (mwc1616_x<<16)+(mwc1616_y&0xffff);
}

但是MWC无法通过TestU01的多项测试。

置换同余生成器（PCG, permuted congruential generator）是LCG的改进，2014产生。它以更小的，快速的代码和小的状态量实现了出色的统计性能。 PCG与传统线性同余生成器在以下三个方面有所不同：

• LCG模数和状态较大，通常是所需输出大小的两倍；
• 它使用2的幂模，这导致使用全周期发生器和无偏输出位特别有效的实现
• 状态不是直接输出，而是状态的最高有效位用于选择按位旋转或移位，将其应用于状态以产生输出。

也正是可变的旋转，消除了2次幂的LCG的低阶比特短时间周期的问题。

PCG系列包括许多变体。虽然实际建议仅使用64位和128位，但核心LCG的定义宽度为8位至128位。如PCG-XSH-RR：

#include <stdint.h>
static uint64_t       state      = 0x4d595df4d0f33173;		// Or something seed-dependent
static uint64_t const multiplier = 6364136223846793005u;
static uint64_t const increment  = 1442695040888963407u;	// Or an arbitrary odd constant

static uint32_t rotr32(uint32_t x, unsigned r)
{ 
        
	return x >> r | x << (-r & 31);
}

uint32_t pcg32(void)
{ 
        
	uint64_t x = state;
	unsigned count = (unsigned)(x >> 59);		// 59 = 64 - 5

	state = x * multiplier + increment;
	x ^= x >> 18;								// 18 = (64 - 27)/2
	return rotr32((uint32_t)(x >> 27), count);	// 27 = 32 - 5
}

void pcg32_init(uint64_t seed)
{ 
        
	state = seed + increment;
	(void)pcg32();
}

Marsaglia建议通过将输出与简单的加法计数器模2^32（Weyl序列）组合，这会使周期增加2^32倍，达到2^192-2^32。

#include <stdint.h>
struct xorwow_state { 
        
  uint32_t a, b, c, d, e;
  uint32_t counter;
}
uint32_t xorwow (struct xorwow_state *state) { 
        
  uint32_t t = state->e;
  uint32_t s = state->a;
  state->e = state->d;
  state->d = state->c;
  state->c = state->b;
  state->b = s;
  t ^= t >> 2;
  t ^= t << 1;
  t ^= s ^ (s << 4);
  state->a = t;
  state->counter += 362437;
  return t + state->counter
}

该生成器的典型案例是 Nvidia CUDA 工具包。

xorshift*采用xorshift生成器，并对其输出应用可逆乘法（以字大小为模），作为非线性变换。

以下具有64位状态的64位生成器的最大周期为2^64-1，仅TestU01 BigCrush的MatrixRank测试失败：

#include <stdint.h>
struct xorshift64s_state { 
        
  uint64_t a;
}
uint64_t xorshift64s (struct xorshift64s_state *state) { 
        
  uint64_t x = state->a;
  x ^= x >> 12;
  x ^= x << 25;
  x ^= x >> 27;
  state->a = x;
  return x * UINT64_C(0x2545F4914F6CDD1D);
}

Vigna 建议使用下面的xorshift1024 *生成器，该生成器具有1024位状态且最大周期为2^1024-1：

#include <stdint.h>

/* The state must be seeded so that there is at least one non-zero element in array */
struct xorshift1024s_state { 
        
	uint64_t array[16];
	int index;
};

uint64_t xorshift1024s(struct xorshift1024s_state *state)
{ 
        
	int index = state->index;
	uint64_t const s = state->array[index++];
	uint64_t t = state->array[index &= 15];
	t ^= t << 31;		// a
	t ^= t >> 11;		// b
	t ^= s ^ (s >> 30);	// c
	state->array[index] = t;
	state->index = index;
	return t * (uint64_t)1181783497276652981;
}

除了使用乘法，还可以使用加法作为更快的非线性变换。这个想法最初是由Saito和Matsumoto（还负责Mersenne Twister）在XSadd生成器中提出的，该生成器基于32位移位将基础xorshift生成器的两个连续输出相加。

但是，其输出的低阶位有一些弱点；当输出字位反转时，它在几次BigCrush测试中失败。为了解决这个问题，Vigna 引入了基于64位移位的xorshift +系列：以下xorshift128 +生成器使用128位状态，最大周期为2^128-1。

#include <stdint.h>

struct xorshift128p_state { 
        
  uint64_t a, b;
};

/* The state must be seeded so that it is not all zero */
uint64_t xorshift128p(struct xorshift128p_state *state)
{ 
        
	uint64_t t = state->a;
	uint64_t const s = state->b;
	state->a = s;
	t ^= t << 23;		// a
	t ^= t >> 17;		// b
	t ^= s ^ (s >> 26);	// c
	state->b = t;
	return t + s;
}

此生成器是通过 BigCrush 测试最快的生成器之一，添加连续输出的一个缺点是，底层的xorshift128生成器是二维分布的，而关联的xorshift128 +生成器只有一维分布的。

xorshift+没有通过BigCrush的反向测试。

xoshiro和xoroshiro是移位寄存器生成器的其他变体，除了移位之外还使用旋转。根据Vigna的说法，与xorshift相比，它们速度更快，并产生更好的质量输出。 此类生成器具有32位和64位整数和浮点输出的变体。对于浮点数，取高53位（对于binary64）或高23位（对于binary32），因为高位的质量比浮点生成器中的低位更好。该算法还包括跳转功能，该功能将状态向前设置若干步（通常为2的幂），以允许许多执行线程从不同的初始状态开始。

xoshiro256**是该系列的通用随机64位数字生成器：

/* Adapted from the code included on Sebastian Vigna's website */

#include <stdint.h>

uint64_t rol64(uint64_t x, int k)
{ 
        
	return (x << k) | (x >> (64 - k));
}

struct xoshiro256ss_state { 
        
	uint64_t s[4];
};

uint64_t xoshiro256ss(struct xoshiro256ss_state *state)
{ 
        
	uint64_t *s = state->s;
	uint64_t const result = rol64(s[1] * 5, 7) * 9;
	uint64_t const t = s[1] << 17;

	s[2] ^= s[0];
	s[3] ^= s[1];
	s[1] ^= s[2];
	s[0] ^= s[3];

	s[2] ^= t;
	s[3] = rol64(s[3], 45);

	return result;
}

xoshiro256+ 比 xoshiro256** 快约15％，但最低的三位具有较低的线性复杂度；因此，应仅通过提取高53位数值并将其用于浮点结果。

#include <stdint.h>

uint64_t rol64(uint64_t x, int k)
{ 
        
	return (x << k) | (x >> (64 - k));
}

struct xoshiro256p_state { 
        
	uint64_t s[4];
};

uint64_t xoshiro256p(struct xoshiro256p_state *state)
{ 
        
	uint64_t (*s)[4] = &state->s;
	uint64_t const result = s[0] + s[3];
	uint64_t const t = s[1] << 17;

	s[2] ^= s[0];
	s[3] ^= s[1];
	s[1] ^= s[2];
	s[0] ^= s[3];

	s[2] ^= t;
	s[3] = rol64(s[3], 45);

	return result;
}

值得注意的是：如果空间有限，则xoroshiro128 **等同于xoshiro256 **，xoroshiro128 +等同于xoshiro256 +。它们具有较小的状态空间，因此对于大规模并行程序没有多大用处。 xoroshiro128 +还表现出对Hamming权重的轻微依赖，在输出5TB数据后会产生故障。 对于32位输出，xoshiro128 **和xoshiro128 +完全等同于xoshiro256 **和xoshiro256 +。

原则上，在设计随机算法的时候，我们会希望初始状态使用尽可能少的内存，执行速度更快，循环周期更长，且随机性质量越高。

理论检验主要是检验随机数的周期，内在结构和关联性。不过光理论检验是不够的，我们还需要经验检验。经验检验性对较简单，且有很多的方法。比如：等分布检验：又称均匀性检验，将[0,1)区间分成k个相等的子区间，落在每个子区间的伪随机数个数应该相等。常用的是x^2检验。序列检验：相继的数偶独立一致分布。计算的是数偶对的出现次数。间隔检验：用来考察在某个范围内序列的出现的间隔长度。扑克检验：考虑5个相继整数组成的n个组，观察其出现模式…

随机性的质量测试主要可以通过工具TestU01，它是一个软件库，以ANSI C语言实现，并提供了用于对统一随机数生成器进行经验统计测试的实用程序的集合。 TestU01直接测试的必须是一个能生成[0,1)之间浮点数的函数或一个能生成32位全随机无符号整数的函数。 它的整套测试被称作BigCrush，是目前最严格的随机数测试套件，它可以很好地检测一个随机数生成器的可靠程度。

如下图所示，只有7中RNG算法通过了测试，详细测试报告可见： https://www.pcg-random.org/pdf/toms-oneill-pcg-family-v1.02.pdf

*优秀的PRNG还有很多，像非常安全且周期极长的Mersenne Twiste（马特赛特旋转演算法），本文暂时先不记录。

js随机数的主要生成方式毫无疑问是通过Math.random()，首先我们来看官方文档（ES6）对这个方法的介绍：

简单翻译来说就是通过在0～1范围内以大致均匀的分布随机或伪随机生成一个0～1范围内的数字。以V8引擎为例，我们来看看Math.random()的具体实现：

直至v4.9.40版本，V8选择的是MWC1616算法，核心内容：

uint32_t state0 = 1;
uint32_t state1 = 2;
uint32_t mwc1616 () {
  state0 = 18030 * (state0 & 0xFFFF) + (state0 >> 16);
  state1 = 30903 * (state1 & 0xFFFF) + (state1 >> 16);
  return state0 << 16 + (state1 & 0xFFFF);
}

MWC1616的计算速度很快，且占用很少的内存，但是它的质量低于标准水平：

• 1.范围限制：这种算法的随机数范围在2^32以内，但是双精度浮点数能代表2^52个0～1的小数；

• 2.结果控制：结果值的上半部分非常依赖于state0；如果初始状态选择不正确，周期长度可能少于4亿；

• 3.无法通过TestU01的多项测试。

因此V8对随机数的算法进行了优化（v4.9.41.0开始，Chrome49开始）：

首先来看Math.random()的源码（v8/src/numbers/math-random.cc）：

// Copyright 2018 the V8 project authors. All rights reserved.
// Use of this source code is governed by a BSD-style license that can be
// found in the LICENSE file.

#include "src/numbers/math-random.h"

#include "src/base/utils/random-number-generator.h"
#include "src/common/assert-scope.h"
#include "src/execution/isolate.h"
#include "src/objects/contexts-inl.h"
#include "src/objects/fixed-array.h"
#include "src/objects/smi.h"

namespace v8 {
namespace internal {

void MathRandom::InitializeContext(Isolate* isolate,
                                   Handle<Context> native_context) {
  Handle<FixedDoubleArray> cache = Handle<FixedDoubleArray>::cast(
      isolate->factory()->NewFixedDoubleArray(kCacheSize));
  for (int i = 0; i < kCacheSize; i++) cache->set(i, 0);
  native_context->set_math_random_cache(*cache);
  Handle<PodArray<State>> pod =
      PodArray<State>::New(isolate, 1, AllocationType::kOld);
  native_context->set_math_random_state(*pod);
  ResetContext(*native_context);
}

void MathRandom::ResetContext(Context native_context) {
  native_context.set_math_random_index(Smi::zero());
  State state = {0, 0};
  PodArray<State>::cast(native_context.math_random_state()).set(0, state);
}

Address MathRandom::RefillCache(Isolate* isolate, Address raw_native_context) {
  Context native_context = Context::cast(Object(raw_native_context));
  DisallowHeapAllocation no_gc;
  PodArray<State> pod =
      PodArray<State>::cast(native_context.math_random_state());
  State state = pod.get(0);
  // Initialize state if not yet initialized. If a fixed random seed was
  // requested, use it to reset our state the first time a script asks for
  // random numbers in this context. This ensures the script sees a consistent
  // sequence.
  if (state.s0 == 0 && state.s1 == 0) {
    uint64_t seed;
    if (FLAG_random_seed != 0) {
      seed = FLAG_random_seed;
    } else {
      isolate->random_number_generator()->NextBytes(&seed, sizeof(seed));
    }
    state.s0 = base::RandomNumberGenerator::MurmurHash3(seed);
    state.s1 = base::RandomNumberGenerator::MurmurHash3(~seed);
    CHECK(state.s0 != 0 || state.s1 != 0);
  }

  FixedDoubleArray cache =
      FixedDoubleArray::cast(native_context.math_random_cache());
  // Create random numbers.
  for (int i = 0; i < kCacheSize; i++) {
    // Generate random numbers using xorshift128+.
    base::RandomNumberGenerator::XorShift128(&state.s0, &state.s1);
    cache.set(i, base::RandomNumberGenerator::ToDouble(state.s0));
  }
  pod.set(0, state);

  Smi new_index = Smi::FromInt(kCacheSize);
  native_context.set_math_random_index(new_index);
  return new_index.ptr();
}

}  // namespace internal
}  // namespace v8

总的来说，此文件与随机数生成相关的主要内容：