sicp-1-2-6-primality-testing

#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.