Hatena::ブログ(Diary)

あどけない話

2010-01-31

素数判定

要約:素数判定に使われるミラーラビン法を解説しながら、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 5632
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 となれば、フェルマーテストを騙すことができる。

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 5631
expmod 562 2 5631

たとえば、67 は二乗しても 1 にはならない。

expmod 67 2 563548

ここで、カーマイケル数 561 に対しても同じ計算をしてみる。

expmod 1 2 5611
expmod 560 2 5611
expmod 67 2 5611

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]

追記

カーマイケル数は、フェルマーテストを必ず騙すことを思い出したので、書き換えた。そう言えば、以前自分で証明したことがあったのだ。^^;

ミラーラビン法は、対偶をとる必要はないと気付いたので、コードを書き換えた。

乗越雅光乗越雅光 2012/10/04 06:50 教えてください。

「素数判定ーあどけない話」のMiller-Rabin法の記事は大変興味深く読ませていただきました。ありがとうございます。
http://d.hatena.ne.jp/kazu-yamamoto/20100131/1264921432の後半にある「同様の考察を続けると、n が素数であるための条件は以下の通り。

a^d ≡ 1 (mod n)、または
ある r (0 ≦ r < s - 1) に対して a^(d × 2^r) ≡ -1 (mod n) 」
で2つの合同式の値はそれぞれ1又はー1となっていますが、双方とも±1でもよいのではないでしょうか?
どちらにしても少なくとも1回は2乗されるのですから。
初学者の質問ですので的外れ化もしれませんが、よろしくご教授根が合います。
以上

kazu-yamamotokazu-yamamoto 2012/10/05 15:54 よいのではないかと言われれば、いいんですけど、条件としては冗長であるということです。たとえばプログラムを書くときに、-1 であるか判定すればいいところを、1 かまたは -1 と判定するのは無駄ですね。

スパム対策のためのダミーです。もし見えても何も入力しないでください
ゲスト


画像認証

トラックバック - http://d.hatena.ne.jp/kazu-yamamoto/20100131/1264921432