N 個の要素セットの中からランダムに整数を1つ取り出すと仮定します。使っているコンピュータに32ビット整数の乱数を生成する機能がある場合、その数をどのように N より小さいインデックスに変換すればいいのでしょうか。例えば、サイズが N のハッシュテーブルがあると仮定します。このような場合でも、ハッシュ値（通常32ビットや64ビットの整数）を N より小さいインデックスに変換する必要があります。このような場合、大抵プログラマは解決の際、 N が2のべき乗であるようにしますが、これは必ずしも理想的とは言えません。

任意の整数 N をできるだけ公平にマッピングしたいとします。2 ³² 個存在する32ビットの整数全てから始める場合、{0, 1 ,…, N – 1}の値域に定められた値に、ちょうど2 ³² / N 個の値をマッピングできるようにするのが理想的です。

残念ながら、2 ³² は N によって割り切れないため、公平にマッピングすることはできません。しかし、完璧にでなくても、floor(2 ³² /N)かceil(2 ³² / N )のどちらかの個数の値が、値域のそれぞれの値にマッピングされるように設定することはできます。

もし、2 ³² より N が小さい場合、上記の方法で完璧に近いマッピングができるでしょう。

これらを解決する際、一般的に使用されるのがモジュロ演算で、 x mod N と表現されます（特記がない限り、コンピュータサイエンスにおいて、モジュロ演算とは除算により求めた余りと定義されています）。

uint32_t reduce(uint32_t x, uint32_t N) {
  return x % N;
}

公平にできたかの確認はどのようにすればいいのでしょうか。 x の値を0から始めてみましょう。すると、 x の値が増えるにつれ、モジュロ演算の値は、0, 1, …, N – 1, 0, 1, …などをとることがわかります。そして、最終的に x の値は(2 ³² – 1)に到達し、0, 1, …, (2 ³² – 1) mod N の値はceil(2 ³² / N )回登場し、その他の値はfloor(2 ³² / N )回登場します。小さい値にバイアスがかかった公平なマッピングとなっています。

これでもいいのでしょうが、モジュロ演算は除算を行い、除算は安価な演算ではありません。乗算より高くなります。32ビットの除算を最近のx64プロセッサで実行すると、命令スループットは1命令毎6サイクルで、レイテンシは26サイクルとなります。対照的に、乗算の場合、1命令毎1サイクルで、レイテンシは3サイクルとなります。

N が事前に分かっていれば、工夫を凝らした方法でモジュロ演算を複数の乗算や他の演算に変身させ、”事前計算”することができます。コンパイル時に N が分かっていれば、コンパイラはこれらの演算を事前に実行してくれます。そうでなければ、ソフトウェアライブラリあるいは独自の公式を使って事前計算することができます。

しかし、もっといい方法が実はあるのです。それは、パフォーマンスに影響せず、簡単に実装でき、モジュロ演算で提供されるマップと同等のものを出してくれます。

まず、 x と N を32ビットの整数と想定し、64ビットの x * N を見てみましょう。( x * N ) div 2 ³² は値域内に収まっており、公平にマップされています。

uint32_t reduce(uint32_t x, uint32_t N) {
  return ((uint64_t) x * (uint64_t) N) >> 32 ;
}

( x * N ) div 2 ³² の演算は64ビットのプロセッサを使えば、非常に高速にできます。この計算は、乗算をした後ビットをシフトしています。最近のIntelプロセッサでは、恐らくレイテンシは約4サイクルで、スループットは最低でもコール上は毎2サイクルだと思います。

では、32ビットのモジュロ演算に比べてどれだけマッピングは早くなるのでしょうか。

これはテストするために、サイズ N の配列のインデックスにランダムに繰り返しアクセスするようベンチマークを実装しました。インデックスはモジュロ演算か上記の方法で取得できます。最近のIntelプロセッサ（Skylake）を使うと、毎アクセスでのCPUサイクルは次のとおりになります。

結果は4倍速となります。悪くないです。

いつもどおり、 私のコードは自由に利用できるようになっています 。

これは何に役に立つのか。例えば、除算を高価にしないため、今まで配列とハッシュテーブルに2のべき乗の容量を無理やり持たせていた場合、上記の高速処理マップを使うことで、パフォーマンスを犠牲にすることなく、任意の容量をサポートすることができます。さらに、乱数の生成も速くなり、高速の乱数発生器を使っている場合は、大いに役に立ちます。

では、マップが公平かをどうのように確認すればいいのでしょうか。

N を掛けることで[0, 2 ³² )の整数値を[0, N * 2 ³² )の N の倍数にマッピングします。2 ³² で割ることで、[0, 2 ³² )にある N の倍数全てを0にマッピングし、[2 ³² , 2 * 2 ³² )にある N の倍数を1にマッピングする……といったことができます。公平にマッピングされているか確認するには、長さ2 ³² の区間での N の倍数の個数をカウントすればいいのです。カウントは、ceil(2 ³² / N )かfloor(2 ³² / N )になるべきです。

例えば、「区間の最初の値が N の倍数」という、区間中に存在する倍数の数が最大であることが明白なシナリオを仮定します。いくつあるでしょう。ceil(2 ³² / N )となります。長さ N の部分区間を設定すると、全ての完全な区間は N の倍数で始まり、剰余がある場合は N の倍数が1つ追加されます。最悪のシナリオでは、区間内で最初 N の倍数が現れるのは N – 1の位置になります。この場合、倍数はfloor(2 ³² / N )個になります。ここでも、なぜこのような結果になるのかを見る場合は、長さ N の部分区間を設定します。全ての完全な区間は N の倍数で終わります。

これで、マップが公平であることを証明できました。

では、もう少し正確にしてみましょう。「 N の倍数の数を最大化できるのは、長さ2 ³² の区間の最初の値が N の倍数である場合である」と前に述べました。その区間の最後には、不完全な長さ2 ³² mod N の区間ができることになります。しかし、 N の倍数が区間の最初に現れるようにするのではなく、インデックス2 ³² mod N の位置に現れるようにすれば、不完全な部分区間ができることはありません。つまり、2 ³² mod N より前に N が現れると、ceil(2 ³² / N )個の倍数が得られ、そうでない場合はfloor(2 ³² / N )個の倍数を得られるのです。

どの結果がfloor(2 ³² / N )の頻度で現れ、どの結果がceil(2 ³² / N )の頻度で現れるのかを知る術はあるのでしょうか。あります。 k という値を出力したと想定します。まず、 N の倍数のうちk 2 ³² 以上のものがどこにあるのか探す必要があります。これはceil(k 2 ³² / N ) N – k 2 ³² にあるのですが、これを2 ³² mod N と比較するだけです。もし、2 ³² mod N よりも小さければ、カウントはceil(2 ³² / N )になり、そうでない場合のカウントはfloor(2 ³² / N )となります。

関連記事

• math/rand: speed up Int31n with multiply/shift instead of modulo (golang issue 16213)

• Agner Fog, Pseudo-Random Number Generators for Vector Processors and Multicore Processors 、Journal of Modern Applied Statistical Methods、2015年。

（更新：Kendall Willetsのコメントを受け、証明をより直感的なものにしました）