SICP 1.2.6 Example: Testing for Primality - プログラミング再入門

整数計算のなかなか難しい話が沢山出て来ますが、調べまくった内容を兎に角メモしておきます。

ノート

例：素数判定

Searching for divisors（約数を探す）

基本的には2から始めて順に割り切れるかを確かめる。問題はいくつまで確かめれば十分なのか。n/2より大きい数字を確かめる必要は無い事は明らか。ここではnが素数でない場合、必ず $sqrt{n}$ 以下に約数を持つと言う性質から $sqrt{n}$ までの数を確かめている。

脚注にある通り、dがnの約数であれば、n/dは整数。d × n / d = nと表せる。dとn/dの二つの整数はちょうど反比例の関係で、どちらかの小さい方の数を考えた時に、その最小値そのものはdとn/dが同じ数の時に最大になる。dとn/dが同じ数の時、d= $sqrt{n}$ なので、n最小の約数は $sqrt{n}$ を超えない事が分かる。

プログラムの基本は：

候補約数（test-divisor）が $sqrt{n}$ より大きかったらn自身が最小の約数。即ちnは素数。
候補約数でnが割り切れたら、その候補約数がnの最小約数。
候補約数をひとつ増やして再帰
候補約数は２から始める。
割り切れるかを判定する為にdivise?と言う補助関数（述語）を使っている。
その上で最小約数が自身と同じであれば素数と判定している。

The Fermat test（フェルマーテスト）

フェルマーの小定理を利用。

フェルマーの小定理: 素数nと、nより小さい任意の正の整数aがあるとき、anをnで割った余りはaをnで割った余りと一致する。

と書いてある。
二つの数が『congruent modulo nである時』とは互いにnで割ったときの余りが等しい時を言う。n進数に直した時の１の位が同じとも言える。

よく見かける定義は

a^n-1≡1 (mod n)　nは素数。aはnとは互いに素な任意の正の整数。

つまり、a^n-1をnで割った余りは常に１。

この二つの形式が同じである事を示すのはそんなに単純な話ではない。
一般の定義から
$a^{n-1} = nb + 1$ （bは未知）
両辺にaを掛けると、
$a^{n} = anb + a$
つまりaⁿをnで割った余りはa × n × b + aをnで割った余りと等しい。

より前の部分はaの倍数なので、a × n × b + aをnで割った余りは、a / nの余りである。

ここではa < nに限定しているので、a / nの余りはaである。

纏めるとaⁿをnで割った余りはaである（nは素数、aはnより小さい任意の正の整数）。と言うのがこの本に書いてあるフェルマーの小定理。aとnが互いに素、即ちaは１より大きい事になるが、aが１でも取り敢えずこの定理は成り立っている。

nが素数でない場合、nより小さい任意の正の整数aでは上記の定理は殆ど成り立たない。仮に成り立ったとしても別の任意の数aを選んで同じ判定を行えばnが素数である確率は更に上がる。ただし素数ではないのに任意のaに対してこの定理が成り立ってしまう数が知られており、それらはカーマイケル数と呼ばれている。

このaⁿ%m計算する所でfast-exptで使ったテクニックで計算量を減らすのだが、その根拠が以下の通り。二つの数xとyをそれぞれmで割る時、
$x=md_{x}+(x%m)$
$y=md_{y}+(y%m)$
と表すとすると、 $x\times y$ を考えると
$x\times y = m^{2}d_{x}d_{y}+md_{x}(y%m)+md_{y}(x%m)+(x%m)(y%m)$
右辺は最後の項を除いてmを含んでいるので（つまりmの倍数なので）、両辺をmで割った余りは
$(x\times y)%m = {(x%m)(y%m)}%m$
つまりxとyの積をmで割った余りは、xをmで割った余りとyをmで割った余りの積を更にmで割った余りに等しい。
これを応用して、expが偶数の場合にxとyを共にbase^exp/2として当てはめて
$(base^{\exp})%m$
$=(base^{\exp/2}\times base^{\exp/2})%m$
$=(base^{\exp/2}%m)^{2}%m$ （expが偶数の場合）
またexpが奇数の場合はx=base^exp-1、y=baseとして当てはめて
$(base^{\exp})%m$
$=(base^{\exp-1}\times base)%m$
$=(base^{\exp-1}%m\times base%m)%m$ （expが奇数の場合）
これをコードのしたのがexpmod。ただしbaseはmで割り切れない事が前提なのか奇数の場合のbase%mは求めずbaseをそのまま使用している。

fermat-testのプログラムを走らせるにあたり、R5RSには乱数を返す関数が定義されていないとの事なので、DrRacketのPLaneTにあるモジュールを利用する。言語設定を「ソースで決める」にしておいて、先頭に下記の２行を追加する。

#lang racket
(require (planet neil/sicp:1:16))

最初に少しダウンロード等が走る。fast-prime?を実際に動かしてみると

> (fast-prime? 2 3)
#t
> (fast-prime? 5 2)
#t
> (fast-prime? 41 2)
#t
> (fast-prime? 97 2)
#t
> (fast-prime? 99 2)
#f
> (fast-prime? 101 2)
#t
> (fast-prime? 103 2)
#t
> (fast-prime? 107 2)
#t
> (fast-prime? 109 2)
#t
> (fast-prime? 111 2)
#f
> (fast-prime? 123 2)
#f
>

またカーマイケル数は素数ではないが、何度繰り返した所で素数として判定される。

> (fast-prime? 561 10)
#t
> (fast-prime? 1105 100)
#t
> (fast-prime? 1729 100)
#t
>

Probablistic Method（確率的方法（確率的アルゴリズム））

フェルマーテストは偽となればその数が素数でない事が確定するが、真であっても素数である事が確定する訳ではない。何度テストをやっても真となる数はテストの回数を増やす程素数である確率が上がって欲しい所であるが、残念ながらそう言う訳でもない。素数ではないのに常にフェルマーの小定理を満たす数が存在する（カーマイケル数。7/10000個程度）。ただその数は十分に少ないので、フェルマーテストは十分実用的である。

フェルマーテストにはバリエーションがある。フェルマーテストの様に任意に選んだa (a < n)とnの関係を調べnが素数かを判定するもの、対照的にnが素数でない限り殆ど満たされない条件を証明する物もある。
従ってひとつのテストを通ると正しい確率は1/2より大きくなり、二つのテストに通ると正しい確率は3/4よりも大きくなる。このテストを繰り返せば正しい確率はどんどん上がる。（と書いてあるが本当か？）

この様に繰り返すほど間違う確率が下がる方法は確率的アルゴリズムと呼ばれる。

Exercise 1.21

find-divisorに出現するsquareだけ自分で定義して：
実行結果

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

Exercise 1.22

与えられた引数nが素数だったとき***と計算に要した時間を表示するプログラム：timed-prime-test。これを使ってある範囲の奇数が素数かどうか判定するプログラムを書く。
問題文にあるruntimeはR5RSには存在しないらしい。そこで言語をRaketにする必要があるが、そうすると今度はRacketがelseが無いifは許していないので若干問題文の関数を書き換える必要がある。

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

特に意味はないがelseに入ったらfalseを返しておく。

runtimeに相当する関数としてはcurrent-millisecondsとcurrent-inexact-millisecondsがあるが、整数値のcurrent-millisecondsでは計算が速すぎて計算量が0としか表示されないので、浮動小数点のcurrent-inexact-millisecondsを使う。ただし、この時間は実際の時刻を取って来るので、他のプロセスやOSが消費する時間も含まれてしまう。current-process-millisecondsと言う関数もあるが、これはやはり単位が大きすぎて結果が全て0となってしまう。

(define (runtime) (current-milliseconds))

計算する関数は与えられた範囲の奇数列を返すodd-sequenceと、その各要素にtimed-prime-testを適用するmapを呼び出す様にsearch-for-primesを定義する。最後にmapが結果のリストを返しそれが表示されてしまうので、それを防ぐ為に余分なnewlineを入れてある。

(define (odd-sequence min max)
  (cond ((even? min) (odd-sequence (+ 1 min) max))
        ((> min max) '())
        (else (cons min (odd-sequence (+ 2 min) max)))))

(define (search-for-primes min max)
  (map timed-prime-test (odd-sequence min max))
  (newline))

あと、R5RSにはsquareも無いので定義しておく。
実行結果：

> (search-for-primes 1000 1030)

1001
1003
1005
1007
1009 *** 0.004150390625
1011
1013 *** 0.004150390625
1015
1017
1019 *** 0.004150390625
1021 *** 0.004150390625
1023
1025
1027
1029
>

これを使って概ねの計算量を測定する。他のプロセス等の影響を見る為に何度かデータを取ってみる。

素数	１回目	２回目	３回目	４回目
1009	0.00390625	0.004150390625	0.004150390625	0.00390625
1013	0.00390625	0.004150390625	0.003173828125	0.00390625
1019	0.00390625	0.00390625	0.00390625	0.00390625

どうも３種類しか結果が出て来ない。周囲の影響は遅くなる方にしか作用しない筈と割り切って考えて最速の値を採用する事にする。ここでは0.003173828125。
同様に

素数	１回目	２回目	３回目	４回目
10007	0.010986328125	0.010009765625	0.010986328125	0.010986328125
10009	0.010009765625	0.010986328125	0.01123046875	0.010009765625
10037	0.010986328125	0.010986328125	0.010986328125	0.011962890625

最速は0.010009765625。

素数	１回目	２回目	３回目	４回目
100003	0.032958984375	0.033935546875	0.032958984375	0.0341796875
100019	0.031982421875	0.031982421875	0.031982421875	0.034912109375
100043	0.031005859375	0.03515625	0.035888671875	0.031005859375

最速は0.031005859375。

素数	１回目	２回目	３回目	４回目
1000003	0.10009765625	0.10986328125	0.1220703125	0.11181640625
1000033	0.097900390625	0.096923828125	0.156005859375	0.09716796875
1000037	0.157958984375	0.113037109375	0.156005859375	0.096923828125

最速は0.096923828125。

1000近傍の素数を計算するのに要する時間に対して10000近傍の場合は√10倍の時間が掛かると予想されているが、

> (* 0.003173828125 (sqrt 10))
0.010036525776901594
> (* 0.003173828125 (sqrt 100))
0.03173828125
> (* 0.003173828125 (sqrt 1000))
0.10036525776901595
>

Do your timing data bear this out?

の問いには、まあまあ1000近くの計算量に対して10000の場合は√10倍、100000の場合は√100=10倍、1000000の場合は√1000と言えそう。

How well do the data for 100,000 and 1,000,000 support the n prediction?

にはどう答えるべきか。誤差の％であれば、

> (* (/ (- 0.010009765625 (* 0.003173828125 (sqrt 10))) (* 0.003173828125 (sqrt 10))) 100)
-0.26662764084341956
> (* (/ (- 0.031005859375 (* 0.003173828125 (sqrt 100))) (* 0.003173828125 (sqrt 100))) 100)
-2.307692307692308
> (* (/ (- 0.096923828125 (* 0.003173828125 (sqrt 1000))) (* 0.003173828125 (sqrt 1000))) 100)
-3.4289053010118025
>

0.3%、-2.3%、-3.4%と言った所。

Is your result compatible with the notion that programs on your machine run in time proportional to the number of steps required for the computation?

これも何と答えるのか。。。1000に対してその10倍、100倍、1000倍の数の計算をして、√10、√100、√1000に対して数パーセント以内の誤差であればO(n)よりは遥かにO(logn)に合致していると言って良いのかな。

Exercise 1.23

最小の約数かどうかを2から順に判定するのに、2が約数ではなかった場合にそれ以降偶数をテストするのは無駄なので奇数のみを返す関数にする。本当はエラストネスの篩の様な感じでそれまでにテストした数の倍数は除いて返す（結果として素数のみを返す事になる筈）関数の方がよりsmallest-divisiorの計算量は減る。ただ素数を返す関数の方がややリソースを食う気がする。
ここでは初期値は2から始めるが、与えられた引数より大きい奇数を返す関数nextを定義する。一般化すると引数が偶数の時は+1、奇数の時は+2する、となるが問題文の通り引数が2の時は3を、それ以外は引数+2を返すと言うここだけの特殊な関数でも良い。

(define (next pre)
  (if (= pre 2)
      3
      (+ pre 2)))

書き換えたfind-divisorはこう：

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

timed-prime-testを実行してみる。

> (timed-prime-test 1009)

1009 *** 0.005859375

どうも、Exercise 1.22で行った奇数のみを連続でテストする時と条件が異なる模様。search-for-primesでテストしてみる。

> (search-for-primes 1000 1030)

1001
1003
1005
1007
1009 *** 0.003173828125
1011
1013 *** 0.00390625
1015
1017
1019 *** 0.003173828125
1021 *** 0.0029296875
1023
1025
1027
1029
>

比較をするためにはこちらの数値の方が良さそう。
測定の制度の問題かこちらの方が『若干』速いだけの様に見える。

オーダー	Exercise 1.22	今回	比率
1000代	0.003173828125	0.0029296875	92%
10000代	0.010009765625	0.0068359375	68%
100000代	0.031005859375	0.02099609375	67%

やはり1000代の数では何か他の影響によるオフセットが大きいのか、10000代以上で漸く速くなったと言えそう。でも時間にして2/3なので1.5倍程度。

how do you explain the fact that it is different from 2?

この説明はなかなか難しい。
smallest-divisorとdvidedes?を呼び出す回数は減る（半減）代わりに

nextの呼び出し
ifによる条件分岐

が増えている。
全体としてはそれ以外の部分もあるので元々オフセットされている。
ただ、1.5倍程度が妥当かどうかは簡単には判断出来ない。
もっと大きな数で試してみると、100000000代で、0.10888671875が0.06396484375になるので59%。100000000代で0.966064453125が0.626220703125になるので65%程度と、オフセットは余り大きくはないのかもしれない。

Exercise 1.24

prime?の代わりにfast-prime?を使用。繰り返しは取り敢えず２としておく。

(define (start-prime-test n start-time)
  (if (fast-prime? n 2)
      (report-prime (- (runtime) start-time))
      #f))

実行してみる。

> (timed-prime-test 1009)

1009 *** 0.01904296875
>

数回実行して最速の値を取る。1000代

素数	計算時間
1009	0.01904296875
1013	0.018798828125
1019	0.019775390625

10000代

素数	計算時間
10007	0.02099609375
10009	0.02099609375
10037	0.02294921875

100000代

素数	計算時間
100003	0.05419921875
100019	0.052978515625
100043	0.048095703125

1000000代

素数	計算時間
1000003	0.06298828125
1000033	0.060791015625
1000037	0.066162109375

how would you expect the time to test primes near 1,000,000 to compare with the time needed to test primes near 1000?

計算量の議論でのlogの底は多分2だが、底が何であれlog(1000000)=log(1000^^2^^)=2log(1000)なので、n=1000000の時の計算量はn=1000の時の計算量の２倍辺りが期待される。

Do your data bear this out?

log10000/log1000=4/3
10000の時の実測値は1000の時の凡そ２倍。4/3倍が期待されていた事を考えると、その更に3/2で1.5倍。

log100000/log1000=5/3
100000の時の実測値は1000の時の凡そ2.5倍。5/3倍が期待されていた事を考えると、やはりその更に3/2で1.5倍。

log1000000/log1000=6/3=2
1000000の時の実測値は1000の時の凡そ３倍。２倍が期待されていた事を考えると、やはりその更に1.5倍。

Can you explain any discrepancy you find?

説明は出来ないが、Θ(log n)と言うのは単純に計算量がlog nに比例していると言うより、計算量の増加がlog nの増加に比例していると言う事か。一定の係数なのでn log nとも違うし。

Exercise 1.25

Alyssaの実装を使ってみる。

> (timed-prime-test 1000037)

1000037 *** 15346.0830078125
>

元々0.066162109375だった事と比較すると圧倒的に遅くなっている。

Alyssaが書いたプログラムでのステップ数はほぼ同じに見える。違うとするとremainderを呼び出す回数とその第１引数の大きさ。元のexpmodの方が呼び出す回数が多い代わりに第１引数の数は高々base²に押さえられる（本当は高々baseに押さえられる筈だが）。

例えば元の実装

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

で5³%4の計算過程を考えてみる。
$5^3%4$
$=(5^2%4\times 5)%4$
$=\left[\left{(5^1%4)^2\right}%4 \times 5\right]%4$
$=\left[\left{(5^0%4\times 5)^2\right}%4\times 5\right]%4$
$=\left[\left{(1\times 5)^2\right}%4\times 5\right]%4$
$=(25%4\times 5)%4$
$=(1\times 5)%4$
$=5%4$
$=1$
一方、Alyssaの実装はremainderを呼び出す回数が１回な代わりに第１引数はbase^expなのでexpが大きいと非常に大きな数となる。フェルマーテストでは大きな素数の判定に使うのでexpが非常に大きく、base^expは天文学的数字である。つまり非常に大きい数字に対するremainderの性能が響いていると思われる。

DrRacketのプロファイラではremainderの実行時間は測ってくれないので、数値は出せないが割り算の効率は一般的に悪く、特にマシンの整数レジスタに入りきらない大きな整数の割り算は更に効率が悪いと思われるので、remainderを呼び出す回数が多くても小さい数の割り算にとどめる方が実行時間が短いと予想出来る。

Exercise 1.26

expが偶数のときexpmodを２回呼んでしまうので、本来ここで半分になる筈の計算量が全く減っていない。

expが偶数のとき次のexpを半分にする代わりにexpを２回呼んぶ。
expが奇数の時は従来通りexpを1回呼んで、次のexpを１減らす。

であればexpmodは結局exp回呼ぶ事になるの。ここのexpが素数判定しようとしている数なのでn回expmodを呼んでいる事になる。つまりΘ(n)オーダーの計算量となっている。

Exercise 1.27

フェルマーテストを乱数に対して行うのではなく単にフェルマーの小定理を確認する述語にする。

(define (fermat-test a n)
  (= (expmod a n n) a))

動作確認

> (fermat-test 3 561)
#t
>

全てのaをテストする為に単純な数列を形成する簡易版のrangeを作る。

(define (range begin end)
  (define (make-list a)
    (if (= a end) '() (cons a (make-list (+ a 1)))))
  (make-list begin))

動作確認：

> (range 2 13)
(2 3 4 5 6 7 8 9 10 11 12)
>

全てのaに対してfermat-testするfermat-test-all-rangeを定義する。

(define (fermat-test-all-range n r)
    (map (lambda (a) (fermat-test a n)) r))

動作確認

> (fermat-test-all-range 13 (range 2 13))
(#t #t #t #t #t #t #t #t #t #t #t)
>

これらを使って素数判定をする述語関数prime?を定義する。ここでは確率的どころかa

(define (prime? n)
  (fold (lambda (res a) (and res a)) #t (fermat-test-all-range n (range 2 n))))

foldを使うのにSRFI-1を必要とする。ここではR6RS用のライブラリが必要。DrRacketでは以下の１文で取り込む。

(#%require (lib "srfi/%3a1"))

動作確認

> (prime? 10)
#f
> (prime? 11)
#t
> (prime? 12)
#f
> (prime? 13)
#t
>

カーマイケル数を素数判定してみる。

> (prime? 561)
#t
> (prime? 1105)
#t
> (prime? 1729)
#t
> (prime? 2465)
#t
> (prime? 2821)
#t
> (prime? 6601)
#t
>

全て素数と判定される。
smallest-divisorで約数を捜すと、

> (smallest-divisor 561)
3
> (smallest-divisor 1105)
5
> (smallest-divisor 1729)
7
> (smallest-divisor 2465)
5
> (smallest-divisor 2821)
7
> (smallest-divisor 6601)
7
>

となり、全て素数ではない事が分かる。

Exercise 1.28

ミラー-ラビン素数判定法。

Theorem, which states that 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.

これはフェルマーの小定理（一般に見る形）。

whenever we perform the squaring step in expmod, we check to see if we have discovered a "nontrivial square root of 1 modulo n," that is, a number not equal to 1 or n - 1 whose square is equal to 1 modulo n.

expmodで２乗の計算をする時に、従来は２乗を計算してnで割った余りを関数値としていた所を、

まず、２乗する値が１かn - 1でない事を確認する
その数を２乗してnで割った余りが１であれば、自明でない平方根が見つかった事になる。ここでは0を返す事にする。
そうでない場合はその余りを関数値とする

(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (letrec ((root (expmod base (/ exp 2) m))
                  (rem (remainder (square root) m)))
           (if (and (and (not (= root 1)) (not (= root (- m 1))))
                    (= rem 1))
               0
               rem)))
        (else
         (remainder (* base (expmod base (- exp 1) m))
                    m))))

(define (prime? n)
  (and (millar-rabin-test n) (millar-rabin-test n)))

素数判定を行うのにmillar-rabin-testを２回呼び出す事にする。

実行結果。Exercise 1.24で既出の素数に対しては#tと判定し、カーマイケル数に対しては#fを返している。ただし6601で１回騙されているが、２回目で正解してた。

> (prime? 1009)
#t
> (prime? 1013)
#t
> (prime? 1019)
#t
> (prime? 1021)
#t
> (prime? 10007)
#t
> (prime? 10009)
#t
> (prime? 10037)
#t
> (prime? 561)
#f
> (prime? 1105)
#f
> (prime? 1729)
#f
> (prime? 2465)
#f
> (prime? 2821)
#f
> (prime? 6601)
#t
> (prime? 6601)
#f
>