Haskellで素数を探そうと思った話

はじめに

SICPを読んでたらフェルマーテストについて書いてあって、ほんとにこれ速いのかやってみようと思った話です。

ソースはfermat_test.hsからダウンロード!

フェルマーテストとは

フェルマーの小定理はご存知ですよね？

$p$を素数のとき、$a$を$p$と互いに素な整数に対して、 $$a ^ p \equiv a \pmod{p}$$ が成り立つ。

これの対偶を考えます。

整数$n$について、$n$と互いに素な整数$a$に対して、 $$ a ^ n \not \equiv a \pmod{n} $$ が成り立つとき、整数$n$は合成数である。

これによって、ある数が素数であることは示せないが合成数であることを速く示すことはできます。

実装

部品を 1 つ 1 つ揃えていきます。

互いに素

最小公倍数を求める関数を考えます。

gcd' :: Integer -> Integer -> Integer
gcd' n 0 = n
gcd' n m = gcd' m (n `mod` m)

これを用いて、

gcd' n m == 1

とすれば$n, m$は互いに素です。 (Haskell には標準で gcd があるのでそれを使っても良い)

$a^{n-1} \pmod{n}$を求める

いわゆる、$a^7 = a * (a^3)^2 = a * (a * a^2)^2$として計算量をへらすやつを使います。途中で$\pmod{n}$を挟みます。

exp_mod :: Integer -> Integer -> Integer
exp_mod a n = iter a (n-1) n
  where
    iter a n base
      | n == 1 = a
      | even n = (iter a (n `div` 2) base) ^ 2 `mod` base
      | otherwise = a * (iter a ((n-1) `div` 2) base) ^ 2 `mod` base

これを用いて、

exp_mod a n

を計算することで、$a^{n-1} \pmod{n}$が計算できる。

素数判定

素数判定を実装すると次のようになります。

is_prime' :: Integer -> Bool
is_prime' n = and $ map (\a -> test n a) [2..(n-1)]
  where
    test n a
      | gcd' n a == 1 = (1 == exp_mod a n)
      | otherwise = False

最小公倍数が 1 でない整数が見つかったとき、すぐに素数でないことがわかります。
互いに素な整数$a$に対して$a^{n-1} \pmod{n}$を求めて$1$と一致するか調べます。一致しなければすぐに素数でないことがわかります。
Haskell は遅延評価なのでandはFalseが見つかった瞬間にFalseを返すように出来てるはずです。

この実装では$n$未満のすべての自然数に対して、フェルマーテストを実行しています。この方法では素数$p$に対しては、$2 \sim (p-1)$で割っているので完全に素数を判定することが出来ています。 [2..(n-1)]の部分をランダムな自然数に変更することで、フェルマーテストになります。(Haskell の乱数めんどくさいのでやりません。)

まとめると

これまでのコードをまとめると

gcd' :: Integer -> Integer -> Integer
gcd' n 0 = n
gcd' n m = gcd' m (n `mod` m)

exp_mod :: Integer -> Integer -> Integer
exp_mod a n = iter a (n-1) n
  where
    iter a n base
      | n == 1 = a
      | even n = (iter a (n `div` 2) base) ^ 2 `mod` base
      | otherwise = a * (iter a ((n-1) `div` 2) base) ^ 2 `mod` base

is_prime' :: Integer -> Bool
is_prime' n = and $ map (\a -> test n a) [2..(n-1)]
  where
    test n a
      | gcd' n a == 1 = (1 == exp_mod a n)
      | otherwise = False

となる。

遊ぶ

Haskell の素数ライブラリを install する。

cabal install primes

$2^{67}-1$を判定

さて、$ 2^{67} - 1 = 147573952589676412927$は素数ではありません。 (参考:メルセンヌ数)

Haskell のライブラリで素数か確認すると

ghci> import Date.Numbers.Primes
ghci> isPrime $ 2^67 - 1

処理が終わりません。($\sqrt{n}$までprimesで割っていく実装とここに書いてある)

しかし、今回つくったis_prime'では

ghci> is_prime' $ 2^67 - 1
False

一瞬で処理が終わり合成数であることがわかります。(すごい)

このアルゴリズムのすごいところは、素因数分解は出来ないけど合成数であることはすぐに分かるところですね。判定したい数$n$の

$$ 2^{n-1}, 3^{n-1} \cdots \pmod{n} $$

を計算するしていって 1 にならないのを一つ見つけるだけで良い。

メルセンヌ数探し

primes' :: [Integer]
primes' = 2: filter is_prime' [3,5..]

is_prime'' :: Integer -> Bool
is_prime'' n = and $ map (\a -> test n a) $ takeWhile (< n) [2,3,5]
  where
    test n a
      | gcd' n a == 1 = (1 == exp_mod a n)
      | otherwise = False

を追加します。この関数は2,3,5の倍数でない数に対して、

$$ 2^{n-1}, 3^{n-1}, 5^{n-1} \pmod{n} $$

を計算します。 (これによってis_prime''でTrueになったとしても素数でない数が存在することになります。)

メルセンヌ数とは、

p を素数とし$M_p = 2^p-1$もまた素数になる場合、素数$p$をメルセンヌ数といいます。

早速探してみましょう。(1 分程度動かしました。)

ghci> filter (\p -> is_prime'' $ 2^p -1) primes'
[2,3,5,7,13,17,19,31,61,89,107,127,521,607,1279,2203,2281,3217,4253,4423

実際のメルセンヌ数は次の通りです。(参考:メルセンヌ数)

p = 2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689, ..

今回探した範囲では完全に一致してますね。(すごい)

ちなみに、Haskell の素数ライブラリで 1 分程度探すと

ghci> filter (\p -> isPrime $ 2^p -1) primes
[2,3,5,7,13,17,19,31

となり、31 までしか探せません。

メルセンヌ数はリュカ・テストで速くきちんとした素数判定ができるらしい。

カーマイケル数

少し気になるのがis_prime''を通過するけど、実際に素数でない数はどんなかずなのか？ですよね。実際に探してみましょう。

ghci> filter (not . isPrime) $ filter is_prime'' [2..]
[1729,2821,6601,8911,15841,29341,41041,46657,52633,63973,75361,101101

となり、10 万以下に 11 個あります。(結構少ない)

これらの数$n$は、$2, 3, 5$と互いに素であり、

$$ 2^{n-1}, 3^{n-1}, 5^{n-1} \equiv 1 \pmod{n} $$

を満たす合成数である。

じゃあ、合成数$n$について、n と互いに素な任意のｍに対して

$$ m^{n-1} \equiv 1 \pmod{n} $$

が成り立つような整数を探したくなります。 (このような合成数$n$をカーマイケル数といいます。)

このような数はフェルマーテストでは弾くことは出来なくて、実際に素因数を見つけて弾くしかありません。

カーマイケル数を実際に探してみましょう。

import Data.Numbers.Primes hiding (primes') -- プログラムの最初に追加

is_carmichael :: Integer -> Bool
is_carmichael n = and $ ((not . isPrime) n):(map (\a -> test n a) [2..(n-1)])
  where
    test n a
      | gcd' n a == 1 = (1 == exp_mod a n)
      | otherwise = True

not . isPrimeを追加して素数でない数から探します。
is_primeをotherwise = Trueとするだけで良い。(フェルマーテストのみを通過できる数を探せる。)

ghci> filter is_carmichael [2..]
[561,1105,1729,2465,2821,6601,8911,10585,15841,29341,41041,46657,52633,62745,63973,75361,101101

10 万以下に 16 個あった。

これらの数はたまたま約数でテストしない限りフェルマーテストで弾くことができない。(憎い)

並列化

Haskell の並列計算については、あんまりわかってないけど適当に実装してみた。

import Control.Parallel.Strategies

これ以降、コンパイルは

ghc -O2 fermat_test.hs -threaded -rtsopts

で行い、実行は

time ./fermat_test +RTS -N4

で行う。

-s をつけると詳細が見れる。
-N4 はスレッド数

実装

filter の部分で並列にする。

parFilter :: (a -> Bool) -> [a] -> Eval [a]
parFilter f [] = return []
parFilter f (a:as) = do
        b <- rpar (f a)
        bs <- parFilter f as
        return (if b then a:bs else bs)

main :: IO ()
main = do
  let ps = takeWhile (< 5000) primes'
      mersenne = runEval $ parFilter (\n -> is_prime'' (2^n - 1)) ps
  print mersenne

rparでタスクが作られるらしい。(タスクの振り分けは Haskell が勝手にやってくれるらしい。)

実行結果は

real    0m7.037s
user    0m27.018s
sys     0m0.313s

となった。

普通のfilterで実行すると

real    0m27.217s
user    0m28.602s
sys     0m3.278s

となり、だいたい 4 倍も速くなった。(すごい)

憎きカーマイケル数を探す

is_carmichaelでは調べたい合成数$n$について、互いに素な全ての整数$m$に対して、

$$ m^{n-1} \equiv 1 \pmod{n} $$

が成り立つものを言うが、もし$m = p_1 \cdot p_2 \cdots$と素因数分解できるとき、

$$ p_1^{n-1} \equiv 1, p_2^{n-1} \equiv 1, \cdots \pmod{n} $$

ならば、自動的に$m^{n-1} \equiv 1 \pmod{n}$が成り立つから、カーマイケル数を求めるには、

合成数$n$について、$n$未満の任意の素数$p$について、

$$ p^{n-1} \equiv 1 \pmod{n} $$

を計算すれば良いことになる。

is_carmichael' :: Integer -> Bool
is_carmichael' n = and $ map (\a -> test n a) $ takeWhile (< n) primes
  where
    test n a
      | gcd' n a == 1 = (1 == exp_mod a n)
      | otherwise = True

main :: IO ()
main = do
  let carmichael = runEval $ parFilter is_carmichael' $ filter (not . isPrime) $ takeWhile (<1000000) [2..]
  print carmichael

並列化を行う際にはrparで生成されるタスクが少ない方が良いので、parFilter に通す前に素数は弾いてしまっています。
上の考察から素数のみでテストするのはis_carmichael'で実装。(これだけで約 4 倍速くなった。)

実行結果は

real    0m5.806s
user    0m10.653s
sys     0m1.572s

となった。普通のfilterで実行すると

real    0m9.437s
user    0m9.912s
sys     0m1.567s

となり、だいたい 2 倍も速くなった。

(計算しかしてないからユーザーモードの実行時間とリアルの実行時間を見るだけで、並列化の効率を見ることができるかも？)

感想

素因数をみつけるという方法以外で合成数かどうか判断できるのがすごいを思った($M_{67}$を判定できたのはすごい)。大きな素数と大きな素数の積についてフェルマーテストを行い合成数であることがわかっても、その約数はぜんぜんわかりませんっていう不思議な状態になる。
1 回フェルマーテストで弾ける合成数の割合とかは気になる。また、どの数でフェルマーテストをすると弾きやすいのかとかが気になる。
並列化してパソコンに負荷をかけるのが楽しかった。(タスクが増えすぎてsegmentation faultしたりした。)
カーマイケル数を効率的に探す方法とか、フェルマーテストの亜種とか、ぜんぜん違う素数判定のアルゴリズムとかを調べたら面白そうだなとおもった。

つづく…

561,1105,1729,2465,2821,6601,8911,10585,15841,29341,41041,46657,52633,62745,63973,75361,101101,115921,126217,162401,172081,188461,252601,278545,294409,314821,334153,340561,399001,410041,449065,488881,512461,530881,552721,656601,658801,670033,748657,825265,838201,852841,997633,1024651,1033669,1050985,1082809,1152271,1193221,1461241,1569457,1615681,1773289,1857241,1909001,2100901,2113921,2433601,2455921,2508013,2531845,2628073,2704801,3057601,3146221,3224065,3581761,3664585,3828001,4335241,4463641,4767841,4903921,4909177,5031181,5049001,5148001,5310721,5444489,5481451,5632705,5968873,6049681,6054985,6189121,6313681,6733693,6840001,6868261,7207201,7519441,7995169,8134561,8341201,8355841,8719309,8719921,8830801,8927101,9439201,9494101,9582145,9585541,9613297,9890881

はじめに#

フェルマーテストとは#

実装#

互いに素#

$a^{n-1} \pmod{n}$を求める#

素数判定#

まとめると#

遊ぶ#

$2^{67}-1$を判定#

メルセンヌ数探し#

カーマイケル数#

並列化#

実装#

憎きカーマイケル数を探す#

感想#