問題3.62-3.69 - nrvctの日記

ストリーム続きです。

3.62

(define (div-series s1 s2)
  (if (= (stream-car s2) 0)
      (error "s2 has 0 constant term" s2)
      (mul-series s1 
                  (scale-stream (invert-unit-series s2) (/ 1 (stream-car s2))))))

> (define tangent-series (div-series sine-series cosine-series))
> (disp-stream-n tangent-series 10)
(0 1 0 1/3 0 2/15 0 17/315 0 62/2835)

3.63

(define (sqrt-stream x)
  (define guesses
    (cons-stream 1.0
                 (stream-map (lambda (guess)
                               (sqrt-improve guess x))
                             guesses)))
  guesses)
(define (sqrt-stream-2 x)
  (cons-stream 1.0
               (stream-map (lambda (guess)
                             (sqrt-improve guess x))
                           (sqrt-stream-2 x))))

> (time (stream-ref (sqrt-stream 2) 100))
cpu time: 4 real time: 2 gc time: 0
1.414213562373095
> (time (stream-ref (sqrt-stream-2 2) 100))
cpu time: 60 real time: 60 gc time: 0
1.414213562373095
> (time (stream-ref (sqrt-stream 2) 1000))
cpu time: 12 real time: 16 gc time: 0
1.414213562373095
> (time (stream-ref (sqrt-stream-2 2) 1000))
cpu time: 7228 real time: 7242 gc time: 1076
1.414213562373095

sqrt-stream-2では手続きの再帰的呼び出しによってストリームを作っている。
つまり新しいストリームを何度も作り出すことになるので、メモ化が働かず同じ値を何度も計算してしまうことになる。よって非効率。
sqrt-streamでは同じストリームguessesを使うので、メモ化が働き既知の値を計算しなおすということがない。メモ化を外すとsqrt-stream-2と同じくらい時間がかかるようになる。

(define-syntax delay
  (syntax-rules ()
    ((_ a) (lambda () a))))
(define (force a) (a))

> (time (stream-ref (sqrt-stream 2) 100))
cpu time: 48 real time: 50 gc time: 0
1.414213562373095
> (time (stream-ref (sqrt-stream-2 2) 100))
cpu time: 48 real time: 50 gc time: 0
1.414213562373095
> (time (stream-ref (sqrt-stream 2) 1000))
cpu time: 5152 real time: 5150 gc time: 248
1.414213562373095
> (time (stream-ref (sqrt-stream-2 2) 1000))
cpu time: 5036 real time: 5044 gc time: 136
1.414213562373095

3.64

(define (stream-limit s tolerance)
  (let ((s0 (stream-car s)) (s1 (stream-car (stream-cdr s))))
    (if (< (abs (- s0 s1)) tolerance)
        s1
        (stream-limit (stream-cdr s) tolerance))))

> (stream-limit (sqrt-stream 2) 0.0001)
1.4142135623746899
> (stream-limit (sqrt-stream 2) 0.000000001)
1.414213562373095

3.65

3つの並びってなんだろ。級数をそのまま計算するのと、オイラーの加速技法を使うのと、さらにタブローの先頭要素を使ったものということかな。

(define ln-series-terms
  (stream-map (lambda (x) (if (even? x) (- 0.0 (/ 1 x)) (/ 1 x))) integers))
(define ln2-1 (partial-sums ln-series-terms))
(define ln2-2 (euler-transform ln2-1))
(define ln2-3 (accelerated-sequence euler-transform ln2-1))

級数そのまま

> (disp-stream-n ln2-1 20)
(1
 0.5
 0.8333333333333334
 0.5833333333333334
 0.7833333333333333
 0.6166666666666667
 ...(略)
 0.6661398242280596
 0.7187714031754281
 0.668771403175428)

一段階オイラー加速

> (disp-stream-n ln2-2 20)
(0.7
 0.6904761904761905
 0.6944444444444444
 0.6924242424242424
 0.6935897435897436
 0.692857142857143
 ...(略)
 0.6931303775344024
 0.6931616470778671
 0.6931346368409873)

タブロー加速

> (disp-stream-n ln2-3 20)
(1
 0.7
 0.6932773109243697
 0.6931488693329254
 0.6931471960735491
 0.6931471806635636
 0.693147180560404
 0.6931471805599446
 0.6931471805599428
 0.6931471805599454
 +nan.0
 +nan.0
 ...(略)
 +nan.0
 +nan.0)

tableau加速は第11項で0割が発生してしまいました。オイラー加速の式を見ると分かりますが、S_{n-1} - 2S_n + S_{n+1}が非常に小さくなると0割が起きます。
たったの第11項で発生するなんて、すごく早く収束してる証拠ですね。

3.66

interleaveを初めに思いついた人すごいなー。
(pairs integers integers)での、対の生成順序について考える。

> (disp-stream-n sp 30)
((1 1)
 (1 2)
 (2 2)
 (1 3)
 (2 3)
 (1 4)
 (3 3)
 (1 5)
 (2 4)
 (1 6)
 (3 4)
 (1 7)
 (2 5)
 (1 8)
 (4 4)
 (1 9)
 (2 6)
 (1 10)
 (3 5)
 (1 11)
 (2 7)
 (1 12)
 (4 5)
 (1 13)
 (2 8)
 (1 14)
 (3 6)
 (1 15)
 (2 9)
 (1 16))

便宜上、(i j)を行列のi行j列として考えることにする。
見れば分かる通り、始めを除いて第1行の要素は1つおきに現れるので、(1 100)は2(100-1)=198回目に現れる。つまりそれ以前の対は197個(といってもこれが本当に正しいかどうかは証明しないとだめだけれど)。
同様に考えてみると、第2行は始めを除いて3つおきに現れていることが分かるし、第3行も始めを除いて7つおきに現れている。第4行は15個おき。
そうすると、1,3,7,15,...の並びを一般化できるか考えてみる。階差をとると2のべき乗になっていることが分かる。
そうすると、第i行の並びが始め以外はm_i個おきであるとしたとき、

$m_i=1+\sum_{k=1}^{i-1}2^k$
となるはず。
次に、i=jとなる(i j)が初めに出てくるのは何番目なのかですが、これは上の実行結果を眺めてると分かりますが、これも1,3,7,15,...の並びになっています。さらに、始めの(i j)と(i j+1)の間も、1,3,7,15,...を一つずらしたもの、つまりi-1の時の始めでないところ、(i-1 j)と(i-1 j+1)の間隔と同じになります。
したがって、一般に対(i j) (i≦j) が現れるのをf(i, j)とすると、
$\begin{eqnarray*}(j-1 > i \geq 2) \hspace{10pt} f(i, j) &=& m_i + (m_{i-1} + 1) + (j-i-1)(m_{i} + 1)\\ &=& 3 + 2^{i-1} + \sum_{k=2}^{i-1} + (j-i-1)(4 + \sum_{k=2}^{i-1}2^k)\\ &=& 3 + 4(j-i-1) + 2^{i-1} + (j-i)\sum_{k=2}^{i-1}2^k \\ (j-1 = i \geq 2) \hspace{10pt} f(i, j) &=& m_i + (m_{i-1} + 1) = 1 + \sum_{k=1}^{i-1}2^k + (2 + \sum_{k=1}^{i-2}2^k) \\ &=& 3 + 2^{i-1} + \sum_{k=2}^{i-1}2^k \\ (j = i \geq 2) \hspace{10pt} f(i, j) &=& m_i = 1 + \sum_{k=1}^{i-1}2^k \\ (j > i=1) \hspace{10pt} f(i, j) &=& 2(j-i) \\ (j=i=1) \hspace{10pt} f(i, j) &=& 1\end{eqnarray*}$
というのが推測されます。証明は数学的帰納法を使えばいいと思いますが、非常に面倒くさそうなのでやりません。
(12/25 11:30追記
2のべき乗の和って等比級数だから簡単に計算できますね。まあもうf作っちゃったしTeX書くの面倒なので上式を変える気はあまりないんですが、ド忘れしてたのはちょっとショックです…＾＾；
)
f(i, j)をscheme上で実装します。

(define expt2-stream (cons-stream 1 (add-streams expt2-stream expt2-stream)))
(define (f i j)
  (cond ((= i 1)
         (if (= i j)
             1
             (* 2 (- j i))))
        ((= i j) 
         (+ 1 (stream-ref (partial-sums (stream-cdr expt2-stream)) (- (- i 1) 1))))
        ((= (- j 1) i)
         (+ 3
            (stream-ref expt2-stream (- i 1))
            (if (< (- i 1) 2)
                0
                (stream-ref (partial-sums (stream-cdr (stream-cdr expt2-stream)))
                            (- (- i 1) 2)))))
        (else
         (+ 3
            (* 4 (- j i 1))
            (stream-ref expt2-stream (- i 1))
            (* (- j i)
               (if (< (- i 1) 2)
                   0
                   (stream-ref (partial-sums (stream-cdr (stream-cdr expt2-stream)))
                               (- (- i 1) 2))))))))

> (f 2 4)
9
> (f 1 6)
10
> (f 3 4)
11
> (f 2 6)
17

上の(pairs integers integers)が生成したものと比べると、fはとりあえず正しい挙動を示しています。
これで、対(1, 100), (99, 100), (100, 100)が何番目に出てくるのかが分かります。

> (f 1 100)
198
> (f 99 100)
950737950171172051122527404031
> (f 100 100)
1267650600228229401496703205375

＾＾；
ということで、iの値が大きくなると、その対が出てくるのは非常に後になるということがわかりました。これ使えるんでしょうか…。

3.67

うーん難しい。

(define (pairs-2 s t)
  (let ((u integers))
    (cons-stream
     (list (stream-car s) (stream-car t))
     (interleave
      (stream-map (lambda (x) (list x (stream-car u)))
                  (stream-cdr t))
      (interleave
       (stream-map (lambda (x) (list (stream-car s) x))
                   (stream-cdr t))
       (pairs-2 (stream-cdr s) (stream-cdr t)))))))
(define count 0)
(define (have-pair? pair stream)
  (set! count (+ count 1))
  (if (equal? pair (stream-car stream))
      (begin (display count)
             (set! count 0))
      (have-pair? pair (stream-cdr stream))))

> (define ps2 (pairs-2 integers integers))
> (have-pair? '(3 8) ps2)
309
> (have-pair? '(8 9) ps2)
54613
> (have-pair? '(10 11) ps2)
873813
> (have-pair? '(2 9) ps2)
109
> (have-pair? '(1 100) ps2)
395
> (have-pair? '(3 8) ps2)
309

pairsの意味を考えて書いてみるとたぶんこんな感じのような気がします。ややこしすぎてあまり自信がありません。

3.68

(define (pairs-3 s t)
  (interleave
   (stream-map (lambda (x) (list (stream-car s) x))
               t)
   (pairs-3 (stream-cdr s) (stream-cdr t))))

> (disp-stream-n (pairs-3 integers integers) 20)

実行すると無限ループに入ります。
たぶんこういうことだと思います。

(pairs-3 integers integers)
=>
(interleave
 (stream-map (lambda (x) (list (stream-car s) x))
             t)
 (interleave
  (stream-map (lambda (x) (list (stream-car (stream-cdr s) x)))
              (stream-cdr t))
  (interleave ...

つまり、interleaveの第2引数に値を渡せないため、無限ループに陥るのだと思います。

3.69

分かりません。頭痛い。
他の人のを参考にしました。pairs使えばいいのね…。

(define (triples s t u)
  (cons-stream
   (list (stream-car s) (stream-car t) (stream-car u))
   (interleave
    (stream-map (lambda (x) (list (stream-car s) (stream-car t) x))
                (stream-cdr u))
    (interleave
     (stream-map (lambda (x) (cons (stream-car s) x))
                 (pairs (stream-cdr t) (stream-cdr u)))
     (triples (stream-cdr s) (stream-cdr t) (stream-cdr u))))))

(define pythagoras
  (stream-filter
   (lambda (triples)
     (= (+ (sqr (car triples)) (sqr (cadr triples)))
        (sqr (caddr triples))))
   (triples integers integers integers)))

これで実行すると、

> (time (disp-stream-n pythagoras 1))
cpu time: 1532 real time: 1534 gc time: 0
((3 4 5))
> (time (disp-stream-n pythagoras 2))
cpu time: 9229 real time: 9234 gc time: 2256
((3 4 5) (6 8 10))

なんとかでます。しかしとてつもなく遅い。
ピタゴラス数の並びは表を(s0,t0,u0)を頂点とした四角錐と見ると、大体真ん中ら辺にあると思われますが、triplesやpairsで生成する並びのストリームは問題3.66で確かめたように、端っこが早くでてきて、真ん中のあたりの要素はとても後に出てきます。これが遅くなる原因だと思われます。
確かめてみます。次のようにカウンター変数を置きます。

(define count 0)
(define pythagoras
  (stream-filter
   (lambda (triples)
     (begin (set! count (+ count 1))
            (= (+ (sqr (car triples)) (sqr (cadr triples)))
               (sqr (caddr triples)))))
   (triples integers integers integers)))

すると

> (time (disp-stream-n pythagoras 1))
cpu time: 1880 real time: 1878 gc time: 320
((3 4 5))
> count
36181
> (set! count 0)
> (time (disp-stream-n pythagoras 2))
cpu time: 10281 real time: 10286 gc time: 3216
((3 4 5) (6 8 10))
> count
159232

で、非常に効率の悪いことをしていると分かります(ちなみに(disp-stream-n pythagoras 3)は3分放置しても動いてました)。
このことから、並びを作る順序はある程度目的に沿って決める必要があると分かります。

ということで、次の問題3.70は並びの目的に沿った生成順序に関するもののようですが、大分進めたので今日はここまでにします。