この記事は2008年に発表されました。
オリジナル記事を読み、私のニュースレターを購読するには、ここ でご覧ください。
SICPのセクション1.3には素数を判定する方法が記載されています。学生たちに広く知られている単純な素数判定法は、nの最小の約数が1以外で自分自身であれば、nは素数であるというものです。
(define (remainder num divisor)
(- num (* divisor (round (/ num divisor)))))
; nの最小の約数を返す
(define (smallest-divisor n)
; nの約数(2からnまで)を探す
(define (find-divisor n test-divisor)
(cond ((< n (sqrt test-divisor)) n)
((= (remainder n test-divisor) 0) test-divisor)
(else (find-divisor n (+ test-divisor 1)))))
(find-divisor n 2))
; nが素数かどうかをチェックする
(define (prime? n)
(= n (smallest-divisor n )))
しかし、この方法は遅く、複雑度はO(n)です。
SICPでは次にフェルマーの小定理が紹介されています：
もしnが素数であり、a<nの任意のaについて、a^n≡a(mod n)である。
残念ながら、逆定理は成り立ちませんが、それでも有用です。nより小さい数aを選び、上記のテストを行います。もしa^n≡a(mod n)が成り立つなら、nは非常に高い確率で素数です。より多くのaをnに対してテストするほど、nが素数である可能性が高くなります。したがって、これは確率的なアルゴリズムであり、十分なテスト回数を行えば、nが素数でない確率は非常に小さくなります。
これにより、フェルマーのテストプログラムが導かれます。a^nの計算効率を心配する必要はありません：
考えてみましょう：
a^n = (a^(n/2))^2; nが偶数の場合
a^n = a * a^(n-1); nが奇数の場合
これは常に成り立つので、ほぼすべての乗算でnが2倍になります。数直線を見てみると、これは基本的に二分探索です。したがって、時間計算量が減少します。これにより、プログラミング言語が提供するpow関数が自分でa*a*a…と書くよりも常に速い理由が説明されます。
しかし、直接組み込みのpow関数を使用すべきではありません。理由は、我々の目標はaのn乗ではなく、aのn乗をmで割った余りがaの余りに等しいかどうかを確認することだからです。簡単に言えば、この式が成り立つかどうかを確認するだけです：a^n mod m = a。もしpow(a, n) % m == aと直接表現すると、pow(a, n)の値が非常に大きくなる可能性があります。したがって、同様の技法を使用する必要があります：
a^n mod m = ((a^(n / 2))^2) mod m; nが偶数の場合
a^n mod m = (a * a^(n-1)) mod m; nが奇数の場合
この方法では、計算された値はmを超えることはなく、非常に大きな数値での計算を避けることができます。したがって、次のプログラムが簡単に得られます：
(define (remainder num divisor)
(- num (* divisor (round (/ num divisor)))))
(define (expmod base exps m)
(define (square n) (* n n))
(cond ((= exps 0) 1)
((even? exps)
(remainder (square (expmod base (/ exps 2) m))
m))
(else
(remainder (* base (expmod base (- exps 1) m))
m))))
; フェルマーのテスト、nをチェックする。もし(a^n) mod n == aなら、ほぼ間違いなくnは素数
(define (fermat-test n)
(define (try-it a)
(= a (expmod a n n)))
;(= a (remainder (expt a n) n)))
(try-it (+ 1 (random (- n 1)))))
; フェルマーのテストで素数を取得する。nが素数かどうかをチェックする。試行回数timesを試す
(define (fermat-prime? n times)
(cond ((= 0 times) #t)
((fermat-test n) (fermat-prime? n (- times 1)))
(else #f)))
しかし、完全ではありません。
フェルマーのテストを欺くことができるカーマイケル数という種類の数が存在します。2から10000の間に、いくつかのこのような数があります：
561, 1105, 1729, 2465, 2821, 6601, 8911
これらの頑固な数は合成数ですが、どのaを使ってもフェルマーのテストにおいて素数のように振る舞うことができます。彼らは数が少ないわけではなく、広く分布しており、100,000,000の範囲内に255のカーマイケル数があります。したがって、カーマイケル数のリストを維持したくない場合、素数判定方法を強化することができます。
幸運なことに、エクササイズ1.28でミラーとラビンがミラー・ラビンテスト方法をもたらしました。この方法はフェルマーの小定理の変形に基づいています：
nが素数なら、a < nの任意のaについて、a^(n-1) ≡ 1(mod n)
MRテストを使用するには、フェルマーのテストの平方ステップで「非自明な平方根1 mod n」を確認する必要があります。もしs = a^2; とs <> 1およびs <> n-1が存在し、s^2 mod m = 1なら、nは素数ではありません。
Matrix67：ミラー・ラビン素数判定法も確率的アルゴリズムであり、基数aに対するミラー・ラビンテストを通過する合成数を基数aに対する強擬素数と呼びます。基数2に対する最初の強擬素数は2047です。基数2および3に対する最初の強擬素数は1,373,653のように大きいです。
したがって、十分なテスト回数があれば、この方法は工業用途において安全性を保証できます。以下は私の質問1.28に対する答えです：
(define (remainder num divisor)
(- num (* divisor (round (/ num divisor)))))
(define (mr-expmod base exps m)
(define (square x) (* x x))
(define (check-square-root x n)
(define (check r)
(if (and (not (= x 1)) (not (= x (- n 1))) (= r 1)) 0 r))
(check (remainder (square x) n)))
(cond ((= exps 0) 1)
((even? exps)
(check-square-root(mr-expmod base (/ exps 2) m) m))
(else
(remainder (* base (mr-expmod base (- exps 1) m))
m))))
; ミラー・ラビンテスト、nをチェックする。もし(a^(n-1)) mod n == 1なら、ほぼ間違いなくnは素数
(define (mr-test n)
(define (try-it a)
(= 1 (mr-expmod a (- n 1) n)))
(try-it (+ 1 (random (- n 1)))))
; ミラー・ラビンテストで素数を取得する。nが素数かどうかをチェックする。試行回数timesを試す
(define (mr-prime? n times)
(cond ((= 0 times) #t)
((mr-test n) (mr-prime? n (- times 1)))
(else #f)))
いくつかのデータをまとめ、[1000, 1100]、[10000, 10100]、[100000, 100100]、および[1000000, 1000100]の範囲でこれらの方法がどのように機能するかを見てみました：
従来の方法 従来の方法の約数2の最適化 ミラー・ラビンテスト方法
[1000, 1100] 1 ms 1 ms 1 ms
[10000, 10100] 13 ms 8 ms 1 ms
[100000, 100100] 142 ms 78 ms 1 ms
[1000000, 1000100] 1337 ms 830 ms 1 ms
これは私のIntel Core 1 1.83GHzマシン（Ubuntu 8.04 + DrScheme）で実行されました。
また、Cバージョンを書き、以下の結果が得られました：
範囲 [2, 10000]：
$ ./prime
Normal: return 0 cost 10ms
Fermat: return 0 cost 24ms
Miller-Rabin: return 0 cost 28ms
範囲 [200000, 1000000]：
$ ./prime
Normal: return 0 cost 3088ms
Fermat: return 0 cost 824ms
Miller-Rabin: return 0 cost 959ms
gprof分析結果は概ね予想通りであり、レポートの一部を以下に示します：
各サンプルは0.01秒としてカウントされます。
% 累積時間 自己時間 自己呼び出し 合計
時間 秒 秒 呼び出し us/呼び出し us/呼び出し 名前
63.55 2.58 2.58 800000 3.23 3.23 is_prime(int)
15.76 3.22 0.64 800005 0.80 0.84 fermat_expmod(long, long, long)
10.10 3.63 0.41 800012 0.51 0.88 mr_expmod(long, long, long)
7.27 3.92 0.29 14813786 0.02 0.02 chk_unusual_square_root(long, long)
1.23 3.98 0.05 main
0.86 4.01 0.04 800000 0.04 0.93 mr_test(long, long)
0.74 4.04 0.03 14813658 0.00 0.00 square(long)
0.25 4.05 0.01 800000 0.01 0.94 mr_is_prime(int)
0.25 4.06 0.01 800000 0.01 0.85 fermat_is_prime(int)
0.00 4.06 0.00 800000 0.00 0.84 fermat_test(long, long)
小さなn値でのパフォーマンスが低い理由は、Cでの通常の再帰処理のオーバーヘッドによるものであると思います。しかし、データサイズが増えると、このオーバーヘッドは比較的無視できるものになります。誰かが上記のアルゴリズムを反復形式（私にはできません :P）や末尾再帰形式（gccは末尾再帰を反復形式に最適化できます）に変換できれば、パフォーマンスはさらに顕著になると予想されます。