Testing for Primes: Scheme
I cannot possibly look at parentheses the same way now. And that’s all thanks to Scheme. I can only see function application. Even when solving the matrix chain multiplication problem, I still had this odd idea that in fact is a function applied to some other function. Somehow that made sense.
But today, let’s talk about primes. How to test if a number is prime and how to generate the stream of primes. And which method works best. And most importantly, why does the method I thought should work best run out of memory?
I’ll try testing for primes in 3 ways: finding relevant divisors (prime?), the Fermat primality test (fermat-prime?) and the Rabin-Miller primality test (rabinmiller-prime? and rabinmiller-true-prime?). After that, we’ll try to generate the stream of prime numbers and see which of these algorithms perform better and also compare them with the tried and tested Erathosthenes’ Sieve method.
Let’s first have a look at everyone’s favorite primality test, search for divisors up to its square root:
(define divides? (lambda (a b) (= (remainder a b) 0))) (define prime? (lambda (n) (cond ((or (= n 1) (= n 0)) #f) ((= n 2) #t) ((even? n) #f) (else (let prime-test ( (d 3) ) (cond ((> (square d) n) #t) ((divides? n d) #f) (else (prime-test (+ d 2)))))))))
This is a really basic test that for a given odd number checks all odd numbers starting from 3 to see whether or not they divide . Obviously, a better approach would have been to do the test only for prime , but there’s no way of figuring that out without first having all prime numbers less than already. I remember implementing this in high-school.
Let’s get to something more interesting, the Fermat test, which is nothing but a restatement of Fermat’s Little Theorem: If is a prime number and , then .
So, everyone, open section 1.2 of the SICP and let’s have a look at exponentiation (fast-exponentiation to be exact). We’ll have this done in logarithmic time and we’ll also have the remainder calculated at each step to avoid the size of the numbers from getting out of hand:
(define expmod (lambda (base exp mod) (remainder (cond ((= exp 0) 1) ((even? exp) (square (expmod base (/ exp 2) mod))) (else (* base (expmod base (- exp 1) mod)))) mod)))
There are some issues with the Fermat primality test. Most importantly, it is not conclusive. We’ll set a fixed value of and test the condition in Fermat’s Little Theorem. If is prime, then the test is true for all values of , but testing all values of a would not be feasible.
(= (expmod a n n) a)
only proves that is probably prime, but there is no way of knowing for sure in only one test. Picking a few random values of and having the test succeed is a pretty good indicator of the primality of the number and is in fact what the PGP program uses to test for primes.
So, our Fermat test would look something like:
(define fermat-prime? (lambda (n) (let fermat-tests ((i 5)) (if (= i 0) #t (let ( (a (+ 1 (random (- n 1)))) ) (if (= (expmod a n n) a) (fermat-tests (- i 1)) #f))))))
I could have specified the number of tests to run as a parameter to fermat-prime? instead of 5. Finding a single failure is proof that the number is not prime. However, there are certain numbers that always verify this condition despite not being prime!
They are called Carmichael numbers and they are not detected by the fermat-prime? test:
1 ]=> (prime? 9746347772161)
1 ]=> (fermat-prime? 9746347772161)
In fact, .
This is where the Miller-Rabin Test comes in. It’s based on a slightly modified version of Fermat’s Little Theorem with a small twist that gets rid of Carmichael numbers as well. Instead of fixing an and checking to see if , it checks to see if .
The key difference is that when performing the squaring step in the expmod function we check to see if the number we’re supposed to be squaring is a nontrivial square root of 1 modulo n.
Let’s explain what that actually means (courtesy of the Wiki article on the Miller-Rabin test). Let be a prime number. Clearly , for any . Those are called the trivial square roots of 1. We’ll prove that if is prime, there can be no others: suppose there is another such that , so divides the product, therefore it must divide one of the factors, meaning either or which means that either or (in the field ).
If we at some point in the exponentiation get a partial result
(s (expmod2 base (/ exp 2) mod)) meaning that and without being either or , then by what we’ve proven above, cannot be prime since we’ve just found a non-trivial square root of 1.
So, we can formulate a function that works with Carmichael numbers too:
(define sqrt-of-1? (lambda (s n) (if (or (= s 1) (= s (- n 1))) #f (= (remainder (* s s) n) 1)))) (define expmod2 (lambda (base exp mod) (remainder (cond ((= exp 0) 1) ((even? exp) (let ((s (expmod2 base (/ exp 2) mod))) (if (sqrt-of-1? s mod) 0 (* s s)))) (else (* base (expmod2 base (- exp 1) mod)))) mod))) (define rabinmiller-prime? (lambda (n) (if (or (= n 0) (= n 1)) #f (let rabinmiller-tests ((i 5)) (if (= i 0) #t (let ( (a (+ 1 (random (- n 1)))) ) (if (= (expmod2 a (- n 1) n) 1) (rabinmiller-tests (- i 1)) #f)))))))
Notice however that we still only test a number 5 times. There is a proven number of times the exponentiation needs to be done in order for the test to be deterministic. Chrid Caldwell’s Primality Proving website lists the following result:
Miller’s Test [Miller76]: If the extended Riemann hypothesis is true, then if n is an a-SPRP for all integers a with 1 < a < 2(log n)2, then n is prime.
So, in fact we could adjust the code to generate the range between 2 and :
(define rabinmiller-true-prime? (lambda (n) (cond ((or (= n 0) (= n 1)) #f) ((= n 2) #t) ((even? n) #f) (else (let rabinmiller-tests ( (test-range (range 2 (min (- n 1) (* 2 (square (floor (log n))))))) ) (if (null? test-range) #t (let ((a (car test-range))) (if (= (expmod2 a (- n 1) n) 1) (rabinmiller-tests (cdr test-range)) #f))))))))
Let’s try some primitive benchmarks. I’m going to try and get the first 5000 prime numbers by filtering the stream of all naturals and benchmarking (subtracting endtime and starttime) for each method above including a sieve approach.
For that, let’s introduce a couple of new constructs for manipulating streams:
(define naturals (let make-naturals ((i 0)) (cons-stream i (make-naturals (+ i 1))))) (define filter-stream (lambda (f s) (let ( (rest (delay (filter-stream f (force (cdr s))))) (head (car s)) ) (if (f head) (cons head rest) (force rest))))) (define take (lambda (n s) (if (= n 0) '() (let ( (rest (take (- n 1) (stream-cdr s))) (head (stream-car s)) ) (cons head rest)))))
The code for the sieve is very basic (done in a lab at school):
(define sieve (lambda (s) (cons (car s) (delay (sieve (filter-stream (lambda (x) (not (= (remainder x (car s)) 0))) s)))))) (define naturals-from-2 (force (cdr (force (cdr naturals))))) (define primes (sieve naturals-from-2))
The results are very interesting:
1 ]=> (benchmark (delay (take 5000 (filter-stream rabinmiller-prime? naturals))))
1 ]=> (benchmark (delay (take 5000 (filter-stream rabinmiller-true-prime? naturals))))
1 ]=> (benchmark (delay (take 5000 (filter-stream prime? naturals))))
1 ]=> (benchmark (delay (take 5000 primes)))
;Aborting!: out of memory
;GC #32: took: 0.37 (51%) CPU time, 0.47 (57%) real time; free: 14882
;GC #33: took: 0.34 (100%) CPU time, 0.36 (99%) real time; free: 14938
Something is seriously wrong with the sieve based approach…
The interesting part is that it works great up about 2390:
1 ]=> (benchmark (delay (take 2390 primes)))
1 ]=> (benchmark (delay (take 2400 primes)))
;Aborting!: out of memory
;GC #58: took: 0.34 (51%) CPU time, 0.34 (50%) real time; free: 13773
;GC #59: took: 0.33 (100%) CPU time, 0.34 (99%) real time; free: 13806
This should be the fastest of all approaches but it just dies.
I’d really like some help here! Volunteers?
Also, it’s interesting to note that while the Rabin-Miller test performs better for testing huge numbers (only one at a time), but only by a small fraction (there is really not a big difference < 1sec for 32416190071, which is prime), both of these tests are ineffective against the simple divisor check fo smaller numbers as can be clearly seen from the benchmarks above. (disclaimer: I know that they are really really imprecise!)
Still, the prime? test beats rabinmiller-prime? by a factor of about 3 and rabinmiller-true-prime? by a factor of over 10.
What’s wrong? Clearly there are too many tests being done. For the biggest prime number I could fine on the Big Primes website, there are a total of 1250 = (* 2 (square (ceiling (log 32416190071)))) being performed.
You can find the code in my pp2011 github repo, here.
Update: I think however that there are major issues with the tests. Restarting Scheme makes a huge difference. Can anyone please confirm? I’m using MIT/GNU Scheme Release 9.0.1 under Mac OS X.
I would have liked to make a point about testing the primality of numbers in an efficient way. It turns out the point is: I suck at it!