Testing for Primality
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