Riddler - Book Proofs

Connect the dots

This week’s Riddler Classic is a problem about connecting dots to create as many non-intersecting polygons as possible. Here is the problem:

Polly Gawn loves to play “connect the dots.” Today, she’s playing a particularly challenging version of the game, which has six unlabeled dots on the page. She would like to connect them so that they form the vertices of a hexagon. To her surprise, she finds that there are many different hexagons she can draw, each with the same six vertices.

What is the greatest possible number of unique hexagons Polly can draw using six points?

(Hint: With four points, that answer is three. That is, Polly can draw up to three quadrilaterals, as long as one of the points lies inside the triangle formed by the other three. Otherwise, Polly would only be able to draw one quadrilateral.)

Extra Credit: What is the greatest possible number of unique heptagons Polly can draw using seven points?

Here is my solution:
[Show Solution]

It turns out this is a well-studied and notoriously difficult problem in combinatorial geometry. While it is tempting to ask the same question for 8, 9, or even $n$ points, we’ll see that the problem gets very difficult very quickly, and the answer for general $n$ is not actually known!

First, we’ll get some terminology out of the way:

A complete graph is a graph with every possible edge drawn in. The complete graph with $n$ nodes is denoted $K_n$.
A Hamiltonian cycle is a path through the graph that finishes where it started and visits every node exactly once. For the graph $K_n$, there are $\frac{1}{2}(n-1)!$ possible Hamiltonian cycles since we can fix the starting node and then there are $n-1$ nodes left to order. We divide by two because each cycle is counted twice (can be traversed forwards or backwards).
A realization of a graph is a particular way of arranging the nodes. Realizations are important because we care about whether edges cross or not; sometimes two different arrangements of the same graph will have different numbers of edge crossings.
The rectilinear crossings of a graph realization $G$ is the number of times edges cross when we connect the nodes with straight lines. We’ll denote this number $\mathrm{rcr}(G)$.
The crossing-free Hamiltonian cycles of a graph realization $G$ is the number of different Hamiltonian cycles of $G$ that do not contain any edges that cross each other. We’ll denote this number $\mathrm{cfhc}(G)$. These are also called “spanning cycles” or “simple polygonalizations” depending on who you ask.

For a simple example, let’s start with $K_4$, the complete graph on 4 nodes. Although there are infinitely many ways to realize this graph (since we can arrange the 4 points in infinitely many ways), rearrangements of the points that preserve how the edges intersect one another are all equivalent for our purpose (this is an example of an equivalence class). So we only need to consider one representative from each class. These representative realizations are called order types. It turns out $K_4$ has two order types:

The first order type has no crossings ($\mathrm{rcr}(G)=0$) and three possible crossing-free Hamiltonian cycles. Here are the cycles:

The second order type has one crossing ($\mathrm{rcr}(G)=1$) and only one possible crossing-free Hamiltonian cycle:

The bad news…

Unfortunately, things get difficult from this point on as we add more nodes:

The number of order types for $K_n$ grows rapidly with $n$. For $n\ge 3$, they are: $\{1, 2, 3, 16, 135, 3315, 158817, 14309547, 2334512907,\dots\}$. This sequence is in OEIS. There is no known general formula.
The minimal rectilinear crossing number for $K_n$ over all possible realizations for $n\ge 3$ is: $\{0, 0, 1, 3, 9, 19, 36, 62, 102, 153,\dots\}$. This sequence is also in OEIS. No known formula known for this one either.
Finally, the minimal number of crossing-free Hamiltonian cycles for $K_n$ for $n \ge 3$ is: $\{1, 3, 8, 29, 92, 339, 1282, 4994,\dots\}$. And this is also in OEIS. You guessed it; no known formula.

These problems are related but different. For example, we might expect realizations with smaller crossing numbers to contain more crossing-free Hamiltonian cycles, since it is easier to avoid intersections when there are fewer of them. But this is not always the case. We might also expect realizations with more symmetry to contain more cycles; also not true in general. What I’m trying to get at is that there are no (known) tricks we can use to reduce this problem to a simpler one. Unfortunately, it appears the only way to solve such problems is to enumerate all order types (even this is difficult!), then enumerate all possible cycles, keeping only the ones that are crossing-free.

Many other variants are just as difficult: computing the crossing number (with curved edges allowed), the number of triangulations, the minimum number of convex polygons needed for a decomposition of the convex hull, and more. These sorts of problems belong to a branch of mathematics called combinatorial geometry. A great deal of research on the topic of order types, crossing numbers, and more, has been conducted by Prof. Oswin Aichholzer and colleagues. If you’re interested in learning more, I recommend checking out his webpage here, which contains an up-to-date database of solutions to many of these and related problems for small-ish $n$.

Visualizing the solutions

I used the point layouts provided on Prof. Aichholzer’s website and wrote some Python code to visualize the results. As a warm-up, I started with the case $K_5$.

The order types are sorted by their rectilinear crossing number. The one with the most crossing-free Hamiltonian cycles is the first one. Here they are:

Six nodes

Here are the 16 order types for $K_6$:

Here are the 29 crossing-free Hamiltonian cycles for the maximal configuration:

Seven nodes

Here are the 135 order types for $K_7$. Interestingly, there are three different realizations for the minimal rectilinear crossing number of $9$, and the most symmetric one (the first one) is not the one with the most crossing-free Hamiltonian cycles!

And here are the 92 crossing-free Hamiltonian cycles for the maximal configuration:

More nodes

The solutions get too large to visualize beyond 7 nodes, but here are the results for $n \leq 10$, courtesy once again of Prof. Aichholzer’s webpage. “CFHC” stands for “crossing-free Hamiltonian cycles”.

Number of nodes	Order types	Max CFHC
3	1	1
4	2	3
5	3	8
6	16	29
7	135	92
8	3,315	339
9	158,817	1,282
10	14,309,547	4,994

The Python code I wrote to produce all the figures is available here.

Dungeons & Dragons

This week’s Riddler Classic is a probability problem about the game Dungeons & Dragons. Here it goes:

When you roll a die “with advantage,” you roll the die twice and keep the higher result. Rolling “with disadvantage” is similar, except you keep the lower result instead. The rules further specify that when a player rolls with both advantage and disadvantage, they cancel out, and the player rolls a single die. Yawn!

There are two other, more mathematically interesting ways that advantage and disadvantage could be combined. First, you could have “advantage of disadvantage,” meaning you roll twice with disadvantage and then keep the higher result. Or, you could have “disadvantage of advantage,” meaning you roll twice with advantage and then keep the lower result. With a fair 20-sided die, which situation produces the highest expected roll: advantage of disadvantage, disadvantage of advantage or rolling a single die?

Extra Credit: Instead of maximizing your expected roll, suppose you need to roll N or better with your 20-sided die. For each value of N, is it better to use advantage of disadvantage, disadvantage of advantage or rolling a single die?

Here is a detailed derivation of the relevant probabilities:
[Show Solution]

We’ll use the following notation:
\begin{align}
f_k &= \mathbf{Prob}(X=k) & &\text{(probability mass function, pmf)} \\
F_k &= \mathbf{Prob}(X\leq k) & &\text{(cumulative mass function, cmf)} \\
\bar F_k &= \mathbf{Prob}(X\geq k) & &\text{(complementary cmf)}
\end{align}where $X$ is the outcome of a single roll of a fair $n$-sided die. Since it’s a fair die, these are:
\[
f_k=\frac{1}{n}, \qquad F_k=\frac{k}{n}, \qquad \bar F_k=\frac{n-k+1}{n}.
\]Let’s tackle the more complicated cases one by one. For the other events of interest, we’ll use similar notation but with superscripts “A” (advantage), “D” (disadvantage), “AD” (advantage of disadvantage), and “DA” (disadvantage of advantage). So, for example, $f^\textrm{DA}_k$ is the pmf for rolling with disadvantage of advantage.

Rolling with advantage

Here, we roll the die twice and take the largest value. We can leverage the following property of the maximum function: $\max(x,y) \le k$ if and only if $x\le k$ and $y\le k$. If $X_1$ and $X_2$ are independent rolls of an $n$-sided die, we have:
\begin{align}
F^{\mathrm{A}}_k &= \mathbf{Prob}(\max(X_1,X_2) \le k) \\
&= \mathbf{Prob}(X_1\le k,X_2\le k) \\
&= \mathbf{Prob}(X_1\le k)\,\mathbf{Prob}(X_2\le k) \\
&= F_k^2 \\
&= \left(\frac{k}{n}\right)^2
\end{align}Finally, we can use the fact that $f^{\mathrm{A}}_k = F^{\mathrm{A}}_k-F^{\mathrm{A}}_{k-1}$ to compute the pmf:
\[
f^{\mathrm{A}}_k = \frac{2k-1}{n^2}
\]

Rolling with disadvantage

Here, we use a similar approach, but with the complementary cmf instead:
\begin{align}
\bar F^{\mathrm{D}}_k &= \mathbf{Prob}(\min(X_1,X_2) \ge k) \\
&= \mathbf{Prob}(X_1\ge k,X_2\ge k) \\
&= \mathbf{Prob}(X_1\ge k)\mathbf{Prob}(X_2\ge k) \\
&= \bar F_k^2 \\
&= \left(\frac{n-k+1}{n}\right)^2
\end{align}Finally, we can use the fact that $f^{\mathrm{D}}_k = \bar F^{\mathrm{D}}_k-\bar F^{\mathrm{D}}_{k+1}$ to compute the pmf:
\[
f^{\mathrm{D}}_k = \frac{2n-2k+1}{n^2}
\]

Rolling with advantage of disadvantage

This case can be computed by iterating the strategies used in the previous two cases. Specifically:
\begin{align}
F^{\mathrm{AD}}_k &= \left( F^{\mathrm{D}}_k \right)^2 \\
&= \left( 1-\bar F^{\mathrm{D}}_{k+1} \right)^2 \\
&= \left( 1-\left(\frac{n-k}{n}\right)^2 \right)^2 \\
&= \frac{k^2(2n-k)^2}{n^4}
\end{align}As before, we can compute the pmf using $f^{\mathrm{AD}}_k = F^{\mathrm{AD}}_k-F^{\mathrm{AD}}_{k-1}$:
\[
f^{\mathrm{AD}}_k = \frac{(2 k-2 n-1) \left(2 k^2-4 k n-2 k+2 n+1\right)}{n^4}
\]

Rolling with disadvantage of advantage

This case can also be computed by iterating the previous cases. Specifically:
\begin{align}
\bar F^{\mathrm{DA}}_k &= \left( \bar F^{\mathrm{A}}_k \right)^2 \\
&= \left( 1-F^{\mathrm{A}}_{k-1} \right)^2 \\
&= \left( 1-\left(\frac{k-1}{n}\right)^2 \right)^2 \\
&= \frac{(k-n-1)^2 (k+n-1)^2}{n^4}
\end{align}As before, we can compute the pmf using $f^{\mathrm{DA}}_k = \bar F^{\mathrm{DA}}_k-\bar F^{\mathrm{DA}}_{k+1}$:
\[
f^{\mathrm{DA}}_k = \frac{(1-2 k) \left(2 k^2-2 k-2 n^2+1\right)}{n^4}
\]

Continuous limit

We can also compute the continuous limit as $n$ tends to infinity. To make sense of this, we’ll imagine that “rolling an $n$-sided die” is replaced by “picking a uniformly chosen random number in $[0,1]$”. So the probability mass functions are replaced by probability densities. The derivations are essentially identical to those above, except the densities for the base case are:
\[
f(z)=1, \qquad F(z)=z, \qquad \bar F_k=1-z.
\]

And here are the results:
[Show Solution]

Table of results

Here is a table containing the probability mass function (pmf), cumulative distribution (cmf), and complimentary cumulative distribution (ccmf) for the five combinations of advantage and disadvantaged dice rolls.
\[
\begin{array}{c|c|c|c}
\text{case} & \text{pmf }f_k & \text{cmf }F_k & \text{ccmf }\bar F_k \\ \hline
\text{A}+\text{D} & \frac{1}{n} & \frac{k}{n} & \frac{n-k+1}{n} \\
\text{A} & \frac{2k-1}{n^2} & \frac{k^2}{n^2} & \frac{(n+k-1)(n-k+1)}{n^2} \\
\text{D} & \frac{2n-2k+1}{n^2} & \frac{k(2n-k)}{n^2} & \frac{(n-k+1)^2}{n^2} \\
\text{A of D} & \frac{(2 k-2 n-1) \left(2 k^2-4 k n-2 k+2 n+1\right)}{n^4} & \frac{k^2(2n-k)^2}{n^4} & \frac{(n-k+1)^2 \left(n^2-k^2+2kn+2k-2 n-1\right)}{n^4}\\
\text{D of A} & \frac{(2k-1) \left(2n^2-2k^2+2k-1\right)}{n^4} & \frac{k^2 \left(2 n^2-k^2\right)}{n^4} & \frac{(n+k-1)^2 (n-k+1)^2}{n^4} \\
\end{array}
\]And here is a table for the continuous case, containing the probability density function (pdf), cumulative distribution (cdf), and complimentary distribution (ccdf).
\[
\begin{array}{c|c|c|c}
\text{case} & \text{pdf }f(z) & \text{cdf }F(z) & \text{ccdf }\bar F(z) \\ \hline
\text{A}+\text{D} & 1 & z & 1-z \\
\text{A} & 2z & z^2 & 1-z^2 \\
\text{D} & 2-2z & z(2-z) & (1-z)^2 \\
\text{A of D} & 4z(1-z)(2-z) & z^2(2-z)^2 & (1-z)^2(1+2z-z^2)\\
\text{D of A} & 4z(1-z^2) & z^2(2-z^2) & (1-z^2)^2 \\
\end{array}
\]

Expected value

We can compute expected values as a function of $n$ by evaluating $\sum_{k=1}^n k\,f_k$ for each of the cases. We can also evaluate the continuous limit by computing $\int_0^1 z\,f(z)\,\mathrm{d}z$ for each case. I rearranged the columns so that they are sorted from largest to smallest.
\[
\begin{array}{c|c}
\text{case} & \text{Expectation} & \text{Continuous limit} \\ \hline
\text{A} & \frac{(n+1) (4 n-1)}{6 n} & \frac{2}{3} \\
\text{D of A} & \frac{(n+1) \left(16 n^3-n^2+n-1\right)}{30 n^3} & \frac{8}{15} \\
\text{A}+\text{D} & \frac{n+1}{2} & \frac{1}{2} \\
\text{A of D} & \frac{(n+1) (2 n+1) \left(7 n^2-3 n+1\right)}{30 n^3} & \frac{7}{15} \\
\text{D} & \frac{(n+1) (2 n+1)}{6 n} & \frac{1}{3} \\
\end{array}
\]Indeed, if we take the Expectation column $E$, which takes on values in $\{1,\dots,n\}$, compute $\frac{E-1}{n-1}$ in order to normalize the range of possible outcomes to $[0,1]$, then let $n\to\infty$, we recover the “Continuous limit” column. Here is a plot of this phenomenon in action.

Ultimately, we conclude that “disadvantage of advantage” is better than “advantage of disadvantage”. This idea that the “minimums of maximums” is always larger than the “maximum of minimums” is actually a much more general notion that applies to arbitrary functions. The result is known as the max-min inequality and it states that for any function $f:Z\times W \to \mathbb{R}$, we have:
\[
\inf_{w\in W}\sup_{z\in Z}f(z,w) \geq \sup_{z\in Z}\inf_{w\in W}f(z,w)
\]

Rolling $k$ or better

If the task is to roll $k$ or better, then we are seeking the probability that our outcome $X$ satisfies $X \ge k$. This is precisely the complimentary cumulative distribution function! Plotting these values for $n=6$, we obtain:

As anticipated, when $k=1$, all probabilities are $1$. The “advantage” roll is always best, while “disadvantage” is always worst. “disadvantage of advantage is second best for $k=1,2,3,4$, but for $k=5,6$ it’s best to roll a single die. We observe a similar result when $n=20$:

This time, “disadvantage of advantage” is second best for $1\le k \le 13$ and rolling a single die is best for $14 \le k \le 20$. In general, we can compute where the transition occurs as a function of $n$ by seeing which value of $k$ leads to $\bar F^{\text{DA}}_k = \bar F_k$. This leads to the equation:
\[
\frac{(n+k-1)^2 (n-k+1)^2}{n^4} = \frac{n-k+1}{n}
\]Solving this equation for $k$ and discarding extraneous solutions, we find the threshold occurs at:
\[
k_0 = 1 + \left( \frac{\sqrt{5}-1}{2} \right) n
\]So if $1\leq k < k0$, we have $\bar F^{\text{DA}}_k \gt \bar F_k$ and if $k_0 \lt k \leq n$, we have $\bar F^{\text{DA}}_k \lt \bar F_k$. Ok, so what happens in the limit? We can plot the continuous complimentary cdfs and obtain:

The pattern is the same, and we can solve for the cross-over point as before, by setting $\bar F^{\text{DA}}(z) = \bar F(z)$. This time, the equation is $(1-z^2)^2=1-z$. Solving for $z$, we obtain
\[
z_0 = \frac{\sqrt{5}-1}{2}
\]This is not surprising, as we could have simply normalized the discrete result as we did for expectations and take the limit $n\to\infty$ to obtain the same result. This number is the inverse of the Golden Ratio, neat!

Flip to freedom

This week’s Riddler Classic is a problem about coin flipping. The text of the original problem is quite long, so I will paraphrase it here:

There are $n$ prisoners, each with access to a random number generator (generates uniform random numbers in $[0,1]$) and a fair coin. Each prisoner is given the opportunity to flip their coin once if they so choose. If all of the flipped coins come up Heads, all prisoners are released. But if any of the flipped coins come up Tails, or if no coins are flipped at all, the prisoners are not released. If the prisoners cannot communicate in any way and do not know who is flipping their coin or not, how can they maximize their chances of being released?

Here is my solution:
[Show Solution]

This problem is a balancing act. If $k$ of the prisoners flip their coins, the chance of being released is
\[
\mathrm{Prob}(\text{win if }k\text{ coins are flipped})=\begin{cases} 0 & k=0\\ \frac{1}{2^k} & k=1,2,\dots,n\end{cases}
\]The best-case scenario is if $k=1$ prisoner flips their coin and the others do not, ensuring the maximum probability of release of $1/2$. But how should the prisoners decide who flips their coin? If they were allowed to communicate ahead of time, they could elect somebody to flip their coin. But since the prisoners cannot communicate so there is no way to plan ahead!

Since each prisoner has access to a random number generator, they can use it to decide whether they should be flipping their coin or not. Here is the strategy. We must decide on a threshold $t$. Each prisoner queries their random number generator and if the number they obtain is less than $t$, they flip their coin. The larger $t$ is, the more likely prisoners are to flip coins. So there is a trade-off between ensuring at least one prisoner flips a coin and ensuring not too many prisoners flip coins. Let’s figure out the best choice of $t$.

Let’s call the probability the prisoners are released $f_n(t)$. It depends on the threshold they choose. Now,
\begin{align*}
f_n(t) &= \sum_{k=1}^n \frac{1}{2^k} \mathrm{Prob}(\text{there are }k\text{ flips}) \\
&= \sum_{k=1}^n \frac{1}{2^k} \binom{n}{k}t^k(1-t)^{n-k} \\
&= \sum_{k=1}^n \binom{n}{k}\left(\frac{t}{2}\right)^k (1-t)^{n-k} \\
&= ( 1-\tfrac{1}{2}t)^n-(1-t)^n
\end{align*} In the last step, we used the Binomial Theorem to evaluate the sum. So our goal is to find $t\in[0,1]$ that minimizes $f_n(t)$. Here is what the function looks like for different values of $n$:

As anticipated, there is an ideal value of the threshold $t$ that maximizes the probability of winning. However, the maximum depends on $n$, the number of prisoners. To find the maximum, we can use calculus: set the derivative of $f_n(t)$ equal to zero:
\[
\frac{\mathrm{d}f_n(t)}{\mathrm{d}t}=n(1-t)^{n-1}-\tfrac{1}{2}n(1-\tfrac{1}{2}t)^{n-1}=0
\]Solving for the optimal threshold $t$, we obtain:

$\displaystyle
t_\text{opt}=1-\frac{1}{2^{n/(n-1)}-1}
$

As a function of $n$, $t_\text{opt}$ is the threshold each prisoner should use. If their random number generator produces a value less than the threshold, they should flip their coin. Otherwise, they should not. As anticipated, in the limit $n\to\infty$, $t_\text{opt}\to 0$ because as the number of prisoners increases, we want each prisoner to be less likely to flip their coin. As an approximation, the first order asymptotic expansion of $t_\text{opt}$ is
\[
t_\text{opt} = \frac{2\log(2)}{n} + O\left(\frac{1}{n^2}\right)\approx\frac{1.386}{n}
\]So this can be used as a quick way to approximate the threshold (in case the prisoners don’t have access to a calculator!).

How about the probability of winning assuming the prisoners use this strategy? This simply amounts to evaluating $f_n(t_\text{opt})$. Doing so, we obtain:

$\displaystyle
\textrm{Prob}(\text{win}) = \frac{1}{\left(2^{n/(n-1)}-1\right)^{n-1}}
$

Here is a table comparing the optimal threshold and corresponding probability of winning for different values of $n$:

Number of prisoners	threshold $t_\text{opt}$	Probability of freedom
1	1	0.5
2	0.666667	0.333333
3	0.453082	0.299119
4	0.342037	0.284842
5	0.274529	0.277001
10	0.13802	0.262709
50	0.0277034	0.252429
100	0.0138574	0.251208

So when $n=1$, the lone prisoner flips every time and the probability of winning is $\frac{1}{2}$. With two prisoners, each prisoner should flip $\frac{2}{3}$ of the time, and the probability of winning is $\frac{1}{3}$. As $n\to\infty$, the threshold goes to zero as described above, and we can also verify that the probability of winning tends to $\frac{1}{4}$ in the limit.

So this strategy ensures that the probability of all prisoners being released is always at least 25%. The only caveat is that in order to compute which threshold to use, the prisoners must know $n$, the total number of prisoners playing the game.

When did the snow start?

This week’s Riddler Classic is a neat calculus problem:

One morning, it starts snowing. The snow falls at a constant rate, and it continues the rest of the day.

At noon, a snowplow begins to clear the road. The more snow there is on the ground, the slower the plow moves. In fact, the plow’s speed is inversely proportional to the depth of the snow — if you were to double the amount of snow on the ground, the plow would move half as fast.

In its first hour on the road, the plow travels 2 miles. In the second hour, the plow travels only 1 mile.

When did it start snowing?

Here is my solution:
[Show Solution]

Flipping your way to victory

This week’s Riddler Classic concerns a paradoxical coin-flipping game:

You have two fair coins, labeled A and B. When you flip coin A, you get 1 point if it comes up heads, but you lose 1 point if it comes up tails. Coin B is worth twice as much — when you flip coin B, you get 2 points if it comes up heads, but you lose 2 points if it comes up tails.

To play the game, you make a total of 100 flips. For each flip, you can choose either coin, and you know the outcomes of all the previous flips. In order to win, you must finish with a positive total score. In your eyes, finishing with 2 points is just as good as finishing with 200 points — any positive score is a win. (By the same token, finishing with 0 or −2 points is just as bad as finishing with −200 points.)

If you optimize your strategy, what percentage of games will you win? (Remember, one game consists of 100 coin flips.)

Extra credit: What if coin A isn’t fair (but coin B is still fair)? That is, if coin A comes up heads with probability p and you optimize your strategy, what percentage of games will you win?

Here is my solution:
[Show Solution]

At first, the problem appears to be asking the impossible. How can flipping fair coins over and over lead to anything other than winning half the time? The key here is that “winning” only requires getting a positive score. If we looked at the “expected total score” after 100 flips, this would be 0 no matter what strategy was used (you can see this by applying linearity of expectation). But since we’re counting number of positive outcomes rather than total score, it’s possible to win more than half of the time on average.

As a simpler example of how this can happen, consider a game where you roll a die. If you get a 1, you lose \$5. But if you get anything else, you win \$1. The expected score in the game is $\tfrac{1}{6}(-5) + \tfrac{5}{6}(+1) = 0$. However, if “winning” simply means obtaining a positive score, then we win $\tfrac{5}{6}$ of the time. The flipping game is similar; our total score can average out to zero, even though it ends up being positive more than half of the time. Practically speaking, this means that when we do lose, we lose big.

Probability of winning assuming optimal play

Let’s suppose we are doing $N$ flips, and we are choosing between:

Coin $A$: wins $+1$ with prob $a$ and wins $-1$ with prob $1-a$.
Coin $B$: wins $+2$ with prob $b$ and wins $-2$ with prob $1-b$.

Our task is to determine the percentage of games that we will win if we optimize our strategy. One way to approach this is using dynamic programming. The basic idea is that our best move at any point in the game should only depend on (1) our current score and (2) how many moves we have left. Therefore, we define:
\[
V_t(s) = \left\{\begin{aligned}
&\text{Probability that we win if we play optimally}\\
&\text{from here on out, given that our current score}\\
&\text{is }s\text{ and we are currently at turn }t\text{ of }N
\end{aligned}\right.
\]We can view each $V_t$ as a vector indexed by $s \in \{\dots,-1,0,1,\dots\}$, which is our current score. The idea is to solve for each $V_t$ by working our way backwards from the final step until we reach $V_0$, which is our probability of winning for each different starting score. Ultimately, the question is asking us to compute $V_0(0)$.

The terminal distribution is simply given by our probability of winning given that we have achieved a given final score. The result in this case is simple, since at that point there are no moves left to make:
\[
V_N(s) = \begin{cases}
0 & \text{if }s \leq 0\\
1 & \text{if }s > 0
\end{cases}
\]At each previous step, we pick the coin that will give us our highest chance of winning. This leads to the recursive formula for $t=n-1,\dots,2,1$:
\[
V_{t-1}(s) = \max\left\{ a V_t(s+1)+(1-a) V_t(s-1),\, b V_t(s+2)+(1-b) V_t(s-2) \right\}
\]I couldn’t find a way to evaluate this expression in closed form, but it is rather straightforward to evaluate it numerically. Here is efficient Python code that evaluates $V_0(0)$ for any choice of $(a,b,N)$.

import numpy as np

# probability of winning if coin A has prob a of +1 and (1-a) of -1,
# coin B has prob b of +2 and (1-b) of -2, and we do N flips in total
def compute_win_prob( a, b, N ):

    # shift the indices so that a score of "s" is at index v[s+2*N]    
    v = np.zeros(4*N+1)
    
    # initialize (time=N)
    v[:2*N] = 0
    v[2*N+1:] = 1
    
    # overwrite results as we recurse backwards to time=0
    # this code is efficient: uses vectorization and overwrites past results
    for t in range(N):
        coin_A = a*v[3:4*N] + (1-a)*v[1:4*N-2]
        coin_B = b*v[4:4*N+1] + (1-b)*v[0:4*N-3]
        v[2:4*N-1] = np.maximum(coin_A, coin_B)
    
    # return the probability of winning assuming we start with a score of 0
    return v[2*N]

When we run this script using $a=b=0.5$ and $N=100$, we obtain a probability of winning of about $0.6403$. We can get the exact answer as well using the following Mathematica code:

ComputeWinProb[a_,b_,n_]:=(
  v = ConstantArray[0,4n+1];
  v[[2n+2;;4n+1]] = 1;
  For[t=0, t<n, t++,
    CoinA = a v[[4;;4n]] + (1-a) v[[2;;4n-2]];
    CoinB = b v[[5;;4n+1]] + (1-b) v[[1;;4n-3]];
    v[[3;;4n-1]] = MapThread[Max,{CoinA,CoinB}];
  ];
  v[[2n+1]])
ComputeWinProb[1/2,1/2,100]

This gives us the exact result:
\[
\text{Prob}(\text{win})=\tfrac{811698796376000066208208781649}{1267650600228229401496703205376} \approx 0.640317
\]

Does doing more flips help?

We can win 64% of the time if we do 100 flips, but What if we did 1000 flips, or more? Could we keep increasing our chances of winning? It doesn’t appear to be the case. As we do more and more flips (increase $N$), the probability of winning increases slightly, but only up to a limit. To see this, I solved the problem for $N$ ranging up to 100,000 flips. Here is what I found:

Given how slowly the function is growing, it appears that there should be a finite limit, and that this limit is around $\frac{2}{3}$ (the last point on the plot is at $0.6658$). But is there another way we can compute this limit?

Limiting behavior for a large number of flips

The optimal strategy for the case $a=b=\tfrac{1}{2}$ is to flip coin A when our score is positive and to flip coin B when our score is negative. How often will we win the game as $N\to\infty$ if we use this strategy?

We can think of the coin-flipping problem in the context of random walks. It is well-known that random walks on the 1D line (with equal probability of moving left or right at each step) will visit every integer an infinite number of times, given enough time. It is also known that the expected number of steps before a walk first returns to its starting point is infinite. Therefore, as $N$ gets large, we can expect that the the score will spend only an infinitesimal fraction of the time near 0.

We can model this process as a Markov chain with states:
\[
\dots,-8,-6,-4,-2,0,1,2,3,4,\dots
\]When our score is $\leq 0$, we will flip coin B and move $\pm 2$. When our score is $\gt 0$, we will flip coin A and move $\pm 1$. Note that the states $-3,-5,-7,\dots$ can never be reached because we only move by $\pm 2$ when our score is negative. Based on the random walk argument above, we can expect the Markov chain to spend most of its time with a large-magnitude score (either on the positive or negative side).

As an approximation, we can imagine that once we hit $+2$, we will spend a large amount of time on the positive side (with a winning score) before we return to $+1$. Likewise, once we hit $0$, we will spend a large amount of time on the negative side (with a losing score) before we return. Since the positive and negative halves of this infinite Markov chain have identical probability distributions, we can expect the time spent on both halves to be equally large. We model this situation by truncating the Markov chain and assigning a transition probability of $\gamma \gg 0$ of remaining in a winning or losing state once we get there. Here is the resulting truncated Markov chain:

The transition matrix for this Markov chain is:
\[
A = \begin{bmatrix}
\gamma & \frac{1}{2} & 0 \\
0 & 0 & 1-\gamma \\
1-\gamma & \frac{1}{2} & \gamma
\end{bmatrix}
\]The steady-state distribution (right eigenvector corresponding to the eigenvalue of $1$) is given by: $\begin{bmatrix} \tfrac{1}{2} & 1-\gamma & 1\end{bmatrix}$. As anticipated, as $\gamma \to 1$, less and less probability is associated with the intermediate state 1. In the limit, the fraction of time spent in the winning state is
\[
\text{Prob}(\text{win as }N\to\infty) = \frac{1}{\tfrac{1}{2}+1} = \frac{2}{3}
\]So this justifies why we found the probability of winning tending to $\tfrac{2}{3}$ in the numerical simulation above.

As an aside, @electricvishnu with some help from @DarkAdonisSA found a formula that computes the probability of winning the game as a function of $N$. This formula is:

$\displaystyle
\text{Prob}(\text{win in }N\text{ steps}) = \frac{2}{3} + \frac{(-1)^N}{2^N}
\left[ \frac{1}{3} + \sum_{k=1}^{N-1} (-1)^k \binom{k}{\lfloor k/2 \rfloor} \right]
$

We don’t have a proof of formula as of yet, but it does produce the exact answer when $N=100$, and it also clearly reveals the limit of $\tfrac{2}{3}$ since the second term tends to zero as $N\to\infty$. There does not appear to be a meaningful way to simplify this formula either.

Changing the probabilities

We can also examine what happens when we try varying the probabilities away from $a=b=\tfrac{1}{2}$. Using the Python function shown above, this is a straightforward exercise. Here is what we obtain when we assume coin B is fair but coin A is not fair:

As expected, when coin A has a low probability, the optimal strategy is to flip coin B a lot more, which just leads to a probability of 0.5 of winning. As coin A’s probability increases, as do our chances of winning. We can also see what happens as we increase either $a$ or $b$ (neither coin is a fair coin). Here is what we obtain:

Once again, I struggled to find closed-form expressions for any of the above cases ($a$ and/or $b$ varying). Even for simpler instances, the solutions get complicated in a hurry. For example, the probability of winning as a function of $a$, with $b=\tfrac{1}{2}$, when the game only lasts 4 flips is:
\[
\text{Prob}(\text{win}) = \begin{cases}
\frac{1}{8}(a+4) & 0 \leq a \leq \tfrac{1}{2} \\
-\frac{1}{8} \left(16 a^4-32 a^3+20 a^2-11 a-1\right) & \tfrac{1}{2} \leq a \leq \gamma \\
-\frac{1}{2} a \left(6 a^3-11 a^2+6 a-3\right) & \gamma \leq a \leq 1
\end{cases}
\]where $\gamma \approx 0.7328$ is one of the roots of $8\gamma^3-4\gamma^2-1=0$. (this was found using the Mathematica function I included above). Each time we increase the number of flips, the number of parts in the piecewise polynomial can potentially double. So it seems highly unlikely that we’ll have a simple closed-form expression by the time we arrive at $N=100$!

Penny Pinching

This week’s Riddler Classic is indeed a classic! Here it goes (paraphrased to make it a bit more general):

The game starts with $n$ pennies, which I then divide into two piles any way I like. Then we alternate taking turns, with you first, until someone wins the game. For each turn, a player may take any number of pennies he or she likes from either pile, or instead take the same number of pennies from both piles. Each player must also take at least one penny every turn. The winner of the game is the one who takes the last penny.

If we both play optimally, what starting numbers of pennies guarantee that you can win the game?

Here is my solution:
[Show Solution]

This is a spin on the classical game called Nim, which dates back to the beginning of the 16th century, but it wasn’t until 1901 when Charles L. Bouton of Harvard came up with the name “Nim” and developed the complete mathematical theory surrounding the game. In classical Nim, there can be many piles of coins, but you can only remove coins from one of the piles during your turn. The “Penny Pinching” variant in this week’s Riddler Classic is known as Wythoff’s Game and the theory surrounding it is quite fascinating. I’ll try to explain it in my own words here.

Each position can be described by a pair $(a,b)$ (the number of pennies in each pile). A position is called “safe” if it’s safe for me to leave the coins in this particular arrangement; no matter what the other player does, they can never win the game. By default, $(0,0)$ is safe, since reaching this position means I won the game.

For example, $(1,2)$ is safe because my opponent’s only legal moves are:

Remove from the first pile to obtain $(0,2)$.
Remove from the second pile to obtain $(1,1)$ or $(1,0)$.
Remove from both piles to obtain $(0,1)$.

In all of these cases, we can win in one move, so $(1,2)$ is a safe position for us. There are two key properties of safe positions:

It is impossible to reach a safe position from another safe position in one move. If this were possible, then our opponent could move to a safe position and then they would be the safe ones!
From any unsafe position, it is possible to move to a safe position in one move. This is because if we start from a safe position, then by definition we must be able to return to a safe position whenever it’s our turn again.

One way to generate safe positions is to draw a lattice. Each cell $(a,b)$ corresponds to a position of the game. We can color each cell green if it is safe and red if it is unsafe. We can do this iteratively, starting by marking $(0,0)$ safe and then proceeding as described in the following diagram:

Here is Python snippet that generates a solution grid of any size:

import numpy as np
N = 51  # change this for a different number of points
A = np.zeros((2*N,2*N),int)
for i in range(2*N):
    for j in range(i+1):
        A[i-j,j] = ~( (A[:i-j,j]==1).any() |
                      (A[i-j,:j]==1).any() |
                      (np.diag(A[i-j-1::-1,j-1::-1])==1).any() )

We can then visualize the grid using the Python code below:

import matplotlib.pyplot as plt
import matplotlib.cm as cm

fig,ax = plt.subplots(figsize=(9,9))
ax.imshow(A[:N,:N], origin='lower', interpolation=None, cmap=cm.Blues)
ax.set_xticks(range(N), minor=True)
ax.set_yticks(range(N), minor=True)
ax.grid(which='minor', color='lightgray')
ax.grid(which='major', color='dimgray')
ax.plot([0,20],[20,0],'r--')
ax.plot([0,30],[30,0],'r--')
ax.set_xlabel('coins in first pile')
ax.set_ylabel('coins in second pile')
ax.set_title("Safe combinations in Pinching Pennies (Wythoff's game)")

The resulting figure (for $N=51$) is:

When the game begins, I get to split the coins any way they like. For example, if there are 20 pennies, then I get to pick any starting position along the lower red dashed diagonal (these are the positions for which the total number of coins in both piles is 20). Since there are no safe positions along this diagonal, I cannot play safe, so I will lose. However, if the initial number of pennies is: 3, 8, 11, 16, 21, 24, 29, 32, 37, 42, 45, 50, 55, 58, 63, 66, 71, 76, 79, 84, 87, 92, 97, 100,… (the diagonals that contain safe positions) then I can split the coins to reach a safe position and therefore win. This sequence of numbers is interesting enough that it got its own name: the Wythoff AB numbers.

The distribution of safe positions seems to form lines on the plot above, but upon closer inspection, there is no obvious pattern in how the points are distributed. However, Wythoff was able to prove that the slopes of these lines are $\varphi$ and $\varphi^{-1}$, where $\varphi= \tfrac{1}{2}(1+\sqrt{5}) \approx 1.618$ is the golden ratio! Not only that, we can generate all safe pairs using the formula:

$\displaystyle
\text{The set of all safe pairs }(a,b)\text{ with }a \le b\text{ is}\quad\left( \bigl\lfloor n \varphi\bigr\rfloor,\, \bigl\lfloor n \varphi^2 \bigr\rfloor \right)\quad\text{for }n=0,1,2,\dots
$

Where $\lfloor x \rfloor$ is the floor function, i.e. the integer you get by rounding down. This generates the steeper line on the plot above. We can get the other half of the points by swapping the two coordinates. This means that the starting player can always win the game when the starting number of pennies is of the form:

$\displaystyle
\bigl\lfloor n \varphi\bigr\rfloor + \bigl\lfloor n \varphi^2 \bigr\rfloor \quad\text{for }n=0,1,2,\dots
$

When we substitute the formula above using for example the Python code below, we obtain Wythoff’s sequence:

import numpy as np
phi = (1+np.sqrt(5))/2
np.array([ int(np.floor(n*phi) + np.floor(n*phi**2)) for n in range(1,118) ])

# OUTPUT:
array([  3,   8,  11,  16,  21,  24,  29,  32,  37,  42,  45,  50,  55,
        58,  63,  66,  71,  76,  79,  84,  87,  92,  97, 100, 105, 110,
       113, 118, 121, 126, 131, 134, 139, 144, 147, 152, 155, 160, 165,
       168, 173, 176, 181, 186, 189, 194, 199, 202, 207, 210, 215, 220,
       223, 228, 231, 236, 241, 244, 249, 254, 257, 262, 265, 270, 275,
       278, 283, 288, 291, 296, 299, 304, 309, 312, 317, 320, 325, 330,
       333, 338, 343, 346, 351, 354, 359, 364, 367, 372, 377, 380, 385,
       388, 393, 398, 401, 406, 409, 414, 419, 422, 427, 432, 435, 440,
       443, 448, 453, 456, 461, 464, 469, 474, 477, 482, 487, 490, 495])

If you’re interested in the original proof of this formula (or a fascinating glimpse into what written mathematics looked like in the early 20th century), here is Wythoff’s original paper from 1907.

Delirious ducks

This week’s Riddler Classic is about random walks on a lattice:

Two delirious ducks are having a difficult time finding each other in their pond. The pond happens to contain a 3×3 grid of rocks.

Every minute, each duck randomly swims, independently of the other duck, from one rock to a neighboring rock in the 3×3 grid — up, down, left or right, but not diagonally. So if a duck is at the middle rock, it will next swim to one of the four side rocks with probability 1/4. From a side rock, it will swim to one of the two adjacent corner rocks or back to the middle rock, each with probability 1/3. And from a corner rock, it will swim to one of the two adjacent side rocks with probability 1/2.

If the ducks both start at the middle rock, then on average, how long will it take until they’re at the same rock again? (Of course, there’s a 1/4 chance that they’ll swim in the same direction after the first minute, in which case it would only take one minute for them to be at the same rock again. But it could take much longer, if they happen to keep missing each other.)

Extra credit: What if there are three or more ducks? If they all start in the middle rock, on average, how long will it take until they are all at the same rock again?

Here is my solution:
[Show Solution]

One duck

We will solve this problem by starting with a simpler one and building up to the full problem. Let’s consider a single duck on the grid of rocks. At each step, the duck moves randomly to one of the other nearby rocks. We can describe this process using a Markov chain. Each state is a different rock, and the directed edges are labeled with the transition probabilities. Here is a diagram of what the Markov chain looks like:

If we number the states in the same way as above, then the transition matrix for this Markov chain is:
\[
P = \begin{bmatrix}
0 & \tfrac12 & 0 & \tfrac12 & 0 & 0 & 0 & 0 & 0 \\
\tfrac13 & 0 & \tfrac13 & 0 & \tfrac13 & 0 & 0 & 0 & 0 \\
0 & \tfrac12 & 0 & 0 & 0 & \tfrac12 & 0 & 0 & 0 \\
\tfrac13 & 0 & 0 & 0 & \tfrac13 & 0 & \tfrac13 & 0 & 0 \\
0 & \tfrac14 & 0 & \tfrac14 & 0 & \tfrac14 & 0 & \tfrac14 & 0 \\
0 & 0 & \tfrac13 & 0 & \tfrac13 & 0 & 0 & 0 & \tfrac13 \\
0 & 0 & 0 & \tfrac12 & 0 & 0 & 0 & \tfrac12 & 0 \\
0 & 0 & 0 & 0 & \tfrac13 & 0 & \tfrac13 & 0 & \tfrac13 \\
0 & 0 & 0 & 0 & 0 & \tfrac12 & 0 & \tfrac12 & 0 \\
\end{bmatrix}
\]The way to interpret this matrix is that $P_{ij}$ is the probability that we will transition from $i$ to $j$. This explains why the rows sum to 1 (the matrix is right-stochastic). We can express this fact neatly using matrix multiplication as: $P \mathbf{1} = \mathbf{0}$, where $\mathbf{0}$ and $\mathbf{1}$ are column vectors of all-zeros and all-ones, repectively. Mathematically, we can propagate probability distributions through this Markov chain via matrix multiplication. If we have some initial distribution $\mathbf{a}$ across the states, for example
\[
\mathbf{a}^\mathsf{T} = \begin{bmatrix} \tfrac12 & \tfrac12 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}
\]i.e. we are equally likely to be in states 1 or 2. Then at the next step, the probability distribution will transition $\mathbf{a}\to\mathbf{b}$ with
\[
\mathbf{b}^\mathsf{T} = \mathbf{a}^\mathsf{T}P = \begin{bmatrix} \tfrac16 & \tfrac14 & \tfrac16 & \tfrac14 & \tfrac16 & 0 & 0 & 0 & 0 \end{bmatrix}
\]So in the next step, we might be in states 1, 2, 3, 4, 5 with the associated probabilities above.

It’s important to distinguish states e.g. $\{1,2,\dots,9\}$ from distributions over states, which are the vectors of probabilities such as the $\mathbf{a}$ and $\mathbf{b}$ used above. In the case where a distribution across states is degenerate, i.e. concentrated entirely on a particular state $s$, then I’ll use the following notation to denote the corresponding distribution over states:
\[
\left( \mathbf{e}_s\right)_i = \begin{cases} 1 & \text{if } i=s \\ 0 & \text{otherwise} \end{cases}
\]So, for example, $\mathbf{e}_2^\mathsf{T} = \begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}$.

Stopping time

We are interested in the expected hitting time, which is the number of steps it will take on average to travel from some initial state $s$ to some target state $t \in \mathcal{T}$. To keep this general, I’m assuming there can be multiple target states, indicated by the set $\mathcal{T}$. Let’s define $\mathbf{q}_s$ to be the expected hitting time to get from state $s$ to any state in our terminal set $\mathcal{T}$. It turns out that $\mathbf{q}$ satisfies the recursive relationship:
\[
\mathbf{q}_s = \begin{cases}
0 & \text{if }s \in \mathcal{T} \\
1 + \sum_{i} P_{si} \mathbf{q}_i & \text{otherwise}
\end{cases}
\]The first case is clear: if we start in the terminal set, then we have already arrived, so the hitting time is zero. If we are outside the terminal set, then the expected hitting time will be $1$ plus the weighted sum of the hitting times of wherever we end up after the next transition. The recursion above is essentially the Bellman equation from dynamic programming.

Define the vector $\mathbf{t}_i = \begin{cases} 0 & \text{if }i \in \mathcal{T}\\ 1 & \text{otherwise} \end{cases}$.
We can rewrite the equation above in concise vector form as:
\[
\mathbf{q} = \textrm{diag}(\mathbf{t}) \left( \mathbf{1} + P \mathbf{q} \right)
\]Using the fact that $\textrm{diag}(\mathbf{t})\mathbf{1} = \mathbf{t}$, we can further simplify and obtain:
\[
\left( I-\textrm{diag}(\mathbf{t}) P \right) \mathbf{q} = \mathbf{t}
\]For the case of one duck, if we set the terminal set to be $\mathcal{T} = \{5\}$, then we substitute $\mathbf{t} = \begin{bmatrix} 1 & 1 & 1 & 1 & 0 & 1 & 1 & 1 & 1 \end{bmatrix}^\mathsf{T}$ above and obtain:
\[
\begin{bmatrix}
1 & -\tfrac12 & 0 & -\tfrac12 & 0 & 0 & 0 & 0 & 0 \\
-\tfrac13 & 1 & -\tfrac13 & 0 & -\tfrac13 & 0 & 0 & 0 & 0 \\
0 & -\tfrac12 & 1 & 0 & 0 & -\tfrac12 & 0 & 0 & 0 \\
-\tfrac13 & 0 & 0 & 1 & -\tfrac13 & 0 & -\tfrac13 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & -\tfrac13 & 0 & -\tfrac13 & 1 & 0 & 0 & -\tfrac13 \\
0 & 0 & 0 & -\tfrac12 & 0 & 0 & 1 & -\tfrac12 & 0 \\
0 & 0 & 0 & 0 & -\tfrac13 & 0 & -\tfrac13 & 1 & -\tfrac13 \\
0 & 0 & 0 & 0 & 0 & -\tfrac12 & 0 & -\tfrac12 & 1 \\
\end{bmatrix}
\begin{bmatrix}\mathbf{q}_1\vphantom{\tfrac12} \\ \mathbf{q}_2\vphantom{\tfrac12} \\ \mathbf{q}_3\vphantom{\tfrac12} \\
\mathbf{q}_4\vphantom{\tfrac12} \\ \mathbf{q}_5\vphantom{\tfrac12} \\ \mathbf{q}_6\vphantom{\tfrac12} \\
\mathbf{q}_7\vphantom{\tfrac12} \\ \mathbf{q}_8\vphantom{\tfrac12} \\ \mathbf{q}_9\vphantom{\tfrac12} \end{bmatrix}
= \begin{bmatrix}
1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\
0\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12}
\end{bmatrix}
\]Inverting this matrix and solving for $\mathbf{q}$, here is the result in array form:
\[
\begin{array}{|c|c|c|} \hline
\mathbf{q}_1 = 6 & \mathbf{q}_2 = 5 & \mathbf{q}_3 = 6 \\ \hline
\mathbf{q}_4 = 5 & \mathbf{q}_5 = 0 & \mathbf{q}_6 = 5 \\ \hline
\mathbf{q}_7 = 6 & \mathbf{q}_8 = 5 & \mathbf{q}_9 = 6 \\ \hline
\end{array}
\]So if we start at node 1, it will take us on average 6 moves to reach node 5. It is by no means obvious that the expected hitting times should be integers, since they represent an average path length over an infinite number of possible paths. They just turn out to be integers in this case. The symmetry observed here also makes sense, and we could leverage it to simplify the problem to a Markov chain with only 3 states, but I’ll leave that as an exercise to the reader! (the reason I didn’t do this from the start is that I want to generalize easily to the multi-duck case, and the symmetries there are more complicated).

What about the variant where we start at the origin but we want to know how much time it will take us to return to the origin? In this case, we simply start counting from our second move, and we add 1 to the answer to account for the move we skipped. If we start at node 5, i.e. our initial distribution is $\mathbf{e}_5^\mathsf{T}$, then our distribution on the second turn is $\mathbf{e}_5^\mathsf{T} P$, and so the expected number of moves to start at 5 and return to 5 is: $\mathbf{e}_5^\mathsf{T} P \mathbf{q} + 1 = 6$.

Two ducks

At first glance, the two-duck version might seem much more challenging than the on-duck version. As we shall see, the problem is essentially the same once we look at it the right way. The key is to imagine a Markov chain where instead of the states being $\{1,2,\dots,9\}$, the states are $\{(1,1), (1,2), \dots, (9,9)\}$. In other words, there are 81 states, consisting of all possible pairs $(s_1,s_2)$, where $s_i \in \{1,2,\dots,9\}$ is the position of duck $i$.

What would be the transition matrix for such a Markov chain? The transition probability $(a_1, a_2) \to (b_1, b_2)$ is simply $P_{a_1,b_1} P_{a_2,b_2}$, i.e. the product of the respective transition probabilities for each duck. This means that if we order our states according to the first duck first (lexicographical ordering), then the transition matrix for this augmented Markov chain will be $P\otimes P$, where $\otimes$ is the Kronecker product.

So we can solve the problem in a very similar fashion. This time, our terminal set consists of all the nodes $(s_1, s_2)$ where $s_1=s_2$. There are 9 such nodes, and we can form the associated $\mathbf{t}$ by evaluating $\mathbf{t} = \mathbf{1}-\textrm{vec}(I)$, where $\mathrm{vec}$ is the vectorized version of the matrix, which you obtain by enumerating the columns. For example:
\[
\mathrm{vec} \begin{bmatrix}a_{11} & a_{12} & a_{13} \\
a_{21} & a_{22} & a_{23} \\
a_{31} & a_{32} & a_{33}
\end{bmatrix} =
\begin{bmatrix}
a_{11} \\ a_{21} \\ a_{31} \\ a_{12} \\ a_{22} \\ a_{32} \\ a_{13} \\ a_{23} \\ a_{33} \end{bmatrix}
\]the new equation for $\mathcal{q}$ is $\bigl( I-\textrm{diag}(\mathbf{t}) (P\otimes P) \bigr) \mathbf{q} = \mathbf{t}$ and the solution to the problem where you start at $(5,5)$ and want to know the expect time it takes for us to reach the terminal set (where the two ducks are on the same rock) can be computed following the same logic as in the one-duck case. That is, $\mathbf{e}_{(5,5)}^\mathsf{T} (P\otimes P) \mathbf{q} + 1$.

A cautionary note

While it may seem like this is straightforward algebra, there is a major issue I neglected to discuss: invertibility. When we solve the equation for $\mathbf{q}$ in the one-duck case, the equation looks like $A\mathbf{q} = \mathbf{t}$, and the solution is $\mathbf{q} = A^{-1}\mathbf{t}$. But what if $A$ isn’t invertible? This corresponds to cases where the Markov chain is not connected. So if you start on one island and the terminal set is on another island, you’ll never get there; the stopping time is infinite. That doesn’t happen with one duck, but it does happen with two ducks! This is because every time a duck moves to a new stone, the parity of the stone changes from odd to even or vice versa. So if one duck started on rock 1 and the other started on rock 2, they would never meet.

From a linear algebra standpoint, it just means that even though the matrix is not invertible, if we restrict our attention to the “island” (read: subspace) where the two ducks have the same parity, then the associated transition matrix for that island will be invertible and everything will be fine. From a practical standpoint, the equation above can still be solved; we just need to replace the inverse by a pseudoinverse, and ignore all the components of $\mathbf{q}$ for which the two components have different parity. More on this later…

Computation

I computed the result using the following Matlab code:

% transition matrix for one duck
P = [ 0 1/2 0 1/2 0  0  0  0  0
     1/3 0 1/3 0 1/3 0  0  0  0
      0 1/2 0  0  0 1/2 0  0  0
     1/3 0  0  0 1/3 0 1/3 0  0
      0 1/4 0 1/4 0 1/4 0 1/4 0
      0  0 1/3 0 1/3 0  0  0 1/3
      0  0  0 1/2 0  0  0 1/2 0
      0  0  0  0 1/3 0 1/3 0 1/3
      0  0  0  0  0 1/2 0 1/2 0 ];

T = kron(P,P);     % transition matrix
t = 1-vec(eye(9)); % terminal states

% starting distribution
s = kron([0 0 0 0 1 0 0 0 0],[0 0 0 0 1 0 0 0 0])';

% compute stopping time as a rational number
stop_time = s'*T*( (eye(81)-diag(t)*T)\t ) + 1
rats(stop_time)

The resulting expected stopping time was $\frac{363}{74}$, or approximately $4.905$ steps. We can also find an approximate answer via simulation. Here is what you get by simulating one million trials:

The staggered decrease of the probabilities is a real effect; it’s not a result of approximation error!

Many ducks

If we have $n$ ducks, we can clearly just generalize the approach used above. This time, the transition matrix will be $\underbrace{P\otimes \cdots \otimes P}_{n\text{ times}}$. The issue here is that our transition matrix will be quite large: $9^n \times 9^n$, to be precise. One way to shrink the number of states is to return to the idea of parity. The reason our transition matrix is large is because it describes all transitions between all possible duck configurations. Since we know all ducks start on the same node and that parity will be preserved, we only need to worry about the subset of states for which all ducks have the same parity. Namely, there will be $5^n$ odd states and $4^n$ even states. So we can restrict ourselves to a smaller $(5^n+4^n)\times(5^n+4^n)$ transition matrix. This is still large, but it’s an improvement. We will also leverage the fact that the transition matrix will be sparse, which will help with computation.

After all this work, we end up with the following result (for up to 6 ducks):

The values are as follows:

Number of ducks	Expected rendez-vous time
2	4.9054
3	18.4360
4	66.7420
5	237.3955
6	825.3364

It’s clear the expected rendez-vous time grows exponentially, but I haven’t found a practical way to approximate it or bound it. Unfortunately, it’s pretty difficult to extend my approach beyond $n=6$ ducks. To give you an idea of scale, our naive transition matrix would have had $9^6=531,441$ rows and columns. After our reduction procedure, we’re down to $5^6+4^6=19,721$. The transition matrix is sparse (about 98.5% of its entries are zeros), but that still leaves about 6 million nonzero entries, which makes computing $\mathbf{q}$ a challenge.

Further reductions are possible if we leverage symmetries. Let’s say we have 2 ducks. If we don’t do any reductions, there are $9^2=81$ states. If we reduce using even/odd parity, we drop to $5^2+4^2 = 41$.

From a stopping time standpoint, the only thing that matters is how many ducks are on each stone; it doesn’t matter which ducks are on which stones. This immediately allows us to reduce $5^n+4^n$ to ${n+4 \choose 4} + {n+3\choose 3}$, which is a dramatic improvement. In the two-duck case, this brings us to $25$.
Deeper symmetries occur as well. For example, a pair of ducks on 1 and 3 is equivalent to a pair of ducks on 7 and 9 by rotational symmetry. In the two-duck case, this brings us to $8$ states. Counting and bookkeeping these states might require some group theory…

Ultimately, these reductions would have a dramatic effect, reducing the scaling from exponential to polynomial. I didn’t implement any of these additional reductions because at this point (up to $n=6$) the exponential trend was clear. I didn’t think that extending the trend further would be particularly enlightening.

Mismatched socks

This week’s Riddler Classic is a problem familiar to many…

I have $n$ pairs of socks in a drawer. Each pair is distinct from another and consists of two matching socks. Alas, I’m negligent when it comes to folding my laundry, and so the socks are not folded into pairs. This morning, fumbling around in the dark, I pull the socks out of the drawer, randomly and one at a time, until I have a matching pair of socks among the ones I’ve removed from the drawer.

On average, how many socks will I pull out of the drawer in order to get my first matching pair?

Here is my solution:
[Show Solution]

Computing the expectation

Define $T$ to be number of turns that the game lasts (this is a random variable). Define $q_k = \textbf{P}(T > k)$ to be the probability that the game has still not ended after $k$ draws. We can compute this explicitly. There are $\binom{2n}{k}$ possible (equally likely) ways of choosing $k$ socks among the $2n$ socks in the drawer. If there is no match, we must have drawn $k$ socks of different colors. There are $\binom{n}{k}$ ways to choose the $k$ colors and $2^k$ ways to pick the socks once the colors are chosen. Putting everything together and performing some algebraic manipulations, we obtain:
\[
q_k = \frac{\binom{n}{k} 2^k}{\binom{2n}{k}}
= \frac{n! \cdot (2n-k)! \cdot 2^k}{(2n)! \cdot (n-k)!}
= \frac{ \binom{2n-k}{n} 2^k }{\binom{2n}{n}}
\]Let’s also define $p_k = \textbf{P}(T=k)$ to be the probability that the game ends on precisely the $k^\text{th}$ draw. Note that $q_{k} = p_{k+1}+p_{k+2}+\cdots+p_{n+1}$, or alternatively, $p_k = q_{k-1}-q_k$. From the definition of expected value:
\[
\textbf{E}(T) = \sum_{k=1}^{n+1} k\,p_k = \sum_{k=1}^{n+1}k(q_{k-1}-q_k) = \sum_{k=0}^n q_k
\]This is a manifestation of the general fact that $\textbf{E}(T) = \sum_{k\ge 0}\textbf{P}(T > k)$. Make the change $k\mapsto n-k$ in the expectation:
\[
\textbf{E}(T) = \frac{\sum_{k=0}^n \binom{2n-k}{n} 2^k}{\binom{2n}{n}}
= \frac{2^n \sum_{k=0}^n \binom{n+k}{k} 2^{-k}}{\binom{2n}{n}}
\]Consider the sum in the numerator. Define:
\[
S_n = \sum_{k=0}^n \binom{n+k}{k} 2^{-k}
\]We can evaluate this sum recursively:
\begin{align}
S_{n+1} &= \sum_{k=0}^{n+1} \binom{n+k+1}{k} 2^{-k} \\
&= \sum_{k=0}^{n+1}\left[ \binom{n+k}{k} + \binom{n+k}{k-1} \right] 2^{-k} \\
&= \sum_{k=0}^{n+1} \binom{n+k}{k}2^{-k} + \sum_{k=0}^{n+1} \binom{n+k}{k-1} 2^{-k} \\
&= S_n + \binom{2n+1}{n+1}2^{-n-1} + \sum_{k=0}^{n} \binom{n+k+1}{k} 2^{-k-1} \\
&= S_n + \tfrac{1}{2}S_{n+1}
\end{align}Solving, we obtain $S_{n+1} = 2S_n$, and therefore conclude that $S_n = 2^n$. Substituting this result back into the formula we had above, we obtain:

$\displaystyle
\textbf{E}(T) = \frac{4^n}{\binom{2n}{n}}
$

To get a feel for how this expression varies with $n$, we can use Stirling’s approximation, which leads us to the approximation:
\[
\textbf{E}(T) \approx \sqrt{\pi n}
\]So the expected number of socks we’ll have to draw grows as the square root of the number of pairs. For example, when we have $n=10$ pairs of socks, $\textbf{E}(T) = \frac{262144}{46189}\approx 5.67546$. Meanwhile, the approximation yields $\sqrt{10\pi} \approx 5.60499$. The approximation gets better as $n$ gets larger. Here is a plot showing the true expected value and the approximation:

Computing the distribution

We have the expectation $\textbf{E}(T) = \frac{4^n}{\binom{2n}{n}} \approx \sqrt{\pi n}$. But what about the distribution of $T$ itself? Does it approach anything interesting as $n$ gets large? First, we need to compute the exact probability mass function, which is nothing other than $p_k$. We can do this from the recursion we derived earlier: $p_k = q_{k-1}-q_k$. This simplifies to:
\[
p_k = \frac{k-1}{n-k+1} \frac{\binom{2n-k}{n} 2^{k-1} }{\binom{2n}{n}}\qquad\text{for }k=1,2,\dots,n+1
\]We can plot these distributions along with their means to see what it looks like. Here are the distributions of $T$ for $n=\{10,20,50,100\}$.

To understand what this distribution looks like as $n$ gets very large, we have to do a bit of asymptotic analysis. This part is a little rough around the edges (read: not particularly rigorous), but it seems to give a reasonable answer. My reasoning was that as $n\to\infty$, $k$ is small compared to $n$ because the mean grows like $\sqrt{n}$. Therefore, we will assume $k \sim \sqrt{n}$.

Using Mathematica, we can look at how $\frac{p_k}{k-1}$ varies for small $k$. Using a series approximation about $k=0$ and $n=\infty$, we find that:
\[
\log\left( \frac{p_k}{k-1} \right) \approx \log(\tfrac{1}{2n}) + O(\tfrac{1}{n}) + \tfrac{5}{4n} k-\tfrac{1}{4n}k^2+O(k^3)
\]Therefore, we conclude that:
\[
p_k \approx \frac{k}{2n} e^{-\tfrac{k^2}{4n}}
\]This means that as $n\to\infty$, the distribution of stopping time $T$ approaches a Rayleigh distribution with parameter $\sigma = \sqrt{2n}$. This also gives the correct mean for large $n$, as the mean of the Rayleigh distribution is $\sigma\sqrt{\frac{\pi}{2}} = \sqrt{\pi n}$. For a comparison, here is the plot for $n=500$ along with the Rayleigh approximation:

Prismatic puzzle

This week’s Riddler Classic is simple to state, but tricky to figure out:

What whole number dimensions can rectangular prisms have so that their volume (in cubic units) is the same as their surface area (in square units)?

Here is my solution:
[Show Solution]

If the side lengths of the prism are $a,b,c$, then the volume being equal to the surface area implies that
\[
abc = 2ab+2ac+2bc.
\]Let’s assume the side lengths are ordered, so $a \ge b \ge c$. We therefore have: $ab \ge ac \ge bc$. So $abc = 2ab+2ac+2bc \leq 6ab$. Simplifying, we obtain $c \leq 6$. We can now consider each case separately:

If $c=6$, our equation reduces to $ab=3a+3b$. Rearranging and factoring, this is $(a-3)(b-3)=9$. Since $c=6$ was assumed to be the smallest side length, each of the factors $(a-3)$ and $(b-3)$ is at least $3$. So the only solution is $(6,6,6)$.
If $c=5$, our equation reduces to $3ab=10a+10b \le 20a$. Therefore, $3b \le 20$ and so $b \le 6$. Since $c\le b$, this means we only have to check $b=5$ and $b=6$. Checking each case, we find the only solution is $(10,5,5)$.
If $c=4$, our equation reduces to $ab=4a+4b$. Rearranging and factoring, this is $(a-4)(b-4)=16$. Each way of factoring $16$ as a product of two factors yields a different solution. These are: $16\times 1$, $8\times 2$, and $4\times 4$. The resulting $(a,b,c)$ triples are: $(20,5,4)$, $(12,6,4)$, $(8,8,4)$.
If $c=3$, our equation reduces to $ab=6a+6b$. Rearranging and factoring, this is $(a-6)(b-6)=36$. Each way of factoring $36$ as a product of two factors yields a different solution. These are: $36\times 1$, $18\times 2$, $12\times 3$, $9\times 4$, $6\times 6$. The resulting triples are: $(42,7,3)$, $(24,8,3)$, $(18,9,3)$, $(15,10,3)$, $(12,12,3)$.
If $c=2$, our equation reduces to $0=4a+4b$, which has no solutions.
If $c=1$, our equation reduces to $0=ab+2a+2b$, which also has no solutions.

And that’s it! Overall, the problem has exactly 10 solutions. They are:
$(6,6,6)$, $(10,5,5)$, $(20,5,4)$, $(12,6,4)$, $(8,8,4)$, $(42,7,3)$, $(24,8,3)$, $(18,9,3)$, $(15,10,3)$, $(12,12,3)$.

Note: I assumed that “whole number” in this context means a positive integer. If we allow for side lengths of zero, then any degenerate prism with side lengths $(a,0,0)$ will have both area and volume equal to zero, so there are infinitely many solutions.

Another note: If there were some way to upper-bound $a$, then we could do an exhaustive numerical search to find all solutions because there would only be finitely many triples $(a,b,c)$ to search. Alas, I have not been able to find an upper bound on $a$. I suspect it might not be possible, since if we set $a\to\infty$, the problem reduces to $bc=2b+2c$, which is the 2D version of the problem (a rectangle whose area equals its perimeter)!

Optimal HORSE

This week’s Riddler Classic is about how to optimally play HORSE — the playground shot-making game. Here is the problem.

Two players have taken to the basketball court for a friendly game of HORSE. The game is played according to its typical playground rules, but here’s how it works, if you’ve never had the pleasure: Alice goes first, taking a shot from wherever she’d like. If the shot goes in, Bob is obligated to try to make the exact same shot. If he misses, he gets the letter H and it’s Alice’s turn to shoot again from wherever she’d like. If he makes the first shot, he doesn’t get a letter but it’s again Alice’s turn to shoot from wherever she’d like. If Alice misses her first shot, she doesn’t get a letter but Bob gets to select any shot he’d like, in an effort to obligate Alice. Every missed obligated shot earns the player another letter in the sequence H-O-R-S-E, and the first player to spell HORSE loses.

Now, Alice and Bob are equally good shooters, and they are both perfectly aware of their skills. That is, they can each select fine-tuned shots such that they have any specific chance they’d like of going in. They could choose to take a 99 percent layup, for example, or a 50 percent midrange jumper, or a 2 percent half-court bomb.

If Alice and Bob are both perfect strategists, what type of shot should Alice take to begin the game?

What types of shots should each player take at each state of the game — a given set of letters and a given player’s turn?

Here is how I solved the problem:
[Show Solution]

Let’s denote the score numerically to make notation a bit easier. A score of $k$ means the player has accrued $k$ letters. First player to reach a score of 5 loses. We’ll consider the game from the point of view of Alice and define $V(a,b)$ to be the probability that Alice will win if it’s currently her turn, she has a score of $a$, and Bob has a score of $b$. This assumes both players are playing optimally.

Suppose it’s Alice’s turn, the score is $(a,b)$, and Alice attempts a shot that succeeds with probability $p$. One of three things can happen:

Alice and Bob both make it. This happens with probability $p^2$. The score remains $(a,b)$, it’s still Alice’s turn, and she wins with probability $V(a,b)$.
Alice makes the shot but Bob misses. This happens with probability $p(1-p)$. The score is now $(a,b+1)$, and it’s still Alice’s turn. So Alice wins with probability $V(a,b+1)$.
Alice misses her shot. This happens with probability $1-p$. In this case, it’s now Bob’s turn and the score is still $(a,b)$. Since both players are playing optimally, by symmetry Bob will use the strategy Alice would use if the score had been $(b,a)$ and therefore win with probability $V(b,a)$. So Alice wins with probability $1-V(b,a)$.

We can therefore write the following recursion:
\[
V(a,b) = p^2 V(a,b) + p(1-p) V(a,b+1) + (1-p) \bigl( 1-V(b,a) \bigr)
\]Alice will choose $p$ in order to maximize $V(a,b)$, her probability of winning. We’ll now work backwards to find $V(a,b)$ for all values of $a$ and $b$.

Computing $V(4,4)$

Setting $a=b=4$, defining $x=V(4,4)$, and using the fact that $V(4,5)=1$ (Alice wins if Bob reaches a score of 5), the recursion becomes:
\[
x = p^2 x + p(1-p) + (1-p)( 1-x)
\]Rearranging, this becomes:
\[
(1-p)(2+p)x = (1-p)(1+p)
\]One solution is to pick $p=1$, which makes the equation vacuously true and we learn nothing about $x$. This corresponds to the case where Alice chooses a shot that goes in 100% of the time. In such a case, Alice would keep shooting forever and the game would never end. So Alice would never lose, but she would also never win. To avoid this case (and to make this solution more realistic), let’s assume $0 \le p \le 1-\epsilon$, where $\epsilon\gt 0$ is a small positive number that characterizes how far from perfect the “easiest” shot is. This leaves us with:
\[
x = \frac{1+p}{2+p}
\]Alice wants to choose $p$ to maximize $x$, her probability of winning. The function above is an increasing function of $p$, and its maximum value is therefore $x=\frac{2-\epsilon}{3-\epsilon}$, which occurs when $p=1-\epsilon$. In short, Alice should shoot her easiest shot. In the limit, her probability of winning will be $\frac{2}{3}$.

Computing $V(a,b)$ in general

Rearranging the recursion and dividing by $p-1$, we obtain:
\[
V(a,b) = \frac{p}{p+1} V(a,b+1) + \frac{1}{p+1} (1-V(b,a))
\]Let’s consider the cases $(a,b)=(3,4)$ and $(a,b)=(4,3)$ where the decisions made are $p$ and $q$ respectively. This leads to the pair of equations:
\begin{align}
V(3,4) &= \frac{p}{p+1} + \frac{1}{p+1} (1-V(4,3)) \\
V(4,3) &= \frac{q}{q+1}V(4,4) + \frac{1}{q+1}(1-V(3,4))
\end{align}Solving these equations, we obtain:
\begin{align}
V(3,4) &= 1-\frac{q}{p q+p+q} V(4,4) \\
V(4,3) &= \left(1-\frac{p}{p q+p+q}\right) V(4,4)
\end{align}Recall the goal is to pick $p$ to maximize $V(3,4)$ and pick $q$ to maximize $V(4,3)$. This happens when $p$ and $q$ are as large as possible. So setting $p=q=1-\epsilon$, we obtain:
\[
V(3,4)=\frac{\epsilon ^2-5 \epsilon +7}{(\epsilon -3)^2},\quad\text{and}\quad
V(4,3)=\frac{(\epsilon -2)^2}{(\epsilon -3)^2}
\]Or, when $\epsilon \to 0$, $V(3,4)=\frac{7}{9}$ and $V(4,3)=\frac{4}{9}$. Ultimately, we can keep repeating this process of solving for pairs of values at a time, until we fill out the whole table. In every case, the same thing happens: Alice should always play as well as she can. The easiest way to solve these equations is to write them out with $p=1-\epsilon$:
\[
V(a,b) = \frac{1-\epsilon}{2-\epsilon} V(a,b+1) + \frac{1}{2-\epsilon}(1-V(b,a))
\]and these are simultaneous linear equations for all choices of $a$ and $b$.

and here is the solution:
[Show Solution]

The optimal strategy is for both Alice and Bob to attempt the surest shots possible no matter what the score is. Using this strategy, if it’s Alice’s turn, her chances of winning are:
\[
\begin{array}{c|cccccc}
& 0 & 1 & 2 & 3 & 4 & 5 \\ \hline
0 & \frac{10802}{19683} & \frac{4241}{6561} & \frac{545}{729} & \frac{617}{729} & \frac{227}{243} & 1 \\
1 & \frac{2984}{6561} & \frac{1216}{2187} & \frac{487}{729} & \frac{191}{243} & \frac{73}{81} & 1 \\
2 & \frac{256}{729} & \frac{328}{729} & \frac{46}{81} & \frac{19}{27} & \frac{23}{27} & 1 \\
3 & \frac{176}{729} & \frac{80}{243} & \frac{4}{9} & \frac{16}{27} & \frac{7}{9} & 1 \\
4 & \frac{32}{243} & \frac{16}{81} & \frac{8}{27} & \frac{4}{9} & \frac{2}{3} & 1 \\
\end{array}
\]Here, the rows correspond to Alice’s current score, and the columns correspond to Bob’s current score. If Alice and Bob can’t play perfectly, and the highest percentage shot they can make goes in with probability $1-\epsilon$, then the probability of winning as a function of $\epsilon$ is given by this table.

What about at the beginning of the game? Assuming Alice goes first and she can attempt shots that go in with probability $0 \le p \le 1-\epsilon$, her chances of winning as a function of $\epsilon$ are:
\[
P(\text{Alice wins}) = \tfrac{(\epsilon -2) \left(\epsilon ^8-20 \epsilon ^7+194 \epsilon ^6-1130 \epsilon ^5+4221 \epsilon ^4-10245 \epsilon ^3+15704 \epsilon ^2-13870 \epsilon +5401\right)}{(\epsilon -3)^9}
\]Here is what the plot looks like as a function of $(1-\epsilon)$.

So, for example, say your easiest layup shot goes in 85% of the time. Then you would substitute $\epsilon = 0.15$ into the formula above (or read it off the plot, and find that if you go first, your chances of winning against an equally matched opponent are about 54.3%. The optimal strategy is for both players to attempt the 85% shot every time.

In the limit of perfect shooting, when players can shoot a shot that goes in 99.999…% of the time, the player that goes first will win on average $\tfrac{10802}{19683}$ of the time, or about 54.88% of the time. Note that playing this way will lead to very long games!

I did not address the version of HORSE with $N$ players, as there are many regional variants of this game with different rules and I wasn’t sure which to look at. It would be interesting to see if the optimal strategy is still for all players to play their best shot all the time!