オンラインSICP読書女子会 #5(1.2.4~1.2.6 最大公約数と素数判定) まとめ

#ladiescpp の主催してる オンラインSICP読書女子会 #4 (1.2.4~1.2.6) - connpassのまとめです〜。

指数計算（前回続き）

練習問題1.19

fib(n+1)からfib(n)がわかるなら、
fib(n+2)もfib(n)から求められるはず！ みたいな問題。

\begin{align}
a_{n+1} &= b_nq + a_nq + a_np \\
b_{n+1} &= b_np + a_nq
\end{align}

とおく。

\begin{align}
a_{n+2} &= (b_np + a_nq)q + (b_nq + a_nq + a_np)q + (b_nq + a_nq + a_np)p \\
b_{n+2} &= (b_np + a_nq)p + (b_nq + a_nq + a_np)q
\end{align}

なので、展開していくと

\begin{align}
\qquad &a_{n+2} \\
&\quad = (b_np + a_nq)q + (b_nq + a_nq + a_np)q + (b_nq + a_nq + a_np)p \\
&\quad = a_n(2pq + q^2) + b_n(2pq + q^2) + a_n(p^2 + q^2)\\
\\
\qquad &b_{n+2}\\
&\quad = (b_np + a_nq)p + (b_nq + a_nq + a_np)q \\
&\quad = a_n(2pq + q^2) + b_n(p^2 + q^2)
\end{align}

最初の式と見比べて

\begin{align}
p' &= p^2 + q^2 \\
q' &= 2pq + q^2
\end{align}

元の関数にこれを代入

(define (fib n) (fib-iter 1 0 0 1 n))
(define (fib-iter a b p q count)
  (cond ((= count 0) b)
    ((even? count)
     (fib-iter a
           b
           (+ (* p p) (* q q))
           (+ (* 2 p q) (* q q))
           (/ count 2)))
    (else (fib-iter
           (+ (* b q) (* a q) (* a p))
           (+ (* b p) (* a q))
           p
           q
           (- count 1)))))

(fib 10) ; 55 あってる！

最大公約数

練習問題1.20

以下@hio-harp さんの解答引用 引用元はこちら

正規順序評価規則の場合

; [入力]
(gcd 206 40)

; [step 1]
; a1=206, b1=40. <env:none>
(if (= b1 0)
    a1
    (gcd b1 (remainder a1 b1)))
; <eval b1> <env:a1=206, b1=40>
b1
40
; pred ==> false.
else-clause
(gcd b1 (remainder a1 b1))

; [step 2]
; a2=b1, b2=(remainder a1 b1) <env:a1=206, b1=40>
(if (= b2 0)
    a2
    (gcd b2 (remainder a2 b2)))
; <eval b2> <env:a1=206, b1=40; a2=b1, b2=(remainder a1 b1)>
b2
(remainder a1 b1)
(remainder 206 40)
; REMAINDER 1 回目.
6 ; <env:a1=206, b1=40, a2=b1, b2=6>
; pred = false
else-clause
(gcd b2 (remainder a2 b2))

; [step 3]
; a3=b2, b3=(remainder a2 b2) <env:a1=206, b1=40, a2=b1, b2=6>
(if (= b3 0)
    a3
    (gcd b3 (remainder a3 b3)))
; <eval b3> <env:a1=206, b1=40; a2=b1, b2=6; a3=b2, b3=(remainder a2 b2)>
b3
(remainder a2 b2)
(remainder b1 6)
(remainder 40 6) ; <env:a1=206, b1=40; a2=40, b2=6; a3=b2, b3=(remainder a2 b2)>
; REMAINDER 2 回目.
4 ; <env:a1=206, b1=40; a2=40, b2=6; a3=b2, b3=4>
; pred = false
else-clause
(gcd b3 (remainder a3 b3))

; [step 4]
; a4=b3, b4=(remainder a3 b3) <env:a1=206, b1=40; a2=40, b2=6; a3=b2, b3=4>
(if (= b4 0)
    a4
    (gcd b4 (remainder a4 b4)))
; <eval b4> <env:a1=206, b1=40; a2=40, b2=6; a3=b2, b3=4; a4=b3, b4=(remainder a3 b3)>
b4
(remainder a3 b3)
(remainder b2 4)
(remainder 6 4) ; <env:a1=206, b1=40; a2=40, b2=6; a3=6, b3=4; a4=b3, b4=(remainder a3 b3)>
; REMAINDER 3 回目.
2 ; <env:a1=206, b1=40; a2=40, b2=6; a3=6, b3=4; a4=b3, b4=2>
; pred = false
else-clause
(gcd b4 (remainder a4 b4))

; [step 5]
; a5=b4, b5=(remainder a4 b4) <env:a1=206, b1=40; a2=40, b2=6; a3=6, b3=4; a4=b3, b4=2>
(if (= b5 0)
    a5
    (gcd b5 (remainder a5 b5)))
; <eval b5> <env:a1=206, b1=40; a2=40, b2=6; a3=6, b3=4; a4=b3, b4=b3; a5=b4, b5=(remainder a4 b4)>
b5
(remainder a4 b4)
(remainder b3 2)
(remainder 4 2) ; <env:a1=206, b1=40; a2=40, b2=6; a3=6, b3=4; a4=2, b4=2; a5=b4, b5=(remainder a4 b4)>
; REMAINDER 4 回目.
0 ; <env:a1=206, b1=40; a2=40, b2=6; a3=6, b3=4; a4=2, b4=2; a5=b4, b5=0>
; pred = true
then-clause
a5
; <eval a5> <env:a1=206, b1=40; a2=40, b2=6; a3=6, b3=4; a4=2, b4=2; a5=b4, b5=0>
a5
b4
2 ; answer.

; remainder 演算は合計4回.

適用順序評価規則の場合

; [入力]
(gcd 206 40)

; [step 1]
; a1=206, b1=40. <env:none>
(if (= b1 0)
    a1
    (gcd b1 (remainder a1 b1)))
; <eval b1> <env:a1=206, b1=40>
b1
40
; pred ==> false.
else-clause
(gcd b1 (remainder a1 b1))
(gcd 40 (remainder 206 40))
; REMAINDER 1 回目.
(gcd 40 6)

; [step 2]
; a2=40, b2=6
(if (= b2 0)
    a2
    (gcd b2 (remainder a2 b2)))
; <eval b2> <env:a2=40, b2=6>
b2
6
; pred ==> false.
else-clause
(gcd b2 (remainder a2 b2))
(gcd 6 (remainder  40 6))
; REMAINDER 2 回目.
(gcd 6 4)

; [step 3]
; a3=6 b3=4
(if (= b3 0)
    a3
    (gcd b3 (remainder a3 b3)))
; <eval b3> <env:a3=6, b3=4>
b3
4
; pred ==> false.
else-clause
(gcd b3 (remainder a3 b3))
(gcd 4 (remainder  6 4))
; REMAINDER 3 回目.
(gcd 4 2)

; [step 4]
; a4=4 b4=2
(if (= b4 0)
    a4
    (gcd b4 (remainder a4 b4)))
; <eval b4> <env:a4=4, b4=2>
b4
2
; pred ==> false.
else-clause
(gcd b4 (remainder a4 b4))
(gcd 2 (remainder  4 2))
; REMAINDER 4 回目.
(gcd 2 0)

; [step 5]
; a5=2 b5=0
(if (= b5 0)
    a5
    (gcd b5 (remainder a5 b5)))
; <eval b5> <env:a5=2, b5=0>
b5
0
; pred ==> true
then-clause
a5
; <eval a5> <env:a5=2, b5=0>
a5
2 ; answer.

; こちらの場合でも remainder 演算は合計4回.

素数判定

フェルマーの小定理 

nが素数, n < a (aは整数)
a^n mod n = a mod n

をaが素数の時に満たす
aが素数じゃなくてもたまに満たす

1回フェルマーテスト(times=1)するのの計算時間は O(log(n))

練習問題1.21

smallest-divisor手続きを使って次の数値の最小の約数を求めよ。199, 1999, 19999.

;;smallest-divisor.scm
(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 (+ test-divisor 1)))))

(define (divides? a b) (= (remainder b a) 0))

(define (prime? n)
  (= n (smallest-divisor n)))

(load "./smallest-divisor.scm")

(smallest-divisor 199) ;199
(smallest-divisor 1999) ;1999
(smallest-divisor 19999) ;7

練習問題1.22

n~mまでの数字が素数かどうかみていって、
素数だったらそれを判定するのに何マイクロ秒かかったかを調べる関数
fermatテストを使わずに、普通にnに対して√nまでで割れるかで素数判定をする。

;; 1.22.scm
(load "./smallest-divisor.scm")

(define (timed-prime-test n)
  (newline)
  (display n)
  (start-prime-test n (runtime)))

(define (start-prime-test n start-time)
  (if (prime? n)
      (report-prime (- (runtime) start-time))))

(define (report-prime elapsed-time)
  (display " *** ")
  (display elapsed-time))

(define (runtime)
  (use srfi-11)
  (let-values (((a b) (sys-gettimeofday)))
    (+ (* a 1000000) b))) ;; gosh にruntime実装がない

(define (search-for-primes n end)
  (cond ((even? n) (search-for-primes (+ n 1) end))
    ((> n end) (display "end"))
    (else (timed-prime-test n) (search-for-primes (+ n 2) end))))

(search-for-primes 1000 1050) ; 1009 1013 1019 ｜6
(search-for-primes 10000 10100) ; 10007 10009 10037 ｜17
(search-for-primes 100000 100100) ; 100003 100019 100043｜ 55 
(search-for-primes 1000000 1000030) ; 395
(search-for-primes 10000000 10000030) ; 797
; √10000          100: 17
; √1000000            1000: 395

O(√n)に対して 妥当なのかどうかは良くわからない
100:17 ≒ 1000: 395 といえるか？

※何度かテストすると、実行時間2倍位計算時間にばらつきがあるので怪しい。

練習問題1.23

1.22では、nに対して1 ~ √nまでで割れるか調べていたが、
2以上の偶数で割っても仕方がない。そこで2以上の偶数は飛ばすようにする。
そうすると2倍早くなるはず

たしかに2倍くらい早くなった

(load "./1.22.scm")

(define (find-divisor n test-divisor)
  (cond ((> (square test-divisor) n) n)
    ((divides? test-divisor n) test-divisor)
    (else (find-divisor  n (+ test-divisor 1)))))

(search-for-primes 1000 1020) ; 1009 1013 1019 ｜9マイクロ秒
(search-for-primes 10000 10138) ; 10007 10009 10037 ｜26
(search-for-primes 100000 100044) ; 100003 100019 100043 ｜65 

(define (next x) (if (= x 2) 3 (+ x 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)))))


(search-for-primes 1000 1020) ; 1009 1013 1019 ｜4
(search-for-primes 10000 10138) ; 10007 10009 10037 ｜11
(search-for-primes 100000 100044) ; 100003 100019 100043 ｜36 

練習問題1.24

フェルマーの小定理を使って素数を求めてみる。

(use srfi-27) ;; randam-integerが入ってるmodule
(load "./1.22.scm")

;;ばらつくので3回ずつ実行
(search-for-primes 1000 1020) ; 6, 8, 6
(search-for-primes 10000 10040) ; 15, 26, 30
(search-for-primes 100000 100043) ; 47, 65, 54
(search-for-primes 1000032 1000033) ; 395, 185, 268
(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)))

(define fermet-time 3)
(define (start-prime-test n start-time)
  (if (fast-prime? n fermet-time)
      (report-prime (- (runtime) start-time))))

(search-for-primes 1000 1020) ; 10, 18, 29 (ふえた！)
(search-for-primes 10000 10040) ; 23, 20, 40
(search-for-primes 100000 100043) ; 29, 27, 24

(define fermet-time 10)
(search-for-primes 1000 1020) ; 21, 15, 37
(search-for-primes 10000 10040) ; 19, 22, 21
(search-for-primes 100000 100043) ; 28, 29, 19
;; あんまり増えない
(search-for-primes 1000032 1000033) ; 74, 52, 37
;; 元のprime?よりは大きくした時にはやい

fermatテストはO(log(n))で増加するので

>>> math.log(10000, 2)
13.28771237954945
>>> math.log(1000000, 2)
19.931568569324174

1.5倍くらいかかるはず。
実際は以下なので、妥当な結果？
1000のとき(21, 15, 37) -> 1000000のとき(74, 52, 37)

もっと大きい値をみてみる

(define fermet-time 10)
(search-for-primes 10000 10040) ; 19, 22, 21 ※log(10000) = 13.287...

(search-for-primes 1000000000 1000000033) ; 162, 111, 119 ※log(1000000000) = 26.57..
;;10000の2倍以上はある

log(n)の値は2倍なのに、実行時間は5倍くらい有るのが気になる・・・。
原因はまだ不明。