Pill splitting

Problems of this sort may seem tedious at first glance, but there is a simple and systematic way to approach them. We can encode the pills in the bottle using polynomials. Specifically, we’ll write $x^n y^m$ when there are $n$ full pills and $m$ half-pills. Since there is an element of randomness with how the pills are picked, we will use a weighted sum of such polynomials to describe a mixture or possible states. For example, we could write:
\[
\tfrac{2}{3} x y + \tfrac{1}{3} y^3
\]to describe a scenario where there is a $\frac{2}{3}$ probability that the bottle contains one full pill and one half-pill, and a $\frac{1}{3}$ probability that the bottle contains three half-pills. Note that the coefficients sum to $1$, so exactly one of the two scenarios must occur.

So what happens when we take our daily dose? We can consider each scenario in turn.

• With probability $\frac{\binom{n}{2}}{\binom{n+m}{2}}$, we will pick two full pills. In this case, we cut one in half and return half to the bottle. This means the bottle now contains $x^{n-2}y^{m+1}$.
• With probability $\frac{nm}{\binom{n+m}{2}}$, we will pick one full pill and one half-pill. In this case, we have our exact dose and the bottle now contains $x^{n-1}y^{m-1}$.
• With probability $\frac{\binom{m}{2}}{\binom{n+m}{2}}$, we will pick two half-pills. In this case, we must go back into the bottle and pick one more pill. There are two sub-cases here:
  • With probability $\frac{n}{n+m-2}$, we pick one of the full pills. So we split it in half, and return half to the bottle. The jar has now lost a total of one full pill and one half-pill, so it contains $x^{n-1}y^{m-1}$.
  • With probability $\frac{m-2}{n+m-2}$, we pick one of the half-pills. So we have our exact dose and the bottle now contains $x^n y^{m-3}$.

Combining all the possible cases, we find that every day, the polynomial describing the contents of the bottle change according to the mapping
\begin{multline}
x^n y^m \mapsto \tfrac{n (n-1)}{(m+n) (m+n-1)}x^{n-2} y^{m+1} + \tfrac{m n (3 m+2 n-5)}{(m+n) (m+n-1) (m+n-2)} x^{n-1} y^{m-1} \\
+ \tfrac{m (m-1) (m-2)}{(m+n) (m+n-1) (m+n-2)} x^n y^{m-3}
\end{multline}
There are two nice things about this formulation. First, edge cases are taken care of automatically. For example, if the bottle contains $n=2$ full pills and $m=2$ half-pills, the third term ($x^n y^{m-3}$) in the expansion is impossible, since the bottle does not contain three half-pills. But the associated coefficient is zero, so this term cancels automatically. Similarly, if $m=0$ (the bottle only contains full pills), two of the coefficients vanish and the third becomes $1$, so we obtain $x^n \mapsto x^{n-2} y$.

Second, when the bottle contains a mixture of possible states (described by a polynomial whose coefficients sum to 1), we can apply the above transformation to each term of the polynomial separately, and the final result will again be a mixture of possible states. This follows because the sum of the coefficients of the mapping above sum to 1 as well.

To compute the mixture of states after 10 days, I wrote a Mathematica script that starts at $x^{15}$ (15 full pills) and applies the above transformation recursively. Here is the Mathematica script:

P = x^15; (* initial polynomial *)
Print[P];
For[iter = 0, iter < 9, iter++,
  Q = 0;
  ML = MonomialList[P];
  For[i = 1, i <= Length[ML], i++,
    term = ML[[i]];
    coef = term /. {x -> 1, y -> 1};
    n = Exponent[ term, x];
    m = Exponent[term, y];
    Q = Q + coef ((m(m-1)(m-2))/((n+m)(n+m-1)(n+m-2)) x^n y^(m-3)
          + (n m(2n+3m-5))/((n+m)(n+m-1)(n+m-2)) x^(n-1) y^(m-1)
          + (n(n-1))/((n+m)(n+m-1)) x^(n-2) y^(m+1));
   ];
   P = Q // Expand;
   Print[P];
]

and here is the output from the script:
\begin{gather}
x^{15} \\
x^{13} y\\
\tfrac{1}{7}x^{12}+\tfrac{6 }{7}x^{11} y^2\\
\tfrac{36 }{91}x^{10} y+\tfrac{55 }{91}x^9 y^3\\
\tfrac{23 }{308}x^9+\tfrac{2385}{4004} x^8 y^2+\tfrac{30}{91} x^7 y^4\\
\tfrac{4 }{13}x^7 y+\tfrac{81}{143} x^6 y^3+\tfrac{18}{143} x^5 y^5\\
\tfrac{335}{4004} x^6+\tfrac{87}{154} x^5 y^2+\tfrac{185}{572} x^4 y^4+\tfrac{4}{143} x^3 y^6\\
\tfrac{22573}{56056} x^4 y+\tfrac{42605}{84084} x^3 y^3+\tfrac{101}{1144} x^2 y^5+\tfrac{1}{429}x y^7\\
\tfrac{44783}{240240} x^3+\tfrac{362603}{560560} x^2 y^2+\tfrac{27185}{168168} x y^4 + \tfrac{61}{12012} y^6\\
\tfrac{80529}{101920} x y+\tfrac{21391}{101920} y^3\\
\end{gather}

Therefore, on the $10^\text{th}$ day, there is a $\tfrac{21391}{101920} \approx 20.988\%$ chance the bottle contains three half-pills, and a $\tfrac{80529}{101920} \approx 79.012\%$ chance the bottle contains one full pill and one half-pill. This is the probability that the correct dosage was taken, which is what is asked in the original problem statement.

Large-pill limit

We saw that on the final day, the probability that the bottle contains one full pill and one half-pill is about $79\%$. A natural question to ask is: what happens to this probability as we vary the number of starting pills? If we start with only 3 pills, then we always have one full pill and one half-pill on the next day. So the probability is $100\%$. Here is a plot of what happens as we vary the number of initial pills:

It appears that as we increase the number of initial pills, the probability of having one full pill and one half-pill on the last day decreases. I don’t have a great explanation for why this happens; if you made it this far and have any ideas, write a comment; I would love to hear it!

Generating functions?

I am well aware that the polynomials I used look a lot like generating functions. I looked for a nice way to express the pill-taking transformation as a differential operator acting on the generating function, but I was thwarted by the denominators. The numerators are easy, e.g.
\[
\tfrac{\mathrm{d}^3}{\mathrm{d}y^3}(x^n y^m) = m(m-1)(m-2) x^n y^{m-3}
\]but I couldn’t find a nice way to deal with the denominators. For the case above, the denominator should be $(n+m)(n+m-1)(n+m-2)$. If anybody found a way to solve this using the full generating function machinery, I would love to hear about it!