Prime Numbers

Prime numbers are natural numbers with exactly two distinct natural number divisors: itself and 1. There are an infinite number of prime numbers, starting with 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, and 97; by definition, 1 is not prime. Prime numbers appear frequently in the study of number theory and have applications in public-key cryptography as well as many other facets of mathematics.

The Greek mathematician Eratosthenes (276BC-194BC), who was born in Cyrene (modern-day Libya) and lived most of his life in Alexandria, Egypt, and who most famously calculated the circumference of the Earth to within an accuracy of 1%, invented a method of calculating the prime numbers: make an ascending list of integers starting from two; then, repeatedly take the first unmarked number on the list, which is prime, and cross off all its multiples, which are not prime. This algorithm has become known as the Sieve of Eratosthenes, since the repeated marking of multiples acts like a sieve that blocks composite numbers and allows only primes to fall through.

Oleg Kiselyov used the same algorithm to create the (primes n) function that he posted to comp.lang.scheme on November 13, 2001, with one simple optimization: since the only even prime number is two, he sieved only odd numbers:

(define (primes n)
  (let* ((max-index (quotient (- n 3) 2))
         (v (make-vector (+ 1 max-index) #t)))
    (let loop ((i 0) (primes '(2)))
      (cond ((< max-index i) (reverse primes))
            ((vector-ref v i)
              (let ((prime (+ i i 3)))
                (do ((j (+ 3 (* 3 i)) (+ j prime)))
                    ((< max-index j))
                  (vector-set! v j #f))
              (loop (+ 1 i) (cons prime primes))))
            (else (loop (+ 1 i) primes))))))

The code is quite straight forward. A vector v is initialized to #t, then the let loop runs through each number, and when it finds a prime, the do loop crosses off all its multiples. The mysterious 3s sprinkled through the code implement the odd-only optimization.

It was another Greek mathematician, Euclid, the "Father of Geometry" (who flourished circa 300BC, and also lived in Alexandria), who proved that there are an infinite number of prime numbers. In his book Elements, the most successful mathematics textbook ever written, he described an algorithm for finding the greatest common divisor of two numbers:

(define (gcd m n)
  (if (zero? n)
      (gcd n (modulo m n))))

This simple algorithm, based on repeated subtraction via the modulo function, is the oldest algorithm still in use today; Donald Knuth (AoCP1) calls it the "grand-daddy" of all algorithms. Although it bears his name, Euclid almost certainly did not invent the algorithm; he gets the credit because he was the first to write it down. Euclid's Algorithm is important to us because it allows us to determine if two numbers are coprime (relatively prime); that is, they have no common divisors except unity, so (gcd m n) is 1.

The standard test for primality of a given number is trial division by all prime numbers up to the square root of the given number. But in the same way that a simple optimization of the Sieve of Eratosthenes is to sieve only odd numbers, coprimes provide a basis for an improved primality test that doesn't require a pre-existing list of prime numbers. The improvement extends the simple odd-number optimization by considering additional prime numbers beyond two. For instance, once all the multiples of two are eliminated, only half the numbers remain. If multiples of three are eliminated from the remaining half, only one-third of the numbers remain. And so on. This makes it faster to determine by trial division if a number is prime, since there are fewer numbers remaining.

One way this is done is by wheel factorization. Consider the image of a wheel rolling along a number line; the wheel's circumference is the width of six numbers, and has six points numbered zero through five. Spokes extend from the wheel at points one and five. The wheel is positioned with the spoke at point five pointing at the five on the number line. As the wheel rolls along the number line, point zero aligns with number six, then the spoke at point one aligns with the number seven, then points two, three, and four align with numbers eight, nine, and ten, the spoke at point five aligns with the number eleven, and so on.

The wheel works by pointing to candidate primes: five, seven, and eleven are all candidate primes, but six, eight, nine and ten are certainly not prime. Candidate primes must still be tested; for instance, 25 will be a candidate prime. But such testing is quicker because so many non-prime numbers are identified by the wheel.

The size of the wheel, the points where spokes are located, and the starting position of the wheel on the number line are all special, and are based on coprimes. The wheel described above is known as a 2,3-wheel because it is based on the starting prime numbers two and three. Its circumference is six because 2×3=6. Spokes are at points one and five because those are the only two numbers coprime to six. The wheel starts with the last spoke at the next prime number after the highest prime size of the wheel. As a further example, a 2,3,5-wheel has a circumference of 30, spokes at 1, 7, 11, 13, 17, 19, 23, and 29, and starts with the 29 spoke positioned at 7 on the number line.

Wheels quickly reach a point of diminishing returns: a 2-wheel reduces the number of candidate primes by 50%, a 2,3-wheel by 67%, and a 2,3,5,7-wheel by 77%, but to get a 90% reduction requires the primes up to 257, a 95% reduction requires the primes up to 75037, and a 96% reduction requires the primes up to 1246379. The following test for primality is based on a 2,3,5,7-wheel, which is a convenient balance between speed and code-size:

(define (prime? n)
  (define (wheel . xs) (set-cdr! (last-pair xs) xs) xs)
  (if (< n 2) #f
      (let loop ((i 2) (z (cons 1 (cons 2 (cons 2 (cons 4 (wheel
          2 4 2 4 6 2 6 4 2 4 6 6 2 6 4 2 6 4 6 8 4 2 4 2 4 8
          6 4 6 2 4 6 2 6 6 4 2 4 6 2 6 4 2 4 2 10 2 10)))))))
        (cond ((< n (* i i)) #t)
              ((zero? (modulo n i)) #f)
              (else (loop (+ i (car z)) (cdr z)))))))

The internal function wheel creates a circular list by tying the last cdr of a list back to its head. The special cases that start the wheel are added to the front of the wheel with cons. The wheel is specified by the distances between the spokes rather than the actual points of the spokes, so the sequence, after starting with 2, 3, 5, 7, and 11, continues with 11+2=13, 13+4=17, 17+2=19, 19+4=23, 23+6=29, 29+2=31, 31+6=37, and so on. A similar technique lets us compute the prime factors of a number:

(define (factors n)
  (define (wheel . xs) (set-cdr! (last-pair xs) xs) xs)
  (if (< n 2) '()
    (let loop ((n n) (i 2) (f '()) (z (cons 1 (cons 2 (cons 2 (cons 4
        (wheel 2 4 2 4 6 2 6 4 2 4 6 6 2 6 4 2 6 4 6 8 4 2 4 2 4 8
               6 4 6 2 4 6 2 6 6 4 2 4 6 2 6 4 2 4 2 10 2 10)))))))
      (cond ((< n (* i i)) (reverse (cons n f)))
            ((zero? (modulo n i))
              (loop (quotient n i) i (cons i f) z))
            (else (loop n (+ i (car z)) f (cdr z)))))))


Closely related to factors are divisors, those numbers which evenly divide a given number. For instance, the divisors of 28 are 1, 2, 4, 7, 14, and 28. To compute the divisors of a number first compute its factors, then calculate the products of the cross-products of all the powers of its factors:

(define (divisors n)
  (sort <
    (map (lambda (x) (apply * x))
      (apply cross
        (map pow-list
          (cluster =
            (factors n)))))))

This requires a few general-purpose utility functions; cross forms the cross-product of a set of lists, cluster groups and counts factors (for instance, cluster (factors 360) returns ((2 . 3) (3 . 2) (5 . 1))), and pow-list generates a list of the powers of a factor-cluster (for instance, (pow-list (2 . 3)) returns (8 4 2 1).

(define (cross . lists)
  (let recur ((lists lists))
    (if (null? lists)
          (let ((tails (recur (cdr lists))))
            (let outer ((list (car lists)) (product '()))
              (if (null? list)
                  (let ((item (car list)) (more (cdr list)))
                    (let inner ((tails tails) (product product))
                      (if (null? tails)
                          (outer more product)
                          (inner (cdr tails) (cons (cons item (car tails)) product))))))))))))

(define (cluster eql? xs)
  (if (null? xs)
      (let loop ((xs (cdr xs)) (prev (car xs)) (k 1) (result '()))
        (cond ((null? xs) (reverse (cons (cons prev k) result)))
              ((eql? (car xs) prev) (loop (cdr xs) prev (+ k 1) result))
              (else (loop (cdr xs) (car xs) 1 (cons (cons prev k) result))))))

(define (pow-list nk)
  (let ((n (car nk)) (k (cdr nk)))
    (let loop ((k k) (x '(1)))
      (if (zero? k)
          (loop (- k 1) (cons (* n (car x)) x))))))

Sometimes only the count of the divisors, or their sum, is needed:

(define (numdiv n)
  (apply *
    (map (lambda (x) (+ x 1))
      (map cdr
        (cluster =
          (factors n))))))

(define (sumdiv n)
  (apply +
    (map (lambda (x) (apply * x))
      (apply cross
        (map pow-list
          (cluster =
            (factors n)))))))

Divisors, numdiv and sumdiv are defined to include the given number; for instance, 28 is a divisor of 28. Sometimes, the need is for only the proper divisors of a number, excluding the number itself. Given the functions above, it is easy to exclude the given number.


Another function related to prime numbers is Leonhard Euler's totient function, which calculates how many numbers less than a given number are coprime to it. For instance, the totient of 30 is 8, and the eight totatives are 1, 7, 11, 13, 17, 19, 23, and 29; you might recognize those numbers as the spokes on a 2,3,5-wheel, which can be calculated as (totatives (* 2 3 5)):

(define (totatives n)
  (let loop ((i n) (ts '()))
    (cond ((zero? i) ts)
          ((= (gcd i n) 1) (loop (- i 1) (cons i ts)))
          (else (loop (- i 1) ts)))))

The totient function calculates the number of totatives of a given number:

(define (totient n)
  (let loop ((fs (factors n)) (prev 0) (ts '()))
    (cond ((null? fs) (apply * ts))
          ((= (car fs) prev) (loop (cdr fs) (car fs) (cons (car fs) ts)))
          (else (loop (cdr fs) (car fs) (cons (- (car fs) 1) ts))))))

Totients are multiplicative, as long as the two multipliers are coprime, so the totient function calculates the prime factors of the given number and multiplies the totients of each grouping of primes. For instance, (factors 12) is (2 2 3), so the totient of 12 is the totient of (* 2 2) times the totient of 3. The totient of a prime number like 3 is one less than the number, and the totient of a power of a prime number like pn is (p-1)·pn-1. The totient function loops over the prime factors of the given number, performing the requisite multiplications.