ミラー・ラビン素数判定法

数学 プログラミング

ミラー・ラビン素数判定法（Miller-Rabin primarity test）はフェルマーテストを改良したものなので、まずはその説明をします。

フェルマーの小定理は、pを素数、aをpと互いに素な自然数とすると、

an-1 ≡ 1 (mod p)

が成り立ちます。例えば、p = 7, a = 2とすると、

26 = 64 ≡ 1 (mod 7)

ということは、自然数pに対してこれが成り立たなければpは合成数ということになります。
しかし、これが成り立てばpは素数ということにはなりません。例えば、a = 2でpが1000までで上が成り立つ合成数があるか調べてみましょう。

def sieve(max_n):
    a = [ True ] * (max_n + 1)
    for p in (n for n in xrange(2, max_n + 1) if a[n]):
        for k in xrange(p * 2, max_n + 1, p):
            a[k] = False
    return filter(lambda n: not a[n], xrange(2, max_n + 1))

def pow_mod(n, e, d):
    if e == 0:
        return 1
    else:
        m = pow_mod(n, e / 2, d)
        return m * m * n % d if e & 1 else m * m % d

def is_psuedo_prime(n, a):
    return pow_mod(a, n - 1, n) == 1

N = 1000
not_primes = sieve(N)
for n in (n for n in not_primes if is_psuedo_prime(n, 2)):
    print n,

341 561 645

この3つが出てきました。これを擬素数と呼びます。10000まで22個あります。
a = 2だけでなく他の底についてもテストを行うと合成数であると判定できることがあります。例えば、

3340 ≡ 56 (mod 341)

となって、341は合成数であることがわかります。
どんなpと素でないaに対しても剰余が1になる合成数があります。

def is_Carmichael(n):
    return all(is_psuedo_prime(n, a)
            for a in xrange(2, n) if gcd(n, a) == 1)

N = 10000
not_primes = sieve(N)
for n in ifilter(is_Carmichael, not_primes):
    print n,

561 1105 1729 2465 2821 6601 8911

これをカーマイケル数（Carmichael number）と言います。すなわち、いくらテストする回数を増やしても誤判定することがあるということです。


これを改良したのがミラー・ラビン法です。
例えばp = 41を考えます。フェルマーの小定理より、

a40 ≡ 1 (mod 41)

指数を2で割り切れるまで割っていきます。a20の剰余は1か-1になります。1とすると、a10の剰余は1か-1になります。1とすると、a5の剰余は1か-1になります。

指数 40  20  10   5
剰余  1┳ 1┳ 1┳ 1
       ┗-1┗-1┗-1

指数が5,10,20のうちどれかの剰余が-1、または5のとき剰余が1になります。これをチェックするとなぜか誤判定が少なくなります。底が2だけなら10000以下で、

2047 3277 4033 4681 8321

これだとたいしたことはないですが、例えば底が2、3でダブルチェックすると、10^6まで出てきません。チェックの回数が多くなればなるほど誤判定が少なくなるらしいです。通常は疑似乱数で底を決めてチェックするそうです。
チェックの回数を増やすと、非常に高速でかつ正確な素数判定を行うことができます。素数で割っていく方法では対象となる自然数が大きくなれば遅くなりますが、この判定法はあまり遅くならないのが特徴です。ですから、非常に大きな、具体的には億を超えるような数の素数性を判定するのに有効です。

from random import seed, randrange

def Miller_Rabin(n):
    if n % 2 == 0:
        return False
    
    def MR_core(d, a, s = 1):
        if s == 1:
            return MR_core(d >> 1, a, 2) == n
        elif s == 2:
            if d & 1:
                m = MR_core(d >> 1, a, 3)
                s = m * m % n * a % n
                return n if s == n - 1 or s == 1 else s
            else:
                m = MR_core(d >> 1, a, 2)
                if m == n:
                    return n
                else:
                    s = m * m % n
                    return n if s == n - 1 else s
        elif s == 3:
            if d == 0:
                return 1
            else:
                m = MR_core(d >> 1, a, 3)
                if d & 1:
                    return m * m % n * a % n
                else:
                    return m * m % n
    
    k = 3
    return all(MR_core(n - 1, randrange(2, n)) for i in range(k))


Project Eulerに必要な用語集
http://d.hatena.ne.jp/inamori/20091225/p1