要約：素数判定に使われるミラーラビン法を解説しながら、Haskell で実装してみる。
フェルマーテスト
大きな数を確実に素数だと判定するには、大変時間がかかるので、実用的には「ほぼ素数だ」と確率的に判定する。確率的な素数判定の代表格がフェルマーテストである。
フェルマーテストには、以下に示すフェルマーの小定理を利用する。
a^p ≡ a (mod p)
a は任意の整数。p は素数である。法 p の下で a を p 乗したものは、a と合同であると言う意味だ。a には制限はないが、特に a を p より小さい整数、0 ≦ a ≦ p - 1 とすれば、a を p 乗して、p で割ると a に戻るとも解釈できる。
最初に見たときは、だからどうしたと思われるかもしれない。しかし、有名なフェルマーの大定理が実用上何の役にも立たないのに対し、フェルマーの小定理はいろんな場面で活躍する。
実際に計算してみよう。準備として、base を ex 乗して m で割った余りをとる関数 expmod を以下のように定義する。この関数は、高速で、しかも使用するメモリーが少ない。
expmod :: Integer -> Integer -> Integer -> Integer expmod base ex m | ex == 0 = 1 | even ex = square (expmod base (ex `div` 2) m) `mod` m | otherwise = (base * (expmod base (ex - 1) m)) `mod` m where square x = x * x
素数の例として 563 を取り上げる。
expmod 2 563 563 → 2 map (\x -> expmod x 563 563) [0..562] == [0..562] → True
すべての素数はフェルマーの小定理を満たす。よって大きな整数 n が素数であるかを判定するには、任意の整数 a に対して n がフェルマーの小定理を満たすか調べればよい。
フェルマーの小定理を使って、n が素数か確率的に判定する fermat という関数は以下のように書ける。
fermat :: Integer -> Integer -> Bool fermat n a = expmod a n n == a
素数に対して実験してみよう。
fermat 563 2 → True fermat 563 3 → True fermat 563 4 → True
素数であると判定された。
次に、合成数 560 について実験してみる。
fermat 560 2 → False fermat 560 3 → False fermat 560 4 → False
合成であると判定された。
乱数を自動的に生成して、テストするようにしてみよう。
import Random import Control.Applicative fermatRandom :: Integer -> IO Bool fermatRandom n = fermat n <$> getRandom 2 (n - 2)
実際に使ってみる。
fermatRandom 563 → True fermatRandom 560 → False
Bool が返ってきているように見えるが、乱数を使っているので IO Bool が返っていることに注意。
疑素数
フェルマーテストを騙す数を疑素数という。たとえば、合成数 341 (11 × 31) は、乱数がたまたま 2 となれば、フェルマーテストを騙すことができる。
>||haskell|
fermat 341 2
→ True
|
つまり、こういうことだ。
1回のテストでは心もとないので、通常は k 回ほどこのテストを繰り返して確率を上げる。
import Control.Monad fermatTest :: Int -> Integer -> IO Bool fermatTest times n = and <$> replicateM times (fermatRandom n)
341 に対して、10 回テストを繰り返してみよう。
fermatTest 10 341 → False
カーマイケル数
疑素数 341 が 2 のときにしかフェイルマーテストを騙せなかったのに対し、疑素数はフェルマーテストを騙す。
最小のカーマイケル数は、561 (3 × 11 × 17) である。実際に試してみよう。
fermat 561 2 → True fermat 561 3 → True fermat 561 4 → True length $ filter (fermat 561) [0..560] → 561
恐ろしいことにすべての場合で、騙していると分かる。
フェルマーの小定理の変形
カーマイケル数に騙されにくくする方法の1つは、変形したフェルマーの小定理を用いることである。
フェルマーの小定理は、以下の通りであった。
a^p ≡ a (mod p)
a が p の倍数でなければ、両辺を a で割ってよいので、以下のように変形できる。(法 p の下では、p の倍数は 0 と同じ意味になるので、これで割ってはいけない。)
a^(p - 1)≡ 1 (mod p)
我々にとっては、0 を取り除いた、1 ≦ a ≦ p - 1 の範囲で考えれば十分である。変形フェルマーの小定理を使ったテストは以下のように実装できる。
fermat' :: Integer -> Integer -> Bool fermat' n a = expmod a (n - 1) n == 1
実際に使ってみよう。
fermat' 561 2 → True fermat' 561 3 → False fermat' 561 4 → True
カーマイケル数 n が変形フェルマーテストを騙す a の条件は、gcd(n,a) = 1 である。
たとえば、561 の場合、3、11、17 を因数に持つ a は、変形フェルマーテストを騙せない。
length $ filter (fermat' 561) [1..561] → 320
フェルマーテストを変形すると、騙されにくくなったことが分かる。
ミラーラビン法のアイディア
変形フェルマーテストの確率をさらに上げるために、自明でない 1 の平方根を探す方法がある。これは、ミラーラビン法と呼ばれている。
まず、1 の自明な平方根とはなんだろうか？ その数 x は二乗したら 1 になるはずである。
x^2 = 1 (mod n)
変形すると
(x - 1)(x + 1) = 0 (mod n)
n が素数であれば、n は (x - 1) も (x + 1) も割り切らない。そこで、x は 1 か -1 であると分かる。これを 1 の自明な平方根という。-1 とは、n - 1 のことである。
合成数の場合は、1 や -1 以外に、二乗して 1 となる整数が存在する。その数を 1 の自明でない平方根と呼ぶ。
意味が分からないと思うので、例を挙げよう。
素数 563 に対して、自明な平方根は、1 と 562 である。
expmod 1 2 563 → 1 expmod 562 2 563 → 1
たとえば、67 は二乗しても 1 にはならない。
expmod 67 2 563 → 548
ここで、カーマイケル数 561 に対しても同じ計算をしてみる。
expmod 1 2 561 → 1 expmod 560 2 561 → 1 expmod 67 2 561 → 1
67 を 2 乗すると 1 になった。これが1の自明でない平方根である。つまり、66 と 68 を素因数分解してみれば、3、11、17 という因数が現れる！
僕が計算したところ、561 に対して1の自明でない平方根は 67,188,254,307,373,494 であった。
ミラーラビン法のアルゴリズム
ミラーラビン法は、変形フェルマーテストを計算する過程で、1の自明でない平方根を見つけるアルゴリズムである。
素数か合成数か判定したい整数 n は奇数である。(n = 2 は自明すぎて対象外。それ以外の偶数は、もちろん合成数である。) よって、n - 1 は偶数であり、因数 2 を持つ。
n - 1 を d × 2^s と表すことにする。d は奇数である。変形フェルマーの小定理より。
a^(d × 2^s) ≡ 1 (mod n)
つまり、n が素数であれば、a を d × 2^s 乗したときは 1。素数は、1 の自明な平方根しか持たないので、a を d × 2^(s-1) 乗したときは、1 か -1 である。
a を d × 2^(s-1) 乗したとき 1 となった場合、a を d × 2^(s-2) 乗したときは、1 か -1 となる。
同様の考察を続けると、n が素数であるための条件は以下の通り。
- a^d ≡ 1 (mod n)、または
- ある r (0 ≦ r ＜ s - 1) に対して a^(d × 2^r) ≡ -1 (mod n)
これを実装すると、以下のようになる。
factor2 :: Integer -> (Int, Integer) factor2 n = factor2' (0,n) where factor2' (s,d) | even d = factor2' (s + 1, d `div` 2) | otherwise = (s,d) millerRabin :: Integer -> Integer -> Bool millerRabin n a = let (s,d) = factor2 (n - 1) aExpD = expmod a d n double x = expmod x 2 n evens = take s $ iterate double aExpD in aExpD == 1 || any (== n - 1) evens millerRabinRandom :: Integer -> IO Bool millerRabinRandom n = millerRabin n <$> getRandom 2 (n - 2) millerRabinTest :: Int -> Integer -> IO Bool millerRabinTest times n = and <$> replicateM times (millerRabinRandom n)
実際に実験してみよう。
length $ filter (millerRabin 561) [1..560] → 10
fermat' では 320 通り騙されていたのに、millerRabin では 10 通りしか騙されない。つまり、310個の場合の計算過程において、1の自明でない平方根(67,188,254,307,373,494)を見つけたことになる。
もちろん素数は、すべての場合に millerRabin テストにパスする。
length $ filter (millerRabin 563) [1..562] → 562
最終実験
carmichaels :: [Integer] carmichaels = [561,1105,1729,2465,2821,6601,8911,10585,15841] primes :: [Integer] primes = [563,1009,1013,1019,10007,10009,10037,100003,100019]
カーマイケル数は、高い確率で fermatTest を騙す。
mapM (fermatTest 10) carmichaels → [True,True,True,True,True,True,True,True,True] mapM (fermatTest 10) primes → [True,True,True,True,True,True,True,True,True]
一方、カーマイケル数は、millerRabinTest を騙すこともあるが、その確率は高くない。
mapM (millerRabinTest 10) carmichaels → [False,False,False,False,False,False,False,False,False] mapM (millerRabinTest 1) carmichaels → [True,False,False,False,False,False,False,False,False] mapM (millerRabinTest 10) primes → [True,True,True,True,True,True,True,True,True]
