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


 [#ladiescpp](http://ladiespp.connpass.com/) の主催してる [オンラインSICP読書女子会  #6 (1.2.6) - connpass](http://ladiespp.connpass.com/event/28089/)のまとめです〜。


@hio-harp さんのまとめ： [オンラインSICP読書女子会 #6](http://qiita.com/hio-harp/items/9955bf218c6c7de9193e)

## 前回までの復習

フェルマーの小定理

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

```math
a^n \equiv a \pmod{n} 
```

```lisp
(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がやってること

偶数の時

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

奇数の時

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

## 1.25 expmodの書き換え

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

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


元のexpmodは以下

```lisp
(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

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

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

```math
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....


```lisp
;; 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 の平方根```は、つまり

```math
m^2 \equiv 1 \pmod{n}
```
となるようなmのこと、という理解でよいだろうか。(*ただしmはn-1でも1でもない)



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


フェルマーの変形

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

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

```math
\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...)

```lisp
;; 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 さんのミラーラビンはこちら](http://qiita.com/hio-harp/items/9955bf218c6c7de9193e#ex-128-miller-rabin-test)