#lang racket
;; 1.21
(define (divides? a b)
(zero? (remainder b a)))
(define (find-divisor n test-divisor)
(cond [(> (sqr test-divisor) n) n]
[(divides? test-divisor n) test-divisor]
[else (find-divisor n (+ test-divisor 1))]))
(define (smallest-divisor n)
(find-divisor n 2))
(smallest-divisor 199)
;; 199
(smallest-divisor 1999)
;; 1999
(smallest-divisor 19999)
;; 7
;; 1.22
(define (my-prime? n)
(= n (smallest-divisor n)))
;; Time my-prime? to verify it is O(sqrt(n)).
(define (timed-prime-test n)
(define (start-prime-test num start-time)
(cond [(my-prime? num)
(printf "~a is prime *** ~a ms ~n"
num
(inexact->exact
(truncate
(- (current-inexact-milliseconds) start-time))))
#t]
[else #f]))
(start-prime-test n (current-inexact-milliseconds)))
(define (find-primes start times)
; find and display times many primes, with timing info, beginning at start
(cond [(zero? times) (newline)]
[(timed-prime-test start) (find-primes (add1 start) (sub1 times))]
[else (find-primes (add1 start) times)]))
;; I had to make the numbers big to get noticeable elapsed times. The starting
;; numbers go up by a factor of 100. So with a sqrt(n) time complexity, the
;; elapsed times should be going up by a factor of 10 each time, which we do see.
(find-primes 10000000000000 3)
;; 10000000000037 is prime *** 115 ms
;; 10000000000051 is prime *** 106 ms
;; 10000000000099 is prime *** 115 ms
(find-primes 1000000000000000 3)
;; 1000000000000037 is prime *** 1201 ms
;; 1000000000000091 is prime *** 1200 ms
;; 1000000000000159 is prime *** 1203 ms
(find-primes 100000000000000000 3)
;; 100000000000000003 is prime *** 12009 ms
;; 100000000000000013 is prime *** 12025 ms
;; 100000000000000019 is prime *** 12140 ms
;; 1.23
;; Redo 1.22 with an improved find-divisor that does not unnecessarily check even
;; numbers greater than 2.
(define (next n) ; return 3 if n is 2 else n + 2
(if (= n 2) 3 (+ n 2)))
(define (new-find-divisor n test-divisor)
; same as find-divisor but uses next instead of incrementing by 1
(cond [(> (sqr test-divisor) n) n]
[(divides? test-divisor n) test-divisor]
[else (new-find-divisor n (next test-divisor))]))
(define (new-smallest-divisor n)
(new-find-divisor n 2))
(define (new-my-prime? n)
(= n (new-smallest-divisor n)))
(define (new-timed-prime-test n)
(define (start-prime-test num start-time)
(cond [(new-my-prime? num)
(printf "~a is prime *** ~a ms ~n"
num
(inexact->exact
(truncate
(- (current-inexact-milliseconds) start-time))))
#t]
[else #f]))
(start-prime-test n (current-inexact-milliseconds)))
(define (new-find-primes start times)
(cond [(zero? times) (newline)]
[(new-timed-prime-test start) (new-find-primes (add1 start) (sub1 times))]
[else (new-find-primes (add1 start) times)]))
(new-find-primes 10000000000000 3)
;; 10000000000037 is prime *** 82 ms
;; 10000000000051 is prime *** 62 ms
;; 10000000000099 is prime *** 79 ms
(new-find-primes 1000000000000000 3)
;; 1000000000000037 is prime *** 748 ms
;; 1000000000000091 is prime *** 748 ms
;; 1000000000000159 is prime *** 770 ms
(new-find-primes 100000000000000000 3)
;; 100000000000000003 is prime *** 7588 ms
;; 100000000000000013 is prime *** 7591 ms
;; 100000000000000019 is prime *** 7603 ms
;; The new elapsed times are a little more than half the old times. I think this
;; reflects the fact that timed-prime-test has some unavoidable printing overhead
;; that the improved my-prime? cannot help. Nonetheless I think this is a
;; significant improvement, not a big-O improvement, but a nice leading
;; coefficient improvement.
;; 1.24
;; Redo 1.22 using fast-prime? instead of my-prime? and check that we get
;; logarithmic time complexity.
(define (expmod base expo m)
(cond [(zero? expo) 1]
[(even? expo)
(remainder (sqr (expmod base (/ expo 2) m)) m)]
[else (remainder (* base (expmod base (sub1 expo) m)) m)]))
(define (fermat-test n)
(define (try-it a)
(= (expmod a n n) a))
(try-it (add1 (random (min 4294967087 ; biggest number accepted by random
(sub1 n))))))
(define (fast-prime? n times)
(cond [(zero? times) #t]
[(fermat-test n) (fast-prime? n (sub1 times))]
[else #f]))
(define (fast-timed-prime-test n)
; same as timed-prime-test but uses fast-prime? instead of my-prime?
(define (start-prime-test num start-time)
(cond [(fast-prime? num 10) ; run the fermat test ten times
(printf "~a is prime *** ~a ms ~n"
num
(inexact->exact
(truncate
(- (current-inexact-milliseconds) start-time))))
#t]
[else #f]))
(start-prime-test n (current-inexact-milliseconds)))
(define (fast-find-primes start times)
(cond [(zero? times) (newline)]
[(fast-timed-prime-test start) (fast-find-primes (add1 start) (sub1 times))]
[else (fast-find-primes (add1 start) times)]))
;; First let's compare with the numbers we used earlier.
(fast-find-primes 10000000000000 3)
;; 10000000000037 is prime *** 0 ms
;; 10000000000051 is prime *** 0 ms
;; 10000000000099 is prime *** 0 ms
(fast-find-primes 1000000000000000 3)
;; 1000000000000037 is prime *** 0 ms
;; 1000000000000091 is prime *** 0 ms
;; 1000000000000159 is prime *** 0 ms
(fast-find-primes 100000000000000000 3)
;; 100000000000000003 is prime *** 0 ms
;; 100000000000000013 is prime *** 0 ms
;; 100000000000000019 is prime *** 0 ms
;; So we got the same primes but this time basically instantaneously.
;; Now we will have to use huge numbers to register any elapsed time at all. With
;; log running time, doubling the number of digits in a number should double the
;; elapsed time.
(fast-find-primes (inexact->exact 1e25) 3)
;; 10000000000000000905969697 is prime *** 2 ms
;; 10000000000000000905969749 is prime *** 1 ms
;; 10000000000000000905969823 is prime *** 1 ms
(fast-find-primes (inexact->exact 1e50) 3)
;; 100000000000000007629769841091887003294964970946821 is prime *** 4 ms
;; 100000000000000007629769841091887003294964970947027 is prime *** 4 ms
;; 100000000000000007629769841091887003294964970947049 is prime *** 4 ms
(fast-find-primes (inexact->exact 1e100) 3)
;; 100000000000000001590289110975991804683608085639452813897813275577478387
;; 72170381060813469985856815251 is prime *** 15 ms
;; 100000000000000001590289110975991804683608085639452813897813275577478387
;; 72170381060813469985856815363 is prime *** 14 ms
;; 100000000000000001590289110975991804683608085639452813897813275577478387
;; 72170381060813469985856815719 is prime *** 17 ms
(fast-find-primes (inexact->exact 1e200) 3)
;; 999999999999999969733122212510361659474503275455023626482417509503468484
;; 355540755341963384047062518680275124159738824081821357343682784846393850
;; 41047239877871023591066789981811181813306167128854888513 is prime *** 69 ms
;; 999999999999999969733122212510361659474503275455023626482417509503468484
;; 355540755341963384047062518680275124159738824081821357343682784846393850
;; 41047239877871023591066789981811181813306167128854888901 is prime *** 73 ms
;; 999999999999999969733122212510361659474503275455023626482417509503468484
;; 355540755341963384047062518680275124159738824081821357343682784846393850
;; 41047239877871023591066789981811181813306167128854890043 is prime *** 71 ms
;; The run times did not double each time, they go up by a factor
;; of 4, but any constant factor increase is still evidence of logarithmic
;; time. So our test is successful.
;; Note how (inexact->exact 1e200) is not 1 followed by 200 zeroes but is rather a
;; 200-digit number beginning with many 9's. I guess that's due to some machine
;; limitation.
;; 1.25
;; Alyssa's expmod code does not go modulo after each squaring operation, but only
;; takes one remainder at the end. Note that when we use expmod in the fermat
;; test, the number n we are testing for primality is the exponent, and this
;; number can be very large. So the intermediate numbers in Alyssa's expmod
;; process will be super huge and probably cause some kind of register overflow
;; error.
;; 1.26
;; expmod is a recursive procedure and it's time complexity is the number of nodes
;; in the recursion tree. By doubling the number of recursive calls, Louis's
;; expmod turns the recursion tree from a single branch with log(n) nodes, into a
;; full binary tree with n nodes.
;; 1.27
(define carmichael '(561 1105 1729 2465 2821 6601))
(define (passes-fermat? n)
(for/and ([a (in-range 2 n)])
(= (expmod a n n) a)))
(passes-fermat? 5)
;; #t
(passes-fermat? 6)
;; #f
(passes-fermat? 7)
;; #t
(passes-fermat? 8)
;; #f
(passes-fermat? 9)
;; #f
(for/and ([n carmichael]) (passes-fermat? n))
;; #t
(for/or ([n carmichael]) (my-prime? n))
;; #f
;; So the carmichael numbers are non-primes that fool the fermat test.
;; 1.28
;; the Miller-Rabin primality test
(define (special-expmod base expo m)
; same as expmod but immediately returns 0 on encountering
; a non-trivial square root of unity
(cond [(zero? expo) 1]
[(even? expo)
(define u (special-expmod base (/ expo 2) m))
(define u-squared (remainder (sqr u) m))
(if (and (= u-squared 1)
(> u 1)
(< u (sub1 m)))
0
u-squared)]
[else (remainder (* base (special-expmod base (sub1 expo) m)) m)]))
(define (miller-rabin-prime? n [times 10]) ; run 10 times by default
(define (try-it a)
(= (special-expmod a (sub1 n) n) 1))
(if (< n 10)
(member n '(2 3 5 7))
(for/and ([i (in-range times)])
; do not test 0, 1, or n - 1
(try-it (+ 2 (random (min 4294967087 (- n 3))))))))
;; easy tests
(miller-rabin-prime? 10)
;; #f
(miller-rabin-prime? 11)
;; #t
(miller-rabin-prime? 12)
;; #f
(miller-rabin-prime? 13)
;; #t
(miller-rabin-prime? 14)
;; #f
(miller-rabin-prime? 15)
;; #f
(miller-rabin-prime? 16)
;; #f
(miller-rabin-prime? 17)
;; #t
;; big primes from earlier problems
(miller-rabin-prime? 10000000000037)
;; #t
(miller-rabin-prime? 10000000000051)
;; #t
(miller-rabin-prime? 10000000000099)
;; #t
(miller-rabin-prime? 1000000000000037)
;; #t
(miller-rabin-prime? 1000000000000091)
;; #t
(miller-rabin-prime? 1000000000000159)
;; #t
(miller-rabin-prime? 100000000000000003)
;; #t
(miller-rabin-prime? 100000000000000013)
;; #t
(miller-rabin-prime? 100000000000000019)
;; #t
;; close but no cigar
(miller-rabin-prime? (* 10000000000037 10000000000051))
;; #f
(miller-rabin-prime? (* 1000000000000091 1000000000000159))
;; #f
(miller-rabin-prime? (* 100000000000000003 100000000000000019))
;; #f
;; the toughies
(for/or ([n carmichael])
(miller-rabin-prime? n))
;; #f
;; lastly let's time it on some pretty big primes
(time (miller-rabin-prime? 100000000000000007629769841091887003294964970946821))
;; cpu time: 16 real time: 9 gc time: 0
;; #t
(time (miller-rabin-prime? 100000000000000007629769841091887003294964970947027))
;; cpu time: 15 real time: 8 gc time: 0
;; #t
(time (miller-rabin-prime? 100000000000000007629769841091887003294964970947049))
;; cpu time: 0 real time: 5 gc time: 0
;; #t
;; So Miller-Rabin looks like a very efficient and trustworthy primality checker.
Comments