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

とりあえずは、ここに書いてありますexpmodを実装してみる。

  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 ;;runtime;;;;;;;;;;;;;;;;;;;;;;
 15 (define (runtime)
 16   (use srfi-11)
 17   (let-values (((a b) (sys-gettimeofday)))
 18               (+ 1000000 b)))
 19 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 20 
 21 ;;1-21;;;;;;;;;;;;;;;;;;;;;;;;;
 22 (define (square n)
 23   (* n n))
 24 
 25 
 26 (define (smallest-divisor n)
 27   (find-divisor n 2))
 28 
 29 (define (find-divisor n test-divisor)
 30   (cond ((> (square test-divisor) n) n)
 31         ((divides?  test-divisor n) test-divisor)
 32         (else (find-divisor n
 33                             (+ test-divisor 1)))))
 34 
 35 
 36 (define (divides? a b) 
 37   (= (remainder b a) 0))
 38   
 39 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 40 
 41 
 42 ;;P28 - P29;;;;;;;;;;;;;;;;;;;;;;;;;
 43 
 44 ;;ここは前までの実装なのでcomment out
 45 ;(define (expmod base exp m) ;base = a, exp =n, m = n
 46 ;  (cond ((= exp 0) 1)
 47 ;        ((even? exp)
 48 ;         (remainder (square (expmod base (/ exp 2) m))
 49 ;                    m))
 50 ;        (else
 51 ;          (remainder (* base (expmod base (- exp 1) m))
 52 ;                     m))))
 53 
 54 ;;問題1.25 expmod実装
 55 (define (fast-expt b n) ;base = a, exp =n, m = n
 56   (cond ((= n 0) 1)
 57         ((even? n) (square (fast-expt b (/ n 2))))
 58         (else (* b (fast-expt b (- n 1))))))
 59 
 60 (define (expmod base exp m)
 61   (remainder (fast-expt base exp) m))
 62 ;;;問題1.25 expmod実装ここまで。
 63 
 64 
 65 
 66 (define (fermat-test n)
 67   (define (try-it a)
 68     (= (expmod a n n) a))
 69   (try-it (+ 1 (random (- n 1)))))
 70 
 71 (define (fast-prime? n times)
 72   (cond ((= times 0) #t)
 73         ((fermat-test n) (fast-prime? n (- times 1)))
 74         (else #f)))
 75 
 76 (define (prime? n) 
 77   (fast-prime? n 5))
 78 
 79 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 80 
 81 
 82 (define (timed-prime-test n)
 83   (newline)
 84   (display n)
 85   (start-prime-test n (runtime)))
 86 
 87 (define (start-prime-test n start-time)
 88   (if (prime? n)
 89     (report-prime (- (runtime) start-time))))
 90 
 91 (define (report-prime elapsed-time)
 92   (display " *** ")
 93   (display elapsed-time))
 94 
 95 (define (search-for-primes-iter num counter)
 96 (cond ((= counter 3))
 97       ((prime? num) (print num) (search-for-primes-iter (+ num 2) (+ counter 1)))
 98       (else (search-for-primes-iter (+ num 2) counter))))
 99 
100 
101 (define (search-for-primes num)
102   (cond ((divides? 2 num) (search-for-primes-iter (+ num 1) 0))
103         (else (search-for-primes-iter num 0))))
104 
105 
106 ;; main
107 (define (main args)
108   (search-for-primes 1000)
109   (search-for-primes 10000)
110   (search-for-primes 100000)
111   (search-for-primes 1000000)
112   (search-for-primes 10000000)
113   (search-for-primes 100000000)
114   (search-for-primes 1000000000)
115   (search-for-primes 10000000000)
116   (search-for-primes 100000000000)
117 
118 
119   (timed-prime-test 1009)
120   (timed-prime-test 1013)
121   (timed-prime-test 1019)
122 
123   (timed-prime-test 10007)
124   (timed-prime-test 10009)
125   (timed-prime-test 10037)
126 
127   (timed-prime-test 100003)
128   (timed-prime-test 100019)
129   (timed-prime-test 100043)
130 
131   (timed-prime-test 1000003)
132   (timed-prime-test 1000033)
133   (timed-prime-test 1000037)
134 
135   (timed-prime-test 10000019)
136   (timed-prime-test 10000079)
137   (timed-prime-test 10000103)
138 
139   (timed-prime-test 100000007)
140   (timed-prime-test 100000037)
141   (timed-prime-test 100000039)
142 
143   (timed-prime-test 1000000007)
144   (timed-prime-test 1000000009)
145   (timed-prime-test 1000000021)
146 
147   (timed-prime-test 10000000019)
148   (timed-prime-test 10000000033)
149   (timed-prime-test 10000000061)
150 
151   (timed-prime-test 100000000003)
152   (timed-prime-test 100000000019)
153   (timed-prime-test 100000000057)
154   (newline)
155  0)

実行。。。　。　　。　　　。　　　　　。　　　　　　。　　　　　　　。

遅すぎる！！
なので、とりあえず、使えないことはわかる。

なぜ遅いか。考えてみた。

問題1.25のexpmod実装だと、べき乗を求めてからremainderを実行しているしている。

その前の実装だとremainderのタイミングが多い。
remainderを逐次実行すると、積算で扱う数が小さくてすむ。
たとえばここでいうなら、通常のexmodであれば、(m-1)^2が最高にでかい乗算になる。
桁数が大きいと平方は、重そうなイメージ
だからだと思う。たぶん。きっと。