1、原理:
C语言中伪随机数生成算法实际上是采用了"线性同余法"。具体的计算如下:
seed = (seed * A + C ) % M
其中A,C,M都是常数(一般会取质数)。当C=0时,叫做乘同余法。
假设我们定义随机数函数
void rand(int &seed)
{
seed = (seed * A + C ) % M;
}
每次调用rand函数都会产生一个随机值赋值给seed,可以看出实际上用rand函数生成的是一个递推的序列,一切值都来源于最初的 seed。所以当初始的seed取一样的时候,得到的序列都相同。
我们称seed为种子,一个伪随机数常用的原则就是M尽可能的大。例如,对于32位的机器来说,选择M=2^31-1=2147483647, A=7^5=16807时可以取得最佳效果。
2、代码实现:
现在我们来看看levelDB里随机数Random类是如何实现的:
在Random类中,A为16807,M为2147483647,C为0;
// Copyright (c) 2011 The LevelDB Authors. All rights reserved.
// Use of this source code is governed by a BSD-style license that can be
// found in the LICENSE file. See the AUTHORS file for names of contributors.
#ifndef STORAGE_LEVELDB_UTIL_RANDOM_H_
#define STORAGE_LEVELDB_UTIL_RANDOM_H_
#include <stdint.h>
namespace leveldb {
// A very simple random number generator. Not especially good at
// generating truly random bits, but good enough for our needs in this
// package.
class Random
{
private:
uint32_t seed_;
public:
// 0x7fffffffu == 2147483647L == 2^31-1 == 01111111 11111111 11111111 11111111
// 表达式s & 0x7fffffffu,确保结果值在[0,2147483647]范围内
explicit Random(uint32_t s) : seed_(s & 0x7fffffffu)
{
// Avoid bad seeds.
if (seed_ == 0 || seed_ == 2147483647L)
{
seed_ = 1;
}
}
// 16807随机数
uint32_t Next()
{
//01111111 11111111 11111111 11111111
static const uint32_t M = 2147483647L; // 2^31-1
//0100 0001 1010 0111
static const uint64_t A = 16807; // bits 14, 8, 7, 5, 2, 1, 0
// We are computing
// seed_ = (seed_ * A) % M, where M = 2^31-1
//
// seed_ must not be zero or M, or else all subsequent computed values
// will be zero or M respectively. For all other values, seed_ will end
// up cycling through every number in [1,M-1]
uint64_t product = seed_ * A;
// Compute (product % M) using the fact that ((x << 31) % M) == x.
seed_ = static_cast<uint32_t>((product >> 31) + (product & M)); // 为什么会这样计算?请看下面我做的证明和分析
// The first reduction may overflow by 1 bit, so we may need to
// repeat. mod == M is not possible; using > allows the faster
// sign-bit-based test.
if (seed_ > M)
{
seed_ -= M;
}
return seed_;
}
// Returns a uniformly distributed value in the range [0..n-1]
// REQUIRES: n > 0
uint32_t Uniform(int n) { return Next() % n; }
// Randomly returns true ~"1/n" of the time, and false otherwise.
// REQUIRES: n > 0
bool OneIn(int n) { return (Next() % n) == 0; }
// Skewed: pick "base" uniformly from range [0,max_log] and then
// return "base" random bits. The effect is to pick a number in the
// range [0,2^max_log-1] with exponential bias towards smaller numbers.
uint32_t Skewed(int max_log)
{
return Uniform(1 << Uniform(max_log + 1));
}
};
} // namespace leveldb
#endif // STORAGE_LEVELDB_UTIL_RANDOM_H_
3、源码注释提出了个问题,为什么会这样计算?通过公式seed = (seed * A + C ) % M可以知道,
需要判断(product % M)是否与static_cast((product >> 31) + (product & M))相等?
4、证明等式(product%M) == (product>>31)+(product&M)成立。注:M等于2^31-1。
证明:
因为product类型是uint64_t,可以将product的二进制从左到右分解成高33位和低31位,假设高33位的值为H,低31位的值为L,
则product相当于高33位向左移动了31位加上低31位,即H<<31+L。
则product等于H*2^31+L。由源码知道product=seed_*A,而seed_和A都小于M,则H肯定小于M。
从而我们可以得到:
表达式左边product%M = (H*2^31+L)%M = (H*M+H+L)%M = H+L。
表达式右边(product>>31) + (product&M) = (H*2^31 +L)>>31+L = (H*2^31+L)/2^31+L = H+L。
注:细心的同学可能发现,低31位的值L可能等于M,那么表达式左边就等于H了,此时表达式右边等于H+M。
回归源码,我们可以看到这个判断条件
seed_ = static_cast<uint32_t>((product >> 31) + (product & M));
if (seed_ > M)
{
seed_ -= M;
}
经过这个判断条件,表达式右边也就等于H了。
综上左边表达式和右边表达式的结果相等,所以等式成立。
5、经过证明4,我们了解到,可以将一个数,转换成二进制,将其从左到右分解成高m位,假设值为H,低n位,假设值为L,
则这个数可表示为H*2^n+L;
上面只是将一个数的二进制从左到右分解成高低两段。
同理我们将一个数的二进制从左到右分解成3段,第一段m位,假设值为H;第二段b位,假设值为C;第三段n位,假设值为L,
则这个数可表示为H*2^(b+n)+C*2^n+L。
同理也可以了解到一个数分解成n段。
这里为什么要说这么一大段,因为这样的分解是很有作用的。
比如,一个非常长的二进制数,现有的数据类型存放不下,怎么办?就是将这个数分解成多段来解决问题。
6、通过上面第4点和第5点的分析,我们来反观A这个数,也就是16807,它的二进制是01 000001 1 01 001 1 1,
如果我们将位置为1的地方分段,则将其分成7段,分别是2位、6位、1位、2位、3位、1位、1位,值都是1,
则A可表示为1*2^14+1*2^8+1*2^7+1*2^5+1*2^2+1*2^1+1*2^0。如果我们将值1改为seed_,也就是得到seed_*A。
7、证明源码注释里面的等式(x<<31)%M == x成立。注:M等于2^31-1,x为小于M大于0的正整数。
证明:
计算表达式左边(x << 31) % M,由于x<<31等于x*2^31,
则(x << 31) % M=(x*2^31)%M=(x + x*(2^31-1))%M=(x + x*M)%M=x%M=x