Skip to content

Testing for Primes: Scheme

5 March, 2011

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 A_k 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 n checks all odd numbers d starting from 3 to see whether or not they divide n. Obviously, a better approach would have been to do the test only for prime d, but there’s no way of figuring that out without first having all prime numbers less than n 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 n is a prime number and a < n,\, a \in \mathbb{N^*}, then a^n \equiv a\, (\mod n).

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 a and test the condition in Fermat’s Little Theorem. If n is prime, then the test is true for all values of a, but testing all n - 1 values of a would not be feasible.

Checking

(= (expmod a n n) a)

only proves that n is probably prime, but there is no way of knowing for sure in only one test. Picking a few random values of a 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)
;Value: #f
1 ]=> (fermat-prime? 9746347772161)
;Value: #t

In fact, 9746347772161 = 7 \times 11 \times 13 \times 17 \times 19 \times 31 \times 37 \times 41 \times 641.

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 a and checking to see if a^n \equiv a \mod n, it checks to see if a^{n-1} \equiv 1\, (\mod n).

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 p be a prime number. Clearly (-1)^2 = 1^2 \equiv 1\, (\mod p), for any p. Those are called the trivial square roots of 1. We’ll prove that if p is prime, there can be no others: suppose there is another x such that x^2 \equiv 1\, ( \mod p ) \Leftrightarrow x^2 - 1 = (x-1)(x+1) \equiv 0\, (\mod p), so p divides the product, therefore it must divide one of the factors, meaning either x - 1 \equiv 0 (\mod p) or x + 1 \equiv 0\, (\mod p) which means that either x = 1 or x = -1 (in the field \mathbb{Z}/p\mathbb{Z}).

If we at some point in the exponentiation get a partial result  (s (expmod2 base (/ exp 2) mod)) meaning that a^n\,(\mod n) = s^2 \times t\,(\mod n) and s^2 \equiv 1 (\mod n) without s being either 1 or -1 = n - 1, then by what we’ve proven above, n 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 2(\log n)^2:

(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))))
;Value: 7.3799999999999955
1 ]=> (benchmark (delay (take 5000 (filter-stream rabinmiller-true-prime? naturals))))
;Value: 97.03999999999999
1 ]=> (benchmark (delay (take 5000 (filter-stream prime? naturals))))
;Value: 1.8300000000000125
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)))
;Value: .00999999999999801
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! 😀

from → fun, questions, scheme
5 Comments leave one →
  1. Mihai Maruseac permalink
    6 March, 2011 8:42 AM

    You can try something to see if it helps. Make a function to intercept each force call and update a global variable (imperative way only). Check how many updates are done in each of the three variants.

    Something like

    (define oldforce force)
(define forces 0)
(define (force x) ((set! forces (+ 1 forces)) (oldforce x)))

    Also, I’ll try to test them on my computer to see what results I have.

    Reply
  2. Mihai Maruseac permalink
    6 March, 2011 9:42 AM

    I couldn’t test your code (no time left, sorry) but I tried my own implementation. See it here[1].

    Also, see the source of inspiration[2]

    [1]: http://swarm.cs.pub.ro/~mihai/facultate/pp/wh-primes.scm
    [2]: http://www.haskell.org/haskellwiki/Prime_numbers

    Reply
  3. danf permalink*
    6 March, 2011 11:42 AM

    I’m trying to run wh-primes.scm but my Scheme won’t load it at all.
    First of all, I needed to replace the unicode lambda with “lambda”, my foldl function is called fold-left but it still dies saying:


    1 ]=> (load "../wh-primes.scm")
    ;Loading "../wh-primes.scm"...
    ;The object 2, passed as the first argument to car, is not the correct type.

    I’m using MIT/GNU Scheme.

    Ok, actually, it loads (and with foldl) with mzcheme (the binary in the bin/ folder of the Racket installation).
    New problems:
    1. It sometimes loads the file and sometimes it doesn’t (although this always fails with MIT Scheme):

    > (load "wh-primes.scm")
    take: expected argument of type ; given (1 #)
    === context ===
    /Applications/Racket v5.0.2/collects/racket/list.rkt:102:0: take
    /Users/Dan/Dropbox/proj/2/pp/wh-primes.scm:7:0: take
    /Users/Dan/Dropbox/proj/2/pp/wh-primes.scm:49:0: next
    /Applications/Racket v5.0.2/collects/racket/private/list.rkt:179:4: foldl
    /Applications/Racket v5.0.2/collects/racket/private/misc.rkt:74:7
    > (load "wh-primes.scm")
    >

    2. It doesn’t have a (runtime) function. Any suggestions?

    Also, why are these so different? 😐

    Reply
  4. Mithrandir permalink
    6 March, 2011 11:52 AM

    Try using drScheme or drRacket. Scheme is known to have spawned multiple dialects and each tool recognizes only a few of them.

    Reply
  5. danf permalink*
    6 March, 2011 12:05 PM

    Okay, so after setting up a symlink and reloading my .profile, I can now use mzscheme instead of mit-scheme in my Emacs 😀

    Also, there seems to be a (time ) function available in mzscheme, so I’ve modified my “benchmark” a bit:


    (define benchmark
    (lambda (f)
    (time (begin
    (force f)
    '()))))

    And now, I can finally say, victory!:

    > (benchmark (delay (take 50000 primes-wh)))
    cpu time: 197 real time: 203 gc time: 160
    ()

    I’ll play with the other tests as well later, but this wheel approach seems to be the best. I also found something about it in a paper (which I can’t immediately find in my history) where it compared different approaches (in Haskell) – it tried trial-division, the sieving approach, the wheel and then tried doing some stuff with maps or heaps instead of lists.
    I’ll post a link as soon as I find it 🙂

    Reply

Leave a comment