64bit符号なし整数を素因数分解するため、事前に32bit符号なし整数で表せる素数をリストアップしておくとします。
方法1: 素数のリストをビットベクトルで保持し、試し割り
32bit符号なし整数で表せる奇数が素数であるかを、0, 1で表した列(ビットベクトル)を事前に作成し、それを素因数分解プログラムのバイナリに埋め込みます。 ビットベクトルのエントリ数は2Gi個で、1Byte=8bitであるとして256MiByteのバイナリが完成します。
32bit符号なし整数で表せる素数は203280221個(203M個)存在するので、この方法を使うと最悪で203M回の除算が発生します。
最近のCPUであっても除算器はパイプライン化されていないので、203M回×除算器の占有サイクル数くらいの時間がかかります。
この方法を使ったとき、手元のマシンでは、4.8秒くらいかかりました。ちょっと時間がかかりすぎな気がします。
方法2: 素数の逆元をファイルから読み込み、判定
m % p == 0
という剰余演算はp
が奇数の定数であれば、p
の264における逆元を用いて最適化可能です。これを事前に計算しておく方式です。一つの素数当たり16Byteの情報が必要なため、バイナリに埋め込むことはできませんでした。
事前に計算した定数をファイル(バイナリファイルで3.03GiB、おそらくメモリ上にキャッシュ済み)から読み込んだ場合、素因数分解には5.2秒かかりました。一つの素数で割れるかを確認するのは1cycleで可能で非常に高速なのですが、読み込みに5.0秒くらいかかっているので、この方式をこれ以上改善することはできなさそうです。
mmap
で読み込んだら、読み込み1秒、計算1秒くらいまで高速化しました。
素因数が小さいとき、読み込みが高速になるのもよい点です。
さすがに事前計算に時間がかかっているだけあって高速な手法です。
方法3: ポラードのρ法
以前に実装した高速なgcd実装を用い、フロイドの循環検出法を用いたオリジナルバージョンを実装しました。
素数であるかわかってない入力が来る場合、ポラードのρ法は素数であるかを確認するのに非常に不向きなので、別途素数判定アルゴリズムを併用することになります。
素数でないとわかっている入力に対し、何らかの因数を見つけるのにかかる時間は、(乱択アルゴリズムなのでぶれが大きいですが)平均5秒、最悪20秒程度でした。
素因数は232ほどの大きさであり、これの平方根(65536)くらいの回数の乱数生成で発見できると見積もるようですが、実際には数千万個の乱数を生成しないと素因数が発見できない例が多々ありました。
事前の計算は不要でプログラムサイズも小さいですが、素数のリストを事前に作成した試し割り法よりは低速だということがわかりました。