今度は, を使って, の真の値を求めるよー.といってもほとんど SICP にやりかたが書いてあるので,ぼくは Scheme を書くだけだ.
まずはこれが上記の総和をとってるやつの各項をストリームにしたもの.なるほどとても直感的だ.
(define (ln2-summands n) (cons-stream (/ 1.0 n) (stream-map - (ln2-summands (+ n 1)))))
次に部分和を取る.
(define ln2-stream (partial-sums (ln2-summands 1)))
今 ln2-summands の第 i 要素を とすると,ストリーム ln2-stream の第 n 要素は になってるよー.ということは とすると の値が求められているはず…!!
dump してみる.
(dump-stream ln2-stream 10) ; 1.0, 0.5, 0.8333333333333333, 0.5833333333333333, 0.7833333333333332, 0.6166666666666666, 0.7595238095238095, 0.6345238095238095, 0.7456349206349207, 0.6456349206349207, done
このように,意外と収束がおそい.しかし世界はよくできていて,オイラーが,連続する二項の符号が逆になっている数列の部分和を早く計算できることを示しているらしくて,こんな風にしたら早い!!
これは SICP 中で euler-transform って関数として定義されている.
(define (euler-transform s) (let ((s0 (stream-ref s 0)) ; Sn-1 (s1 (stream-ref s 1)) ; Sn (s2 (stream-ref s 2))) ; Sn+1 (cons-stream (- s2 (/ (square (- s2 s1)) (+ s0 (* -2 s1) s2))) (euler-transform (stream-cdr s)))))
dump してみるよー
(dump-stream (euler-transform ln2-stream) 10) ; 0.7, 0.6904761904761905, 0.6944444444444444, 0.6924242424242424, 0.6935897435897436, 0.6928571428571428, 0.6933473389355742, 0.6930033416875522, 0.6932539682539683, 0.6930657506744464, done
収束はやい!!
なんか次は tableau だか,そういうデータ構造が紹介されている.このへんよくわからん!!
(define (make-tableau transform s) (cons-stream s (make-tableau transform (transform s)))) (define (accelerated-sequence transform s) (stream-map stream-car (make-tableau transform s)))
これつかったらもっと収束早い!!
(dump-stream (accelerated-sequence euler-transform ln2-stream) 10) ; 1.0, 0.7, 0.6932773109243697, 0.6931488693329254, 0.6931471960735491, 0.6931471806635636, 0.6931471805604039, 0.6931471805599445, 0.6931471805599427, 0.6931471805599454, done