Determine If Number Is Prime

The SICP book is still being read...

Now I am going to introduce some algorithms that allow us to say whether the number is prime or not. Again we will start from the simplest one.

Take numbers 1 by 1 starting from 2 to √n and check, if it divides n. The complexity of the algorithm is Θ(√n).

(define (prime? n)
  (define (divides? k)
    (= (remainder n k) 0))
  (define (find-divisor n test)
    (cond ((> (sqr test) n) n)
          ((divides? test) test)
          (else (find-divisor n (+ test 1)))))

  (= n (find-divisor n 2)))

Fermat Primality Test

There is a test for primality with Θ(log(n)) complexity. It is based on Fermat little theorem.

an ≡ a (mod n), where 1 ≤ a < n

Which means, that if n is prime, then for any integer a, an − a will be evenly divisible by n.

;this procedure calculates baseexp mod m without using big numbers
(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp) (remainder (sqr (expmod base (/ exp 2) m)) m))
        (else (remainder (* base (expmod base (- exp 1) m)) m))))

(define (fermat-test n)
  (define (try-it a)
    (= (expmod a n n) a))
  (try-it (+ 1 (random (- n 1)))))

Unfortunately there are some issues with this test - it is probabilistic. If the test fails, n is certainly composite. But if it passes, it does not guarantee that n is prime, although it shows high probability for n to be. So we can run as many tests as we need to make the mistake probability as low as possible.

(define (fermat-prime? n times)
  (cond ((= times 0) true)
        ((fermat-test n) (fermat-prime? n (- times 1)))
        (else false)))

Unfortunately that is not all. There are some numbers that pass the test (for every a!), without being prime. These are called Carmichael numbers. So to determine, if n is prime, we need to check if it passes the Fermat test and is not a Carmichael number, or use a similar Miller–Rabin primality test, that cannot be tricked.

Miller–Rabin Primality Test

It starts from an alternate form of Fermat's Little Theorem:

an-1 ≡ 1 (mod n), where 1 ≤ a < n

The difference is that whenever we perform the squaring step in expmod, we check to see if we have discovered a "nontrivial square root of 1 modulo p", that is, a number not equal to 1 or n - 1 whose square is equal to 1 mod n. It is possible to prove that if such a nontrivial square root of 1 exists, then n is not prime. It is also possible to prove that if n is an odd number that is not prime, then, for at least half the numbers a < n, computing an-1 in this way will reveal a nontrivial square root of 1 mod n.

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

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

(define (expmod base exp m)
  (define (!= x y) (not (= x y)))

  (cond ((= exp 0) 1)
        ((and (!= base 1) 
              (!= base (- m 1)) 
              (= (remainder( sqr base) m) 1)) 0)
        ((even? exp) (remainder (sqr (expmod base (/ exp 2) m)) m))
        (else (remainder (* base (expmod base (- exp 1) m)) m))))