ある値より小さい乱数を得る - arc4random_uniform(3)

1つ前の記事に関連して、OpenBSDのarc4random(3)やarc4random_uniform(3)について調べていたのでメモ。

BSD系のlibcを用いているOS（OpenBSD・FreeBSD・OS X・iOSなど）で提供されている*1arc4random(3)は、[0, UINT_MAX)の範囲でランダムな値を返します。これを用いて、特定の範囲の乱数、ここでは、ある値より小さい乱数を得ることを考えます。

剰余を計算する

たとえば[0, 5)の範囲の整数である乱数を得たい場合に、お手軽に

int random_number = arc4random () % 5;

と剰余を計算する方法があります。しかし、この方法では生成される乱数に偏りができてしまいます。これはarc4random(3)自体の問題ではありません（後述）。

[0, 10)の範囲の整数である乱数を生成する関数rand10 ()があるとします。
この関数を用いて、3より小さい整数である乱数（[0, 3): 0・1・2のいずれか）を生成するとき、上記の方法を用いてみます。

rand10 () rand10 () mod 3
0 0
1 1
2 2
3 0
4 1
5 2
6 0
7 1
8 2
9 0

rand10 ()が、0〜9のいずれかの整数をランダムに偏りなく返したとしても、3の剰余を計算すると、0の出現率は4/10、1の出現率は3/10、2の出現率は3/10と、0が1や2よりも多く出現します。つまり、得られた乱数に偏りができてしまいます。

OpenBSDのarc4random(3)のmanpageによると、剰余を用いたときに偏りが出るこの現象を"modulo bias"と呼ぶそうです。

なお、剰余計算では下位のビットしか使われないため、線形合同法のような下位ビットの質が悪い乱数生成器を用いた場合にはさらに偏りができてしまいます。

0以上1未満の範囲にしてn倍する

[0, 5)の範囲の整数である乱数を得たい場合に、

int random_number = ((double) arc4random () / UINT32_MAX+1.0) * 5;

のように、0以上1未満の範囲の実数にした上で5倍する方法があります。

ステップを分解してみます。

  1. uint32_tの範囲の乱数を生成する: [0, UINT32_MAX]
  2. UINT32_MAX+1.0で割る: [0.0, 1.0)
  3. 5倍する: [0.0, 5.0)
  4. 小数点以下を切り捨てて整数にする

この方法でも、得られる乱数に偏りがあります。

再びrand10 ()を用いて、3より小さい整数である乱数（[0, 3): 0・1・2のいずれか）を生成する場合を考えてみます。

rand10 () rand10 () / 10.0 rand10 () / 10.0 * 3 floor (rand10 () / 10.0 * 3)
0 0.0 0.0 0
1 0.1 0.3 0
2 0.2 0.6 0
3 0.3 0.9 0
4 0.4 1.2 1
5 0.5 1.5 1
6 0.6 1.8 1
7 0.7 2.1 2
8 0.8 2.4 2
9 0.9 2.7 2

1・2は3回ですが、0は4回出現しており、偏りができています。
この偏りはrand10 ()が生成する値の範囲を広くすれば相対的に小さくなると考えられますが、それでも一様であるとはいえなくなります。

なお、この方法では上位のビットを使うため、下位ビットの質が悪い乱数生成器を用いた場合でも剰余を用いた方法のような偏りは生じません。

arc4random_uniform(3)を使う

arc4random(3)とともに、ある値より小さい乱数を返すarc4random_uniform(3)という関数も提供されています。

arc4random_uniform(3)は引数として生成される乱数のupper boundを取ります。そして、0以上で指定された値より小さい整数である乱数を返します。

arc4random_uniform(3)も内部ではarc4random(3)を呼び出して乱数を生成し、先に説明した剰余を用いる方法で、指定された値より小さい乱数にしています。しかし、modulo biasの影響がないように、乱数をあらかじめ選別しています。

選別する方法ですが、ソースコードの説明には、[0, 232 mod upper bound)の範囲から外れる乱数が生成されるまで乱数の生成を繰り返すことで、[232 mod upper bound, 232)の範囲の乱数を得る、とあります。説明だと何のことだかさっぱり分かりませんが、mod upper boundされた状態でUINT_MAXから数えて[0, upper bound）の列がぴったりn回出現するように、範囲を制限しているようです。

イメージとしては、0からUINT_MAXまでの数字が書かれた、厚みのないセロハンテープのような細長い紙を筒に最後まで巻いて、内側の切れ端と外側の切れ端が合わない余った部分を内側から切り落とす感じになると思います。

下からではなく上から切っていく方法でもよいと思いますが、この場合の選別する範囲を算出するには、mod upper boundが0になる最大の値を探す必要があります。これはfloor (232/upper bound) * upper boundで計算できます。しかし、232というuint32_tに格納できない値が出てきている点がいけてないです。

範囲から外れた値がarc4random(3)から返ってきた場合には再度arc4random(3)を呼び出している*2ため、運が悪いと何度もarc4random(3)を呼び出すことになりそうです。最良の場合の一例（upper bound = 2; 232 mod 2 = 0）では範囲外となるのは[0, 0)であり範囲から外れる値はありません。最悪の場合（upper bound = 231+1; 232 mod (231+1) = 2147483647）では範囲外となるのは[0, 2147483647)と約50%の確率ですが、よほど運が悪くない限りは多くの試行を必要としないと思われます。

再びrand10 ()を用いて、3より小さい整数である乱数（[0, 3): 0・1・2のいずれか）を生成することを考えてみます。232 mod 3 = 1ですので、範囲チェックとして[0, 1)であればその値を不採用とします。

rand10 () 範囲チェック rand10 () % 3
0 × -
1 ○ 1
2 ○ 2
3 ○ 0
4 ○ 1
5 ○ 2
6 ○ 0
7 ○ 1
8 ○ 2
9 ○ 0

このように、0・1・2が均等に分布していることが分かります。

先に述べましたが、arc4random_uniform(3)も、内部ではarc4random(3)を呼び出しています*3。arc4random(3)の返す乱数が良くないというよりも、その結果を加工する方法が良くないということです。

なおこの方法は乱数に対して剰余計算を行っているため、下位ビットの質が悪い乱数生成器と組み合わせることはできません。

*1:Linuxでも、たとえばDebianではlibbsdライブラリとして提供されています

*2:ソースコードでもfor (;;)を用いた無限ループになっています

*3:libc/crypt/arc4random_uniform.c