Scheme
SICP

オンラインSICP読書女子会 #6(1.2.6 素数判定-後編) まとめ

More than 1 year has passed since last update.

#ladiescpp の主催してる オンラインSICP読書女子会 #6 (1.2.6) - connpassのまとめです〜。

@hio-harp さんのまとめ: オンラインSICP読書女子会 #6

前回までの復習

フェルマーの小定理

nが素数で、nがaよりも小さい任意の正の整数であるとき,
以下が成り立つ

a^n \equiv a \pmod{n} 
(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))))


expmodがやってること

偶数の時

\frac{base^{exp}}{m} \equiv (\frac{base^\frac{exp}{2}}{m})^2 \pmod{m} 

奇数の時

\frac{base^{exp}}{m} \pmod{m} \equiv (\frac{base^{exp-1} * base}{m})^2 \pmod{m} (奇数の時)

1.25 expmodの書き換え

expmodを以下のように書いてはだめか

(define (expmod base exp m) (remainder (fast-expt 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))))

fast-expt: 1.2.4

(define (fast-expt b n)
  (cond ((= n 0) 1)
    ((even? n) (square (fast-expt b (/ n 2))))
    (else (* b (fast-expt b (- n 1))))))

提案された方法では、一度以下を計算してから、その余りを求めている

base^{exp}

この方法だと、
一度上で計算された数値が巨大になるため、
例えばbase = 10^9, exp = 10^9 - 1の時、
従来のexpmodであれば計算可能だが、
提案方法では一度桁が9 ^ (10^9)の数字を生成する必要があり、
整数値としてもてない。(intやlongの値に収まらない)

1.26 expmodの計算量

元のexpmodの計算量は O(log n)
squareの部分によってO(n) -> O(log n)になっていっている。

毎回試行回数を半分にしていったのに、
提案方法だと、毎回2倍呼ぶようにしてしまったから、(相殺して)
結局試行回数はO(n)になってしまう。

1.27 カーマイケル数

カーマイケル数が実際にフェルマーテストを騙すことを実験せよ。

カーマイケル数

フェルマーテストを通過してしまう合成数
561, 1105, 1729, 2465, 2821....

;; base ^ exp % m を求める関数
(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)))) 

;; all
(define (fermat-test n m)
  (= (expmod m n n) m))


;; fermat-all-testを使う
(define (prime? n times)
  (cond ((= times 0)
     (display "Passed test. Prime number ")
     (display n))
    ((fermat-test n times) (prime? n (- times 1)))
    (else
     (display "Failed:")
     (display n))))

(define (fermat-prime? n)
  (prime? n (- n 1))
  (newline))

;; test
(fermat-prime? 4) 
(fermat-prime? 24) 
(fermat-prime? 11) 
(fermat-prime? 19) 

;; カーマイケル数
(fermat-prime? 561)

;;Failed:4
;;Failed:24
;;Passed test. Prime number 11
;;Passed test. Prime number 19
;;Passed test. Prime number 561 <- 通過してしまう

1.28 ミラーラビンテスト

騙すことができないフェルマーテストの変種のひとつにミラー-ラビンテスト (Miller-Rabin test) がある。
これは、次のようなフェルマーの小定理の別の形から始める。
n が素数で、かつ a が n 以下の任意の正の整数 であれば、
a の (n − 1) 乗は法 n に関して 1 と合同であるというものだ。

Miller-Rabin テストで数値 n の素数性をテストするには、a < n である適当な a を選び、expmod 手続きを使って a の (n − 1) 乗の n を法とする剰余を求める。
しかし、expmod 手続きの中で二乗を行うステップで、毎回 “法 n に関する自明でない 1 の平方根” を見つけたかチェックする。これは、1 とも n − 1 とも等しくない数値で、その二乗が法 n に関して 1 に等しい数値である。
そのような自明でない 1 の平方根が存在すれば n は素数ではないということは証明できる。
また、もし n が素数でない奇数であれば、少なくとも a < n であるaのうち半分は、
このように an−1 を演算すると、法 n に関する自明でない 1 の平方根が現れるということも証明できる 
(このため、Miller-Rabin テストは騙すことができない)。 
expmod 手続きを変更し、自明でない 1 の平方根を見つけるとシグ ナルを送るようにして、それを使って fermat-test と同じような 手続きとして Miller-Rabin テストを実装せよ。
既知の素数、非素数を使ってその手続きをチェックせよ。

日本語が難しい・・・。 法 n に関する自明でない 1 の平方根は、つまり

m^2 \equiv 1 \pmod{n}

となるようなmのこと、という理解でよいだろうか。(*ただしmはn-1でも1でもない)

nが素数で、かつaがn以下の任意の正の整数なら以下が成り立つ

フェルマーの変形

a^{n-1} \equiv 1 \pmod{n}

こうしてくとフェルマーテストと同じことがわかる

\begin{align}
a^{n-1} \equiv 1 \pmod{n} \\
a^{n-1}*a \equiv a \pmod{n}\\
a^{n} \equiv a \pmod{n} 
\end{align}

(1.27の関数を上書きするようにしたので、ミラーラビンテストなのに名前がfermet...)

;; 1.28.scm
(load "./1.27.scm")
(print "== ミラーラビン==")
(define (expmod base exp m)
  (cond ((= exp 0) 1)
    ((even? exp) (sq-check (expmod base (/ exp 2) m) m))
    (else (remainder (* base (expmod base (- exp 1) m)) m))))

(define (sq-check n m)
  (define sq_mod (remainder (* n n) m))
  (if (= sq_mod 1)
      (cond ((= n 1) 1)
        ((= n (- m 1)) 1)
        (else 0))
      sq_mod))


;; test
(fermat-prime? 4) ; #f
(fermat-prime? 24) ; #f
(fermat-prime? 11) ; #t
(fermat-prime? 19) ; #t

;; カーマイケル数
(fermat-prime? 561) ; #t <- 通過してしまう

;; Failed:4
;; Failed:24
;; Passed test. Prime number 11
;; Passed test. Prime number 19
;; Passed test. Prime number 561
;; == ミラーラビン==
;; Failed:4
;; Failed:24
;; Passed test. Prime number 11
;; Passed test. Prime number 19
;; Failed:561

もっと詳しい @hio-harp さんのミラーラビンはこちら