計算機プログラムの構造と解釈 第二版 P29 問題1.28

[追記20090226]
これソースコードのところ、間違ってます。

教科書の以下のくだりに目をつけた。

つまり1でもn-1でもないがその二乗がnを法として1になる数がなかったか調べる。そういうnを法とした1の自明でない平方根があれば、nが素数でないことが証明できる

ぎりぎりわかりそうだ。
上記の転載部分の前半を考える。
つまり、ある数の二乗がnを法として1になればいい。

a^2 ≡ 1 (mod n)
a は 1 < a < n-1
のとき、nは素数でないということがわかった。

つまりだ、
wikipedia 合同式
http://ja.wikipedia.org/wiki/%E5%90%88%E5%90%8C%E5%BC%8F

二つの整数 a, b は、a − b が自然数 n で割り切れるとき、n を法として合同（n をほうとしてごうどう、congruent modulo n ）であるといい、

ってかいてあるから、
(a^2 -1) % m
が0になるかどうか確かめればいい。

あとは、書いてあるとおりで、expmodでその数があったら、0を返すようにすれば何とかなる。

素数とか、Carmichael数とかでもテストしてみたけどたぶんいけてる。
もう、能率とか考えずに泥臭いコードを書いてみました。

と、書いてる間にミス発覚、変更点1にand条件がいくつか必要。

以下ソースコード。

  1 #!/usr/local/bin/gosh
  2 ;; -*- coding: utf-8 -*-
  3 
  4 (use ggc.debug.trace)
  5 (use math.mt-random)
  6 
  7 ;;random;;;;;;;;;;;;;;;;;;;;;;
  8 (define mt (make <mersenne-twister> :seed (sys-time)))
  9 (define (random n)
 10   (mt-random-integer mt n))
 11 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 12 
 13 
 14 (define (square n)
 15   (* n n))
 16 
 17 
 18 
 19 (define (expmod base exp m)
 20   (cond ((= exp 0) 1)
 21         ;;ここ変更点1
 22         ((= exp 2)
 23          (cond ((= (remainder (- (square base) 1) m) 0) 0)
 24                (else
 25                  (remainder (square (expmod base (/ exp 2) m))
 26                     m))))
 27        ((even? exp) 
 28          (remainder (square (expmod base (/ exp 2) m))
 29                     m))
 30         (else
 31           (remainder (* base (expmod base (- exp 1) m))
 32                      m))))
 33 
 34 (define (miller-rabin-test n)
 35   (define (try-it a)
 36     ;;ここ変更点2
 37     (and (= (expmod a n n) a) (not(= (expmod a n n) 0))))
 38   (try-it (+ 1 (random (- n 1)))))
 39 
 40 (define (fast-prime? n times)
 41   (cond ((= times 0) #t)
 42         ((miller-rabin-test n) (fast-prime? n (- times 1)))
 43         (else #f)))
 44 
 45 (define (prime? n)
 46   (fast-prime? n 5))
 47 
 48 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 49 
 50 
 51 ;; main
 52 (define (main args)
 53 
 54   (display "cmiller-rabin-test 素数とすでに確認された数が素数かチェック")(newline)
 55   (display "1009:")(print (prime? 1009))
 56   (display "1013:")(print (prime? 1013))
 57   (display "1019:")(print (prime? 1019))
 58   (display "10007:")(print (prime? 10007))
 59   (display "10009:")(print (prime? 10009))
 60   (display "10037:")(print (prime? 10037))
 61   (newline)
 62 
 63   (display "偶数と素数でない奇数を素数かチェック")(newline)
 64   (display "9:")(print (prime? 9))
 65   (display "20:")(print (prime? 20))
 66   (display "21:")(print (prime? 21))
 67   (display "49:")(print (prime? 49))
 68   (display "81:")(print (prime? 81))
 69   (display "100:")(print (prime? 100))
 70   (newline)
 71 
 72   (display "cmiller-rabin-test charmichaelナンバーが、素数かチェック")(newline)
 73   (display "テストを行う。素数性が確認されれば#t")(newline)
 74   (display "561:")(print (prime? 561))
 75   (display "1105:")(print (prime? 1105))
 76   (display "1729:")(print (prime? 1729))
 77   (display "2465:")(print (prime? 2465))
 78   (display "2821:")(print (prime? 2821))
 79 
 80   (newline)
 81  0)