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

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

フェルマーの小定理

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

p を素数とし、ap の倍数でない整数(ap は互いに素)とするときに、

a^{p-1} \equiv 1 \pmod{p}

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

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

  • a = 2、 p = 7 –> 27-1 Mod 7 = 1
  • a = 3、 p = 7 –> 37-1 Mod 7 = 1
  • a = 4、 p = 7 –> 47-1 Mod 7 = 1
  • a = 5、 p = 7 –> 57-1 Mod 7 = 1
  • a = 6、 p = 7 –> 67-1 Mod 7 = 1
  • a = 8、 p = 7 –> 87-1 Mod 7 = 1
  • a = 9、 p = 7 –> 97-1 Mod 7 = 1

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

  • a = 2、 p = 11 –> 211-1 Mod 11 = 1
  • a = 3、 p = 11 –> 311-1 Mod 11 = 1
  • a = 4、 p = 11 –> 411-1 Mod 11 = 1
  • a = 5、 p = 11 –> 511-1 Mod 11 = 1
  • a = 6、 p = 11 –> 611-1 Mod 11 = 1
  • a = 7、 p = 11 –> 711-1 Mod 11 = 1
  • a = 8、 p = 11 –> 811-1 Mod 11 = 1
  • a = 9、 p = 11 –> 911-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をよろしくお願いします。