SICP読書録 #9 1.2.6 素数判定 - 練習問題 1.26〜1.28

hiroshi-manabe様の日本語訳を使わせてもらいます。
練習問題等はGaucheで実行。

前回はこちら

フェルマーテストに関する問題が続く。コードを再掲しとこう。

.scm

(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)))

練習問題 1.26

expmod手続きを以下のように書いたら、計算時間の増加オーダーが$\Theta(\log{n})$から$\Theta(n)$になってしまった。なぜか？

.scm

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

expmodを2回呼んでるのが、オリジナルと異なる点。

呼び出しの階層そのものは $\log_2{n}$ 回程度だけど、そのたびにexpmodの呼ばれる数は2回 > 4回 > 8回と、2倍ずつ増えていく。というわけで、増加オーダーは

 \Theta(\log_2{n^2})=\Theta(n)

となる。

練習問題 1.27

この手続きは、確実な素数判定を行うものではない。
テスト回数timesを増やすことで判定の制度を高めることができるのだが・・・カーマイケル数という数は素数でないのに変わらず、このテストをすり抜けてしまう。

カーマイケル数を最小のものからいくつか挙げると、561, 1105, 1729, 2465, 2821, 6601。
実際試してみよう。

ランダムな数ではなく、$a<n$となる全ての正の整数に対してテストするために以下の手続きを用意する。

.scm

(define (carmichael-test n)
  (define (iter a)
    (cond ((= a n) #t)
          ((= (expmod a n n) a) (iter (+ a 1)))
          (else #f)))
  (iter 1))

実行結果。

.scm

gosh> (carmichael-test 561)
# t
gosh> (carmichael-test 1105)
# t
gosh> (carmichael-test 1729)
# t
gosh> (carmichael-test 2465)
# t
gosh> (carmichael-test 2821)
# t
gosh> (carmichael-test 6601)
# t

練習問題 1.28

カーマイケル数に騙されないテスト方法がある。

これはミラー-ラビンテストと呼ばれ、2乗を行うステップで、1でもn-1でもなく、しかも法mに対して1と合同な数の平方根（つまり2乗してmで割った余りが1になる数)が見つかったら素数ではない、とするものだ。

.scm

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

(define (fermat-test n)
  (define (try-it a)
    (= (mr-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)))

実行結果。もう騙されないぞ。

.scm

gosh> (fast-prime? 561 100)
# f
gosh> (fast-prime? 1105 100)
# f
gosh> (fast-prime? 1729 100)
# f
gosh> (fast-prime? 2465 100)
# f
gosh> (fast-prime? 2821 100)
# f
gosh> (fast-prime? 6601 100)
# f

他の数に対してもテストしなければ。
前回作成した、指定した数から順に素数判定する手続きを以下のように修正し、試してみる。

.scm

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

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

実行すると、前回のと同じ結果になった。

.scm

gosh> (search-for-prime (expt 10 10))
10000000019 *** 0.005072
10000000033 *** 0.006049
10000000061 *** 0.004751
3
gosh> (search-for-prime (expt 10 11))
100000000003 *** 0.013715
100000000019 *** 0.006009
100000000057 *** 0.006815
3
gosh> (search-for-prime (expt 10 12))
1000000000039 *** 0.005402
1000000000061 *** 0.007774
1000000000063 *** 0.003905
3