# 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 [[3]]s 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)
m
(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)))))))
## Divisors

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)
'(())
(reverse
(let ((tails (recur (cdr lists))))
(let outer ((list (car lists)) (product '()))
(if (null? list)
product
(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)
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)
x
(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.
## Totients

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 *p*^{n} is (*p*-1)·*p*^{n}-1. The totient
function loops over the prime factors of the given number, performing the requisite
multiplications.