ウォンツテック

そでやまのーと

SICP 問題を解こう12

久々にSICPの問題を解きたくなったのでやってみる。
問題 1.23
問題1.22の素数テストを改良して検査時間が1/2倍になったかどうかを調べてなってなかったら理由を説明せよという問題。
改良の仕方は素数のテストに\sqrt{n}までの2以上の整数を調べていたのに対し、2と3以上の奇数を調べるという方法に変更というもの。
とりあえずコードを書いて検査時間を調べてみる。

(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 (prime? x)
                                    (- count 1)
                                    count)))))
(define (timed-prime-test n)
  (newline)
  (display n)
  (start-prime-test n (get-internal-run-time)))
(define (start-prime-test n start-time)
  (if (prime? n)
      (report-prime (- (get-internal-run-time) start-time))))
(define (report-prime elapsed-time)
  (display " *** ")
  (display elapsed-time))

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

値1,000,000,000以上の素数を3つ求めるテストで問題1.22の物とほぼ同じ。
違うのは

  (define (next n)
    (if (= n 2)
        3
        (+ n 2)))

の部分。実際に掛かった時間は問題1.22の物だと16〜17単位時間(秒じゃない)で問題1.23の物だと10〜11単位時間。よって大体1.6〜1.7倍程度。
値10,000,000,000以上の3つの素数でもテストしてみると

問題1.22 問題1.23
1,000,000,000 16〜17 10〜11
10,000,000,000 56〜57 35〜36

といった感じで大体1.6倍程度。
ではこれがなぜ2倍じゃないのかってのが問題。

実際に偶数の検査に掛かる時間をT_e,奇数の検査に掛かる時間をT_oとするとこのテストによると (T_e + T_o)/T_o = 1.6 くらいになってる。
ってことはT_eの検査時間がT_oより短いって事。
そんじゃあこの検査に掛かる主要部分はどこか見てみると

    ((divides? test-divisor n) test-divisor)

の部分であり、これは素数nがtest-divisorで割れるかどうかを見ている。
当然T_eの時はtest-divisorが偶数である。
ここでnは素数であるから奇数なのでこのdivides?(実質の処理は(remainder n test-divisor))は即座に結果が返ってくるはずである。
T_oはどうかというとdivides?の処理は実際に割って見ないとわからない。
remainderの実装に依存するが、このschemeの処理系においてはT_e < T_oであり、求める倍率が2倍ではなく1.6倍程度となる。
※1.6〜1.7倍は本質的な数値ではなく単純に10^12〜10^13のオーダーにおけるT_oに掛かる時間が表出してきている値だと思う(T_eはおそらく一定なので)