隨機數漫談

apocelipes發表於2024-07-03

隨機數對程式設計來說很重要,今天就從幾方面探討下一些常見的隨機數相關的問題。

本文只討論整數相關的隨機數,另外需要你對機率論有最基本的瞭解(至少知道古典概型是什麼)。

本文索引

  • 如何從rand7生成rand5
    • go標準庫的做法
  • 從rand5生成rand7
  • 充分利用每一個bit
  • 帶有權重的隨機數
  • 隨機數種子
    • golang和python3選擇了那些源作為種子
  • 總結

如何從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 % nCx % 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,採用了更高效的演算法。

新演算法依舊基於拒絕取樣,但對取樣的範圍進行了變更,具體是這樣的:

  1. 依然用機率均等的rand64生成一個隨機整數x
  2. 現在把x*n,這樣生成的值的範圍是[0, n·2⁶⁴)
  3. 因為是對已有範圍的等比擴大,所以x*n[0, n·2⁶⁴)依舊是均勻分佈的(不過要注意,範圍擴充套件了,但樣本的總量不變還是2⁶⁴個)
  4. [0, n·2⁶⁴)可以均等分成n個範圍,分別是[0, 2⁶⁴), [2⁶⁴, 2*2⁶⁴), [2*2⁶⁴, 3*2⁶⁴) ... [(n-1)*2⁶⁴, n*2⁶⁴)
  5. 這樣每一個均等分割的範圍整除以2⁶⁴就可以得到一個整數k,k一定在[0, n)
  6. 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

沒錯,正好產生了均等分佈的0到24的整數。很神奇吧,其實想明白為什麼不難。我們先看5 * rand5,這樣或產生0、5、10、15、20這五個數字,我們要想有機會生成連續的整數,就一定需要把缺少的1到4,11到14這些數字補上。這時候正巧一個+ rand5就可以把這些缺的空洞全部填上。當然用進位來理解會更簡單。

總結: n * randn + randn可以產生連續的範圍在[0, n*n)的均勻分佈的整數。注意這裡沒有結合律,因為randn每次的結果都是不一樣的

這個樣本空間是遠超rand7的要求的,因此現在問題回到了第一節:如何從rand25生成rand7?現在大家都知道了:

func rand7() int {
    x := 5*rand5() + rand5()
    max := 25 - 25%7
    for x >= max {
        x = 5*rand5() + rand5()
    }
    return x % 7
}

你也可以改寫成上一節說的更高效的做法。

充分利用每一個bit

rand.Uint64()返回的隨機數有足足64bits,而我們通常不需要這麼大的隨機數,舉個例子,假如我們只需要0到15的隨機數,這個數字只需要4bits就夠了,如果用rand.Uint64N(16)來生成,我們會浪費掉60bits的資料。

為啥這麼說?因為rand.Uint64()保證能機率均等的生成[0, 2⁶⁴)範圍內的整數,這其實說明了兩件事,第一是這個隨機數的每一個bit也都是隨機的,這個很明顯;第二個是每個bits是0還是1的機率也是均等的,這個可以靠列表加歸納法得出。我們來看看第二點這麼證明。

先假設我們有一個能均勻生成[0, 8)範圍內數字的隨機數發生器,現在我們看看所有可能的情況:

生成的數字 二進位制表示 從左到右第一位 從左到右第二位 從左到右第三位
0 000 0 0 0
1 001 0 0 1
2 010 0 1 0
3 011 0 1 1
4 100 1 0 0
5 101 1 0 1
6 110 1 1 0
7 111 1 1 1

不難注意到,三個bit各有八種輸出,0佔四種,1佔剩下四種,0和1的機率均等。這個結論能推廣到任意的[0, 2^n)範圍上。

同樣,基於這個結論,我們還能得到這樣一個結論,任意連續的n個bit,都能產生均勻分佈在[0, 2^n)上的隨機數,這個證明太簡單了,所以我們集中注意力就行了。

現在回頭看看為什麼我說會浪費60bits,因為根據上面的結論,我們64位的隨機整數完全可以按每四位劃分一次,這樣可以分成16組,而每組正好能產生[0, 16)範圍內的隨機數,且機率同樣的均等的。也就是說一次rand.Uint64()理論上應該可以產生16個我們需要的隨機數,但實際上我們只生成了一個。這裡就有很大的提升空間了。

怎麼做呢?我們需要一點位運算,把64位整數分組四位一組:

n := rand.Uint64()
mask := uint64(0xf) // 0b1111
for i := range 10 {
	a[i] = int(n & mask) // 只讓最右邊四位有效(書寫順序的右邊,這裡不討論大小端了因為說明起來太麻煩)
	n >>= 4 // 把剛剛使用過的四位去掉
}

程式碼很簡單,下面看看效能測試:

// 不要這樣寫程式碼
// 我這麼做是為了避免記憶體分配和回收會對測試產生不必要的雜音
func getRand(a []int) {
	if len(a) < 10 {
		panic("length wrong")
	}
	for i := range 10 {
		a[i] = int(rand.Int32N(16))
	}
}

// 不要這樣寫程式碼
// 我這麼做是為了避免記憶體分配和回收會對測試產生不必要的雜音
func getRandSplit(a []int) {
	if len(a) < 10 {
		panic("length wrong")
	}
	n := rand.Uint64()
	mask := uint64(0xf)
	for i := range 10 {
		// 這裡不需要mod
		a[i] = int(n & mask)
		n >>= 4
	}
}

func BenchmarkGetRand(b *testing.B) {
	var a [10]int
	for range b.N {
		getRand(a[:])
	}
}

func BenchmarkGetRandSplit(b *testing.B) {
	var a [10]int
	for range b.N {
		getRandSplit(a[:])
	}
}

測試結果:

goos: windows
goarch: amd64
pkg: benchrand
cpu: Intel(R) Core(TM) i5-10200H CPU @ 2.40GHz
BenchmarkGetRand-8        	15623799	        79.31 ns/op	       0 B/op	       0 allocs/op
BenchmarkGetRandSplit-8   	100000000	        11.18 ns/op	       0 B/op	       0 allocs/op

充分利用每一個bit之後我們的效能提升了整整6倍

到目前為止還不錯,如果你不在乎生成的隨機數的機率分佈或者你只想生成[0, 2^n)範圍的隨機數且這個n可以整除64,那麼可以直接跳到下一節繼續看了。

接著往下看的人肯定是希望不管在什麼範圍內都能生成機率均勻的隨機數且儘量多利用已生成的隨機bits的。但事情往往不盡人意,比如,[0, 13)[0, 7)就是兩個例子。前者右邊界不是2的冪,後者雖然是2的冪但3不能整除64。

我們先從簡單的問題開始解決,[0, 13)。受先表示數字12至少得用4個bit,4能整除64,所以我們還可以每4個連續的bit分割成一組,但這時機率分佈均勻的條件是滿足不了的。無法保證的原因很簡單,和我們第一節裡說的“rand7生成rand5”的情況一樣,每個均勻分割出來的一組連續的bits的組合裡有我們不需要的樣本存在。處理這個情況的方法在第一節裡已經有表述了,那就是拒絕取樣。確定拒絕的範圍也使用第一節說到的辦法,注意到每一組bits能生成[0, 16)的隨機數,再考慮到13本身是素數,這裡只需要簡單地把≥13的樣本全部剔除即可。

所以程式碼變成了下面這樣:

func getRandSplit(a []int) {
	if len(a) < 10 {
		panic("length wrong")
	}
	mask := uint64(0xf)
	count := 0
	for {
		n := rand.Uint64()
		for i := 0; i < 16; i++ {
			sample := int(n & mask)
			n >>= 4
			// 不符合要求後直接跳到下一組去
			if sample >= 13 {
				continue
			}
			// 這裡也不需要mod
			a[count] = sample
			count++
			if count >= 10 {
				return
			}
		}
	}
}

如果一組不滿足取樣要求,我們就跳過直接去下一組,因此有可能16組裡無法獲得足夠的隨機數,因此我們得重新獲取一次64位的隨機數,然後再次進入分割計算。這麼做會對效能產生一點點負面影響,但依舊很快:

goos: linux
goarch: amd64
pkg: benchrand
cpu: Intel(R) Core(TM) i5-10200H CPU @ 2.40GHz
BenchmarkGetRand-8              16242730                72.22 ns/op            0 B/op          0 allocs/op
BenchmarkGetRandSplit-8         37794038                31.81 ns/op            0 B/op          0 allocs/op

這時候效能就只提升了一倍。

上面那種情況還是最簡單的,但[0, 7)就不好辦了。首先表示6需要至少3個bit,而3不能整除64,其次6也不是2的冪。這個怎麼處理呢?

有兩個辦法,核心思想都是拒絕取樣,但拒絕的範圍有區別。

第一個想法是,既然3不能整除64,那我們選個能被3整除的,這裡是63,也就是說超過2⁶³-1的樣本全部丟棄,然後把符合要求的樣本按每連續的3bits進行分割。這樣我們先保證了3bits分割出來的每一組都能均等的生成[0, 8)範圍內的隨機整數。現在問題轉化成了“rand8怎麼生成rand7”,這題我們會做而且做了好多回了,最終程式碼會是這樣:

func getRandSplit(a []int) {
	if len(a) < 10 {
		panic("length wrong")
	}
	// 注意,mask現在只要三位也就是0b111了
	mask := uint64(0x7)
	count := 0
	for {
		n := rand.Uint64()
		// 先拒絕大於2⁶³-1的樣本
		if n > 1<<63-1 {
			continue
		}
		for i := 0; i < 21; i++ {
			sample := int(n & mask)
			n >>= 3
			// 一組bits如果組合出來的數大於等於7也拒絕取樣
			if sample > 6 {
				continue
			}
			// 這裡是不用mod的,因為產生的sample本身只會在0-6之間
			a[count] = sample
			count++
			if count >= 10 {
				return
			}
		}
	}
}

程式碼變得很長也很複雜,而且需要兩步拒絕取樣,相對的我們一次也能分割出21組,比4bits的時候多了5組,所以難說效能下降還是不變,因此我們看看測試:

goos: linux
goarch: amd64
pkg: benchrand
cpu: Intel(R) Core(TM) i5-10200H CPU @ 2.40GHz
BenchmarkGetRand-8              16500700                73.77 ns/op            0 B/op          0 allocs/op
BenchmarkGetRandSplit-8         31098928                39.54 ns/op            0 B/op          0 allocs/op

確實慢了一些,但總體上還是提速了85%。

第二種想法只需要一步拒絕取樣,既然3不能整除64,那麼就找到一個離3最近的可以整除64且大於3的整數。在這裡我們可以直接注意到4符合條件,實際開發中如果要找到任意符合條件的數,可以依賴一下線性探測。現在我們按連續的4位把64位隨機整數分割,這樣分割出來的每一組可以生成均勻分佈在[0, 16)上的整數。然後問題變成了“從rand16生成rand7”。程式碼這樣寫:

func getRandSplit2(a []int) {
	if len(a) < 10 {
		panic("length wrong")
	}
	mask := uint64(0xf)
	count := 0
	for {
		n := rand.Uint64()
		for i := 0; i < 16; i++ {
			sample := int(n & mask)
			n >>= 4
			if sample > 13 {
				continue
			}
			// mod不能漏了,因為我們會產生大於等於7的結果
			a[count] = sample % 7
			count++
			if count >= 10 {
				return
			}
		}
	}
}

程式碼簡單了,也只需要一步拒絕取樣就行,但問題在於每一組的生成範圍變大導致我們不得不使用取模操作。看看效能怎麼樣:

goos: linux
goarch: amd64
pkg: benchrand
cpu: Intel(R) Core(TM) i5-10200H CPU @ 2.40GHz
BenchmarkGetRand-8              16451838                75.86 ns/op            0 B/op          0 allocs/op
BenchmarkGetRandSplit-8         30802065                39.15 ns/op            0 B/op          0 allocs/op
BenchmarkGetRandSplit2-8        38995390                30.75 ns/op            0 B/op          0 allocs/op

想法2比想法1快了近20%,看來兩步拒絕取樣成了硬傷,不僅僅是因為多獲取幾次64位隨機數更慢,多出來的一個if還可能會影響分支預測,即便最後我們多了5組可以取樣也無濟於事了。

所以,當你需要的隨機數範圍比較有限的時候,充分利用每一個bit是非常理想的效能提升手段。

帶有權重的隨機數

討論了半天機率均勻分佈的情況,但業務中還有一種常見場景:一組樣本進行取樣,既要有隨機性,又要樣本之間在統計上儘量滿足某些比例關係。

這個場景我想大家最先想到的應該是抽獎。是的沒錯,這是帶權重隨機數的常見應用。但還有一個場景,一個負載均衡器連線著一組權重不同的伺服器硬體,現在想要儘量按權重來分配連結,這時候帶權重隨機數就有用武之地了。

假設我們有樣本(1, 2, 3, 4)四個整數,然後要按1:2:3:4的比例來隨機生成樣本,該怎麼做呢?

按比例,我們可以得到1,2,3,4生成的機率是0.1,0.2,0.3,0.4,這些機率加起來是一定等於1的,所以我們不妨來想象有一個數軸上的0到1的區間,然後我們把這些比例“塞進”數軸裡:

0.0   0.1   0.2   0.3   0.4   0.5   0.6   0.7   0.8   0.9   1.0
|-----|-----|-----|-----|-----|-----|-----|-----|-----|-----|

|_______________________|
            4           
                        |_____|
                           1  
                              |_________________|
                                       3        
                                                |___________|
                                                      2

我故意打亂了順序,實際上順序並不影響結果。每個樣本可以有一個數軸上的範圍,範圍的長度之間也是符合比重的,因此當存在一個可以均勻生成[0, 1)之間所有實數的隨機數生成器時,這個生成器生成的數落在哪個範圍裡,我們就選擇生成這個範圍對應的樣本,舉個例子,如果生成的實數落在[0.0, 0.4)這個區間裡,那麼就生成樣本“4”,如果落在[0.8, 1.0)這個區間,就生成“2”。這樣帶權重的隨機數就生成了。這個看上面那個圖還是挺好理解的,我就不費筆墨去證明了。

但我不是很想用這個方法。為啥呢,因為你看到了,區間的左邊是閉合的,這意味著要做浮點數的等值比較,雖然很簡單,但我很懶不想多寫,而且浮點數沒法精確表示所有情況的比例導致我們區間的兩端都有精度損失,這就需要考慮誤差率,雖然通常這點精度損失帶來的誤差可以忽略不記(尤其是在統計意義上)但只是考慮到這點我就渾身難受了。

所以我要找一種不需要浮點數的解決方案。想出來其實不難,假設我們有0到9的整數,正好10個,現在樣本“1”按比例需要佔用其中1個樣本,樣本“4”按比例要佔用其中四個樣本,現在我們獲得了一個能均勻生成0到9的整數的隨機數生成器,那隻要生成的隨機數正好是樣本“4”佔用的那幾個隨機數我們就生成“4”,生成的隨機數是樣本“1”佔用的那就生成“1”。可以看到只要佔夠一定數量的不同的樣本,那麼我們一樣能生成帶權重的隨機數。

下面有幾個問題,一是樣本總數怎麼確定,這個簡單,每個比例當成整數相加即可,比如1:2:3:4就是1+2+3+4=102:2:3:5就是2+2+3+5=12,依此類推。如果比例是實數呢?2.3 : 3.4 : 4.7怎麼辦?這就要用到比例的性質了,等比擴大後比例不變,所以我們每個實數都乘以10,然後去掉小數點後的0全部當成整數,所以23+34+47=104,理論上任意比例都能這麼整,不過整數最終有大小限制的,你總不能生成個隨機數還用big.Int吧,所以注意總和別超過整數範圍限制。

二是樣本的範圍怎麼算,雖然我們只需要不相同的滿足1裡總數的離散的點就行,但為了方便計算我們還是選擇連續的整數比較好,所以範圍限定為[0, sum-1]。這樣我們能直接利用rand.Uint64N()來生成需要的隨機數生成器。

最後我們只要讓樣本按比例隨機佔領一些連續的整數就行了。而且我們只需要記錄右邊界就夠了,我們從範圍是[0, n]的第一個樣本開始比較,如果生成器給出的隨機數小於等於某個右邊界,那它一定落在邊界代表的樣本上(因為是從最小的邊界開始比較的,所以隨機數必然不可能落在前一個範圍的樣本上)。

其實就是把連續的實數換成了離散的整數點罷了,換湯不換藥。

搞明白思路後程式碼寫起來就是快:

type WeightRandom[T comparable] struct {
	entries       []entry[T]
	upperBoundary uint64
}

type entry[T comparable] struct {
	value T
	end   uint64
}

首先是資料結構,WeightRandom是我們的帶權重隨機數生成器,upperBoundary是樣本數量的總和,entries則是各個樣本和樣本佔領的連續整數的右邊界。

接著來看構造WeightRandom物件的方法:

func NewWeightRandom[T comparable](rndMap map[T]uint64) *WeightRandom[T] {
	var lowerBoundary uint64
	entries := make([]entry[T], 0, len(rndMap))
	for k, v := range rndMap {
		if v == 0 {
			panic("weight cannot be zero")
		}
		if lowerBoundary+v < lowerBoundary {
			panic("overflow")
		}
		lowerBoundary += v
		entries = append(entries, entry[T]{
			value: k,
			end:   lowerBoundary,
		})
	}
	slices.SortFunc(entries, func(a, b entry[T]) int {
		return cmp.Compare(a.end, b.end)
	})

	if len(entries) == 0 {
		panic("no valid sample")
	}

	return &WeightRandom[T]{
		entries:       entries,
		upperBoundary: lowerBoundary,
	}
}

lowerBoundary用來統計有多少樣本,我們最終選擇了左閉右開的區間,這樣方便算。rndMap的key是樣本,value則是比例。當樣本的範圍計算並儲存結束之後,我們需要按照右邊界從小到大排序這些樣本,因為後面的查詢範圍到樣本的對應關係需要右邊界滿足從小到大的順序。

最後是查詢函式:

func (w *WeightRandom[T]) RandomGet() T {
	x := rand.Uint64N(w.upperBoundary)
	for i := range w.entries {
		if x < w.entries[i].end {
			return w.entries[i].value
		}
	}
	panic("not possible")
}

查詢時先生成一個範圍在[0, upperBoundary)之間的隨機數,然後我們從最小的邊界開始逐一比較,一旦發現比自己大的邊界,那麼就說明需要生成邊界對應的樣本。底部那句panic如字面意思,理論上是執行不到的,但go不知道,我們要麼返回個空值要麼panic,考慮到能走到這裡那說明我們的程式或者標準庫的程式碼有重大bug,panic了去除錯才是比較好的辦法。

根據upperBoundary的大小,實際上我們還能複用上一節充分利用每一個bit的辦法,不需要每次都生成新的隨機數,等分割出來的分組都消耗完了再生成,這樣可以大大加速這個函式。不過為了通用性和儘量簡化程式碼,我就不這樣寫了。

最後附加一個用例:

func main() {
	w := NewWeightRandom(map[string]uint64{
		"a": 15,
		"b": 30,
		"c": 45,
		"d": 60,
	})
	m := make(map[string]int, 4)
	const limit = 100_0000_0000
	for range limit {
		m[w.RandomGet()]++
	}
	for k, v := range m {
		fmt.Printf("key: %s, count: %d, p: %g\n", k, v, float64(v)/limit)
	}
}

我們按權重生成“abcd”四個字母,比例是15:30:45:60,簡化一下就是1:2:3:4,所以理論上機率應該接近10%,20%,30%和40%。不過統計上的機率總是有點誤差的,只要大致趨勢接近於這個比例就行了。我們執行100億次來看看結果:

$ go run main.go

key: b, count: 2000011606, p: 0.2000011606
key: a, count: 1000058297, p: 0.1000058297
key: d, count: 3999943022, p: 0.3999943022
key: c, count: 2999987075, p: 0.2999987075

非常符合預期。作為一項最佳化措施,我們可以利用類似二分查詢的辦法來定位樣本,因為右邊界本身是有序的,這可以顯著改善在有大量樣本時的效能表現。不過如果你的樣本不超過10個的話我覺得現在這樣的線性查詢就足夠了。

怎麼說呢,簡單粗暴但有效解決了問題。只是和浮點數的方案比還是有缺點的:

  1. 因為整數有範圍限制,這導致了樣本總量會被限制在一個比浮點數方案更小的範圍內
  2. 同樣因為整數大小限制,比例的選擇範圍會比浮點數更小

因為限制更少,所以在通用的工具庫裡用浮點數方案的人更多,但業務場景和通用工具庫是不一樣的,很多時候選整數方案也沒啥問題,最終你應該選擇一個符合業務需求並且自己和同事都看得懂的方案。

至於效能怎麼樣,浮點數方案的查詢過程和整數方案是一樣的,效能需要做一次完整的測試才能看出孰高孰低,我不好在這憑空幻想。當然測試我就不做了,我偷個懶。

隨機數種子

“種子”其實是指一些偽隨機數生成演算法需要的一些初始狀態,這些演算法會根據初始狀態來生成一系列有隨機性的數值序列。

所以相同的種子通常會帶來相同的序列,這時候雖然每個序列都有隨機性,但兩個序列之間是沒有隨機性的——有了其中一個序列就可以精準預測另一個序列的排列。具體表現很可能會是你編寫了一個遊戲,其中敵人會隨機採取一些行動,然而因為種子沒設定好,導致每次見到這個敵人的時候它都會採取一模一樣的行動步驟,這樣的遊戲是極其無聊的。

不僅如此,產生相同的數值序列後還會帶來無數安全問題。

種子通常只用設定一次,並且在程式第一次需要隨機數的地方設定——理想情況是這樣的,然而總是有倒黴蛋忘記了這一點,於是隨機演算法經常只能使用預設的低質量且相同的種子。所以比較現代的語言,比如go1.22和python3都選擇了在程式剛開始執行的時候幫你自動設定種子。

此外我們還得擔心一下種子的質量。

“種子”的值需要每次獲取都不一樣,符合要求的常見種子來源有以下幾種:

  1. 使用系統時間。這個確實每次獲取都不一樣,獲取速度也很快,是很多開發者和庫的預設選項。但是系統時間是可以修改的,而且世界上還有閏秒這個麻煩東西,你意外得到相同的系統時間的機率其實不低。
  2. 使用隨機裝置,這些是軟體模擬出來的隨機數生成器,以Linux為例,有randomurandom兩種,其中random以作業系統及外部環境的噪音為資料來源產生隨機值,要求不高時我們可以認為這是“真隨機”,當然它的生成速率是比較低的;另一個是urandom,它會用外部噪音生成一個初始狀態,然後基於這個狀態用偽隨機數演算法快速產生隨機值,因此它的生成速率高但質量低。一般使用urandom生成種子是夠用的。
  3. 使用產生真實隨機數的硬體,原理和軟體模擬的類似,和大多數軟體實現轉硬體實現後會效能提升不同,TRNG反而可能會有更低的效能,但相比軟體實現它可以更精細地收集環境裡的雜音從而生成實驗證明的不可預測的高質量的隨機數。常見的TRNG生成速率從數百k/s到數兆/s的都有,但像/dev/random通常速率可以有數兆到數十兆/s。除了效能上的不穩定,不是所有裝置都包含TRNG,這導致了適用面受限,所以直接用它的人不多。不過很多依賴高質量隨機數的場景就不得不考慮TRNG了。
  4. 利用地址空間佈局隨機化。現代作業系統在載入程式後都會給程式的記憶體地址加點隨機的偏移量,所以程式每次執行的時候獲取的變數地址基本都是不同的。這個是成本極低的獲取隨機值的方法,幾乎不花費執行時代價,谷歌的abseil庫裡就有很多用這個方法獲取隨機數種子的程式碼。然而,使用這個方法的前提是你的系統要支援地址空間佈局隨機化,其次系統加的隨機偏移量質量要尚可,這兩個我們都控制不了,我們只能相信常用作業系統都做到這幾點了。另外,高許可權的程式和使用者始終能把一些程式碼寫進有固定地址的地方,雖然這種操作正變得越來越難,但還不是完全不可能,所以需要高質量種子的時候這個方案通常不會被考慮(另一個原因是有的系統可以配置關閉隨機化佈局甚至根本不支援隨機化)。
  5. auxiliary vector。Linux特有的,可以透過/proc/<pid>/auxv或者glibc的函式getauxval來獲取。這個陣列包含一系列作業系統和當前程序的資訊,全部是作業系統在程式載入時寫入的,Windows有類似的東西。這些資料中有些是不變的比如硬體資訊和平臺資訊,有些則是完全隨機的,比如其中有程式的入口地址和vDSO的地址,這些因為ASLR的緣故都是完全隨機的,另外auxv裡還有專門的隨機值欄位,這些資訊加一起比單純依賴ASLR能帶來更強的不可預測。

原則就是儘量讓預測結果的難度增加,最好是能做到完全不可預測。

那作為開發者我用啥呢?一般來說系統時間是夠用了,自己寫著完或者做些簡單工具可以用這個,不過要記住系統時間是可修改的不可靠的。如果是庫或者對隨機性依賴比較重的比如遊戲,/dev/urandom是個比較理想的選擇。追求極致效能並且對種子質量要求沒那麼高時,像谷歌那樣利用ASLR帶來的隨機值也是可以的。

實在有選擇困難症的話,我們來看看別人是怎麼做的。

golang和python3選擇了那些源作為種子

golang實際上是先利用auxv,如果系統不支援,就回退到從urandom之類的隨機裝置讀取隨機值,這個也出問題了就使用系統時間:

// runtime/rand.go

// OS-specific startup can set startupRand if the OS passes
// random data to the process at startup time.
// For example Linux passes 16 bytes in the auxv vector.
var startupRand []byte

func randinit() {
	lock(&globalRand.lock)
	if globalRand.init {
		fatal("randinit twice")
	}

	seed := &globalRand.seed
	// 檢視是否有auxv資訊被系統寫入
	if startupRand != nil {
		for i, c := range startupRand {
			seed[i%len(seed)] ^= c
		}
		clear(startupRand)
		startupRand = nil
	} else {
		// 先從urandom讀取
		if readRandom(seed[:]) != len(seed) {
			// readRandom should never fail, but if it does we'd rather
			// not make Go binaries completely unusable, so make up
			// some random data based on the current time.
			readRandomFailed = true
			readTimeRandom(seed[:])
		}
	}
	globalRand.state.Init(*seed)
	clear(seed[:])
	globalRand.init = true
	unlock(&globalRand.lock)
}

這個是全域性函式的設定,go還能自己建立rand.Source,這個的種子只能顯式傳進去,這時候傳什麼go就沒法管了,靈活的同時犧牲了一定的安全性。

Python3則是先讀取urandom,失敗後會結合系統時間加當前程序pid來生成種子,這樣比單使用系統時間要強:

// https://github.com/python/cpython/blob/main/Modules/_randommodule.c#L293
static int
random_seed(RandomObject *self, PyObject *arg)
{
    int result = -1;  /* guilty until proved innocent */
    PyObject *n = NULL;
    uint32_t *key = NULL;
    size_t bits, keyused;
    int res;

    // 引數為空的時候
	if (arg == NULL || arg == Py_None) {
       if (random_seed_urandom(self) < 0) {
            PyErr_Clear();

            /* Reading system entropy failed, fall back on the worst entropy:
               use the current time and process identifier. */
            if (random_seed_time_pid(self) < 0) {
                return -1;
            }
        }
        return 0;
    }

    // 引數不為空的時候根據引數生成種子

Done:
    Py_XDECREF(n);
    PyMem_Free(key);
    return result;
}

然後這個函式會被Random物件的__init__方法呼叫,如果初始化一個Random物件但不傳seed引數,那麼就會進行預設設定。而random模組裡所有的方法其實都是由一個全域性的Random()物件提供的,因為沒傳seed進去,所以程式碼裡會自動設定seed:

# https://github.com/python/cpython/blob/main/Lib/random.py#L924
_inst = Random()
seed = _inst.seed
random = _inst.random
uniform = _inst.uniform
triangular = _inst.triangular
randint = _inst.randint
choice = _inst.choice
randrange = _inst.randrange
sample = _inst.sample
shuffle = _inst.shuffle
choices = _inst.choices
normalvariate = _inst.normalvariate
lognormvariate = _inst.lognormvariate
expovariate = _inst.expovariate
vonmisesvariate = _inst.vonmisesvariate
gammavariate = _inst.gammavariate
gauss = _inst.gauss
betavariate = _inst.betavariate
binomialvariate = _inst.binomialvariate
paretovariate = _inst.paretovariate
weibullvariate = _inst.weibullvariate
getstate = _inst.getstate
setstate = _inst.setstate
getrandbits = _inst.getrandbits
randbytes = _inst.randbytes

這樣python就防止了使用者忘記設定seed的問題。

總結

關於隨機數要講的暫時就這麼多了,除了一丁點數值演算法之外都是些比較淺顯易懂的東西。

機率和統計對程式設計的影響是很大的,所以我覺得與其花時間看最近比較火的微分方程,稍微抽出點時間看看機率統計對自己的幫助可能更大。

最後,其實標準庫還有各種第三方庫已經貼心準備了幾乎全套功能了,看懂文件就能放心用,而且我也更推薦用這些庫,開源的且久經檢驗的程式碼始終是比自己閉門造車來的強。

參考

https://stackoverflow.com/questions/18394733/generating-a-random-number-between-1-7-by-rand5

https://en.wikipedia.org/wiki/Rejection_sampling

https://www.linuxquestions.org/questions/linux-kernel-70/what-does-proc-pid-auxv-mean-exactly-4175421876/

相關文章