如何从rand7生成rand5
首先是和某个知名算法题相反的问题。
这里给定一个可以概率均等地生成0到6的随机数生成器,要求用这个生成器创造出能概率均等地生成0到4的随机数生成器。
有人可能会立刻给出这样的答案:
func rand5() int { |
return rand7() % 5 |
} |
然而这个答案只满足了输出的范围在0到4,不满足概率均等,所以不正确。这种时候列表法的作用就显现出来了:
| rand7的输出 | rand7 % 5 |
|---|---|
| 0 | 0 |
| 1 | 1 |
| 2 | 2 |
| 3 | 3 |
| 4 | 4 |
| 5 | 0 |
| 6 | 1 |
发现问题了吗,0和1出现了两次,它们出现的概率是其他数字的两倍,因此概率分布并不均等。
通过列表法我们其实也发现了这个问题真正的解法:除掉5和6的话剩下的输出不仅符合要求概率也是均等的,所以代码会变成这样:
func rand5() int { |
n := rand7() |
for n >= 5 { |
n = rand7() |
} |
return n |
} |
上面的代码其实就是随机采样算法中非常重要的一种:拒绝采样。同样上面的rand7生成rand5也可以归类成一大类问题:给定一组满足规律或者特征是g(x)的样本,现在需要从这些样本中筛选出或者生成另一组满足特征是f(x)的样本。解决这类问题的算法很多,而拒绝采样是比较直观的:判断g(x)的样本是否符合要求,不符合的就排除取下一个样本,符合要求的就归类到满足f(x)的样本集合中。
按照这个角度来看,上面的问题可以划分为几个基本元素:
- g(x):rand7
- f(x): rand5
- 需要拒绝的样本:大于等于5的整数
拒绝采样在大多数时间都能获得理想的结果,但还有采样率需要注意。采样率就是g(x)的样本中有多少可以被接受,采样率太低的时候意味着算法的效率也会非常低。所以我们简单算算rand5的采样率,是七分之五,大约70%。这个概率不大不小,勉强合适。
当采样率过低的时候要么得改变拒绝采样的标准或范围,要么就不能再使用拒绝采样了。
go标准库的做法
标准库里当然不会有rand5和rand7,但它提供了一个叫Int63n的函数,它解决的问题是如何从均匀分布在[0, 2⁶⁴)范围上的随机整数中生成均匀分布在范围[0, n)上的随机整数。换句话说,虽然范围不一样了,但还是同一个问题。
我们肯定不能像上面那样把大于等于n的样本全丢了,因为2⁶⁴包含至少1844京(1E16)个整数样本,这会带来低得无法接受的采样率。
但因为我们是用mod来选择范围内的随机数的,因此我们可以选择n的倍数,这个证明太简单了,列表法加归纳就行。或者还可以这么想,有一个整数常数C,x % n和Cx % n能产生的输出的种类和它们的数量都是完全相同的,所以如果样本能均匀分布在[0, n)的范围上,那么范围[0, C·n]只不过是[0, n)重复了C次,所以样本在每一段上都是均匀分布的,整个合起来的区间上也是均匀的。
有了常数C,这样使得我们可以尽可能地让更多的样本被采样,这样能降低重试的次数。
C其实也很好选择,取一个2⁶⁴内的n的最大的倍数就行,如果C本身能被2⁶⁴整除,那么C就是2⁶⁴/n。
所以标准库是这样写的:
func (r *Rand) Int63n(n int64) int64 { |
if n <= 0 { |
panic("invalid argument to Int63n") |
} |
if n&(n-1) == 0 { // n is power of two, can mask |
return r.Int63() & (n - 1) |
} |
max := int64((1 << 63) - 1 - (1<<63)%uint64(n)) |
v := r.Int63() |
for v > max { |
v = r.Int63() |
} |
return v % n |
} |
代码还是很简单的,超过C·n的样本全部拒绝采样,剩下的样本就能保证在mod n的时候获得分布均匀的随机整数了。
采样率是多少?我们可以利用拒绝率来反推,这里拒绝率还挺好算的,就是(1<<63)%uint64(n),算下来拒绝了少的时候是百亿分之一,多的时候是数千万分之一——都很低,基本上大多数时间最多两次重试就能获得想要的结果了。
但作为标准库,它的性能还不够,尤其是go的编译优化非常弱的现实下,更需要高效的算法弥补。问题出在哪?首先不是采样率,这个采样率是足够的,问题出在它需要两次64位除法,除法运算相比其他运算比如右移要慢很多,何况还是两次,别的语言中的标准库采用的算法只需要0到1次除法就够了。
好在go提供了math/rand/v2,采用了更高效的算法。
新算法依旧基于拒绝采样,但对采样的范围进行了变更,具体是这样的:
- 依然用概率均等的rand64生成一个随机整数x
- 现在把
x*n,这样生成的值的范围是[0, n·2⁶⁴) - 因为是对已有范围的等比扩大,所以
x*n在[0, n·2⁶⁴)依旧是均匀分布的(不过要注意,范围扩展了,但样本的总量不变还是2⁶⁴个) [0, n·2⁶⁴)可以均等分成n个范围,分别是[0, 2⁶⁴),[2⁶⁴, 2*2⁶⁴),[2*2⁶⁴, 3*2⁶⁴)...[(n-1)*2⁶⁴, n*2⁶⁴)- 这样每一个均等分割的范围整除以2⁶⁴就可以得到一个整数k,k一定在
[0, n)内 - k可以当作符合要求的结果,而整除以2⁶⁴实际上可以转换成位运算,这样除法运算可以减少一次。
新算法有几个问题,第一个是x*n在大多数情况下会超过2⁶⁴,但这不用担心,因为go提供了高性能128位整数运算。
第二个是x*n虽然在[0, n·2⁶⁴)均匀分布,但我们怎么保证在均等分割的每个[(k-1)*2⁶⁴, k*2⁶⁴)上也是均等分布的呢?
答案是如果只有上面写的六个步骤,我们保证不了。原因是因为要想保证x*n均匀分布在每个[(k-1)*2⁶⁴, k*2⁶⁴)上,我们就要保证x本身要均匀分布在[(k-1)*(2⁶⁴/n), k*(2⁶⁴/n))上,换人话说,就是把2⁶⁴分割成n份,每份里的样本数量都要一致。因为我们的样本都是整数而不是实数,所以动动脚趾就能想到很多数是不能整除2⁶⁴的,因此会留下“余数”。但我们的新算法实际上假设了x均匀分布在2⁶⁴分割出来的均等的范围内。不能整除的情况下意味着即使按最均匀的办法分割,也会存在一部分范围比其他的范围多几个样本或者少几个样本,会多还是会少取决与你对2⁶⁴/N取整的方式。
但这问题不大,通常分段直接的数量差异对概率产生的误差非常小,比如我们把n取6,按尽可能均匀的分割,就存在4个分段比剩下的分段里的样本总数多1个,但每个分段的样本数量都有超过3E18个,多一个还是多两个带来的影响几乎可以忽略不计。
然而标准库最重要的是要保证结果的正确性,即使可能性是3E18分之一,它依旧不是0,函数的实现是不正确的,更何况根据n的选择,n越大分段的样本数量越少,分段之间数量差异带来的影响就会越来越大,总有一个n能让结果的误差大到无法忽略。
问题其实也好解决,因为我们知道始终会有一些分组的样本是多余的,我们只要保证分组里的样本数量一致就行,不需要关心具体剔除的样本是什么。假设我们采用向下取整的办法,那么会存在一些理论上应该在分段k上的样本跑到k+1的分组上,这些样本通常分布在分段的起始位置上,我们可以把这些样本拒绝采样,这样比较容易实现。这些样本乘以n之后会落在[k*2⁶⁴, k*2⁶⁴+(2⁶⁴%n))上。
剔除这些样本后,我们就能保证x*n在每个[(k-1)*(2⁶⁴/n), k*(2⁶⁴/n))上都是均匀分布的了。
思路理解了看代码也就不难了:
func (r *Rand) uint64n(n uint64) uint64 { |
if is32bit && uint64(uint32(n)) == n { |
return uint64(r.uint32n(uint32(n))) |
} |
if n&(n-1) == 0 { // n is power of two, can mask |
return r.Uint64() & (n - 1) |
} |
hi, lo := bits.Mul64(r.Uint64(), n) |
if lo < n { |
thresh := -n % n // 2⁶⁴ % n 的简化形式 |
for lo < thresh { |
hi, lo = bits.Mul64(r.Uint64(), n) |
} |
} |
return hi |
} |
精髓在于利用(x*n) >> 64来避免了x % n带来的除法运算。而且新算法不用一开始就算余数,因此运气好的时候可以一次除法都不做。
还有一个小疑问,128位乘法够了吗?肯定够了,因为n最大也只能取到2⁶⁴,这意味这x*n的范围最大也只到[0, 2⁶⁴·2⁶⁴),128位乘法刚好够用。
最后做下性能测试,标准库里已经提供了,在我的10代i5上旧算法一次调用需要18ns,新算法只需要5ns,两者使用的随机数发生器是一样的,因此可以说新算法快了3倍,提升还是很可观的。
从rand5生成rand7
上一节讨论了从更大的样本空间里筛选出特定特征的子集,这一节我们反过来:从范围更小的样本空间里派生出有某一特征的超集。
同时,这也是一道常见的中等难度的算法题。
首先要考虑的是如何把受限的样本空间尽量扩张。上一节我们用乘法来扩展了样本分布的范围,然而乘法尤其是乘以常数是没法增加样本数量的,因此这个做法只能pass。加法可以平移样本的范围,但也不能增加样本总量,而且我们需要样本空间是[0, x)平移之后起点都变了,因此也不行。
那剩下的可行的也最稳定的办法是rand5() * rand5()。它像乘法一样能扩张样本的范围,同时因为不是乘以常数因此还有机会产生新的样本。我们列个表看看:
| rand5 | rand5 | rand5 * rand5 |
|---|---|---|
| 0 | 0 | 0 |
| 0 | 1 | 0 |
| 0 | 2 | 0 |
| 0 | 3 | 0 |
| 0 | 4 | 0 |
| 1 | 0 | 0 |
| 1 | 1 | 1 |
| 1 | 2 | 2 |
| 1 | 3 | 3 |
| 1 | 4 | 4 |
| 2 | 0 | 0 |
| 2 | 1 | 2 |
| 2 | 2 | 4 |
| 2 | 3 | 6 |
| 2 | 4 | 8 |
| 3 | 0 | 0 |
| 3 | 1 | 3 |
| 3 | 2 | 6 |
| 3 | 3 | 9 |
| 3 | 4 | 12 |
| 4 | 0 | 0 |
| 4 | 1 | 4 |
| 4 | 2 | 8 |
| 4 | 3 | 12 |
| 4 | 4 | 16 |
确实有新样本出现了,但不够连续,比如没有7和10。因此这条路是不通的。
这时候就要上原汁原味的拒绝采样算法了,我们使用5 * rand5 + rand5:
| rand5 | rand5 | 5 * rand5 + rand5 |
|---|---|---|
| 0 | 0 | 0 |
| 0 | 1 | 1 |
| 0 | 2 | 2 |
| 0 | 3 | 3 |
| 0 | 4 | 4 |
| 1 | 0 | 5 |
| 1 | 1 | 6 |
| 1 | 2 | 7 |
| 1 | 3 | 8 |
| 1 | 4 | 9 |
| 2 | 0 | 10 |
| 2 | 1 | 11 |
| 2 | 2 | 12 |
| 2 | 3 | 13 |
| 2 | 4 | 14 |
| 3 | 0 | 15 |
| 3 | 1 | 16 |
| 3 | 2 | 17 |
| 3 | 3 | 18 |
| 3 | 4 | 19 |
| 4 | 0 | 20 |
| 4 | 1 | 21 |
| 4 | 2 | 22 |
| 4 | 3 | 23 |
| 4 | 4 | 24 |