素数判定アルゴリズム

The Haskell Road to Logic, Maths, and Programming( http://fldit-www.cs.uni-dortmund.de/~peter/PS07/HR.pdf )のp3~p8にある、素数判定プログラムの解説に対して感じた疑問を書き留めておく。

以下、\(n\)は常に2以上の整数、すなわち素数か合成数かのいずれか一方であるとする。

教科書では実装に先立ち、まず\(LD(n)\)という関数を定義し、2つの命題を証明している。\(LD(n)\)とは、「\(n\)の約数(1以外)のうち最小のもの」である。

\(n\) 2 3 4 5 6 7 8 9 10 11 12 13 14 ……
\(LD(n)\) 2 3 2 5 2 7 2 3 2 11 2 13 2 ……

この\(LD(n)\)に対して、以下の1,2が成り立つ。
1. \(n\)が素数か合成数かにかかわらず、\(LD(n)\)は素数である。
2. \(n\)が素数でない(合成数である)とき、\([LD(n)]^2\leq n\)である。

1は、表の下段には素数だけが並ぶことを主張している。この命題を素数判定にどのように活用するのだろうかと疑問を覚えたが、後のページを読んでも、どうも使っている気配がない。実際には「\(n\)が素数である ⇔ \(LD(n) = n\)」という同値関係を用いている。「命題1を用いれば、左向きの矢印が言える」と考えることもできるが、そもそもこの同値関係は素数の定義と\(LD(n)\)の定義から直接言えることである。

教科書での実装をいったん忘れて、とりあえず上の関係を使って素朴に素数判定プログラムを書けば

ldf k n | rem n k == 0 = k
        | otherwise    = ldf (k+1) n
prime n = ldf 2 n == n

となる*1。この関数ldfは、「nの約数のうち、k以上で最小のもの」ということができ、その\(k=2\)の場合が\(LD(n)\)である。

このアルゴリズムでは、nを割り切る数が見つかるまで、2から始めて最悪でnまで試さないといけない。しかし、もしもnが合成数であり、\(n=1\times n\)以外に\(n=a\times b\)(ただし\(1\lt a\leq b \lt n\))という2整数の積に分解する方法があるなら、\(a^2\leq ab=n\)だから、\(\sqrt{n}\)以下の範囲でnを割り切る数が見つかるはずである。これは同時に、「その範囲で見つからなければ諦めて探索を切り上げ、\(LD(n)=n\)である(したがって\(n\)は素数である)と結論してよい」ということを意味する。このような「便法」を提供してくれるのが命題2にほかならず、これを盛り込んで実装し直すと

ldf k n | k^2 > n      = n
        | rem n k == 0 = k
        | otherwise    = ldf (k+1) n
prime n = ldf 2 n == n

となる。これでアルゴリズムの効率は向上したが、もはやldfは「nの約数のうち、k以上で最小のもの」を返す関数ではない*2*3。例えば「ldf 6 35」は35を返すが、「35の約数のうち、6以上で最小のもの」は7である。「ldf 2 n」という形で呼ぶ限り\(LD(n)\)を返すことに変わりはないが、ldfだけを切り出して、一般のkに対して「これは何をする関数であるか」という問にクリアカットに答えることは難しくなってしまった。

どうしてこんな、微妙に筋の悪い話に迷い込んでしまったのか考え直すと、そもそも\(LD(n)\)という関数を立てること自体に問題があるようにも思える。見方によっては、この\(LD(n)\)は\(n\)が合成数か素数かによってずいぶん意味合いが異なる。\(n\)が合成数のときの\(LD(n)\)とは、前述の\(n=a\times b\)(ただし\(1\lt a\leq b\lt n\))の\(a\)、すなわち「大きくないほう」に入ることのできる値のひとつ(最小値)である。いっぽう\(n\)が素数のときの\(LD(n)\)は\(n\)だが、これは\(n=1\times n\)という唯一の分解方法における、2数のうちの「大きいほう」である。\(LD(n)\)の値で\(n\)の素数判定をしたいのに、\(n\)が素数かどうかで\(LD(n)\)の意味合いがこんなにも変わるのは、なんだか先に勝負が付いてしまっているような気持ち悪さがある。そう考えると、\(n\)が素数のときの\(LD(n)\)の値は、\(n\)というよりもむしろNULL的なものを想像したほうがよいのかもしれない。するとldfは「\(k\)以上\(√n\)以下の範囲に\(n\)の約数が存在すればその最小値を、存在しなければNULLを返す」関数と言うことができる。これを\(k=2\)で呼び出せば、合成数のときは\(a\)の最小値が、素数のときはNULLが返ってくる。

ここまで考えると、素数判定に関する限り「ldfがNULLを返すかどうか」だけが問題であり、\(a\)の最小値の具体的な値はどうでもよい。そう割り切るなら、

nodiv k n | k^2 > n      = True
          | rem n k == 0 = False
          | otherwise    = nodiv (k+1) n
prime n = nodiv 2 n

と書くこともできる。ここでのnodivは、「k以上√n以下の範囲にnの約数が存在しない」ことを表す。

*1:ただし、教科書にあるような関数divides,ldの定義や不適切なnへの処理は省略した。以下同様。

*2:教科書では、この効率化されたldfに対して「nの約数のうち、k以上で最小のものを表す」と誤って書かれていたが、のちに「nがk未満に約数(1以外)を持たない限りにおいて、nの約数(1以外)で最小のものを表す」と訂正されている( http://homepages.cwi.nl/~jve/HR/errata.pdf )。訂正後の記述は確かに正しいけれども、きわめて理解しにくいものとなっている。「限りにおいて」って、それ以外のときは何を意味するねんと。

*3:正確には教科書での実装ではrem n k == 0の判定をk^2 > nの判定よりも先に行なっている。そのため、例えばldf 6 35は35を返すのにldf 7 35は(35ではなく)7を返すという仕様になっている。いずれにせよ、著者による訂正の記述はそのまま成り立ち、素数判定としての挙動にも影響を与えない。どうせldf 2 35の形でしか使わないから、再帰呼び出しはldf 6 35で終了し、ldf 7 35がどんな値を返そうと知ったこっちゃないのである。