C++随机数之线性同余发生器
Posted -飞鹤-
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了C++随机数之线性同余发生器相关的知识，希望对你有一定的参考价值。
1. 随机数
随机数，也即“随机选择的数”，是在一个有限数集上的一个一致分布的随机序列。随机数在许多方面有应用，如仿真、抽样、数值分析、计算机程序、娱乐等方面都有所应用。计算机用确定算法计算出来的随机数系列是伪随机数，并不是真正的随机数，但是其具有随机数的统计特征，如均匀性、独立性等。
在上世纪初，冯诺伊曼建议用平方，然后取结果中间的数据作为随机数(平方取中法)。结果看起来修似乎是随机的，但是很多人对此质疑。实际上这并不是好的随机方法，尤其是针对一些短循环序列。
2. 线性同余序列
2.1. 概念
后来，有人发明了线性同余方程，用此方程生成的一组序列，随机效果不错。线性同余方程生成的序列，被称为线性同余序列。线性同余序列在一个周期内重复出现。并且方程中有四个魔数，但是这里有四个魔数，对随机的效果影响巨大。线性同余发生器(Linear congruential generator,LCG)用于生成随机序列。
LCG定义：
X
n
+
1
=
(
a
X
n
+
c
)
m
o
d
m
X_{n+1}=(aX_{n}+c)\\mod\\;m
Xn+1=(aXn+c)modm
- X是随机数序列
- m,0<mm,0<m，模
- a,0<a<ma,0<a<m，乘子
- c,0≤c<mc,0≤c<m，增量，也叫做偏移量
- X0 ,0≤X0<mX0,0≤X0<m，开始值，通常叫做“种子”seed
当c=0时候，LCG就变换成了乘法同余发生器，multiplicative congruential generator (MCG)， c≠0时叫做混合同余发生器，mixed congruential generator，MCG。
2.2. 示例
上面的线性同余方程，选定a、X、c和m四个魔数，可以生成线性同余序列，线性同余序列在一个周期内存是随机的，独立的，一个周期完成之后，循环进行。
当m=9,a=2,c=0,seed=1时，得：
X
0
=
s
e
e
d
=
1
X
2
=
(
a
∗
X
1
+
c
)
m
o
d
m
=
(
2
×
1
+
0
)
%
9
=
2
X
3
=
(
a
∗
X
2
+
c
)
m
o
d
m
=
(
2
×
2
+
0
)
%
9
=
4
X
4
=
(
a
∗
X
3
+
c
)
m
o
d
m
=
(
2
×
4
+
0
)
%
9
=
8
X
5
=
(
a
∗
X
4
+
c
)
m
o
d
m
=
(
2
×
8
+
0
)
%
9
=
7
X
6
=
(
a
∗
X
5
+
c
)
m
o
d
m
=
(
2
×
7
+
0
)
%
9
=
5
X
7
=
(
a
∗
X
6
+
c
)
m
o
d
m
=
(
2
×
5
+
0
)
%
9
=
1
\\begin{aligned} &X_0 = seed = 1\\\\ &X_2 = (a∗X_1+c)\\mod\\;m=(2×1+0)\\%9=2\\\\ &X_3 = (a∗X_2+c)\\mod\\;m=(2×2+0)\\%9=4\\\\ &X_4 = (a∗X_3+c)\\mod\\;m=(2×4+0)\\%9=8\\\\ &X_5 = (a∗X_4+c)\\mod\\;m=(2×8+0)\\%9=7\\\\ &X_6 = (a∗X_5+c)\\mod\\;m=(2×7+0)\\%9=5\\\\ &X_7 = (a∗X_6+c)\\mod\\;m=(2×5+0)\\%9=1\\\\ \\end{aligned}
X0=seed=1X2=(a∗X1+c)modm=(2×1+0)%9=2X3=(a∗X2+c)modm=(2×2+0)%9=4X4=(a∗X3+c)modm=(2×4+0)%9=8X5=(a∗X4+c)modm=(2×8+0)%9=7X6=(a∗X5+c)modm=(2×7+0)%9=5X7=(a∗X6+c)modm=(2×5+0)%9=1
数据示图如下：
如上图所示，m、a、c取值不同，其循环的周期长度是不同的，并且随机的效果也是不同的。
2.3. 最大周期长度
m、a、c取值不同，其循环的周期长度是不同的。为了使得线性同余序列的周期长度达到最大m，需要满足以下3个条件：
- m and c are relatively prime（互为质数，即两个数没有除1之外的公约数）
- a-1 is divisible by all prime factors of m（m的质因数（可以整除m的所有质数），都能整除a-1，也即a-1是m的所有质因数的倍数）
- a-1 is divisible by 4 if m is divisible by 4 (如果m被4整除，a-1也被4整除，即m是4的倍数，a-1也要是4的倍数）
详细证明过程，感觉兴趣的可以看《计算机程序设计艺术》第二卷的随机数章节，有详细证明过程的介绍。
-
是质数
/*判断是否为质数*/ bool prime(int N) { if (N < 2) return false; if (N == 2) return true; for (int i = 2; i < std::sqrt(N*1.0); i++) { if (N%i == 0) //i能被j整除，N不是质数 { return false; } } return true; } -
互质
bool isrp(int a, int b) { if(a <=0 || b<=0 || a==b) { // 互质整数不能小于或等于0 return false; } else if(a==1 || b==1) { // 两个正整数中，只有其中一个数值为1，两个正整数为互质数 return true; } else { // 求出两个正整数的最大公约数 while(1) { int t = a%b; if(t == 0) { break; } else { a = b; b = t; } } if( b > 1){ // 如果最大公约数大于1，表示两个正整数不互质 return false; }else{ // 如果最大公约数等于1,表示两个正整数互质 return true; } } } -
所有质因数
class QualityFactor { public: vector<int> m_vInt; public: void QFContract(long a) //用短除法对合数进行分解 { while(a>1) { for(int i=2;i<=a;i++) { if(a%i==0) //短除法 { a = a/i; if (std::find(m_vInt.begin(), m_vInt.end(), i) == m_vInt.end()) m_vInt.push_back(i); break; } } } } }; -
所有质因数整除
int MaxCommondMultiple(const vector<int>& _vInt, BOOL _bFour, int _nSrc) { int nCommonMultiple = 1; while (nCommonMultiple < _nSrc) { for (auto it = _vInt.begin(); it != _vInt.end(); ++it) { nCommonMultiple*= *it; if (_bFour && (nCommonMultiple%4 != 0)) { continue; } BOOL bFind = TRUE; for (auto it = _vInt.begin(); it != _vInt.end(); ++it) { if (nCommonMultiple % *it != 0) { bFind = FALSE; break; } } if (bFind) { return nCommonMultiple; } } } return 0; }- 随机数函数
static int seed = 0; int random(int a, int b, int m) { seed = (seed * a + b) % m; return seed; }
输入最大周期长度m，即可以根据上述函数生成相应的a和c。
2.4. 其他
2.4.1. 标准库rand
m、a、c、X是一级相关数，选定m之后，同样会有很多组对应的a、c、X。不同的数据，对应的随机数分布是不一样的。随机效果也是大大不同。理论上m选择的范围越方，其效果越好。而VS实现的std::rand()即是采用线性同余方程实现的。
*/
int __cdecl rand (
void
)
{
_ptiddata ptd = _getptd();
return( ((ptd->_holdrand = ptd->_holdrand * 214013L
+ 2531011L) >> 16) & 0x7fff );
}
其周期较小，所以随机效果不是特别好。另外其中的几个常数，是经常大量实验选择的，是随机效果相对较好的一组数。
2.4.2. 自定义rand
什么时候需要自定义rand()呢？例如我要取0-10之间的随机数，并且每次不能相同，且可以取完一个周期。一般情况下我们会使用“洗牌”算法，即在数据序列中，依次取出一个数与后面的随机取的数交换，即可以达到效果。标准库提供了std::shuffle实现此功能。
std::vector<int> v = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
std::shuffle (v.begin (), v.end (), std::default_random_engine (0));
大多数时候，我们使用上面的方法，可以得到效果非常好的随机序列。但是如果我需要一个0-0xFFFFFFFF之间的随机序列呢？此时就不方便使用构造一个巨大的数据序列之后，再调用洗牌函数打乱数据。因为数据序列的内存消耗巨大。
此时，就需要使用一个自定义周期长度为0xFFFFFFF的rand()函数了。即，我们选定m=0xFFFFFF，根据三个条件调用相应的函数求出相应的m、c、X。但是因为参数的可能性太多，随机性差异巨大，需要经过大量实验，才能选择一组较好的数据组合供使用。
示例代码
2.4.3. 最佳随机数
线性同余发生器一直是随机数生成的主要方式，但是其随机效果不是特别好。C++11引入了梅森旋转演算法的随机数引擎，是目前效果最好的随机算法。std::default_random_engine默认即使用梅森旋转演算法。也可以直接使用梅森旋转演算法引擎std::mt19937。
以上是关于C++随机数之线性同余发生器的主要内容，如果未能解决你的问题，请参考以下文章