Bounded Counting Knapsack

Counting knapsack problems look like this: given $n$ items with respective integer weights $w_1,\dots,w_n$, how many ways can we form combinations of these items that have total weight $s$?

There are three basic variations of counting knapsack.

1. In the “unbounded” variation we allow combinations to contain repeated items. For instance with 3 items $a,b,c$ with weights 1,3,5 we can form a combination of $3a, 2b, 1c$ with total weight 16. The partition function $p(s)$ in number theory is given by solving this problem for $n=s$ and $w_i=i$ for all $i$.
2. In the “binary” or “0/1” variation, we can only use each item once. The weights need not be distinct, so if the $w_i$ are 1,1,3,5 respectively then we can form a combination 1,1,5 that sums to 7, but we cannot form a combination 1,3,3.
3. In the "bounded" variation, we are provided with integers $c_1, \dots, c_n$ and we are allowed to use up to $c_i$ instances of item $i$ in each combination. This is a generalisation of the two variations above, where in the unbounded variation we have $c_i = \infty$ for all $i$ and in the 0/1 variation we have $c_i = 1$ for all $i$.

The unbounded and 0/1 variations are classic combinatorics examples that can be solved in $O(ns)$ time and $O(n)$ space using basic dynamic programming. Denote $[i,k]$ to be the number of combinations of items $1, \dots, i$ with total weight $k$, so that the solution to the problem is given by $[n,s]$. Then the unbounded variation has recurrence $$[i,k] = [i-1,k] + [i,k-w_i]$$ This recurrence has the interpretation "to form a combination of weight $k$ using items $1, \dots, i$ we can either add the $i$th item to a combination of weight $k-w_i$ formed by items $1, \dots, i$, or we can just take any existing combination of weight $k$ formed by items $1, \dots, i-1$. The 0/1 variation has the very similar recurrence $$[i,k] = [i-1,k] + [i-1,k-w_i]$$ with an analogous interpretation; the only difference is that we are prevented from using item $i$ more than once in a combination. The $O(n)$ space bound is given by observing that in both cases, we can evaluate $[i, k]$ by referring only to $[j,k]$ for $j \geq i-1$, so any $[j,k]$ with $j<i-1$ can be discarded. In all variations we have base cases $[i,k] = 1$ if $i=k=0$ and $[i,k] = 0$ otherwise.

For the bounded variation things get more interesting. We might expect a solution to look like the solution to the solution for the unbounded variation, except that the unbounded solution has no "memory" for how many times it has used an item, so it has no inbuilt way of preventing us from exceeding our allowance for an item. The 0/1 variation, on the other hand, generalises naturally to the bounded case. We just replace the binary decision in its recurrence – "do we use item $i$ or not?" – by a more general decision "how many instances of item $i$ do we use?". The resulting recurrence is

$$[i,k] = \sum_{q=0}^{c_i} [i-1, k-qw_i]$$

When $c_i = 1$, this is exactly equivalent to the 0/1 variation recurrence. However, when $c_i = \infty$, this looks quite different to the recurrence for the unbounded case. Why? It is instructive to think hard about the answer to this question.

What is the time complexity of this solution? It is not $O(ns)$ because clearly the recurrence is no longer $O(1)$. Overall, for each item $i$ we sum over $c_i$ values $s$ times, so we could say it is $O(s\sum_i c_i)$. In the worst case our items have small $w_i$ and $c_i \approx s$, so $\sum_i c_i$ starts to look like $ns$ and we're looking at $O(s^2n)$.

At this point we observe that all $w_i$ are distinct: if we had $w_i = w_j$ for $i \neq j$, we could replace items $i,j$ with a new item $i+j$ with weight $w_{i+j} = w_i + w_j$ and $c_{i+j} = c_i + c_j$. For any $i,k$ we have $[i,k-w_i] = 0$ whenever $k-w_i \leq 0$, so for any fixed $w_i$ we are summing over at most $k/w_i$ terms. Using the distinctness of $w_i$ we have

$$\sum_{i=1}^n \frac{k}{w_i} \leq \sum_{i=1}^n \frac{k}{i} \leq k \log n$$

What are we saying here? If we fix a value $k$ we can evaluate $[i,k]$ across all $i$ in $O(k \log n)$. $k$ is bounded by $s$ and we must evaluate this for each $k$, so the algorithm in fact runs in $O(s^2 \log n)$. Of course, if the $c_i$ are known a priori to be small, this is a very loose bound and it's more useful to think of the algorithm as running in $O(n \sum_i c_i)$. Either way it's a far cry from the $O(ns)$ solutions for the other variations.

Generally we speed up an algorithm by noticing and getting rid of redundant computation, and the case is no different here. The critical observation is best conveyed through example. Fix $w_i = 3, c_i=4$. We have

$$\begin{aligned} [i,k] &= [i-1,k] + [i-1,k-3] + [i-1,k-6] + [i-1,k-9] + [i-1,k-12]\\ [i,k+1] &= [i-1, k+1] + [i-1, k-2] + [i-1, k-5] + [i-1, k-8] + [i-1, k-11] \\ [i,k+2] &= [i-1, k+2] + [i-1, k-1] + [i-1, k-4] + [i-1, k-7] + [i-1, k-10] \end{aligned}$$

No repeated work yet. But see what happens for $[i,k+3]$:

$$[i,k+3] = [i-1, k+3] + [i-1, k] + [i-1, k-3] + [i-1, k-6] + [i-1, k-9]$$

All except the first element of this summation have been seen before in the solution to $[i,k]$! In fact we can write

$$[i,k+3] = [i-1,k+3] + [i,k] - [i,k-12]$$

What is happening here is that to solve $[i,k]$ we are use $[i,k-w_i]$ as a cumulative array that sums values across $[i-1, \dots]$. Unfixing $w_i, c_i$ we arrive at the $O(1)$ recurrence

$$[i,k] = [i-1, k] + [i,k-w_i] - [i,k-(c_i+1)w_i]$$

To me this is a remarkable result. In some informal sense the recurrence appears to "lie between" the recurrences for the 0/1 and unbounded variations, in that it starts with $[i-1, k] + [i, k-w_i]$ like the unbounded variation, but then it corrects itself using $[i-1,k-(c_i+1)w_i]$ to avoid using item $i$ too many times. It makes sense if I squint, but I don't think I could have come up with that correction procedure through combinatorial reasoning. The recurrence seems to say "we make a combination of items up to $i$ with weight $k$ either by just taking a combination of weight $k$ using items up to $i-1$, or by adding item $i$ to a combination $k-w_i$ of items up to $i$, but if it's possible that we used too many of item $i$, uncount the number of combinations that contain too many uses of that item." This gives us a $O(ns)$ time, $O(n)$ space solution to the bounded counting knapsack problem.

In the next entry we will consider another counting knapsack variation which generalises the unbounded counting knapsack problem into another "dimension", with beautiful results.

Some reasons for studying mathematics

A different sort of post today.

I've never had much interest in doing maths. I'm fundamentally incurious about mathematical objects, I don't feel joy or excitement when I prove a theorem, and I rarely enjoy the process of solving mathematical problems. The typical questions that I encounter in my course of mathematical study don't appeal to my intellectual hunger in the same way that, say, a computer science problem can take over a corner of my brain for weeks. This is how I know that mathematics is not for me. So why did I choose to study it?

When I studied algorithms in high school I learned a bunch of graph search techniques and dynamic programming and computational geometry etc, but more importantly I developed a sort of mental toolbox for reasoning about systems and processes. Items in the toolbox could also be described as "problem-solving techniques" or "approaches" or "paradigms" or "meta-algorithms"; they can feel abstract and fuzzy and difficult to describe, but certainly seem to exist as distinct concepts in my brain so that I can manifest them at will. They're obviously useful in the context of algorithmic problem solving, but they're also applicable to problems in totally unrelated fields. There are hundreds of tools in my toolbox of various levels of abstraction and usefulness. I'll have a go at identifying some of the tools in my algorithms toolbox:

• Thinking in terms of growth rates
• The notion of "greediness"
• "Sliding window" methods: we can solve a sequence of overlapping problems by solving the first problem and efficiently transforming that solution into a solution for the subsequent problem, and so on.
• Divide and conquer
• Laziness vs eagerness
• Invariants
• "X reduces to Y", meaning in a precise sense "X is at least as hard as Y, and we can't solve X without solving Y"

This mental toolbox represents the greatest benefit I received from studying algorithmic problem solving.

The study of mathematics has given me a whole new set of tools to add to my toolbox. As in the examples above, these tools tools are used constantly in the study of maths, but they have little to do with mathematical problems themselves. These mathematical tools are often far more insightful and valuable to me than the problems that they help to solve, and I can apply them to problems outside of mathematics. On its own, this is a good enough reason for me to study maths (it's not the only reason – see below). In this entry I will list and briefly describe some of the tools I've added to my toolbox in the past few years of maths education. Most of these tools can be expressed as statements about systems or structures. I may update this list over time as new tools occur to me. All tools are stated informally.

• Action to structure: we can study actions on a structure by treating the actions as first-class objects themselves. (for instance studying symmetries via the symmetry group. Reverse of "Structure to action".)
• Basis: in some cases, the behaviour of a system is wholly described by its behaviour on a small set of representative inputs.
• Continuity: the property of being able to induce arbitrarily small changes in the output of a system via correspondingly small changes in the input.
• Contradiction (proof by): one can assert a proposition by supposing that it doesn't hold and subsequently taking a logical sequence of steps to arrive at an impossibility.
• Convexity: the property of "not having dents"; not necessarily in a geometric setting.
• Density arguments: if a property holds for a subset $A$ of a set $X$ and any element of $X$ is the limit of elements in $A$, then we may be able to extend the property to all elements of $X$. (This is useful when $A$ consists of "simple" examples where the property is easily explored. The obvious example is the treatment of simple functions in integration theory.)
• Diagonalisation: if we have a countable number of instances from some class, we can generate a new instance of that class simply by ensuring that it differs from each previous instance at exactly one point.
• Differentiation: we can learn about a system by observing its instantaneous change.
• Duality: Two seemingly disparate domains sometimes pair up conveniently, so that we can study items in one domain by looking at the corresponding item in the other domain.
• Eigen-: we can learn about a system by studying the inputs on which it is particularly well-behaved.
• Exchange arguments: we can show that a state is optimal by taking any other state and showing that it doesn't get any worse if we bring it a step closer to the proposed optimal state.
• Graphs: complex systems and structures often have natural representations as networks.
• Integration: we can learn about a system by observing its cumulative effect.
• Isomorphism: rather than speaking about the equality of two objects, it's often more meaningful to consider equivalence with respect to some structural property.
• Lifts: methods that transform solutions for a problem in one domain to solutions in a broader domain, hence "lifting" them.
• Limits: properties that hold as we get close to something may or may not hold once we get there, depending on the attributes of the surrounding space.
• Locality vs globality: properties that hold at a micro level often don't hold at the macro level and vice versa.
• Orthogonality: a stronger or more precise form of independence. If A, B are orthogonal then we can change A in any way without affecting B and vice versa.
• Pigeonhole principle: if you're placing balls into boxes with more balls than boxes, at least one box will have more than one ball. (And in reverse: if there are less balls than boxes, at least one box will be empty.)
• Probabilistic method: informally, if a property has non-zero chance of holding for a randomly-selected element of a set, then the property holds for at least one element of the set. If the mean of the values in a set is $x$ then at least one of the values is at least $x$ and at least one of the values is at most $x$.
• Projection: we can learn about a complex structure by studying what it looks from a fixed perspective. I'm stunned at how often the notion of projection appears in mathematics and how abstract it can become. We can project in the literal geometric sense – as in projective geometry, or projecting a vector onto another – but we can also project the integers onto the set {0,1} via the equivalence relation of parity, or we can project a vector space to another vector space via a linear map, or we can project a sigma-algebra to a sub-sigma-algebra via conditional expectation. Projection is closely linked to the notion of homomorphism – so much so that I wish our study of homomorphisms had begun with an explicit acknowledgement that homomorphisms are precisely the maps that behave "projectively".
• Representatives: we can sometimes study the interactions of complex structures (often equivalence classes) by instead considering the interactions of individual representatives from those complex structures. (e.g. dealing with cosets via coset representatives.)
• Structure to action: we can learn about a structure by constructing an action of that structure on some class of objects and studying the action instead. (for instance studying a group via its group action on a set. The reverse of "Action to structure".)
• Substructure: we can learn about a structure by studying the smaller structures embedded within it. (see subgroup, subspace, submodule, subgraph, etc.)
• Symmetry: we can learn about a structure by studying the ways in which we can transform the structure without changing the way it looks.
• Two-player game: some notions can be understood in terms of a two-player game oppositional game. For instance epsilon-delta proofs interpreted as a game between an epsilon-picker and a delta-picker.

Why else do I study maths if I fundamentally don't care for it? I find maths to be an exceptional form of brain training. Most of this training happens via two mechanisms:

1. Structuring and summarising information: as part of my maths education I'm expected to be able to regurgitate in exam conditions something like fifty theorems and proofs per course. Each proof can consist of pages of dense equation manipulation or tricky argumentation, which is far too much to learn by rote. Instead, I need to memorise the core ideas of each theorem and proof in summarised form, but in a way where I can elaborate on those summaries as much as required until I arrive at their full unsummarised forms. I also need to develop a rich mental map of the relationships between the objects of study. Both of these processes develop my reasoning abilities.
2. Coping with abstraction: the objects of study in pure mathematics become increasingly abstract throughout a course of study, but we still need to reason about them. So we develop skills and techniques for coping with abstraction. For instance, we can fall back to geometric intuition (vector space = "space with arrows that can be added to each other and scaled"; orbit-stabiliser theorem = "sending a fixed vertex of a polyhedron to all other vertices then counting symmetries which fix that vertex in place"), as long as we keep in mind the limitations of that geometric intuition. Or we can comprehend abstraction in terms of sufficiently distinct representative examples whose generalisation "spans" the full abstraction (for instance it's difficult to understand abstract measure theory in full generality, so we develop most of our intuition by comparing and contrasting canonical examples: Lebesgue measure, counting measure, Hausdorff measure, Dirac measure, ...). If we think of an abstract concept as a space of instances of examples ("measure theory" as the study of the set of measures), this technique corresponds to choosing a set of examples at different points on the boundary of the space that help us to reason about any instance inside the boundary. I suspect that this skill is exercised more in mathematical study than in any other area of study.

Studying maths for me feels a lot like going to the gym. My interest in number theory or complex analysis is about as deep as my interest in push ups or squats. To me their value is as means of self-improvement rather than as ends in themselves. I deeply respect my lecturers and mathematically-minded peers, but as my study has progressed and the mathematical objects in consideration have become more contrived, I increasingly struggle to understand the motivation behind a career in research mathematics. But of course many academic mathematicians say the exact same thing about careers in tech or finance. Different strokes for different folks...

Developing intuition for Noether's isomorphism theorems

In abstract algebra, Noether's isomorphism theorems are extremely general statements about algebraic objects and the relationships between them. I've studied the theorems in separate courses on groups, rings, and modules over commutative rings (including vector spaces), so I've spent plenty of time trying to grasp their essence. In this post I share some of the intuition that I've developed. If nothing else, this intuition helps to remember the statement of the theorems. I will present the theorems in the context of abelian groups, where the statements are especially straight-forward, in order to avoid clouding the intuition with details about normal subgroups and so on. I will number the theorems 1, 2, 3, 4.

Theorem 1. If $f$ is a homomorphism between abelian groups $A, B$ then $$\frac{A}{\text{Ker } f} \cong \text{Im } f$$

This is the core theorem from which theorems 2 and 3 easily follow. By taking the quotient of $A$ with the kernel of $f$, we are forming a group that looks like $A$ except where we've combined the elements that together map to the same element in the codomain. I think of this as forcing $f$ to become injective by taking all the sets of elements where it fails to be injective and combining those elements into a single element. Since our version of $f$ is now injective, and it maps onto its image (by definition of image), then obviously it's a bijection onto its image. Moreover, the structure-preserving nature of the homomorphism means that structure is preserved through the bijection – i.e. there is an isomorphism between the LHS and RHS.

To me, the statement becomes almost tautological under this interpretation. It seems to say "if we take a homomorphism and disregard all the points where it fails to be injective, then it becomes injective".

Theorem 2. If $A, B$ are subgroups of an abelian group $G$, then $$\frac{A+B}{A} \cong \frac{B}{A \cap B}$$

I think of this statement as a broad generalisation of a basic result from number theory: if $a, b$ are positive integers then $a/\gcd(a,b) = \mathrm{lcm}(a,b)/b$ (this rearranges into the more familiar form $ab/\gcd(a,b) = \mathrm{lcm}(a,b)$. If we draw the diamond-shaped partial order of the "divides" relation with the positive integers $a, b, \gcd(a,b), \mathrm{lcm}(a,b)$ and then we draw the partial order of the "subgroup" relation with groups $A, B, A \cap B, A + B$, we see that the number-theoretic statement on the "divides" partial order looks very similar to what the second isomorphism theorem says about the "subgroup" partial order.

Indeed, $A \cap B$ as the "greatest common subgroup" of $A, B$ and $A + B$ is the "smallest common supergroup" of $A, B$. We know that the quotient operation feels a lot like division, so it seems natural to me to relate the quotient operation with greatest common subgroups and smallest common supergroups to the division operation with greatest common divisors and least common multiples.

The number theory statement can in fact be derived from the second isomorphism theorem applied to the $(\mathbb Z, +)$, so we're not saying anything groundbreaking. The benefit of this view of things is that you probably already have a good intuition for the properties of GCDs and LCMs, so it's just a matter of translating that intuition to the world of abstract algebra.

This isn't a complete intuition. In particular it's not clear under this intuition why we get an isomorphism between the LHS and RHS rather than just a bijection. But it's a useful starting point, or at the very least a mnemonic device.

Theorem 3. If $A \leq B \leq G$ are abelian groups, then $$\frac{G/A}{B/A} \cong \frac{G}{B}$$

This is probably the easiest theorem to remember out of the three. Quotients feel like fractions, and indeed we can "multiply both sides" of the fraction to cancel out quotients of quotients. Unfortunately I don't yet have a good intuition for the preservation of structure.

It's worth keeping in mind that the proofs of the three theorems are all in the "follow your nose" style. For the first theorem we want to establish an isomorphism between the LHS and RHS; it turns out that the most natural possible mapping from the LHS to the RHS turns out to be this isomorphism. For the remaining theorems we again pick the most natural mappings between the LHS and RHS (it'll be easiest to take the mapping from the RHS to the LHS). Again this mapping will happen to be a homomorphism, and the second and third theorems fall out immediately from applying the first theorem to this homomorphism. This is quite remarkable and it means that even if you haven't developed a fully-formed intuition for the theorems, it's easy to convince yourself rigorously that they're correct.

Theorem 4. Let $P$ be a subgroup of an abelian group G. Then there exists a bijective correspondence$$\left\{ \text{subgroups of } \frac{G}{P} \right\} \longleftrightarrow \left\{ \text{subgroups of G that include P}\right\}$$

The right-to-left map looks like this: if $Q \leq P$ is a subgroup satisfying $P \subseteq Q$, then $P$ is also a subgroup of $Q$, so we can form the quotient $Q/P$. Since $Q \leq G$ it would seem natural that $Q/P \leq G/P$. Indeed this is the case, so $Q/P$ is in the LHS. In other words, we map from the right to the left simply by taking the quotient of the subgroup with $P$.

For the left-to-right map, notice that any subgroup of $G/P$ has elements of the form $g+P$ where $g$ is a coset representative. Denote the set of all such representatives by $H$. We can show that $H$ is a subgroup of $G$ (it inherits its structure from $G/P$). Moreover, the representatives of the identity element of the subgroup of $G/P$ are the elements of $P$ itself, so the elements of $P$ are contained within $H$. Thus $P \subseteq H \leq G$ so $H$ lands in the RHS. I visualise this map as "unpacking" a subgroup of the quotient into its representatives, and it makes intuitive sense to me that (1) these representatives form a subgroup of the quotient numerator, and (2) the identity of the quotient (i.e. the coset $0+P$) is in the subgroup and unpacks to the set $P$.

Some work has to be done to prove that these maps are inverses of each other, but if we think of the maps above as packing elements into a quotient and then unpacking the quotient back to elements, the fact that the maps are inverses of each other shouldn't come as a surprise.

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