C++随机数之线性同余发生器

1. 随机数

随机数，也即“随机选择的数”，是在一个有限数集上的一个一致分布的随机序列。随机数在许多方面有应用，如仿真、抽样、数值分析、计算机程序、娱乐等方面都有所应用。计算机用确定算法计算出来的随机数系列是伪随机数，并不是真正的随机数，但是其具有随机数的统计特征，如均匀性、独立性等。
在上世纪初，冯诺伊曼建议用平方，然后取结果中间的数据作为随机数(平方取中法)。结果看起来修似乎是随机的，但是很多人对此质疑。实际上这并不是好的随机方法，尤其是针对一些短循环序列。

2. 线性同余序列

2.1. 概念

后来，有人发明了线性同余方程，用此方程生成的一组序列，随机效果不错。线性同余方程生成的序列，被称为线性同余序列。线性同余序列在一个周期内重复出现。并且方程中有四个魔数(常数)，但是这四个魔数，对随机的效果影响巨大。线性同余发生器(Linear congruential generator,LCG)用于生成随机序列。
LCG定义：
$$X_{n+1}=(a{X}_n+c)\mathrm{m}\mathrm{o}\mathrm{d}m$$

• X是随机数序列
• m，0<m，模
• a，0<a<m，乘子
• c，0≤c<m，增量，也叫做偏移量
• X₀ , 0≤X₀<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个条件：

1. m and c are relatively prime（互为质数，即两个数没有除1之外的公约数）
2. a-1 is divisible by all prime factors of m（m的质因数（可以整除m的所有质数），都能整除a-1，也即a-1是m的所有质因数的倍数）
3. 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。