Library CheckerのFactorizeという問題でfastestを取ったので、そのコードの解説です（速度ランキング）。 Factorizeは、 $10^{18}$ までの整数が高々100個与えられ、それらを素因数分解せよ、という問題です。

概要

上位陣（13ms～15ms）の解法では、ECM（楕円曲線を用いた素因数分解法）が使われているようですが、よくわからないのでこれには手を出しません。 ECMを用いていない提出の中で最速の提出はhqztrue氏の16msであり、20%以上の差があります。この提出では、ポラードのロー法という手軽に実装できるアルゴリズムが用いられていました。これを効率化してみます。

ポラードのロー法には、乱数を連鎖的に生成する部分が存在します。この部分がレイテンシに支配されていて並列性に乏しいという問題点に着目しました。これに対し、異なる二つの系列の乱数列を同時に生成する手法を導入します。これにより、プログラムの並列度をスループットに支配されるところまで上げることができます。並列化のオーバーヘッドがかなりあり、二倍速くなるというわけではないですが、それでも12msと上位陣の解法を上回る速度を達成できました。

基礎知識

純粋にアルゴリズムの基礎知識です。 CPUの特性を生かした高速化に関する記述はほぼないので、高速化テクニックが知りたい場合は「本題：私の提出コードの解説」まで読み飛ばしても多分大丈夫です。ポラードのロー法を“貼って”いるけれど、何をやっているのかあまりよく知らない、という場合は読んでみると面白いかもしれません。

ポラードのロー法

ポラードのロー法は、合成数 $n$ を入力すると、その非自明な因子（ $n$ の $1$ や $n$ でない約数） $a$ を返す乱択アルゴリズム（ラスベガス法）です。ここで、返ってくる $a$ は素数とは限らない点に注意してください。とはいえ、素数判定と組み合わせ、素数でなければ再帰的にロー法を適用する、とすることで素因数分解を行うことができます。

乱数を使って素因数分解するというのはどういうことでしょうか。単純な方法として、 $[0, n)$ の範囲の整数を一様ランダムに生成し、それと $n$ の最大公約数（greatest common divisor, GCD）を求める、という方法があります。これを繰り返していれば、そのうち $a$ が返ってくることがあるでしょう。しかし、これには平均で $a$ に比例した回数の乱数生成が必要です。

誕生日のパラドックスを考えると、次のような方法で乱数生成回数を減らせます。 $[0, n)$ の範囲の整数をランダムに $\sqrt a$ 個程度生成し、それら同士の差を計算します。ここで計算された「差」一つ一つに対して、 $n$ との間の最大公約数を求めれば、乱数生成は $\sqrt a$ に比例した回数で済みます。しかし、乱数を二つ取り出す組み合わせの数を考えると、「差」は $a/2$ 個ほどあるので、計算量はやはり $a$ に比例してしまいます。なんとか無意味な組み合わせの「差」を求めることを回避できないでしょうか。

これに対してはうまい方法があり、それがポラードのロー法です。ここまでは一様ランダムな真の乱数を考えていましたが、そうではなく、あえて疑似乱数を使うというものです。例えば、なんらかの整数 $C$ に依存した疑似乱数系列 $x_i \mapsto x_i^2 + C \bmod n$ がよく用いられます。‘

「差」と $n$ の間の最大公約数を求めたときに $a$ が返ってくるのは、「差」を取った二つの乱数が $\bmod a$ で一致していた場合なので、乱数の $a$ での剰余に着目します。いきなり $a$ での剰余に着目してもいいのですが、それだと当たり前のことを言っているだけのようにも感じてしまうかもしれません。そこで、あえて任意の $K$ での剰余の場合にどうなるかを考え、その後 $K = a$ という特殊な状況ではどのように性質が変わるかを見てみます。

さて、この疑似乱数系列を $\mod\! K$ で見た $y_i = x_i \bmod K$ を考えます。この時、 $y_i \mapsto y_{i+1}$ という変換は一意に決まる（写像になる）でしょうか？ $K$ がなんでもない数であれば、これは一意に決まるとは限りません。 $x_i$ は $y_i$ よりも多くの情報を持っているので、 $y_i$ が決まっても $y_{i+1}$ には複数の可能性があります。

例で考えてみます。 $x_i \mapsto x_i^2 + 1\bmod 15$ という疑似乱数生成式は、 $x_1 = 1$ に対して1→2→5→11→2→5→11→……という疑似乱数系列を生み出します。ここで、 $K = 6$ で割った余りを考えると、1→2→5→5→2→5→5→……となっています。つまり、5の次が5だったり2だったりしています。 $y_i \mapsto y_{i+1}$ という変換が一意でないということです。

しかし、 $K$ が $n$ の約数の場合は状況が違います。たしかに $x_i$ は $y_i$ より多くの情報を持っていますが、それは $y_i$ と $y_{i+1}$ の関係には影響を及ぼしません。なぜなら、 $y_i \mapsto y_i^2 + C \bmod K$ が成り立つからです（ $\bmod rs$ で成り立つ式は $\bmod r$ でも成り立つ、という性質を使いました）。よって、 $K$ が $n$ の約数の場合は変換が一意です。変換が一意なので、系列中に同じ値が出現すると、以降は全く同じ系列が繰り返されることになります（上記例では5→5となったからといって、その後も5→5→5→5→……となるわけではないのと対比してみてください）。以下では、系列が繰り返される周期を $P$ と表記します。

話を戻すと、乱数 $x_i$ と乱数 $x_j$ の差 $x_i - x_j$ を求め $n$ との最大公約数を求める手順で、無意味な組み合わせを回避したいという話でした。無意味な組み合わせとは、 $x_i - x_j$ が $a$ の倍数にならないことが明らかな組み合わせのことです。これは、 $K=a$ で考えれば、 $y_i$ と $y_j$ が必ず異なることが明らかな組み合わせのことです。ポラードのロー法では、疑似乱数系列に周期があることをうまく活用し、無駄な組み合わせを回避します。例えば、 $y_{i+1}$ と $y_{i+3}$ が異なっていれば、 $y_i$ と $y_{i+2}$ も異なっています（対偶を考えれば明らかです）。よって、 $x_{i+1}-x_{i+3}$ と $n$ の最大公約数が1だった場合、 $x_i-x_{i+2}$ と $n$ の最大公約数を求める必要はありません。

この考え方、つまり添え字の差が同じペアは調べても無駄、という事実を活用すると、調べるべきペアの個数を $P$ 個程度に減らすことができます。 $x_i$ と $x_{i+1}$ の差と $n$ の最大公約数を求める、 $x_i$ と $x_{i+2}$ の差、 $x_i$ と $x_{i+3}$ の差、 $x_i$ と $x_{i+4}$ の差、……とP回繰り返せばよいからです（※）。乱数を $P$ 個生成した時、ペアは $P^2/2$ 個程度ありますが、調べる必要があるのは $P$ 個程度でしかない、ということです。

周期 $P$ は $K=a$ 以下*1であることが明らかですが、 $y_{i+1} = x_{i+1} \bmod K = y_i^2 + C \bmod K$ という疑似乱数生成式を使った場合、周期 $P$ は $\sqrt K=\sqrt a$ 程度になることが経験的に知られています。 $n$ の任意の因子 $c$ について、 $\sqrt c$ 個程度のペアを調べると $c$ が見つかるということなので、基本的には小さい因子から見つかっていきます。 $n$ の最大の素因子はこの方法で見つけなくていい（ $n$ を見つかった因子で割るということを繰り返して残った数が素数かを判定すればよい）ので、 $n$ の二番目に大きな素因子の大きさの正の平方根に比例するくらいの時間がかかるとみてよいでしょう。 $n$ の二番目に大きな素因子の大きさは高々 $\sqrt n$ なので、 $n^\frac14$ に比例した回数の乱数生成により素因数分解できることが経験的に知られています。

※実際には、 $i$ をどう選ぶかという問題があります。 $y_1$ が循環に含まれるとは限らないからです（例えば、上で紹介した1→2→5→11→2→5→11→……の例では、1は循環に含まれていません）。これは、適宜 $i$ 自体を進めていく、という方法で解決可能です。ポラードが提案した進め方（二回に一回進める）よりも乱数生成回数が少なくなる進め方（ $2^t$ 回に一回、系列を $2^{t-1}$ だけ進める）をブレントが提案しています。

参考

以下の記事に、アニメーションを用いたわかりやすい解説があります。

ポラード・ロー素因数分解法について #競技プログラミング - Qiita

モンゴメリ乗算

$\mod\!n$ の計算は、普通にやると除算命令が必要です。除算命令は一般に遅いのでどうにか回避したいです。

モンゴメリ乗算は、バレット還元と同様に、いくらかの事前計算を行うことで、剰余演算を除算命令無しに行う方法です。モンゴメリ乗算では、モンゴメリ表現と呼ばれる数値表現で値を保持したまま計算が進みます。モンゴメリ表現は、 $\bmod n$ の代表元 $X$ の代わりに、 $2^{64}$ などの切りのいい数を掛けた $2^{64}X \bmod n$ を使います。

事前計算で求めるものは、モンゴメリ表現から通常の表現に戻すときに使う値 $R$ です（他にも用途ごとに他の数を事前計算することがあります）。これは $\bmod 2^{64}$ での $n$ での逆元であり、 $nR \equiv 1 \mod 2^{64}$ を解けば求まります。この $R$ を用いると、乗算を以下のように行うことができます。

元の表現	モンゴメリ表現
$X$	$A = 2^{64}X \bmod n$
$Y$	$B = 2^{64}Y \bmod n$
$XY$	$C = 2^{64}XY\bmod n$ $\;\;\;\:= 2^{-64}AB\bmod n$

さて、この $2^{-64}AB\bmod n$ はどのように求めたらよいでしょうか？ここで $\bmod n$ が発生していたら意味がありません。これを事前計算された $R$ を使って回避するのが、モンゴメリ乗算です。 $Z \equiv 2^{-64}D \mod n$ となる $Z$ は次のように計算できます。

$Z = \frac{D}{2^{64}} - \frac{(DR \bmod 2^{64})n}{2^{64}}$

この手順には $2^{64}$ での除算や $2^{64}$ での剰余を求める部分がありますが、64bitCPUであれば実質ノーコストである点に注意します。

なぜこれでいいのでしょうか？まず、 $Z$ が整数であることから確認します。

$D - (DR \bmod 2^{64})n \equiv D - DRn \equiv D - D\times1 \equiv 0 \mod 2^{64}$

なので、 $D - (DR \bmod 2^{64})n$ は $2^{64}$ の倍数であり、 $2^{64}$ で割っても整数です。また、 $D$ や $(DR \bmod 2^{64})n$ は $2^{64}$ の倍数とは限りませんが、その差が $2^{64}$ の倍数であるということは、下位64bitは同じになっています。よって、もし $\frac{D}{2^{64}}$ や $\frac{(DR \bmod 2^{64})n}{2^{64}}$ を切り捨て除算で計算しても、結果は同じになります。 $Z = \frac{D - (DR \bmod 2^{64})n}{2^{64}}$ のように計算すると $D - (DR \bmod 2^{64})n$ の下位64bitの計算が一見必要そう（下位64bitは除算で消えるにせよ、ボローが生じるかを計算しないといけなさそう）ですが、実はそれは不要、ということを意味しています。

つぎに、 $Z \equiv 2^{-64}D\mod n$ であることを確認します。

$Z = 2^{-64}D - 2^{-64}(DR \bmod 2^{64})n \equiv 2^{-64}D \mod n$

大丈夫そうです。

このように、 $2^{64}$ などの切りのいい数を掛けた表現で持つことで、多少の事前計算と引き換えに剰余演算を除去するのが、モンゴメリ乗算です。

なお、単に $2^{64}$ を掛けた表現で持っているだけなので、加減算は普通の剰余環でのものと同じようにして実行可能です。

ポラードのロー法では非常に多くの剰余演算が使われるため、高速化のためには事実上必須のテクニックです。以下では、ポラードのロー法といったとき、暗黙でこのテクニックが適用されているものとします。

素数判定

ポラードのロー法で素因数分解するためには、素数であるかを（素因数分解によらずに）判定する必要があります。いろいろな高速な手法が提案されていますが、素因数分解にかかる時間が圧倒的なため、素数判定の高速化に躍起になる必要はありません。今回はミラーラビン素数判定法を愚直に実装したものを使いました。ミラーラビン素数判定法は確率的判定法ですが、 $2^{64}$ までの整数の素数判定に対しては、特定の7つの底を使えば間違いなく判定できることが知られています。

Binary GCD

Binary GCD（Stein's algorithm）は、最大公約数を求めるアルゴリズムです。二進法を採用している計算機（今日のほぼすべての計算機）に適しており、高速に動作します。

最大公約数を求めるアルゴリズムとして有名なユークリッドの互除法は、剰余演算が含まれるため、あまり高速ではありません。アルゴリズムの進行とともに入力依存で法が変化するため、剰余演算を事前計算で高速化することも不可能です。 Binary GCDは除算を含まないため、高速に動作します。

以前の記事（最大公約数をもっと高速に求める番外編（その2） - よーる）では、C言語コードを工夫して書けばclangが最適なコードを吐くようにできることを紹介しました（インラインアセンブリに手を出さなくて済むということです）。当時はgccだと最適なコードとならない問題がありましたが、現在のgcc（ジャッジサーバーで使われているのはgcc 13.2）ではほぼ最適なコードが出るようになっているようです。

60bit整数同士の最大公約数を求める場合を実験してみると、このコードのループは平均で40回程度実行されるようです。 Binary GCD (GNU MP 6.1.2)によれば1bitあたり0.68回の反復が必要、と書かれているのでこれと合致します。ループ一周は5cycleで実行可能*2なため、最大公約数を求めるのにかかる時間は200cycle程度です。

ポラードのロー法で、最大公約数をまとめて計算するテクニック

ポラードのロー法を愚直に実装すると、「乱数を生成し、差を計算し、それと $n$ の最大公約数を求める」というのを繰り返すコードになると思います。しかし、これはかなり遅くなります。先ほど示したように、最大公約数を求めるのにかかる時間は200cycle程度です。これは乱数の生成にかかる時間が10cycleや20cycleであるのに対して10倍も大きいです。これではアルゴリズムの実行時間の大半が最大公約数を求めることに費やされてしまいます。

これには面白い解決策があります。 $n$ との間の最大公約数を求める数を、1つの差とするのではなく、 $M$ 個の差を全部掛け合わせたもの（の $n$ での剰余）とするのです。これにより、最大公約数を求めるのにかかる時間を $1/M$ にすることができます。ただし、全部掛け合わせると $\bmod n$ で $0$ になってしまう場合（≒複数の因数が $M$ 個の中で同時に発見された場合）の対処が追加で必要になります。

この手法は非常に大きな性能向上につながるため、高速化のためには事実上必須のテクニックです。以下では、ポラードのロー法といったとき、暗黙でこのテクニックが適用されているものとします。

本題：私の提出コードの解説

ボトルネックの解析

ポラードのロー法のボトルネックは、疑似乱数を順次作っていく部分にあります。 $x_i$ のモンゴメリ表現から $x_{i+1} = x_i^2 + C \bmod n$ のモンゴメリ表現を求めるには、モンゴメリ乗算が1回必要です。モンゴメリ乗算には、以下が必要です。

64bit×64bit=128bitの乗算3回
- フルの128bitが必要なのが1回
- 下位64bitだけ得られれば良いのが1回
- 上位64bitだけ得られれば良いのが1回

これらの乗算はいずれも、IntelのSandy Bridge以降のプロセッサでは3cycleで実行可能です。これらの乗算は直列に依存しているため、これだけで9cycleが必要です。さらに一回の加減算が必須で、あわせて10cycleのレイテンシが避けがたいです。

差を求めて掛け合わせる部分でもモンゴメリ乗算が1回必要なので、合計で一つの乱数生成当たり6乗算必要（ブレントの改良方式の場合）ということになります。乗算器はパイプライン化されているので、10cycleの間にスループットとして10回の乗算が可能であり、たった6回しか乗算を行わないのは非常にもったいないです。このように、疑似乱数を順次作っていく部分がクリティカルパスとなっており、並列に実行できる命令が十分にないため、CPUはかなりの時間遊んでいます（とはいっても平均IPCが2を超えるような命令レベル並列性が豊富なプログラムなのですが、最近のCPUの最大IPCは4以上あるので、それには届かないということです）。

ボトルネックの解決策

乗算器が余っていることに着目し、並列に二つの疑似乱数系列を生成していくようにすることで、高速化が期待できます。このようにすると、ループ中に12乗算含まれる形になります。 IntelのCPUは1cycleにスループットとして1回の乗算が可能であり、このループの逆数スループットは12cycleとなります。これは疑似乱数生成のレイテンシは10cycleよりも長いので、理想的にはスループットに支配されたプログラムとなります。ただし、逆数スループット12cycleはクリティカルパスレイテンシ10cycleとあまり差がないため、あまりにも雑な命令スケジューリングが行われるとレイテンシが見えるようになってしまう点に注意する必要があります。普通のCPUの命令スケジューリングでは古い命令が優先的に選択されることを考えると、掛け合わせる部分の命令列よりも十分前に疑似乱数を生成する部分の命令列があれば良さそうです。この目的のためには、ソフトウェアパイプライン化を施せばよいです。

実は、以下の条件を満たす場合はソフトウェアパイプライン化を施さなくてもそこそこの命令列が生成されていました。

clangでコンパイルする
掛け合わせる部分のコードが実行されるのが、乱数生成直後ではない（例えば、差を計算するのに時間がかかるなど）

gccでコンパイルすると遅いのはなんでだろうとずっと悩んでいたのですが、clangを使っていても差を計算する部分を軽量化（後述）したときにかえって遅くなったことから、命令スケジューリングに問題があることに気づくことができました。ソフトウェアパイプライン化を施したコードは、上記条件を満たすときのclangの出力コードよりも5%程度高速であり、これが単独fastestを取る詰めの一手となったようです。

（あまり効果のない）軽量化

モンゴメリリダクションしない

普通の感覚では、差を掛け合わせた数と $n$ の最大公約数を求めるとき、まずは差を掛け合わせた数のモンゴメリ表現qに対し元の表現Qを求め（この手順をモンゴメリリダクションといいます）、その後Qとnの最大公約数を求めると思います。しかし、qはQを $2^{64}$ 倍した数を $\bmod n$ で考えているだけなので、nの奇数の約数aについて、Qがaの倍数であればqもaの倍数です。よって、元の表現Qを求めなくても、qとnの最大公約数を求めれば目的を達成できます。

正規化しない引き算

普通の感覚では、モンゴメリ表現xとモンゴメリ表現zの差を求める際、x < z ? x - z + mod : x - zなどと求めるはずです。しかし、モンゴメリ表現は $[0, n)$ の範囲に正規化されている必要はないため、x < zかにかかわらず常にx - z + modとしても大丈夫です。また、xはかなり長い間固定なので、x + modを計算しておけば、差を求めるのはsub命令一つでできるようになります。

なお、これにより掛け合わせる部分のコードが実行されるのが乱数生成直後になってしまいかえって遅くなってしまったことが、命令スケジューリングに問題があることに気づけたきっかけでした。

謎の数のモンゴメリ表現を使う

普通の感覚では、 $x_{i+1} = x_i^2 + C \bmod n$ の計算にモンゴメリ乗算を使う場合、 $C$ のモンゴメリ表現を求める必要があるように思われます。しかし、 $C$ の具体的な値は重要ではなく、 $C$ には任意の整数が使えます。よって、「モンゴメリ表現が1になるような謎の定数 $C$ を使った疑似乱数生成器」を使っているんだと考えることで、モンゴメリ表現に変換するコードを記述しなくて済むようになります。この定数 $C$ は $n$ に依存するためアルゴリズムが正しく動作するのか怪しいところがありますが、少なくとも $2^{40}$ までの合成数の非自明な因数を見つけることは可能なようです（Intel Core i9-12900Kにて24スレッド並列で1日くらいかけて確認）。

この技術は、あらゆるところに使われています。乱数のシードZ1, Z2はモンゴメリ表現である必要がありますが、いかにも普通の表現のような顔をして1などと書かれています。掛け合わせていくときの初期値q1, q2は、普通は1のモンゴメリ表現を使いますが、ここでもわざと普通の表現と混同して1などと書いています。

これらは、プログラム中の生存変数を減らすことを目的に行っています。ただ、頑張って減らしても結局gccが発生させるスピルを取り除くことはできませんでした（clangだと生存変数が多くてもスピルも無駄なmov命令もほぼなかったのに……）。最初はgccの生成するコードが遅い原因をスピルなのではないかと疑っていた（命令スケジューリングに問題があることに気づいていなかった）ため、生存変数を減らす努力をしていたのですが、結構無意味な努力だったようです。

疑似乱数生成器を固定する

一般に一つの疑似乱数生成器だけでは素因数分解ができない合成数があるため、疑似乱数生成器を[n]( uint64_t x ) { return __uint128_t(x) * x % n + 1; }などとハードコーディングするのは正しくない（無限ループしうる）ことを以前に指摘しました（間違ったポラードのロー法の使い方 - よーる）。しかし、今回は二つの疑似乱数生成器を並列で使っているため、このような問題は多分起こりません（確証なし）。よって、プログラムの命令数を減らすために、疑似乱数生成器をハードコーディングすることができそうです。プログラムの命令数を減らすことが高速化につながっているのかはよくわかっていません。

ポラードの方法（フロイドの循環検出法）でもブレントの方法でもない循環検出法

そんな大層なものではなく、ブレントの方法の亜種です。

周期を $P$ として、循環を検出する（ $n$ との最大公約数が1でない「差」を見てくる）のに必要な乱数生成回数は、フロイドの循環検出法では $3P$ 回程度、ブレントの改良方式では $2P$ 回程度です（正確には周期だけでなく循環に達するまでの反復回数（ρのしっぽの部分の長さ）を考えなければなりません。これに関するわかりやすい説明が前掲のポラード・ロー素因数分解法について #競技プログラミング - Qiitaに載っています）。このことからブレントの改良方式が優れているとの考えが主流ですが、これは必ずしも正しくありません。二つの手法における計算量を確認してみます。

フロイドの循環検出法では、一周当たり20cycle（乱数生成を直列に2回やるレイテンシ。もう一系列の乱数生成のレイテンシは並列実行で隠蔽可能）で乗算が12回（3乗算×乱数生成3回＋3乗算×掛け合わせ1回）です。循環検出までにかかる周回数は $P$ 回程度です。

ブレントの改良方式では、一周当たり10cycle（乱数生成を1回やるレイテンシ）で乗算が6回（3乗算×乱数生成1回＋3乗算×掛け合わせ1回）です。巡回検出までにかかる周回数は $2P$ 回程度です。

このようにしてみると、どちらが特に優れているなどとは言えなさそうに見えてくると思います。

ブレントの改良はもう一つあり、差を取る乱数の片側（提出コードにおけるx側）を変更した直後は掛け合わせを省くという最適化です。前掲のポラード・ロー素因数分解法について #競技プログラミング - Qiitaによれば、

もし r/2 < k < 3r/4 の間に約数が検出されるとすると、残りの 3r/4 ≦ k < r の間でも再度同じ約数が検出されるため

という説明がされています。これはアルゴリズムの正当性（および反復回数を上から抑える議論）としては正しいものですが、近年のCPUにおいて高速化することの説明としては正しくなさそうです。ブレントの改良方式では、一周当たり10cycleで乗算が6回なので、ここから掛け合わせを省いて乗算を3回に減らしても特に高速化しないからです（クリティカルパス上に無い命令が減ってもレイテンシに支配されているプログラムは高速化しません）。逆にこの“スキップ”により、せっかく乱数を生成しているにもかかわらず因数を見つけるチャンスを失っていることになります。この機会損失と最大公約数計算の手間のどちらが大きいかは簡単には分かりませんが、単純に「仕事量が減っているから高速化するはず」という雑な議論では済まなさそうです。

私の導入した二つの疑似乱数系列の並列生成を導入した場合、以下のようになります。

掛け合わせを省かない場合、一周当たり12乗算なのでスループットで支配されて12cycleかかる。よって一周は12cycle
掛け合わせを省く場合、一周当たり6乗算なのでスループットは十分で、レイテンシの10cycleがみえる。よって一周は10cycle

つまり、掛け合わせを省くことで多少は高速化するようです。しかし、実験してみるとこのような先ほどよりも良い条件にも関わらず、素因数分解の速度の向上は見られませんでした。やはり、“スキップ”中の因子発見チャンスを捨てている損失がある程度存在するようです。

フロイドの循環検出法を用いたものも実験してみましたが、ほぼ同じかやや遅いくらいの性能になりました。原因は正確には調査できていないですが、どうも「全部掛け合わせると $\bmod n$ で $0$ になってしまう場合」という運の悪いパターンの発生率が上がっているようでした。ソフトウェアパイプライン化がうまくいっていないなどの原因も考えられます。

以上の実験により、提出コードは“スキップ”を省いたブレントの改良方式を使っています。

パラメータの調整

差を512個掛け合わせてから最大公約数を求める

差を掛け合わせる個数 $M$ を大きくすることには、以下のような得失があります。

最大公約数を求める回数が減るため、高速化する
「全部掛け合わせると $\bmod n$ で $0$ になってしまう場合」の発生率が上がる（発生するとかなり低速になる）

実験を行った結果、 $M=512$ 程度を選ぶと実行時間が最小化されるという結論に達しました。

最悪実行時間に最適化

ブレントの改良方式では、差を取る乱数のうち固定側（提出コードにおけるx側）を更新するのは乱数生成の1, 2, 4, 8, ...回目にするのが普通です。しかしこれは増やしていくという部分が重要なのであって、最初は1回で更新、などとする必要はありません。逆にこのようにすると、頻繁に最大公約数を求める手順が挟まるため、大きな合成数の素因数分解では不利そうです。よって、提出コードでは、乱数生成の512, 1024, 2048, 4096, ...回目に更新することにしました。これにより、小さい合成数に対する速度がかなり低下しますが、（少なくともこの問題に対しては）小さい合成数に対する速度が上がっても得にならないので、このような最適化を行いました。そもそも、小さい合成数の素因数分解を速く行いたいのなら、ロー法ではなく試し割り法の方が向いているという事実もあります。

その他

決定的素数判定法のコードに関する疑惑

本題と関係ありませんが、ミラーラビン素数判定法でif( n <= a ) { return true; }とやるのはまずいという話を見つけました（ミラーラビン素数判定法についてメモのコメント）。これは一般には多分そうなのですが、提出コードで使われている64bit整数に対する決定的素数判定法（競技プログラミングでよく使われる方法）では問題ないはずです。まず、n == 2の場合は素数、そうでなくn % 2 == 0の場合は素数ないと判定する部分は問題ありません。その後、まずはa = 2でミラーラビン素数判定が行われます。底2での強擬素数（ミラーラビン素数判定が誤判定する合成数）は最小でも2047なので、このテストが通っている時点でif( n <= 7 ) { return true; }やif( n <= 61 ) { return true; }は正当です。七つの底を使う方では、そもそもそちらが実行される条件がn >= 4'759'123'141なので、if( n <= 1'795'265'022 ) { return true; }などの条件が成立することは決してありません。

プログラムの実行時間の比率と乗算器利用効率

おおよそ以下のような感じでした。

素数判定1.5%
最大公約数を求める部分4.3%
最大公約数を求め終わった時の分岐予測ミスペナルティ0.4%
乱数生成部分93.1%
その他0.7%

乱数生成部分は一周あたり12cycleで512回繰り返されるので6144cycle、最大公約数の計算は200cycle程度なので、比率のオーダー感としてはあっています。分岐予測ミスペナルティは20cycleくらいなので、最大公約数の計算にかかる時間との比率もあっています。乱数生成部分で乗算器を使った回数をカウントし、それを実行にかかった時間で割って得た乗算器利用効率は、86.3%でした。乱数生成部分が占める時間は93.1%なので、92.7%ほどは乗算器を有効利用できているということになりそうです。

なぜ二並列なのに二倍速くならないのか

まず、二並列にした時点でレイテンシ支配ではなくスループット支配に変わっています。レイテンシは一周当たり10cycle、逆数スループットは12cycle（つまり一系列当たり6cycle）なので、この時点で1.67倍までしか速くならないことになります。

また、ロー法で循環を検出するまでにかかる平均周回数が半分になるわけではないのも原因のようです。ロー法で循環を検出するまでにかかる周回数が指数分布に従っているとすれば、（仕事を並列化するのではなく）二系列同時にやって早く巡回を検出したほうを採用、という手順で平均周回数が半分になって、二倍高速化するはずです。しかし、ロー法で巡回を検出するまでにかかる周回数は指数分布から外れているらしく、二系列同時にやっても平均周回数は半分にはならず、2/3程度にしかならないようです。分布の標準偏差は平均値の半分程度（指数分布なら平均値と同じくらいのはず）であり、短い巡回が見つかる可能性がかなり低いことを示しています。非巡回部分（ρのしっぽの部分）に原因があるのかもしれません。また、歪度は1程度です（これは短い巡回には1という自明な下限があるにもかかわらず長い巡回には上限がないことが原因で、普通です。指数分布なら2あるはずなので、やはり短い巡回が少ないことを意味しているようです）。尖度は1（四次中央モーメントで書くと4）程度であり、裾の重い分布であることを示しています。

これらから計算すると、提出コードは二並列にしているにもかかわらず実行時間は $2/3\times12/10 = 0.8$ 倍程度、あるいは実行速度は1.25倍程度、にしかならないとわかります。