オンラインSICP読書女子会 #6

オンラインSICP読書女子会 #6 (1.2.6), くーむ氏のまとめ

前回: オンラインSICP読書女子会 #5 (1.2.4~1.2.6 最大公約数と素数判定, 練習問題1.19 .. 練習問題1.24), くーむ氏のまとめ, 自分のまとめ

ex-1.21 (復習) smallest-divisor

smallest-divisor 手続きを使って次の数値の最小の約数を求めよ。199, 1999, 19999.

(define (smallest-divisor n)
	(find-divisor n 2))
(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 (divides? a b)
	(= (remainder b a) 0))
(define (square a) (* a a))

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

おまけ:

; gosh> (/ 19999 7)
; 2857
; gosh> (smallest-divisor 2857)
; 2857
; gosh> (* 7 2857)
; 19999

19999 = 7 * 2857 でこの2つは素数の模様.

ex-1.22 (復習) 素数判定の処理時間, 総当りの場合

(timed-prime-test 手続きをつかって) 指定した範囲の連続した奇数について素数判定を行う手続き search-for-primes を書け。
その手続きを使って、1000, 10,000, 100,000 より大きな素数をそれぞれ 3 つ見つけよ。
判定アルゴリズムは Θ(√n) の増加オーダーを持っているので、10,000 あたりの判定には 1000 あたりの √10 倍程度の時間がかかるはずである。あなたの計測データはこれを裏付けているだろうか。
100,000 や 1,000,000 のデータは、Θ(√n) という予想はどれだけ当たっているだろうか。
あなたの結果は、演算に必要なステップ数に比例して実行時間が増えるという概念に矛盾していないだろうか。

(define (ex-1.22)
	(search-for-primes 1000)
	(search-for-primes 10000)
	(search-for-primes 100000)
)


(define (search-for-primes n)
	(define n2 (if (odd? n) n (+ n 1)))
	(search-for-primes-iter n2 3)
)
(define (search-for-primes-iter n rest)
	(if (timed-prime-test n)
		(if (= rest 1)
			(undefined) ; exit loop.
			(search-for-primes-iter (+ n 2) (- rest 1)))
		(search-for-primes-iter (+ n 2) rest)
	)
)


; 速すぎて計測できないので判定は2万回実行.
(define (prime? n)
	(prime-loop n 20000))
(define (prime-loop n loop)
	(if
		(= n (smallest-divisor n))
		(if (= loop 1)
			#t
			(prime-loop n (- loop 1)))
		#f))


; gauche には runtime がなかったので自前実装.
(use srfi-19)
(define (runtime)
	(define (to-microseconds t)
		(+
			(* (time-second t) 1000 1000)
			(/ (time-nanosecond t) 1000)))
	(to-microseconds (current-time)))


; timed-prime-test はテキストの記載から若干変更している.
(define (timed-prime-test n)
	(start-prime-test n (runtime)))
(define (start-prime-test n start-time)
	(if (prime? n)
		(report-prime n (- (runtime) start-time))
		#f))
(define (report-prime n elapsed-time)
	(display n)
	(display " *** ")
	(display elapsed-time)
	(newline)
	#t)


(define (smallest-divisor n)
	(find-divisor n 2))
(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 (divides? a b)
	(= (remainder b a) 0))
(define (square a) (* a a))

実行結果:

1009 *** 266833
1013 *** 262664
1019 *** 260001
10007 *** 853538
10009 *** 844711
10037 *** 853837
100003 *** 2657958
100019 *** 2655992
100043 *** 2670254

それぞれの中間値を代表値として比較:

start    prime   elapsed
  1,000    1013  0.262664
 10,000   10007  0.853538  (3.24 倍; sqrt(10) + 2.7%)
100,000  100003  2.657958  (3.11 倍; sqrt(10) - 1.5%) ( 10.11 倍)

おおよそ想定通りの増分となっているといえる.

ex-1.23 (復習)

smallest-divisor の最適化 (1)

(define (smallest-divisor n)
	(find-divisor 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 (next test-divisor)
	(if (= test-divisor 2)
		3
		(+ test-divisor 2)))

実行結果:

1009 *** 168980
1013 *** 171566
1019 *** 172254
10007 *** 513457
10009 *** 511730
10037 *** 525784
100003 *** 1620600
100019 *** 1615102
100043 *** 1657809

最適化結果:

start    prime   変更前    変更後    比率
  1,000    1013  0.262664  0.171566  65.3%
 10,000   10007  0.853538  0.513457  60.1%
100,000  100003  2.657958  1.620600  60.9%

60% 程度にとどまっていて半分までは達しなかった.

smallest-divisor の最適化 (2)

next 手続きにおける定数分の処理が考えられるので,
test-divisor の値で判定するのではなく,
smalltest-divisor で最初に偶数判定を行い,
next では (+ test-divisor 2) のみを行うように修正を行った:

(define (smallest-divisor n)
	(if
		(even? n)
		2
		(find-divisor n 3)))


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

実行結果:

1009 *** 137096
1013 *** 147942
1019 *** 137724
10007 *** 430376
10009 *** 428306
10037 *** 422682
100003 *** 1337870
100019 *** 1342409
100043 *** 1622971

start    prime   変更前    変更後 (2) 比率 (2)
  1,000    1013  0.262664  0.137724   52.4%
 10,000   10007  0.853538  0.428306   50.1%
100,000  100003  2.657958  1.342409   50.5%

ほぼ半分に達したといえそう.

ex-1.24 (復習) fermat-test

練習問題 1.22の timed-prime-test 手続きを fast-prime?(フェルマー法) を使うように修正し、練習問題 1.22で見つけた 12 個の素数をそれぞれテストせよ。フェルマーテストは Θ(log n) で増加するが、1,000,000 に近い素数をテストするのにかかる時間は 1000 に近い素数をテストするのに必要な時間と比べてどの程度になると予想できるか。実際とのずれがあれば、それを説明できるだろうか。

log(1,000,000) / log(1,000) = 2 であることから2倍程度と予測される.

ex-1.23 のコードを用いて 1,000,000 に近い素数を抽出:

(search-for-primes 1000000)
; 1000003 *** 4222554
; 1000033 *** 4231522
; 1000037 *** 4243597

この値織り込んで実装:

(use srfi-27) ; random relatives.
; (random-inreger n) ==> a ∈ [0,n-1]


(define (ex-1.24)
	(timed-prime-test    1009)
	(timed-prime-test    1013)
	(timed-prime-test    1019)
	(timed-prime-test   10007)
	(timed-prime-test   10009)
	(timed-prime-test   10037)
	(timed-prime-test  100003)
	(timed-prime-test  100019)
	(timed-prime-test  100043)
	(timed-prime-test 1000003)
	(timed-prime-test 1000033)
	(timed-prime-test 1000037)
)


; ちょうどいいので判定をこれまでと同数の2万回実行.
; 今回は素数に対して判定しているので必ず2万回実行し, 真が返る.
(define (prime? n)
	(fast-prime? n 20000))


(define (expmod base exp m)
	(cond ((= exp 0) 1) ; base**0 ==> 1, よって (base**0)%m = 1%m = 1
		((even? exp) (remainder (square (expmod base (/ exp 2) m)) m))
		(else (remainder (* base (expmod base (- exp 1) m)) m))))


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


(define (fast-prime? n times)
	(cond
		((= times 0) #t)
		((fermat-test n) (fast-prime? n (- times 1))) ; pass at this time.
		(else #f) ; certainly not prime.
		))


; gauche には runtime がなかったので自前実装.
(use srfi-19)
(define (runtime)
	(define (to-microseconds t)
		(+
			(* (time-second t) 1000 1000)
			(/ (time-nanosecond t) 1000)))
	(to-microseconds (current-time)))


; timed-prime-test はテキストの記載から若干変更している.
(define (timed-prime-test n)
	(start-prime-test n (runtime)))
(define (start-prime-test n start-time)
	(if (prime? n)
		(report-prime n (- (runtime) start-time))
		#f))
(define (report-prime n elapsed-time)
	(display n)
	(display " *** ")
	(display elapsed-time)
	(newline)
	#t)

実行:

(ex-1.24)
; 1009 *** 194461
; 1013 *** 207848
; 1019 *** 213064
; 10007 *** 262958
; 10009 *** 256652
; 10037 *** 261537
; 100003 *** 310040
; 100019 *** 310115
; 100043 *** 322753
; 1000003 *** 356723
; 1000033 *** 361114
; 1000037 *** 361577

比率:

;     1,013 *** 0.207848
; 1,000,033 *** 0.361114
; 比率 = 173.7%

予想より速く, 70% の追加で判定は終了した.
処理回数を 10 万回に増やしても 80% 程度であり, 2倍は必要としなかった.

ずれの理由は不明...；ω；

ex-1.25 畳み込まないシンプルな expmod とした場合

(define (expmod base exp m)
	(remainder (fast-expt base exp) m))

この場合, 一旦 base の exp 乗という値を完全に構築する必要があり,
例えば素数 1009 の場合で 3000 桁以上, 1万bit 級の数値を構築する必要が生じてしまう.
このような桁の多い数値に対する演算は処理に非常に時間がかかってしまうため望ましくない.

お手本にある expmod の場合であれば演算対象となるのは判定対象としている数値の自乗までであるため, 素数 1009 に対してせいぜい 20 bit 程度の作業領域で事足りる.

ex-1.26 展開された 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))))

fast-prime? では, (even? exp) 時の expmod の呼び出し回数を半分にすることで
$\Theta(\log_2(n))$ すなわち $\Theta(\log(n))$ の計算量を実現していた.
しかしこの記述では, expmod の呼び出しが明示的に2回行われてしまうため,
半分にしたかった糸が完全に廃棄されてしまうことになる.
さらには prime? 手続きの場合では $\sqrt{n}$ での終端条件により
$\Theta(\sqrt{n})$ の計算量となっていたが,
こちらのロジックではそれとは異なる終端条件となっている.
結果, どちらの計算量抑制のロジックも働かず, $\Theta(n)$ という計算量に増大してしまっている.

ex-1.27 Carmichael numbers

(define (ex-1.27)
	(test-carmichael-number  561)
	(test-carmichael-number 1105)
	(test-carmichael-number 1729)
	(test-carmichael-number 2465)
	(test-carmichael-number 2821)
	(test-carmichael-number 6601)
	(display "-- non-primes --")
	(newline)
	(test-carmichael-number 39)
	(test-carmichael-number 19999)
	;(test-carmichael-number (* 2821 6601))
	(display "-- primes --")
	(newline)
	(test-carmichael-number 199)
	(test-carmichael-number 1999)
	(test-carmichael-number    1009)
	(test-carmichael-number    1013)
	(test-carmichael-number    1019)
	(test-carmichael-number   10007)
	(test-carmichael-number   10009)
	(test-carmichael-number   10037)
	(test-carmichael-number  100003)
	(test-carmichael-number  100019)
	(test-carmichael-number  100043)
	;(test-carmichael-number 1000003)
	;(test-carmichael-number 1000033)
	;(test-carmichael-number 1000037)
)


(define (test-carmichael-number n)
	(display n)
	(display " ... ")
	(test-iter n 2 0 0)
	(newline))


(define (test-iter n a succ fail)
	(cond
		((= n a) (result n succ fail))
		((= (expmod a n n) a)
			(test-iter n (+ a 1) (+ succ 1) fail))
		(else
			(test-iter n (+ a 1) succ (+ fail 1)))))


(define (result n succ fail)
	(display " succ ")
	(display succ)
	(display " fail ")
	(display fail))


(define (expmod base exp m)
	(cond ((= exp 0) 1) ; base**0 ==> 1, よって (base**0)%m = 1%m = 1
		((even? exp) (remainder (square (expmod base (/ exp 2) m)) m))
		(else (remainder (* base (expmod base (- exp 1) m)) m))))


; gosh> (ex-1.27)
; 561 ...  succ 559 fail 0
; 1105 ...  succ 1103 fail 0
; 1729 ...  succ 1727 fail 0
; 2465 ...  succ 2463 fail 0
; 2821 ...  succ 2819 fail 0
; 6601 ...  succ 6599 fail 0
; -- non-primes --
; 39 ...  succ 7 fail 30
; 19999 ...  succ 47 fail 19950
; -- primes --
; 199 ...  succ 197 fail 0
; 1999 ...  succ 1997 fail 0
; 1009 ...  succ 1007 fail 0
; 1013 ...  succ 1011 fail 0
; 1019 ...  succ 1017 fail 0
; 10007 ...  succ 10005 fail 0
; 10009 ...  succ 10007 fail 0
; 10037 ...  succ 10035 fail 0
; 100003 ...  succ 100001 fail 0
; 100019 ...  succ 100017 fail 0
; 100043 ...  succ 100041 fail 0
; #<undef>

ex-1.28 Miller-Rabin test

フェルマーの小定理: n が素数で、a が n より小さい任意の正の整
数であるとき、a の n 乗は法 n に関して a と合同である

フェルマーの小定理の別の形:
n が素数で、かつ a が n 以下の任意の正の整数であれば、a の (n − 1) 乗は法 n に関して 1 と合同である

a は n より小さい正の整数の誤りと思われる:

if n is a prime number and a is any positive integer less than n, then a raised to the (n − 1)-st power is congruent to 1 modulo n.

「n が素数で, a が n より小さい任意の正の整数であるとき」まではどちらの形式でも同一。
別形式では $a^{n-1} \equiv 1 \pmod{n}$ とある.
これに $a$ を乗ずると, $a ^n \equiv a \pmod{n}$ となる.
($a < n$ であるため, $a \mod n = a$ を維持する)
これは前述の形式と同一である.
(逆向きも言えるのかまでは調べていない)

Miller-Rabin テストで数値 n の素数性をテストするには、 a < n である適当な a を選び、expmod 手続きを使って a の (n − 1) 乗の n を法とする剰余を求める。

これは別形式にある式であり, n が素数であればその値は 1 となる.

しかし、expmod 手続きの中で二乗を行うステップで、毎回 “法 n に関する自明でない 1 の平方根” を見つけたかチェックする。

$(square x)$ の結果が 1 となるなら, $x$ は 1 の平方根である.
通常の自然数等であれば 1 の平方根は 1 であるが,
例えば $\mod 8$ の元であれば $3 * 3 = 9 \equiv 1 \pmod{8}$ であり, 「3 は法 8 における 1 の平方根である」と言える.

(mod 8)
0^2 = 0      (0 は 1 の平方根? no)
1^2 = 1      (1 は 1 の平方根? YES)
2^2 = 4      (2 は 1 の平方根? no)
3^2 = 9  = 1 (3 は 1 の平方根? YES)
4^2 = 16 = 0 (4 は 1 の平方根? no)
5^2 = 25 = 1 (5 は 1 の平方根? YES)
6^2 = 36 = 4 (6 は 1 の平方根? no)
7^2 = 49 = 1 (7 は 1 の平方根? YES)
8^2 = 64 = 0 (8 は 1 の平方根? no)

$1$ と $n-1$ は, 法 n においては常に $1$ の平方根である.
$1$ については自明なので省略.
$(n-1)$ においては, $(n-1)^2 = n^2 - 2n + 1 = n(n - 2) + 1$
$n(n - 2)$ は法 n において $0$ であるため, 全体としては常に $1$ となる.

$m$ が法 n において $1$ の平方根であれば, $n-m$ もまた法 n において $1$ の平方根.
前提より $m^2 \equiv 1 \pmod{n}$
$(n-m)^2 = n(n - 2m) + m^2 \equiv m^2 = 1$

$(square x)$ は expmod 中で普通に算出して利用しているので,
そこに1つロジックはさんであげればよい.

(define (ex-1.28-1)
	(test-carmichael-number  561)
	(test-carmichael-number 1105)
	(test-carmichael-number 1729)
	(test-carmichael-number 2465)
	(test-carmichael-number 2821)
	(test-carmichael-number 6601)
	(display "-- non-primes --")
	(newline)
	(test-carmichael-number 39)
	(test-carmichael-number 19999)
	;(test-carmichael-number (* 2821 6601))
	(display "-- primes --")
	(newline)
	(test-carmichael-number 199)
	(test-carmichael-number 1999)
	(test-carmichael-number    1009)
	(test-carmichael-number    1013)
	(test-carmichael-number    1019)
	(test-carmichael-number   10007)
	(test-carmichael-number   10009)
	(test-carmichael-number   10037)
	(test-carmichael-number  100003)
	(test-carmichael-number  100019)
	(test-carmichael-number  100043)
	;(test-carmichael-number 1000003)
	;(test-carmichael-number 1000033)
	;(test-carmichael-number 1000037)
)


(define (test-carmichael-number n)
	(display n)
	(display " ... ")
	(test-iter n 2 0 0)
	(newline))


(define (test-iter n a succ fail)
	(cond
		((= n a) (result n succ fail))
		((= (expmod a n n) a)
			(test-iter n (+ a 1) (+ succ 1) fail))
		(else
			(test-iter n (+ a 1) succ (+ fail 1)))))


(define (result n succ fail)
	(display " succ ")
	(display succ)
	(display " fail ")
	(display fail))


(define (expmod base exp m)
	(cond ((= exp 0) 1) ; base**0 ==> 1, よって (base**0)%m = 1%m = 1
		((even? exp) (expmod-2 (expmod base (/ exp 2) m) m))
		(else (remainder (* base (expmod base (- exp 1) m)) m))))


(define (expmod-2 prev m)
	(define sq (remainder (square prev) m))
	; Test if prev is sqrt of 1.
	; The value sq is square of prev in mod m.
	; If sq is 1, then prev is sqrt of 1
	(if
		(= sq 1)
		(cond
			((= prev 1      ) 1) ; trivial.
			((= prev (- m 1)) 1) ; trivial.
			(else 0)) ; nontrivial sqrt of 1. signal to expmod procedure.
		sq)) ; not sqrt of 1.

; gosh> (ex-1.28-1)
; 561 ...  succ 249 fail 310
; 1105 ...  succ 365 fail 738
; 1729 ...  succ 593 fail 1134
; 2465 ...  succ 741 fail 1722
; 2821 ...  succ 929 fail 1890
; 6601 ...  succ 1649 fail 4950
; -- non-primes --
; 39 ...  succ 5 fail 32
; 19999 ...  succ 29 fail 19968
; -- primes --
; 199 ...  succ 197 fail 0
; 1999 ...  succ 1997 fail 0
; 1009 ...  succ 1007 fail 0
; 1013 ...  succ 1011 fail 0
; 1019 ...  succ 1017 fail 0
; 10007 ...  succ 10005 fail 0
; 10009 ...  succ 10007 fail 0
; 10037 ...  succ 10035 fail 0
; 100003 ...  succ 100001 fail 0
; 100019 ...  succ 100017 fail 0
; 100043 ...  succ 100041 fail 0
; #<undef>