From Mathematics to Generic Programming (2011)

---

Counted ranges are easier for the compiler, because it knows the number of iterations, called the trip count, in advance. This allows the compiler to make certain optimizations, such as loop unrolling or software pipelining.

Exercise 11.6. Why don’t we provide the following interface?

Click here to view code image

pair<I0, I1> swap_ranges_n(I0 first0, I1 first1, N0 n0, N1 n1)

11.3 Rotation

One of the most important algorithms you’ve probably never heard of is rotate. It’s a fundamental tool used behind the scenes in many common computing applications, such as buffer manipulation in text editors. It implements the mathematical operation rotation.

Definition 11.5. A permutation of n elements by k where k ≥ 0:

(k mod n, k + 1 mod n, ..., k + n – 2 mod n, k + n – 1 mod n)

is called an n by k rotation.

If you imagine all n elements laid out in a circle, we’re shifting each one “clockwise” by k.

At first glance, it might seem that rotation could be implemented with a modular shift, taking the beginning and end of the range, together with the amount to shift, as arguments. However, doing modular arithmetic on every operation is quite expensive. Also, it turns out that rotation is equivalent to exchanging different length blocks, a task that is extremely useful for many applications. Viewed this way, it is convenient to present rotation with three iterators: f,m, and l, where [f,m) and [m,l) are valid ranges.⁴ Rotation then interchanges ranges [f,m) and [m,l). If the client wants to rotate k positions in a range [f,l), then it should pass a value of m equal to l – k. As an example, if we want to do a rotation with k = 5 on a sequence specified by the range [0, 7), we choose m = l – k = 7 – 5 = 2:

⁴ The names f,m, and l are meant to be mnemonics for first, middle, and last.

0 1 2 3 4 5 6
f m l

which produces this result:

2 3 4 5 6 0 1

Essentially, we’re moving each item 5 spaces to the right (and wrapping around when we run off the end).

Exercise 11.7. Prove that if we do rotate(f, m, l), it performs distance(f, l) by distance(m, l) rotation.

An important rotate algorithm was developed by David Gries, a professor at Cornell University, together with IBM scientist Harlan Mills:

Click here to view code image

template <ForwardIterator I>
void gries_mills_rotate(I f, I m, I l) {
// u = distance(f, m) && v = distance(m, l)
if (f == m || m == l) return; // u == 0 || v == 0
pair<I, I> p = swap_ranges(f, m, m, l);
while(p.first != m || p.second != l) {
if (p.first == m) { // u < v
f = m; m = p.second; // v = v - u
} else { // v < u
f = p.first; // u = u - v
}
p = swap_ranges(f, m, m, l);
}
return; // u == v
}

The algorithm first uses the regular swap_ranges function to swap as many elements as we can—as many elements as there are in the shorter range. The if test checks whether we’ve exhausted the first range or the second range. Depending on the result, we reset the start positions of f and m. Then we do another swap, and repeat the process until neither range has remaining elements.

This is easier to follow with an example. Let’s look at how a range is transformed, and how the iterators move, as the algorithm runs:

Start:

0 1 2 3 4 5 6
f m l

Swap [0, 1] and [2, 3]. Have we exhausted both ranges? No, only the first one, so we set f = m and m = p.second, which points to the first element in the sequence that hasn’t yet been moved:

2 3 0 1 4 5 6
f m l

Swap [0, 1] and [4, 5]. Have we exhausted both ranges? No, only the first one, so again we set f = m and m = p.second:

2 3 4 5 0 1 6
f m l

Swap [0] with [6]. Have we exhausted both ranges? No, only the second one this time, so we set f = p.first:

2 3 4 5 6 1 0
f m l

Swap [1] with [0]. Have we exhausted both ranges? Yes, so we’re done.

2 3 4 5 6 0 1
f m
l

Now look at the comments in the gries_mills_rotate code (shown in boldface). We call u the length of the first range [f,m), and v the length of the second range [m,l). We can observe something remarkable: the annotations are performing our familiar subtractive GCD! At the end of the algorithm, u = v = GCD of the lengths of the initial two ranges.

Exercise 11.8. If you examine swap_ranges, you will see that the algorithm does unnecessary iterator comparisons. Rewrite the algorithm so that no unnecessary iterator comparisons are done.

It turns out that many applications benefit if the rotate algorithm returns a new middle: a position where the first element moved. If rotate returns this new middle, then rotate(f, rotate(f, m, l), l) is an identity permutation. First, we need the following “auxiliary rotate” algorithm:

Click here to view code image

template <ForwardIterator I>
void rotate_unguarded(I f, I m, I l) {
// assert(f != m && m != l)
pair<I, I> p = swap_ranges(f, m, m, l);
while (p.first != m || p.second != l) {
f = p.first;
if (m == f) m = p.second;
p = swap_ranges(f, m, m, l);
}
}

The central loop is the same as in the Gries-Mills algorithm, just written differently. (We could have written it this way before but wanted to make the u and v computations clearer.)

We need to find m′—the element whose distance to the last is the same as m’s distance from the first. It’s the value returned by the first call to swap_ranges. To get it back, we can embed a call to rotate_unguarded in our final version of rotate, which works as long as we have forward iterators. As we’ll explain shortly, the forward_iterator_tag type in the argument list will help us invoke this function only in this case:

Click here to view code image

template <ForwardIterator I>
I rotate(I f, I m, I l, std::forward_iterator_tag) {
if (f == m) return l;
if (m == l) return f;
pair<I, I> p = swap_ranges(f, m, m, l);
while (p.first != m || p.second != l) {
if (p.second == l) {
rotate_unguarded(p.first, m, l);
return p.first;
}
f = m;
m = p.second;
p = swap_ranges(f, m, m, l);
}
return m;
}

How much work does this algorithm do? Until the last iteration of the main loop, every swap puts one element in the right place and moves another element out of the way. But in the final call to swap_ranges, the two ranges are the same length, so every swap puts both elements it is swapping into their final positions. In essence, we get an extra move for free on every swap. The total number of swaps is the total number of elements n, minus the number of free swaps we saved in the last step. How many swaps did we save in the last step? The length of the ranges at the end, as we saw earlier, is gcd(n – k,k) = gcd(n,k), where n = distance(f,l) and k = distance(m,l). So the total number of swaps is n – gcd(n,k). Also, since each swap takes three assignments (tmp = a; a = b; b = tmp), the total number of assignments is 3 (n – gcd(n,k)).

11.4 Using Cycles

Can we find a faster rotate algorithm? We can if we exploit the fact that rotations, like any permutations, have cycles. Consider a rotation of k = 2 for n = 6 elements:

The item in position 0 ends up in position 2, 2 ends up in position 4, and 4 ends up back in position 0. These three elements form a cycle. Similarly, item 1 ends up in 3, 3 in 5, and 5 back in 1, forming another cycle. So this rotation has two cycles. Recall from Section 11.1 that we can perform any permutation in n–u+v assignments, where n is the number of elements, u is the number of trivial cycles, and v is the number of nontrivial cycles. Since we normally don’t have any trivial cycles, we need n + v assignments.

Exercise 11.9. Prove that if a rotation of n elements has a trivial cycle, then it has n trivial cycles. (In other words, a rotation either moves all elements or no elements.)

It turns out that the number of cycles is gcd(k,n), so we should be able to do the rotation in n + gcd(k,n) assignments,⁵ instead of the 3(n – gcd(n,k)) we needed for the Gries-Mills algorithm. Furthermore, in practice GCD is very small; in fact, it is 1 (that is, there is only one cycle) about 60% of the time. So a rotate algorithm that exploits cycles always does fewer assignments.

⁵ See Elements of Programming, Section 10.4, for the proof.

There is one catch: Gries-Mills only required moving one step forward; it works even for singly linked lists. But if we want to take advantage of cycles, we need to be able to do long jumps. Such an algorithm requires stronger requirements on the iterators—namely, the ability to do random access.

To create our new rotate function, we’ll first write a helper function that moves every element in a cycle to its next position. But instead of saying “which position does the item in position x move to,” we’ll say “in which position do we find the item that’s going to move to position x.” Even though these two operations are symmetric mathematically, it turns out that the latter is more efficient, since it needs only one saved temporary variable per cycle, instead of one for every item that needs to be moved (except the last).

Here’s our helper function:

Click here to view code image

template <ForwardIterator I, Transformation F>
void rotate_cycle_from(I i, F from) {
ValueType<I> tmp = *i;
I start = i;
for (I j = from(i); j != start; j = from(j)) {
*i = *j;
i = j;
}
*i = tmp;
}

Note that we’re using the ValueType type function we defined near the end of Section 10.8.

How does rotate_cycle_from know which position an item comes from? That information will be encapsulated in a function object from that we pass in as an argument. You can think of from(i) as “compute where the element moving into position i comes from.”

The function object we’re going to pass to rotate_cycle_from will be an instance of rotate_transform:

Click here to view code image

template <RandomAccessIterator I>
struct rotate_transform {
DifferenceType<I> plus;
DifferenceType<I> minus;
I m1;

rotate_transform(I f, I m, I l) :
plus(m - f), minus(m - l), m1(f + (l - m)) {}
// m1 separates items moving forward and backward

I operator()(I i) const {
return i + ((i < m1) ? plus : minus);
}
};

The idea is that even though we are conceptually “rotating” elements, in practice some items move forward and some move backward (because the rotation caused them to wrap around the end of our range). When rotate_transform is instantiated for a given set of ranges, it precomputes (1) how much to move forward for items that should move forward, (2) how much to move backward for things that move backward, and (3) what the crossover point is for deciding when to move forward and when to move backward.

Now we can write the cycle-exploiting version of algorithm for rotation, which is a variation of the algorithm discovered by Fletcher and Silver in 1965:

Click here to view code image

template <RandomAccessIterator I>
I rotate(I f, I m, I l, std::random_access_iterator_tag) {
if (f == m) return l;
if (m == l) return f;
DifferenceType<I> cycles = gcd(m - f, l - m);
rotate_transform<I> rotator(f, m, l);
while (cycles-- > 0) rotate_cycle_from(f + cycles, rotator);
return rotator.m1;
}

After handling some trivial boundary cases, the algorithm first computes the number of cycles (the GCD) and constructs a rotate_transform object. Then it calls rotate_cycle_from to shift all the elements along each cycle, and repeats this for every cycle.

Let’s look at an example. Consider the rotation k = 2 for n = 6 elements that we used at the beginning of this section. For simplicity, we’ll assume that our values are integers stored in an array:

0 1 2 3 4 5

We also assume our iterators are integer offsets in an array, starting at 0. (Be careful to distinguish between the values at a position and the position itself.) To perform a k = 2 rotation, we’ll need to pass the three iterators f = 0, m = 4, and l = 6:

0 1 2 3 4 5
f m l

The boundary cases of our new rotate algorithm don’t apply, so the first thing it does is compute the number of cycles, which is equal to gcd(m – f, l – m) = gcd(4, 2) = 2. Then it constructs the rotator object, initializing its state variables as follows:

plus ← m – f = 4 – 0 = 4
minus ← m – l = 4 – 6 = – 2
m1 ← f + (l – m) = 0 + (6 – 4) = 2

The main loop of the function rotates all elements of a cycle, then moves on to the next cycle. Let us see what happens when rotate_cycle_from is called.

Initially, we pass f + d = 0 + 2 = 2 as the first argument. So inside the function, i = 2. We save the value at position 2, which is also 2, to our tmp variable and set start to our starting position of 2.

Now we go through a loop as long as a new variable, j, is not equal to start. Each time through the loop, we are going to set j by using the rotator function object that we passed in through the variable from. Basically, all that object does is add the stored values plus or minus to its argument, depending on whether the argument is less than the stored value m1. For example, if we call from(0), it will return 0 + 4, or 4, since 0 is less than 2. If we call from(4), it will return 4 + (–2), or 2, since 4 is not less than 2.

Here’s how the values in our array change as we go through the loop in rotate_cycle_from:

i ← 2, j ← from(2) = 0

0 1 2 3 4 5
j i

*i ← *j

0 1 0 3 4 5
j i

i ← j = 0, j ← from(0) = 4

0 1 0 3 4 5
i j

*i ← *j

4 1 0 3 4 5
i j

i ← j = 4, j ← from (4) = 2 which is start, so loop ends

4 1 0 3 4 5
j i

*i ← tmp

4 1 0 3 2 5
j i

This completes the first call to rotate_cycle_from in the while loop of our rotate function.

Exercise 11.10. Continue to trace the preceding example until the rotate function finishes.

Notice that the signatures of this rotate function and the previous one differ by the type of the last argument. In the next section, we’ll write the wrapper that lets the fastest implementation for a given situation be automatically invoked.