1.2.6 素数判定

(define (smallest-divisor n) (find-divisor n 2))
(define (find-divisor n test-divisor)
  (cond ((> (square test-divisor) n) n)
        ((divisor? test-divisor n) test-divisor)
        (else (find-divisor n (+ test-divisor 1)))))
(define (divisor? a b) (= (remainder b a) 0))
(define (prime? n)
  (= n (smallest-divisor n)))

最も小さい約数を探す手続きsmallest-divisorが$n$をそのまま返して来たら、$n$は素数ということになる。

$n$の約数は$\sqrt{n}$まで探して見つからなければ打切り。よってこのアルゴリズムの時間計算量の増加オーダーは$\Theta(\sqrt{n})$。これじゃ大きな数に対しては使えない

フェルマーテスト

$n$が素数ならば、$a$が$n$より小さい正の整数であるとき、$a^n$を$n$で割った余りは$a$と等しい。

等しくなければ$n$は素数ではない。だが等しい場合は、素数である可能性が高い。このテストを繰り返し行えば、素数かどうかの確信を高めることができる。

(use srfi-27)

(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (remainder
           (square (expmod base (/ exp 2) m))
           m))
        (else
          (remainder
            (* base (expmod base (- exp 1) m))
            m))))

(define (fermat-test n)
  (define (try-it a)
    (= (expmod a n n) a))
  (try-it (+ 1 (random-integer (- n 1)))))

(define (fast-prime? n times)
  (cond ((= times 0) #t)
        ((fermat-test n) (fast-prime? n (- times 1)))
        (else #f)))

この方法は、確実な結果を返すわけではない。素数でないにも関わらず、フェルマーテストをすり抜けるカーマイケル数というのが知られている。だが桁が増えるとそうした数は非常にまれ。

テストを繰り返すことで、間違いの可能性を好きなだけ小さくできる。そうしたアルゴリズムを確率的アルゴリズムという。RSAはこれを応用している。世の中の役に立ってるのだ。

練習問題1.21

smallest-divisorを動かしてみるというだけの問題。

gosh> (smallest-divisor 199)
199
gosh> (smallest-divisor 1999)
1999
gosh> (smallest-divisor 19999)
7

199、1999は素数だった。

練習問題 1.22

手続きprime?の計算時間の増加オーダーは$\Theta(\sqrt{n})$。ほんとかどうか実際試してみよう。

指定した範囲の連続した奇数について素数判定する手続きを作れという問題だったけど、ちょっとアレンジして、3つの奇数が見つかったら終了としてみた。

計測にはgauche.timeを利用。

.scm

(use gauche.time)
(define (if-prime-report n)
  (let ((tc (make <real-time-counter>)))
    (if (with-time-counter tc (prime? n))
      (begin (format #t "~D *** ~D\n" n (time-counter-value tc)) 1)
      0)))

(define (search-for-prime n)
  (define (next n)
    (+ n (if (= (remaider n 2) 0) 1 2)))
  (define (search n c)
    (if (= c 3)
      c
      (search (next n) (+ c (if-prime-report n)))))
  (search n 0))

1000、10000、100000に対して試せと書いてあったけど、数が小さすぎて測定不能だった。$10^{10}$、$10^{11}$、$10^{12}$に対しそれぞれ実行。

.scm

gosh> (search-for-prime (expt 10 10))
10000000019 *** 0.021751
10000000033 *** 0.020621
10000000061 *** 0.02121
3
gosh> (search-for-prime (expt 10 11))
100000000003 *** 0.072056
100000000019 *** 0.068461
100000000057 *** 0.065637
3
gosh> (search-for-prime (expt 10 12))
1000000000039 *** 0.196225
1000000000061 *** 0.191329
1000000000063 *** 0.199522
3

$\sqrt{10}=3.162...$なので、まず予測通りと言える。

練習問題 1.23

find-divisorは無駄が多い。2以外の偶数についてはチェックしなくていいじゃないか。

.scm

(define (smallest-divisor n) (find-divisor n 2))
(define (find-divisor n test-divisor)
  (define (next divisor) (+ divisor (if (= divisor 2) 1 2)))
  (cond ((> (square test-divisor) n) n)
        ((divisor? test-divisor n) test-divisor)
        (else (find-divisor n (next test-divisor)))))
(define (divisor? a b) (= (remainder b a) 0))
(define (prime? n)
  (= n (smallest-divisor n)))

この方法ならば、さっきの倍の速度で実行できるんじゃないか。やってみよう。

.scm

gosh> (search-for-prime (expt 10 10))
10000000019 *** 0.01474
10000000033 *** 0.011265
10000000061 *** 0.012591
3
gosh> (search-for-prime (expt 10 11))
100000000003 *** 0.037631
100000000019 *** 0.042453
100000000057 *** 0.031286
3
gosh> (search-for-prime (expt 10 12))
1000000000039 *** 0.113101
1000000000061 *** 0.102724
1000000000063 *** 0.10425
3

倍、ではないな。せいぜい1.7倍。nextの呼び出しの分だけ時間がかかるため、と思われる。

練習問題 1.24

先ほど見つけた素数に対し、fast-prime?でテストしてみる。
テスト回数は適当に、1000回としてみた。

.scm

(use gauche.time)
(define (prime-report n)
  (let ((tc (make <real-time-counter>)))
    (if (with-time-counter tc (fast-prime? n 1000))
        (begin (format #t "~D *** ~D\n" n (time-counter-value tc)) #t)
        #f)))

この手続きの計算時間の増加オーダーは$\Theta(\log{n})$なので、
$10^{10}$付近の素数のテストにかかる時間を基準とすると、
$10^{11}$付近の素数については $\frac{\log{10^{11}}}{\log{10^{10}}} = 1.1$倍、
$10^{12}$付近の素数については $\frac{\log{10^{12}}}{\log{10^{10}}} = 1.2$倍、
...
の計算時間がかかることが予想される。

以下、実行結果。

.scm

gosh> (prime-report 10000000019)
10000000019 *** 0.031272
# t
gosh> (prime-report 10000000033)
10000000033 *** 0.02699
# t
gosh> (prime-report 10000000061)
10000000061 *** 0.029298
# t
gosh> (prime-report 100000000003)
100000000003 *** 0.03167
# t
gosh> (prime-report 100000000019)
100000000019 *** 0.03104
# t
gosh> (prime-report 100000000057)
100000000057 *** 0.0336
# t
gosh> (prime-report 1000000000039)
1000000000039 *** 0.0322
# t
gosh> (prime-report 1000000000061)
1000000000061 *** 0.03624
# t
gosh> (prime-report 1000000000063)
1000000000063 *** 0.032391
# t

$10^{10}$付近の平均は0.02919。
$10^{11}$付近の平均は0.03210。増加率は1.0997。
$10^{12}$付近の平均は0.03361。増加率は1.1514。

おまけでもうひとつ上の桁を。

.scm

gosh> (prime-report 10000000000037)
10000000000037 *** 0.038639
# t
gosh> (prime-report 10000000000051)
10000000000051 *** 0.037296
# t
gosh> (prime-report 10000000000099)
10000000000099 *** 0.036793
# t

$10^{13}$付近の平均は0.03757。増加率は1.2870。

練習問題 1.25

expmod手続きを以下のように実装してはダメなのか？

.scm

(define (expmod base exp m)
  (remainder (expt base exp) m))

これでは、$a^n$を計算してしまう。例えば$10000$という値の素数判定を行うなら、$a^{10000}$を計算しようとするのだ。これはまずい。

元のexpmodを再掲する。

.scm

(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (remainder
           (square (expmod base (/ exp 2) m))
           m))
        (else
          (remainder
            (* base (expmod base (- exp 1) m))
            m))))

（原書の注46によると）これは、$xy \div m$の余りを求めるには、$x \div m$の余りと、$y \div m$の余りを掛けた値を、さらに$m$で割ってその余りを計算すればよいという事実に基づいている。

$n$が偶数なら$a^n$を $a^{n/2} \times a^{n/2}$ に分割し、
奇数なら $a^{n-1} \times a$ に分割している。