Riddler - Book Proofs

Quarter flipping on a spinning table

This week’s Riddler Classic is a simple-sounding problem about statistics. But it’s not so simple! Here is a paraphrased version of the problem:

There is a square table with a quarter on each corner. The table is behind a curtain and thus out of your view. Your goal is to get all of the quarters to be heads up — if at any time all of the quarters are heads up, you will immediately be told and win.

The only way you can affect the quarters is to tell the person behind the curtain to flip over as many quarters as you would like and in the corners you specify. (For example, “Flip over the top left quarter and bottom right quarter,” or, “Flip over all of the quarters.”) Flipping over a quarter will always change it from heads to tails or tails to heads. However, after each command, the table is spun randomly to a new orientation (that you don’t know), and you must give another instruction before it is spun again.

Can you find a series of steps that guarantees you will have all of the quarters heads up in a finite number of moves?

Here is my solution:
[Show Solution]

Every time a move is made, the table spins around randomly. This sort of randomization actually simplifies the game considerably by reducing the number of distinct positions and moves possible:

Although there are $2^4 = 16$ different ways to arrange quarters in combinations of Heads (H) and Tails (T) on the four corners of the table, if one position can be reached from another via rotation, then those two positions can be treated as one, since there is no way to distinguish them after the table spin. There are in fact only 6 distinct positions:
- all Heads: (H,H,H,H). This is the winning position
- all Tails: (T,T,T,T).
- opposite Heads: (H,T,H,T). There are 2 configurations for this position.
- one Head: (H,T,T,T). There are 4 configurations for this position.
- one Tail: (T,H,H,H). There are 4 configurations for this position.
- two at a time: (H,H,T,T). There are 4 configurations for this position.
Although there are also $2^4-1=15$ different moves we can make (minus one because we exclude the move where you flip nothing), again the spinning makes several of these moves equivalent. There are really only five different possible moves:
- flip one quarter (ONE)
- flip two opposite quarters (OPP)
- flip two adjacent quarters (ADJ)
- flip three quarters (TRI)
- flip all the quarters (ALL)

Given this information, we can look at all the possible transitions between states that can occur if we make a certain move from any given position. Here is a diagram showing all possible transitions. I excluded the “flip three quarters” move because this move turns out not to be necessary.

For example, if we “flip two opposite quarters”, the (H,H,T,T) position always remains unchanged. Meanwhile, the (H,T,H,T) could lead us to the winning positions (H,H,H,H), or it could lead us to the (T,T,T,T) position.

Winning strategy

A key observation: the “one-off” positions (H,T,T,T) and (T,H,H,H), which I colored in purple, will only ever transition between each other if we flip two quarters or all quarters. We can use this to our advantage: we will begin by only using these moves to try to win. If we cannot win, then we conclude that we must be in one of the “one-off” positions. Essentially, the strategy consists of making moves that successively test possible states. Every time we test a state, if we are correct we reach (H,H,H,H) and the game stops. Otherwise, the game continues and we can eliminate that state from the list of possible states. We continue in this fashion until there are no states left. Here are the logical steps:

Are we in (T,T,T,T)? Try (ALL) to find out.
If the game continues, we must not be in (T,T,T,T).
Are we in (H,T,H,T)? Do (OPP) to escape.
If the game continues, we moved to (T,T,T,T).
Repeat step 1 to test if this is true.
If the game continues, we must not be in (H,T,H,T).
Are we in (H,H,T,T)? Do (ADJ) to escape.
If the game continues, we moved to (T,T,T,T) or (H,T,H,T).
Repeat steps 1-2 to test if this is true.
If the game continues, we must not be in (H,H,T,T).
By elimination, we must be in (T,H,H,H) or (H,T,T,T). Do (ONE) to escape.
If the game continues, we moved to (T,T,T,T), (H,T,H,T), or (H,H,T,T).
Repeat steps 1-3 to test if this is true.

If we “unroll” the recipe above, the sequence guaranteed to win is:
ALL,OPP,ALL,ADJ,ALL,OPP,ALL,ONE,ALL,OPP,ALL,ADJ,ALL,OPP,ALL

These 15 moves are guaranteed to find (H,H,H,H) at some point along the way, no matter what position we start in, and no matter how unlucky we get with the spins of the table.

Tank counting

This week’s Riddler Classic is a simple-sounding problem about statistics. But it’s not so simple! Here is a paraphrased version of the problem:

There are $n$ tanks, labeled $\{1,2,\dots,n\}$. We are presented with a sample of $k$ tanks, chosen uniformly at random, and we observe that the smallest label among the sample is $22$ and the largest label is $114$. Both $n$ and $k$ are unknown. What is our best estimate of $n$?

Here is a short introduction on parameter estimation, for the curious:
[Show Solution]

This problem falls under the purview of parametric statistics, and it’s worth doing a short introduction for those of you not familiar, because this comes up a lot in statistics. The basic setup is that there is a “population” of values that follows some assumed probability distribution that depends on a set of fixed parameters. The goal is to estimate the parameters of this distribution by observing a random “sample” from the population.

For example, our population could be the heights of all people living in the US. Let’s say this population follows a normal distribution, but we don’t know the mean or variance parameters of the distribution. We’d like to estimate the mean, because we want to know the average height of all people in the US. It’s impractical to measure everybody in the US, so instead we survey 100 people at random and obtain their heights. We might then compute the mean of our sample and use this as an estimate for the mean of the population. This is an example of an estimator; a function that maps samples to a parameter estimate.

What are important properties of an estimator? What makes a “good” estimator? If we fix the size of the sample, but take lots of different random samples, we will obtain a set of sample means. They might all be different due to the variability in our samples. Two important concepts are:

The bias of an estimator: this is the difference between the true population mean and the mean of all our sample means. Ideally, we want bias to be as small as possible. For our example, the sample mean estimator has zero bias, so we say that it is unbiased.
The variance of an estimator: this is the sample variance of our set of sample means. Roughly, this is how much variability we can expect our parameter estimates to have. Clearly, we would also like the variance to be as small as possible.

The notions of bias and variance are analogous to the notions of accuracy and precision. Would you rather be accurate or precise? It turns out it’s generally impossible to have both. In parametric statistics, this is referred to as the bias-variance trade-off. We can’t make both bias and variance zero. We can only make one small at the expense of the other. Continuing our average height example, we could have used a much simpler estimator: always estimate six feet as the average height! This estimator has no variance at all because it always estimates the same thing, but it is surely biased, because the true population mean likely isn’t six feet.

Although this was a rather simple example, I hope it illustrates that it doesn’t really make sense to talk about the “best estimator” unless we are specific about our notion of “best”. Sometimes, variance can be reduced dramatically and it only costs us a little bit of bias. This is one of the main benefits of using regularization in machine learning. It’s always a trade-off, and the notion of “best” depends on the context and the application.

With that out of the way, let’s get back to counting tanks!

Here is my solution to the problem:
[Show Solution]

The parameters of our distribution are $n$ (the number of tanks) and $k$, the size of the group. The random variables we get to observe are the minimum and maximum labels, which I’ll call $x$ and $y$, respectively. The probability mass function is
\[
P_{n,k}(x,y) = \frac{y-x-1 \choose k-2}{n \choose k}
\]because there are $y-x-1 \choose k-2$ ways of picking as set of $k$ distinct integers whose minimum and maximum values are $x$ and $y$ respectively, and there are $n\choose k$ total ways of picking $k$ distinct integers out of a set of $n$. We can verify that this is a legitimate probability mass function by checking that
\[
\sum_{x=1}^{n-k+1} \sum_{y=x+k-1}^{n} P_{n,k}(x,y) = 1
\]The reason for the strange bounds in the sums is to account for the fact that not all minimum and maximum values are admissible. For example, if we pick $k=3$ different numbers among $\{1,2,3,4,5\}$, the only possible values for the minimum $x$ are $\{1,2,3\}$.

For this problem, our “sample” is the single pair $(x,y)$ that we observe, and our task is to estimate $n$. Unlike our population height example from above, it’s no longer clear what sort of estimator we should use here. I’ll present some options.

Maximum likelihood estimator

One possible way to estimate $n$ is to use the maximize likelihood estimator (MLE). The likelihood is what you get when you fix $(x,y)$ to be what you observed, and you view $P_{n,k}(x,y)$ as a function of $(n,k)$. This tells you which values of the parameters were likelier given the data you observed. For any fixed $k,x,y$, it’s clear that we can maximize $P_{n,k}(x,y)$ by making $n$ as small as possible, so the MLE is to choose $\hat n = y$. What can we say about bias and variance for the MLE? We can directly compute these quantities using the $P_{n,k}(x,y)$ above:
\[
\text{MLE:}\,(\hat n = y)\, \left\{
\begin{aligned}
\text{Bias} &= -\frac{n-k}{k+1} \\
\text{Variance} &= \frac{k (n+1) (n-k)}{(k+1)^2 (k+2)}
\end{aligned}
\right.
\]As we can see, the estimator has bias. This makes sense: our estimate is the largest number we saw, so we will always be under-estimating the true value of $n$.

Minimum-variance unbiased estimator

Another way to estimate $n$ is to restrict our attention to unbiased estimators. Since we only get to observe $x$ and $y$, a reasonable thing to try is to compute the expectations $\mathbb{E} x$ and $\mathbb{E} y$. Doing so, we obtain:
\[
\mathbb{E} x = \frac{n+1}{k+1},\qquad
\mathbb{E} y = \frac{k(n+1)}{k+1}
\]We can eliminate $k$ by adding these together! Indeed, we find that:
\[
n = \mathbb{E}( x+y-1 )
\]This tells us that if we use $\hat n = x+y-1$ as our estimator, it will be unbiased! We can compute bias and variance as before, and we obtain:
\[
\text{MVUE:}\,(\hat n = x+y-1)\, \left\{
\begin{aligned}
\text{Bias} &= 0 \\
\text{Variance} &= \frac{2 (n+1) (n-k)}{(k+1) (k+2)}
\end{aligned}
\right.
\]Using some more advanced statistics, it’s actually possible to prove that among all unbiased estimators, this is the one with the smallest variance. Not bad! But note that this estimator still has about twice the variance of MLE estimator. Another manifestation of the bias-variance trade-off.

Solution

The notion of “best estimate” is an ambiguous one, so we have to be clear about what we mean by “best”. The maximum-likelihood estimator (MLE) is to estimate $\hat n = y$, the largest label observed. For the example in the problem, we would pick $\hat n = 114$. Although this is the likeliest $n$, it is a biased estimate: it will always under-estimate the true $n$.

Another option is to look for an estimator that is unbiased. That is, if we had a large number of independent samples of groups of tanks and we averaged the estimates across groups, we would on average obtain the true $n$. The minimum-variance unbiased estimator (MVUE) is to pick $\hat n = x+y-1$ (the sum of largest and smallest labels observed, minus one). For the example in the problem, this yields $\hat n = 135$. While the MVUE is unbiased, it has about twice the variance as the MLE, which means it produces more highly variable estimates.

This higher variability is unavoidable; if we want to have an unbiased estimator, we must be willing to give up some variance! This is analogous to the notion of precision vs accuracy. The MVUE above is actually the estimator with the smallest variance among all unbiased estimators, so if you’re looking for an unbiased estimator, this is the “best” one!

Simulation

To demonstrate that the MVUE estimator works as intended, I created the following simulation: I set $n=100$ and $k=15$, and then picked $k$ tanks at random among the $n$, recorded the minimum and maximum labels $x$ and $y$, and computed the MVUE estimator for $n$, which is $\hat n = x+y-1$. I repeated the process $10^6$ times, and recorded all $\hat n$ values. Based on the calculations above, the mean and standard deviation of this set of $\hat n$ estimates should be:
\[
\mathbb{E}(\hat n) = n = 100,\quad\text{and}\quad
\mathrm{Var}(\hat n) = \sqrt{\tfrac{2(n+1)(n-k)}{(k+1)(k+2)}} = 7.9451
\]The values I obtained were $100.0015$ and $7.9279$, respectively. Looks like a match! Here is a histogram of the results of the simulation:

So it looks like the MVUE estimator performs as expected. Note that even if $k$ is completely unknown to us or if every random sample of tanks we observe has a different size $k$, the mean will still be $n$. Changing $k$ only affects our confidence in our estimates. A larger $k$ means we can be more confident (lower variance).

Thanos snaps

This week’s Riddler Express problem is a nod to the recently released finale to the Avengers saga. No spoilers!

Thanos, the all-powerful supervillain, can snap his fingers and destroy half of all the beings in the universe.

But what if there were 63 Thanoses, each snapping his fingers one after the other? Out of 7.5 billion people on Earth, how many can we expect would survive?

If there were N Thanoses, what would the survival fraction be?

Here is my solution:
[Show Solution]

First things first: If there are $S$ Thanos snaps and the probability of survival after one snap is $p$, then the expected fraction of survivors is $p^S$. In the case where the population is $7.5$ billion and the probability of survival is $p=0.5$, if there are $S=63$ snaps, we can calculate the expected number of survivors to be:
\[
\frac{7.5\times 10^9}{2^{63}} = 8.13\times 10^{-10}
\]This number is so small that there will almost certainly be no survivors.

A more interesting interpretation: But what if Thanos were not immune to the snap? After all, when Thanos snaps his fingers, half of all the beings in the universe are destroyed! This should include Thanos! So when the first Thanos snaps his fingers, on average half of the remaining Thanoses will be destroyed before having the opportunity to snap their fingers.

This means that $S$, the total number of snaps, is a random variable. We can compute its probability mass function recursively. Specifically, let’s define $P_N(s)$ to be the probability that there will be $S=s$ snaps given that there are initially $N$ Thanoses. Let’s start with some easy cases.

If there are no Thanoses, there will be no snaps. So $P_0(0) = 1$.
If there is at least one Thanos, there will be at least one snap. So $P_N(0) = 0$ for $N=1,2,\dots$.

The recursion for larger values of $s$ is given by:
\begin{align}
P_N(s) &= \sum_{k=0}^{N-1} \mathrm{Prob}(k\text{ Thanoses survive first snap}) P_k(s-1) \\
&= \sum_{k=0}^{N-1} \binom{N-1}{k} p^{k}(1-p)^{N-1-k} P_{k}(s-1)
\end{align}This follows because if there are $s$ snaps remaining, then during the first one, the remaining $N-1$ Thanoses are at risk of destruction. The number of surviving Thanoses follows a binomial distribution with $N-1$ trials and probability $p$, which explains why the probability that $k$ survive the first snap is equal to $\binom{N-1}{k} p^{k}(1-p)^{N-1-k}$. This recursion can be evaluated numerically. Here is what the distribution of snaps looks like for $p=0.5$ and $N=63$:

So there will most likely be $5$ or $6$ snaps. Once we know how many snaps there are, the effect on the population of the universe is the same as in the first interpretation. So the expected fraction of survivors in the universe will be $\mathbb{E}(p^S)$, or:
\[
\mathbb{E}(\text{fraction of survivors if we have }N\text{ Thanoses}) = \sum_{s=0}^N p^s P_{N}(s)
\]Here is a plot of the surviving fraction as a function of $N$ when $p=0.5$:

When there are 63 Thanoses, the expected survival rate is about 2.22%. So with an initial population of 7.5 billion, we will be left with, on average, about 166.5 million.

If you’re curious about the shape of this graph, we can approximate it like this: With $N$ Thanoses and $p=0.5$, there are roughly $\log_2(N+1)$ flips that occur. If exactly that many flips occur, the fraction of survivors will be roughly $(0.5)^{\log_2(N+1)} = \tfrac{1}{N+1}$. This tells us that as $N$ gets large, the fraction of survivors will scale roughly like $\frac{1}{N}$, where $N$ is the number of Thanoses. This is a crude approximation, but it captures the correct asymptotic behavior.

If you’re interested in the (Python) code I wrote to generate this solution (and the figures), you can find it on github here.

Gift card puzzle

Here is a puzzle from the Riddler about gift cards:

You’ve won two gift cards, each loaded with 50 free drinks from your favorite coffee shop. The cards look identical, and because you’re not one for record-keeping, you randomly pick one of the cards to pay with each time you get a drink. One day, the clerk tells you that he can’t accept the card you presented to him because it doesn’t have any drink credits left on it.

What is the probability that the other card still has free drinks on it? How many free drinks can you expect are still available?

Here is my solution:
[Show Solution]

Let’s suppose each card starts with $n$ drinks. There are three ways the sequence of buying drinks can end:

The first card gets maxed out and the second card has $k$ drinks remaining with $k\ge 1$. This means we purchased a total of $2n-k$ drinks and $n$ of them ended up on the first card. Then, we tried to purchase one additional drink on the first card and the buying stopped. The probability of this occurring is: $(\tfrac{1}{2})^{2n-k+1}\binom{2n-k}{n}$.
both cards get maxed out and then we try to buy one more drink on either card. This means we purchased a total of $2n$ drinks and $n$ of them ended up on the first card. The probability of this occurring is $(\tfrac{1}{2})^{2n}\binom{2n}{n}$.
The second card gets maxed out and the first card has $k$ drinks remaining with $k\ge 1$. This is exactly analogous to the first case, and the probability of this occurring is again: $(\tfrac{1}{2})^{2n-k+1}\binom{2n-k}{n}$.

Putting all of this together, we end up with a tidy formula for $p(n,k)$, the probability that when the game ends, the other card has exactly $k$ drinks remaining on it:

$\displaystyle
p(n,k) = \frac{1}{2^{2n-k}} \binom{2n-k}{n}
$

To double-check that this is a valid probability mass function, we should have $\sum_{k=0}^n p(n,k) = 1$. This is indeed the case, but it’s rather challenging to verify. Here is a link to the WolframAlpha computation. Here is plot of the probability mass function for $n=50$.

We can now answer the first question: the probability that the other card still has drinks on it is one minus the probability that there are no drinks left, i.e. $1-p(n,0)$. This evaluates to:

$\displaystyle
\mathrm{Prob}\Bigl(\begin{smallmatrix}\text{there are drinks left}\\\text{on the other card}\end{smallmatrix}\Bigr)
= 1-\frac{1}{4^n} \binom{2n}{n} \approx 1-\frac{1}{\sqrt{\pi n}}
$

This makes sense; the only way there can be no drinks left is if we purchased $2n$ drinks, and we were lucky to use each card exactly $n$ times. We would expect this probability to get smaller as $n$ gets larger. The approximation I used above comes from Stirling-like approximations for binomial coefficients (see here). When $n=50$ as in the problem statement, the probability evaluates to $92.04\%$. Here is a plot showing how the probability of there being drinks on the other card slowly tends to $1$ as $n$ gets larger:

Now let’s address the second question: what is the expected number of drinks on the other card? We will compute the expected value directly from the probability mass function: $E_n = \sum_{k=0}^n k\, p(n,k)$. We find:

$\displaystyle
\mathbb{E}\Bigl(\begin{smallmatrix}\text{number of drinks left}\\\text{on the other card}\end{smallmatrix}\Bigr)
= \frac{2n+1}{4^n}\binom{2n}{n}-1 \approx \frac{2n+1}{\sqrt{\pi n}}-1
$

Again, the derivation is challenging so I omitted it, but here is the WolframAlpha link. Here is a plot of this function as it changes with $n$ along with the approximation (which is quite good!)

When $n=50$ as in the problem statement, the expected number of drinks remaining on the other card is $7.0385$. This corresponds to the expected value of the probability distribution plotted in the first image.

Elf music

This holiday-themed Riddler problem is about probability:

In Santa’s workshop, elves make toys during a shift each day. On the overhead radio, Christmas music plays, with a program randomly selecting songs from a large playlist.

During any given shift, the elves hear 100 songs. A cranky elf named Cranky has taken to throwing snowballs at everyone if he hears the same song twice. This has happened during about half of the shifts. One day, a mathematically inclined elf named Mathy tires of Cranky’s sodden outbursts. So Mathy decides to use what he knows to figure out how large Santa’s playlist actually is.

Help Mathy out: How large is Santa’s playlist?

Here is my solution:
[Show Solution]

Settlers in a circle

In this Riddler problem, the goal is to spread out settlements in a circle so that they are as far apart as possible:

Antisocial settlers are building houses on a prairie that’s a perfect circle with a radius of 1 mile. Each settler wants to live as far apart from his or her nearest neighbor as possible. To accomplish that, the settlers will overcome their antisocial behavior and work together so that the average distance between each settler and his or her nearest neighbor is as large as possible.

At first, there were slated to be seven settlers. Arranging that was easy enough: One will build his house in the center of the circle, while the other six will form a regular hexagon along its circumference. Every settler will be exactly 1 mile from his nearest neighbor, so the average distance is 1 mile.

However, at the last minute, one settler cancels his move to the prairie altogether (he’s really antisocial). That leaves six settlers. Does that mean the settlers can live further away from each other than they would have if there were seven settlers? Where will the six settlers ultimately build their houses, and what’s the maximum average distance between nearest neighbors?

Here is my solution:
[Show Solution]

I approached this problem from a modeling and optimization perspective. If there are $N$ settlers and we imagine the coordinates of the settlers are $(x_i,y_i)$ for $i=1,\dots,N$ and $d_i$ is the distance between the $i^\text{th}$ settler and its nearest neighbor, then we can model the problem as follows:
\[
\begin{aligned}\underset{x,y,d\,\in\,\mathbb{R}^N}{\text{maximize}} \quad & \frac{1}{N}\sum_{i=1}^N d_i \\
\text{such that} \quad & x_i^2 + y_i^2 \le 1 &&\text{for }i=1,\dots,N\\
& (x_i-x_j)^2 + (y_i-y_j)^2 \ge d_i^2 &&\text{for }i,j=1,\dots,N\text{ and }j\ne i
\end{aligned}
\]The objective is to minimize the average distance between each settler and its nearest neighbor (the average of the $d_i$). The first constraint says that each settler must be within the circle (a distance of at most $1$ from the origin). The second constraint says that the $i^\text{th}$ settler is a distance at least $d_i$ from each of the other settlers.

The game is then to find $(x_i,y_i,d_i)$ that satisfy these constraints and maximize the objective. Broadly speaking, this is an example of a mathematical optimization problem. Unfortunately, the problem as stated is nonconvex problem because of the form of the second constraint. This means that the problem may have local maximizers that are not globally optimal. There are two ways forward:

First, we could use local search: start with randomly placed settlers, wiggle their positions in ways such that the objective continues to increase, and then stop once we can no longer improve the objective. Examples of such approaches include hill-climbing and interior-point methods. In general, such approaches only find local maxima, so we should try many random starting positions to ensure we find the best possible answer.

Second, we could use global search: these are methods that attempt to find a globally optimal solution by considering many configurations simultaneously and adding random perturbations that can “kick” a solution out of a local maximum. Examples include simulated annealing and particle swarms. These approaches tend to be much slower than local search, but they have a much better chance of finding a global optimum.

Ultimately, I used local search with several random initializations. Here is how the maximum average distance scales with the number of settlements.

The bar plot also shows how well we would do if we distributed the settlements evenly on the circumference (“all on the border”) and if we put one in the middle and the rest evenly distributed on the circumference (“one in the center”). We can also examine what the settlement distributions actually look like. Here they are:

The case $N=6$ is particularly interesting because it has a settlement that is almost in the center, but not quite. Indeed, the average minimum distance to neighbors for this case is $1.002292$, so we benefit ever so slightly by placing one settlement off-center.

We can solve for the $N=6$ case analytically as well, but the solution isn’t pretty (fair warning!) It turns out the settlements are distributed throughout the circle in the following way:

Here, $x$ is the distance between the origin and the point $A$. The other distances satisfy $z \lt y \lt x+1$, and the average distance to the nearest neighbor is therefore $\tfrac{1}{6}(1+x+2y+3z)$. Note that once we set $x$, this uniquely determines the positions of all the other points. After some messy trigonometry, we obtain the following equations relating $x$, $y$, $z$ and two auxiliary variables $p$ and $q$:
\begin{align}
y^2&=1+x-x^2-x^3\\
z^2&=1-x^2+2x\left(pq-\sqrt{1-p^2}\sqrt{1-q^2}\right) \\
p&= \tfrac{1}{2}(1-2x-x^2)\\
q&= \tfrac{1}{2}(1-x+x^2+x^3)
\end{align}We can eliminate $\{y,z,p,q\}$ and solve for the average distance $\bar d = \frac{1+x+2y+3z}{6}$ as a function of $x$ alone:
\[
\bar d = \frac{x+1}{12} \left(2+4 \sqrt{1-x}+3 \sqrt{2}\sqrt{1-x} \sqrt{\left(x^3+2 x^2-x+2\right)-x \sqrt{x+3} \sqrt{x^3+x^2-x+3}} \right)
\]We can then plot this function of $x$, and we get the following curve:

Interestingly, this function reaches its maximum just a bit after $x=0$. Zooming in, we can see more clearly:

The maximum occurs at about $x=0.0555108$ and the maximum value is $1.00229$, just as we found before. If we try to solve for this analytically by taking the derivative of this function, setting it equal to zero, and eliminating rationals, we obtain the following result… That the optimal $x$ a root of the following $30^\text{th}$ order polynomial:
\begin{align}
x^{30}&+\tfrac{152}{9} x^{29}+\tfrac{9259}{81} x^{28}+\tfrac{6851}{18} x^{27}+\tfrac{383363 }{648}x^{26}+\tfrac{1999495} {5832}x^{25}\\
&+\tfrac{43478263}{52488} x^{24}+\tfrac{75726373}{23328} x^{23}+\tfrac{3282369305}{1679616} x^{22}-\tfrac{31642768891}{7558272} x^{21}\\
&+\tfrac{367775655235}{136048896} x^{20}+\tfrac{392027905801}{34012224} x^{19}-\tfrac{2016793647847}{136048896} x^{18}-\tfrac{1258157703385}{68024448} x^{17}\\
&+\tfrac{4545130566235}{136048896} x^{16}-\tfrac{22059424877}{17006112} x^{15}-\tfrac{8448216227365}{136048896} x^{14}+\tfrac{2878062707167}{68024448} x^{13}\\
&+\tfrac{7398046956073}{136048896} x^{12}-\tfrac{312655708903}{3779136} x^{11}+\tfrac{17266856539}{1679616} x^{10}+\tfrac{63670102441}{839808} x^9\\
&-\tfrac{29827845749}{559872} x^8-\tfrac{541401565}{23328} x^7+\tfrac{619158287}{23328} x^6+\tfrac{6168305}{2592} x^5\\
&-\tfrac{156285}{32} x^4-\tfrac{5153}{36} x^3+\tfrac{3923}{12} x^2+\tfrac{45}{2} x-\tfrac{35}{16}
\end{align}Yuck! Here is what you get if you plot this polynomial; as expected, there is a root at $x=0.0555108$.

Alice and Bob fall in love

In this interesting Riddler problem, we’re dealing with a possibly unbounded sequence of… children? Here it goes:

As you may know, having one child, let alone many, is a lot of work. But Alice and Bob realized children require less of their parents’ time as they grow older. They figured out that the work involved in having a child equals one divided by the age of the child in years. (Yes, that means the work is infinite for a child right after they are born. That may be true.)

Anyhow, since having a new child is a lot of work, Alice and Bob don’t want to have another child until the total work required by all their other children is 1 or less. Suppose they have their first child at time T=0. When T=1, their only child is turns 1, so the work involved is 1, and so they have their second child. After roughly another 1.61 years, their children are roughly 1.61 and 2.61, the work required has dropped back down to 1, and so they have their third child. And so on.

(Feel free to ignore twins, deaths, the real-world inability to decide exactly when you have a child, and so on.)

Five questions: Does it make sense for Alice and Bob to have an infinite number of children? Does the time between kids increase as they have more and more kids? What can we say about when they have their Nth child — can we predict it with a formula? Does the size of their brood over time show asymptotic behavior? If so, what are its bounds?

Here is an explanation of my derivation:
[Show Solution]

Let’s call $t_k$ the time at which Alice and Bob have their $k^\text{th}$ child. The first child is born at time zero, so $t_1 = 0$. Since the work due to a child is one divided by their age, the work due to child $k$ at time $t$ is given by:
\[
\text{work at time }t = \begin{cases}
\frac{1}{t-t_k} & \text{if }t > t_k \\
0 & \text{otherwise}
\end{cases}
\]We also know that each new child is born once the total work due to all previous children equals $1$. Therefore, we have the following equations:
\begin{align}
t_2\text{ satsifies:} && \frac{1}{t_2-t_1} &= 1 \\
t_3\text{ satisfies:} && \frac{1}{t_3-t_2} + \frac{1}{t_3-t_1} &= 1 \\
t_4\text{ satisfies:} && \frac{1}{t_4-t_3} + \frac{1}{t_4-t_2} + \frac{1}{t_4-t_1} &= 1 \\
\cdots\qquad \\
t_k\text{ satisfies:} && \sum_{i=1}^{k-1} \frac{1}{t_k-t_i} &= 1
\end{align}and so on. We also know that $t_1 \lt t_2 \lt t_3 \lt \ldots$ since the children are born sequentially. The nice thing about these equation is that they can be solved sequentially. We know that $t_1=0$. Substituting this into the first equation yields $t_2=1$. Substituting these values into the second equation yields:
\[
\frac{1}{t_3-1}+\frac{1}{t_3} = 1
\]Rearranging yields $t_3^2-3t_3+1=0$. This equation has two solutions, but we take the one that satisfies $t_3 \gt t_2$, and we obtain $t_3=\tfrac{3+\sqrt{5}}{2} \approx 2.618$. We can continue in this fashion, solving each subsequent equation by substituting the previous values… but things get complicated in a hurry. The next equation becomes:
\[
t_4^3-\left(\tfrac{11+\sqrt{5}}{2}\right)t_4^2+\left(\tfrac{13+3\sqrt{5}}{2}\right)t_4-\left(\tfrac{3+\sqrt{5}}{2}\right)=0
\]Once again, we find that there is a unique solution that satisfies $t_4 \gt t_3$. This solution is approximately $t_4\approx 4.5993$. There doesn’t appear to be any automatic way to continue this process. Every time we try to compute the next $t_k$, we have to find the root of an even more complicated polynomial!

We’ll take a numerical approach, but let me first address one issue: how can we be sure that these ever-more-complicated polynomials always have a root that satisfies our constraints? For example, what if the polynomial for $t_4$ had no roots larger than $t_3$? What if it had many such roots? Thankfully, there is a simple way we can guarantee that there is always exactly one admissible root to each polynomial by leveraging the intermediate value theorem. Specifically, the function
\[
f_k(t) = \frac{1}{t-t_1} + \frac{1}{t-t_2} + \cdots + \frac{1}{t-t_{k-1}}
\]with $t \gt t_{k-1} \gt \ldots \gt t_1$ is a continuous and strictly decreasing function. Moreover, $\lim_{t\to {t_{k-1}}^+}f_k(t) = +\infty$ and $\lim_{t\to\infty} f_k(t) = 0$. Therefore, by the intermediate value theorem, there must be some point $t_k \in (t_{k-1},\infty)$ such that $f_k(t_k)=1$. The strict monoticity of $f_k$ ($f_k$ is strictly decreasing) implies that this $t_k$ must be unique.

The root may be unique, but how can we find it efficiently? Here, we can leverage a further property of $f_k$; namely convexity. The fact that $f_k$ is convex (and smooth) means that we can leverage some efficient methods from convex analysis that make use of derivative information. There are many ways to get the answer from here… What I ended up doing was to use the L-BFGS algorithm on the objective $(f_k(t)-1)^2$. My implementation in Julia using the Optim.jl package took about 40 seconds to evaluate birth dates of the first 1,000 children.

Note: there are many ways to formulate this as an optimization problem. One way is to note that solving $f_k(t)=1$ is equivalent to extremizing the antiderivative of $f_k(t)-1$. In other words, we could solve:
\[
t_k = \underset{t \gt t_{k-1}}{\text{arg min}} \left( t-\sum_{i=1}^{k-1}\log(t-t_i) \right)
= \underset{t \gt t_{k-1}}{\text{arg max}}\,\, e^{-t}\prod_{i=1}^{k-1}(t-t_i)
\]Incidentally, this last expression for $t_k$ makes it clear that there should be a unique $t_k$. When $t\gt t_{k-1}$, the product is a polynomial that is strictly increasing, and it’s multiplied by $e^{-t}$, which will drive it to zero eventually. There must therefore be a unique maximizer. I tried implementing both of these approaches and both turned out to be slower than directly minimizing $(f_k(t)-1)^2$.

If you’re just interested in the answers to the questions, here they are:
[Show Solution]

Does it make sense for Alice and Bob to have an infinite number of children? Biological limitations aside, yes! We proved that for each $t_k$, there exists a unique $t_{k+1} \gt t_k$ that satisfies the constraints of the problem. So by induction, Alice and Bob can have as many children as they like.
Does the time between kids increase as they have more and more kids? The answer is yes, and I will provide a “proof by picture”. After finding $t_k$ for $k=1,2,\dots$, we know the birth dates of each child. Then it remains to compute $(t_{k+1}-t_k)$ for each $k$ (the time between births). Here is the plot

The first point on the left says that we wait $1$ year after child $1$ is born, we wait about $1.6$ years after child $2$ is born, and so on. Plotted on a log scale as above, we see the curve is almost a perfect straight line. Here is a formula we can use to approximate the number of years to wait:

$\displaystyle
\text{years to wait after child }k\text{ is born} \,\approx\, 1 + 2.19\cdot \log_{10} k
$

So yes, the time between children increases (without bound), but it does so slowly. By the time Alice and Bob have their 1,000th child, they will be waiting approximately 7.5 years between kids.
What can we say about when they have their $N^\text{th}$ child? Is there a formula we can write down? I couldn’t find an analytic expression, but I did find a pretty good approximation. First, let’s make a plot to see what we’re working with. Here is what the birth years is like for the first 19 children

We can see here that the growth rate is slightly faster than linear. Since the time between children fit so well to a curve of the form $a+b\log x$, it makes sense to fit the birth dates themselves to the integral of this form. In an effort to fit the simplest possible curve, I did a bit of experimenting and found that the form $a x + x \log x$ does a really good job. I found the equation of best fit using linear regression and put the equation in the plot above. If we extend the time horizon, we get a similar result:

This is a pretty good fit; the plot actually depicts the true data with the fit overlaid on top. The fit is so good that you can’t even see the data underneath! So for large $N$, a good approximation is:

$\displaystyle
\text{year when child }N\text{ is born} \,\approx\,-0.315N + N\log N
$
The final question is whether the size of the brood over time shows asymptotic behavior. We already plotted the time of birth vs child number, and found the relationship to be approximately: $y \approx -0.315 N + N\log N$. So here, we need to look at the inverse relationship and ask what $N$ does as a function of $y$. Since $y$ grows a little faster than linearly, this means $N$ will grow a little more slowly than linearly.

Beer pong

This interesting twist on the game of Beer Pong appeared on the Riddler blog. Here it goes:

The balls are numbered 1 through N. There is also a group of N cups, labeled 1 through N, each of which can hold an unlimited number of ping-pong balls. The game is played in rounds. A round is composed of two phases: throwing and pruning.

During the throwing phase, the player takes balls randomly, one at a time, from the infinite supply and tosses them at the cups. The throwing phase is over when every cup contains at least one ping-pong ball. Next comes the pruning phase. During this phase the player goes through all the balls in each cup and removes any ball whose number does not match the containing cup. Every ball drawn has a uniformly random number, every ball lands in a uniformly random cup, and every throw lands in some cup. The game is over when, after a round is completed, there are no empty cups.

How many rounds would you expect to need to play to finish this game? How many balls would you expect to need to draw and throw to finish this game?

Here is my solution:
[Show Solution]

We’ll address the second question first: How many balls would you expect to need to draw and throw to finish the game? For this question, the rounds don’t matter at all. The game will end once there is at least one correct ball in each cup, and the pruning phase at the end of each round only removes incorrect balls so we can effectively neglect it.

Each time a ball is thrown, there is a $1/N$ chance it will land in the correct cup. This is a classical coupon collector problem (see my write-up on the card collection completion problem for more background). The solution rests on the fact that if some trial is successful with probability $p$, then we must perform $1/p$ trials on average before we obtain our first success. The probability of landing our first correct ball is $1/N$ so it will take an average of $N$ tosses to do so. The probability of landing our next new correct ball (in a different cup) is $\tfrac{1}{N}\cdot \tfrac{N-1}{N}$. This takes on average $\tfrac{N^2}{N-1}$ more tosses. Continuing in this fashion and summing up all of the intervals until each cup is filled, the total expected number of tosses is:

$\displaystyle
T_\text{tot} = N^2 \sum_{k=1}^N \frac{1}{k} \approx N^2(\log N + \gamma)
$

where $\gamma \approx 0.5772$ is the Euler-Mascheroni constant.

Now let’s turn our attention to the first question: How many rounds would you expect to play to finish this game? This is a much more complicated question, so to illustrate the approach I’ll solve it for the case $N=3$. We’ll begin by breaking each round into a number of stages, which track the progress of cups as they are filled. In the diagram below, each cup is either empty (no ball), is red (incorrect ball), or is green (correct ball). We can also figure out all the transition probabilities and write them on the arrows.

In this diagram, self-loops have whatever probability is required so that the outgoing arrows sum to $1$. This is an example of a Markov chain. We start on the left (all cups empty), and we work our way to the right (all cups full). Here is what the transition matrix looks like for $N=3$ when we arrange the states from top to bottom and left to right:
\[
A =
\begin{bmatrix}
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
2/3 & 2/9 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
1/3 & 1/9 & 1/3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 4/9 & 0 & 4/9 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 2/9 & 4/9 & 2/9 & 5/9 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 2/9 & 0 & 1/9 & 2/3 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 2/9 & 0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1/9 & 2/9 & 0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 1/9 & 2/9 & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 1/9 & 0 & 0 & 0 & 1
\end{bmatrix}
\begin{array}{c}
(k=0,j=0)\\
(k=1,j=0)\\
(k=1,j=1)\\
(k=2,j=0)\\
(k=2,j=1)\\
(k=2,j=2)\\
(k=3,j=0)\\
(k=3,j=1)\\
(k=3,j=2)\\
(k=3,j=3)\\
\end{array}
\]Note that the columns sum to $1$ since all outgoing probabilities must sum to $1$. The question is: what is the probability that we’ll end up in each of the possible ending point? (the four last states). This is called the limiting distribution. In our case, we can compute it easily. First, we can eliminate the self-loops by normalizing. This works because we don’t actually care how much time is spent at each stage. See below for an illustration of this normalization.

Here is what our transition matrix looks like after normalization
\[
A_\text{norm} =
\begin{bmatrix}
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
2/3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
1/3 & 1/7 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 4/7 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 2/7 & 2/3 & 2/5 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 1/3 & 0 & 1/4 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 2/5 & 0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1/5 & 1/2 & 0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 1/4 & 2/3 & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 1/3 & 0 & 0 & 0 & 1
\end{bmatrix}
\begin{array}{c}
(k=0,j=0)\\
(k=1,j=0)\\
(k=1,j=1)\\
(k=2,j=0)\\
(k=2,j=1)\\
(k=2,j=2)\\
(k=3,j=0)\\
(k=3,j=1)\\
(k=3,j=2)\\
(k=3,j=3)\\
\end{array}
\]Once again, the columns sum to $1$, but this time only the absorbing states have self-loops. Finding the limiting distribution simply amounts to finding the limit $\lim_{t\to\infty} A_\text{norm}^t$. Since the upper-left block of $A_\text{norm}$ is nilpotent (the diagonal is zero), powers of the matrix will eventually be constant. It suffices to evaluate $A_{\text{norm}}^{2N}$. Extracting the last N+1 rows and the columns corresponding to the states $\{1,3,6,\dots,\tfrac{(N+1)(N+2)}{2}\}$ (the states where we have exactly $k$ correct balls in cups, for $k=0,1,\dots,N$), we obtain the transition matrix for one round of the game. Here is the result we obtain in the case $N=3$:
\[
P = \begin{bmatrix}
16/105 & 0 & 0 & 0 \\
41/105 & 1/3 & 0 & 0 \\
5/14 & 1/2 & 2/3 & 0 \\
1/10 & 1/6 & 1/3 & 1
\end{bmatrix}
\begin{array}{c}
(k=0)\\
(k=1)\\
(k=2)\\
(k=3)\\
\end{array}
\]In other words, $P_{ij}$ is the probability that if we currently have $j$ correct balls in cups, we will have $i$ after the next round is over. From here, our next task is to compute the expected number of rounds it will take to travel from $k=0$ to $k=N$. This can be computed as follows: if we call $E_k$ the expected number of rounds starting from $k$, then we have:
\begin{align}
E_0 &= 1 + P_{00} E_0 + P_{10}E_1 + \cdots + P_{N0}E_N \\
E_1 &= 1 + P_{01} E_0 + P_{11}E_1 + \cdots + P_{N1}E_N \\
&\vdots \\
E_{N-1} &= 1 + P_{0,N-1} E_0 + P_{1,N-1}E_1 + \cdots + P_{N,N-1} E_N \\
E_N &= 0
\end{align}Or, in matrix-vector notation: $E = \mathbb{1} – e_N + P^\mathsf{T} E$. This can be solved via simple matrix inversion. We can then extract $E_0$ from the solution, which is the expected number of rounds starting with all empty cups. As far as I can tell, there does not appear to be an easy way to extract an analytic formula, but all the steps above are very easy to compute numerically. Here are plots of the results for small values of $N$:

I also overlayed a result obtained via Monte-Carlo simulation using $1000$ simulated games. This shows that the analytical results agree well with simulation. Here are further numerical results, this time extended up to $N=100$:

It appears that the number of rounds required grows approximately linearly with $N$ while the number of balls required grows quadratically. This latter result is to be expected, since we know the formula in this case is approximately $N^2 \log N$.

Hand sort

A card-rearranging problem on the Riddler blog. Here it goes:

You play so many card games that you’ve developed a very specific organizational obsession. When you’re dealt your hand, you want to organize it such that the cards of a given suit are grouped together and, if possible, such that no suited groups of the same color are adjacent. (Numbers don’t matter to you.) Moreover, when you receive your randomly ordered hand, you want to achieve this organization with a single motion, moving only one adjacent block of cards to some other position in your hand, maintaining the original order of that block and other cards, except for that one move.

Suppose you’re playing pitch, in which a hand has six cards. What are the odds that you can accomplish your obsessive goal? What about for another game, where a hand has N cards, somewhere between 1 and 13?

Here is my solution:
[Show Solution]

I was unable to find an analytic or closed-form expression for the solution to this problem, so I’ll present an old fashioned brute-force solution. The idea is simple:

Enumerate all possible hands
For each hand, enumerate all possible moves that consist of moving some group of adjacent cards to somewhere else in the hand.
If any such move succeeds in sorting the hand, mark the hand as “good”.
Once we’re done, count all “good” hands and divide that by the total number of hands

Simple, right? I’d like to point out is that this is not an approximate solution; there is no simulation involved. I’m counting all the possible scenarios, so this approach produces the exact answer. The only problem is that there are a lot of hands to check, and the bigger the hands get, the longer it takes to check them, since there are more possible moves. Successfully solving this problem requires coding it up in a way that a computer can find the solution in a reasonable amount of time. In order to make this work, I used several tricks to reduce the computation required:

We don’t have to enumerate all possible hands because only the suits matter, but we do have to be careful about how we count. We begin with a full deck of cards (52 distinct cards). Let’s suppose we want to count hands of size 2. There are $52\cdot 51 = 2652$ possible hands. But since numbers don’t matter, we can remove the number and just look at suits. There are then 16 possible hands:
\[
\begin{array}{cccc}
♠♠ & ♠♣ & ♠\color{red}{♥} & ♠\color{red}{♦} \\
♣♠ & ♣♣ & ♣\color{red}{♥} & ♣\color{red}{♦} \\
\color{red}{♥}♠ & \color{red}{♥}♣ & \color{red}{♥}\color{red}{♥} & \color{red}{♥}\color{red}{♦} \\
\color{red}{♦}♠ & \color{red}{♦}♣ & \color{red}{♦}\color{red}{♥} & \color{red}{♦}\color{red}{♦}
\end{array}
\]suppose for example that the hands with different suits but same color are not sortable. There are four such hands: $\{ ♠♣, ♣♠, \color{red}{♥}\color{red}{♦}, \color{red}{♦}\color{red}{♥} \}$. Then we might conclude that $\tfrac{4}{16} = 0.25$ of hands are not sortable. But this would be wrong! (thanks to Guy Moore for pointing this out). The reason is that these 16 hands don’t all occur with equal probability. There are $13\cdot 13 = 169$ ways of picking two cards of different suits, but only $13\cdot 12 = 156$ ways of picking two cards of the same suit! So each type of card should be weighted by its likelihood of occurrence. This means the probability of picking an unsortable hand should really be:
\[
\frac{4\cdot 169}{4\cdot 156 + 12\cdot 169} = \frac{13}{51} \approx 0.2549
\]So we can collapse the large number of possible hands to this more manageable size as long as we compensate by appropriately weighting the different items. In combinatorics, this collapsed group of cards (where we care about the order and the suit, but not the number of the card) is an example of a multiset permutation.

Note: If we don’t weigh the hands according to likelihood (i.e. assume all hands are equally likely), this is equivalent to assuming the deck has infinitely many cards (but the same number of cards in each suit).
Adjacent cards of the same suit can be collapsed to a single card. This is because there is never any benefit to moving a block of cards if it splits an existing contiguous block. So for example:
\[
\{♠,♠,\color{red}{♥},\color{red}{♥},\color{red}{♦}\} \implies \{♠,\color{red}{♥},\color{red}{♦}\}.
\]This greatly reduces the computations because if two different hands collapse to the same reduced hand, we only need to do the work for one of them! In computer programming, this technique of storing the values of computations so that you don’t end up computing the same thing many times is called memoization.
If a hand consists of $N$ cards (after we collapse adjacent suit repetitions as described above), then there are $\binom{N+1}{3}$ possible moves we can make. This is because every move is characterized by three locations. I’ll illustrate this with an example. Say our hand looks like this: $\{♠,\color{red}{♥},\color{red}{♦},♠,\color{red}{♥},♣\}$ and we want to sort it by moving cards 4 and 5 and inserting them between cards 1 and 2. In other words, we want to perform:
\[
\{♠,\color{red}{♥},\color{red}{♦},(♠,\color{red}{♥}),♣\} \implies \{♠,(♠,\color{red}{♥}),\color{red}{♥},\color{red}{♦},♣\}
\implies \{♠,\color{red}{♥},\color{red}{♦},♣\}
\]where we did one final collapse at the end. This transformation can be represented by inserting three “bars” between cards as follows. Then the swap simply consists of exchanging the cards between the two pairs of bars! Here is the diagram:
\[
\{\,♠\,\,|\underbrace{\color{red}{♥}\,\,\color{red}{♦}}_\text{swap this}|\underbrace{♠\,\,\color{red}{♥}}_\text{with this}|\,\,♣\,\}
\implies \{♠, ♠, \color{red}{♥},\color{red}{♥},\color{red}{♦},♣\}
\]Each triplet of bars corresponds to a possible swap, and we can insert bars in $N+1$ possible spots (in between each card and on either end).
If the hand has 8 or more cards in it after we perform the collapse, then it’s impossible to sort it in one move. This is because our final hand can have at most 4 groups of cards in it, and a single move can cause at most 3 collapses. With 8 cards, we can get down to 5 but not 4. Here is an example of a 7-card hand that can be sorted in one move:
\[
\{\,♠\,|\,\color{red}{♥}\,♣\,|\,♠\,\color{red}{♥}\,|\,♣\,\color{red}{♦}\,\}
\implies \{\,♠\,♠\,\color{red}{♥}\,\color{red}{♥}\,♣\,♣\,\color{red}{♦}\,\}
\implies \{\,♠\,\color{red}{♥}\,♣\,\color{red}{♦}\,\}
\]

Numerical results

There was a bit of confusion as to how to interpret the meaning of the requirement “the cards of a given suit are grouped together and, if possible, such that no suited groups of the same color are adjacent”. One way to interpret this statement is that we only require cards of a different suit but same color to be separated if it’s actually possible to separate them. For example, the two-card hand $\{\color{red}{♥}\,\color{red}{♦}\}$ counts as being sorted because there are no spades or clubs that can be placed between the heart and the diamond. In this scenario, here are the results:

Note that hands of up to size $N=4$ are always sortable in one move. The exact solution for the case $N=6$ using this interpretation is $\tfrac{51083}{83895} \approx 60.8892\%$. Another way to interpret the requirement is that cards of a different suit but same color can never be adjacent. This means that hands like $\{\color{red}{♥}\,\color{red}{♦}\}$ are never sortable no matter what you do. In this case, we get the result:

It’s interesting to note that the probability of a “sortable” hand in this interpretation actually reaches a maximum when $N=4$ and then decreases afterwards. The exact solution for the case $N=6$ using this interpretation is $\tfrac{1735996}{2936325} \approx 59.1214\%$.

I wrote my code in Julia, and you can view my notebook here. The computation gets slower in an exponential fashion as we increase $N$. The cases $N \leq 8$ take on the order of seconds, but computation time quadruples every time $N$ increases by 1. The last case, $N=13$, took about 40 minutes! While there is a significant amount of computation to do, there is also significant computational overhead due to using exact arithmetic instead of making approximations. For this, I used Julia’s BigInt datatype. So the solutions I obtained were exact. For example, the probability of a sortable hand for $N=13$ (using the first interpretation) turns out to be:
\[
\tfrac{10954472929065768960}{3954242643911239680000} = \tfrac{30785713171}{11112737293000} \approx 0.27703\%
\]Using approximations the whole way leads to a computation time of less than 10 minutes for the most complicated case.

Tether your goat!

A geometry problem from the Riddler blog. Here it goes:

A farmer owns a circular field with radius R. If he ties up his goat to the fence that runs along the edge of the field, how long does the goat’s tether need to be so that the goat can graze on exactly half of the field, by area?

Here is my solution:
[Show Solution]

There are many ways to solve this problem, so I figured I would pick an unorthodox way: using calculus! Consider the diagram below.

The circular field is centered at $O$. The goat is tethered at $C$ and the tether has length $r$. Let $B$ be the point diametrically opposed to $C$ and let $A$ be the farthest point along the perimeter that the goat can reach. Let $\theta$ be the angle $BOA$, as shown. If $r$ changes a little bit, by some amount $dr$, then $\theta$ will change by some amount $d\theta$. Likewise, the fraction of the total area reachable by the goat will change by $dA$. Let’s calculate how these quantities are related.

First, note that $dA$ is simply the area between the two circular arcs centered at $C$ and passing through $A$ and $A’$, respectively. We must also divide by the total area of the field in order to get a ratio. As $dr\to 0$, we have:
\[
dA = \frac{\tfrac{1}{2} \theta\, d(r^2)}{\pi R^2}
\]Our next task is to find out how $d\theta$ is related to $d(r^2)$. Look at the triangle $OAC$. By the law of cosines, we have:
$r^2 = 2R^2(1 + \cos(\theta))$. Taking differentials, we obtain:
\[
d(r^2) = -2R^2 \sin(\theta)\,d\theta
\]Substituting back into our expression for $dA$, we obtain:
\[
dA = -\tfrac{1}{\pi} \theta\,\sin(\theta)\,d\theta
\]We can now integrate! Note that $A=1$ when $\theta=0$. Therefore, we have:
\begin{align}
A(\theta) &= 1 + \int_{0}^\theta \tfrac{1}{\pi} t\, \sin(t)\,dt \\
&= \frac{\pi + \theta\,\cos(\theta)-\sin(\theta)}{\pi}
\end{align}This tells us the ratio of areas as a function of $\theta$. If we want this as a function of $r$ instead, we can use the law of cosines again to eliminate $\theta$ and write $A$ as a function of $r$. This yields:
\[
A(\rho) = 1+\frac{\left(\frac{\rho ^2}{2}-1\right) \cos ^{-1}\left(\frac{\rho ^2}{2}-1\right)}{\pi }-\frac{\rho \sqrt{4-\rho ^2}}{2\pi}
\]where I made the substitution $\rho := r/R$. We can plot this quantity as a function of $\rho$, and we get:

As we might expect, when $r=0$, $A=0$. And the area ratio grows monotonically with $r$ until we reach $r=2R$, at which point the goat can reach the farthest point (which is point $B$) and therefore can reach the entire field. There is no closed-form expression for the exact tether length that leads to a particular area ratio, since that would require solving the above $A(\rho)$ for $\rho$. However, we can easily solve the equation numerically. Doing so, we obtain:

$\displaystyle
\frac{\text{length of tether}}{\text{radius of field}} \approx 1.15873
$

Finally, here is a diagram of what the field and tether look like when the area ratio is exactly 1/2.