From Mathematics to Generic Programming (2011)

---

Unfortunately, Pythagoras and his followers kept their work secret, so none of their writings survive. However, we know from contemporaries what some of his discoveries were. Some of these come from a first-century book called Introduction to Arithmetic by Nicomachus of Gerasa. These included observations about geometric properties of numbers; they associated numbers with particular shapes.

Triangular numbers, for example, which are formed by stacking rows representing the first n integers, are those that formed the following geometric pattern:

Oblong numbers are those that look like this:

It is easy to see that the nth oblong number is represented by an n × (n + 1) rectangle:

_n = n(n + 1)

It’s also clear geometrically that each oblong number is twice its corresponding triangular number. Since we already know that triangular numbers are the sum of the first n integers, we have

So the geometric representation gives us the formula for the sum of the first n integers:

Another geometric observation is that the sequence of odd numbers forms the shape of what the Greeks called gnomons (the Greek word for a carpenter’s square; a gnomon is also the part of a sundial that casts the shadow):

Combining the first n gnomons creates a familiar shape—a square:

This picture also gives us a formula for the sum of the first n odd numbers:

Exercise 3.1. Find a geometric proof for the following: take any triangular number, multiply it by 8, and add 1. The result is a square number. (This problem comes from Plutarch’s Platonic Questions.)

3.2 Sifting Primes

Pythagoreans also observed that some numbers could not be made into any nontrivial rectangular shape (a shape where both sides of the rectangle are greater than 1). These are what we now call prime numbers—numbers that are not products of smaller numbers:

2, 3, 5, 7, 11, 13,...

(“Numbers” for the Greeks were always whole numbers.) Some of the earliest observations about primes come from Euclid. While he is usually associated with geometry, several books of Euclid’s Elements actually discuss what we now call number theory. One of his results is this theorem:

Theorem 3.1 (Euclid VII, 32): Any number is either prime or divisible by some prime.

The proof, which uses a technique called “impossibility of infinite descent,” goes like this:¹

¹ Euclid’s proof of VII, 32 actually relies on his proposition VII, 31 (any composite number is divisible by some prime), which contains the reasoning shown here.

Proof. Consider a number A. If it is prime, then we are done. If it is composite (i.e., nonprime), then it must be divisible by some smaller number B. If B is prime, we are done (because if A is divisible by B and B is prime, then A is divisible by a prime). If B is composite, then it must be divisible by some smaller number C, and so on. Eventually, we will find a prime or, as Euclid remarks in his proof of the previous proposition, “an infinite sequence of numbers will divide the number, each of which is less than the other; and this is impossible.”

This Euclidean principle that any descending sequence of natural numbers terminates is equivalent to the induction axiom of natural numbers, which we will encounter in Chapter 9.

* * *

Another result, which some consider the most beautiful theorem in mathematics, is the fact that there are infinitely many primes:

Theorem 3.2 (Euclid IX, 20): For any sequence of primes {p₁,...,p_n}, there is a prime p not in the sequence.

Proof. Consider the number

where p_i is the ith prime in the sequence. Because of the way we constructed q, we know it is not divisible by any p_i. Then either q is prime, in which case it is itself a prime not in the sequence, or q is divisible by some new prime, which by definition is not in the sequence. Therefore, there are infinitely many primes.

One of the best-known techniques for finding primes is the Sieve of Eratosthenes. Eratosthenes was a 3rd-century Greek mathematician who is remembered in part for his amazingly accurate measurement of the circumference of the Earth. Metaphorically, the idea of Eratosthenes’ sieve is to “sift” all the numbers so that the nonprimes “fall through” the sieve and the primes remain at the end. The actual procedure is to start with a list of all the candidate numbers and then cross out the ones known not to be primes (since they are multiples of primes found so far); whatever is left are the primes. Today the Sieve of Eratosthenes is often shown starting with all positive integers up to a given number, but Eratosthenes already knew that even numbers were not prime, so he didn’t bother to include them.

Following Eratosthenes’ convention, we’ll also include only odd numbers, so our sieve will find primes greater than 2. Each value in the sieve is a candidate prime up to whatever value we care about. If we want to find primes up to a maximum of m = 53, our sieve initially looks like this:

In each iteration, we take the first number (which must be a prime) and cross out all the multiples except itself that have not previously been crossed out. We’ll highlight the numbers being crossed out in the current iteration by boxing them. Here’s what the sieve looks like after we cross out the multiples of 3:

Next we cross out the multiples of 5 that have not yet been crossed out:

And then the remaining multiples of 7:

We need to repeat this process until we’ve crossed out all the multiples of factors less than or equal to , where m is the highest candidate we’re considering.

In our example, m = 53, so we are done. All the numbers that have not been crossed out are primes:

Before we write our implementation of the algorithm, we’ll make a few observations. Let’s go back to what the sieve looked like in the middle of the process (say, when we were crossing out multiples of 5) and add some information—namely, the index, or position in the list, of each candidate being considered:

Notice that when we’re considering multiples of factor 5, the step size—the number of entries between two numbers being crossed out, such as 25 and 35—is 5, the same as the factor. Another way to say this is that the difference between the indexes of any two candidates being crossed out in a given iteration is the same as the factor being used. Also, since the list of candidates contains only odd numbers, the difference between two values is twice as much as the difference between two indexes. So the difference between two numbers being crossed out in a given iteration (e.g., between 25 and 35) is twice the step size or, equivalently, twice the factor being used. You’ll see that this pattern holds for all the factors we considered in our example as well.

Finally, we observe that the first number crossed out in each iteration is the square of the prime factor being used. That is, when we’re crossing out multiples of 5, the first one that wasn’t previously crossed out is 25. This is because all the other multiples were already accounted for by previous primes.

3.3 Implementing and Optimizing the Code

At first glance it seems like our algorithm will need to maintain two arrays: one containing the candidate numbers we’re sifting—the “values”—and another containing Boolean flags indicating whether the corresponding number is still there or has been crossed out. However, after a bit of thought it becomes clear that we don’t actually need to store the values at all. Most of the values (namely, all the nonprimes) are never used. When we do need a value, we can compute it from its position; we know that the first value is 3 and that each successive value is 2 more than the previous one, so the ith value is 2i + 3.

So our implementation will store just the Boolean flags in the sieve, using true for prime and false for composite. We call the process of “crossing out” nonprimes marking the sieve. Here’s a function we’ll use to mark all the nonprimes for a given factor:

Click here to view code image

template <RandomAccessIterator I, Integer N>
void mark_sieve(I first, I last, N factor) {
// assert(first != last)
*first = false;
while (last - first > factor) {
first = first + factor;
*first = false;
}
}

We are using the convention of “declaring” our template arguments with a description of their requirements. We will discuss these requirements, known as concepts, in detail later on in Chapter 10; for now, readers can consult Appendix C as a reference. (If you are not familiar with C++ templates, these are also explained in this appendix.)

As we’ll see shortly, we’ll call this function with first pointing to the Boolean value corresponding to the first “uncrossed-out” multiple of factor, which as we saw is always factor’s square. For last, we’ll follow the STL convention of passing an iterator that points just past the last element in our table, so that last - first is the number of elements.

* * *

Before we see how to sift, we observe the following sifting lemmas:

• The square of the smallest prime factor of a composite number c is less than or equal to c.

• Any composite number less than p² is sifted by (i.e., crossed out as a multiple of) a prime less than p.

• When sifting by p, start marking at p².

• If we want to sift numbers up to m, stop sifting when p² > m.

We will use the following formulas in our computation:

step between multiple k and multiple k + 1 of value at i:

index of square of value at i:

We can now make our first attempt at implementing the sieve:

Click here to view code image

template <RandomAccessIterator I, Integer N>
void sift0(I first, N n) {
std::fill(first, first + n, true);
N i(0);
N index_square(3);
while (index_square < n) {
// invariant: index_square = 2i^2 + 6i + 3
if (first[i]) { // if candidate is prime
mark_sieve(first + index_square,
first + n, // last
i + i + 3); // factor
}
++i;
index_square = 2*i*(i + 3) + 3;
}
}

It might seem that we should pass in a reference to a data structure containing the Boolean sequence, since the sieve works only if we sift the whole thing. But by instead passing an iterator to the beginning of the range, together with its length, we don’t constrain which kind of data structure to use. The data could be in an STL container or in a block of memory; we don’t need to know. Note that we use the size of the table n rather than the maximum value to sift m.

The variable index_square is the index of the first value we want to mark—that is, the square of the current factor. One thing we notice is that we’re computing the factor we use to mark the sieve (i + i + 3) and other quantities (shown in slanted text) every time through the loop. We can hoist common subexpressions out of the loop; the changes are shown in bold:

Click here to view code image

template <RandomAccessIterator I, Integer N>
void sift1(I first, N n) {
I last = first + n;
std::fill(first, last, true);
N i(0);
N index_square(3);
N factor(3);
while (index_square < n) {
// invariant: index_square = 2i^2 + 6i + 3,
// factor = 2i + 3
if (first[i]) {
mark_sieve(first + index_square, last, factor);
}
++i;
factor = i + i + 3;
index_square = 2*i*(i + 3) + 3;
}
}

The astute reader will notice that the factor computation is actually slightly worse than before, since it happens every time through the loop, not just on iterations when the if test is true. However, we shall see later why making factor a separate variable makes sense. A bigger issue is that we still have a relatively expensive operation—the computation of index_square, which involves two multiplications. So we will take a cue from compiler optimization and use a technique known as strength reduction, which was designed to replace more expensive operations like multiplication with equivalent code that uses less expensive operations like addition.² If a compiler can do this automatically, we can certainly do it manually.

² While multiplication is not necessarily slower than addition on modern processors, the general technique can still lead to using fewer operations.

Let’s look at these computations in more detail. Suppose we replaced

Click here to view code image

factor = i + i + 3;
index_square = 3 + 2*i*(i + 3);

with

factor += δ_factor;

index_square += δ_{index_square};

where δ_factor and δ_{index_square} are the differences between successive (ith and i + 1st) values of factor and index_square, respectively:

δ_factor is easy; the variables cancel and we get the constant 2. But how did we simplify the expression for δ_{index_square}? We observe that by rearranging the terms, we can express it using something we already have, factor(i), and something we need to compute anyway, factor(i + 1). (When you know you need to compute multiple quantities, it’s useful to see if one can be computed in terms of another. This might allow you to do less work.)

With these substitutions, we get our final version of sift; again, our improvements are shown in bold:

Click here to view code image

template <RandomAccessIterator I, Integer N>
void sift(I first, N n) {
I last = first + n;
std::fill(first, last, true);
N i(0);
N index_square(3);
N factor(3);
while (index_square < n) {
// invariant: index_square = 2i^2 + 6i + 3,
// factor = 2i + 3
if (first[i]) {
mark_sieve(first + index_square, last, factor);
}
++i;
index_square += factor;
factor += N(2);
index_square += factor;
}
}

Exercise 3.2. Time the sieve using different data sizes: bit (using std::vector<bool>), uint8_t, uint16_t, uint32_t, uint64_t.

Exercise 3.3. Using the sieve, graph the function

π(n) = number of primes < n

for n up to 10⁷ and find its analytic approximation.

We call primes that read the same backward and forward palindromic primes. Here we’ve highlighted the ones up to 1000:

Interestingly, there are no palindromic primes between 1000 and 2000:

Exercise 3.4. Are there palindromic primes > 1000? What is the reason for the lack of them in the interval [1000, 2000]? What happens if we change our base to 16? To an arbitrary n?

3.4 Perfect Numbers

As we saw in Section 3.1, the ancient Greeks were interested in all sorts of properties of numbers. One idea they came up with was that of a perfect number—a number that is the sum of its proper divisors.³ They knew of four perfect numbers:

³ A proper divisor of a number n is a divisor of n other than n itself.

Perfect numbers were believed to be related to nature and the structure of the universe. For example, the number 28 was the number of days in the lunar cycle.

What the Greeks really wanted to know was whether there was a way to predict other perfect numbers. They looked at the prime factorizations of the perfect numbers they knew:

and noticed the following pattern:

The result of this expression is perfect when the the second term is prime. It was Euclid who presented the proof of this fact around 300 BC.

Theorem 3.3 (Euclid IX, 36):