an:sicp: 1章

ラベル 1章の投稿を表示しています。すべての投稿を表示

2010-12-01

2010-11-10

問題1.26

オリジナルのexpmodにおいて、(evenp exp)が成り立つexpmodの呼び出しの数をaとすると、

$2^a \le n < 2^{a + 1}$

なので、

$a \le \log_2n < a + 1$

(oddp exp)が成り立つexpmodの呼び出しの数をbとすると、

$0 \le b \le a + 1$

であるから、

$0 \le b \le \log_2n + 1$

(zerop exp)が成り立つexpmodの呼び出しは一回限りなので、expmodの呼び出しの合計は、

$\log_2n < a + b + 1 \le 2\log_2n + 2$

であり、プロセスはΘ(log n)。

これに対し、Louisのexpmodでは、まず、(evenp exp)が成り立つexpmodの呼び出しの数は、

$\sum_{i=0}^{a - 1}2^i$

(evenp exp)が成り立つごとにプロセスが二方向に分岐するため。

次に、(oddp exp)が成り立つexpmodの呼び出しの数はn。プロセスの過程でexpが1/2になるごとに、(oddp exp)が成り立つ呼び出しの数が2倍になり、結果的にnと等しくなる。

最後に、(zerop exp)が成り立つexpmodの呼び出しの回数は2^a。本来プロセスの中で一回だけのものが、最終的に(evenp exp)が成り立つexpmodの呼び出しの数だけ2倍になる。

これより、Louisのexpmodでのexpmodの呼び出しの合計は、

$\sum_{i=0}^{a - 1}2^i + n + 2^a$

となり、nに比例し、プロセスはΘ(n)。

問題1.25

結果の同一性に関しては正しいが、計算の効率で劣る。

例えば、2¹⁰⁰を3で割った剰余を計算する場合、Alyssaのexpmodでは、

> (expmod-alyssa 2 100 3)
0> Calling (REMAINDER 1267650600228229401496703205376 3) 
<0 REMAINDER returned 1
1
> (time (iter (for i below 10000000) (expmod-alyssa 2 100 3)))
(ITER (FOR I BELOW 10000000) (EXPMOD-ALYSSA 2 100 3)) took 11,578 milliseconds (11.578 seconds) to run 
                    with 2 available CPU cores.
During that period, 11,406 milliseconds (11.406 seconds) were spent in user mode
                    78 milliseconds (0.078 seconds) were spent in system mode
860 milliseconds (0.860 seconds) was spent in GC.
 800,000,032 bytes of memory allocated.
NIL

となり、remainderで2¹⁰⁰ = 1267650600228229401496703205376に対する計算をしなければならない。対して、本来のexpmodでは、

> (expmod 2 100 3)
0> Calling (REMAINDER 2 3) 
<0 REMAINDER returned 2
0> Calling (REMAINDER 4 3) 
<0 REMAINDER returned 1
0> Calling (REMAINDER 2 3) 
<0 REMAINDER returned 2
0> Calling (REMAINDER 4 3) 
<0 REMAINDER returned 1
0> Calling (REMAINDER 1 3) 
<0 REMAINDER returned 1
0> Calling (REMAINDER 1 3) 
<0 REMAINDER returned 1
0> Calling (REMAINDER 2 3) 
<0 REMAINDER returned 2
0> Calling (REMAINDER 4 3) 
<0 REMAINDER returned 1
0> Calling (REMAINDER 1 3) 
<0 REMAINDER returned 1
1
> (time (iter (for i below 10000000) (expmod 2 100 3)))
(ITER (FOR I BELOW 10000000) (EXPMOD 2 100 3)) took 10,000 milliseconds (10.000 seconds) to run 
                    with 2 available CPU cores.
During that period, 10,000 milliseconds (10.000 seconds) were spent in user mode
                    0 milliseconds (0.000 seconds) were spent in system mode
 32 bytes of memory allocated.
NIL

といったように、最大でも4に対する計算で済む。実（被除数）に対するremainderのステップ数の増加の程度が大きいほど、Alyssaのexpmodは速度において大きな問題を抱えることになる。また、一般的に計算機では大きい数の計算は遅いため、その意味でも不利になる。

例外として、法が極端に大きい場合、大きな数の計算の遅さから、Alyssaのexpmodの方が効率が良くなる場合がある。

> (time (iter (for i below 10000000) (expmod 2 2 (expt 2 100))))
(ITER (FOR I BELOW 10000000) (EXPMOD 2 2 (EXPT 2 100))) took 1,781 milliseconds (1.781 seconds) to run 
                    with 2 available CPU cores.
During that period, 1,781 milliseconds (1.781 seconds) were spent in user mode
                    0 milliseconds (0.000 seconds) were spent in system mode
 32 bytes of memory allocated.
NIL
> (time (iter (for i below 10000000) (expmod-alyssa 2 2 (expt 2 100))))
(ITER (FOR I BELOW 10000000) (EXPMOD-ALYSSA 2 2 (EXPT 2 100))) took 1,297 milliseconds (1.297 seconds) to run 
                    with 2 available CPU cores.
During that period, 1,297 milliseconds (1.297 seconds) were spent in user mode
                    0 milliseconds (0.000 seconds) were spent in system mode
 32 bytes of memory allocated.
NIL

2010-11-03

問題1.24

(eval-when (:compile-toplevel :load-toplevel :execute)
  (ql:quickload "cffi"))

(cffi:defcvar "errno" :int)

(cffi:defcfun "strerror" :string
  "Map the error number to a locale-dependent error message string."
  (errnum :int))

(cffi:define-foreign-library librt
  (t (:default "librt")))

(cffi:use-foreign-library librt)

(cffi:defctype clockid-t :int)
(cffi:defctype time-t :long)

(cffi:defcstruct timespec
  (tv_sec time-t)
  (tv_nsec :long))

(cffi:defcfun (%clock-gettime "clock_gettime") :int
  (clock-id clockid-t)
  (tp (:pointer timespec)))

(defconstant +clock-realtime+ 0)
(defconstant +clock-monotonic+ 1)
(defconstant +clock-process-cputime+ 2)
(defconstant +clock-thread-cputime+ 3)

(defun clock-gettime (clock-id)
  "Return the current value for the specified clock."
  (let ((id (ccase clock-id
              ((:realtime) +clock-realtime+)
              ((:monotonic) +clock-monotonic+)
              ((:process-cputime) +clock-thread-cputime+)
              ((:thread-cputime) +clock-thread-cputime+))))
    (cffi:with-foreign-object (tp 'timespec)
      (unless (zerop (%clock-gettime id tp))
        (error "Foreign function clock-gettime faild: ~a" (strerror *errno*)))
      (values (cffi:foreign-slot-value tp 'timespec 'tv_sec)
              (cffi:foreign-slot-value tp 'timespec 'tv_nsec)))))

(defmacro runtime ()
  `(multiple-value-bind (sec nsec)
       (clock-gettime :process-cputime)
     (+ (* sec (expt 10 9)) nsec)))

(defun expmod (base exp m)
  (cond ((zerop exp) 1)
        ((evenp exp)
         (mod (expt (expmod base (/ exp 2) m) 2) m))
        (t (mod (* base (expmod base (1- exp) m)) m))))

(defun fermat-test (n)
  (flet ((try-it (a)
           (= (expmod a n n) a)))
    (try-it (1+ (random (1- n))))))

(defun fast-primep (n times)
  (cond ((zerop times) t)
        ((fermat-test n) (fast-primep n (1- times)))
        (t nil)))

(defun report-prime (elapsed-time)
  (princ " *** ")
  (princ elapsed-time))

(defun start-prime-test (n start-time)
  (when (fast-primep n 1)
    (report-prime (- (runtime) start-time))))

(defun timed-prime-test (n)
  (terpri)
  (princ n)
  (start-prime-test n (runtime)))

(defun search-for-primes (start end)
  (when (oddp start)
    (timed-prime-test start))
  (when (< start end)
    (search-for-primes (1+ start) end)))

Fermatテストはθ(log n)でステップが増加するが、

loga(1000000) = loga(1000^2) = 2loga(1000)

なので、1000000付近の素数をテストするときは、1000付近のときの二倍前後のステップ数と予想され、時間もステップ数に比例すると思われる。

> (search-for-primes 1000 1050)

1001 *** 10851
1003
1005
1007
1009 *** 2683
1011
1013 *** 2683
1015
1017
1019 *** 2695
1021 *** 2761
1023
1025
1027
1029
1031 *** 2635
1033 *** 2635
1035
1037
1039 *** 2581
1041
1043
1045
1047
1049 *** 2557
NIL
> (search-for-primes 1000000 1000050)

1000001
1000003 *** 7510
1000005
1000007
1000009
1000011
1000013
1000015
1000017
1000019
1000021
1000023
1000025
1000027
1000029
1000031
1000033 *** 7707
1000035
1000037 *** 7264
1000039 *** 7552
1000041
1000043
1000045
1000047
1000049
NIL

おおむね二倍から三倍の間の値となるため、予想通りと言っても差し支えない範囲と思われる。

問題1.23

(defun next-divisor (n)
  (if (= n 2)
      3
      (+ n 2)))

(defun find-divisor (n test-divisor)
  (cond ((> (expt test-divisor 2) n) n)
        ((dividesp test-divisor n) test-divisor)
        (t (find-divisor n (next-divisor test-divisor)))))

問題1.22のバージョンでの計測結果は、

> (mapc #'(lambda (n)
            (timed-prime-test n))
        '(1009 1013 1019
          10007 10009 10037
          100003 100019 100043
          1000003 1000033 1000037))

1009 *** 18605
1013 *** 10713
1019 *** 10994
10007 *** 25079
10009 *** 25055
10037 *** 25061
100003 *** 73716
100019 *** 74464
100043 *** 74554
1000003 *** 225458
1000033 *** 229998
1000037 *** 225387
(1009 1013 1019 10007 10009 10037 100003 100019 100043 1000003 1000033 1000037)

nextを使ったバージョンの計測結果は、

> (mapc #'(lambda (n)
            (timed-prime-test n))
        '(1009 1013 1019
          10007 10009 10037
          100003 100019 100043
          1000003 1000033 1000037))

1009 *** 12869
1013 *** 6695
1019 *** 6336
10007 *** 13707
10009 *** 14019
10037 *** 13785
100003 *** 39834
100019 *** 38745
100043 *** 38475
1000003 *** 112975
1000033 *** 112921
1000037 *** 112813
(1009 1013 1019 10007 10009 10037 100003 100019 100043 1000003 1000033 1000037)

となり、速度の比は平均で約1.8になった。

ステップが半分になった代わりに、ステップごとにnextが呼ばれ、その中でnと2の比較が行われるようになったため、単純に速度の比は2にならない。

2010-10-30

問題1.22

(eval-when (:compile-toplevel :load-toplevel :execute)
  (ql:quickload "cffi"))

(cffi:defcvar "errno" :int)

(cffi:defcfun "strerror" :string
  "Map the error number to a locale-dependent error message string."
  (errnum :int))

(cffi:define-foreign-library librt
  (t (:default "librt")))

(cffi:use-foreign-library librt)

(cffi:defctype clockid-t :int)
(cffi:defctype time-t :long)

(cffi:defcstruct timespec
  (tv_sec time-t)
  (tv_nsec :long))

(cffi:defcfun (%clock-gettime "clock_gettime") :int
  (clock-id clockid-t)
  (tp (:pointer timespec)))

(defconstant +clock-realtime+ 0)
(defconstant +clock-monotonic+ 1)
(defconstant +clock-process-cputime+ 2)
(defconstant +clock-thread-cputime+ 3)

(defun clock-gettime (clock-id)
  "Return the current value for the specified clock."
  (let ((id (ccase clock-id
              ((:realtime) +clock-realtime+)
              ((:monotonic) +clock-monotonic+)
              ((:process-cputime) +clock-thread-cputime+)
              ((:thread-cputime) +clock-thread-cputime+))))
    (cffi:with-foreign-object (tp 'timespec)
      (unless (zerop (%clock-gettime id tp))
        (error "Foreign function clock-gettime faild: ~a" (strerror *errno*)))
      (values (cffi:foreign-slot-value tp 'timespec 'tv_sec)
              (cffi:foreign-slot-value tp 'timespec 'tv_nsec)))))

(defmacro runtime ()
  `(multiple-value-bind (sec nsec)
       (clock-gettime :process-cputime)
     (+ (* sec (expt 10 9)) nsec)))

(defun dividesp (a b)
  (zerop (mod b a)))

(defun find-divisor (n test-divisor)
  (cond ((> (expt test-divisor 2) n) n)
        ((dividesp test-divisor n) test-divisor)
        (t (find-divisor n (1+ test-divisor)))))

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

(defun primep (n)
  (= n (smallest-divisor n)))

(defun report-prime (elapsed-time)
  (princ " *** ")
  (princ elapsed-time))

(defun start-prime-test (n start-time)
  (when (primep n)
    (report-prime (- (runtime) start-time))))

(defun timed-prime-test (n)
  (terpri)
  (princ n)
  (start-prime-test n (runtime)))

(defun search-for-primes (start end)
  (when (oddp start)
    (timed-prime-test start))
  (when (< start end)
    (search-for-primes (1+ start) end)))

runtimeがないので、clock_gettimeを使って定義する。精度はナノ秒。また、定義が自明なため、奇数の判別には組み込みのoddpを使っている。

> (search-for-primes 1000 1050)

1001
1003
1005
1007
1009 *** 9114
1011
1013 *** 8887
1015
1017
1019 *** 9054
1021 *** 9078
1023
1025
1027
1029
1031 *** 9275
1033 *** 9431
1035
1037
1039 *** 9162
1041
1043
1045
1047
1049 *** 9066
NIL
> (search-for-primes 10000 10050)

10001
10003
10005
10007 *** 23444
10009 *** 23803
10011
10013
10015
10017
10019
10021
10023
10025
10027
10029
10031
10033
10035
10037 *** 23983
10039 *** 23486
10041
10043
10045
10047
10049
NIL
> (search-for-primes 100000 100050)

100001
100003 *** 73907
100005
100007
100009
100011
100013
100015
100017
100019 *** 73650
100021
100023
100025
100027
100029
100031
100033
100035
100037
100039
100041
100043 *** 73279
100045
100047
100049 *** 73584
NIL
> (search-for-primes 1000000 1000050)

1000001
1000003 *** 217447
1000005
1000007
1000009
1000011
1000013
1000015
1000017
1000019
1000021
1000023
1000025
1000027
1000029
1000031
1000033 *** 216698
1000035
1000037 *** 216829
1000039 *** 227579
1000041
1000043
1000045
1000047
1000049
NIL

おおむねそれらしい値になっている。純粋な処理時間以外にも、手続きを呼び出すオーバーヘッドなども含まれるので、多少の誤差は許容範囲内と思われる。

問題1.21

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

(defun find-divisor (n test-divisor)
  (cond ((> (expt test-divisor 2) n) n)
        ((dividesp test-divisor n) test-divisor)
        (t (find-divisor n (1+ test-divisor)))))

(defun dividesp (a b)
  (zerop (mod b a)))

(smallest-divisor 199)
;; => 199
(smallest-divisor 1999)
;; => 1999
(smallest-divisor 19999)
;; => 7

2010-10-27

2010-10-23

問題1.19

定義により、(Tpq)^2は

a←(bp + aq)q + (bq + aq + ap)q + (bq + aq + ap)p
b←(bp + aq)p + (bq + aq + ap)q

であり、これがTp'q'と等しいのであれば、

bq' + aq' + ap' = (bp + aq)q + (bq + aq + ap)q + (bq + aq + ap)p
bp' + aq' = (bp + aq)p + (bq + aq + ap)q

が成り立つ。この連立方程式を解くと、

aq' = (bp + aq)p + (bq + aq + ap)q - bp'
 q' = (bp^2 + bq^2 - bp') / a + q^2 + 2pq

b{(bp^2 + bq^2 - bp') / a + q^2 + 2pq} + a{(bp^2 + bq^2 - bp') / a + q^2 + 2pq} + ap'
  = (bp + aq)q + (bq + aq + ap)q + (bq + aq + ap)p
(a^2 - b^2 - ab)p' = (a^2 - b^2 - ab)p^2 + (a^2 - b^2 - ab)q^2
                p' = p^2 + q^2

b(p^2 + q^2) + aq' = (bp + aq)p + (bq + aq + ap)q
 bp^2 + bq^2 + aq' = bp^2 + apq + bq^2 + aq^2 + apq
                   = bp^2 + aq^2 + bq^2 + 2apq
               aq' = aq^2 + 2apq
                q' = q^2 + 2pq

となり、p' = p^2 + q^2、q' = q^2 + 2pqのとき、(Tpq)^2 = Tp'q'は成立する。

手続きfibの定義は以下の通り。

(5am:test exercise-1.19
  (5am:is (= 0 (fib 0)))
  (5am:is (= 1 (fib 1)))
  (5am:is (= 5 (fib 5)))
  (5am:is (= (fib 10)
             (+ (fib 8) (fib 9)))))

(defun square (x)
  (* x x))

(defun fib-iter (a b p q count)
  (cond ((zerop count) b)
        ((evenp count)
         (fib-iter a
                   b
                   (+ (square p) (square q))
                   (+ (square q) (* 2 p q))
                   (/ count 2)))
        (t (fib-iter (+ (* b q) (* a q) (* a p))
                     (+ (* b p) (* a q))
                     p
                     q
                     (1- count)))))

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

2010-10-21

問題1.18

(5am:test exercise-1.18
  (5am:is (= 0 (multi 1 0)))
  (5am:is (= 6 (multi 2 3)))
  (5am:is (= (multi 4 5)
             (multi 5 4))))

(defun double (x)
  (* x 2))

(defun halve (x)
  (/ x 2))

(defun multi-iter (a b r)
  (cond ((zerop b) r)
        ((evenp b) (multi-iter (double a) (halve b) r))
        (t (multi-iter a (1- b) (+ r a)))))

(defun multi (a b)
  (multi-iter a b 0))

問題1.17

(5am:test exercise-1.17
  (5am:is (= 0 (multi 1 0)))
  (5am:is (= 6 (multi 2 3)))
  (5am:is (= (multi 4 5)
             (multi 5 4))))

(defun double (x)
  (* x 2))

(defun halve (x)
  (/ x 2))

(defun multi (a b)
  (cond ((zerop b) 0)
        ((evenp b) (double (multi a (halve b))))
        (t (+ a (multi a (1- b))))))

2010-10-18

問題1.16

(5am:test exercise-1.16
  (5am:is (= 1 (fast-expt 2 0)))
  (5am:is (= 2 (fast-expt 2 1)))
  (5am:is (= 8 (fast-expt 2 3)))
  (5am:is (= (* 2 (fast-expt 2 9))
             (fast-expt 2 10))))

(defun square (x)
  (* x x))

(defun fast-expt-iter (b n a)
  (cond ((zerop n) a)
        ((evenp n) (fast-expt-iter (square b) (/ n 2) a))
        (t (fast-expt-iter b (1- n) (* a b)))))

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

2010-10-17

問題1.15

aについて。(sine 12.15)でpが呼ばれるのは5回。以下traceの出力。

> (sine 12.15)
0> Calling (P 0.049999997) 
<0 P returned 0.1495
0> Calling (P 0.1495) 
<0 P returned 0.43513453
0> Calling (P 0.43513453) 
<0 P returned 0.9758465
0> Calling (P 0.9758465) 
<0 P returned -0.7895632
0> Calling (P -0.7895632) 
<0 P returned -0.39980328
-0.39980328

bについて。ステップ数はsineが呼ばれる回数に比例する。ここで、sineが呼ばれる回数をnとすると、

a / 3^n ≦ 0.1 < a / 3^(n - 1)

という関係が成り立つ。これは、

a / 3^n ≦ 0.1 < a / 3^(n - 1)
a / 0.1 ≦ 3^n < 3 * a / 0.1
log3(a / 0.1) ≦ n < log3(3 * a / 0.1)
log3(a) - log3(0.1) ≦ n < log3(a) - log3(0.1) + 1

よって、nは3を底とするaの対数関数であることが分かる。sineではスペースとステップ数は単純に比例するため、増加の程度は共にΘ(log a)。

問題1.14

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

プロセスが使うスペースについて。木構造再帰的プロセスに必要なスペースは、木構造の最大の深さに比例する。1セントだけで両替した場合が最も木構造が深く、両替する金額に比例するため、増加の程度はΘ(n)。

プロセスのステップ数について。両替する金額をnとする。1セントだけで両替した場合、nに比例する深さのステップの枝が一本、1セントと5セントで両替した場合、1セントでの両替に分岐する枝が、nに比例する本数存在する。1セントと5セントと10セントの場合、1セントと5セントの両替に分岐する枝が同様となる。このように、ステップ数は、nを底、硬貨の種類を指数とする冪乗に比例する。つまり、五種類の硬貨での増加の程度はΘ(n^5)。

2010-10-14

問題1.13

ψ = (1 - √5) / 2
f(n) = (φ^n - ψ^n) / √5

と仮定する。

f(0) = 0        ― (1)
f(1) = 1        ― (2)

となり、

f(n - 1) + f(n - 2) = (φ^(n - 1) - ψ^(n - 1)) / √5 + (φ^(n - 2) - ψ^(n - 2)) / √5
                    = (φ^(n - 2) * (φ + 1) - ψ^(n - 2) * (ψ + 1)) / √5     ― (3)

となる。ここで、

φ^2 = ((1 + √5) / 2)^2 = (3 + √5) / 2 = (1 + √5) / 2 + 1
    = φ + 1
ψ^2 = ((1 - √5) / 2)^2 = (3 - √5) / 2 = (1 - √5) / 2 + 1
    = ψ + 1

であるので、それぞれ(3)に代入すると、

f(n - 1) + f(n - 2) = (φ^(n - 2) * φ^2 - ψ^(n - 2) * ψ^2) / √5
                    = (φ^n - ψ^n) / √5
                    = f(n)      ― (4)

Fib(n)の定義は

         0                              n = 0
Fib(n) = 1                              n = 1
         Fib(n - 1) + Fib(n - 2)        otherwise.

であるから、(1)、(2)、(4)より

Fib(n) = f(n)   ― (5)

Fib(n)がφ^n / √5に最も近い整数であるなら、

Fib(n) - 0.5 < φ^n / √5 < Fib(n) + 0.5        ― (6)

が成り立つ。(5)より、

               f(n) - 0.5 < φ^n / √5 < f(n) + 0.5
   (φ^n - ψ^n) / √5 - 0.5 < φ^n / √5 < (φ^n - ψ^n) / √5 + 0.5
φ^n / √5 - ψ^n / √5 - 0.5 < φ^n / √5 < φ^n / √5 - ψ^n / √5 + 0.5      ― (7)

となる。ここで、

ψ^n / √5 = ((1 - √5) / 2)^n / √5

だが、

2 < √5 < 3

であるため、

|ψ^n / √5| < 0.5

は自明であり、(7)および(6)が成り立つ。よって、Fib(n)はφ^n / √5に最も近い整数である。

（証明終）

2010-10-13

問題1.12

(5am:test exercise-1.12
  (5am:is (= 1 (pascals-triangle 1 1)))
  (5am:is (= 6 (pascals-triangle 3 5)))
  (5am:is (= (+ (pascals-triangle 4 9)
                (pascals-triangle 5 9))
             (pascals-triangle 5 10)))
  (5am:is-false (pascals-triangle 3 2))) ; invalid coord

(defun pascals-triangle (x y)
  (cond ((> x y) nil)
        ((or (= x 1) (= x y)) 1)
        (t (+ (pascals-triangle (1- x) (1- y))
              (pascals-triangle x (1- y))))))

問題1.11

(5am:test exercise-1.11
  (5am:is (zerop (f 0)))
  (5am:is (= 2 (f 2)))
  (5am:is (= (+ (f 2)
                (* 2 (f 1))
                (* 3 (f 0)))
             (f 3))))

;; recursive
(defun f (n)
  (if (< n 3)
      n
      (+ (f (1- n))
         (* 2 (f (- n 2)))
         (* 3 (f (- n 3))))))

;; (iter (for i below 10) (collect (f i)))
;; => (0 1 2 4 11 25 59 142 335 796)

;; f(n) = f(n - 1) + 2f(n - 2) + 3f(n - 3)
;; f(3) = f(2) + 2f(1) + 3f(0)

;; 0   1   2
;; a + b + c

;; a <- b
;; b <- c
;; c <- 3a + 2b + c

;; iterative
(defun f (n)
  (if (< n 3)
      n
      (f-iter 0 1 2 n)))

(defun f-iter (a b c n)
  (if (< n 3)
      c
      (f-iter b c (+ (* 3 a) (* 2 b) c) (1- n))))

2010-08-08

問題1.10

(A 1 10)
(A 0 (A 1 9))
(A 0 (A 0 (A 1 8)))
(A 0 (A 0 (A 0 (A 1 7))))
(A 0 (A 0 (A 0 (A 0 (A 1 6)))))
(A 0 (A 0 (A 0 (A 0 (A 0 (A 1 5))))))
(A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 1 4)))))))
(A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 1 3))))))))
(A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 1 2)))))))))
(A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 1 1))))))))))
(A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 2)))))))))
(A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 4))))))))
(A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 8)))))))
(A 0 (A 0 (A 0 (A 0 (A 0 (A 0 16))))))
(A 0 (A 0 (A 0 (A 0 (A 0 32)))))
(A 0 (A 0 (A 0 (A 0 64))))
(A 0 (A 0 (A 0 128)))
(A 0 (A 0 256))
(A 0 512)
1024

(A 2 4)
(A 1 (A 2 3))
(A 1 (A 1 (A 2 2)))
(A 1 (A 1 (A 1 (A 2 1))))
(A 1 (A 1 (A 1 2)))
(A 1 (A 1 (A 0 (A 1 1))))
(A 1 (A 1 (A 0 2)))
(A 1 (A 1 4))
(A 1 (A 0 (A 1 3)))
(A 1 (A 0 (A 0 (A 1 2))))
(A 1 (A 0 (A 0 (A 0 (A 1 1)))))
(A 1 (A 0 (A 0 (A 0 2))))
(A 1 (A 0 (A 0 4)))
(A 1 (A 0 8))
(A 1 16) ; 2^16
; 中略
65536

(A 3 3)
(A 2 (A 3 2))
(A 2 (A 2 (A 3 1)))
(A 2 (A 2 2))
(A 2 (A 1 (A 2 1)))
(A 2 (A 1 2))
(A 2 (A 0 (A 1 1)))
(A 2 (A 0 2))
(A 2 4)
65536

fは2n、gは2ⁿ、hは2↑↑n。

2010-08-07

問題1.9

(defun + (a b)
  (if (= a 0)
      b
      (inc (+ (dec a) b))))

(+ 4 5)
(inc (+ 3 5))
(inc (inc (+ 2 5)))
(inc (inc (inc (+ 1 5))))
(inc (inc (inc (inc (+ 0 5)))))
(inc (inc (inc (inc 5))))
(inc (inc (inc 6)))
(inc (inc 7))
(inc 8)
9

となり、再帰的プロセス。

(defun + (a b)
  (if (= a 0)
      b
      (+ (dec a) (inc b))))

(+ 4 5)
(+ 3 6)
(+ 2 7)
(+ 1 8)
(+ 0 9)
9

となり、反復的プロセス。

問題1.8

(defun good-enough? (guess x)
  (< (abs (- (expt guess 3) x)) 0.001))

(defun improve (guess x)
  (/ (+ (/ x (expt guess 2)) (* guess 2)) 3))

(defun cbrt-iter (guess x)
  (if (good-enough? guess x)
      guess
      (cbrt-iter (improve guess x) x)))

(defun cbrt (x)
  (cbrt-iter 1.0 x))

二乗や三乗を求める関数の定義は自明なので、exptで。

登録: 投稿 (Atom)

an:sicp