[SICP] 問題 1.29 : Simpsonの公式

Simpsonの公式は上に示したのより更に正確な数値積分法である.
Simpsonの公式を使えば, aからbまでの関数fの積分は, 偶整数nに対し, h = (b - a)/nまたy_k = f(a + kh)として

(h/3)[y₀ + 4y₁ + 2y₂ + 4y₃ + 2y₄ + ・・・ + 2y_n-2 + 4y_n-1 + y_n]

で近似出来る. (nの増加は近似の精度を増加する.)
引数としてf, a, bとnをとり, Simpson公式を使って計算した積分値を返す手続きを定義せよ. その手続きを使って(n = 100とn = 1000で)0から1までcubeを積分し, またその結果を上のintegral手続きの結果と比較せよ.

Simpsonの公式の大カッコの中身をじっと見つめると, 次のようなパターンが見えてきます.

(y₀ + 4y₁ + y₂) + (y₂ + 4y₃ + y₄) + ・・・ + (y_n-4 + 4y_n-3 + y_n-2) + (y_n-2 + 4y_n-1 + y_n)

y₂に注目してください. 公式では係数が2となっていますが, それを2つに分けています.

つまり, 大カッコの中は３つの項をひとまとめにするとスッキリと書けます. 手続きを定義します.

(require (lib "racket/trace.ss"))

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

(define (simpson f a b n)
  (define h (/ (- b a) n))
  (define (y k) (f (+ a (* k h))))
  (define (term k) (+ (y k) (* 4.0 (y (+ k 1))) (y (+ k 2))))
  (define (next k) (+ k 2))
  (* (/ h 3.0)
     (sum term 0 next (- n 2))))

0から1までのcubeを積分してみます.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (integral cube 0 1 0.01)
0.24998750000000042
> (integral cube 0 1 0.001)
0.249999875000001
> (simpson cube 0 1 100)
0.2500000000000001
> (simpson cube 0 1 1000)
0.2500000000000001
> 

計算結果を見ると精度が全く異なっています.

[SICP] 問題 1.17 : 対数的ステップ数の乗算手続き（再帰的プロセス）

本節のべき乗アルゴリズムは, 乗算の繰返しによるべき乗の実行に基づいていた. 同様に整数の乗算を加算の繰返しで実行出来る. 次の乗算手続きは(この言語には加算はあるが, 乗算はないと仮定する)expt手続きに似たものである:

(define (* a b)
  (if (= b 0)
      0
      (+ a (* a (- b 1)))))

このアルゴリズムはbに線形のステップ数をとる. 加算の他に整数を二倍するdouble演算と(偶数の)整数を2で割るhalve演算もあるとしよう. これらを使い, fast-exptと類似の対数的ステップ数の乗算手続きを設計せよ.

p.25のfast-exptを参考に手続きを定義します.

(require (lib "racket/trace.ss"))

(define (double x) (* x 2))

(define (halve x) (/ x 2))

(define (fast-mult a b)
  (cond ((< b 0) (* -1 (fast-mult a (abs b))))
        ((= b 0) 0)
        ((= b 1) a)
        ((even? b) (double (fast-mult a (halve b))))
        (else (+ a (fast-mult a (- b 1))))))

(trace fast-mult)

fast-exptと同じような構造をしています. 負の数をかけた場合と０をかけた場合の扱いを追加しています.

a×bの計算について, 奇数と偶数の場合に分けて, 計算をどのように進めるかを考えます.

1) bが奇数の場合, b = m+1とする. 次の式が成立する.

a×(m+1) = a + (a×m)

2) nが偶数の場合, b = 2mとする. 次の式が成立する.

a×(2m) = double(a×m)

実行結果は次のとおりです.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (fast-mult 3 45)
>(fast-mult 3 45)
> (fast-mult 3 44)
> >(fast-mult 3 22)
> > (fast-mult 3 11)
> > >(fast-mult 3 10)
> > > (fast-mult 3 5)
> > > >(fast-mult 3 4)
> > > > (fast-mult 3 2)
> > > > >(fast-mult 3 1)
< < < < <3
< < < < 6
< < < <12
< < < 15
< < <30
< < 33
< <66
< 132
<135
135
> (fast-mult 4 0)
>(fast-mult 4 0)
<0
0
> (fast-mult 5 -8)
>(fast-mult 5 -8)
> (fast-mult 5 8)
> >(fast-mult 5 4)
> > (fast-mult 5 2)
> > >(fast-mult 5 1)
< < <5
< < 10
< <20
< 40
<-40
-40
> 

トレースの結果から、対数的ステップ数で計算していることが分かります.

[SICP] 問題 1.16 : べき乗を求める反復的アルゴリズム

fast-exptのように, 逐次平方を使い, 対数的ステップ数の反復的べき乗プロセスを生成する手続きを設計せよ.

(ヒント: (b^n/2)² = (b²)^n/2に注意し, 指数nと底bの他にもう一つの状態変数aを用意する.
状態の移変りで積abⁿが不変であるように状態遷移を定義する.
プロセス開始時にaを1とし, プロセス終了時にaが結果になるようにする.
一般に, 状態の移変りに不変のままの不変量(invariant quantity)を定義する技法は, 反復的アルゴリズムの設計に有用である.)

問題文のヒントを参考に手続きを定義します.

(require (lib "racket/trace.ss"))

(define (square x) (* x x))

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

(define (fast-expt b n)
  (fast-expt-iter b n 1))

(trace fast-expt-iter)

p×bⁿを考え, pに結果をため込むようにします. 奇数と偶数の場合に分けて, pとbがどのように変化するかを考えます.

1) nが奇数の場合, n = m+1とする. 次の式が成立する.

p×b^m+1 = (p×b)×b^m

右辺のカッコ内の(p×b)が次のステップのpとなる.

2) nが偶数の場合, n = 2mとする. 次の式が成立する.

p×b^2m = p×(b²)^m

右辺のカッコ内の(b²)が次のステップのbとなる.

実行結果は次のとおりです.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (fast-expt 2 8)
>(fast-expt-iter 2 8 1)
>(fast-expt-iter 4 4 1)
>(fast-expt-iter 16 2 1)
>(fast-expt-iter 256 1 1)
>(fast-expt-iter 256 0 256)
<256
256
> 

トレースの結果から, 反復的プロセスを生成しており, 対数的ステップ数で計算していることが分かります.

[SICP] 問題 1.15 : 正弦の計算

(ラジアンで表す)角度の正弦は, xが十分小さい時, sin x ≒ xの近似と, 正弦の引数の値を小さくするための三角関係式

sin(x) = 3sin(x/3) - 4sin³(x/3)

を使って計算出来る. (この問題のためには, 角度の大きさが0.1ラジアンより大きくなければ「十分小さい」と考える.) この方法は次の手続きに採用してある:

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

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

(define (sine angle)
   (if (not (> (abs angle) 0.1))
       angle
       (p (sine (/ angle 3.0)))))

a. (sine 12.15)の評価で, 手続きpは何回作用させられたか.

b. (sine a)の評価で, 手続きsineの生成するプロセスが使うスペースとステップ数の増加の程度は(aの関数として)何か.

問題文の手続きを定義します. 手続きpを作用させた回数を調べるためにトレースを使います.

(require (lib "racket/trace.ss"))

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

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

(define (sine angle)
  (if (not (> (abs angle) 0.1))
      angle
      (p (sine (/ angle 3.0)))))

(trace p)

計算結果は次のとおりです.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (sin 12.5)
-0.06632189735120068
> (sine 12.5)
>(p 0.051440329218107005)
<0.153776521096694
>(p 0.153776521096694)
<0.4467840153484446
>(p 0.4467840153484446)
<0.9836111719853284
>(p 0.9836111719853284)
<-0.8557060643295382
>(p -0.8557060643295382)
<-0.060813768577286265
-0.060813768577286265
> 

トレースの結果から, pを作用させた回数は５回です.

手続きの定義から, aの値が0.1より大きい時は手続きを再帰的に呼び出し, 次のaの値は1/3になっています. 「十分小さい」と判断する値をe, ステップ数をnとすると次の関係が成り立ちます.

(a/3)n	<	e
a/e	<	3n
log3(a/e)	<	n

このことから, ステップ数の増加の程度はlog aであるといえます. また, 手続きは末尾再帰になっていないため, スペースの増加の程度もステップの増加の程度と同じです.

[筋トレ] 2014.10.20の週

今週の記録です。

10/21（火）74.05kg
　スクワット　　　２０、２０　　（5kgの鉄アレイ）
　リアレイズ　　　１５、　８　　（5kgの鉄アレイ）
　レッグレイズ　　１０、　７

10/22（水）73.85kg
　プッシュアップ　１７、　７　　（膝をついた状態）
　ナロースタンス　　０　　　　　（両手の間を拳一つ開けた状態）
　アームカール　　　６、　４　　（10kgのダンベル）

10/24（金）75.40kg
　スクワット　　　２０、２０　　（5kgの鉄アレイ）
　リアレイズ　　　１２、　７　　（5kgの鉄アレイ）
　レッグレイズ　　１３、　５

10/26（日）75.65kg
　プッシュアップ　１６、　４　　（膝をついた状態）
　ナロースタンス　　２、　３　　（両手の間を拳一つ開けた状態）
　アームカール　　　６、　４　　（10kgのダンベル）

スクワットとアームカールに使うウエイトを増やしました。スクワットは20回できてしまうので、ウエイトを増やしても良さそうです。

リアレイズについては、正しいフォームを意識したため、回数が減っています。数をこなせば良いというものではないはずです。

トレーニングは、朝にしており、体重は夜、風呂を上がった後で量っています。金曜日は飲み会だったので、体重は多めです。

更に、土日の夕食がおでんだったので、体重は増える一方です(~_~;)

[SICP] 問題 1.14 : 両替の計算が生成するプロセス

1.2.2節のcount-change手続きが11セントの両替の場合に生成するプロセスを示す木構造を描け. 両替の金額の増加につれ, このプロセスが使うスペースとステップ数の増加の程度は何か.

図を書きます. 木の節点に金額と両替に使えるコインを括弧付きで描いています.

 11 -----------------  -39
(50, 25, 10, 5, 1)     (50, 25, 10, 5, 1)
  |
  |
 11 -----------------  -14
(25, 10, 5, 1)         (25, 10, 5, 1)
  |
  |
 11 --------- 1 ------  -9
(10, 5, 1)  (10, 5, 1) (10, 5, 1)
  |           |
  |           |
  |           1 ------  -4
  |          (5, 1)     (5, 1)
  |           |
  |           |
  |           1  ------  0
  |          (1)        (1)
  |
 11 --------  6 -------  1 ----  -4
 (5, 1)      (5, 1)     (5, 1)   (5, 1)
  |           |          |
  |           |          |
  |           |          1 -----  0
  |           |         (1)      (1)
  |           |          |
  |           |          |
  |           |          1
  |           |         ( )
  |           |
  |           6 -- 5 -- 4 -- 3 -- 2 -- 1 -- 0
  |          (1)  (1)  (1)  (1)  (1)  (1)  (1)
  |           |    |    |    |    |    |
  |           |    |    |    |    |    |
  |           6    5    4    3    2    1
  |          ( )  ( )  ( )  ( )  ( )  ( )
  | 
 11 -- 10 -- 9 -- 8 -- 7 -- 6 -- 5 -- 4 -- 3 -- 2 -- 1 -- 0
 (1)   (1)  (1)  (1)  (1)  (1)  (1)  (1)  (1)  (1)  (1)  (1)
  |     |    |    |    |    |    |    |    |    |    |
  |     |    |    |    |    |    |    |    |    |    |
 11    10    9    8    7    6    5    4    3    2    1
 ( )   ( )  ( )  ( )  ( )  ( )  ( )  ( )  ( )  ( )  ( )

このプロセスが使うスペースとステップ数の増加の程度については後で考えます.

[SICP] 問題 1.13 : フィボナッチ数を求める代数方程式

Φ= (1+√5)/2としてFib(n)がΦⁿ/√5に最も近い整数であることを証明せよ. ヒント: Ψ= (1-√5)/2とする. 帰納法とFibonacci数の定義(1.2.2節参照)を用い, Fib(n)=(Φⁿ-Ψⁿ)/√5を証明せよ.

まず, Fib(n)=(Φⁿ-Ψⁿ)/√5を証明する.

(1) n=1のとき

Fib(1)	=	{Φ1-Ψ1} / √5
	=	{(1+√5)/2 - (1-√5)/2} / √5
	=	2√5 / √5
	=	1

(2) n=2のとき

Φ2	=	{(1+√5)/2}2
	=	{1 + 2√5 + 5} / 4
	=	{2 + 2√5 + 4} / 4
	=	{1 + √5} / 2 + 1
	=	Φ + 1

Ψ2	=	{(1-√5)/2}2
	=	{1 - 2√5 + 5} / 4
	=	{2 - 2√5 + 4} / 4
	=	{1 - √5} / 2 + 1
	=	Ψ + 1

Fib(2)	=	{Φ2-Ψ2} / √5
	=	{(Φ + 1) - (Ψ + 1)} / √5
	=	{Φ - Ψ} / √5
	=	1

(3) n=k, k+1のときFib(n)=(Φⁿ-Ψⁿ)/√5が成立すると仮定する. n=k+2のとき

Fib(k+2)	=	Fib(k+1) + Fib(k)
	=	(Φk+1-Ψk+1)/√5 + (Φk-Ψk)/√5
	=	{(Φk+1+Φk) - (Ψk+1+Ψk)} / √5
	=	{Φk(Φ+1) - Ψk(Ψ+1)} / √5
	=	{ΦkΦ2 - ΨkΨ2} / √5
	=	{Φk+2 - Ψk+2} / √5

(4) 以上より, 1以上の整数nについて, Fib(n)=(Φⁿ-Ψⁿ)/√5が成立する.

次に, Ψ/√5の値を計算する.

Ψ/√5 = (1-√5)/2 ≒ -0.2764

Ψ/√5 < 1/2 であるから, 1以上の整数nについて, Ψⁿ/√5 < 1/2 である.
よって、Fib(n)がΦⁿ/√5に最も近い整数であるといえる.