SICPの問題を解こう13

※Guileだと時間計測にまともな機能が無いのでGaucheに変更。Gaucheだとnanosecondまで計測出来るみたいだしいい感じ。

問題 1.24

素数のテストにFermatテストを使う
Fermatテストの骨格となる以下のコードを見てみる

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

これはbaseのexp乗をmで割った余りを求めてる。
途中ちと直感で納得いかない部分があったので考えてみた。
簡単のためa^bをmで割った余りをa^b(mod m)としておくと上記の偶数部分は

base^exp(mod m) = {base^(exp/2)(mod m)}^2(mod m)

となってる。ん？ほんとにそうかなと思ってちょっと計算

base^(exp/2) = m*k + s ( k,sは整数で s

となるので確かに正しかった。
では本題のコード

(use srfi-27)

(define (search-for-primes n)
  (if (even? n)
      (search-for-primes-iter (+ n 1) 3)
      (search-for-primes-iter n 3)))
(define (search-for-primes-iter x count)
  (if (> count 0)
      (begin
        (timed-prime-test x)
        (search-for-primes-iter (+ x 2)
                                (if (fast-prime? x 10)
                                    (- count 1)
                                    count)))))

(define (timed-prime-test n)
  (start-prime-test n (runtime)))
(define (start-prime-test n start-time)
  (if (fast-prime? n 10)
      (report-prime n (- (runtime) start-time))))
(define (report-prime n elapsed-time)
  (newline)
  (display n)
  (display " *** ")
  (display elapsed-time))

(define (runtime) (time->seconds (current-time)))

(define (fast-prime? n times)
  (cond ((= times 0) #t)
        ((fermat-test n) (fast-prime? n (- times 1)))
        (else #f)))
(define (fermat-test n)
  (define (try-it a)
    (= (expmod a n n) a))
  (try-it (+ 1 (random-integer (- n 1)))))
(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 (square n) (* n n))

(display (search-for-primes 1000))
(display (search-for-primes 1000000))
(display (search-for-primes 1000000000))
(display (search-for-primes 1000000000000))

計測結果

値	通常版の秒数	Fermat版の秒数
1009	0.0129e-4	1.359e-4
1000003	3.55e-4	5.20e-4
1000000007	283.7e-4	8.1e-4
1000000000039	1.9	0.00187

10^6の5.20e-4と10^12の0.00187を比べるとlog(10^12)/log(10^6)=2だから大体2倍くらいの値になってる。ただ10^3と10^6は2倍ではない。。（これは別の要因が関係してるかな？）

問題 1.25

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

baseのexp乗をmで割った余りを求めるのに上記の方法を使った場合、高速テストとして使えるかという問題。

実際にテストをしてみると遅くて使えない。
では理由を考える。

Fermatテストの肝の部分は「ランダムなn以下の整数aのn乗をnで割った余りがaをnで割った余りと等しい」という部分であり、用はa^nをnで割るという処理がこの素数テストの大半の処理時間を占める。ではこの値の２つの求め方を比較してみる。

  (= (expmod a n n) a))
...
(define (expmod base exp m)
  (remainder (fast-expt base exp) m))

  (= (expmod a n n) a))
...
(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))))

前者はn以下のランダムな整数aに対してa^nを全て求めてからnで割っている。a^n自体はfast-exptを使っているのでlog(n)のオーダーで済むがa^nをnで割るのが定時間ではなくなる。(おそらくnのオーダー）
一方後者の方はaをnで割る計算(定時間)をlog(n)回繰り返しているのでlog(n)のオーダーとなる。

問題 1.26

 (remainder (square (expmod base (/ exp 2) m))

正しい処理ではこの部分でlog(n)回になっているが、正しくない処理

 (remainder (* (expmod base (/ exp 2) m)
               (expmod base (/ exp 2) m))

ではexp/2として処理回数が1/2回になっても結局それを2回やっているので回数は減っていない。