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

  • 0
    いいね
  • 0
    コメント

    #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 さんのミラーラビンはこちら