モジュロ演算の替わりとなる高速処理

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

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

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

もし、232より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の値は(232 – 1)に到達し、0, 1, …, (232 – 1) mod Nの値はceil(232/N)回登場し、その他の値はfloor(232/N)回登場します。小さい値にバイアスがかかった公平なマッピングとなっています。

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

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

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

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

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

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

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

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

モジュロ演算 高速処理
8.1 2.2

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

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

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

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

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

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

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

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

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

関連記事

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