Help us understand the problem. What is going on with this article?

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

More than 3 years have 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 さんのミラーラビンはこちら

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away