爆速な素数判定の怪

本記事では，プログラムの高速化を競うプラットフォームである highload.fun における，素数判定の高速化を競う問題で，驚異的な速さを誇る解答が現れたことについて紹介します．

なお，highload.fun では他の人の解答が見られないため，その爆速な解答の中身について直接触れることはしません（というかできません）．そのため，本記事は，爆速な解答が現れたという事実や事の経緯をまとめるにとどまる，いわば謎を謎のまま終わらせる記事であることを先に記しておきます．

highload.fun について

著者（@matsuoka-601）は，highload.fun という，プログラムの高速化を競うプラットフォームで遊ぶことを趣味としています．

highload.fun とは，「解くこと自体は簡単だが，高速化の余地が大いにある」ようなタスクについて，どれだけ高速化できるかを競うコンペティションサイトです．たとえば，以下のようなタスクがあります．

• ASCII エンコードされた 50,000,000 個の u32 を足し合わせた値を求める（Parse integers）
• 250,000,000 個の u32 をフォーマットして，フォーマット結果のチェックサムを求める（Format integers）
• 30,000,000 個の文字列の重複除去を行い，ユニークな文字列の個数を求める（Unique Strings）
• 15,000,000 個の UUID をソートし，出力する（Sort UUIDs）

これらのタスクは，いずれも解くこと自体は簡単であり，正解を計算するだけならたかだか 10 行くらいのコードを書けば事足ります．しかしながら，順位はプログラムの実行時間で決定されるため，上位に食い込めるかどうかは，その解答をどれだけ高速化できるかにかかっています．

以下は，もっとも挑戦者数の多い Parse integers のランキングです．

この問題では，「ns単位で測った実行時間を 10,000 で割った値」がスコアになります．つまり，1 位の解答は約 29ms で 50,000,000 個の u32 をパースしていることになります．一方，ナイーブな C++ の解答は約 13 秒（=13,000 ms） かかるので，1 位の解答はナイーブな解答と比べて約 450 倍も高速化してあることになります．

ちなみに，この問題のトップ層の解答は，ジャッジサーバのメモリの Sequential Read の帯域幅にかなり近いものになっています（約 17GB/s 前後）．すなわち，ハードウェアの制約上，これ以上高速化できないところまできているということです．

解答は一部の問題を除きシングルコアで実行されるため，上位陣は SIMD の高速化バトルになることが多いようです．ジャッジサーバの CPU は Haswell というひと昔前のチップですが，AVX2 を使うことができるのでそれなりに面白い戦いになります．また，期限のないコンペティションなので，1 位の記録が数年越しに破られることもしばしばあり，それも醍醐味の一つになっています．

Task: Sum of Prime Numbers

highload.fun のタスクの一つに，Sum of prime numbers というものがあります．このタスクも，以下のような非常にシンプルなものです．

1,000,000 個の u32 が与えられる．そのうち素数であるようなものの和を求めよ．（筆者注 : 与えられる整数は [0, 2^32) の範囲の一様乱数のようです）

要は素数判定の高速化を競う問題です。この問題の順位表のトップは，長らく（1 年くらい？）15ms 前後で推移してきており，数か月前に 13ms の解答が出て少し話題になりました．

つい先日までの順位表（Parse integers とはスコアの計算方法が異なります）

そのため，highload.fun のコミュニティの中でも，「このあたりのスコアが限界値だろう」という認識があったように思います．

zielaj 現る

しかし，つい先日，いままでの全ての記録を過去のものにしてしまうようなとんでもない解答が，zielaj というユーザーによって提出されました．

2 位の解法よりも 4 倍近く速い，3.6ms という記録です．

highload.fun コミュニティの反応

この解答の出現に highload.fun コミュニティは大いに沸き立ち，どのようなアルゴリズムを用いているのかを予想する議論が交わされました．

この問題に対するこれまでの上位陣の解法は，まず小さい素数で試し割りをして，素数でないものを除外し，残った数に対して 32bit 整数向けの高速なミラーラビン素数判定を行う，というものだと思われます（少なくとも，自分の 15ms の解答はそのようなものです）．ミラーラビン素数判定は結構重いので，最初に試し割りでできる限り素数の候補を削っておくのが重要になります．私の解答の試し割りは AVX2 でバキバキに高速化してあり，ボトルネックはミラーラビン素数判定の部分です．

しかし，この方針の解答では，3.6ms というスコアは達成不可能だという予想が大方であったため，全く違う手法を用いているのではないかという予想が立ちました．その予想として真っ先に上がったのが，コンパイル時に素数か否かのテーブルを構築し，実行時にはそのテーブルを参照して素数判定をするという手法です．

訳：わぁ，zielaj はコンパイル時に素数（かどうかのテーブル）を前計算することに成功したのかもしれない．

しかし，正しく素数判定をするためには，32bit のすべての符号なし整数について，それが素数かどうかをコンパイル時に前計算する必要があります（小さい素数の倍数を削って計算量やメモリ消費をある程度減らすことはできますが）．このコストは，エラトステネスの篩を高速化してもなお非常に大きなコストであり，ジャッジサーバのコンパイル時間制限に引っかかってしまいます．実際，何人かのユーザーがこの方法を試したみたいですが，いずれも失敗に終わりました．

そんな中，zielaj が自分の解法について発言しました．

意訳 : 私も前計算する良い方法が見つけられなかったので，私の解法は実行時の素数判定を行っています．

コンパイル時の前計算は行わず，実行時に素数判定をしているとのことです．この発言により，ますます謎は深まりました．実行時の素数判定をする場合には，基本的にはミラーラビン素数判定やそれに類するアルゴリズムを使う必要があるはずですが，それらのアルゴリズムをここまで高速化することが到底不可能であるように思われるためです．

最後に

ここまでが現在までの経緯であり，記事はここで終わります（尻切れトンボですが・・・）．自分も zielaj の記録に迫るための突破口を模索していますが，まだ見つけられていない状況です．コンペティションは無期限につづいているので，高速化の腕に覚えのある方は zielaj の記録に挑戦してみてはいかがでしょうか．

記事の内容について質問があったりした場合には，お気軽に著者（@matsuoka-601）までどうぞ．