SICP 読書ノート#38 - 3.5.3 ストリームパラダイムの開発 (pp.198-205)
前半はある一定の値に収束していく無限ストリーム(数列)の収束を加速させる手法や、対のストリームについて。数学的なトピックが中心だと思うので、バッサリ飛ばします。早く4章に行きたいのも大きいけど…
後半の「信号としてのストリーム」の節はストリームの工学的応用で気になるので解いて行きたいと思います。
(2015/09/06追記) その後「対の無限ストリーム」が§4.4のquery systemの内部実装で登場したので戻ってきて解きました。結果的にストリームの理解が深まったのでよかったかな
反復をストリームとして形式化する
まずはsqrt-streamを写経して動かす。わりとすぐに収束しているのがわかる。
racket@> (map (lambda (k) (stream-ref (sqrt-stream 2) k)) (enumerate-interval 0 9)) '(1.0 1.5 1.4166666666666665 1.4142156862745097 1.4142135623746899 1.414213562373095 1.414213562373095 1.414213562373095 1.414213562373095 1.414213562373095)
一方、pi-streamはすぐには収束しない。
racket@> (map (lambda (k) (stream-ref pi-stream k)) (enumerate-interval 0 9)) '(4.0 2.666666666666667 3.466666666666667 2.8952380952380956 3.3396825396825403 2.9760461760461765 3.2837384837384844 3.017071817071818 3.2523659347188767 3.0418396189294032)
しかし、オイラーによる方法やタブローを使って収束を加速させることができる。
racket@> (map (lambda (k) (stream-ref (euler-transform pi-stream) k)) (enumerate-interval 0 9)) '(3.166666666666667 3.1333333333333337 3.1452380952380956 3.13968253968254 3.1427128427128435 3.1408813408813416 3.142071817071818 3.1412548236077655 3.1418396189294033 3.141406718496503) racket@> (map (lambda (k) (stream-ref (accelerated-sequence euler-transform pi-stream) k)) (enumerate-interval 0 9)) '(4.0 3.166666666666667 3.142105263157895 3.141599357319005 3.1415927140337785 3.1415926539752927 3.1415926535911765 3.141592653589778 3.1415926535897953 3.141592653589795)
数学的な証明は置いておいて、感じはつかめたので練習問題へ。
問題 3.63
2つのsqrt-streamをguessの推移がわかるように実装する。
(define (sqrt-stream-1 x) (define guesses (cons-stream 1.0 (stream-map (lambda (guess) (display (format "guess=~a ~%" guess)) (sqrt-improve guess x)) guesses))) guesses) (define (sqrt-stream-2 x) (cons-stream 1.0 (stream-map (lambda (guess) (display (format "guess=~a ~%" guess)) (sqrt-improve guess x)) (sqrt-stream-2 x))))
テスト。後者の例はこれまで結果が蓄積されないのでsqrt-streamが繰り返し呼ばれる。
racket@> (stream-ref (sqrt-stream-1 2) 5) guess=1.0 guess=1.5 guess=1.4166666666666665 guess=1.4142156862745097 guess=1.4142135623746899 1.414213562373095 racket@> (stream-ref (sqrt-stream-2 2) 5) guess=1.0 guess=1.0 guess=1.5 guess=1.0 guess=1.5 guess=1.4166666666666665 guess=1.0 guess=1.5 guess=1.4166666666666665 guess=1.4142156862745097 guess=1.0 guess=1.5 guess=1.4166666666666665 guess=1.4142156862745097 guess=1.4142135623746899 1.414213562373095
ただし、メモ化が行われない場合は、前者も後者と同じ結果となる。
問題 3.64
\( s_{i} \) と \( s_{i+1} \) の差分がtoleranceを下回った時点で \( s_{i} \) を返すような手続きを実装すればよい。
(define (sqrt x tolerance) (stream-limit (sqrt-stream x) tolerance)) (define (stream-limit s tolerance) (define (iter s count) (let ((s0 (stream-ref s 0)) (s1 (stream-ref s 1))) (if (< (abs (- s0 s1)) tolerance) (cons s1 count) (iter (stream-cdr s) (+ count 1))))) (iter s 0))
テスト。sqrtはわりと早く収束しているのがわかる。
racket@> (sqrt 2 0.01) '(1.4142156862745097 . 2) racket@> (sqrt 2 0.00001) '(1.4142135623746899 . 3)
問題 3.65
pi-streamに倣って素直に実装するだけ。
(define (ln2-sum n) (cons-stream (/ 1.0 n) (stream-map - (ln2-sum (+ n 1))))) (define ln2-stream (partial-sums (ln2-sum 1)))
テスト。タブローによる加速が速い！
racket@> (log 2) 0.6931471805599453 racket@> (map (lambda (n) (stream-ref ln2-stream n)) (enumerate-interval 0 10)) '(1.0 0.5 0.8333333333333333 0.5833333333333333 0.7833333333333332 0.6166666666666666 0.7595238095238095 0.6345238095238095 0.7456349206349207 0.6456349206349207 0.7365440115440116) racket@> (map (lambda (n) (stream-ref (euler-transform ln2-stream) n)) (enumerate-interval 0 10)) '(0.7 0.6904761904761905 0.6944444444444444 0.6924242424242424 0.6935897435897436 0.6928571428571428 0.6933473389355742 0.6930033416875522 0.6932539682539683 0.6930657506744464 0.6932106782106783) racket@> (map (lambda (n) (stream-ref (accelerated-sequence euler-transform ln2-stream) n)) (enumerate-interval 0 10)) '(1.0 0.7 0.6932773109243697 0.6931488693329254 0.6931471960735491 0.6931471806635636 0.6931471805604039 0.6931471805599445 0.6931471805599427 0.6931471805599454 +nan.0)
対の無限ストリーム
interleaveの発想がおもしろい。こんなの思いつかないわ。
問題 3.66
まずはpairsを動かしてみる。integersは§3.5.2のものを流用。
racket@> (map (lambda (k) (let ((p (stream-ref (pairs integers integers) k))) (cons k (list p)))) (enumerate-interval 0 25)) '((0 (1 1)) (1 (1 2)) (2 (2 2)) (3 (1 3)) (4 (2 3)) (5 (1 4)) (6 (3 3)) (7 (1 5)) (8 (2 4)) (9 (1 6)) (10 (3 4)) (11 (1 7)) (12 (2 5)) (13 (1 8)) (14 (4 4)) (15 (1 9)) (16 (2 6)) (17 (1 10)) (18 (3 5)) (19 (1 11)) (20 (2 7)) (21 (1 12)) (22 (4 5)) (23 (1 13)) (24 (2 8)) (25 (1 14)))
\( \left( S_{i}, T_{j} \right) \) と \( k \) の関係について \( i \) を固定して見た時の\( k \) と \( T_{j} \) に関係性がありそう。
- \( k \) と \( \left( S_{0}, T_{j} \right) \)
racket@> (filter (lambda (x) (= (caadr x) 1)) (map (lambda (k) (let ((p (stream-ref (pairs integers integers) k))) (cons k (list p)))) (enumerate-interval 0 10))) '((0 (1 1)) (1 (1 2)) (3 (1 3)) (5 (1 4)) (7 (1 5)) (9 (1 6)))
- \( k \) と \( \left( S_{1}, T_{j} \right) \)
racket@> (filter (lambda (x) (= (caadr x) 2)) (map (lambda (k) (let ((p (stream-ref (pairs integers integers) k))) (cons k (list p)))) (enumerate-interval 0 20))) '((2 (2 2)) (4 (2 3)) (8 (2 4)) (12 (2 5)) (16 (2 6)) (20 (2 7)))
- \( k \) と \( \left( S_{2}, T_{j} \right) \)
racket@> (filter (lambda (x) (= (caadr x) 3)) (map (lambda (k) (let ((p (stream-ref (pairs integers integers) k))) (cons k (list p)))) (enumerate-interval 0 30))) '((6 (3 3)) (10 (3 4)) (18 (3 5)) (26 (3 6)))
- \( k \) と \( \left( S_{3}, T_{j} \right) \)
racket@> (filter (lambda (x) (= (caadr x) 4)) (map (lambda (k) (let ((p (stream-ref (pairs integers integers) k))) (cons k (list p)))) (enumerate-interval 0 60))) '((14 (4 4)) (22 (4 5)) (38 (4 6)) (54 (4 7)))
\( \left( S_{i}, T_{j} \right) \) の時の \( k \) を \( k_{ij} \) とするとこんな感じ。
- \( i = 0 \) かつ \( j = 0 \) の場合、\( k_{ij} = 0 \)
- \( i = j \) の場合、\( k_{ij} = k_{(i-1)(j-1)} + 2^{i} \)
- \( i = j + 1 \) の場合、\( k_{ij} = k_{i(j-1)} + 2^{i} \)
- \( i \geq j + 2 \) の場合、\( k_{ij} = k_{i(j-1)} + 2^{i+1} \)
これをそのまま実装。\( \left( i, j \right) = \left( S_{i-1}, T_{j-1} \right) \) なので1を引くことを忘れずに。
(define (pairs-index i j) (letrec ((iter (lambda (i j) (cond ((> i j) (error "unexpected index " i j)) ((and (= i 0) (= j 0)) 0) ((= i j) (+ (iter (- i 1) (- j 1)) (expt 2 i))) ((= i (- j 1)) (+ (iter i (- j 1)) (expt 2 i))) (else (+ (iter i (- j 1)) (expt 2 (+ i 1)))))))) (iter (- i 1) (- j 1))))
実行結果。
racket@> (pairs-index 1 100) 197 racket@> (pairs-index 99 100) 950737950171172051122527404030 racket@> (pairs-index 100 100) 1267650600228229401496703205374
答え合わせ。
racket@> (stream-ref (pairs integers integers) 197) => '(1 100) racket@> (stream-ref (pairs integers integers) 950737950171172051122527404030) => いくら待っても返らない...
問題 3.67
上手く説明できないけど縦軸でも伸びていくように実装すればよい。
(define (pairs-ex s t) (cons-stream (list (stream-car s) (stream-car t)) (interleave (interleave (stream-map (lambda (x) (list (stream-car s) x)) (stream-cdr t)) (stream-map (lambda (x) (list x (stream-car t))) (stream-cdr s))) (pairs (stream-cdr s) (stream-cdr t)))))
テスト。
racket@> (map (lambda (k) (let ((p (stream-ref (pairs-ex integers integers) k))) (cons k (list p)))) (enumerate-interval 0 25)) '((0 (1 1)) (1 (1 2)) (2 (2 2)) (3 (2 1)) (4 (2 3)) (5 (1 3)) (6 (3 3)) (7 (3 1)) (8 (2 4)) (9 (1 4)) (10 (3 4)) (11 (4 1)) (12 (2 5)) (13 (1 5)) (14 (4 4)) (15 (5 1)) (16 (2 6)) (17 (1 6)) (18 (3 5)) (19 (6 1)) (20 (2 7)) (21 (1 7)) (22 (4 5)) (23 (7 1)) (24 (2 8)) (25 (1 8)))
問題 3.68
interleaveの呼び出しが無限に続き処理が返らないから。
問題 3.69
これも上手く説明できないけど、(pairs)がs,tの2次元の軸で伸びていくのに対し、さらにuという新たな軸を追加して3次元で伸びていくイメージ。
(define (triples s t u) (cons-stream (list (stream-car s) (stream-car t) (stream-car u)) (interleave (stream-map (lambda (x) (flatten (list (stream-car u) x))) (pairs (stream-cdr s) (stream-cdr t))) (triples (stream-cdr s) (stream-cdr t) (stream-cdr u)))))
テスト。
racket@> (map (lambda (k) (let ((p (stream-ref (triples integers integers integers) k))) (cons k (list p)))) (enumerate-interval 0 25)) '((0 (1 1 1)) (1 (1 2 2)) (2 (2 2 2)) (3 (1 2 3)) (4 (2 3 3)) (5 (1 3 3)) (6 (3 3 3)) (7 (1 2 4)) (8 (2 3 4)) (9 (1 3 4)) (10 (3 4 4)) (11 (1 2 5)) (12 (2 4 4)) (13 (1 4 4)) (14 (4 4 4)) (15 (1 2 6)) (16 (2 3 5)) (17 (1 3 5)) (18 (3 4 5)) (19 (1 2 7)) (20 (2 4 5)) (21 (1 4 5)) (22 (4 5 5)) (23 (1 2 8)) (24 (2 3 6)) (25 (1 3 6)))
このストリームに対しピタゴラスの定理を満たす要素をフィルタリング。
(define pythagoras (stream-filter (lambda (triple) (let ((x (car triple)) (y (cadr triple)) (z (caddr triple))) (= (+ (expt x 2) (expt y 2)) (expt z 2)))) (triples integers integers integers)))
結果を確認。
racket@> (map (lambda (k) (stream-ref pythagoras k)) (enumerate-interval 0 3)) '((3 4 5) (6 8 10) (5 12 13) (9 12 15))
問題 3.70
weight手続きが最初はよくわからなかったが、比較する際に要素の重み付けを行う手続きということみたい。
§3.5.2のmergeを参考に実装。重さが小さい方が前にくる。
(define (merge-weighted s1 s2 weight) (cond ((stream-null? s1) s2) ((stream-null? s2) s1) (else (let* ((s1-car (stream-car s1)) (s2-car (stream-car s2)) (w1 (weight s1-car)) (w2 (weight s2-car))) (if (<= w1 w2) (cons-stream s1-car (merge-weighted (stream-cdr s1) s2 weight)) (cons-stream s2-car (merge-weighted s1 (stream-cdr s2) weight)))))))
weight-pairsはmerge-weightedによる重みづけ元に順序づけされる。
(define (weight-pairs s t weight) (cons-stream (list (stream-car s) (stream-car t)) (merge-weighted (stream-map (lambda (x) (list (stream-car s) x)) (stream-cdr t)) (weight-pairs (stream-cdr s) (stream-cdr t) weight) weight)))
a.
2つの自然数の和が小さい順に並んだ対のストリームを作る。
(define p1 (weight-pairs integers integers (lambda (pair) (+ (car pair) (cadr pair)))))
結果を確認。
racket@> (map (lambda (n) (stream-ref p1 n)) (enumerate-interval 0 24)) '((1 1) (1 2) (1 3) (2 2) (1 4) (2 3) (1 5) (2 4) (3 3) (1 6) (2 5) (3 4) (1 7) (2 6) (3 5) (4 4) (1 8) (2 7) (3 6) (4 5) (1 9) (2 8) (3 7) (4 6) (5 5))
b.
(define (divisible? n) (or (eq? (remainder n 2) 0) (eq? (remainder n 3) 0) (eq? (remainder n 5) 0))) (define p2 (stream-filter (lambda (pair) (and (not (divisible? (car pair))) (not (divisible? (cadr pair))))) (weight-pairs integers integers (lambda (pair) (let ((i (car pair)) (j (cadr pair))) (+ (* 2 i) (* 3 j) (* 5 i j)))))))
結果を確認。いまいちピンとこないけど。
racket@> (map (lambda (n) (stream-ref p2 n)) (enumerate-interval 0 24)) '((1 1) (1 7) (1 11) (1 13) (1 17) (1 19) (1 23) (1 29) (1 31) (7 7) (1 37) (1 41) (1 43) (1 47) (1 49) (1 53) (7 11) (1 59) (1 61) (7 13) (1 67) (1 71) (1 73) (1 77) (1 79))
問題 3.71
Ramanujan数をググッてみたけどそのエピソードがおもしろい。
1918年2月ごろ、ラマヌジャンは療養所に入っており、見舞いに来たハーディは次のようなことを言った。
「乗ってきたタクシーのナンバーは1729だった。さして特徴のない数字だったよ」
これを聞いたラマヌジャンは、すぐさま次のように言った。
「そんなことはありません。とても興味深い数字です。それは2通りの2つの立方数の和で表せる最小の数です」
実は、1729は次のように表すことができる。
1729 = 123 + 13 = 103 + 93
すなわち、1729が「A = B3 + C3 = D3 + E3」という形で表すことのできる数 A のうち最小のものであることを、ラマヌジャンは即座に指摘したのである。
(中略)
この逸話のため、1729は俗にハーディ・ラマヌジャン数やタクシー数などと呼ばれており、スタートレックやフューチュラマなどのSFや、ハッカー文化の文脈では「一見すると特に意味のない数」のような文脈でこの数が使われていることがある。
まずは対の3乗の和の手続き。
(define (cube n) (* n n n)) (define (sum-of-cube pair) (+ (cube (car pair)) (cube (cadr pair))))
3乗の和の小さい順でソートした対の並びは(weight-pairs integers integers sum-of-cube)と書けるので、それの対に対し \( i_{k}^{3} + j_{k}^{3} = i_{k+1}^{3} + j_{k+1}^{3} \) を満たす対を見つければよい。
(define (ramanujan-filter s) (let* ((s1 (stream-ref s 0)) (s2 (stream-ref s 1)) (w1 (sum-of-cube s1)) (w2 (sum-of-cube s2))) (if (= w1 w2) (cons-stream (list w1 s1 s2) (ramanujan-filter (stream-cdr s))) (ramanujan-filter (stream-cdr s))))) (define ramanujan-numbers (ramanujan-filter (weight-pairs integers integers sum-of-cube)))
実行結果を確認。できてるっぽい。
racket@> (map (lambda (k) (stream-ref ramanujan-numbers k)) (enumerate-interval 0 4)) '((1729 (1 12) (9 10)) (4104 (2 16) (9 15)) (13832 (2 24) (18 20)) (20683 (10 27) (19 24)) (32832 (4 32) (18 30)))
問題 3.72
問題文がわかりにくだけで、やりたいことは問題 3.71とほとんどいっしょ。
(define (sum-of-square pair) (+ (square (car pair)) (square (cadr pair)))) (define (sum-of-squares-filter s) (let* ((s1 (stream-ref s 0)) (s2 (stream-ref s 1)) (s3 (stream-ref s 2)) (w1 (sum-of-square s1)) (w2 (sum-of-square s2)) (w3 (sum-of-square s3))) (if (= w1 w2 w3) (cons-stream (list w1 s1 s2 s3) (sum-of-squares-filter (stream-cdr s))) (sum-of-squares-filter (stream-cdr s))))) (define sum-of-square-numbers (sum-of-squares-filter (weight-pairs integers integers sum-of-square)))
ストリームを確認。まあこんなこもんかなといった感じ。
racket@> (map (lambda (k) (stream-ref sum-of-square-numbers k)) (enumerate-interval 0 15)) '((325 (1 18) (6 17) (10 15)) (425 (5 20) (8 19) (13 16)) (650 (5 25) (11 23) (17 19)) (725 (7 26) (10 25) (14 23)) (845 (2 29) (13 26) (19 22)) (850 (3 29) (11 27) (15 25)) (925 (5 30) (14 27) (21 22)) (1025 (1 32) (8 31) (20 25)) (1105 (4 33) (9 32) (12 31)) (1105 (9 32) (12 31) (23 24)) (1250 (5 35) (17 31) (25 25)) (1300 (2 36) (12 34) (20 30)) (1325 (10 35) (13 34) (22 29)) (1445 (1 38) (17 34) (22 31)) (1450 (9 37) (15 35) (19 33)) (1525 (2 39) (9 38) (25 30)))
信号としてのストリーム
離散値の積分は以下のように表現できる。
\( S_{i} = C + \sum_{j = 1}^{i} x_{i} dt \)
\( x_{i} \) をストリームintegrandとすると積分器は以下のように実装される。
(define (integral integrand initial-value dt) (define int (cons-stream initial-value (add-streams (scale-stream integrand dt) int))) int)
問題 3.73
テキストの図のRC回路をほんとそのまま実装すればよい。
(define (RC R C dt) (define (proc integrand v0) (add-streams (scale-stream integrand R) (integrand (scale-stream integrand (/ 1 C)) v0 dt))) proc)
確かめる方法が思いつかないので、テストはパスで。
問題 3.74
(define zero-crossings (stream-map sign-change-detector sense-data <expression>))
の<expression>を補完せよ。もうここまでヒントがあったら解く前からほとんど答えはわかるけど。
それにしてもEvaの洞察がすごい。こんな上司かっこいいな。
まずはAlyssaが実装したmake-zero-crossingsからsign-change-detectorをリーバスエンジニアリング。
(define (sign-change-detector x last) (cond ((and (< x 0) (> last 0)) -1) ((and (> x 0) (< last 0)) 1) (else 0)))
次にsense-dataを定義する。(有限だけど)
(define sense-data (list->stream (list 1 2 1.5 1 0.5 -0.1 -2 -3 -2 -0.5 0.2 3 4)))
またzero-crossingsの結果は、
1 2 1.5 1 0.5 -0.1 -2 -3 -2 -0.5 0.2 3 4 ... ? ? ? ? ? ? ? ? ? ? ? ? ? ... ↓ 0 0 0 0 0 -1 0 0 0 0 1 0 0 ...
となればよいので、zero-crossingsは以下のように定義される。
(define zero-crossings (stream-map sign-change-detector sense-data (cons-stream 0 sense-data)))
テスト。OK。
racket@> (map (lambda (i) (stream-ref zero-crossings i)) (enumerate-interval 0 12)) => '(0 0 0 0 0 -1 0 0 0 0 1 0 0)
問題 3.75
Louisの実装はmake-zero-crossingsが再帰で実行される際のlast-valueが前のavptとなっているため、avptが時刻tとt-1の平均ではない点にバグがある。
そこでavptが正しく算出されるように修正する。
(define (make-zero-crossings input-stream last-value last-avpt) (let ((avpt (/ (+ (stream-car input-stream) last-value) 2))) (cons-stream (sign-change-detector avpt last-avpt) (make-zero-crossings (stream-cdr input-stream) (stream-car input-stream) avpt))))
テスト。たぶん合ってそう。
racket@> (map (lambda (i) (stream-ref (make-zero-crossings sense-data 0 0) i)) (enumerate-interval 0 12)) => '(0 0 0 0 0 0 -1 0 0 0 0 1 0)
問題 3.76
問題3.75の実装では、入力信号を平滑化する処理とゼロ交差を検出する処理が分離されていないので、個別のフィルタ関数としてリファクタリングする。
平滑化する手続きsmoothは、
(define (average x y) (/ (+ x y) 2)) (define (smooth input-stream) (stream-map average input-stream (cons-stream 0 input-stream)))
ゼロ交差の検出結果を返す手続きmake-zero-crossingsは、問題 3.74を応用すればよくて、
(define (make-zero-crossings input-stream) (stream-map sign-change-detector input-stream (cons-stream 0 input-stream)))
テスト。
racket@> (map (lambda (i) (stream-ref (smooth sense-data) i)) (enumerate-interval 0 12)) => '(1/2 3/2 1.75 1.25 0.75 0.2 -1.05 -5/2 -5/2 -1.25 -0.15 1.6 7/2) racket@> (map (lambda (i) (stream-ref (make-zero-crossings (smooth sense-data)) i)) (enumerate-interval 0 12)) => '(0 0 0 0 0 0 -1 0 0 0 0 1 0)
楽しいぜ。
次回は「§3.5.4 ストリームと遅延評価」から。
