Permalink
Browse files

added 1.21 - 1.28

  • Loading branch information...
1 parent c6acfe9 commit a5ec58a634b342838f082ea9e98d588cc7e8a71e @sarabander committed Jul 3, 2011
Showing with 312 additions and 0 deletions.
  1. +24 −0 1.2/Ex1.21.scm
  2. +54 −0 1.2/Ex1.22.scm
  3. +23 −0 1.2/Ex1.23.scm
  4. +59 −0 1.2/Ex1.24.scm
  5. +52 −0 1.2/Ex1.25.scm
  6. +4 −0 1.2/Ex1.26.scm
  7. +16 −0 1.2/Ex1.27.scm
  8. +80 −0 1.2/Ex1.28.scm
View
24 1.2/Ex1.21.scm
@@ -0,0 +1,24 @@
+(define (square x) (* x x))
+
+(define (smallest-divisor n)
+ (find-divisor n 2))
+
+(define (find-divisor n test-divisor)
+ (cond ((> (square test-divisor) n) n)
+ ((divides? test-divisor n) test-divisor)
+ (else (find-divisor n (+ test-divisor 1)))))
+
+(define (divides? a b)
+ (= (remainder b a) 0))
+
+(define (prime? n)
+ (= n (smallest-divisor n)))
+
+(smallest-divisor 199) ; 199
+(smallest-divisor 1999) ; 1999
+(smallest-divisor 19999) ; 7
+
+; more tests
+(smallest-divisor 561) ; 3
+(smallest-divisor 241) ; 241
+(smallest-divisor 169) ; 13
View
54 1.2/Ex1.22.scm
@@ -0,0 +1,54 @@
+;; Needs the procedures of 1.21
+
+;; Racket needs this
+(define (runtime) (current-milliseconds))
+
+(define (timed-prime-test n)
+ (start-prime-test n (runtime)))
+
+;; The procedures differ from book originals
+(define (start-prime-test n start-time)
+ (if (prime? n)
+ (report-prime n (- (runtime) start-time))
+ false))
+
+(define (report-prime n elapsed-time)
+ (newline)
+ (display n)
+ (display " *** ")
+ (display elapsed-time)
+ (newline))
+
+(timed-prime-test 242)
+(time (timed-prime-test 4398042316799))
+(time (timed-prime-test 1125899839733759))
+
+(define (search-for-primes starting-with how-many)
+ (cond ((zero? how-many) (void))
+ ((odd? starting-with)
+ (if (timed-prime-test starting-with)
+ (search-for-primes (+ 2 starting-with) (- how-many 1))
+ (search-for-primes (+ 2 starting-with) how-many)))
+ (else (search-for-primes (+ 1 starting-with) how-many))))
+
+(search-for-primes 1000 3) ; 1009, 1013, 1019
+(search-for-primes 10000 3) ; 10007, 10009, 10037
+(search-for-primes 100000 3) ; 100003, 100019, 100043
+(search-for-primes 1000000 3) ; 1000003, 1000033, 1000037
+
+;; search finishes faster than the runtime granularity (1 ms)
+;; need to use larger numbers:
+
+(search-for-primes 1000000000 3) ; 1000000007, 1000000009, 1000000021
+;; average time 16 ms
+(search-for-primes 10000000000 3) ; 10000000019, 10000000033, 10000000061
+;; average time 42 ms
+(search-for-primes 100000000000 3) ; 100000000003, 100000000019, 100000000057
+;; average time 85 ms
+(search-for-primes 1000000000000 3); 1000000000039, 1000000000061, 1000000000063
+;; average time 200 ms
+
+(* 30 (sqrt 10))
+
+;; Little less than (sqrt 10) growth in time when input grows 10 times.
+;; Seems to support the proportionality of number of steps and time.
View
23 1.2/Ex1.23.scm
@@ -0,0 +1,23 @@
+;; Depends on 1.21 and 1.22
+
+;; modified from 1.21 to use next
+(define (find-divisor n test-divisor)
+ (cond ((> (square test-divisor) n) n)
+ ((divides? test-divisor n) test-divisor)
+ (else (find-divisor n (next test-divisor)))))
+
+(define (next input)
+ (if (= input 2)
+ 3
+ (+ 2 input)))
+
+(search-for-primes 1000000000 3) ; 1000000007, 1000000009, 1000000021
+;; average time 9 ms
+(search-for-primes 10000000000 3) ; 10000000019, 10000000033, 10000000061
+;; average time 20 ms
+(search-for-primes 100000000000 3) ; 100000000003, 100000000019, 100000000057
+;; average time 40 ms
+(search-for-primes 1000000000000 3); 1000000000039, 1000000000061, 1000000000063
+;; average time 115 ms
+
+;; Yes, runs about twice as fast as in 1.22
View
59 1.2/Ex1.24.scm
@@ -0,0 +1,59 @@
+;; Depends on 1.22
+
+;; How many times to run fermat-test
+(define tries 100000)
+
+;; Modified to use fast-prime?
+(define (start-prime-test n start-time)
+ (if (fast-prime? n tries)
+ (report-prime n (- (runtime) start-time))
+ false))
+
+(define (fast-prime? n times)
+ (cond ((= times 0) true)
+ ((fermat-test n) (fast-prime? n (- times 1)))
+ (else false)))
+
+(define (fermat-test n)
+ (define (try-it a)
+ (= (expmod a n n) a))
+ (try-it (+ 1 (random (- n 1)))))
+
+;; To count how many times expmod divides exponent by two
+(define halvings 0)
+;; and decrements it
+(define decrements 0)
+
+(define (expmod base exp m)
+ (cond ((= exp 0) 1)
+ ((even? exp)
+ (and (set! halvings (+ halvings 1))
+ (remainder (square (expmod base (/ exp 2) m))
+ m)))
+ (else
+ (and (set! decrements (+ decrements 1))
+ (remainder (* base (expmod base (- exp 1) m))
+ m)))))
+
+(search-for-primes 1000 3)
+;; average time 540 ms
+(search-for-primes 1000000 3)
+;; average time 970 ms
+
+;; Ratio of execution times should be:
+(/ (log 1000000) (log 1000)) ; 2
+;; experiment almost confirms this: ratio of (/ 970 540 1.0) is 1.8,
+;; the small difference probably comes from this:
+;;
+;; We could guess that when exponentiating with successive squaring,
+;; the process is less efficient with small exponents than with large ones.
+;; With small exponents the algorithm probably decremets the exponent
+;; relatively more often than it divides by 2. Let's test the hypothesis:
+
+;; when using 1000 and 1000000
+(/ halvings decrements 1.0) ; 1.13 and 2.04, respectively
+
+;; indeed, our guess is confirmed!
+
+(set! halvings 0)
+(set! decrements 0)
View
52 1.2/Ex1.25.scm
@@ -0,0 +1,52 @@
+(define (expmod base exp m)
+ (remainder (fast-expt base exp) m))
+
+;; recursive process
+(define (fast-expt b n)
+ (cond ((= n 0) 1)
+ ((even? n) (square (fast-expt b (/ n 2))))
+ (else (* b (fast-expt b (- n 1))))))
+
+;; iterative process from Ex 1.16
+(define (fast-expt b n)
+ (fast-expt-iter 1 b n))
+
+(define (fast-expt-iter a b n)
+ (cond ((= n 0) a)
+ ((even? n) (fast-expt-iter a (square b) (/ n 2)))
+ ((odd? n) (fast-expt-iter (* a b) b (- n 1)))))
+
+;; In principle, we could use fast-expt, but it makes the calculation
+;; terribly slow. I would guess it's because of the procedure-calling
+;; overhead (environment creation) which is magnified by 100000 trials.
+;; Better to use inline expressions here.
+
+(search-for-primes 1000 3)
+;; average time 7600 ms
+
+;; This hypothesis is supported by the observation that using the
+;; iterative fast-expt makes search-for-primes even slower, because
+;; of the additional call to fast-expt-iter, which makes bigger
+;; environment with 3 parameters.
+
+(search-for-primes 1000 3)
+;; average time 8900 ms
+
+;; If we call fast-expt-iter directly, we would expect faster execution:
+
+(define (expmod base exp m)
+ (remainder (fast-expt-iter 1 base exp) m))
+
+(search-for-primes 1000 3)
+;; average time 8900 ms
+
+;; No, it isn't any faster. This would suggest the overhead comes mainly
+;; from environment creation.
+
+;; ...
+
+;; After some pondering, I would disregard much of the above explanation.
+;; The main reason of slowness is that fast-expt produces a huge number,
+;; and operations with those are slow. The original expmod keeps the
+;; numbers small by repeatedly taking remainders before squaring or
+;; multiplying. Alyssa's expmod calculates the remainder only once.
View
4 1.2/Ex1.26.scm
@@ -0,0 +1,4 @@
+;; The multiplication expression evaluates expmod twice, resulting in binary
+;; process tree, which grows exponentially as O(2^n). This cancels the
+;; advantage we had with O(log n) algorithm, because 2^(log n) = n.
+
View
16 1.2/Ex1.27.scm
@@ -0,0 +1,16 @@
+;; Depends on 1.21 and 1.24
+
+(define (fermat-prime? n)
+ (define (fast-prime n a)
+ (cond ((= a 0) true)
+ ((fermat-test n a) (fast-prime n (- a 1)))
+ (else false)))
+ (fast-prime n (- n 1)))
+
+(define (fermat-test n a)
+ (define (try-it a)
+ (= (expmod a n n) a))
+ (try-it a))
+
+(map fermat-prime? '(561 1105 1729 2465 2821 6601)) ; Carmichael no-s fool it,
+(map prime? '(561 1105 1729 2465 2821 6601)) ; but not this
View
80 1.2/Ex1.28.scm
@@ -0,0 +1,80 @@
+;; Depends on 1.24
+
+;; How many times to run prime-test
+(define tries 1000)
+
+;; Probabilistic test
+(define (miller-rabin-prime? n)
+ (define (fast-prime n times)
+ (cond ((= times 0) true)
+ ((prime-test n) (fast-prime n (- times 1)))
+ (else false)))
+ (fast-prime n tries))
+
+(define (prime-test n)
+ (define (try-it a)
+ (= (expmod a (- n 1) n) 1))
+ (try-it (+ 1 (random (- n 1)))))
+
+;; ;; Exhaustive test
+;; (define (miller-rabin-prime? n)
+;; (define (fast-prime n a)
+;; (cond ((= a 0) true)
+;; ((prime-test n a) (fast-prime n (- a 1)))
+;; (else false)))
+;; (fast-prime n (- n 1)))
+
+;; (define (prime-test n a)
+;; (define (try-it a)
+;; (= (expmod a (- n 1) n) 1))
+;; (try-it a))
+
+;; Modified to signal the discovery of nontrivial square root of 1
+;; by returning 0
+(define (expmod base exp m)
+ (cond ((= exp 0) 1)
+ ((even? exp)
+ (let* ((aexp (expmod base (/ exp 2) m))
+ (rem (remainder (square aexp)
+ m)))
+ (if (= rem 1)
+ (if (or (= aexp 1) (= aexp (- m 1)))
+ 1
+ 0) ; nontrivial sqrt of 1 discovered!
+ rem)))
+ (else (remainder (* base (expmod base (- exp 1) m))
+ m))))
+
+;; Tests
+;; -----
+
+;; Carmichael numbers
+(map miller-rabin-prime? '(561 1105 1729 2465 2821 6601 8911 10585 15841))
+; all false
+
+;; Primes from Wikipedia:
+
+;; Centered decagonal primes
+(map miller-rabin-prime?
+ '(11 31 61 101 151 211 281 661 911 1051 1201 1361 1531)) ; all true
+
+;; Dihedral primes
+(map miller-rabin-prime?
+ '(2 5 11 101 181 1181 1811 18181 108881 110881 118081 120121 121021))
+; all true
+
+;; Fibonacci primes
+(map miller-rabin-prime?
+ '(2 3 5 13 89 233 1597 28657 514229)) ; all true
+
+;; Mersenne primes
+(map miller-rabin-prime?
+ '(3 7 31 127 8191 131071 524287)) ; all true
+
+;; Primes # 301-320
+(map miller-rabin-prime?
+ '(1993 1997 1999 2003 2011 2017 2027 2029 2039 2053
+ 2063 2069 2081 2083 2087 2089 2099 2111 2129)) ; all true
+
+;; Non-primes
+(map miller-rabin-prime? '(6 9 46 78 27 93)) ; all false

0 comments on commit a5ec58a

Please sign in to comment.