\documentclass[12pt,reqno]{amsart}%{article}
%\usepackage{a4ams}
\usepackage{amsfonts,amsmath,amssymb}
\usepackage{rotating}
\DeclareMathOperator{\dilog}{dilog}
\topmargin 0 pt                        
\textheight 46\baselineskip     
\advance\textheight by \topskip
%\setlength{\parindent}{0pt}            
%\setlength{\parskip}{2pt plus 2pt minus 1pt}     
\setlength{\textwidth}{155mm}         
\setlength{\oddsidemargin}{5.6mm}     
\setlength{\evensidemargin}{5.6mm}   

\def\al{\alpha}

\newtheorem{theorem}{Theorem}

\newtheorem{lemma}{Lemma}

\newtheorem{prop}{Proposition}
\numberwithin{equation}{section}

\def\P{{\Bbb {P}}}
\def\t{\vartheta}\def\tt{\bar{\vartheta}}
\def\v{\textsf{value}}
\def\mm{\widehat m}
\def\ma{\widehat a}
\def\Card{\mathit{Card}}
% other useful commands
\newcommand{\Leaf}{\qedsymbol}

\title[The pills problem revisited]
{The pills problem revisited}
\author[C.~Brennan]{Charlotte A.~C.~Brennan}

\address{C.~Brennan\\
The John Knopfmacher Centre for Applicable Analysis and Number Theory\\
 School of Mathematics\\
University of the Witwatersrand, P.~O. Wits\\ 
2050 Johannesburg, South Africa}
\email{charlotte@maths.wits.ac.za}

\author[H.~Prodinger]{Helmut Prodinger$^\dag{}$}
\thanks{$^\dag{}$Supported by NRF Grant 2053748}
\address{H.~Prodinger\\
The John Knopfmacher Centre for Applicable Analysis and Number Theory\\
 School of Mathematics\\
University of the Witwatersrand, P.~O. Wits\\ 
2050 Johannesburg, South Africa}
\email{helmut@maths.wits.ac.za}


\date{May 7, 2003}
\subjclass{Primary: 05A15; Secondary: 60C05}
\keywords{Pills problem, recursions, harmonic numbers, generating functions}
\begin{document}
\begin{abstract}
We revisit the pills problem proposed by Knuth and McCarthy.
In a bottle there are $m$ large pills and $n$ small pills. 
The large pill is equivalent to two small pills. 
Every day a person chooses a pill at random. 
If a small pill is chosen, it is eaten up, if a large pill is chosen it is broken 
into two halves, one half is eaten and the other half which is now considered 
to be a small pill is returned to the bottle. How many pills are left, on 
average, when the last large pill has disappeared? We show how to compute the moments,
in particular the variance, and then generalize the problem in various ways.
\end{abstract}


\maketitle

\section{Introduction}

The following problem was  proposed by Donald E. Knuth and John McCarthy
in \cite[p.~264]{KnMcC91}; the solution appeared in \cite[p.~684]{KnMcC92}.

\textsl{``In a bottle there are $m$ large pills and $n$ small pills. 
The large pill is equivalent to two small pills. 
Every day a person chooses a pill at random. 
If a small pill is chosen, it is eaten up, if a large pill is chosen it is broken 
into two halves, one half is eaten and the other half which is now considered 
to be a small pill is returned to the bottle.''}

The proposed problem  was to find the expected number of small pills remaining when 
there are no more large pills left in the bottle.

We will revisit this problem by showing how one can derive all the moments, at
least in principle. We compute them up to the third moment; in particular we
get the variance.

Then we generalize the problem in various ways: We assume two types of pills, 
but this time the large one is equivalent to $p$ of the small ones. 
When it is chosen, one unit is eaten, and the remaining $p-1$ small ones
are put back into the bottle.

We also consider the model of three types of pills. The large ones are equivalent
to 3 small ones, and the middle ones to 2 small ones. The simpler model is 
as follows:
when the large one is chosen, one unit is eaten, and the rest goes back to the bottle
in the form of 2 small units. The more complicated model is the one
where what remains from
a large pill goes back into the bottle in the form of a middle sized pill.

The methodology is as follows: Recursions are set up that describe the process.
Often, the solutions are found by inspection. For this, computer software is available in the form of \textsc{Gfun}, a Maple package described in \cite{SaZi94}. Once the answers are known, the proofs are a routine verification. One just has to check that the recursion and the
initial conditions are satisfied.

However, this approach has its limitations when the problems become too large,
and then the solutions had to be found in a different way. By inspection again, 
one can make a meaningful \textit{ansatz} and compute the coefficients 
involved in it by unwinding the recursions which they satisfy.


\section{The original problem revisited}

It is meaningful to introduce
 generating functions $f_{m,n}(u)$
where the coefficient of $u^k$ is the probability that $k$ small pills are left in
the bottle after the last large pill has been selected. The following recursion
is self-explanatory:
%
\begin{equation*}
f_{m,n}(u)=\frac{m}{m+n}f_{m-1,n+1}(u)+\frac{n}{m+n}f_{m,n-1}(u)\quad
\textrm{for $m\ge1$ \quad and}\quad
f_{0,n}(u)=u^n.
\end{equation*}
Quantities with negative indices are defined to be zero. 

We cannot hope to solve this system explicitly but it can be used to compute
moments by differentiations with respect to $u$, and then setting $u=1$.
Since the parameter $u$ appears explicitly only in the initial conditions, we
can compute the $s$th moment by using the same recursion and replacing the initial
conditions by $f_{0,n}=n^s$.

Let us first revisit the expected value $e_{m,n}:=f^{\prime}_{m,n}(1)$.
 \[
e_{m,n}=\frac{m}{m+n}\,e_{m-1,n+1}+\frac{n}{m+n}\,e_{m,n-1}
\quad
\textrm{for $m\ge1$ \quad and}\quad
e_{0,n}=n.\]

One can see quickly that
$e_{m,n}=\dfrac{n}{m+1}+H_m$,
where $H_m$ is the harmonic number defined as
\[
H_m=1+\frac{1}{2}+\frac{1}{3}+\cdots+\frac{1}{m}.
\]
A general reference for the harmonic numbers and many of their properties 
is \cite{gkp94}.
As mentioned before, to check that this is indeed the solution is a routine task.

Now let us move to the second moments which we denote by $e^{(2)}_{m,n}$.
We have the recursion
\[
e^{(2)}_{m,n}=\frac{m}{m+n}\,e^{(2)}_{m-1,n+1}+\frac{n}{m+n}\,e^{(2)}_{m,n-1}
\quad\textrm{and}\quad
e^{(2)}_{0,n}=n^2.
\]

If one computes a few of these values one sees some patterns evolving. 
For fixed $m$, the numbers $e^{(2)}_{m,n}$ are given by a quadratic polynomial in
$n$. Maple's \textsc{interpolate} command is useful to find them. Here is
a list of the first few instances:
%
\begin{align*}
e^{(2)}_{0,n}&=n^2,\\
e^{(2)}_{1,n}&=\frac{1}{3}n^2+\frac{7}{6}n+1,\\
e^{(2)}_{2,n}&=\frac{1}{6}n^2+\frac{23}{18}n+\frac{5}{2},\\
e^{(2)}_{3,n}&=\frac{1}{10}n^2+\frac{37}{30}n+\frac{71}{18},\\
e^{(2)}_{4,n}&=\frac{1}{15}n^2+\frac{29}{25}n+\frac{95}{18},\\
e^{(2)}_{5,n}&=\frac{1}{21}n^2+\frac{38}{35}n+\frac{2927}{450}, \ \text{\&c.}\\
\end{align*}

The coefficients of $n^2$ are clearly given by
$$
\frac{2}{(m+1)(m+2)}.
$$

It is not so obvious to find the general law of the coefficients of $n$.
However, here, \textsc{Gfun} is helpful. We form the list
\begin{equation*}
L=\bigg[
\frac76,{\frac {23}{18}},{\frac {37}{30}},{\frac {29}{25}},{\frac {38}{
35}},{\frac {997}{980}},{\frac {4819}{5040}},{\frac {5113}{5670}},{
\frac {59141}{69300}},{\frac {30883}{38115}}
,\dots
\bigg]
\end{equation*}
of these coefficients.
Using a few intermediate steps that we are not going to describe here, we find
a differential equation for the generating function $f(z)$ of the coefficients
of $L$:
\begin{align*}
0&=(180{z}^{3}-256{z}^{2}+92z)f(z)+
(360{z}^{4}-680{z}^{3}+412{z}^{2}-84z)f^{\prime}(z)\\
&+(150{z}^{5}-332{z}^{4}+238{z}^{3}-56{z}^{2})
f^{\prime\prime}(z)+\left (15{z}^{6}-37{z}^{5}+29{z}^{4
}-7{z}^{3}\right )f^{\prime\prime\prime}(z),\\
& f(0)=\frac76,\quad f^{\prime}(0)={\frac {23}{
18}},\quad f^{\prime\prime}(0)={\frac {37}{15}}.
\end{align*}
The solution is
\[f(z)=\frac{2\log^2 (1-z)}{z^2}-\frac{2\log (1-z)}{z^3}
+\frac{5\log (1-z)}{z^2}+\frac{4\dilog (1-z)}{z^2}-\frac{2}{z^2}.\]
Let us recall the $\dilog$ function:
\begin{equation*}
\dilog(1-z)=\sum_{n\ge1}\frac{z^n}{n^2}.
\end{equation*}
As a check, we compute the series expansion
\[
f(z)=\frac{7}{6}+\frac{23}{18}z+\frac{37}{30}z^2+\frac{29}{25}z^3+\frac{38}{35}z^4+
\frac{997}{980}z^5+\frac{4819}{5040}z^6+\frac{5113}{5670}z^7+\frac{59141}{69300}z^8
+O(z^9).
\]
We can thus find the coefficient of the linear term in $e^{(2)}_{m,n}$
as $[z^{m-1}]f(z)$.
The relevant formul{\ae}, that we need are
$$
[z^n]\log (1-z)=-\frac{1}{n},\quad
[z^n]{\log^2 (1-z)}=\frac{2}{n}H_{n-1},\quad
[z^n]\dilog(1-z)=\frac{1}{n^2}.
$$
Eventually we get
$$
[z^{m-1}]f(z)=\frac{2}{m+2}-\frac{5}{m+1}+\frac{4H_{m+1}}{m+1}.
$$

A similar procedure will now be applied to the sequence of the constants
$$\bigg
[1,\frac52,{\frac {71}{18}},{\frac {95}{18}},{\frac {2927}{450}},{\frac {
3437}{450}},{\frac {191633}{22050}},{\frac {1706629}{176400}},{\frac {
16826581}{1587600}},{\frac {18210313}{1587600}},\dots\bigg],
$$
leading to the generating function
\begin{align*}
f(z)&=\frac{3\log (1-z)}{z}+\frac{4\dilog(1-z)}{z}+\frac{2\log^2 (1-z)}{z}\\
&+\frac{3\log (1-z)}{1-z}+\frac{4\dilog(1-z)}{1-z}+\frac{2\log^2 (1-z)}{1-z}.
\end{align*}
As before we seek the coefficient of $z^{m-1}$.
We need additional formul\ae, see \cite{grkn82},
\[
[z^n]\frac{\dilog(1-z)}{1-z}=H^{(2)}_n,\quad \textrm{where}\quad H^{(r)}_m
=\sum_{j=1}^m\frac{1}{j^r},
\]
and
\[[z^n]\frac{\log^2 (1-z)}{1-z}=2{\sum_{p=2}^n\frac{1}{p}H_{p-1}}
=H^2_n-H^{(2)}_n.\]
This leads to
$$
[z^{m-1}]f(z)=-3H_m+2H^2_m+2H^{(2)}_m.
$$
Hence we found the second moment
$$
e^{(2)}_{m,n}=\dfrac{2}{(m+1)(m+2)}n^2+
\left[\dfrac{4H_{m+1}}{m+1}+\dfrac{2}{m+2}-\dfrac{5}{m+1}\right]n
-3H_m+2H^2_m+2H^{(2)}_m;
$$
the correctness of the formula was then checked subsequently by induction.

This gives us the variance $v_{m,n}=e^{(2)}_{m,n}-e^{2}_{m,n}$
as a corollary:
$$
v_{m,n}=\dfrac{m}{(m+1)^2(m+2)}n^2+\left[\dfrac{2H_m}{m+1}-\dfrac{m(3m+7)}
{(m+1)^2(m+2)}\right]n-3H_m+H^2_m+2H^{(2)}_m.
$$

Let us now discuss the third moment. We tried a similar approach as before, but
the order of the 
differential equations was too high, so we came up with  something different.

It is a meaningful guess to set
\[e^{(3)}_{m,n}:=a_mn^3+b_mn^2+c_mn+d_m.\]
Plugging that into the  recurrence
$$
e^{(3)}_{m,n}=\frac{m}{m+n}e^{(3)}_{m-1,n+1}+\frac{n}{m+n}e^{(3)}_{m,n-1}
\quad\text{and}\quad e^{(3)}_{0,n}:=n^3,
$$
we find
\begin{align*}
a_mn^3+b_mn^2+c_mn+d_m&=
\frac{m}{m+n}\left[a_{m-1}(n+1)^3+b_{m-1}(n+1)^2+c_{m-1}(n+1)+d_{m-1}\right]\\
&+\frac{n}{m+n}\left[a_m(n-1)^3+b_m(n-1)^2+c_m(n-1)+d_m\right].
\end{align*}
We  multiply this by $m+n$ and then compare 
coefficients of $n^i$, $i=0,\dots,4$:
\begin{alignat*}{2}
(m+3)a_m&=ma_{m-1},\qquad &a_0&=1,\\
(m+2)b_m&=3ma_{m-1}+mb_{m-1}+3a_m,\qquad &b_0&=0,\\
(m+1)c_m&=3ma_{m-1}+2mb_{m-1}+mc_{m-1}-a_m+b_m,\qquad &c_0&=0,\\
d_m&=a_{m-1}+b_{m-1}+c_{m-1}+d_{m-1},\qquad &d_0&=0.
\end{alignat*}
These are all first order recursions and can be solved one after the
other by iteration. The $a_m$'s are the easiest ones:
\[a_m=\frac{6}{(m+1)(m+2)(m+3)}.\]
Plugging that in we find the equation for the $b_m$'s:
\[(m+1)(m+2)b_m=m(m+1)b_{m-1}+\frac{18(m+4)}{(m+2)(m+3)},\]
and eventually
\[b_m=\frac{18}{(m+1)(m+2)}\left(H_m-\frac{m(7m^2+36m+41)}{6(m+1)(m+2)(m+3)}\right).\]\\
The equation for the $c_m$'s is then
\begin{align*}
(m+1)c_m&=mc_{m-1}+\left(\frac{54}{m+1}-\frac{18}{m+2}\right)H_m
\\&+\frac{54}{(m+1)^2}
-\frac{93}{m+1}-\frac{18}{(m+2)^2}+\frac{63}{m+2}-\frac{12}{m+3}.
\end{align*}
Iterating this we find
\begin{align*}
(m+1)c_m&=\sum_{j=1}^m\frac{54}{j+1}H_j-\sum_{j=1}^m\frac{18}{j+2}H_j
+\sum_{j=1}^m\frac{54}{(j+1)^2}-\sum_{j=1}^m\frac{93}{j+1}\\
&-\sum_{j=1}^m\frac{18}{(j+2)^2}+\sum_{j=1}^m\frac{63}{j+2}-\sum_{j=1}\frac{12}{j+3}.
\end{align*}
Using formul\ae{} such as
\[\sum_{j=1}^m\frac{1}{j+1}H_j=\frac{1}{2}H^{2}_m-\frac{1}{2}H^{(2)}_m+\frac{1}{m+1}H_m\]
and
\[\sum_{j=1}^m\frac{1}{j+2}H_j=
\frac{1}{2}H_m^2-\frac{1}{2}H^{(2)}_m+\frac{2m+3}{(m+1)(m+2)}H_m-\frac{m}{m+1}\]
we get eventually
\begin{align*}
c_m&=\frac{18}{m+1}\left(H^2_m+H^{(2)}_m\right)
-\frac{6(7m^2+18m+5)}{(m+1)^2(m+2)}H_m\\&
+\frac{m(7m^4+42m^3+40m^2-150m-227)}{(m+1)^3(m+2)^2(m+3)}.
\end{align*}
The equation for the $d_m$'s becomes after a few simplifications
\[
d_m=d_{m-1}-\frac{42}{m}H_{m-1}+\frac{36}{m^2}H_{m-1}+\frac{18}{m}H^2_{m-1}+\frac{18}{m}H^{(2)}_{m-1}
+\frac{36}{m^3}-\frac{42}{m^2}+\frac{7}{m}
\]
and furthermore
\begin{align*}
d_m
&=-42\sum_{j=1}^m\frac{1}{j}H_j+18\sum_{j=1}^m\frac{1}{j}\left[H^2_j+H^{(2)}_j\right]+7H_m.\\
\end{align*}
We use
$$
\sum_{j=1}^m\frac{1}{j}H_j=\frac{1}{2}\left(H^2_m+H^{(2)}_m\right)\ \text{and}\
\sum_{j=1}^m\frac{1}{j}\left[H^2_j+H^{(2)}_j\right]=\frac{1}{3}
\left[H^3_m+3H_mH^{(2)}_m+2H^{(3)}_m\right],
$$ 
which, once they are known, are not difficult to prove by induction.
(Carsten Schneider's page \textsc{Sigma} \cite{Sc01} is helpful for such
calculations.)
%
So finally
\begin{align*}
d_m=7H_m-21H^2_m+6H^3_m-21H^{(2)}_m+12H^{(3)}_m+18H_mH^{(2)}_m.
\end{align*}

\noindent The third moment is therefore\\
\begin{align*}
e^{(3)}_{m,n}&=\frac{6}{(m+1)(m+2)(m+3)}n^3\\
&+\frac{18}{(m+1)(m+2)}\left(H_m-\frac{m(7m^2+36m+41)}
        {6(m+1)(m+2)(m+3)}\right)n^2\\
&+\bigg[\frac{18}{m+1}\left(H^2_m+H^{(2)}_m\right)-\frac{6(7m^2+18m+5)}{(m+1)^2(m+2)}H_m
\\&
\hspace*{4cm}+\frac{m(7m^4+42m^3+40m^2-150m-227)}{(m+1)^3(m+2)^2(m+3)}\bigg]n\\
&+7H_m-21H^2_m+6H^3_m-21H^{(2)}_m+12H^{(3)}_m+18H_mH^{(2)}_m.
\end{align*}


\section{The large pill equals $p$  small pills}


We shall now look at a variation of the previous problem.
Suppose that we have $m$ large pills and $n$ small pills, 
but this time a large pill is equivalent to $p$ small pills.
If a large pill is chosen, it is broken into $p$ parts, 
one part is swallowed and the other $p-1$ parts
 are returned to the bottle. If a small pill is chosen it is simply swallowed.
The question is again the average number of small pills when the last large
pill has been chosen. 

Let us start with the instance $p=3$.
 The recurrence relation for this average given that we start with 
$m$ resp. $n$ pills, is 
\[e_{m,n}=\frac{m}{m+n}e_{m-1,n+2}+\frac{n}{m+n}e_{m,n-1}\quad\text{
for $m\ge1$, \ \ and}\quad
e_{0,n}=n.
\]
One quickly finds that 
$e_{m,n}=a_mn+b_m$. The first few values of $a_m$ for $m\ge1$
are given by
\[
\frac{1}{2},\frac{1}{3},\frac{1}{4},\frac{1}{5},\frac{1}{6},\frac{1}{7},\frac{1}{8},
\frac{1}{9},\frac{1}{10},\frac{1}{11},
\]
whence $a_m=1/(m+1)$; the $b_m$'s are a bit more complicated:
\[
2,3,\frac{11}{3},\frac{25}{6},\frac{137}{30},\frac{49}{10},\frac{363}{70},
\frac{761}{140},\frac{7129}{1260},\frac{7381}{1260}.
\]
However, with the help of \textsc{Gfun}, we have no trouble to find
the formula
$$
e_{m,n}=\frac{n}{m+1}+2H_m.
$$
Now we turn to the general case, with the recursion
\[e_{m,n}=\frac{m}{m+n}e_{m-1,n+p-1}+\frac{n}{m+n}e_{m,n-1}.
\]
Working out a few cases, one sees the formula
\[e_{m,n}=\frac{n}{m+1}+(p-1)H_m,\]
which subsequently can be proved by induction. We leave such a proof to
the interested reader. 

\section{Three types of pills}



Let us now look at a similar pill problem as before, 
but this time it involves three types of
pills. In a bottle there are $m$ large pills, $n$ medium pills and $p$ small pills. 
The large pill is equivalent to three small pills, the medium pill to two small 
pills. 

\subsection*{The simple model}
A pill is chosen at random and if it is a large pill, 
it is broken into three parts, one part is swallowed and the 
remaining two parts which are equivalent to two small pills are 
returned to the bottle. If the chosen pill is a medium pill, 
it is broken into two parts, one part is swallowed and the 
remaining part which is equivalent to one small pill is returned to the bottle.






We ask ourselves three questions:
What is the expected number of small pills left over when

\begin{itemize}
\item[i)] there are no more large pills,
\item[ii)] there are no more medium pills,
\item[iii)] there are no more large pills and no more medium pills?
\end{itemize}

The recursion
\[e_{m,n,p}:=\frac{m}{m+n+p}e_{m-1,n,p+2}+\frac{n}{m+n+p}e_{m,n-1,p+1}+
\frac{p}{m+n+p}e_{m,n,p-1}
\]
is the same for all three instances, only the initial conditions and the
range of its validity differ.

In the case i) it holds for $m\ge1$, $n,p\ge0$, and $e_{0,n,p}=p$.

In the case ii) it holds for $n\ge1$, $m,p\ge0$, and $e_{m,0,p}=p$.

In the case iii) it holds for $m,n,p\ge0$, $(m,n)\neq (0,0)$ and $e_{0,0,p}=p$.

Let us start with case i). We proceed as before and calculate $e_{m,n,p}$ for
$m=1,2,\dots$, with the results
\begin{align*}
e_{1,n,p}&=\frac{p}{2}+\frac{n}{4}+2,\\ 
e_{2,n,p}&=\frac{p}{3}+\frac{5n}{18}+3,\\ 
e_{3,n,p}&=\frac{p}{4}+\frac{13n}{48}+\frac{11}3,\\ 
e_{4,n,p}&=\frac{p}{5}+\frac{77n}{300}+\frac{25}{6},\\
e_{5,n,p}&=\frac{p}{6}+\frac{29n}{120}+\frac{137}{30}.
\end{align*}

It is not too hard to guess that $e_{m,n,p}=\frac{p}{m+1}+a_mn+2H_m$, and we
will find the sequence $a_m$ by plugging that into the recursion 
\[
e_{m,n,p}=\frac{m}{m+n+p}e_{m-1,n,p+2}+\frac{n}{m+n+p}e_{m,n-1,p+1}
+\frac{p}{m+n+p}e_{m,n,p-1}.
\]
From this we find eventually that
\[
(m+1)a_m=ma_{m-1}+\frac{1}{m+1}.
\]
The solution is again found by iteration:
\[
a_m=\frac{H_{m+1}-1}{m+1}.
\]

So our final result for the expected number of small pills left over
 when there are no more large pills is
\[e_{m,n,p}=\frac{p}{m+1}+\frac{H_{m+1}-1}{m+1}\,n+2H_m.
\]

In the case ii) the computations are very similar and the result is
\[e_{m,n,p}=\frac{p}{n+1}+\frac{2(H_{n+1}-1)}{n+1}m+H_n.
\]

The case iii) is a bit more involved.
First we note that $e_{0,n,p}=\frac{p}{n+1}+H_n$, which is just the
original pill problem.

We do a few more computations:
\begin{align*}
e_{1,n,p}&=\frac{p}{n+2}+\frac{n+2}{n+1}H_{n+1},\\
e_{2,n,p}&=\frac{p}{n+3}+\frac{n+4}{n+2}H_{n+2},\\
e_{3,n,p}&=\frac{p}{n+4}+\frac{n+6}{n+3}H_{n+3},\\
\end{align*}
from which the general formula can be seen:
\[e_{m,n,p}=\frac{p}{m+n+1}+\frac{2m+n}{m+n}H_{m+n};
\]
for $m=n=0$, it must be interpreted as $e_{0,0,p}=p$.

\subsection*{The difficult model}

We still have the three types of pills, the large, the medium, and the small. 
The way in which they are broken up and swallowed is different. This time,
 if a large pill is chosen, it is broken into two parts in the ratio $1:2$, 
the small part which is equivalent to a small pill is swallowed and the remaining
 part which is equivalent to a medium pill is returned to the bottle. However 
if a medium pill is chosen it is broken into two equal parts, both equivalent to a 
small pill, one part is swallowed and the other part is returned as a small pill 
into the bottle. If a small pill is chosen it is simply swallowed.

We look again at the expected number of small pills left over
but this time only for two 
cases, since the third case makes no sense in this setting:

\begin{itemize}
\item[i)] when there are no more large pills,
\item[ii)] when there are no large pills and no medium pills left.
\end{itemize}

We start with case i). The recursion  is
\[e_{m,n,p}:=\frac{m}{m+n+p}e_{m-1,n+1,p}
+\frac{n}{m+n+p}e_{m,n-1,p+1}+\frac{p}{m+n+p}e_{m,n,p-1},\]
for $m\ge1$ and  the initial conditions $e_{0,n,p}=p$.

Again, working out a few cases for $m=1,2,\dots$, one arrives
at
\[
e_{m,n,p}=\frac{p}{m+1}+\frac{H_{m+1}-1}{m+1}\,n+a_m,
\]
and we will determine $a_m$ by plugging that into the recursion, which leads to
\[
ma_m=ma_{m-1}+H_{m}-1\quad\text{and}\quad a_0=0.
\]
This has the solution (again by iteration)
\begin{align*}
a_m=\sum_{i=1}^m\left(\frac{H_i}{i}-\frac{1}{i}\right)
=\frac{H^2_m+H^{(2)}_m}{2}-H_m.
\end{align*}
%
Hence our final answer is
\[
e_{m,n,p}=\frac{p}{m+1}+\frac{H_{m+1}-1}{m+1}\,n+\frac{{{H}_m^2+{H}_m^{(2)}}}{2}-H_m.
\]

Now we turn to the second case, counting the number of small pills remaining when
the other ones have disappeared.
In this case the recursion holds for all $(m,n)\neq(0,0)$, and
$e_{0,0,p}=p$.
%
It turns out that 
\begin{equation*}
e_{m,n,p}=\lambda_{m,n}p+\mu_{m,n}.
\end{equation*}
Using that in the recursion, one finds
\begin{equation*}
(m+n+1)\lambda_{m,n}=m\lambda_{m-1,n+1}+n\lambda_{m,n-1}
\quad\text{for $m\ge1$\quad and}\quad \lambda_{0,n}=\frac1{n+1}.
\end{equation*}
Multiplying that with $(m+n)^{\underline{m}}:=(m+n)\dots(n+1)$ we get
\begin{equation*}
(m+n+1)^{\underline{m+1}}\,\lambda_{m,n}=(m+n)^{\underline{m+1}}\,\lambda_{m,n-1}
+m(m+n)^{\underline{m}}\,\lambda_{m-1,n+1}.
\end{equation*}
Iterating this we find, noticing that $\lambda_{m,0}=\frac{m}{m+1}\lambda_{m-1,1}$, 
\begin{equation*}
\lambda_{m,n}=\frac m{(m+n+1)^{\underline{m+1}}}
\sum_{k=1}^{n+1}(k+m-1)^{\underline{m}}\,\lambda_{m-1,k}.
\end{equation*}
The recursion looks nicer in terms of
$\Lambda_{m,n}:=(m+n)^{\underline{m+1}}\lambda_{m,n}$:
\begin{equation*}
\Lambda_{m,n}=\frac{m\,n}{m+n+1}
\sum_{k=1}^{n+1}\Lambda_{m-1,k}.
\end{equation*}
One could now iterate this with respect to $m$ and thus write $\Lambda_{m,n}$
as an $m$-fold iterated sum.



Let us compute the first three $\lambda_{m,n}$'s:\\
For $m=1$ we get
\begin{align*}
\lambda_{1,n}
&=\frac1{(n+2)(n+1)}\sum_{k=1}^{n+1}k\lambda_{0,k}\\
&=\frac1{(n+2)(n+1)}\bigg[n+2-H_{n+2}\bigg]=
\frac1{n+1} -\frac{H_{n+2}}{(n+2)(n+1)}.
\end{align*}
				%
For $m=2$ we obtain
\begin{align*}
\lambda_{2,n}&=
\frac{2}{(n+1)(n+2)(n+3)}\sum_{k=1}^{n+1}\left(k-\frac{k}{k+2}H_{k+2}\right)\\
&=\frac{1}{n+1}-\frac{2(n+4)H_{n+3}}{(n+1)(n+2)(n+3)}+\frac{2\left[H_{n+3}^2+H^{(2)}_{n+3}\right]
}{(n+1)(n+2)(n+3)}.
\end{align*}
For $m=3$ we have
\begin{multline*}
\lambda_{3,n}=\frac{1}{n+1}+\\\frac{3\left[(2n+13)(H_{n+4}^2+H_{n+4}^{(2)})-(n+4)(n+5)H_{n+4}-2H_{n+4}^3-6H_{n+4}H_{n+4}^{(2)}-4H_{n+4}^{(3)}\right]}{(n+1)(n+2)(n+3)(n+4)}.
\end{multline*}
The results look rather cumbersome, so we stop here for the $\lambda$'s.\\

Now for the $\mu$'s, if we use $e_{m,n,p}=\lambda_{m,n}p+\mu_{m,n},$ in the recursion, we find
\begin{equation*}
(m+n)\mu_{m,n}=n\mu_{m,n-1}+m\mu_{m-1,n+1}+n\lambda_{m,n-1}
\quad\text{for $n\ge1$\quad and}\quad \mu_{0,n}=H_n.
\end{equation*}
Multiplying that with $\dfrac{(m+n-1)!}{n!}$ we get
\[\frac{(m+n)!}{n!}\mu_{m,n}=\frac{(m+n-1)!}{(n-1)!}\mu_{m,n-1}+\frac{m(m+n-1)!}{n!}\mu_{m-1,n+1}+\frac{(m+n-1)!}{(n-1)!}\lambda_{m,n-1}.\]
Iterating this and noticing that $\mu_{m,0}=\mu_{m-1,1}$, we obtain 
\[\mu_{m,n}=\frac{n!}{(m+n)!}\left[\sum_{i=0}^n\frac{m(m+i-1)!}{i!}\mu_{m-1,i+1}+\sum_{i=1}^n\frac{(m+i-1)!}{(i-1)!}\lambda_{m,i-1}\right].\]
Let us compute the first two values of $\mu_{m,n}$.\\
For $m=1$ we have
\[\mu_{1,n}=\frac{1}{n+1}\left[(n+2)H_{n+1}-\frac{H^2_{n+1}+H^{(2)}_{n+1}}{2}\right].\]
For $m=2$ we have
\[\mu_{2,n}=\frac{(n+2)(n+3)H_{n+2}-(n+5)(H_{n+2}^2+H_{n+2}^{(2)})+H_{n+2}^3+3H_{n+2}H_{n+2}^{(2)}+2H_{n+2}^{(3)}}{(n+1)(n+2)}.\]
Again the results are not attractive so we stop here.\\
Table 1 shows a few values of $e_{m,n,p}$ for that instance.
\begin{scriptsize}
{\renewcommand{\arraystretch}{3}
\begin{sidewaystable}[p!]
\caption{The first few values of $e_{m,n,p}.$}
\begin{tabular}{|c||c|c|c|c|c|}%|c|*3{p{8.5em}|}}
\hline 
$m \raisebox{-0.02cm}{\text{\begin{large}{${\setminus}$}\end{large}}} n$&0&1&2&3\\
\hline\hline
0&
$p$&
$\frac{1}{2}p+1$&
$\frac{1}{3}p+\frac{3}{2}$
&$\frac{1}{4}p+\frac{11}{6}$\\
\hline
1&
$\frac{1}{4}p+1$&
$\frac{7}{36}p+\frac{11}{8}$&
$\frac{23}{144}p+\frac{179}{108}$
&$\frac{163}{1200}p+\frac{1085}{576}$\\
\hline
2&
$\frac{7}{54}p+\frac{11}{8}$&
$\frac{97}{864}p+\frac{347}{216}$&
$\frac{5359}{54000}p+\frac{2075}{1152}$
&$\frac{9619}{108000}p+\frac{1417813}{720000}$\\
\hline
3&
$\frac{97}{1152}p+\frac{347}{216}$&
$\frac{54997}{720000}p+\frac{12259}{6912}$&
${\frac{151187}{2160000}}p+\frac{10376087}{5400000}$
&$\frac{111775621}{1728720000}p+\frac{44369339}{21600000}$\\
\hline
4&
$\frac{54997}{900000}p+\frac{12259}{6912}$&
$\frac{460463}{8100000}p+\frac{41129339}{21600000}$&
$\frac{7241046271}{136136700000}p+\frac{78657551}{38880000}$
&$\frac{290480276101}{5808499200000}p+\frac{16255908076501}{7623655200000}$\\
\hline
5&
$\frac{460463}{9720000}p+\frac{41129339}{21600000}$&
$\frac{51185279267}{1143548280000}p+\frac{390968681}{194400000}$&
$\frac{24851973140539}{585496719360000}p+\frac{112626727006711}{53365586400000}$
&$\frac{17235071397413921}{426827108413440000}p+\frac{120391962039184219}{54646360473600000}$\\
\hline
6&
$\frac{51185279267}{1334139660000}p+\frac{390968681}{194400000}$&
$\frac{200170674968477}{5464636047360000}p+\frac{336486120012803}{160096759200000}$&
$\frac{628506421179609521}{17926738553364480000}p+\frac{204870652213546079}{93679475097600000}$&
$\frac{10048104718893889789}{298778975889408000000}p+
\frac{9754683417101073727043}{4302417252807475200000}$\\
\hline
%7&
%$\frac{200170674968477}{6245298339840000}p+\frac{336486120012803}{160096759200000}$&
%$\frac{11369422537341929933}{368778621669212160000}p+\frac{2859481756726972261}{131151265136640%0000}$&
%$\frac{1369430873074088345537}{46097327708651520000000}p+\frac{65481847688598369594119}{290413%16456450457600000}$\\
%\hline
\end{tabular}
\end{sidewaystable}}
\end	{scriptsize}
\newpage



\bibliographystyle{plain}
\begin{thebibliography}{1}

\bibitem{KnMcC92}
T.~Hesterberg et~al.
\newblock Problems and solutions, {E}3429.
\newblock {\em American Mathematical Monthly}, 99(7):684, 92.

\bibitem{gkp94}
R.~L. Graham, D.~E. Knuth, and O.~Patashnik.
\newblock {\em Concrete Mathematics}.
\newblock Addison-Wesley, 1994.

\bibitem{grkn82}
D.~Greene and D.~E. Knuth.
\newblock {\em Mathematics for the Analysis of Algorithms}.
\newblock Boston: Birkhauser, 1982.

\bibitem{KnMcC91}
D.~E. Knuth and J.~McCarthy.
\newblock Problem {E}3429: Big pills and little pills.
\newblock {\em American Mathematical Monthly}, 98(3):264, 91.

\bibitem{SaZi94}
B.~Salvy and P.~Zimmermann.
\newblock Gfun: a {M}aple package for the manipulation of generating and
  holonomic functions in one variable.
\newblock {\em ACM Transactions on Mathematical Software}, 20(2):163--177,
  1994.

\bibitem{Sc01}
C.~Schneider.
\newblock {\em Symbolic Summation in Difference Fields}.
\newblock PhD thesis, RISC, J Kepler University, Linz, 2001.

\end{thebibliography}


%\bibliography{charlbib}


\end{document}


