ある範囲に収まる乱数を得るために剰余（モジュロ）演算を書くとき、レビューするときに意識すること

はじめに

ある乱数生成器が N 個のセットのなかからランダムに一つを返すとき、その返り値をそれよりも小さな範囲に収まるようにしてから利用したい、という要件にたまに出会います。例えば、[0, 2^32) の範囲内の乱数を生成する乱数生成器を利用できる環境で、サイコロの目をランダムに計算するには、何らかの方法を使って [0, 6) の範囲の乱数に収める必要があります。このような getrandom(2) や /dev/urandom を使った乱数生成器の例以外にも、例えば Int64 のユーザー属性値を入力にしてユーザーを 10 種類に均等に分類したいという類の要件を過去にレビューしたこともあります。

ある値域をより小さい値域にマップするために、よく利用されるのは剰余（モジュロ）演算です。乱数生成器の例でいえば、その返り値を X とすると、 X % 6 を計算すれば結果は [0, 6) に収まります。

このモジュロ演算の使い方は、モジュロバイアス{modulo bias} と呼ばれる問題があることが知られています。ネットを検索しても仕事で書くときレビューするときにさっと参照しやすいページが見つけられなかったので、書くことにします。このエントリでは、モジュロバイアスとは何か、バイアスが発生しないケースとは何か、バイアスが発生する場合その程度をどう計算して評価するか、バイアスが発生するとしてそれを回避する具体的な実装はどうするべきか、説明します。

モジュロバイアスとは

より小さい値域にマップしたときに、期待しているのはその値の範囲でランダムな値を得ることです。例えば、[0, 2^32) の範囲の乱数 X を [0, 6) にマップするとき、その得られた値は 0, 1, 2, 3, 4, 5 のどれかになり、出現確率はそれぞれ同じ 1/6 であることが期待されます。最終的に欲しいのは [0, 6) の乱数だからです。

残念ながら、このとき X % 6 は 1/6 の均等な（uniform）出現確率でランダムな値を生成しません。わずかに 5 が出る確率が他（0, 1, 2, 3, 4）よりも少ないです。この均等でない出現確率（= 偏り）のことをモジュロバイアスと呼んでいます。

モジュロバイアスは、入力となる値域の最大数（上の例で言えば 2^32 のこと）と出力となる値域の最大数（上の例で言えば 6）の具体的な組み合わせ、そしてこれが発生するシチュエーションによって問題になる場合とならない場合が実務上はありえます。一律に許容できないわけではないので、理解した上で具体的なケースに適用する必要があります。

ここからは、わかりやすさのために、入力となる値域を小さくした例で説明していきます。

入力となる乱数の範囲: [0, 2^12=4096)
出力したい乱数の範囲: [0, 20)

この例では、17, 18, 19 が [0, 16] よりも出にくいバイアスがあります。

モジュロバイアスが発生しないケース

TL;DR: X % n = 0（X は入力の最大数、n は出力の最大数）が成り立つ場合はモジュロバイアスは発生しません。つまり、上の例で言えば、出力したい乱数の範囲が [0, 16) であれば 4096 % 16 = 0 が成り立つので均等な乱数を得られます。

モジュロ演算を使うと必ずバイアスが発生するわけではありません。それを確認するために、まずは {0, 1,..., 4095} までの値を順番に X % 20 した場合を考えてみます。{0, 1,...,19} までは均等にマップできます。その後 {20, 21,..., 4079} まで続けても均等なままです。残りの {4080, 4081,..., 4095} は {0, 1,..., 15} にマップされます。結果としてマップされなかった {16, 17, 18, 19} は出現頻度が少ないことがわかります。

それでは、{0, 1,..., 4095} までの値を順番に X % 16 した場合を考えてみます。{0, 1,...,15} までは均等にマップできます。 {16, 17,..., 4079} までも同様です。残りの {4080, 4081,..., 4096} は {0, 1,..., 15} にマップされます。余りは発生しないので出現頻度に偏りはありません。

別の言い方をすると、{0, 1,..., 4095} までの値を順番に X % 16 したとき、{0, 1,..., 15} のどれかにマップされる個数はそれぞれ 4096/16 で求められます。つまり、均等に 256 個ずつ振り分けられます。一方で、X % 20 のときは 4096/20 は割り切れないので、{0, 1,...,15} は ceil(4096/20) で 205 個ずつ、{16, 17, 18, 19} は floor(4096/20) で 204 個ずつ振り分けられます（合計すると 205 * 16 + 204 * 4 = 4096）。

モジュロバイアスが発生することがわかった場合は、それがそのアプリケーションにおいて許容できる程度かを判断する必要があります。

モジュロバイアスの程度を計算する方法

TL;DR: (X % n) / X * 100（X は入力の最大数、n は出力の最大数）で発生頻度に最大でどれくらい差が出るのかを計算できます。上の例で言えば、(4096 % 20) / 4096 * 100 は 0.390625 なので、[0, 16) と [16, 20) の間におよそ 0.4 % の偏りが発生します。具体的なイメージとしては、1000 回試行するうちの 4 回は不当な結果が出る程度と受け取れると思います。

{0, 1,..., 4095} までの値を順番に X % 20 した場合に偏りが出るのは、最後の {4080, 4081,...,4095} だけです。つまり、{16, 17, 18, 19) の発生頻度に比べて、{0, 1,..., 15} は 1 回多く出ます。

実験してみます。(rand(0...4096)) % 20 を 1 億回実行してそれぞれの結果とその頻度（ヒストグラム）を算出します。

rs = []
100_000_000.times { rs.push(((rand(0...4096)) % 20).round) }
rs.group_by{ |v| v }.each{ |k, v| puts "#{k},#{v.size}\n" }.map { }

その表とグラフが以下です。理想的にはそれぞれ 500 万回ずつ得ることですが、{16, 17, 18, 19} はそれよりも 2 万回前後明らかに出ていません。一方でそれ以外が当然多く出てます。

この程度のバイアスであれば許容できるかもしれませんが、(X % n) / X * 100 の値が大きくなる場合にはバイアスが許容できないほど大きくなります。上記の例では (4096 % 20) / 4096 * 100 は 0.390625 と小さいですが、20 を 2049(2^11+1) に変えるとバイアスは最大になります。このとき、(4096 % 2049) / 4096 * 100 は 49.9755859375 でほぼ 50 % です。これが意味するのは {0, 1,..., 2047} がそれぞれ 2 回ずつ発生するときに 2048 は 1 回しか発生しないということです。

バイアスが 50 % というのは 100 回試行するうちの 50 回は不当な結果が出る程度と受け取れるので許容できないと思います。また、セキュリティや基盤に関わる実装であれば多少のバイアスでも念のため避けたいというケースもあると思います。多くのプログラミング言語では、getrandom(2) や /dev/urandom を使って int32 や int64 の乱数を生成しつつ、Go であれば func (r *Rand) Int31n(n int32) int32、Ruby であれば rand(max) → number、Swift であれば func next<T>(upperBound: T) -> T のように、より小さな範囲を指定した乱数取得をサポートしています。これらの関数はモジュロバイアスを含まないように実装されています。

すでに用意された安全な関数がない場合、自分でモジュロバイアスを回避する実装を行う必要があります。

モジュロバイアスを回避する実装例

具体的な値によって取りうる回避策が異なるので、簡単なものから順に場合分けしつつ説明します。
ここでは、代表的な例として二つ挙げます。以下の X は入力の最大数、n は出力の最大数を指します。

1. X % n = 0 になるように X または n を調整する
2. 入力となる乱数を実際に生成した後に、それがバイアスになり得る値の範囲であれば生成し直す

X % n = 0 になるように X または n を調整する

もっともシンプルな解決方法ですが、当然 X または n を調整できない要件の場合には採用できません。

X が 2 の累乗の場合は、n がそれよりも小さい 2 の累乗であれば割り切れます。例えば、2^12（4096） は 2^1（2） でも 2^2（4） でも 2^11（2048） でも割り切れます。

このような X も n も 2 の累乗であれば除算を伴うモジュロ演算を使わなくても、(n - 1) をビットマスクとして使えます。パフォーマンスを意識した実装方法として目にする機会があるかもしれないので紹介しておきます。Go の Int31n 関数では X は 2^32 なので、n int32 が 2 の累乗であれば、すぐに ビット演算 AND を適用して結果を返しています。例えば、4096 & 15 は二進数だと 1000000000000 & 1111 なので 0 です。4095（111111111111） & 1111 は 15 です。左の数字を 1 ずつ減らしていくと結果も 1 ずつ減って 0（= 4080 & 15） になったらまた 15（= 4079 & 15） に戻ります。ちなみに、n&(n-1) == 0 は例えば n=16 に当てはめると 10000&1111 になるので成立します。

// Int31n returns, as an int32, a non-negative pseudo-random number in [0,n).
// It panics if n <= 0.
func (r *Rand) Int31n(n int32) int32 {
...
	if n&(n-1) == 0 { // n is power of two, can mask
		return r.Int31() & (n - 1)
	}
...

バイアスになり得る乱数を生成したら、再度生成し直す

入力となる乱数を生成し直すことが許容できる場合には、これが採用されます。X=4096, n=20 の場合で言えば、入力となる乱数が 4080 以上であれば捨てて、それ未満が生成されるまで繰り返すということです。n=2049 になると半分捨てることになるので無駄が多くなります。

以下は Go の Int31n 関数の実装続きです。読み解くために X は 4096(2^12) で n は 20 とすると、max は (1 << 12) - 1 - (1 << 12)%20) で要するに 4096 - 1 - 16 = 4079 です。その後は 4079 以下の乱数が生成されるまでループに入って、抜けた後モジュロ演算を適用して結果を返します。

// Int31n returns, as an int32, a non-negative pseudo-random number in [0,n).
// It panics if n <= 0.
func (r *Rand) Int31n(n int32) int32 {
...
	max := int32((1 << 31) - 1 - (1<<31)%uint32(n))
	v := r.Int31()
	for v > max {
		v = r.Int31()
	}
	return v % n
}

まとめ

乱数にモジュロ演算を適用してより範囲を絞った乱数を得ようとするとき、モジュロバイアスが発生する可能性があります。モジュロバイアスが許容できない場合は、回避する実装を検討する必要があります。乱数に関しては、プログラミング言語の標準関数などで範囲指定可能な関数が用意されていることが多いので、バイアスを回避するように実装されていることを念のため確認の上、それを積極的に使うのが良いでしょう。

このエントリ内の多くの例は、もとになる乱数を 2^12 と小さい値にしています。そのため過度な印象を与えるかもしれないので最後に補足しておきます。多くの乱数生成器が返す int32 や int64 で発生するバイアスは計算すれば分かる通り、特殊な用途でない限り無視できる程度がほとんどなので、過剰に気を使う必要はありません。

参考資料

• A fast alternative to the modulo reduction – Daniel Lemire's blog