カテゴリー : Source

SICP Reading #10 高階関数な練習問題を解く

最近

研究室を決める時期だったりそういう感じで皆殺伐としていたり、課題ヤラで結構忙しくなってきた時期。SICPを読んでる時間は楽しい(ならもっと時間取れよという)。

元気に問題解いていきます。

問題1.29

前回実装した近似式による積分の実装integralよりも、正確な積分近似のSimpsonの公式を高階関数で実装しろというもの。前回実装したintegralは四角形で領域を近似していたけど、Simpsonの公式では積分区間を分割し、各3点を2次関数で近似するといもの(用は、弧を用いた近似)。

大学での並列分散システムでMPIを利用したSimpsonの公式の実装をした矢先、SICPで問題として取り上げられててちょっぴり嬉しかったり。

以下、実装

(define (cube x) (* x x x))

;;級数の和を抽象化
(define (sum term a next b)
  (if (> a b)
      0
      (+ (term a)
         (sum term (next a) next b))))

;;四角形による積分の近似
(define (integral f a b dx)
  (define (add-dx x) (+ x dx))
  (* (sum f (+ a (/ dx 2.0)) add-dx b) dx))

;;Simpsonの公式
(define (simpson f a b n)
  (define h (/ (- b a) n))
  (define (y k) (f (+ a (* k h))))
  (define (term k)
    (cond ((or (= k 0) (= k n)) (y k))
          ((even? k) (* 2 (y k)))
          (else (* 4 (y k)))))
  (define (next s) (+ s 1))
  (exact->inexact (* (/ h 3) (sum term 0 next n))))

テスト

(add-load-path ".")
(load "q-1-29")

(simpson  cube 0 1 10) ;;gosh> 0.25
(simpson  cube 0 1 2) ;;gosh> 0.25
;;(* (/ (/ 1 2) 3)
;;   (+ (* 4 (cube (/ 1 2)))
;;    (cube 1)))

(integral cube 0 1 0.0001)
;gosh> 0.24999999874993412

(use math.const)

(simpson  sin 0 pi 1000)
;gosh> 2.000000000001081
(integral sin 0 pi 0.0001)
;gosh> 2.0000000008052607

integralは、でおよそ3万の領域に分割して近似しているけど、Simpsonの公式の方では1000の領域にしか分割してないけど、より正確な値が出ている。

問題1.30

先の問題で、sumは再帰的な実装をされていた。反復的にsumを実装しろという問題。

(add-load-path ".")
(load "q-1-29")

;;sumのoverride(iterative)
(define (sum term a next b)
  (define (iter a result)
    (if (> a b) result
        (iter (next a) (+ result (term a)))))
  (iter a 0))

なんか

反復的な実装のほうがシックリくるな。

問題1.31

級数の和を求める手続きはsumで実装。級数の積を求める手続きproductを実装しろという問題。
高階関数で級数を連続的に扱う構造を抽象化していたおかげで、演算子を変更するだけでOK。

反復的に実装

;;級数の積を求める(iterative)
(define (product term a next b)
  (let iter ((a a) (result 1))
    (if (> a b) result
        (iter (next a) (* result (term a))))))

(define (factorial n)
  (product (lambda (n) n) 1
           (lambda (n) (+ n 1)) n))

(define (square n) (* n n))
;;piの近似を求める(引数sは精度,大きいほど正確)
(define (calc-pi s)
  (define (next n) (+ n 1))
  (* 4 (/ (* (product (lambda (n) (* n 2)) 1 next s)
             (product (lambda (n) (* (+ n 1) 2)) 1 next s))
          (square
           (product (lambda (n) (+ (* n 2) 1)) 1 next s)))))

問題文を読んでサクッと反復的に実装したら、小問bでしっかり再帰/反復両方で実装しろとのこと。以下が再帰的な実装。

;;級数の積を求める(factorial)
(define (product term a next b)
  (if (> a b) 1
      (* (term a) (product term (next a) next b))))

テスト

(add-load-path ".")
(load "q-1-31")

(use math.const)

pi
;gosh> 3.141592653589793

(exact->inexact (calc-pi 100))
;gosh> 3.1493784731686008
(exact->inexact (calc-pi 1000))
;gosh> 3.142377365093878
(exact->inexact (calc-pi 10000))
;gosh> 3.1416711865344635

おぉ、PIに収束してる。でも、なんか収束が遅いな。

追記

id:yamanetoshiさんとの勉強会してる最中に

(define (calc-pi s)
  (define (next n) (+ n 1))
  (* 4 (/ (* (product (lambda (n) (* n 2)) 1 next s)
             (product (lambda (n) (* (+ n 1) 2)) 1 next s))
          (square
           (product (lambda (n) (+ (* n 2) 1)) 1 next s)))))

よりもうちょっとましな実装を思いついた。

(define (calc-pi s)
  (define (next n) (+ n 1))
  (/ (/ (square (product (lambda (n) (* n 2)) 1 next (+ s 1))) (+ s 1))
     (square
      (product (lambda (n) (+ (* n 2) 1)) 1 next s))))

こんな感じ。こっちの実装を公式に持ってた方がきれいな公式になると思うけどなー。とりあえずsquareと2つのproductのみで実装できたので満足。(product使わなくても良いならもっと速い実装はできるっぽいけど、本末転倒)

問題1.32

sumとproductをさらに

(accumulate combiner null-value term a next b)

みたいな関数で抽象化せよという問題。

たしかに、演算とその単位元をさらに抽象かすれば級数の和も積も抽象化できる。
以下、実装

;;null-valueが演算の単位元
(define (accumulate combiner null-value term a next b)
  (if (> a b) null-value
      (combiner (term a)
                (accumulate combiner null-value term (next a) next b))))

(define (accumulate combiner null-value term a next b)
  (let iter ((result (term a)) (a (next a)))
    (if (> a b) result
        (iter (combiner result (term a)) (next a)))))

テスト

(add-load-path ".")
(load "q-1-31")
(load "q-1-32")

(define (inc n) (+ n 1))
(define (square n) (* n n))
(define (cube n) (* n n n))
(define (identity n) n)

(define (sum term a next b)
  (accumulate + 0 term a next b))

(sum identity 1 inc 10)
;(accumulate + 0 (lambda (n) n) 1 (lambda (n) (+ n 1)) 10)

(sum cube 1 inc 5)
;gosh> 225
(apply + (map cube '(1 2 3 4 5)))
;gosh> 225

(define (product term a next b)
  (accumulate * 1 term a next b))

(exact->inexact (calc-pi 1000))
;gosh> 3.142377365093878

問題1.33

フィルターを用いたaccumulateの実装。

これだと、nextで「次の素数を求める手続き」というやや難しい手続きを実装しなくても、filterの方で「素数か?」という手続きを、nextで「1つ加算」という手続きを渡すと素数に対する演算を実装できる。filter便利。

(add-load-path ".")
;;prime?手続きをload
(load "q-1-23")

;;フィルター装備なaccumulate
(define (filterd-accumulate combiner null-value filter term a next b)
  (let iter ((result null-value) (a a))
    (cond ((> a b) result)
          ((filter a)
           (iter (combiner result (term a)) (next a)))
          (else (iter result (next a))))))

簡潔簡潔。部品を組み立ててる感じがもの凄い楽しい。

素数の二乗和

(add-load-path ".")
(load "q-1-33")
(use srfi-1)

(define (square n) (* n n))
(define (inc n) (+ n 1))
(define (identity n) n)
(define (sum-of-square a b) (+ a (square b)))

;;区間a,bの素数の二乗の和
(define (sum-of-square-of-primes a b)
  (filterd-accumulate sum-of-square 0 prime? identity a inc b))

(use srfi-1)
(define (square n) (* n n))

(sum-of-square-of-primes 1 10)
;gosh> 87

;;こんな感じで書き換え可
(define (sum-of-square-of-primes a b)
  (apply + (map square (filter prime? (iota (+ (- b a) 1) a)))))

互いに素な整数の積

;;nと互いに素でnより小さい正の整数の積
;;(i < n で,GCD(i, n) = 1なる全てのi)の積
(define (product-of-mutually-prime-of n)
  (filterd-accumulate * 1 (mutually-prime? n) identity 1 inc n))

;;互いに素か判定する関数を返す
(define (mutually-prime? b)
  (lambda (a) (= (GCD a b) 1)))

;;最大公約数
(define (GCD a b)
  (if (= b 0) a
      (gcd b (remainder a b))))

(product-of-mutually-prime-of 10)
;gosh> 189
(filter (mutually-prime? 10) (iota 10 1))
;gosh> (1 3 7 9)
(apply * '(1 3 7 9))

SICP Reading #9 高階関数

ようやっと1.3章までたどりついたー。ペース遅いのは気にしない:-)

高階手続きというプログラミングっぽいお話。

(define (cube x) (* x x x))

のように、三乗の概念を定義しておくことで言語のレベルに縛られずよく使うパターンの抽象を直接扱えるようにするのは強力なプログラミング言語の要請の一つ。

抽象化!さくっ、と実装しちゃうとついパターンに偏ったプログラミングになるけど、手続きとしてまとめておくとあとあと便利だよね。という。(コードも読みやすいし)

高階手続き

手続きに渡せるのが数値だけの場合、抽象の能力がかなり狭くなる。で、手続きに手続きを与える手続き(また、手続きを返してくれたりする手続き)が構成できると、抽象の幅が広がる。
C言語では、関数ポインタとして手続きを渡せたりしますね。うん。(qsortとかおいしい例)

級数の和

;;aからbまでの整数の和を計算
(define (sum-integers a b)
  (if (> a b)
      0
      (+ a (sum-integers (+ a 1) b))))

(sum-integers 1 10)
;;gosh> 55

;;与えられた範囲の整数の三乗の和を計算
(define (cube x) (* x x x))
(define (sum-cubes a b)
  (if (> a b)
      0
      (+ (cube a) (sum-cubes (+ a 1) b))))

(sum-cubes 1 10)
;;gosh> 3025

;;級数項の並びの和
(define (pi-sum a b)
  (if (> a b)
      0
      (+ (/ 1.0 (* a (+ a 2))) (pi-sum (+ a 4) b))))

(pi-sum 1 10)
;;gosh> 0.372005772005772
;;(/ pi 8.0)に収束

のように、級数の和として同じようなパターンを持つ計算を、共通部を抜き出し以下のように抽象化。

(define (<name> a b)
  (if (> a b)
      0
      (+ (<term> a)
         (<name> (<next> a) b))))

termが、各項を求める手続き。nextが次項を求めるる手続き
これは、級数の総和を求める数学の

をプログラミングで表現した感じ。

この抽象をsumと定義して、先ほどの特殊例を高階関数で書き直す。

;;級数の和を抽象化
(define (sum term a next b)
  (if (> a b)
      0
      (+ (term a)
         (sum term (next a) next b))))

;;次項:一つ進める
(define (inc n) (+ n 1))

;;高階関数で書き直し
;;三乗の和
(define (sum-cube a b)
  (sum cube a inc b))

;;項=値
(define (identify x) x)
;;和
(define (sum-integers a b)
  (sum identify a inc b))

;;級数の並びの和
(define (pi-sum a b)
  (define (pi-term x)
    (/ 1.0 (* x (+ x 2))))
  (define (pi-next x)
    (+ x 4))
  (sum pi-term a pi-next b))

;;piの近似
(* 8 (pi-sum 1 10000))

積分近似

級数の和を定義いたことで、積分の近似式の部品としても使用できる。
以下な積分の近似式をschemeで、

;;積分の近似も定義できる!
(define (integral f a b dx)
  (define (add-dx x) (+ x dx))
  (* (sum f (+ a (/ dx 2.0)) add-dx b) dx))

(integral cube 0 1 0.001)
;;gosh> 0.249999875000001
(integral cube 0 1 0.0001)
;;gosh> 0.24999999874993412

単純な近似式を用いてるとはいえ、積分がこんなに簡単なコードで表現できるのは嬉しい。
高階関数すごい。

問題は、、勉強会の後でまとめる!

SICP Reading #8 素数性テスト (Fermat,Miller-Rabin)

素数性テストな問題の連発。

問題1.21

1でない最小の除数を求める手続きを実装。

;;除数の探索
(define (square n) (* n 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 (+ test-divisor 1)))))

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

(map smallest-divisor  '(199 1999 19999))
;;gosh> (199 1999 7)

結果から、199,1999は素数で、19999は7を約数に持つことが分かる。

問題1.22

素数性テストにかかる時間を実測するとい問題。
SICPでは、現在時刻を取得するruntime手続きなるものが使用されていた。

runtimeは?

gaucheでのruntime。 n-oohiraのSICP日記さんでちょっとまとめられてたり。

一番なにやってるか分かりやすい sys-gettimeofday からのマイクロ秒変換での方法を使う事に。情報感謝。

;;gauche用runtime(現在時刻)関数
;;sys-gettimeofdayは(values 秒数 マイクロ秒)を返す
(define (runtime)
  (use srfi-11)
  (let-values (((a b) (sys-gettimeofday)))
              (+ (* a 1000000) b)))
;;秒を1000000倍してマイクロ基準の秒数に

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

;;n以上の、k個の素数リストを返す
(define (search-for-primes n k)
  (if (even? n) (set! n (+ n 1)))
  (let iter ((n n) (k k) (lst '()))
    (cond ((= k 0) lst)
          ((prime? n) (iter (+ n 2) (- k 1) (append lst (list n))))
          (else (iter (+ n 2) k lst)))))

;;O(sqrt(n))な実装
(define (prime? n)
  (if (< n 2) #f
      (let iter ((n n) (i 2) (l (floor->exact (sqrt n))))
        (cond ((> i l) #t)
              ((= 0 (remainder n i)) #f)
              (else (iter n (+ i 1) l))))))

テスト

(add-load-path ".")
(load "q-1-22")

;;一応、テスト
((lambda (micro)
   (let1 s (runtime)
         (sys-nanosleep (* 1000 micro))
         (- (runtime) s))) 100000)
;;だいたい正しいね

(map (lambda (n) (for-each timed-prime-test (search-for-primes n 3)))
     '(1000 10000 100000 1000000 10000000 100000000 1000000000))
;;gosh>
1009 *** 17
1013 *** 10
1019 *** 10
10007 *** 23
10009 *** 24
10037 *** 24
100003 *** 66
100019 *** 66
100043 *** 65
1000003 *** 226
1000033 *** 198
1000037 *** 202
10000019 *** 616
10000079 *** 614
10000103 *** 648
100000007 *** 1939
100000037 *** 1939
100000039 *** 1791
1000000007 *** 11712
1000000009 *** 12113
1000000021 *** 8847
;;だいたい(sqrt 10)倍ずつ #=> 3.1622776601683795

問題1.23

prime?手続きをもうちょっと賢くした実装。
全てのn以下な整数について「割り切れないか」をチェックするのではなく、全てのn以下な奇数について調べれば良いという実装。(nが偶数の場合は、必ず2で割り切れる)

(add-load-path ".")
(load "q-1-22")

;;O(sqrt(n))な実装
;;もうちょっと細かくいうとsqrt(n)/2ステップぐらい
(define (prime? n)
  (cond ((= n 2) #t)
        ((or (< n 2) (= 0 (remainder n 2))) #f)
        (else
         (let iter ((n n) (i 3) (l (floor->exact (sqrt n))))
           (cond ((> i l) #t)
                 ((= 0 (remainder n i)) #f)
                 (else (iter n (+ i 2) l)))))))

テスト

(add-load-path ".")
(load "q-1-23")

(map (lambda (n) (for-each timed-prime-test (search-for-primes n 3)))
     '(1000 10000 100000 1000000 10000000 100000000 1000000000))
;;だいたい時間が(q-1-22)の半分位に
1009 *** 7
1013 *** 5
1019 *** 4
10007 *** 10
10009 *** 10
10037 *** 10
100003 *** 28
100019 *** 27
100043 *** 28
1000003 *** 83
1000033 *** 80
1000037 *** 83
10000019 *** 247
10000079 *** 314
10000103 *** 271
100000007 *** 851
100000037 *** 848
100000039 *** 897
1000000007 *** 7138
1000000009 *** 4428
1000000021 *** 4269

問題1.24

確率的手法Fermatテストを利用した fast-prime? での時間測定。

(add-load-path ".")
(load "q-1-22")

;;fast-prime?で上書き
(define (start-prime-test n start-time)
  (if (fast-prime? n 3) ;;3はてけとー
      (report-prime (- (runtime) start-time))))

;;素数の確率的判定
;;この判定を騙す素数でない自然数をCarmichael数という
;;'(561 1105 1729 2465 2821 6601)

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

;;メルセンヌツイスターを使用(いつか自分で実装)
(use srfi-27)
(define (fermat-test n)
  (if (< n 2) #f
      (let1 a (+ (random-integer (- n 1)) 1)
            (= (expmod a n n) a))))

;;高速な素数判定(確率的)
(define (fast-prime? n times)
  (cond ((= times 0) #t)
        ((fermat-test n) (fast-prime? n (- times 1)))
        (else #f)))

テスト

(add-load-path ".")
(load "q-1-24")

(map (lambda (n) (for-each timed-prime-test (search-for-primes n 3)))
     '(1000 10000 100000 1000000 10000000 100000000 1000000000))
;;めちゃんこ早い
1009 *** 49
1013 *** 44
1019 *** 31
10007 *** 36
10009 *** 36
10037 *** 36
100003 *** 81
100019 *** 111
100043 *** 93
1000003 *** 126
1000033 *** 120
1000037 *** 128
10000019 *** 154
10000079 *** 152
10000103 *** 170
100000007 *** 172
100000037 *** 178
100000039 *** 196
1000000007 *** 188
1000000009 *** 161
1000000021 *** 161

1-22,1-23の実装との速度比較

(use srfi-1)
(map (lambda (n) (for-each timed-prime-test (search-for-primes n 1)))
  (map (lambda (n) (expt 2 n))(iota 40 1)))

こんな感じでそれぞれ2^n (n 1 ~ 40) < な最小の素数に対してそれぞれテスト。
かかった時間をgnuplotでグラフに。

prime?とその改良版,及びfast-prime?の時間比較

prime?とその改良版,及びfast-prime?の時間比較

改良する前のprime?(q22)と、奇数のみをチェックする改良したソレ(q23)のグラフが奇麗に1/2になっている。そしてFermatテストを利用した確率的なfast-prime?(q24)はそれらに比べると時間はほぼ無視できる位小さい。これがO(sqrt(n))とO(log(n))の違い。

上では全く見えていないfast-prime?のみをグラフにするとこんな感じ。

fast-prime?の増加具合

fast-prime?の増加具合

うーん、指数的なステップの増加に対して直線的な増加。log。

問題1-25

expmodを、再帰的ではなく超簡単に実装した場合、計算が遅くなるのはなぜか?という問題。
実装は以下な感じ

(add-load-path ".")
(load "q-1-24")

;;めちゃくちゃ遅い!!
(define (fast-expt n exp)
  (cond ((= exp 0) 1)
        ((even? exp)
         (square (fast-expt n (/ exp 2))))
        (else (* n (fast-expt n (- exp 1))))))

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

このexpmodと再帰的なexpmodの一番の違いは、扱う値の大きさ。
引数としてはおなじ(expmod 10 97 97)とかでも、再帰的なexpmodでは再帰のたびにremainderを適用しているので、評価する値が巨大にならない。

だけど、この実装では(expmod 10 97 97)は実際に10^97をもとめているので、扱う値が巨大になって時間がかかってしまう。

つまり、再帰で部分的にremainderで剰余を求めていかないと、引数がどんどん巨大に。#=> 処理が重くなる><

問題1.26

指数が偶数の時、square演算を用いて呼び出しを半分にしているけど、square演算ではなくて単純に*で2乗してはだめか?という問題。
実装は以下

(add-load-path ".")
(load "q-1-24")

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

;;ものすごく遅い
(map (lambda (n) (for-each timed-prime-test (search-for-primes n 3)))
     '(1000 10000 100000))

これって、偶数の時は(/ exp 2)で半分にした呼び出しを2回求めているので、結局呼び出し回数は1/2*2=1となってしまう。なので、オーダはO(log(n)^2)となりO(n)

じつはこの問題以前に、「こっちの実装の方が簡単じゃん」と思って、いかな実装でq-1-24を渡航としてgaucheが停止してしまったり:-)

;;簡単に実装(exptも組み込みを使用)
(define (expmod base exp m)
  (remainder (expt base exp) m))

しかも、どうやら単純に対象が自然数な累乗演算は、組み込みのexptよりも自前で実装したfast-expの方が速い感じ。 (exptは負数やfloatにも対応してるけど)

問題1.27

Fermatテストをつねに欺くCarmichael数の存在を証明する問題。
Carmichael数の性質はよく分かってないらしい。

全てのn < cとなるnに対し、(= (expmod n c c) n)が常に真となるかのチェック。

(add-load-path ".")
(load "q-1-24")

;;a < nで全て pass するか?
;;するなら、素数かcarmichael数
(define (test carmichael)
  (let iter ((c carmichael) (i 1))
    (cond ((= i c) #t)
          ((= (expmod i c c) i) (iter c (+ i 1)))
          (else #f))))

;;Carmichae数?
(define Carmichaels '(561 1105 1729 2465 2821 6601))
(map test Carmichaels)
;;gosh> (#t #t #t #t #t)

;;素数?
(map prime? Carmichaels)
;;gosh> (#f #f #f #f #f )

問題1.28

Carmichaels数にだまされないテスト。Fermatテストが疑い深くなったver.

Miller-Rabinテストといって、expmodの再帰の過程で2乗を行なう段階で毎回チェックを施す。
チェックは(expmod a c c)の時、1でもc-1でもなく,かつn^2のcを法とした剰余演算が1を満たす値が現れるか。とうもの。

チェックのルールは証明を見ないと良くわからないけど、ここはスルー:-) (また今度!)
説明によると、cが素数でない場合は上のルールを満たすa(< c)が少なくとも半分以上は存在するということ。

以下、実装

;;素数の確率的判定
;;だまされないFermatテスト
;;'(561 1105 1729 2465 2821 6601)

;;メルセンヌツイスターを使用(いつか自分で実装)
(define (miller-rabin-test n)
  (if (< n 2) #f
      (let1 a (+ (random-integer (- n 1)) 1)
            (= (expmod a n n) a))))

(use srfi-27)
(define (square n) (* n n))
(define (expmod base exp m)
  (cond ((= 0 exp) 1)
        ((even? exp)
         (let* ((ret (expmod base (/ exp 2) m))
                (val (remainder (square ret) m)))
           (if (and (not (= ret 1)) (not (= ret (- m 1))) (= val 1))
               0 val)))
        (else
         (remainder (* base (expmod base (- exp 1) m)) m))))

;;高速な素数判定(確率的)
(define (fast-prime? n times)
  (cond ((= times 0) #t)
        ((miller-rabin-test n) (fast-prime? n (- times 1)))
        (else #f)))

ただ、Miller-RabinテストもCarmichaels数にだまされないと言ってもそれも確率的な話
全ての数に対してチェックが通れば本物の素数だけど、それだと確率的手法のメリットがない。

テスト

(add-load-path ".")
(load "q-1-23")
(load "q-1-28")

(define Carmichaels '(561 1105 1729 2465 2821 6601))
;;判定の可能性 <= 1/2^(テスト回数 - 1) (1/1024)
(map (lambda (n) (fast-prime? n 10)) Carmichaels)
;;gosh> (#f #f #f #f #f #f)

(use srfi-1)
;;本物の素数列
(define primes (filter prime? (iota 1000 1)))
(length primes)
;;gosh> 168

;;miller-rabinテストと比較
(equal? primes (filter (lambda (n) (fast-prime? n 10)) (iota 1000 1)))
;;gosh> #t
;;だまされてないね!

テストの思考回数を10でためしてみると、1000以下の素数判定はパス。
この場合、素数でない数がテストをパスする確率は、1/2^10となる(かな?)。

SICP reading #6 [1.2.4] べき乗

wp-emacsという最強の助っ人が現れたので、ブログの投稿作業自体が楽しくなってきた。
wp-emacsについては別なエントリにあげるとして、とりあえずsicpの問題とくよ!

q1.16

sicp勉強会で、問題を読んで即実装に取りかかって小一時間ハマってしまった問題。
逐次平方で対数的ステップ数の反復的な羃乗プロセスを実装しろとのこと。(長!

逐次平方で対数的ステップ数っていうのは、例題であったように2の羃乗の指数の場合指数を1/2した結果を2乗すれば乗算のステップが半分になるという上手い方法。
普通に再帰で実装すれば羃乗は以下のように

(define (expt b n)
  (if (= n 0)
      1
      (* b (expt b (- n 1)))))

これは、ステップ/スペース数ともにO(n)な実装。

逐次平方を用いた再帰的羃乗プロセスは

(define (square n) (* n n))
(define (fast-expt-recr b n)
  (cond ((= n 0) 1)
        ((even? n) (square (fast-expt b (/ n 2))))
        (else (* b (fast-expt b (- n 1))))))

これだと、指数が2の乗数の場合に最小の計算量が見込め、O(log(n))となる。

ただ、再帰的な実装のためスペース数もO(log(n))となる。
で、q1.16ではそれを反復的に実装しろとのこと。

ハマった

問題を読んで感覚的に実装すると、以下のようなコードに

(define (fast-expt-iter b n)
  (define (iter b n a)
    (cond ((= n 0) a)
          ((even? n) (iter b (/ n 2) (square (* b a))))
          (else (iter b (- n 1) (* b a)))))
  (iter b n 1))

指数nが偶数の時引数としてのnを1/2にするのは正しいけれど底bと結果aが間違ってる。
結果を一回きり2乗して指数を半分にしているため、実際の羃乗よりはるかに小さい値となる。

で、逐次平方で反復的な羃乗のプロセスの生成具合を書き下してみると、指数を半分にしたときに底を2倍すればよいことに気付いた(簡単だよね!最初っから書いておけば。。)

以下、正しい実装。

(define (fast-expt-iter b n)
  (define (iter b n a)
    (cond ((= n 0) a) ;;偶数なら底を2乗 and 指数を/2
          ((even? n) (iter (square b) (/ n 2) a))
          (else (iter b (- n 1) (* b a)))))
  (iter b n 1))

以下、動作確認

(use srfi-1)
(map (lambda (n) (print 2 "^" n "=" (fast-expt-iter 2 n))) (iota 11 0))
;;実行結果
gosh> 2^0=1
2^1=2
2^2=4
2^3=8
2^4=16
2^5=32
2^6=64
2^7=128
2^8=256
2^9=512
2^10=1024
(#<undef> #<undef> #<undef> #<undef> #<undef> #<undef> #<undef> #<undef> #<undef> #<unde\f> #<undef>)
gosh>

上では(というか今ままで全部)、実装と同一ファイルに直接テストコードを埋め込んで確認している。
勉強会を一緒にしてるid:yamanetoshiさんからgauche.testの使い方を教えてもらったので、使ってみた。

gauche.testでのテスト

(use gauche.test)
(add-load-path ".")
(load "q-1-16")

(test-start "expt")
(test-section "expt")
(test* "2^0"
       1
       (fast-expt-iter 2 0))
(test* "2^1"
       2
       (fast-expt-iter 2 1))
(test* "2^2"
       4
       (fast-expt-iter 2 2))
(test* "2^3"
       8
       (fast-expt-iter 2 3))
(test* "2^4"
       16
       (fast-expt-iter 2 4))
(test* "2^5"
       32
       (fast-expt-iter 2 5))
(test* "2^6"
       64
       (fast-expt-iter 2 6))
(test* "2^7"
       128
       (fast-expt-iter 2 7))
(test* "2^8"
       256
       (fast-expt-iter 2 8))
(test* "2^9"
       512
       (fast-expt-iter 2 9))
(test* "2^10"
       1024
       (fast-expt-iter 2 10))

(test-end)
で、実行
SHINYA% gosh q-1-16-test.scm
Testing expt ...
<expt>-------------------------------------------------------------------------
test 2^0, expects 1 ==> ok
test 2^1, expects 2 ==> ok
test 2^2, expects 4 ==> ok
test 2^3, expects 8 ==> ok
test 2^4, expects 16 ==> ok
test 2^5, expects 32 ==> ok
test 2^6, expects 64 ==> ok
test 2^7, expects 128 ==> ok
test 2^8, expects 256 ==> ok
test 2^9, expects 512 ==> ok
test 2^10, expects 1024 ==> ok
passed.

うぉぉ、楽しい!passed!

q1.17

羃乗の次は乗算を対数ステップ数で実装しろという問題。
羃乗の場合は指数を1/2して平方をとっていたから逐次平方、ということはこの場合は逐次2倍?

以下、線形/対数オーダーなそれぞれの実装

;;線形オーダーな実装
(define (* a b)
  (if (= b 0)
      0
      (+ a (* a (- b 1)))))

(define (double n) (+ n n))
(define (halve  n) (/ n 2))
;;対数オーダーな実装
(define (fast-* a b)
  (cond ((= b 0) 0)
        ((even? b) (fast-* (double a) (halve b)))
        (else (+ a (fast-* a (- b 1))))))

gauche profiler

さっきの問題でもスルーしてたけど、実際に対数的かどうかgaucheのprofilerを使って確認してみた。

;;profiler便利
(profiler-reset)
(profiler-start)
(fast-* 2 0)
(profiler-show)
;;以下、実行結果(抜粋)
gosh> Profiler statistics (total 0 samples, 0.0 seconds)
                                                    num    time/    total
Name                                                calls  call(ms) samples
---------------------------------------------------+------+-------+-----------
fast-*                                                  12  0.0000     0(  0%)
even?                                                   11  0.0000     0(  0%)
double                                                  10  0.0000     0(  0%)
halve                                                   10  0.0000     0(  0%)

実際にはprofiler-showでの結果はgaucheの内部っぽい関数(compiled-code-emit1i!とか)も拾ってくるけど、ここでは自分で実装した関数の結果のみ表示。

そのまま実装すると1024回の加算が必要だけどfast-*の呼び出しは12回。これは、最初の1回目の呼び出しとさらにbが0になったとき0を返す分の呼び出しが1回あるので10+2=12。
bが1になった時そのままaを返せば1回分呼び出しが減るけど、(fast-* 2 0)のとき0を返してほしかったのでこう実装しました。

q1.18

q1.16,q1,17の復習問題。ここまで復習すると逐次平方的な考えはバッチリ。

以下、実装

(define (double n) (+ n n))
(define (halve  n) (/ n 2))

(define (fast-* a b)
 (let i ((a a) (b b) (r 0))
   (cond ((= b 0) r)
         ((even? b) (i (double a) (halve b) r))
         (else (i a (- b 1) (+ r a))))))

試験ではmapの手続きにさらにmapを仕込んで9×9の掛け算を生成。

試験(test & profile)

;;試験
(let ((1to9 '(1 2 3 4 5 6 7 8 9)))
  (map (lambda (a)
         (map (lambda (b) (if (= (fast-* a b) (* a b))
                              (* a b) 'error))
              1to9)) 1to9))
;;even?の数で把握
(profiler-reset)
(profiler-start)
(fast-* 10 2048)
(profiler-show)
;;以下、実行結果(test)
gosh> ((1 2 3 4 5 6 7 8 9) (2 4 6 8 10 12 14 16 18) (3 6 9 12 15 18 21 24 27) (4 8 12 16\
 20 24 28 32 36) (5 10 15 20 25 30 35 40 45) (6 12 18 24 30 36 42 48 54) (7 14 21 28 35 \
42 49 56 63) (8 16 24 32 40 48 56 64 72) (9 18 27 36 45 54 63 72 81))
;;profile
gosh> Profiler statistics (total 0 samples, 0.0 seconds)
                                                    num    time/    total
Name                                                calls  call(ms) samples
---------------------------------------------------+------+-------+-----------
even?                                                   12  0.0000     0(  0%)
double                                                  11  0.0000     0(  0%)
halve                                                   11  0.0000     0(  0%)

うーん、素晴らしい。emacs & scheme最高!

q1.19

フィボナッチ数を求めるプロセスを、変換Tの特殊な場合としてTを

みたいに定義している。
p=0,q=1の場合は確かに

となるね。Tに名前はついてるのかな?(誰か教えて

この問題ではT_p_qによる2回変換をT_p’_q’で表せと問うている。実際に計算を行なってp’,q’をp,qで表せば良いので、結構簡単。

手書きで計算したので、画像up。(文字が汚いのは仕様です。

手計算で

手計算で

計算結果をそのままコードに落とし込むと

(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))
                   (* q (+ q (* 2 p)))
                   (/ count 2)))
        (else (fib-iter (+ (* b q) (* a q) (* a p))
                        (+ (* b p) (* a q))
                        p
                        q
                        (- count 1)))))

(map fib '(1 2 3 4 5 6 7 8 9 10))

(profiler-reset)
(profiler-start)
(fib 2048)
(profiler-show)

終了!

まとめ?

ここ最近、emacsの操作方法が色々分かってきたので作業効率が異常なまでに良くなってきている。
wp-emacsとかtwittering-modeとか、emacsは便利なプラグインもさることながら基本操作が半端無く充実しているなー。いずれエントリにもまとめてみたいかも。オライリーのemacs本読もっと。

入門 GNU Emacs 第3版
計算機プログラムの構造と解釈

SICP reading #5 [1.2.3] 増加の程度

前々回の勉強会で宿題となった問題1.14の解答。(時間空け過ぎ!

q1.14

(cc 11 5)が使うステップ/スペース数を求める問題。
ステップ数は再帰の数、スペース数は再帰呼び出しのツリーが最大に展開された時の数。
とりあえず、ccの中身

;;両替の計算
(define (first-denomination kinds-of-coins)
  (cond ((= kinds-of-coins 1) 1)
        ((= kinds-of-coins 2) 5)
        ((= kinds-of-coins 3) 10)
        ((= kinds-of-coins 4) 25)
        ((= kinds-of-coins 5) 50)))

;;両替手続き
(define (cc amount kinds-of-coins)
  (cond ((= amount 0) 1)
        ((or (< amount 0) (= kinds-of-coins 0)) 0)
        (else (+ (cc
                  (- amount
                     (first-denomination kinds-of-coins))
                  kinds-of-coins)
                 (cc amount (- kinds-of-coins 1))))))
(cc 11 5) ;;問題

とりあえず、(cc 11 5)が再帰呼び出し行なってく手続きを手書きでツリーに。

(cc 11 5)の生成プロセス

(cc 11 5)の生成プロセス


見にくいとかのツッコミはおいといて、書いていく事でこのツリーの特徴がいくつか分かった。

  • ノードの左にはコインの数だけ伸びていく
  • ノードの右にはamountの数だけ(に比例てして)伸びていく
  • 全て1 centで両替するプロセスが一番深い枝に
  • (cc 11 5)は55ノード(呼び出し)

amountがマイナス or コインが0の葉は、amountが増えたところでそれ以上呼び出しは行なわれない。(死んだ枝)。
ここで、amountがマイナスの葉は、最初の呼び出しのamountが有る程度増えると伸びていく。
全ての節は必ず2つの子をもつ(葉以外は子が2つ)ので、安易にステップ数はO(n^2)かと思ったけど、死んだ枝の割合も計算してオーダーを求めないと答えが求まらないっぽい。
あ、というか生きた枝(amountが1増えると子が2つ増える葉を含むパス)っていうのは、そのまんまccで求める両替できるパターンの数になることに気付く。(= amount 0)な枝。
たとえば、(cc 11 5)の場合は答えは4なので、(cc 12 5)が呼び出された場合は、この4つの葉に子が2つ増える。

単純に考えてみると、ccの返す数(両替のパターン数)がnに比例の場合cc全体のステップ数は(ccの返す値)^2でO(n^2)かと思ったけど、問題はamountがマイナスな枝(> amount 0)。
こいつは、amountがある程度増えるとまた生きた枝として伸びていく(しかも、コインが1以上)。

強行

色々ごちゃごちゃしてきたので、2分木でありamountが十分大きいと死んだ枝の割合とかは無視できると思うということでステップ数はO(n^2)と予測。

で、実際にccの呼び出し回数(ツリーのノードの数)を計算する関数を以下のように定義して

(define (ccount amount koc)
  (if (or (= koc 0) (not (> amount 0)))
      1 ;;葉
      (+ 1 ;;自分自身(節)
         (ccount amount (- koc 1)) ;;左右の枝
         (ccount (- amount (first-denomination koc)) koc)))))

gnuplotで呼び出し回数のグラフを作成した。

(use srfi-1) ;;フォーマットして出力
(map (lambda (n) (print n " " (ccount n 5))) (iota 10 50 50))

出力したグラフは、以下。

グラフ - ccのステップ数の増加

グラフ - ccのステップ数の増加


おお。すばらしく指数的増加。
この形からn^2でない可能性もあるけど(n^x, x > 1というのは言える?)、二分木ということでO(n^2)と結論。

スペースは単純で、ツリーの最大に展開されるのはamountはそのままでコインが1に(ルートからの左端の葉の親)から、右端の子までのパスなので、単純にamount + 5(コインの数)に比例。ということでO(n)

q1.15

q1.14に比べると、スイーツな問題。
xが十分小さい時、sinは以下のように近似できる(ぽい)

この近似式は、、sinのマクローリン展開からxが小さいという条件でごにょごにょした式かな?そういことにしておこう :-)

とりあえず、問題にのってる実装。

(define (cube n) (* n n n))
(define (p x) (- (* 3 x) (* 4 (cube x))))

;;angleを3で割っていき、0.1より小さくなった時点が
;;再帰の枝となる(0.1を返す)
(define (sine angle)
  (if (not (> (abs angle) 0.1))
      angle
      (p (sine (/ angle 3.0)))))

単純な関数作るのも、schemeなら倍楽しい。(謎
なんで関数名がsineなのかちょっと分からなかったけど、組み込みでsinが定義されているからでした。

実際にこの近似式を実行してみると、

;;組み込み
(sin 12.15)  ;;gosh> -0.4044438228491401
;;近似の利用
(sine 12.15) ;;gosh> -0.39980345741334

おー、掛け算/乗算/引き算の組み合わせでsinがちゃんと近似されてる。すごい。

問題は、(sine 12.15)がプロセスpを生成する回数とsineのスペース/ステップ数

12.15で予備だした場合、3で割った値で再帰的に呼び出すので、、
;;12.15->4.05->1.3499…->0.4499…->0.15->0.0499..(枝)
といった具合のツリーになる。pの生成は再帰の1回目からなのでこの場合のpの生成回数は5回。

sineの作るツリーは枝が1方向に伸びていくだけなので、ステップ数とスペース数は等しい。
オーダーは、3で割って再帰的に呼び出すということでO(log(3,n))というのが大体分かる。

TeXではめんどくさいので手書きで式らしきものを書いた、一応。

引数がaの場合、呼び出し回数nは

オーダーの導出

オーダーの導出


という感じで導出できる。(字が汚いのは、画像の乱れではありません。仕様です。)
オーダーはaが十分大きくなっていた場合に増加の具合を表すもので、この式の定数log(3,10)は省略してO(log(3,n))となる。 参考:アルゴリズムと計算量

 
Better Tag Cloud