文章目录

1. 伪随机数
2. 模运算
3. 乘同余法随机数生成器
- (1) 原理
- (2) 程序实现
4. 混合同余法
5. mt19937

1. 伪随机数

Treap、跳跃表和随机快速排序等需要用到随机数，我们要有一种方法来生成它。不过，真正的随机性在计算机中是稀缺的，如果想要在实际应用中使用生成的随机数，太慢是不行的。一般来说，产生伪随机数 pseudorandom number 就足够了。伪随机数是指看上去像是随机的数。随机数有很多的统计性质，而伪随机数满足其中的大部分。然而，生成伪随机数也不是那么容易的。

以抛一枚硬币为例，随机生成 0 正面或者 1 负面。一种做法是借用系统时钟，将时间记作整数，一个从某个起始时刻开始计数的秒数，用其最低的二进制位。如果需要的是随机数序列，这种方法不够理想——1s1s1s 是一个长的时间段，程序运行过程中时钟可能没有变化。即使时间用微秒 ususus 计量，所生成的数的序列也远不是随机的。

我们真正需要的是随机数的序列 random number sequence 。这些数应该独立出现——如果一枚硬币抛出后是正面，那么下次抛出时出现正面还是反面应该是等可能的。

2. 模运算

在介绍随机数的生成方法之前，先了解一下模运算。我们常常用 % 取模运算符进行简单的取余数操作，但这不代表我们对模运算了解得有多么深刻。

如果 A−BA-BA−B 可以被 NNN 整除，即 AAA 与 BBB 模 NNN 同余 congruent ，记作 A≡B(modN)A\equiv B(mod\ N)A≡B(mod N) 。直观的看法是：无论 AAA 或者 BBB 被 NNN 除，所得到的余数都是相同的。例子有：
81≡61≡1(mod10)81 \equiv 61\equiv 1(mod\ 10)81≡61≡1(mod 10)

模运算有着很好的性质，它和等号一样，若 A≡B(modN)A \equiv B(mod\ N)A≡B(mod N) ，则 A+C≡B+C(modN)A+C\equiv B +C(mod\ N)A+C≡B+C(mod N) 以及 AD≡BD(modN)AD \equiv BD(mod\ N)AD≡BD(mod N) 。还有很多定理适用于模运算，有些需要数论的知识，这里不多涉及。

3. 乘同余法随机数生成器

(1) 原理

产生随机数最简单的方法是线性同余生成器，更准确的说是乘同余生成器，它于1951年被Lehmer首先提出(1951年，emmm)。数 x1,x2,…x_1,x_2,\dotsx1,x2,… 的生成满足：
xi+1=AximodMx_{i+1} = Ax_i\ \ mod\ \ Mxi+1=Axi mod M

为了开始序列的生成，我们必须给出 x0x_0x0 的某个值，也就是我们平常说的随机数种子 seed 。如果给出的 x0=0x_0=0x0=0 ，这个序列就远不是随机的。

但是如果 A,MA, MA,M 被正确选择，那么任何其他的 x0(1≤x0<M)x_0(1\le x_0 \lt M)x0(1≤x0<M) 都是同等有效的。更加有意思的是，如果 MMM 是素数，那么 xix_ixi 就绝不会是 000 。同样举一个例子，M=11,A=7,x0=1M = 11, A = 7, x_0 = 1M=11,A=7,x0=1 ，生成的序列如下：
7,5,2,3,10,4,6,9,8,1,7,5,2,3,10⋯7, 5, 2, 3, 10, 4, 6, 9, 8, 1, 7, 5, 2, 3, 10\cdots7,5,2,3,10,4,6,9,8,1,7,5,2,3,10⋯

我们还可以注意到一个特点：在 M−1=10M - 1= 10M−1=10 个数后，序列将会出现重复，或者说，这个序列的周期为 M−1M - 1M−1 。而我们，必须使得生成的随机数序列的周期尽可能大（鸽巢定理）。我们还知道的是：如果 MMM 是素数，那么总是存在某些 AAA 的值，能够得到全周期 M−1M - 1M−1，但是也有些 AAA 的值得不到这样的周期。如 M=11,A=5,x0=1M = 11, A = 5, x_0 = 1M=11,A=5,x0=1 ，其产生的序列有一个短周期 555 ：
5,3,4,9,1,5,3,4,…5, 3, 4, 9, 1, 5, 3, 4, \dots5,3,4,9,1,5,3,4,…

综上所述，如果我们选择一个很大的素数作为 MMM ，同时选择一个合适的 AAA 值，我们将得到一个周期很长的序列——Lehmer的建议是使用素数 M=232−1=2147483647M = 2^{32}- 1 = 2147483647M=232−1=2147483647 ，针对它，A=48271A = 48271A=48271 是给出全周期生成器的许多值的一个。

(2) 程序实现

这个程序实现起来很简单，类变量 state 保存 xxx 序列的当前值。不过调试随机数程序时，最好是置 x0=1x_0 = 1x0=1 ，这使得总是出现相同的随机序列。使用随机数程序时，可以用系统时钟的值，也可以让用户输入一个值作为 seed 。

另外，返回一个位于 (0,1)(0,1)(0,1) 的随机实数也是常见的；可以转换类变量为 doubledoubledouble 然后除以 MMM 得到。从而，任意闭区间 [α,β][\alpha, \beta][α,β] 内的随机数都可以通过规范化得到。

代码如下，不过有错误：

static const int A = 48271;
static const int M = 2147483647;class Random {private:int state;
public://禁止类构造函数的隐式类型转换//使用传入的初值构造stateexplicit Random(int initialValue = 1) {if (initialValue < 0)initialValue += M;//保证state为正数state = initialValue;if (state <= 0)state = 1;}//返回一个伪随机整数,改变内部state//有bugint randomInt() {return state = (A * state) % M;}//返回一个伪随机实数,(0,1)之间,同时改变内部statedouble random0_1() {return (double) randomInt() / M;}
};

问题在于：乘法可能溢出。这将影响计算的结果和伪随机性。不过经过改造，我们可以让这个过程中所有的计算在32位上进行而不溢出。计算 M/AM / AM/A 的商和余数，并分别定义为 Q,RQ, RQ,R 。此时，Q=44488,R=3399,R<QQ = 44488, R = 3399, R < QQ=44488,R=3399,R<Q ：
xi+1=AximodM=Axi−M⌊AxiM⌋=Axi−M⌊xiQ⌋+M⌊xiQ⌋−M⌊AxiM⌋=Axi−M⌊xiQ⌋+M(⌊xiQ⌋−⌊AxiM⌋)xi=Q⌊xiQ⌋+ximodQ\begin{aligned} x_{i + 1} &= Ax_i\ mod\ M \\ &= Ax_i - M \lfloor {Ax_i \over M }\rfloor \\ &= Ax_i - M \lfloor {x_i \over Q }\rfloor + M \lfloor {x_i \over Q }\rfloor - M \lfloor {Ax_i \over M }\rfloor \\ &= Ax_i - M \lfloor {x_i \over Q }\rfloor + M \Big( \lfloor {x_i \over Q }\rfloor - \lfloor {Ax_i \over M }\rfloor \Big)\\ x_i &= Q \lfloor {x_i \over Q}\rfloor + x_i\ mod\ Q \end{aligned} xi+1xi=Axi mod M=Axi−M⌊MAxi⌋=Axi−M⌊Qxi⌋+M⌊Qxi⌋−M⌊MAxi⌋=Axi−M⌊Qxi⌋+M(⌊Qxi⌋−⌊MAxi⌋)=Q⌊Qxi⌋+xi mod Q
所以有：
xi+1=A(Q⌊xiQ⌋+ximodQ)−M⌊xiQ⌋+M(⌊xiQ⌋−⌊AxiM⌋)=(AQ−M)⌊xiQ⌋+A(ximodQ)+M(⌊xiQ⌋−⌊AxiM⌋)\begin{aligned} x_{i + 1} &= A\Big( Q \lfloor {x_i \over Q} \rfloor + x_i\ mod\ Q\Big) - M \lfloor {x_i \over Q }\rfloor + M \Big( \lfloor {x_i \over Q }\rfloor - \lfloor {Ax_i \over M }\rfloor \Big)\\ &= (AQ - M) \lfloor {x_i \over Q }\rfloor + A(x_i\ mod\ Q) + M \Big( \lfloor {x_i \over Q }\rfloor - \lfloor {Ax_i \over M }\rfloor \Big)\\ \end{aligned} xi+1=A(Q⌊Qxi⌋+xi mod Q)−M⌊Qxi⌋+M(⌊Qxi⌋−⌊MAxi⌋)=(AQ−M)⌊Qxi⌋+A(xi mod Q)+M(⌊Qxi⌋−⌊MAxi⌋)
由于 M=AQ+RM = AQ + RM=AQ+R ，所以 AQ−M=−RAQ - M = -RAQ−M=−R ，得到：
xi+1=A(ximodQ)−R⌊xiQ⌋+M(⌊xiQ⌋−⌊AxiM⌋)=A(ximodQ)−R⌊xiQ⌋+Mδ(xi)δ(xi)=0or1\begin{aligned} x_{i + 1} &= A(x_i\ mod\ Q) -R\lfloor {x_i \over Q }\rfloor + M \Big( \lfloor {x_i \over Q }\rfloor - \lfloor {Ax_i \over M }\rfloor \Big)\\ &= A(x_i\ mod\ Q) -R\lfloor {x_i \over Q }\rfloor + M \delta(x_i) \\ \delta(x_i) &= 0\ or\ 1 \end{aligned} xi+1δ(xi)=A(xi mod Q)−R⌊Qxi⌋+M(⌊Qxi⌋−⌊MAxi⌋)=A(xi mod Q)−R⌊Qxi⌋+Mδ(xi)=0 or 1
验证表明，因为 R<QR < QR<Q ，所以所有余项均可计算没有溢出；另外，当 A(ximodQ)−R⌊xiQ⌋A(x_i\ mod\ Q) -R\lfloor {x_i \over Q }\rfloorA(xi mod Q)−R⌊Qxi⌋ 小于 000 时（⌊xiQ⌋,⌊AxiM⌋\lfloor {x_i \over Q }\rfloor,\lfloor {Ax_i \over M }\rfloor⌊Qxi⌋,⌊MAxi⌋ 都是整数，它们的差不是 000 就是 111 ），δ(xi)=1\delta(x_i) = 1δ(xi)=1 。所以，我们可以通过简单的测试确定 δ(xi)\delta(x_i)δ(xi) 的值。

不溢出的最终程序如下：

static const int A = 48271;
static const int M = 2147483647;
static const int Q = M / A;
static const int R = M % A;class Random {private:int state;
public://禁止类构造函数的隐式类型转换//使用传入的初值构造stateexplicit Random(int initialValue = 1) {if (initialValue < 0)initialValue += M;//保证state为正数state = initialValue;if (state <= 0)state = 1;}//返回一个伪随机整数,改变内部state//有bugint randomInt() {//return state = (A * state) % M;int tempState = A * (state % Q) - R * (state / Q);if (tempState >= 0) state = tempState;else state = tempState + M;return state;}//返回一个(0,1)之间的伪随机实数,同时改变内部statedouble random0_1() {return (double)randomInt() / M;}//返回一个闭区间[low, high]的随机整数,改变内部状态int randomInt(int low, int high) {double partitionSize = (double)M / (high - low + 1);return (int)(randomInt() / partitionSize) + low;}
};

4. 混合同余法

更复杂一点的，我们可以添加一个常数，得到下面的递推公式：
xi+1=(Axi+C)modM\displaystyle x_{i+1} = (Ax_i + C) \bmod Mxi+1=(Axi+C)modM

其中：

MMM 是模数，0<M\displaystyle 0 < M0<M ；(为了达到预期的随机效果，一般希望这个值稍稍大一点，且最好是素数)
AAA 乘数，1≤A<M\displaystyle {1 \le A \lt M}1≤A<M ；(实用起见最好取大于等于 222 的值)
CCC 增量，0<C<M\displaystyle {0\lt C\lt M}0<C<M ；(当 C=0C=0C=0 时，随机数生成过程比不等于零时要稍微快些。它的限制可能缩短这个序列的周期长度，但也有可能得到一个相当长的周期。当 C=0C=0C=0 时被称为乘同余法，C≠0C\ne 0C=0 称为混合同余法。为了一般性，建议选择采用混合同余法)
x0x_0x0 开始值，0≤x0<M\displaystyle 0 \le x_0 \lt M0≤x0<M 。

由这四个整数定义的混合同余序列得到最大周期长度 M−1M - 1M−1 的条件如下：

C,MC,MC,M 互为素数；
对于整除 MMM 的每个素数 PPP ，B=A−1B = A - 1B=A−1 是 PPP 的倍数；
如果 MMM 是 444 的倍数，那么 BBB 也是 444 的倍数。

更多随机数生成算法看这篇文章：戳这里。

5. mt19937

实际应用中，我们不需要手写随机数生成器，C++给我们提供了一个相当强大的随机数算法——mt19937 。

mt19937 是C++11中加入的新特性，作为一种随机数算法，用法与 rand() 函数类似，但是速度快、周期长——它的命名就来自周期长度 219937−12^{19937}-1219937−1 。写在程序中，这个函数的随机范围大概在 [INT_MIN,INT_MAX] 之间。它的用法非常简单：

#include<random>
#include<ctime>
std::mt19937 rnd(time(0));int main() {printf("%lld\n", rnd());return 0;
}

【算法学习】随机化算法随机数生成器和mt19937相关推荐

算法学习--排序算法--插入排序
算法学习--排序算法--插入排序插入排序算法代码实现插入排序算法插入排序(Insertion sort)是一种简单直观且稳定的排序算法.如果有一个已经有序的数据序列,要求在这个已经排好的数据序 ...
随机化算法-数值随机化算法
随机数随机数在随机化算法设计中扮演着十分重要的角色.在现实计算机上无法产生真正的随机数,因此在随机化算法中使用的随机数都是一定程度上随机的,即伪随机数. 线性同余法是产生伪随机数的最常用的方法.由线 ...
php算法学习,php算法学习之动态规划
动态规划程序设计是对解最优化问题的一种途径.一种方法,最终问题的最优解可以通过前面子问题的最优解推导出来. 对于动态规划这个算法,自己学习的还不是很透彻,简单的总结自己学习的感受是: 动态规划思想中融 ...
[算法学习]模拟退火算法（SA）、遗传算法（GA）、布谷鸟算法（CS）、人工蜂群算法（ABC）学习笔记---附MATLAB注释代码
目录 1.模拟退火算法(Simulated Annealing,SA) 1.1 本质: 1.2 算法思想 1.3 SA流程图 1.4 模拟退火过程 1.5 SA解决TSP问题 1.6 SA改进方向 1 ...
算法学习之算法的引入
一．算法的起始 1.第一次尝试如果 a+b+c=1000,且 a^2+b^2=c^2(a,b,c 为自然数),如何求出所有a.b.c可能的组合? (可以考虑到百钱白鸡) 枚举法 # 注意是三重循环 ...
【C语言】算法学习·KMP算法
KMP算法(全称Knuth-Morris-Pratt字符串查找算法,由三位发明者的姓氏命名)是可以在文本串s中快速查找模式串p的一种算法. 要想知道KMP算法是如何减少字符串查找的时间复杂度的,我们不 ...
【C语言】算法学习·回溯算法
目录一.全排列问题二.N 皇后问题三.最后总结回溯算法基本框架解决一个回溯问题,实际上就是一个决策树的遍历过程.你只需要思考 3 个问题: 1.路径:也就是已经做出的选择. 2.选择列表:也 ...
C++算法学习(贪心算法)
贪心算法 1.目标 2.方法 3.例题 [122. 买卖股票的最佳时机 II](https://leetcode-cn.com/problems/best-time-to-buy-and-sell-s ...
C++算法学习(动态规划算法)
动态规划算法 1.目标 2.方法 3.过程 4.例题 (1)[力扣:5. 最长回文子串](https://leetcode-cn.com/problems/longest-palindromic-su ...

【算法学习】随机化算法随机数生成器和mt19937

文章目录

1. 伪随机数

2. 模运算

3. 乘同余法随机数生成器

(1) 原理

(2) 程序实现

4. 混合同余法

5. mt19937

【算法学习】随机化算法随机数生成器和mt19937相关推荐

最新文章

热门文章

【算法学习】随机化算法 随机数生成器和mt19937

文章目录

1. 伪随机数

2. 模运算

3. 乘同余法随机数生成器

(1) 原理

(2) 程序实现

4. 混合同余法

5. mt19937

【算法学习】随机化算法 随机数生成器和mt19937相关推荐

最新文章

热门文章

【算法学习】随机化算法随机数生成器和mt19937

【算法学习】随机化算法随机数生成器和mt19937相关推荐