an:sicp: 2010-11

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にならない。