SICP 問題 1.28 （Miller-Rabinテスト）

【問題】

騙されないFermatテストの変形の一つは、Miller-Rabinテスト（Miller-Rabin test）（Miller 1976: Rabin 1980）という。
これは n が素数であり、a を n より小さい任意の整数とすると、a の (n-1) 乗は n を法として 1 と合同であるという、Fermatの小定理のもう一つの方から始める。
Miller-Rabinテストで数nの素数性をテストするには a < n をランダムにとり、expmod手続きを使って n を法とした a の (n-1) 乗を作る。
しかし、expmodで二乗計算しながら「nを法とした1の自明でない平方根」、つまり1でもn-1でもないがその二乗がnを法として1になる数がなかったか調べる。そういうnを法とした1の自明でない平方根があれば、nが素数でないことは証明できる。
また、nが素数でない奇数なら、a
つまり1でもn-1でもないがその二乗がnを法として1になる数がなかったか調べる。ある判定対象となる数をxとして、上述の文言を数式化すると、

x ≠ 1

x ≠ n-1

x^2 modulo n ＝ 1

と表現できるかな。こいつを手続き化してみよう。

(define (miller-rabin-test x n)
  (if (and (not (= x 1))
	   (not (= x (- n 1)))
	   (= (remainder (* x x) n) 1))
      #t
      #f))

あ、これif使わなくてもandの判定結果をそのまま返せばいいのか。

(define (miller-rabin-test x n)
  (and (not (= x 1))
       (not (= x (- n 1)))
       (= (remainder (* x x) n) 1)))

これをexpmod手続き中のどこかで呼び出す感じ？
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))))

このなかで、判定ケースをもう一つふやしてやればいいんでね？

;お馴染み、冪剰余を算出する関数
(define (expmod base exp m)
  (cond ((= exp 0) 1)
	((miller-rabin-test base m) 1) ;★ここで新しく定義した判定を追加
        ((even? exp)
         (remainder (square (expmod base (/ exp 2) m))
                    m))
        (else
         (remainder (* base (expmod base (- exp 1) m))
                    m))))

じゃ、例のCarmichael数をチェックしてみよう。

gosh> (map realy-prime? '(561 1105 1729 2465 2821 6601))

(#f #f #f #f #f #f)

gosh> 

おおお！ちゃんと非素数判定されてる！
一応他の数もチェックしてみようかね。
指定した区間の数の素数を抽出する手続きを書いてみる。

(define (primes from to answer)
  (if (>= from to)
      (reverse answer)
      (if (realy-prime? from)
	  (primes (+ from 1) to (cons from answer))
	  (primes (+ from 1) to answer))))

読み込ませて次の式を評価。

gosh> (primes 1 561 '())

(1 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557)

gosh> 

ちゃんと561は「素数ではない」と除外されてる。
一応、この結果を一番信頼性のあるsamllest-divisorを使ったprime?手続きでチェック！

gosh> (map prime? (primes 1 561 '()))

(#t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t #t)

gosh> 

いいんでない？