2017/11/29

[SICP] 問題 1.34 : (f f)の評価

手続き

(define (f g)
  (g 2))

を定義したとする. その時


(f square)
4

(f (lambda (z) (* z (+ z 1))))
6

解釈系に組合せ(f f)を(意地悪く)評価させるとどうなるか. 説明せよ. 

問題文の手続きを定義して実行してみます. 手続きは次の通り.
(define (f g)
  (g 2))

(define (square x) (* x x))
実行結果は次のとおりです.
ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (f square)
4
> (f (lambda (z) (* z (+ z 1))))
6
> (f f)
. . application: not a procedure;
 expected a procedure that can be applied to arguments
  given: 2
  arguments...:
   2
> 
(f f)を評価するとエラーが生じています. その理由を考えます.
手続きfの定義から(f f)を評価した場合の置き換えは次のようになります.
  (f f)
= (f 2)
= (2 2)
「2」は手続きではないのでエラーが生じます.

2017/11/26

[SICP] 1.2.6 : expmodでbaseのexp乗をmで割った余りが求まる理由

baseのexp乗をmで割った余りを求める手続きexpmodは次のような定義になります.
(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (remainder (square (expmod base (/ exp 2) m))
                    m))
        (else
         (remainder (* base (expmod base (- exp 1) m))
                    m))))
法演算について馴染みがないと, なぜ, この手続で余りが求まるか分からないと思います. ここで簡単に計算の仕組みを説明します.
まず, 整数aをmで割った余りをr, 整数bをmで割った余りをtとします. すると整数aとbは次のように表現できます.
   a = qm + r
   b = sm + t
このaとbの和をmで割った余りを求めます.
   a + b = (qm + r) + (sm + t)
         = (q + s)m + (r + t)
(q + s)mはmで割り切れるので, a + bをmで割った余りはr+tをmで割った余りと等しくなります.
aにある整数cを掛けた値をmで割った余りも求めてます.
   ca = cqm + cr
caをmで割った余りはcrをmで割った余りと等しくなります.
続けて, aとbの積をmで割った余りを求めます.
   ab = (qm + r)(sm + t)
      = qsm^2 + (qt + rs)m + rt
abをmで割った余りはrtをmで割った余りと等しくなります.
積についての結果から, aのx乗をmで割った余りはrのx乗をmで割った余りと等しくなります.
ここで, 手続きの定義に戻り, baseのexp乗をmで割った余りを考えます.
1)expが偶数のときを考えます.
baseの(exp/2)乗をmで割った余りをrとします. 上で計算した結果から, rの2乗をmで割った余りは, baseのexp乗をmで割った余りと等しくなります.
2) expが奇数のときを考えます.
baseの(exp-1)乗をmで割った余りをrとします. 上で計算した結果から, base*rをmで割った余りは, baseのexp乗をmで割った余りと等しくなります.
これらの結果から, 手続きexpmodによりbaseのexp乗をmで割った余りを求められることが分かります.

2017/11/21

[SICP] 問題 1.27 : カーマイケル数

脚注47のCarmichael数が本当にFermatテストをだますことを示せ. それには整数nをとり, a < nなるすべてのaで, anがnを法としてaの合同になるかどうか見る手続きを書き, Carmichael数でその手続きを使ってみる.
a < nを満たす正の整数aについて, expmodを用いてanがnを法としてaと合同になるか調べます. 手続きは次のとおりです.
(define (square x) (* x x))

(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (remainder (square (expmod base (/ exp 2) m))
                    m))
        (else
         (remainder (* base (expmod base (- exp 1) m))
                    m))))

(define (pass-fermat-test? n)
  (define (loop a)
    (or (= a n)
        (and (= (expmod a n n) a)
             (loop (+ a 1)))))
  (loop 1))
評価してみます.
ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (pass-fermat-test? 561)
#t
> (pass-fermat-test? 1105)
#t
> (pass-fermat-test? 1729)
#t
> (pass-fermat-test? 2465)
#t
> (pass-fermat-test? 2821)
#t
> (pass-fermat-test? 6601)
#t
> (pass-fermat-test? 6605)
#f
> (pass-fermat-test? 6607)
#t
> 
脚注のカーマイケル数は, フェルマー・テストをパスしています. 間違えて, 必ず#tを返す手続きを実装していては困りますので, 合成数6605では#fが返ってくることを確認しています.

2017/11/20

[SICP] 問題 1.26 : expmodの計算に必要なステップ数

Louis Reasonerは問題1.24がなかなかうまくいかなかった. 彼の fast-prime?はprime?よりずっと遅いようであった. Louisは友人のEva Lu Atorに助けを求めた. Louisのプログラムを見ていると, squareを呼び出す代りに乗算を陽に使うように変っていた.
(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (remainder (* (expmod base (/ exp 2) m)
                       (expmod base (/ exp 2) m))
                    m))
        (else
         (remainder (* base (expmod base (- exp 1) m))
                    m))))
「違いが分らん」とLouis, 「分る」とEvaがいった. 「手続きをこのように書き替えたので, Θ(log n)のプロセスをΘ(n)のプロセスにしてしまった.」 説明せよ.
引数expの値に注目して, 手続きが生成するプロセスを図示してみます. まず, squareを用いた場合.
13 -- 12 -- 6 -- 3 -- 2 -- 1 -- 0
まず, squareを用いず, 乗算を用いた場合.
13 -- 12 -- 6 -- 3 -- 2 -- 1 -- 0
       |    |         |
       |    |         + -- 1 -- 0
       |    |
       |    + -- 3 -- 2 -- 1 -- 0
       |              |
       |              + -- 1 -- 0
       |
       + -- 6 -- 3 -- 2 -- 1 -- 0
            |         |
            |         + -- 1 -- 0
            |
            + -- 3 -- 2 -- 1 -- 0
                      |
                      + -- 1 -- 0
このように, squareを用いずに乗算を用いた場合, expの値が偶数の場合, expmodを2回呼び出しています. そのため, Θ(log n)のプロセスをΘ(n)のプロセスにしてしまい, 計算に時間がかかるようになります.

2017/11/18

[SICP] 問題 1.25 : expmodの必要性

Alyssa P. Hackerはexpmodを書くのに多くの余分なことをやったと不満で, 結局べき乗の計算法は知っているから
(define (expmod base exp m)
  (remainder (fast-expt base exp) m))
と単純に書ける筈だといった. これは正しいか. これも高速素数テストと同じに使えるか, 説明せよ.
テキストp.29の脚注48にあるように, ここでは200桁の整数が素数であるかどうかを問題としています. expmodの引数であるbaseやexpに200桁の整数が与えられるため, fast-exptで計算することは困難であり現実的ではありません. そのため, p.28のexpmodの定義が必要となります.

2017/11/17

[SICP] 問題 1.24 : フェルマー・テスト

問題1.22のtimed-prime-test手続きをfast-prime?(Fermat法参照) を使うように変え, その問題で見つけた12個の素数をテストせよ. FermatテストはΘ(log n)の増加なので, 1,000,000近くの素数をテストする時間を1000近くの素数をテストする時間と比べ, どの位と予想するか. データはその予想を支持するか. 違いがあれば説明出来るか.
フェルマーの小定理を用いた素数性のテストをする手続きを定義します. 手続きrandomの引数の範囲に制限があるため, 少し小細工をしています.
(define (square x) (* x x))

(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (remainder (square (expmod base (/ exp 2) m))
                    m))
        (else
         (remainder (* base (expmod base (- exp 1) m))
                    m))))  

(define (rnd n)
  (remainder (* (random (remainder n 4294967087))
                (random (remainder n 4294967087)))
             n))

(define (fermat-test n)
  (define (try-it a)
    (= (expmod a n n) a))
  (try-it (+ 1 (rnd (- n 1)))))

(define (fast-prime? n times)
  (cond ((= times 0) #t)
        ((fermat-test n) (fast-prime? n (- times 1)))
        (else #f)))

;; 素数を求めるために必要な時間を計測する
(define (timed-prime-test n)
  (newline)
  (display n)
  (start-prime-test n (current-milliseconds)))

(define (start-prime-test n start-time)
  (if (fast-prime? n 10000)
      (report-prime (- (current-milliseconds) start-time))))

(define (report-prime elapsed-time)
  (display " *** ")
  (display elapsed-time))

(define (search-for-primes start end)
  (cond ((< start end)
         (timed-prime-test start)
         (search-for-primes (+ start 1) end))
        (else 
         (newline)
         (display "end"))))
 
素数を求めてみます.
ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (search-for-primes 100000000000000 100000000000100)
100000000000031 *** 325
100000000000067 *** 344
100000000000097 *** 312
end
> (search-for-primes 1000000000000000 1000000000000200)
1000000000000037 *** 342
1000000000000091 *** 351
1000000000000159 *** 340
end
> (search-for-primes 10000000000000000 10000000000000100)
10000000000000061 *** 390
10000000000000069 *** 375
10000000000000079 *** 366
end
> 
この程度のデータでは, 何も言えそうにありません.

2017/11/15

[SICP] 問題 1.23 : 素数性を判定する数を奇数に限定

本節の初めのsmallest-divisorは多くの不要な計算をする: 数が2で割り切れるか調べた後は, より大きい偶数で割り切れるか調べるのは無意味である. test-divisorが使う値は, 2, 3, 4, 5, 6, ...ではなく, 2, 3, 5, 7, 9, ...であるべきだ. この変更を実装するため, 入力が2なら3 を返し, そうでなければ入力に2足したものを返す手続きnextを定義せよ. smallest-divisorを(+ test-divisor 1)ではなく, (next test-divisor)を使うように修正せよ. このように修正した smallest-divisorにしたtimed-prime-testで, 問題1.22で見つけた12 個の素数をテストせよ. この修正はテストのステップを半分にしたので, 二倍速く走ることが期待される. この期待は確められるか. そうでなければ, 得られた二つのアルゴリズムの速度の比はどうであったか. それが2と違う事情を説明せよ.
問題文にあるように手続きを修正します.
(define (square x) (* x x))

(define (divides? a b)
  (= (remainder b a) 0))

;; 問題文により追加
(define (next n)
  (if (= n 2)
      3
      (+ n 2)))

(define (find-divisor n test-divisor)
  (cond ((> (square test-divisor) n) n)
        ((divides? test-divisor n) test-divisor)
        (else (find-divisor n (next test-divisor)))))

(define (smallest-divisor n)
  (find-divisor n 2))

(define (prime? n)
  (= n (smallest-divisor n)))

;; 素数を求めるために必要な時間を計測する
(define (timed-prime-test n)
  (newline)
  (display n)
  (start-prime-test n (current-milliseconds)))

(define (start-prime-test n start-time)
  (if (prime? n)
      (report-prime (- (current-milliseconds) start-time))))

(define (report-prime elapsed-time)
  (display " *** ")
  (display elapsed-time))

(define (search-for-primes start end)
  (cond ((< start end)
         (timed-prime-test start)
         (search-for-primes (+ start 1) end))
        (else 
         (newline)
         (display "end"))))
 

問題 1.22で求めた12個の素数を判定するための時間を計測します.
ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (search-for-primes 100000000000000 100000000000100)
100000000000031 *** 624  (939)
100000000000067 *** 571  (1397)
100000000000097 *** 598  (949)
end
> (search-for-primes 1000000000000000 1000000000000200)
1000000000000037 *** 1907  (2891)
1000000000000091 *** 1965  (2928)
1000000000000159 *** 1877  (2919)
end
> (search-for-primes 10000000000000000 10000000000000100)
10000000000000061 *** 5925  (9149)
10000000000000069 *** 5887  (9166)
10000000000000079 *** 5887  (9174)
end
> 
カッコ内に問題 1.22で計測した時間を書いています. 問題 1.22の結果と見比べると, 時間は短くなっていますが, 半分とまでは言えません. 手続きnextでnが2であるかどうかを判定しているからでしょう.

[SICP] 問題 1.22 : 素数を求めるために必要な計算量

殆んどのLispの実装には基本手続きruntimeがあり, システムが走った時間を(例えばマイクロ秒で)示す整数を返す. 次のtimed-prime-test手続きは整数nで呼び出されると, nを印字し, nが素数かどうかを調べ, nが素数なら手続きは三つの星印とテストに要した時間を印字する.

(define (timed-prime-test n)
  (newline)
  (display n)
  (start-prime-test n (runtime)))

(define (start-prime-test n start-time)
  (if (prime? n)
      (report-prime (- (runtime) start-time))))

(define (report-prime elapsed-time)
  (display " *** ")
  (display elapsed-time))

この手続きを使い, 指定範囲の連続する奇整数について素数性を調べる手続きsearch-for-primesを書け. その手続きで, それぞれ1,000, 10,000, 100,000, 1,000,000より大きい, 小さい方から三つの素数を見つけよ. 素数をテストする時間に注意せよ. テストのアルゴリズムはΘ(√n)の増加の程度だから, 10,000前後の素数のテストは1,000前後の素数のテストの√10倍かかると考えよ. 時間のデータはこれを支持するか. 100,000, 1,000,000のデータは√n予想のとおりだろうか. 結果はプログラムが計算に必要なステップ数に比例した時間で走るという考えに合っているか.
手続きruntimeの代わりにcurrent-millisecondsを使います. 手続きは次のようになります.
(define (square x) (* x x))

(define (divides? a b)
  (= (remainder b a) 0))

(define (find-divisor n test-divisor)
  (cond ((> (square test-divisor) n) n)
        ((divides? test-divisor n) test-divisor)
        (else (find-divisor n (+ test-divisor 1)))))

(define (smallest-divisor n)
  (find-divisor n 2))

(define (prime? n)
  (= n (smallest-divisor n)))

;; 素数を求めるために必要な時間を計測する
(define (timed-prime-test n)
  (newline)
  (display n)
  (start-prime-test n (current-milliseconds)))

(define (start-prime-test n start-time)
  (if (prime? n)
      (report-prime (- (current-milliseconds) start-time))))

(define (report-prime elapsed-time)
  (display " *** ")
  (display elapsed-time))

(define (search-for-primes start end)
  (cond ((< start end)
         (timed-prime-test start)
         (search-for-primes (+ start 1) end))
        (else 
         (newline)
         (display "end"))))
問題文では, 素数を求める範囲を1000から始めていますが, 時間を計測できないため, もっと大きな数字から始めます. 測定結果は次のとおりです. 見やすいように, 不要な結果を削除しています.
ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (search-for-primes 100000000000 100000000100)
100000000003 *** 45
100000000019 *** 29
100000000057 *** 28
100000000063 *** 29
end
> (search-for-primes 1000000000000 1000000000100)
1000000000039 *** 97
1000000000061 *** 90
1000000000063 *** 90
end
> (search-for-primes 10000000000000 10000000000100)
10000000000037 *** 293
10000000000051 *** 302
10000000000099 *** 307
end
> (search-for-primes 100000000000000 100000000000100)
100000000000031 *** 939
100000000000067 *** 1397
100000000000097 *** 949
end
> (search-for-primes 1000000000000000 1000000000000200)
1000000000000037 *** 2891
1000000000000091 *** 2928
1000000000000159 *** 2919
end
> (search-for-primes 10000000000000000 10000000000000100)
10000000000000061 *** 9149
10000000000000069 *** 9166
10000000000000079 *** 9174
end
結果を見ると, nを1桁増やすと時間は約3倍になっており、 nを2桁増やすと時間は約10倍になっています. このことから, プログラムが計算に必要なステップ数に比例した時間で走るという考えに合っていると言えます.
問題文には, search-for-primesは「指定範囲の連続する奇整数について」とありますが, この条件を満たしていません. 実害はないのでこのままとします.

2017/11/12

[SICP] 問題 1.21 : 最小除数の探索

smallest-divisor手続きを使い, 次の数の最小除数を見つけよ: 199, 1999, 19999.
テキストの手続きを入力します.
(define (square x) (* x x))

(define (divides? a b)
  (= (remainder b a) 0))

(define (find-divisor n test-divisor)
  (cond ((> (square test-divisor) n) n)
        ((divides? test-divisor n) test-divisor)
        (else (find-divisor n (+ test-divisor 1)))))

(define (smallest-divisor n)
  (find-divisor n 2))

(define (prime? n)
  (= n (smallest-divisor n)))
評価してみます.
ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (smallest-divisor 199)
199
> (smallest-divisor 1999)
1999
> (smallest-divisor 19999)
7
> 

2017/11/10

[SICP] 問題 1.19 : Fibonacci数を対数的ステップ数で計算する

Fibonacci数を対数的ステップ数で計算するうまいアルゴリズムがある. 1.2.2 節のfib-iterプロセスで使った状態変数aとbの変換: a ← a + bとb ← aに注意しよう. この変換をTと呼ぶ. 1と0から始め, Tを繰り返してn回作用させると, Fib(n + 1)とFib(n)の対が出来る. いいかえれば, Fibonacci数は対(1, 0)にTn, つまり変換Tのn乗を作用させると得られる. さて, Tpqは対(a, b)をa ← bq + aq + apとb ← bp + aqに従って変換するものとし, Tを変換族Tpqのp = 0とq = 1の特例としよう. 変換Tpqを二回使うとその効果は同じ形式の変換Tp'q'を一回使ったのと同じになることを示し, p'とq'をp, qを使って表せ. これで変換を二乗する方法が得られ, fast-exptのように逐次平方を使い, Tnを計算することが出来る. これらをすべてまとめ, 対数的ステップ数の以下の手続きを完成せよ.
(define (fib n)
  (fib-iter 1 0 0 1 n))

(define (fib-iter a b p q count)
  (cond ((= count 0) b)
        ((even? count)
         (fib-iter a
                   b
                   <??>      ; p'を計算
                   <??>      ; q'を計算
                   (/ count 2)))
        (else (fib-iter (+ (* b q) (* a q) (* a p))
                        (+ (* b p) (* a q))
                        p
                        q
                        (- count 1)))))
変換Tpqは次のような行列として表現できます.
Tpq= [ p+qq ]
qp
変換Tp'q'はTpq2となります. Tpq2を計算し整理すると次のようになります.
Tp'q'= [ (p2+q2) + (2pq+q2)2pq+q2 ]
2pq+q2 p2+q2
変換Tp'q'とTpqを見比べると, p'とq'はpとqを用いて次のように表せます.
p'=p2+q2
q'=2pq+q2
上で求めたp'とq'を使って手続きを定義すると次のようになります. 計算結果を検証するために, p.22の手続きもfib0として定義します.
(require (lib "racket/trace.ss"))

(define (fib0 n)
  (define (iter a b count)
    (if (= count 0)
        b
        (iter (+ a b) a (- count 1))))
  (iter 1 0 n))

(define (fib-iter a b p q count)
  (cond ((= count 0) b)
        ((even? count)
         (fib-iter a
                   b
                   (+ (* p p) (* q q))
                   (+ (* 2 p q) (* q q))
                   (/ count 2)))
        (else
         (fib-iter (+ (* b q) (* a q) (* a p))
                   (+ (* b p) (* a q))
                   p
                   q
                   (- count 1)))))

(define (fib n)
  (fib-iter 1 0 0 1 n))

;;(trace fib-iter)
計算結果が正しいことを確認します.
ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (map fib0 '(0 1 2 3 4 5 6 7))
(0 1 1 2 3 5 8 13)
> (map fib '(0 1 2 3 4 5 6 7))
(0 1 1 2 3 5 8 13)
> 
トレースしてみます.
ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (fib 30)
>(fib-iter 1 0 0 1 30)
>(fib-iter 1 0 1 1 15)
>(fib-iter 2 1 1 1 14)
>(fib-iter 2 1 2 3 7)
>(fib-iter 13 8 2 3 6)
>(fib-iter 13 8 13 21 3)
>(fib-iter 610 377 13 21 2)
>(fib-iter 610 377 610 987 1)
>(fib-iter 1346269 832040 610 987 0)
<832040
832040
> 
トレースの結果を見ると, 対数的ステップ数で計算していることが分かります.

2017/11/09

[読書] 東京と神戸に核ミサイルが落ちたとき所沢と大阪はどうなる(兵頭二十八 著)


東京と神戸に核ミサイルが落ちたとき所沢と大阪はどうなる (講談社+α新書)

物騒なタイトルの本をジュンク堂で見かけました。「何故、神戸?」と疑問に思いましたが、三菱と川崎が潜水艦を造っていることに気づきました。自分の予想が当たっているか確かめるために購入して読んでみました。

予想は当たっていましたが、次の問題は、「何時、核攻撃を受けるのか?」になります。

この本では、中国共産党が日本を核攻撃するときのシナリオを検討しています。核攻撃する時期として2つを想定しています。
  1. 中国共産党がアメリカとの戦争を開始した直後
  2. 中国共産党が戦争に負け、崩壊する直前
戦争開始直後では、在日米軍の最大の基地である横須賀が目標となります。目的は、中国近辺での米軍の活動を抑えるためです。

中国共産党が崩壊する直前では、戦後、東アジア地域における日本の発言力を抑えるため、目標は、首都東京、核関連施設のある茨城県東海村と青森県六ケ所村、潜水艦を建造している神戸、などとなります。

また、今話題の北朝鮮が狙うとしたら、小牧基地、千歳基地が目標になると理由と共に記述しています。

この本では、核攻撃を受けた後に予想される被害についても書いてあります。私は神戸より西に住んでいるので直接的な被害は無さそうです。ただし、通勤にJRを使用しているので自宅まで歩いて帰ることになりそうです。また、社会は大混乱に陥るでしょうから、2週間程度は補給無しで生きるための準備が必要になります。

核攻撃だけでなく、大規模な自然災害(例えば南海トラフ)が発生した時のための準備をしておくことが、今の私にできる対策です。

[SICP] 問題 1.18 : 対数的ステップ数の乗算手続き(反復的プロセス)

問題1.16, 1.17の結果を使い, 加算, 二倍, 二分による, 対数的ステップ数の, 二つの整数の乗算の反復的プロセスを生成する手続きを工夫せよ.
問題 1.16の結果を参考にして手続きを定義します. 0や負の数をかけた場合の扱いは, ループに入る前に判定してます.
(require (lib "racket/trace.ss"))

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

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

(define (fast-mult-iter a b p)
  (cond ((= b 0) p)
        ((even? b) (fast-mult-iter (double a) (halve b) p))
        (else (fast-mult-iter a (- b 1) (+ a p)))))

(define (fast-mult a b)
  (cond ((< b 0) (* -1 (fast-mult-iter a (abs b) 0)))
        ((= b 0) 0)
        (else (fast-mult-iter a b 0))))

(trace fast-mult-iter)
p+a×bを考え, pに計算の途中経過をため込むことにします. pを0から始めます. p + a×bの計算について, 奇数と偶数の場合に分けて, 計算をどのように進めるかを考えます.
1) bが奇数の場合, b = m+1とする. 次の式が成立する.
p + a×(m+1) = (p+a) + (a×m)
(p+a)が次のpとなります.
2) nが偶数の場合, b = 2mとする. 次の式が成立する.
P + a×(2m) = p + double(a)×m
double(a)が次のaとなります.
実行結果は次のとおりです.
ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (fast-mult 3 45)
>(fast-mult-iter 3 45 0)
>(fast-mult-iter 3 44 3)
>(fast-mult-iter 6 22 3)
>(fast-mult-iter 12 11 3)
>(fast-mult-iter 12 10 15)
>(fast-mult-iter 24 5 15)
>(fast-mult-iter 24 4 39)
>(fast-mult-iter 48 2 39)
>(fast-mult-iter 96 1 39)
>(fast-mult-iter 96 0 135)
<135
135
> (fast-mult 4 0)
0
> (fast-mult 5 -8)
>(fast-mult-iter 5 8 0)
>(fast-mult-iter 10 4 0)
>(fast-mult-iter 20 2 0)
>(fast-mult-iter 40 1 0)
>(fast-mult-iter 40 0 40)
<40
-40
> 
トレースの結果から, 反復的プロセスを生成していることが分かります.

2017/11/07

[SICP] 問題 1.29 : Simpsonの公式

Simpsonの公式は上に示したのより更に正確な数値積分法である.
Simpsonの公式を使えば, aからbまでの関数fの積分は, 偶整数nに対し, h = (b - a)/nまたyk = f(a + kh)として
(h/3)[y0 + 4y1 + 2y2 + 4y3 + 2y4 + ・・・ + 2yn-2 + 4yn-1 + yn]
で近似出来る. (nの増加は近似の精度を増加する.)
引数としてf, a, bとnをとり, Simpson公式を使って計算した積分値を返す手続きを定義せよ. その手続きを使って(n = 100とn = 1000で)0から1までcubeを積分し, またその結果を上のintegral手続きの結果と比較せよ.
Simpsonの公式の大カッコの中身をじっと見つめると, 次のようなパターンが見えてきます.
(y0 + 4y1 + y2) + (y2 + 4y3 + y4) + ・・・ + (yn-4 + 4yn-3 + yn-2) + (yn-2 + 4yn-1 + yn)
y2に注目してください. 公式では係数が2となっていますが, それを2つに分けています.
つまり, 大カッコの中は3つの項をひとまとめにするとスッキリと書けます. 手続きを定義します.
(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
> 
計算結果を見ると精度が全く異なっています.

ヒューマン・リソース・マシーン 入社41年目−並べ替えよ

目次 1)課題 0を終端とした文字列がいくつか流れてきます。各文字列に対してソート(並べ替え)を行い、小さい順(昇順)に右側へ運んでください。 2)状況の確認 この問題では, 予めコードが入っています. このコードを実行して, 何をするコードなのか確かめます.  左のコンベアから...