あけましておめでとうございます

パッケージJava製品開発担当の大です。本年もよろしくお願いします。

さて、今年は平成23年、西暦2011年です。23といえば素数ですが、なんとなく2011も素数みたいな気がしますね。しかし、素数列の最初のほういくつかは覚えてても、4桁になるとちょっと覚えてません。そこで、今年のプログラム書初めは素数判定をやることにします。

フェルマーの小定理

素数の性質についての定理のひとつとして、フェルマーの小定理というのがあります。

p を素数とし、a を p の倍数でない整数（a と p は互いに素）とするときに、

すなわち a を p - 1 乗したものを p で割ったあまりは 1 になるというもの。

ちょっと試してみましょうか。素数7を例にとってみると、

• a = 2、 p = 7 --> 2^7-1 Mod 7 = 1
• a = 3、 p = 7 --> 3^7-1 Mod 7 = 1
• a = 4、 p = 7 --> 4^7-1 Mod 7 = 1
• a = 5、 p = 7 --> 5^7-1 Mod 7 = 1
• a = 6、 p = 7 --> 6^7-1 Mod 7 = 1
• a = 8、 p = 7 --> 8^7-1 Mod 7 = 1
• a = 9、 p = 7 --> 9^7-1 Mod 7 = 1
• ...

素数11を例にとってみると、

• a = 2、 p = 11 --> 2^11-1 Mod 11 = 1
• a = 3、 p = 11 --> 3^11-1 Mod 11 = 1
• a = 4、 p = 11 --> 4^11-1 Mod 11 = 1
• a = 5、 p = 11 --> 5^11-1 Mod 11 = 1
• a = 6、 p = 11 --> 6^11-1 Mod 11 = 1
• a = 7、 p = 11 --> 7^11-1 Mod 11 = 1
• a = 8、 p = 11 --> 8^11-1 Mod 11 = 1
• a = 9、 p = 11 --> 9^11-1 Mod 11 = 1
• ...

なるほど、たしかに余りは必ず1になります。

フェルマーテスト

この定理の対偶をとると、フェルマーテストというものになります。

これは、素数判定の高速なアルゴリズムとして知られています（ただし、確実ではないというところがミソですが、それは後述）。

SICPに出てくるフェルマーテストのコードをPostScriptで写経してみると、こんな感じになりました（Wikipediaの記述に合わせて若干変えてます）。

• is-prime.ps

realtime srand

% n m random
% n以上m未満の範囲の乱数を発生させます
/random {
    0 dict begin
    /m exch def
    /n exch def
    rand m n sub mod n add
    end
} def

% base exp mod expmod
% base^expをmで割った余りを取得します
/expmod {
    0 dict begin
    /m exch def
    /exp exch def
    /base exch def
    exp 0 eq {
        1
    }{
        base exp
        exp 2 mod 0 eq {
            2 idiv m expmod dup
        } {
            1 sub m expmod base
        } ifelse
        mul m mod
    } ifelse
    end
} def

% n fermat-test
% nに対してフェルマーテストを行います
/fermat-test {
    0 dict begin
    /n exch def
    /a 2 n random def
    a n 1 sub n expmod 1 eq
    end
} def

% n times fast-prime?
% nが素数かどうかをtimes回フェルマーテストを実行して判定します
/fast-prime? {
    0 dict begin
    /times exch def
    /n exch def
    times 0 eq {
        true
    } {
        n fermat-test {
            n times 1 sub fast-prime?
        } {
            false
        } ifelse
    } ifelse
    end
} def

% n is-prime?
% nが素数かどうかを判定します
/is-prime? {
    10 fast-prime?
} def

動作確認

実際に動かして確認してみましょう。そうですね、とりあえず3～100までの自然数をすべて素数かどうか判定して、素数だけ抜き出してみましょうか。

$ gs -dNODISPLAY -q  is-prime.ps
GS> [3 1 100 {dup is-prime? not {pop} if} for] ==
[3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97]

判定できてるっぽいですね。では、今回の目的、2011が素数かどうか判定してみます。

GS> 2011 is-prime? ==
true

素数だそうです！たぶん！

確率的素数判定法

「たぶん」といいました。実は、この「フェルマーテスト」は確率的素数判定法と呼ばれる素数判定法で、ある数が「素数である」と確実に判定することはできないのです。これは、フェルマーテストの意味を考えればわかります。自然数 n が合成数であるための十分条件ではありますが、必要条件ではないのですね。なので、上のコードでは、フェルマーテストを10回、適当なaを決めて実行し、「毎回素数と判定されたなら素数」とすることで、多少確率を上げています。

しかしそれでも、フェルマーテストを高い確率で騙してしまう合成数もあります。カーマイケル数と呼ばれるものがそうですが、それについてはまた別の機会に。。。

それでは、たぶん素数な2011年もHOSをよろしくお願いします。