結城浩のSICP日記 RSSフィード

2006-05-20

調和級数(2) 調和級数(2) - 結城浩のSICP日記 を含むブックマーク

級数SICPの索引で調べてみたら問題が見つかった。

問題3.55: 部分和partial-sumsを定義する。(p.195)

そうか、数列を定義して、それを元に部分和を作るのは正しいよな。

(use util.stream)

(define (integers-starting-from n)
  (stream-cons n
    (integers-starting-from (+ n 1))))

(define naturals (integers-starting-from 1))

;問題 3.55
;------------------
(define (partial-sums s)
  (define (s-map proc s)
    (stream-cons
      (proc (stream-car s))
      (s-map proc (stream-cdr s))))
  (stream-cons
    (stream-car s)
    (s-map (lambda (n) (+ (stream-car s) n))
           (partial-sums (stream-cdr s)))))
;------------------

(stream->list (stream-take (partial-sums naturals) 10))
;=> (1 3 6 10 15 21 28 36 45 55)

(define (harmonic)
  (define (s n)
    (stream-cons (/ 1 n) (s (+ n 1))))
  (partial-sums (stream-cons 0 (s 1))))

さっきの表も作ろうか。

(let loop ((n 0) (s (harmonic)))
  (if (< n 10)
      (begin
        (display n)
        (display " : Hn = ")
        (print (stream-car s))
        (loop (+ n 1) (stream-cdr s)))))

実行結果です。

0 : Hn = 0
1 : Hn = 1
2 : Hn = 1.5
3 : Hn = 1.8333333333333333
4 : Hn = 2.083333333333333
5 : Hn = 2.283333333333333
6 : Hn = 2.45
7 : Hn = 2.5928571428571425
8 : Hn = 2.7178571428571425
9 : Hn = 2.828968253968254

調和級数(1) 調和級数(1) - 結城浩のSICP日記 を含むブックマーク

テトラちゃんとハーモニック・ナンバーより、調和数の話。

(use util.stream)

(define (harmonic)
  (let loop ((n 1) (s 0))
    (stream-cons
      s
      (loop (+ n 1) (+ s (/ 1 n))))))

(stream-car (harmonic))
;=> 0

(stream-ref (harmonic) 10)
;=> 2.9289682539682538

表にしてみましょうか。

(let loop ((n 0) (s (harmonic)))
  (if (< n 10)
      (begin
        (display n)
        (display " : Hn = ")
        (print (stream-car s))
        (loop (+ n 1) (stream-cdr s)))))

実行結果です。

0 : Hn = 0
1 : Hn = 1
2 : Hn = 1.5
3 : Hn = 1.8333333333333333
4 : Hn = 2.083333333333333
5 : Hn = 2.283333333333333
6 : Hn = 2.4499999999999997
7 : Hn = 2.5928571428571425
8 : Hn = 2.7178571428571425
9 : Hn = 2.8289682539682537

ついでに、\sum_{k=1}^{n} \frac{1}{k} - \log nの表も作ってみましょう。

(let loop ((n 0) (s (harmonic)))
  (if (< n 10)
      (begin
        (display n)
        (display " : Hn - log n = ")
        (print (- (stream-car s) (log n)))
        (loop (+ n 1) (stream-cdr s)))))

実行結果です。おっと、logに0を渡してしまった。

0 : Hn - log n = #i1/0
1 : Hn - log n = 1.0
2 : Hn - log n = 0.8068528194400547
3 : Hn - log n = 0.7347210446652235
4 : Hn - log n = 0.6970389722134425
5 : Hn - log n = 0.6738954208992329
6 : Hn - log n = 0.6582405307719448
7 : Hn - log n = 0.6469469938018293
8 : Hn - log n = 0.6384156011773068
9 : Hn - log n = 0.6317436766320341

大きな数でのHn - log nを求めてみます。

(define (euler-mascheroni n)
        (- (stream-ref (harmonic) n) (log n)))

(euler-mascheroni 100)
;=> 0.5822073316515288

(euler-mascheroni 1000)
;=> 0.5777155815682065

(euler-mascheroni 10000)
;=> 0.577265664068161

(euler-mascheroni 100000)
;=> 0.5772206648930531

(euler-mascheroni 1000000)
;=> 0.5772161649005074

……ってここまでやってきて気づいたんだけれど、ループで計算するんだったら、別にストリームにする必要ないじゃん。がーん。

トラックバック - http://sicp.g.hatena.ne.jp/hyuki/20060520