Permutations as Sets of Cycle Graphs (and why we care)

This post describes how to transform difficult problems involving permutations into easier problems involving graphs, then lists some problems that require this transformation and gives their solutions. In the process it demonstrates general techniques that may be useful in the face of such problems.

Introduction to permutations

A permutation of a set is a particular ordering of all the elements of the set. If we take the letters a, b, c, d, e, then c, d, b, e, a is a valid permutation, as is a, b, c, d, e or e, d, c, a, b or any other arrangement of these letters. a, c, e is not a permutation as it doesn't use all the elements, and a, a, b, c, d, e is also not a permutation because it uses an element more than once.

A permutation is a bijection from a set to itself, and permuting a sequence corresponds to applying the bijection element-wise to the sequence. If the set is $S = \{a, b, c, d, e\}$ and the permutation $f ~\colon~ S \to S$ is defined by $f(a) = b, f(b) = a, f(c) = e, f(d) = c, f(e) = d$, then applying this permutation to $a, b, c, d, e$ gives the sequence $f(a), f(b), \dots, f(e) = b, a, e, c, d$. The image of $f$ is also its domain so if $x \in S$, then $f(f(x)) \in S$, $f(f(f(x))) \in S$, and so on. In our permutation, $f(a) = b$ and $f(b) = a$, so $f(f(a)) = a$ and $f(f(b)) = b$. We say that $a$ is in a cycle of length 2 because it takes two applications of the permutation to $a$ to get back to $a$. $b$ is also in a cycle of length 2, and in fact we say that $a$ and $b$ are in the same cycle because to get back to $a$ by repeated applications of the permutation, we encounter $b$ and vice versa. In the same vein $c$ is in a cycle of length 3 along with $d$ and $e$. Every element of the set appears in exactly one cycle (note that it is possible to have a cycle of length 1) and cycles don't interact with each other. We say that any permutation can be written as a product of disjoint cycles.

Permutations can be represented as the union of a set of cycle graphs

Let's leave the world of algebra aside for a second and talk about graph theory. A directed graph is a set of vertices $V$ and edges $E$ where each edge maps from one vertex to another . A functional graph is a directed graph in which there is exactly one outgoing edge from each vertex. A permutation has a natural representation as a functional graph: let $V$ be the set on which the permutation is defined, and let there be an edge from $v_1$ to $v_2$ iff $f(v_1) = v_2$. The resulting graph is a specific class of functional graph where each vertex has exactly one incoming edge; in fact, the graph is the union of a set of cycle graphs, with each cycle graph corresponding to a disjoint cycle in the permutation.

Why do we care? This transformation allows us to transform problems about permutations into problems about graphs, and we know how to solve problems about graphs. In the rest of the post, I'll describe three different problems of increasing difficulty and show how to use this transformation to solve them. Without this transformation, these problems would be extremely difficult or impossible. The first problem describes a functional graph directly and the following two problems describe permutations.

Problem 1: a game involving a tennis ball

You are playing a game with $n-1$ other people and one tennis ball. You all sit in a circle and you hold the tennis ball in your hand. On the count of three, two things happen simultaneously:

You yell out a number $k$ between 2 and $n$ inclusive.

Everyone in the circle points at exactly one other person in the circle. You can assume that all participants except for you choose someone in a uniformly random fashion. (By a simple symmetry argument, this means that it doesn't matter who you point to, so we can assume you also point uniformly randomly)

You then throw the ball to the person you're pointing to, and they pass it to the person they're pointing to, and so on, so that the ball is passed a total of $k$ times. Whoever holds the ball at the end of $k$ passes loses the game.

What value of $k$ should you yell to minimise your chance of losing?

The underlying graph structure is immediately evident from the problem description: each player is a vertex and they point to each other to form a functional graph. The graph doesn't necessarily correspond to a permutation because it's possible for a vertex to have no incoming edges.

The solution follows from this line of reasoning: suppose you yell "2" and you lose. This means the ball was passed to whoever you pointed to ("Jane") and Jane returned it to you: you are in a cycle of size 2. Indeed, if you yell "2" you lose iff are in a cycle of size 2. Since each pointing arrangement (i.e. "valid graph" i.e. $n$-vertex functional graph with no cycles of length 1) is equally likely (proof left as exercise to the reader, but it's intuitive that if each person points uniformly randomly at someone else then all valid graphs are equally likely), we conclude that the probability of you losing is exactly the probability of randomly pulling a graph where you are in a cycle of size 2 from the uniform probability distribution over all valid graphs. This probability is equal to the number of valid graphs where you are in a cycle of size 2 divided by the total number of valid graphs.

Suppose you yelled "4" instead. If your graph happens to include you in a cycle of size 4, you will lose; but if you are in a cycle of size 2, you will also lose! And the number of valid graphs where you are in a cycle of size 2 OR a cycle of size 4 is obviously bigger than the number of valid graphs where you are only in a cycle of size 2. It follows that yelling "4" is always a worse choice than yelling 2. The same argument can be made for yelling 6, 8, and so on. In fact, the same argument can be made for yelling any composite number; it's always better to yell one of its prime factors instead.

In fact, this restricts us to yelling only primes. We just have to choose which prime. Some back-of-the-envelope combinatorial calculations show us the probability of being in a cycle of size 2 is $1/(n-1)$ which is larger than the probability of being in a cycle of size 3 which is larger than the probability of being in a cycle of length 5 and so on. The probability of being a member of a large cycle is smaller than the probability of being in a small cycle. We conclude that the solution to the problem is to choose $k$ to be the largest prime no greater than $n$.

The most interesting observation from this solution actually comes from the last few calculations that I skimmed: it turns out that large cycles in random functional graphs or large disjoint cycles in random permutations are extremely unlikely. In fact, the chance that you are in a cycle of length $n$ is $(n-1)!/(n-1)^{n-1}$. When $n=100$, this value in the order of $10^{-44}$.

Problem 2: 100 Prisoners Mk2

This is a variation on the classic 100 Prisoners problem.

The director of a prison offers 100 death row prisoners, who are numbered from 1 to 100, a last chance. A room contains a cupboard with 100 boxes in a row. The director randomly puts one prisoner's number in each closed box. The prisoners enter the room, one after another. Each prisoner may open and look into 50 boxes in any order. The boxes are closed again afterwards. If, during this search, every prisoner finds his number in one of the boxes, all prisoners are pardoned. If just one prisoner does not find his number, all prisoners die. Before the first prisoner enters the room, the prisoners may discuss strategy — but may not communicate once the first prisoner enters to look in the drawers.

The prisoners have also bribed a guard such that, before the prisoners commence the exercise, the guard goes into the room and is allowed to look into all the boxs and then swap the numbers in two boxes if the prisoners' strategy calls for it. The guard has no further communication with the prisoners after entering the room.

There is a strategy allowing all prisoners to survive. Find it.

The solution is to observe that the boxes can be numbered themselves from 1 to 100, in which case box $a$ containing the number of prisoner $b$ represents a permutation $f$ for which $f(a) = b$. Now suppose a prisoner $a$ enters the room and opens box $a$. $a$ points to $b$, which in turn points to $c$, and so on; since he is effectively re-applying the permutation (or walking the cycle graph, however you want to interpret it), he will eventually return to box $a$. But returning to box $a$ means opening a drawer containing prisoner number $a$ inside, which is the terminating condition.

The number of boxes the prisoner opens is exactly the size of the cycle in the graph containing his box. But what if this cycle is larger than size 50? Well, there can only be one such cycle in a graph of 100 vertices, and this is where the bribed guard plays his role! The guard first enters the room and determines if the graph corresponding to the arrangement of numbers contains a cycle of size 51 or more. If so, he simply breaks it into two smaller cycles by swapping two numbers in boxes. The resultant graph has no cycles larger than 50 vertices. Now all prisoners can follow this same strategy and all will survive. Brilliant!

The observation to take away from this problem is the same underlying graph/permutation can be interacted with in different ways (in this case it's walked along in different ways) but we take advantage of the fact that the graph/permutation itself is constant.

Problem 3: repeated number in array

Although the description for this problem is simple, it's an extremely difficult problem with a beautiful solution. It supposedly took Don Knuth 24 hours to solve; I discussed it with my team of Google engineers for over an hour and we concluded that it didn't have a solution. We were wrong.

Unlike the last problem, this problem benefits immensely from reducing it to a graph problem. However, now that you've seen the general strategy to solve problems of this type, I recommend you to attempt to solve it yourself before reading the solution.

Given an array of $n$ elements ranging from $1$ to $n-1$ in which exactly one element appears more than once, find the element in linear time and constant space without modifying the array.

Examples of arrays: [1, 2, 3, 4, 5, 5] or [4, 4, 4, 4, 1] or [6, 5, 2, 3, 5, 4, 1].

There's a trivial linear time and space solution using a flag array, and a fairly simple $O(n \log n)$ constant space solution using a modified binary search, but neither fulfils the requirements of the problem.

To solve it, notice that the first $n-1$ elements of the array are a permutation of the integers $1\dots n-1$, but where zero or the elements have been 'overwritten' with the integer $x$ in this range. If we try to construct the graph representing this almost-permutation, we get a functional graph where the vertex $V_x$ representing $x$ potentially has several incoming edges, and other vertices (the vertices representing overwritten array elements) now have no incoming edges. All components but one in the graph are cycle graphs, and the remaining component is rather messy: it comprises a cycle graph containing $V_x$ with a bunch of path graphs that enter the cycle at $V_x$. Note that the answer to the problem is $x$.

So far, we have ignored the final array element $k$ (represented by vertex $V_n$). $V_n$ has no incoming edges (since an incoming edge would mean there is an element $n$ somewhere in the array, which isn't possible because elements only range up to $n-1$) and one outgoing edge to vertex $V_k$ in the messy component of the functional graph described above. If we follow the edge from $V_n$ we go to $V_k$, and $V_k$ is either on a path graph to the messy component's cycle or it's on the cycle itself. So $V_n$ is the start of a path graph that ends at the messy component's cycle at vertex $V_x$.

Ignoring the vertices not in this path graph or the cycle it leads to, we end up with a rho-shaped graph, also sometimes (misleadingly) called a "linked list with a cycle", and if we find the length of the cycle, we can derive the length of the path graph preceding the cycle and hence find the vertex $V_x$ giving us the value $x$.

But we have efficient algorithms for finding cycles and their lengths! In particular, the hare and tortoise algorithm solves it in linear time and constant space.

In summary, the solution to the problem is to run a linear-time constant-space cycle detection algorithm on the implicit functional graph formed by the array beginning at its last element, and output the index in the array that corresponds to the vertex which starts the cycle.

The take-away here is that even problems that involve not-quite-permutations can be solved by relating the properties of the not-quite-permutation to the properties of permutations.

Conclusion

Problems that require interpreting an array as a permutation, or a permutation as a graph, are fairly rare. However, they often don't have alternative solutions, and if you haven't seen these tricks and techniques before, they can prove to be unsolvable.

Furthest 'Close' Pair in a Set

Another great algorithms question I heard recently:

Given a set $S$ of real numbers, call a pair of numbers $x,y \in S$ close if $x$ and $y$ appear consecutively in a sorted representation of $S$ (i.e. $x,y$ are close iff $S \cap (x,y) = \emptyset$). Find a close pair $x,y \in S$ that maximises $|x-y|$.

The trivial $O(n \log n)$ solution is to sort $S$ using any $O(n \log n)$ comparison sort, iterate through each pair of elements in $S$, and pick the pair whose absolute difference is greatest. However, there is a simple and elegant $O(n)$ solution.

If $S$ contained integers bounded in a reasonably-sized interval, we could modify the trivial solution described above to use a $O(n)$ sort such as radix sort. Unfortunately the elements of $S$ are arbitrary real numbers so no such property holds. However, we will try to use a bucket sort approach anyway. The first step is to decide on a bucketing scheme, and this is where the first "aha!" observation in the solution comes into play:

Claim: a solution $x,y$ to this problem satisfies $|x-y| \geq \frac{\max(S) - \min(S)}{|S|-1}$.

Where does that right-hand-side expression come from? To me it represents a sort of discrete mean value theorem. There are $|S|-1$ close pairs in $S$: the number of consecutive pair of elements in the sorted representation of $S$. Each of these pairs represents a disjoint interval in $\left[\min(S), \max(S)\right]$, and together they partition this entire interval. Say we divide the entire contiguous interval covered by $S$ into $|S|-1$ disjoint uniform partitions of length $L$. In any other partitioning of this interval (including the partitioning denoted by the close pairs in $S$), at least one partition must have length at least $L$ (similarly, at least one partition must have length at most $L$). This is intuitive and easy to visualise, but for the sticklers there's an easy proof by contradiction: suppose that the close pair $x,y \in S$ maximising $|x-y|$ has $|x-y| < \frac{\max(S) - \min(S)}{|S|-1}$. Then the size of the union of all the $|S|-1$ close pair intervals in $S$ is less than $\max(S) - \min(S)$. But the close pair intervals partition $\left[\min(S), \max(S)\right]$ so the size of the union of all close pair intervals must equal $\left|\left[\min(S), \max(S)\right]\right| = \max(S) - \min(S)$.

We will proceed with our bucketing using $|S|-1$ buckets each of size $\frac{\max(S) - \min(S)}{|S|-1} = M$. The first bucket contains elements $s \in S$ where $s < \min(S) + M$, the second bucket contains $s \in S$ where $\min(S) + M \leq s < \min(S) + 2M$, etc. The bucketing process is a linear time operation because for each $s \in S$, a simple mathematical expression tells us which bucket will contain $s$, and we can put it in a bucket in constant time. (Buckets would probably be implemented as linked lists or arrays).

Now what? Well, we've asserted that the elements in the solution pair have absolute distance at least $M$, so they can't both be in the same bucket. We only have to consider pairs of elements in different buckets. Also, if a bucket contains more than 2 elements, we can ignore all but the minimum and maximum in that bucket: no element $x$ in the bucket can form a close pair with an element $y$ of another bucket, because either the minimum or the maximum in the bucket will lie in between $x$ and $y$. So we will proceed as follows: for each bucket, find the minimum in that bucket and the maximum in that bucket, and delete all other elements of the bucket. This is a linear-time operation as every element in $S$ is processed at most twice (once to check if it's a maximum or minimum in its bucket, and potentially once more to erase it from its bucket).

Now we have $|S|-1 = O(n)$ buckets, each with at most two elements. We can simply iterate through each bucket in increasing order and build a list of the elements in all buckets, by adding first the minimum then the maximum in each bucket. This list has $O(n)$ elements, more importantly it's sorted, and most importantly, both the elements in the solution pair are guaranteed to be in the list. To finish off, we simply iterate through each consecutive pair of the list and pick the pair with the greatest difference.

In my opinion this is too much of an "aha!" problem to give a useful signal in an interview (either you see the $M$ trick or you don't), but it's fascinating nonetheless.

Edit: a colleague of mine found a similar solution that follows easier reasoning. Use $|S|+1$ buckets of equal size across the range $\left[\min(S), \max(S)\right]$. Since there are only $|S|$ elements, by [an alternative formulation of] the pigeonhole principle at least one of the buckets must be empty, so the difference between the closest pair is at least one bucket-width. This allows us to ignore all elements that aren't the min or max in a bucket. We now proceed as before.

Don't bound random integers with %

This will be a brief remark that will be the first in a series of posts about random number generation. Rather than discussing generation processes themselves (I'm no where near qualified enough to do that), I'll cover some notions of randomness, tests for randomness, and the role of those tests in software engineering.

You want to generate a random integer between 0 and $n$ inclusive. The usual advice is to do this:

num = random() % (n+1)

Where random() is some standard library function that outputs uniformly distributed random integers. This works. For most purposes, it works well. For cases where a uniform distribution is required, it's dangerous.

The problem is that random() outputs an integer less than some specific maximum. This maximum varies based on your programming language, so we'll refer to it as RAND_MAX – the name of the macro in C. On $k$-bit machines, RAND_MAX is often $2^{k-1} - 1$.

Assume for a second that RAND_MAX = 7, and you want a random number between 0 and 2 inclusive. You do num = random() % 3. Which outcomes of random() give num = 0? 0, 3, 6. What about num = 1? 1, 4, 7. What about num = 2? 2 and 5. Do you see the problem? Recall that the outputs of random() are uniformly randomly distributed, so the probability of num = 0 is $\frac37$, the probability that num = 1 is $\frac37$, but the probability that num = 2 is only $\frac27$. num is thus not uniformly distributed.

Obviously this effect is far less noticeable for larger RAND_MAX values, but it's present nonetheless. In cases where precision is a high priority, this bias is concerning.

So what's the right way? We can do:

do { num = random() } while (num > n)

but the expected number of iterations of this loop grows in inverse proportion to the size of $n$. Better to create a value $n'$ that's the lowest power of 2 not less than $n$, and instead do:

do { num = random() % n' } while (num > n)

This solution takes advantage of the fact that RAND_MAX is one less than a power of 2, so num % n' is uniformly distributed. For instance, in our previous example $n' = 4$, so num = random() % n' gives num = 0 when random() is 0 or 4, num = 1 when random() is 1 or 5, num = 2 when random() is 2 or 6, and num = 3 when random() is 3 or 7. This means that each possible outcome of num occurs with equal probability! However, if num = 3, we have to iterate again. Fortunately $n' < 2n$ so the expected number of iterations is less than 2 – a fair compromise for true uniform randomness.

A final word: I assume that your implementation of random() yields uniformly random integers. Standard library implementations of random() often do not satisfy this. To ensure uniform randomness, use a good random number generator.

In a future post, I'll discuss how to test the conformity of a black-box random number generator to the probability distribution it claims to choose from.

Counting Integers with a Particular Binary Property

The Hamming weight or popcount of $n$ is the number of 1s in the binary representation of $n$.

Question: how many positive integers less than $N$ have a popcount divisible by 3?

The naive $O(N \log N)$ solution is to loop through each integer $n < N$ and count the number of bits that are set. This is slow.

Don't fret! There is an $O(\log^2 N)$ solution.

Suppose that $N$ is a power of 2. Then its binary representation is 1 followed by $\lg(N)$ 0s. We can count the number of positive $n\leq N$ with $\mathrm{popcount}(n)$ by simply doing $C(\lg(N), 3) + C(\lg(N), 6) + ... + C(\lg(N), 3k)$ while $3k \leq \lg(N)$. For instance, when $N=16$, $C(4,3) = 4$ which is correct.

Next, suppose that $\mathrm{popcount}(N) = 2$, (i.e. $N$ has exactly two bits set, i.e. N$ $is the sum of a pair of distinct powers of 2). E.g. say $N = 20$ (binary 10100). The solution that worked for $N=16$ will still work. Now in a sense we've "considered" that power of 2. Then consider the next power of 2 that makes up the number, 4. 4 has 2 0's following, so we use the same principle as before! However, we remember that the bit representing 16 has already been set, so we're no longer looking for popcounts that are multiples of 3; instead, we're looking for popcounts which give a remainder of 2 modulo 3. Mathematically, we do $C(\lg(4), 2) + ... + C(\lg(4), 3k-1)$ while $3k-1 \leq 4$. In this case, we have only $C(2,2) = 1$. This corresponds to $n = 19$ (binary 10011). So for $N=20$, the answer is $4+1=5$.

OK now we generalise the pattern.

from euler import popcount, binom
from numpy import log2

def solve(N):
    ans = 0
    num_considered = 0
    p = 2**60
    oldN = N
    while p >= 1:
        if N & p:
            if (num_considered + 1) % 3 == 0:
                ans += 1
            mult = 3 - (num_considered % 3)
            while mult <= log2(p):
                ans += binom(log2(p), mult)
                mult += 3
            num_considered += 1
        p //= 2
    return ans

I'm not even sure that this is $O(\log^2 N)$... there are two nested loops that each run $O(\log n)$ times, and inside the nested loop we call a binomial function which in theory is $O(nk)$, and we call it with $n,k \leq \log N$ producing an extra $O(\log^2 N)$ multiplier. But due to the ranges of the arguments we use for the binomial function, I'm pretty sure each call is amortized constant time. Maybe.

euler is my Project Euler helper functions module which is available on my Github.

P vs NP for Dummies

I seem to type up a variation of this post every few months. I'm saving this one for the permanent record.

P vs NP for Dummies

There are some yes/no problems which are easy. If I give you a list on numbers and ask you to tell me if the number 5 is in that list, you can solve it easily by looking at each number in the list. In the worst case, when the number 5 is at the very end of the list, this takes about n steps/instructions.

There are yes/no problems which are slightly harder -- for example, I might ask you to tell me if 12*14=168. When you multiply two n-digit numbers using the method you learnt in school, you need about n² steps.

Both of these problems require some number of steps that's bounded by a polynomial of the input size; i.e. for each of these problems there is a fixed k so that the number of steps required to solve an instance of that problem of size n it is always bounded by n^k . This class of problems is called P which is a terrible pseudo-abbreviation for 'polynomial time on a deterministic Turing machine' (a deterministic Turing machine is a standard computer, more or less).

Here's another yes/no problem: given a map of a country with n cities, determine whether there's a route that can be taken that visits every city exactly once, and whose total distance is less that a given distance d. If we are told that the answer is yes and we are given the corresponding route, we can easily check its total distance by simply following the route and adding up all the distances between cities. This verification process takes about n steps since you're adding n-1 numbers together. The set of yes/no problems whose solutions can be verified in polynomial time (given some extra information, in this case the route itself) is called NP. This specific problem is a special case called an NP-complete problem, meaning that it is in a sense 'harder' than all other problems in NP.

One might conjecture that, because you can verify a solution to this problem in a polynomial number of steps, then perhaps you can also solve it from scratch in polynomial time? Obviously all problems in P are also in NP, since if we can verify a solution to any problem in P in polynomial time by simply solving that problem and checking it against the given solution. Since this problem is 'harder' than all the other NP problems, if we can solve this one in polynomial time, then we can solve all the other NP problems in polynomial time too! This would mean that all NP problems are also P problems, in other words P = NP.

So you set out to try to find a polynomial time method to produce a solution to this problem. But you can't. No one can. We don't have a polynomial time method to solve for this problem. In fact we don't have a polynomial time method for any of the dozens of known NP-complete problems. These problems range from finding routes across maps, to packing items efficiently into a backpack, to optimally playing candy crush. On the flip side though, no one's managed to prove that a polynomial time method doesn't exist. Such a proof would demonstrate that P =/= NP.

So there you are. The problem of P vs NP is the problem of proving either that P=NP (this can be done simply by finding a polynomial time method to solve problem above) or that P =/= NP (this intuitively seems harder to prove).

For the route-finding problem described above (it's called the travelling salesman problem), the best currently known method is a small optimisation on "brute force every permutation of cities to visit". This method bound the number of steps by about n² * 2^n. If you can get rid of the exponential term in that expression, you get a million dollar prize and immediately render most of the field of cryptography obsolete.

An Original String Manipulation Problem

This is an original string manipulation problem that I wrote. It was included on the 2014 Australian Invitational Informatics Olympiad. In my opinion it would be a question of moderate difficulty on an IOI exam.

You are given a critical string of length $C$ in some arbitrarily large alphabet $\Sigma$ and some text of length $T$ in $\Sigma$ ($C \leq T \leq 10^6$). You are guaranteed that each character in $\Sigma$ appears at most once in the critical string. You are then instructed to perform $U$ update operations; each operation is of the form "set the character at position $i$ in the text to be $x$" ($0 \leq i < T, x \in \Sigma, U \leq 10^6$) . After each update operation, output the maximum number of consecutive occurrences ("runs") of the critical string in the text. The intended time complexity is linearithmic.

The difficulty in this question arises from the number of different sub-algorithms/techniques (three or four of them) that need to be applied to construct a complete solution. What follows is my solution.
Before we can track runs of the critical word within the text, we first need the ability to detect the formation of individual instances of the critical word. Intuitively, some update operations make us closer to a critical word; e.g. if the critical word is "badger", then an update operation "booger" $\to$ "bodger" gets us closer to "badger". Inversely, "bodger" $\to$ "booger" gets us further from the formation of the critical word. Specifically, some $C$-length substring of the text gets closer/further if its edit distance to the critical string (the number of characters that differ across both strings) decreases/increases respectively.
The following observation drives our first technique.

---

Observation 1. Each update operation results in a decrease in edit distance to at most one instance of a $C$-length substring of the text, and an increase in edit distance to most one instance of a $C$-length substring in the text.
This is a result of the character-distinctness property of the critical string. Say we set a character in the text to 'a'. If the critical string does not contain 'a', then no text substring's edit distance could have decreased. If the critical string contains one occurrence of 'a' at position $p_a$, then only the text substring whose $p_a$th character has become $a$ has a reduction in edit distance. The critical string cannot contain more than one occurrence of 'a' so we need not consider more cases.
Similarly, if we are replacing the character 'b' with 'a', then an increase an edit distance occurs only for the substring whose $p_b$th character is replaced, or for no strings if 'b' does not appear in the critical string.

---

Sub-algorithm 1. By maintaining a map of character to critical string position (i.e. $c \to p_c$) we can find which $C$-length substring gets closer (if any) and which $C$-length substring gets further (if any) in constant time.
If we are updating position $i$ in the text from character $c$ to character $d$, the substring that gets closer begins at position $i - p_c$ in the text and the substring that gets further begins at position $i - p_d$ (all positions are 0-based by the way).

---

Sub-algorithm 2. Store, for each position $i \in 0..T-C+1$ in the text, the edit distance of the $C$-length substring beginning at position $i$ in the text. After each update operation, make the zero, one, or two necessary modifications to these values (as indicated by sub-algorithm 1). After replacing the character $c$ at position $i$ with $d$, we can determine in constant time if a critical word instance has been formed or un-formed by querying the edit distance at $i-p_c$ and $i-p_d$. This is all constant time.

---

Sub-algorithm 2 allow us to detect formations and un-formations of the critical word in the text as we go, in $O(1)$ time for each update. The next step is to detect runs of the critical word in the text. To do this, we introduce a complex data structure into the equation.
I'll begin by posing a more straight-forward problem. We have a binary string. We then perform a series of updates on the string. In each update we set or unset a bit in the string. After each update, we need to output the length of the longest run of 1's in the string.
In fact, we can implement a range tree that handles these exact update and query operations. I won't go into the implementation details, but I will mention what each tree node stores:

• The longest run in the interval contained by the node
• The size of the run that is incident to the left side of the node's interval
• The size of the run that is incident to the right side of the node's interval

We'll call this the longest-run tree. Using this data structure (in fact, $C$ different instances of it!), we get our third technique:

---

Sub-algorithm 3. Instantiate $C$ longest-run trees, where each represents a different 'offset' of a critical string run within the text. The first 'offset' has the first letter of the critical string at indices 0, $C, 2C, ...$ within the text, the second 'offset' has the first letter at indices $1, C+1, 2C+1, ...$, etc.
Each time we detect the completion of a critical word within the text (using sub-algorithm 2), let the position of the beginning of the critical word be $p_0$. Take the $(p_0 \bmod C)$'th longest-run tree and "set the bit" in the tree at index $\left\lfloor \frac{p_0}{C} \right\rfloor$. Do the same for any critical words that we unformed -- find their offset's tree, find the position in the tree, and unset the bit.
Now if we query that tree, we discover the longest run of consecutive critical words in that 'offset'. Setting and unsetting bits is $O(\lg T)$ and querying is also $O(\lg T)$.

---

We are almost done! After making an update, we can query the longest run for a particular offset of the critical word in logarithmic time. But how do we know if this offset contains the maximum-length run?
For this, we introduce yet another data structure (an easier one this time). We use a max-first priority queue.

---

Sub-algorithm 4. Initialise a priority queue of integer pairs where each pair is of the form (maximum run-length in tree, tree ID), sorted in descending order by the former. Each time we "set a bit" in a tree, query the maximum run-length of that tree. Push onto the priority queue the maximum run-length in that tree as well as the tree ID.
Now when it comes time to output the maximum run-length of the critical string, query the priority queue. This will give you the maximum run-length across all trees, as well as the tree in which that maximum run-length occurs. However, we can't necessarily trust this value! What if we'd pushed it to the queue a long time ago, and subsequent updates have removed that run of the critical word? To cater for this, we must now query the tree with the tree ID that we just popped to ensure its run length is truly the longest. If it's not, we ignore this priority queue element, pop it off and try the next one.
You can also use a set as your priority queue, and erase elements when you unset bits rather than the lazy deletion we're using.

---

Each time we do an update, we push something to the queue. We might set a bit up to $T$ times, so potentially we could have $T$ elements on the queue (though this is really unlikely). A queue with $T$ elements has $O(\lg T)$-time operations. Each time we do an update, we potentially pop several things off the queue, but since the total number of elements on the queue is bounded, this amortizes to one pop per update. So the time complexity of each update is the time complexity of a constant number of range tree operations followed by an amortized constant number of priority queue operations. $O(\lg T + \lg U) = O(\log TU)$.
Since we do this for each update, we end up with $O(U \log TU)$ as required.

---

The implementation is not as horrible as it seems! Mine was 110 lines of C++.

#include <cstdio>
#include <cstring>
#include <map>
#include <algorithm>
#include <queue>
#include <cassert>
using namespace std;
#define MAX_N 1000005
#define MAX_W MAX_N
#define ROOT 0
#define LC(n) (2*(n)+1)
#define RC(n) (2*(n)+2)

struct node {
    int left, right, mid;
};

node *trees[MAX_N];
map<int, int> subtract;
int counters[MAX_N];
int N, W, Q, word[MAX_W], searchStr[MAX_N];
priority_queue<pair<int, int> > maxConsPQ; // <max cons in tree, tree id>
int LAST;
int total;

void put(node t[], int cur, int u, int ns, int ne, bool v) {
    assert(ns <= u && u <= ne);
    if (ns == u && u == ne) {
        t[cur].left = t[cur].right = t[cur].mid = v;
    } else {
        if (u <= (ns+ne) / 2) {
            put(t, LC(cur), u, ns, (ns+ne)/2, v);
        } else {
            put(t, RC(cur), u, (ns + ne)/2 + 1, ne, v);
        }
        int half = (ne - ns + 1) / 2, mid = 0;
        mid = max(mid, t[LC(cur)].right + t[RC(cur)].left);
        mid = max(mid, max(t[LC(cur)].mid, t[RC(cur)].mid));
        mid = max(mid, max(t[LC(cur)].left, t[RC(cur)].right));
        t[cur].mid = mid;
        t[cur].left = t[LC(cur)].left;
        if (t[LC(cur)].left == half) {
            t[cur].left += t[RC(cur)].left;
        }
        t[cur].right = t[RC(cur)].right;
        if (t[RC(cur)].right == half) {
            t[cur].right += t[LC(cur)].right;
        }
    }
}

void change(int pos, int nc, bool isInitial) {
    int oc = searchStr[pos];
    if (!isInitial && subtract.find(oc) != subtract.end() && pos - subtract[oc] >= 0) {
        int po = pos - subtract[oc];
        if (W == counters[po]) {
            total--;
            put(trees[po % W], ROOT, po / W, 0, LAST, false);
            maxConsPQ.push(make_pair(trees[po%W][ROOT].mid, po % W));
        }
        counters[po]--;
        assert(counters[po] >= 0);
    }
    searchStr[pos] = nc;
    if (subtract.find(nc) != subtract.end() && pos - subtract[nc] >= 0) {
        int pn = pos - subtract[nc];
        counters[pn] ++;
        assert(counters[pn] <= W);
        if (W == counters[pn]) {
            total++;
            put(trees[pn % W], ROOT, pn / W, 0, LAST, true);
            maxConsPQ.push(make_pair(trees[pn%W][ROOT].mid, pn % W));
        }
    }
}

int query() {
    while (!maxConsPQ.empty() && trees[maxConsPQ.top().second][ROOT].mid != maxConsPQ.top().first) {
        maxConsPQ.pop();
    }
    return maxConsPQ.empty() ? 0 : maxConsPQ.top().first;
}

int main() {
    scanf("%d %d", &W, &N);
    int treeSize = 1;
    while (treeSize < N/W) treeSize <<= 1;
    LAST = treeSize - 1;

    for (int w = 0; w < W; w++) {
        scanf("%d", &word[w]);
        trees[w] = (node*) malloc(sizeof(node) * 2*treeSize);
        memset(trees[w], 0, sizeof(node) * 2*treeSize);
        subtract[word[w]] = w;
    }
    for (int n = 0; n < N; n++) {
        scanf("%d", &searchStr[n]);
        change(n, searchStr[n], true);
    }

    scanf("%d", &Q);
    for (int q = 0; q < Q; q++) {
        int i, x;
        scanf("%d %d", &i, &x);
        change(i-1, x, false);
        printf("%d %d\n", total, query());
    }

    return 0;
}

NFA to DFA

I've started reading a theory of computation textbook in my spare time. In the first few chapters, I re-learnt the algorithm for converting a non-deterministic finite automaton to a deterministic finite automaton that recognises the same language. I say "relearnt" because I actually learnt the same algorithm 2 years ago at the National Computer Science Summer School.

2 years ago, implementing the algorithm was one of the most frustrating (to debug and to reason about) programming tasks I'd ever done. If you're not familiar with the algorithm, it essentially involves breadth-first searching over a graph -- except it's more like a metagraph because each node represents a subset of the states (i.e. nodes) in the NFA (which itself is a graph). So in this metagraph, $A \xrightarrow{c} B$ iff the $B$'s state subset is the union of all the states reachable by following a $c$-transition from any of the states in $A$'s state subset. There's some extra stuff to handle $\epsilon$-transitions too.

Debugging this algorithm was as hard as it sounds. Here's the code from that troublesome day.

EPS = '~'

class NFANode(object):
    def __init__(self, is_final):
        self.is_final = is_final
        self.edges = set()

    def add(self, toks, next):
        for t in toks:
            self.edges.add((t, next))

    def EC(self):
        nodes = {self}
        for tok, next in self.edges:
            if tok == EPS:
                nodes.update(next.EC())
        return nodes


class DFANode(object):
    def __init__(self):
        self.edges = []
        self.states = set()
        
    def is_final(self):
        return any(k.is_final for k in self.states)
    
    def EC(self):
        nodes = set()
        for nfaNode in self.states:
            nodes.add(nfaNode)
            for tok, next in nfaNode.edges:
                if tok == EPS:
                    nodes.update(next.EC())
        dfaNode = DFANode()
        dfaNode.states = nodes
        return dfaNode

    def M(self, queryTok):
        nodes = set()
        for nfaNode in self.states:
            for tok, next in nfaNode.edges:
                if tok == queryTok:
                    nodes.add(next)
        dfaNode = DFANode()
        dfaNode.states = nodes
        return dfaNode


class NFA(object):
    def __init__(self, start):
        self.start = start

    def get_toks(self):
        tok_set = set()
        self.accepts(EPS*1000, tok_set)
        return tok_set.difference({EPS})
    
    def accepts(self, query, all_toks=None):
        """
        >>> n1, n2, n3 = NFANode(True), NFANode(False), NFANode(False)
        >>> n1.add('b', n2)
        >>> n1.add(EPS, n3)
        >>> n2.add('a', n2)
        >>> n2.add('ab', n3) 
        >>> n3.add('a', n1)
        >>> nfa = NFA(n1)
        >>> assert nfa.accepts('aaa')
        >>> assert not nfa.accepts('bb')
        >>> assert nfa.accepts('abababababababababbababababa')
        >>> assert nfa.accepts('baa')
        >>> assert nfa.accepts('baaaaaa')
        >>> assert nfa.accepts('aa')
        >>> assert nfa.accepts('baba')
        >>> assert not nfa.accepts('baababaaaaaaaaaaaaaaaaaab')
        >>> assert nfa.accepts('baababaaaaaaaaaaaaaaaaaaba')
        """
        q = [(0, self.start)]
        seen = set()
        
        while q:            
            pos, node = q.pop(0)
            if (pos, node) in seen:
                continue
            seen.add((pos, node))

            if pos == len(query):
                if node.is_final:
                    return True
            else:
                for (tok, next) in node.edges:
                    if isinstance(all_toks, set):
                        all_toks.add(tok)
                    if tok == EPS or (tok == query[pos] and pos < len(query)):
                        q.append((pos+(tok!=EPS), next))
        return False


class DFA(NFA):
    def __init__(self, nfaToConvert):
        self.nfa = nfaToConvert
        self.all_toks = self.nfa.get_toks()
        self.start = DFANode()
        self.start.states.add(self.nfa.start)
        # mapping of states tuple to DFANode
        self.cache = {}
        self._convert()

    def _convert(self):
        self.start = self.start.EC()
        self.cache[tuple(self.start.states)] = self.start
        
        q = [self.start]
        seen = set()

        while q:
            states = q.pop(0)
            for tok in self.all_toks:
                newDFANode = self.retrieve(states.M(tok).EC())
                if newDFANode.states:
                    states.edges.append((tok, newDFANode))
                    if tuple(newDFANode.states) not in seen:
                        seen.add(tuple(newDFANode.states))
                        q.append(newDFANode)
    
    def retrieve(self, node):
        if tuple(node.states) not in self.cache:
            self.cache[tuple(node.states)] = node
        return self.cache[tuple(node.states)]

    def accepts(self, query):
        """
        >>> n1, n2, n3 = NFANode(True), NFANode(False), NFANode(False)
        >>> n1.add('b', n2)
        >>> n1.add(EPS, n3)
        >>> n2.add('a', n2)
        >>> n2.add('ab', n3) 
        >>> n3.add('a', n1)
        >>> nfa = NFA(n1)
        >>> dfa = DFA(nfa)
        >>> assert dfa.accepts('aaa')
        >>> assert not dfa.accepts('bb')
        >>> assert dfa.accepts('abababababababababbababababa')
        >>> assert dfa.accepts('baa')
        >>> assert dfa.accepts('baaaaaa')
        >>> assert dfa.accepts('aa')
        >>> assert dfa.accepts('baba')
        >>> assert not dfa.accepts('baababaaaaaaaaaaaaaaaaaab')
        >>> assert dfa.accepts('baababaaaaaaaaaaaaaaaaaaba')
        """
        q = [(0, self.start)]
        seen = set()

        while q:            
            pos, node = q.pop(0)
            if (pos, node) in seen:
                continue
            seen.add((pos, node))

            if pos == len(query):
                if node.is_final():
                    return True
            else:
                for (tok, next) in node.edges:
                    if tok == query[pos] and pos < len(query):
                        q.append((pos+(tok!=EPS), next))
        return False


if __name__ == '__main__':
    import doctest
    doctest.testmod()

Combinatorial Knapsack Family

This cute little set of algorithmic triplets solve like 60% of all dynamic programming problems.

(created by me. correctness not guaranteed)

Tutorial on MaxScore (Information Retrieval)

In the Information Retrieval & Web Search course I'm taking, we look at an algorithm for document-at-a-time scoring called MaxScore (also written as max_score, maxscore or max score in the literature). Now MaxScore is a pretty cool algorithm, but unfortunately you wouldn't know that, because to learn MaxScore you have to piece together the ideas behind the technique from half a dozen hand-waved explanations in various academic papers.

This is a great loss to the recreational programming community, so I've written (to my knowledge) the first thorough and in-depth explanation of MaxScore. It includes a general outline, examples, and implementation details, as well as a detailed explanation and proof of why it works. Hopefully it will be useful to others.

• Link here
• Alternative link

Proving a Bound on Quad-tree Queries

A while ago, I made the claim on this blog that "Quad-tree Range Queries are O(sqrt(S)) Time". Today, I thought I'd prove that this is true. Technically, I'll prove a much more formal statement, and then hand-wave a bit to give the above. We will prove that:

The tightest attainable upper-bound on the size of the minimum set of quad-tree nodes whose union is exactly the span of a given rectangle in an $N\times N$ space is $N$.

To prove this, we will consider the general statement "the tightest attainable upper-bound on the number of quad-tree nodes whose union spans a rectangle in an $N\times N$ space is $K$". Then we will show that $K \leq N$ and furthermore that $N \leq K$, allowing us to conclude that $N=K$.

First, we show that $K \leq N$. To do this, it is sufficient to show that there is at least one rectangle that requires a minimum set of nodes of size at least $N$. We can easily do this by choosing any rectangle of height 1 and width $N$. Since each quad-tree node spans a square, we will need each node's corresponding square to have a height of exactly 1 (so it doesn't exceed the height of our rectangle). But since they're squares, their width must also be 1 – so we'd need $N$ of these nodes to describe the rectangle. Thus our minimum set is of size $N$, so $K \leq N$.

Now we will show that $N \leq K$. Intuitively this seems more difficult because we have to prove that for all rectangles, we can find a minimum spanning set of nodes of size at most $N$. We will show this result by using the steps of the standard recursive algorithm for querying a rectangle.

Recall that the algorithm to operate on a quad-tree range looks like:

FUNCTION RangeOperate(node, query_rect)
  IF node.rect is completely contained within query_rect:
    DO OPERATION ON node // This node has been 'selected' by the query
  ELSE IF node.rect is completely disjoint from query_rect:
    RETURN // Nothing more to do
  ELSE:
    RangeOperate(node.quadrant_1, query_rect)
    RangeOperate(node.quadrant_2, query_rect)
    RangeOperate(node.quadrant_3, query_rect)
    RangeOperate(node.quadrant_4, query_rect)

We seek an upper bound on the number of nodes that will be operated upon. Each time RangeOperate recurses, it shifts its attention from an $n\times n$ node to an $\frac{n}{2} \times \frac{n}{2}$ node. Because the function can initially be called on a rectangle of size $N\times N$, and since it can't go any deeper than $1\times 1$ nodes, the function can recurse to a maximum depth of $\log_2 N$. Keep this in mind, we'll need it later.

We make at most 4 recursive sub-calls per call of RangeOperate. Let's list all the possible types of rectangles that query_rect could be, and for each one, consider what these sub-calls will do.

• If query_rect completely contains or is completely disjoint from node.rect, we don't recurse at all. Call this the termination condition.
• If query_rect lies within only one of the four quadrants of node, then all but one of the four sub-calls will encounter the termination condition. The other will continue recursing.
• If query_rect lies within two of the four quadrants of node, then two of the four sub-calls will encounter the termination condition. The other two will continue recursing.
• If query_rect lies within all four quadrants of node, none of the four sub-calls will terminate. However, the query_rect for each of these sub-calls will overlap the sub-call's node.rect at one of its corners. Let's use subnode for a sub-call's node rect; query_rect still refers to the same rectangle. Without loss of generality, let's assume that query_rect overlaps the top-right corner of subnode.rect (the three other cases are rotations of this case).
  • If query_rect only overlaps with quadrant 1 of subnode, then only one sub-call of this sub-call (the sub-sub-call if you will) will not meet the termination condition.
  • If query_rect overlaps with quadrants 1 and 2 of subnode or quadrants 1 and 3 of subnode, then only these two sub-sub-calls will not meet their termination condition.
  • If query_rect overlaps with all four quadrants of subnode, then quadrant 1 of subnode must be completely within the query_rect; so this sub-sub-call meets the termination condition. The sub-sub-calls on quadrants 2 and 3 have query_rect overlapping two corners of subsubnode.rect, and thinking about this for a moment will convince you that at least two of the sub-sub-sub-calls will then meet their termination condition.

By visualising the cases, you'll quickly convince yourself that it is impossible for query_rect to lie within exactly three of node's quadrants, so this list is exhaustive. What was the point of all this? Well, we established that after recursing down three levels from any call of RangeOperate, at most 8 sub-sub-sub-calls won't have met their termination conditions. We can amortize this to prove that at most 2 sub-calls will remain active from any given call of RangeOperate – and since our recursion tree has depth $\log_2 N$ and branching factor 2, the total number of nodes operated on is at most $2^{\log_2 N} = N$.

Of course, there are other interesting questions raised by this exercise:

• Do these results generalise to higher dimensions? Oct-trees and such?
• What can we say about quad-trees that span non-square spaces? (Hint: if the space spanned is $W\times H$, we can adapt what we've done already to show $K=\max(W, H)$.
• We've shown that the upper-bound on the minimum set size is $N$, and we have (implicitly) shown that the standard quad-tree querying algorithm never queries a set of nodes larger than $N$ by proving the $N \leq K$ result using the steps of that same algorithm. But what guarantee do we have that this quad-tree querying algorithm will find a minimum set for any rectangle? For instance, say we had two algorithms, the standard query algorithm $S$ and another query algorithm $P$. We have established that the number of nodes queried by each of these algorithms is at most $N$, but suppose that $S$ queries on average $N/2$ nodes and $Q$ queries on average $N/3$ nodes. Can we prove that $S$ does in fact select a minimum set for any rectangle query?

Simple In-Place Swap

Beginner programmers are often given the task of swapping the values of two variables. First, they will try to assign to each variable the value of the other.

x := y
y := x

Of course, this doesn't behave as expected -- the first assignment overwrites the original value of $x$. The standard solution to this issue is to use a temporary variable to duplicate the value of one of the input variables.

temp := x
x := y
y := temp

For the more mathematically-inclined student, there's an extension to this puzzle: swap the variables' values without storing intermediate values.

The trick is to notice that you can hold the original values of $x$ and $y$ implicitly as simple arithmetic results. For example, let $x := x + y$. We haven't lost the original value of $x$ -- we still have it 'stored' as $x - y$. Once this observation is discovered, some trivial arithmetic on a few small examples is sufficient to find the solution.

x := x + y
y := x - y
x := x - y

To verify that this algorithm actually works, we can let $y'$ be the final value of $y$ and $x'$ be the final value of $x$. By substituting previous values into the above equations, we see that $y' = (x + y) - y = x$ and $x' = (x + y) - ((x + y) - y) = y$, confirming that the values have indeed been swapped.

This is a particularly good task to give to children who are learning to program. It introduces to them the idea that programming isn't just about menial bit-twiddling, and reinforces the importance of analytical problem solving in computer science. It's also very approachable yet still conceptually interesting.

Capacity Scaling for Fast Max Flow

In my previous post, I presented an implementation of the Ford-Fulkerson method (termed 'method' and not 'algorithm' because it encompasses a number of different algorithms) for finding the maximum flow through a network.

The Ford-Fulkerson method involves repeatedly finding so-called augmenting paths in the residual network, then saturating these paths, terminating only when there are no augmenting paths left. The order in which these paths are found is dependent on the specific algorithm employed. The Edmonds-Karp algorithm, for instance, uses a breadth-first search to find the shortest augmenting path. Another algorithm employs a priority queue to find the maximum-flow augmenting path.

There is another approach, however, which accomplishes the same goal but attacks the problem from a different direction: the capacity scaling algorithm. Capacity-scaling for max flow has woefully little literature on it. I found mentions of it in CS lecture slides from several universities, and also in a TopCoder algorithms tutorial. Apparently it's presented in Network Flows: Theory, Algorithms, and Applications, but I'd neven read about it before. In practice, capacity scaling is just as simple to conceptualise, just as easy to implement, and quite a bit faster than traditional Ford-Fulkerson.

Ford-Fulkerson first finds an augmenting path, then saturates that path. Capacity scaling turns this principle on its head: start by assuming that an augmenting path exists with residual flow at least $r$, then find such a path and add $r$ to its flow. (this only works for integral edge weights by the way!)

In particular, we will begin by choosing a value of $r$ such that $r$ is the largest power of 2 that is less than the maximum edge weight in the graph. We will then push $r$ units of flow out of the source node, and hope that all $r$ units reach the sink. If they do, we try push another $r$ units out. Otherwise, set $r := r/2$ and repeat. Iterate while $r \geq 1$.

max_flow = 0
r = 2 ^ max_edge_weight
while r >= 1:
  do:
    push r units from source to sink
    max_flow += r
  while (successful)
  r /= 2
output max_flow

The algorithm should be fairly intuitive. It's clear that it finds the maximum flow, because it terminates only when it cannot push a single unit of flow i.e. there are no augmenting paths left. Furthermore, it pushes flow down an edge at most $\left \lfloor \lg U \right \rfloor$ times, where $U$ is the maximum edge weight in the graph: because each pass considers a different power of 2, one can interpret the algorithm as effectively 'setting each bit' of the flow amount for each edge.

Now to derive the time complexity. First, I want to make a claim that will help us out. I try to avoid proofs, but this one is important so I'll sketch it out.

EDIT: I'm no longer convinced of this proof's correctness. The problem is that by pushing negative flow, you can effectively undo work you did on a previous iteration. I'll leave it here anyway for readers to pick holes in.

Claim: For a graph $G(V,E)$ the algorithm retains a single value of $r$ for at most $|E|$ iterations of the flow-pushing routine.

Proof: Consider the initial case $r=2^{ \left \lfloor \lg U \right \rfloor}, U=\max \; \{ W(e) \mid e \in E \} $, where by design, no edge has residual flow more than twice $r$. Each time we successfully push flow, we decrease the residual flow of at least one edge $e$ by $2^{ \left \lfloor \lg U \right \rfloor}$. Now $W(e)-2^{ \left \lfloor \lg U \right \rfloor} \leq U - 2^{ \left \lfloor \lg U \right \rfloor}$ (since $U$ is chosen to be the greatest of all $W(e)$'s, and furthermore $U - 2^{ \left \lfloor \lg U \right \rfloor} \leq 2^{ \left \lfloor \lg U \right \rfloor - 1}$. So $e$ now has weight $2^{ \left \lfloor \lg U \right \rfloor - 1} < 2^{ \left \lfloor \lg U \right \rfloor}$ so we cannot push $r$ more units of flow down $e$. This process can only occur once for each edge, so the value of r is retained for at most $|E|$ iterations. Following this, we let $r:=r/2$ and classify edges as follows: those that we augmented on the previous iteration (which now have residual flow not more than twice $r$), those that we could not traverse on the previous iteration (because their values were not more than twice $r$), and those which were unreachable on the previous iteration (either because they're not in the connected component or because they can only be reached by first traversing an edge of residual flow not more than twice $r$ -- and hence can be thought of as having the residual flow of this earlier edge). Observe that this scenario is identical to our initial state, in which all reachable edges have $W(e) \leq 2r$, so by induction the algorithm retains a single value of $r$ until the base case $r < 1$.

Since there are $\left \lfloor \lg U \right \rfloor$ values of $r$ to trial, at most $|E|$ repetitions of the push-routine for each value of $r$, and each push-routine conducts an $O(E)$ search (BFS works well), we arrive at a final complexity of $O(E^2 \log U)$. This outperforms Edmonds-Karp's $O(VE^2)$ and the traditional iterated DFS both in theory and in practice.

Code will follow shortly.