brian_hassinger

software developer

Testing for Primality

February 11, 2013

Section 1.2.5 is short, and right on topic with 1.2.6 so for this post I’ve joined them together.

Exercise 1.20

Illustrate the process shape for (gcd 206 40) for both normal-order evaluation and applicative-order.

To refresh some memory, in normal-order evaluation, all operators and operands are expanded until they become a primitive expression. In applicative-order all operators are expanded once and then reduced. In both forms of evaluation the special-form if is evaluated the same. Where the condition is evaluated first, then the appropriate consequent or alternative.

;; For brevity, use % instead of remainder
(define (% a b) (remainder a b))

;; Normal-order evaluation
(gcd 206 40)
(if (= 40 0) a (gcd b (% a b)))
(gcd 40 (% 206 40))
(if (= (% 206 40) 0) a (gcd b (% a b)))
(gcd (% 206 40) (% 40 (% 206 40)))
(if (= (% 40 (% 206 40)) 0) a (gcd b (% a b)))
(gcd (% 40 (% 206 40)) (% (% 206 40) (% 40 (% 206 40))))
(if (= (% (% 206 40) (% 40 (% 206 40))) 0) a (gcd b (% a b)))
(gcd (% (% 206 40) (% 40 (% 206 40)))
     (% (% 40 (% 206 40)) (% (% 206 40) (% 40 (% 206 40)))))
(if (= (% (% 40 (% 206 40)) (% (% 206 40) (% 40 (% 206 40)))) 0) a (gcd b (% a b)))
(% (% 206 40) (% 40 (% 206 40)))
2

;; Applicative-order evaluation
(gcd 206 40)
(if (= 40 0) a (gcd b (% a b)))
(gcd 40 (% 206 40))
(gcd 40 6)
(if (= 6 0) a (gcd b (% a b)))
(gcd 6 (% 40 6))
(gcd 6 4)
(if (= 4 0) a (gcd b (% a b)))
(gcd 4 (% 6 4))
(gcd 4 2)
(if (= 2 0) a (gcd b (% a b)))
(gcd 2 (% 4 2))
(gcd 2 0)
(if (= 0 0) a (gcd b (% a b)))
2

In the normal-order evaluation of the (gcd 206 40) process the remainder operation is performed 14 times before the process meets the conditions to finish evaluation. Then the process still needs to reduce the answer and performs 4 more remainder operations for a total of 18 calls to the remainder procedure. The applicative-order evaluation calls the remainder procedure 4 times.

Exercise 1.21

Find the smallest divisor of the following numbers: 199, 1999, 19999.

Simple enough, setup the procedure (given in the book) and run it.

(smallest-divisor 199)
;Value: 199
(smallest-divisor 1999)
;Value: 1999
(smallest-divisor 19999)
;Value: 7

The results show that 199 and 1999 are both prime numbers. However 19999 is a factor of 7.

Exercise 1.22

Show how long it takes to find primes within a given domain.

Starting to get a little more interesting here. Almost all the steps are setup, the only thing left is to iterate over a problem domain.

(define (search-for-primes min max)
  (cond ((> min max) (newline) (display "done"))
        ((even? min) (search-for-primes (+ min 1) max))
        (else (timed-prime-test min) (search-for-primes (+ min 2) max))))

Ok great! Stop when the min is larger than the max, if the min is even make it odd, and finally test for prime and try the next iteration. I didn’t like the printing non-prime numbers so I changed the output a little bit.

(search-for-primes 1000 1019)
1009 *** 0.
1013 *** 0.
1019 *** 0.
done

(search-for-primes 10000 10037)
10007 *** 0.
10009 *** 0.
10037 *** 0.
done

(search-for-primes 100000 100043)
100003 *** 0.
100019 *** 0.
100043 *** 0.
done

(search-for-primes 1000000 1000037)
1000003 *** 9.999999999999995e-3
1000033 *** 1.0000000000000009e-2
1000037 *** 9.999999999999995e-3
done

The most interesting bit is the gap between CPUs in the 80’s and today (2013). When I run the search for the primes in the domains asked for, my computer can hardly tell it’s thinking. I’m not even on a powerful computer, these tests were run on a netbook with an Intel Atom 1.8Ghz processor. The scheme community suggests upping the anti with a problem domain from 1e9 to 1e12.

(search-for-primes 1000000000 1000000021)
1000000007 *** .25
1000000009 *** .23999999999999994
1000000021 *** .25
done

(search-for-primes 10000000000 10000000061)
10000000019 *** .8299999999999998
10000000033 *** .8300000000000001
10000000061 *** .81
done

(search-for-primes 100000000000 100000000057)
100000000003 *** 2.67
100000000019 *** 2.6900000000000004
100000000057 *** 2.66
done

(search-for-primes 1000000000000 1000000000063)
1000000000039 *** 8.499999999999998
1000000000061 *** 8.5
1000000000063 *** 8.48
done

Ok. 8 seconds. tsk tsk Atom processor.

In accordance with the book’s expectations, each exponential step is roughly (sqrt 10) times the computation time as the previous exponent. This is obvious because the times are different by a factor of ~3 which is ~= (sqrt 10).

Exercise 1.23

Improve the speed of th timed-tests by creating a next procedure to reduce the number of calculations on unnecessary numbers.

(define (next input)
  (if (= input 2) 3 (+ input 2)))

(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)))))

Running the same scaled tests from the last exercise:

(search-for-primes 1000000000 1000000021)
1000000007 *** .16
1000000009 *** .15999999999999998
1000000021 *** .15000000000000002
done

(search-for-primes 10000000000 10000000061)
10000000019 *** .52
10000000033 *** .5
10000000061 *** .5
done

(search-for-primes 100000000000 100000000057)
100000000003 *** 1.6
100000000019 *** 1.6
100000000057 *** 1.62
done

(search-for-primes 1000000000000 1000000000063)
1000000000039 *** 5.09
1000000000061 *** 5.099999999999998
1000000000063 *** 5.07
done

The results are improved, but they are not half the speed of the first tests. The ratio for each is roughly 2/3 speed increase.

Exercise 1.24

Modify the timed-test to use the Fermat method.

I would expect prime numbers of any size using the fast-prime? procedure to take logarithmic time to complete. When running the tests the differences are minimal and increase slightly, suggesting logarithmic growth.

(search-for-primes 1000000000 1000000021)
1000000007 *** 0.
1000000009 *** 9.999999999999995e-3
1000000021 *** 0.
done

(search-for-primes 10000000000 10000000061)
10000000019 *** 9.999999999999995e-3
10000000033 *** 9.999999999999995e-3
10000000061 *** 9.999999999999995e-3
done

(search-for-primes 100000000000 100000000057)
100000000003 *** 0.
100000000019 *** 9.999999999999995e-3
100000000057 *** 1.0000000000000009e-2
done

(search-for-primes 1000000000000 1000000000063)
1000000000039 *** 1.0000000000000009e-2
1000000000061 *** 1.0000000000000009e-2
1000000000063 *** 9.999999999999981e-3
done

Exercise 1.25

This procedure does work to calculate the expmod

(expmod 10 10 10)
;; 0
(expmod 2 10 10)
;; 4

However when run with the fast-prime? procedure the process takes much longer to compute than the other solution. The key to understanding why this takes longer is hidden in the evaluation. Let’s use (expmod 2 5 5)

Remember the definition of fast-expt is:

(define (fast-expt b n)
  (cond ((= n 0) 1)
        ((even? n) (square (fast-expt b (/ n 2))))
        (else (* b (fast-expt b (- n 1))))))

For Alyssa P. Hacker’s expmod the expansion:

(expmod 2 5 5)
(remainder (fast-expt 2 3) 5)
(remainder (* 2 (fast-expt 2 4)) 5)
(remainder (* 2 (square (fast-expt 2 2))) 5)
(remainder (* 2 (square (square (fast-expt 2 1)))) 5)
(remainder (* 2 (square (square (* 2 (fast-expt 2 0))))) 5)
(remainder (* 2 (square (square (* 2 1)))) 5)
(remainder (* 2 (square (square 2))) 5)
(remainder (* 2 (square 4)) 5)
(remainder (* 2 16) 5)
(remainder 32 5)
2

And the faster expmod definition used in the book:

(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (remainder (square (expmod base (/ exp 2) m))
                    m))
        (else
         (remainder (* base (expmod base (- exp 1) m))
                    m))))

With the expansion:

(expmod 2 5 5)
(remainder (* 2 (expmod 2 4 5)) 5)
(remainder (* 2 (remainder (square (expmod 2 2 5)) 5)) 5)
(remainder (* 2 (remainder (square (remainder (square (expmod 2 1 5)) 5)) 5)) 5)
(remainder (* 2 (remainder (square (remainder (square (remainder (* 2 (expmod 2 0 5)) 5)) 5)) 5)) 5)
(remainder (* 2 (remainder (square (remainder (square (remainder (* 2 1) 5)) 5)) 5)) 5)
(remainder (* 2 (remainder (square (remainder (square (remainder 2 5)) 5)) 5)) 5)
(remainder (* 2 (remainder (square (remainder (square 2) 5)) 5)) 5)
(remainder (* 2 (remainder (square (remainder 4 5)) 5)) 5)
(remainder (* 2 (remainder (square 4) 5)) 5)
(remainder (* 2 (remainder 16 5)) 5)
(remainder (* 2 1) 5)
(remainder 2 5)
2

The original expmod requires a larger space compared to Alyssa P. Hacker’s expmod because the remainder is calculated at every step. This reduction is important in larger numbers (e.g 1e9) because the square procedure has a smaller number to calculate at every step. In Alyssa’s procedure the number to square will grow quickly making each step take longer for the scheme compiler to complete.

Exercise 1.26

The reason the square procedure is able to compute faster in the expmod procedure is shown in the process tree below.

;; Process Tree
;; Using (%) for (remainder)
(expmod 2 3 4)
(% (* 2 (expmod 2 2 4)) 4)
(% (* 2 (% (* (expmod 2 1 4) (expmod 2 1 4)) 4)) 4) ;; Excessive Operations
(% (* 2 (% (* (% (* 2 (expmod 2 0 4)) 4) (% (* 2 (expmod 2 0 4)) 4)) 4)) 4)
(% (* 2 (% (* (% (* 2 1) 4) (% (* 2 1) 4)) 4)) 4)
(% (* 2 (% (* (% 2 4) (% 2 4)) 4)) 4)
(% (* 2 (% (* 2 2) 4)) 4)
(% (* 2 (% 4 4)) 4)
(% (* 2 0) 4)
(% 0 4)
0

On third step of the process (expmod 2 3 4) the primitive procedure * is given two identical processes. Both of these will produce the same result, but the compiler does not know this and must compute each one separately. The square procedure is able to reuse the result of evaluating its operator once, thus not having to do twice as much work.

Exercise 1.27

Testing to see if the Carmichael numbers can fool the fermat test.

(define (extended-fermat-test n)
  (define (iter a)
    (cond ((> a n) #t)
          ((= (expmod a n n)
              (remainder a n)) (iter (+ a 1)))
          (else #f)))
  (iter 2))

;; Test - Carmichael Numbers
(extended-fermat-test 561)
(extended-fermat-test 1105)
(extended-fermat-test 1729)
(extended-fermat-test 2465)
(extended-fermat-test 2821)
(extended-fermat-test 6601)
;; These numbers pass the fermat-test, but are not prime.

The extended-fermat-test procedure will perform the fermat-test on every number (starting with 2) lower than the given n. The procedure will return #false if given a non-prime number and return #true if given a prime. When given a Carmichael number, the fermat-test is fooled into thinking it has found a prime number, when it has not.

;; Test - Carmichael Numbers
(extended-fermat-test 561)
;Value: #t
(extended-fermat-test 1105)
;Value: #t
(extended-fermat-test 1729)
;Value: #t
(extended-fermat-test 2465)
;Value: #t
(extended-fermat-test 2821)
;Value: #t
(extended-fermat-test 6601)
;Value: #t

Exercise 1.28

Using the miller-rabin test to check for prime numbers.

(define (expmod base exp m)
  (define (trivial? num root)
    (if (and (not (= 1 num))
             (not (= (- m 1) num))
             (= 1 root))
        0
        root))
  (cond ((= exp 0) 1)
        ((even? exp)
         (let ((num (expmod base (/ exp 2) m)))
           (trivial? num (remainder (square num) m))))
        (else
          (remainder (* base (expmod base (- exp 1) m))
                    m))))

(define (miller-rabin n)
  (define (try-it a)
    (= (expmod a (- n 1) n) 1))
  (try-it (+ 1 (random (- n 1)))))

The miller-rabin procedure will return #t if the test determines n to be prime or #f if the test finds a non-trivial root of n. This test chooses a random number a < n and tries to find a non-trivial root. Because of this random selection, some non-prime numbers are reported as prime. To correctly determine if a number is prime this test should be performed multiple times.

;; Test Runner

(define (fast-prime? n times)
  (cond ((= times 0) true)
        ((miller-rabin n) (fast-prime? n (- times 1)))
        (else false)))

(define (test n)
  (newline)
  (display n)
  (display ": ")
  (if (fast-prime? n 100)
    (display "Prime")
    (display "Non-prime")))

;; Tests

;; Mix primes and non-primes
(test 5)
(test 6)
(test 7)
(test 9)

;; Carmichael Numbers
(test 561)
(test 1105)
(test 1729)
(test 2465)
(test 2821)
(test 6601)

With the results:

;; Mix primes and non-primes
5: Prime
6: Non-prime
7: Prime
9: Non-prime

;; Carmichael Numbers
561: Non-prime
1105: Non-prime
1729: Non-prime
2465: Non-prime
2821: Non-prime
6601: Non-prime