**ポラード・ロー法(Pollard's rho algorithm)**は非自明な約数を高速に発見するアルゴリズムです。
ただし、このアルゴリズムは必ず約数を発見できるわけではなく、失敗する可能性があります。
今回は、このアルゴリズムをHaskellで実装して、その速度を計測してみました。
また、素朴な試し割り法による約数の探索についても速度を計測して、ポラード・ロー法の速さと比較してみようと思います。
実装
次のように実装しました。
-- ポラード・ロー法
rho :: Integer -> Maybe Integer
rho n
| n < 4 = Nothing
| otherwise = loop 2 2 1 -- フロイドの循環検出法で循環が検出されるまで試し割りする
where
random = \x -> mod (x*x + 1) n --疑似乱数
loop x0 y0 d
| d == n = Nothing -- 非自明な約数がない、もしくは発見できなかった
| d /= 1 = Just d -- d == gcd (x0-y0) n なので d は n の約数
| otherwise = loop x1 y1 (gcd (x1-y1) n)
where
x1 = random x0
y1 = random (random y0)
疑似乱数を用いて約数を探すのですが、約数を発見する前に疑似乱数が循環してしまう可能性があります。
そこで、ポラード・ロー法ではフロイドの循環検出法を用いて疑似乱数の循環を検出した場合は探索を終了します。
なので、ポラード・ロー法では非自明な約数を発見する前に探索が打ち切られる可能性があります。
速度の比較
次にポラード・ロー法の速度を測ってみます。
まず、GCHiオプション+sで経過時間などを表示するよう設定します。
*Main> :set +s
ここでは$100000980001501(=10000019\times10000079)$の約数を探索することにします。
*Main> rho 100000980001501
Just 10000079
(0.02 secs, 10,298,096 bytes)
ちゃんと約数($10000079$)を探すことができました。
経過時間は$0.02$秒と表示されています。
次に、素朴な試し割り法での速度を計測してみます。
試し割りのプログラムは以下のように実装しました。
-- trial division
trial_div :: Integer -> Maybe Integer
trial_div n = head' (filter ((0 ==) . mod n) (takeWhile (\d -> d*d <= n) (2:[3,5..])))
where
head' xs
| xs == [] = Nothing
| otherwise = Just (head xs)
$2$と奇数からなるリスト2:[3,5..]のうち$\sqrt n$以下のものについて、$n$を割り切るもののみをフィルタして、その結果の先頭を返します。
つまり、$n$を割り切る非自明な約数のうち、最小のものを求めています。
$\sqrt n$以下で存在しない場合、$n$を割り切る非自明な約数はないのでNothingを返します。
ポラード・ロー法のときと同様に時間を計測します。
*Main> trial_div 100000980001501
Just 10000019
(1.86 secs, 1,160,342,328 bytes)
約数($10000019$)を見つけることができましたが、$1.86$秒かかりました。
ポラード・ロー法の方が数十倍速い結果になりました。
非自明な約数を発見できない例
先程も述べたように、ポラード・ロー法では非自明な約数を発見できない場合があります。
どの程度発見できないのか、確認してみます。
以下は、今回の方法で非自明な約数を見つけることに失敗した合成数の例です。
$$4,8,16,21,22,25,32,44,62,64,88,95,\dots$$
このような合成数の数はどの程度あるのでしょうか？
$n$以下の合成数のうち、ポラード・ロー法で非自明な約数を発見できた数を計算してみます。
以下のように、合成数を数えるプログラムと、ポラード・ロー法の成功率を計算するプログラムを用意しました。
count_composite :: Integer -> (Integer -> Maybe Integer) -> Int
count_composite n algorithm = length (filter ((Nothing /=) . algorithm) [0..n])
rate_success :: Integer -> Double
rate_success n = (f rho) / (f trial_div) where f = fromIntegral . (count_composite n)
結果は以下のとおりです。
*Main> rate_success 10
0.6
*Main> rate_success 100
0.8378378378378378
*Main> rate_success 1000
0.9626955475330926
*Main> rate_success 10000
0.9905359179019384
*Main> rate_success 100000
0.9973895826650592
$100000$以下の合成数の$99%$以上について、非自明な約数を発見することができています。
まとめ
- ポラード・ロー法は素朴な試し割り法より速い。
- ポラード・ロー法は失敗する可能性がある。