\documentclass[a4paper,reqno,12pt]{amsart}

\usepackage{latexsym}
%\usepackage{emlines2}




\marginparwidth 0pt \oddsidemargin -0.4 in \evensidemargin -0.4 in
\marginparsep 0pt \topmargin -0.3 in \textwidth 7 in \textheight 9
in

\hfuzz=30pt
\def\j{\mathbf{j}}

\newtheorem{theorem}{Theorem}%[section]
\newtheorem{proposition}[theorem]{Proposition}%[section]
\newtheorem{lemma}[theorem]{Lemma}%[section]
\newtheorem{cor}[theorem]{Corollary}
\renewcommand{\thecor}{}
\newtheorem{ncor}[theorem]{Corollary}

\newtheorem{df}{Definition}[section]
\newtheorem{ex}{Example}
\renewcommand{\theex}{}

\newcommand\U{{\mathcal{U}}}
\newcommand\E{{\mathbb{E}}}
\newcommand\V{{\mathbb{V}}}
\newcommand\T{{\mathcal{T}}}
\renewcommand\P{{\mathcal{P}}}
\renewcommand\a{{\bold{a}}}
\newcommand\p{{\bold{p}}}
\newcommand\mcirc{{\bigcirc}}

\newcommand\C{{\mathcal C}}

\newcommand\bl{\operatorname{L}}

\newcommand{\floor}[2]{\genfrac{\lfloor}{\rfloor}{}{}{#1}{#2}}
\newcommand{\ceil}[2]{\genfrac{\lceil}{\rceil}{}{}{#1}{#2}}

\renewcommand{\thefootnote}{\fnsymbol{footnote}}

\setcounter{topnumber}{5}
\renewcommand{\topfraction}{1.0}
\renewcommand{\textfraction}{.0}


%\renewcommand{\baselinestretch}{1.8}


\begin{document}

\title[Binary search tree recursions with harmonic toll functions]
{Binary search tree recursions with harmonic toll functions}
\author[A.~Panholzer and H.~Prodinger]
{Alois Panholzer$^{\dag}$ and Helmut Prodinger$^{\ddag}$}

\address{${}^{\dag}$Institut f\"ur Algebra und Computermathematik\\
Technische Universit{\"a}t Wien\\ Wiedner Hauptstra{\ss}e 8--10\\ A-
1040 Wien\\ Austria.}
\email{Alois.Panholzer@tuwien.ac.at}\vskip0pt
\address{${}^{\ddag}$The John Knopfmacher
Centre for Applicable Analysis and Number Theory\\ Mathematics
Department\\ University of the Witwatersrand\\ P.O. Wits\\ 2050
Johannesburg\\ South Africa.}
\email{helmut@gauss.cam.wits.ac.za\vskip0pt {\it WWW-address{\rm:
}}\texttt{http://www.wits.ac.za/helmut/index.htm}}

\date{January 18, 20001}
\begin{abstract}
We study recursions that are traditional in the context
of binary search trees and quicksort with
the (non--standard) toll function $H_n$.
\end{abstract}

\maketitle

\rightline{{\sl When an explicit formula can be}}
\rightline{{\sl  found there is nothing to beat it.}}
\rightline{{\sl In the rush to abstraction and }}
\rightline{{\sl generalization we often forget this.}}
\rightline{{\sf Richard Askey}}


\section{Introduction}

Recursions like
\begin{equation*}
a_n=t_n+\frac{1}{n}\sum_{k=0}^{n-1}(a_k+a_{n-1-k}),
\qquad n\ge1,\  a_0=1
\end{equation*}
with $t_n=1$ or $t_n=n$ are very common in analysis of algorithms
and data structures, especially in the context of binary search trees and
quicksort, see \cite{SeFl96}. The extra quantities $t_n$ are usually
called the {\it toll function}.

In the recent paper \cite{GrPr00} however,
quantities $t_n$ that have a logarithmic behaviour, occur.
(The paper \cite{Fill96} also studies logarithmic costs, but the
results are disjoint from ours.)
To be more specific, a typical recursion from
\cite{GrPr00} is

\begin{equation*}   
T_n=1+S_n+\frac{1}{n}\sum_{k=0}^{n-1}(T_k+T_{n-1-k}),
\end{equation*} 
where the quantities $S_n$ are themselves solutions of 
divide--and--conquer recursions, and turn out to have logarithmic 
behaviour.

Although such recursions might look exotic at a first glance, they
really are not, and appear naturally when e.~g. a 
(recursive) divide--and--conquer strategy has extra costs at each
level of recursion which are again given by a (different) recursion.
Therefore, the recursions of interest  occur typically when one
recursive algorithm calls another (usually simpler) recursive algorithm.


Thus, in this paper, we undertake a thorough study of such recursions,
where we consider instead of $t_n=\log n$  quantities like $t_n=H_n$, with
 harmonic numbers $H_n=\sum_{1\le k\le n}\frac1k$, which are much more
manageable and represent the discrete analogue of the logarithm.

We refrain from investigating more exotic toll functions like $\frac{n}{\sqrt{H_n}}$,
since they are not likely to occur as costs of (auxiliary) algorithms, used
as subroutines in recursive algorithms, and since the class of toll functions
that we consider is  a reasonably large one where one can still hope
for closed form solutions.

This paper can be seen both as a tuturial, showing the interested
reader methods to attack such recursions, but also as a 
{\em catalogue
of solutions}, as explained e.~g. in \cite{GrKn82}:
Since solutions are linear, with a suitably large collection of
basic solutions one can get ones own solution by constructing
an appropriate linear combination.







A good way of imagining what we are doing is by thinking of binary
search trees, where each node is weighted by a harmonic number, indexed by
the size of the subtree having the particular node as a root. The basic
parameter of interest is then the total weight of a binary search tree.
(The sum of the weights of all nodes in the binary search tree.)

\section{Total weight of a binary search tree\label{sect2}}

Let $f_n(v)$ be the  probability generating function of the
parameter ``total weight of a binary search tree,"
i.~e. $[v^k]f_n(v)$ is the probability that a random binary search
tree of size $n$ has weight $k$. Observe carefully that this is {\it not\/}
a polynomial, as the exponents will be in general some rational numbers.
That is, however, not a major drawback, although it is harder to set up
multivariate generating functions. In this section we are interested in
the average $a_n=f_n'(1)$ and variance $b_n+a_n-a_n^2$, with
$b_n=f^{\prime\prime}_n(1)$.

Up to a certain level of complexity, one can resort to elementary methods,
as we will demonstrate here. In later sections, however, the appropriate
methodology is to work with differential equations.

The recursion
\begin{equation*}
f_n(v)=v^{H_n}\frac{1}{n}\sum_{k=0}^{n-1}f_k(v)f_{n-1-k}(v),
\qquad n\ge1,\ f_0(v)=1
\end{equation*}
is quite straightforward; compare \cite{SeFl96} for very similar
recursions.

From that we find the prototype of a binary search tree recursion with a
harmonic toll function, viz.

\begin{equation*}
a_n={H_n}+\frac{2}{n}\sum_{k=0}^{n-1}a_k,
\qquad n\ge1,\ a_0=0.
\end{equation*}
One can find the solution $a_n=2n-H_n$  for instance by inspection,
which is not really a method but works well in practice when everything
is sufficiently simple. Typically, a routine proof by induction would
justify the answer. 
Now set $b_n=f_n^{\prime\prime}(1)$,
then
\begin{equation*}
b_n=H_n(H_n-1)+2H_n(a_n-H_n)+\frac{2}{n}\sum_{k=0}^{n-1}b_k
+\frac{2}{n}\sum_{k=0}^{n-1}a_ka_{n-1-k},\qquad n\ge1,\ b_0=0.
\end{equation*}
Plugging in the known values and simplifying (we spare the reader with the
simple but boring details) we find

\begin{align*}
nb_n-(n+1)b_{n-1}&=
4n^2
-H_{n}^2-3H_n+\frac{2}{n}H_n
-2H_{n}^{(2)}+\frac{1}{n^2}
-1,
\end{align*}
where also harmonic numbers of the second order occur.
The technique is explained in many textbooks on analysis of 
algorithms, e.~g. \cite{SeFl96}. The first step is to multiply by
$n$, so that the difference of two such equations, once indexed
by $n$ and once by $n-1$, makes the sum $\sum_{k=0}^{n-1}b_k$
disappear. 




The recursion may be rewritten as


\begin{align*}
\frac{b_n}{n+1}-\frac{b_{n-1}}{n}&=
4-\frac{4}{n+1}
-\frac{H_{n}^2}{n}+\frac{H_{n+1}^2}{n+1}
-\frac{2H_{n+1}}{(n+1)^2}
+\frac{2H_n}{n^2}-\frac{5H_n}{n}
+\frac{5H_{n+1}}{n+1}\\&-\frac{5}{(n+1)^2}
-\frac{2H_{n}^{(2)}}{n}+
\frac{2H_{n+1}^{(2)}}{n+1}-\frac{1}{(n+1)^3}
+\frac{1}{n^3}
-\frac{1}{n^2},
\end{align*}
which can now be summed (a telescoping sum), with the result

\begin{align*}
\frac{b_n}{n+1}&=
4n-4H_n-\frac{4}{n+1}+4
+
\frac{H_{n}^2}{n+1}
+
5\frac{H_{n}}{n+1}-
6H_{n}^{(2)}
+2\frac{H_{n}^{(2)}}{n+1},
\end{align*}
or

\begin{align*}
b_n&=
4n(n+2)-(4n-1)H_n
+
H_{n}^2-
2(3n+2)H_{n}^{(2)}.
\end{align*}
In the simplification process of these sums, the most important identity
is
\begin{equation*}
\sum_{k=1}^n\binom km H_k=\binom{n+1}{m+1}\Big(
H_{n+1}-\frac1{m+1}\Big),
\end{equation*} which can be found in \cite{Knuth68} and other textbooks. 



Thus the variance turns out to be
\begin{align*}
b_n+a_n-a_n^2&=
4n(n+2)-(4n-1)H_n
+
H_{n}^2-
2(3n+2)H_{n}^{(2)}\\&+2n-H_n-4n^2+4nH_n-H_n^2\\
&=
10n
-2(3n+2)H_{n}^{(2)}.
\end{align*}
Asymptotically, the variance is given by
\begin{equation*}
(10-\pi^2)n=0 .130395596\cdot n.
\end{equation*}

\section{Refined results a la quickselect (=Hoare's Find algorithm)}

Traditionally, one  associates the numbers $1,2,\dots,n$ with the nodes of
the binary search tree by traversing it in {\it inorder}. This reflects
also the intuition that smaller elements are to the left, and larger
elements  to the right of the root. If one uses the quicksort strategy
only to {\em find\/} element $j$, then the recursion of interest is



\begin{equation*}
P_{n,j}(v)=v^{*}\frac{1}{n}\bigg[
\sum_{k=1}^{j-1}P_{n-k,j-k}(v)+1+\sum_{k=j+1}^nP_{k-1,j}(v)
\bigg],\quad P_{1,1}(v)=v,
\end{equation*}
where $*$ stands for 1 in the instance ``number of passes,'' and for
$n-1$ in the instance ``number of comparisons." See for instance
\cite{KiPr98} and the references therein.

In the context of our study of harmonic recursions, it is thus naturally
to introduce the probability generating functions


\begin{equation*}
P_{n,j}(v)=v^{H_n}\frac{1}{n}\bigg[
\sum_{k=1}^{j-1}P_{n-k,j-k}(v)+1+\sum_{k=j+1}^nP_{k-1,j}(v)
\bigg],\quad P_{1,1}(v)=v.
\end{equation*}
The coefficient of $v^k$ in it is the probability that in a random binary
search tree of size $n$ the node $j$ is connected to the root via the path
$r_1,r_2,\dots,r_s=j$, and the sum of the weights of the nodes
$r_1,r_2,\dots,r_s $ is $k$.


The expectations are again given by $a_{n,j}=P_{n,j}'(1)$, or
\begin{equation*}
a_{n,j}={H_n}+\frac{1}{n}\bigg[
\sum_{k=1}^{j-1}a_{n-k,j-k}+\sum_{k=j+1}^na_{k-1,j}
\bigg],\quad a_{1,1}=1.
\end{equation*}
Although even here is it possible to solve it by elementary means, it is
the right time to introduce generating functions,
 \begin{equation*}
F(z,u)=\sum_{1\le j\le n}a_{n,j}z^nu^j. 
\end{equation*}
In order to translate the recursion into a differential equation, one must
first compute \begin{equation*} \sum_{1\le j\le n}nH_nz^{n-1}u^j=
\frac{u}{1-u}\big(S(z) -S(zu)\big) \end{equation*} with

\begin{equation*}
S(z)=\left(\frac{1}{1-z}\log\Big(\frac{1}{1-z}\Big)\right)'
=\frac{1}{(1-z)^2}\log\Big(\frac{1}{1-z}\Big)+\frac{1}{(1-z)^2}.
\end{equation*}
Hence we obtain the differential equation


\begin{equation*}
\frac{\partial}{\partial z}F(z,u)=
\frac{u}{1-u}\big(S(z)
-S(zu)\big)
+\bigg(\frac{u}{1-zu}+\frac{1}{1-z}\bigg)F(z,u),
\qquad F(0,u)=0
\end{equation*}
which is well within the range of Maple, and we obtain


\begin{align*}
F(z,u)&=\frac{1}{2(1-z)(1-zu)}\log^2\frac{1}{1-zu}
+\frac{u}{2(1-z)(1-zu)}\log^2\frac{1}{1-z}\\&
+\frac{2-u-zu}{(1-z)(1-zu)(1-u)}\log\frac{1}{1-zu}
+\frac{u(1+zu-2u)}{(1-z)(1-zu)(1-u)}\log\frac{1}{1-z}\\&
-\frac{2zu}{(1-zu)(1-z)}.
\end{align*}
Reading off coefficients is a (tedious) routine operation, with the result







\begin{align*}
a_{n,j}=\frac{1}{2}H_{j}^2+\frac{1}{2}H_{n+1-j}^2
-\frac{1}{2}H_{j}^{(2)}-\frac{1}{2}H_{n+1-j}^{(2)}
-H_n+2H_j+2H_{n+1-j}-2.
\end{align*}


One should observe the symmetry of the results (obvious from the
definition): for each term with $j$, there is also a corresponding
term with $n+1-j$.

To gain some asymptotic information, it is customary in this
context to set $j=\rho n$, with a fixed ratio $0<\rho<1$.
Then $a_{n,j}\sim \log^2n$.

Observe also that
\begin{align*}
\sum_{j=1}^na_{n,j}&=\sum_{j=1}^n
\Big(H_{j}^2-H_{j}^{(2)}\Big)-nH_n+4\sum_{j=1}^nH_j-2n\\
&=(n+1)\Big(H_{n}^2-H_{n}^{(2)}\Big)+2n-2nH_n-nH_n+4(n+1)H_n-
4n-2n\\
&=(n+1)\Big(H_{n}^2-H_{n}^{(2)}\Big)+nH_n+4H_n-4n.
\end{align*}

Now we turn to  the second factorial moments, which are given
by $b_{n,j} = P^{\prime\prime}_{n,j}(1)$.
We get here the recurrence
\begin{equation*}
b_{n,j} = H_{n} (H_{n}-1) + 2 \frac{H_{n}}{n} \left[ \sum_{k=1}^{j-1} a_{n-k, j-k} 
+ \sum_{k=j+1}^{n} a_{k-1,j}\right] + \frac{1}{n} \left[ \sum_{k=1}^{j-1} b_{n-k, j-k} 
+ \sum_{k=j+1}^{n} b_{k-1,j}\right], \enspace b_{1,1} = 0,
\end{equation*}
which simplifies by the above recursion for the $a_{n,j}$'s to
\begin{align*}
b_{n,j}&=-3H_n^2-5H_n+H_nH_{j}^2+H_nH_{n+1-j}^2-
H_nH^{(2)}_{j}-H_nH^{(2)}_{n+1-j}+4H_nH_j+4H_nH_{n+1-j}\\
&+\frac{1}{n}\sum_{k=1}^{j-1}b_{n-k,j-k}
+\frac{1}{n}\sum_{k=j+1}^{n}b_{k-1,j}.
\end{align*}

The solution can be obtained by the following considerations.
The recurrence
\begin{equation*}
   n b_{n,j} = \sum_{k=1}^{j-1} b_{n-k,j-k} + \sum_{k=j+1}^{n} b_{k-1,j} +
X_{n,j} \, , 
\end{equation*} 
where $X_{n,j}$ is an arbitrary toll function
with  ordinary generating function $X(z,u)$, leads to the differential
equation 
\begin{equation*}
   \frac{\partial}{\partial z} F(z,u) -
\left(\frac{1}{1-z}+\frac{u}{1-uz}\right)F(z,u) =
   \frac{X(z,u)}{z} \, ,
\end{equation*}
which has the solution
\begin{equation*}
   F(z,u) = \frac{1}{(1-z)(1-uz)}\int_{t=0}^{z} \frac{(1-t)(1-ut)}{t}
   X(t,u)
dt \, .
\end{equation*}
(In this case, Maple cannot do that, since the integral is
no longer  expressible in elementary functions.)

Extracting coefficients gives after some manipulations
\begin{equation*}
   b_{n,j} =
\sum_{k=1}^{j}\sum_{r=k}^{n-j+k}\frac{X_{r,k}}{r}-\sum_{k=1}^{j}\sum_{r=k+
1} ^{n-j+k}\frac{X_{r-1,k}}{r}
   - \sum_{k=2}^{j}\sum_{r=k}^{n-j+k}\frac{X_{r-1,k-1}}{r} +
\sum_{k=2}^{j}\sum_{r=k+1}^{n-j+k}\frac{X_{r-2,k-1}}{r} \, .
\end{equation*}
 In our case, the $X_{n,j}$'s are  given by
 \begin{align*}
    X_{n,j} &= -3 n H_{n}^{2} - 5 n H_{n} + n H_{n} H_{j}^2 + n H_{n}
H_{n+1-j}^2 \\&-
n H_{n} H^{(2)}_{j} - n H_{n} H^{(2)}_{n+1-j} + 4 n H_{n} H_{j} + 4 n
H_{n} H_{n+1-j} \, .
 \end{align*}
 The following observation is also of some interest:
If two toll functions $X_{n,j}$ resp.\ $\bar{X}_{n,j}$ are given by
\begin{equation*}
   X_{n,j} = C_{n} D_{j} \quad \text{and} \quad \bar{X}_{n,j} = C_{n}
D_{n+1-j} \, , \end{equation*} then for $b_{n,j}$ resp.\ $\bar{b}_{n,j}$
we have \begin{equation*}
   \bar{b}_{n,j} = b_{n,n+1-j} \, .
\end{equation*}

So the problem is reduced to  double summations for the terms $n H_{n}$, $n H_{n}^{2}$,
$n H_{n} H_{j}$, $n H_{n} H_{j}^{(2)}$ and $n H_{n} H_{j}^{2}$.
Because of  cancellations that occur
 a simple rewriting of $b_{n,j}$ is advantageous with
respect to an easier manipulation 
of the terms involved:
\begin{equation*}
   b_{n,j} = \sum_{k=1}^{j} \sum_{r=k+1}^{n-j+k} \frac{1}{r} (X_{r,k} - X_{r-1,k})
   - \sum_{k=2}^{j} \sum_{r=k+1}^{n-j+k} \frac{1}{r} (X_{r-1,k-1}-X_{r-2,k-1})
   + \sum_{k=2}^{j} \frac{1}{r} (X_{k,k} - X_{k-1,k-1}) + X_{1,1}.
\end{equation*}
For all the terms noted above, we can substitute the inner sums 
by expressions involving harmonic
numbers, and simplify the remaining sum  using  formul{\ae} like
\begin{align*}
   \sum_{k=1}^{n} \frac{H_{k}}{k} & = \frac{1}{2} H_{n}^{2} + \frac{1}{2} H_{n}^{(2)} \, , \\
   \sum_{k=1}^{n} \frac{H_{k}^{(2)}}{k^{2}} & = \frac{1}{2} (H_{n}^{(2)})^{2} + 
   \frac{1}{2} H_{n}^{(4)} \, , \\
   \sum_{k=1}^{n} \frac{H_{k}^{2}}{k} & = \frac{1}{3} H_{n}^{3} + H_{n} H_{n}^{(2)} + \frac{2}{3} H_{n}^{(3)} 
   - \sum_{k=1}^{n} \frac{H_{k}^{(2)}}{k} \, , \displaybreak[1] \\
   \sum_{k=1}^{j} \frac{H_{n-j+k}}{k} & + \sum_{k=1}^{j} \frac{H_{k}}{n-j+k} = 
   H_{n} H_{j} + \frac{1}{n-j} \left( H_{j} + H_{n-j} - H_{n}\right) \, , \\
   \sum_{k=1}^{j} \frac{H_{k}^{2}}{n-j+k} & = \sum_{k=1}^{j} \frac{H_{n-j+k}}{k^{2}} 
   - 2 \sum_{k=1}^{j} \frac{H_{k} H_{n-j+k}}{k} + \frac{2}{n-j} \sum_{k=1}^{j}
   \frac{H_{n-j+k}}{k} \\
   & + \frac{H_{j}^{2}}{n-j} + H_{n} H_{j}^{2} -
   \frac{2 H_{n} H_{j}}{n-j} - \frac{1}{(n-j)^{2}} \left( H_{j} + H_{n-j} - H_{n} \right) \, .
\end{align*}
Eventually we find
\begin{align*}
b_{n,j} & = \sum_{k=1}^{j}\frac{H_{k} H_{n-j+k}^{2}}{k} -
\sum_{k=1}^{j}\frac{H_{k}H_{n-j+k}^{(2)}}{k}
+4 \sum_{k=1}^{j}\frac{H_{k}H_{n-j+k}}{k} -
\sum_{k=1}^{j}\frac{H_{n-j+k}^{2}}{k^{2}}
+ \sum_{k=1}^{j}\frac{H_{n-j+k}^{(2)}}{k^{2}} \\
& + 2 \sum_{k=1}^{j}\frac{H_{n-j+k}^{2}}{k}
- 4 \sum_{k=1}^{j}\frac{H_{n-j+k}}{k^{2}} - 2
\sum_{k=1}^{j}\frac{H_{n-j+k}^{(2)}}{k}
+ 8 \sum_{k=1}^{j}\frac{H_{n-j+k}}{k} \displaybreak[1] \\
& + \sum_{k=1}^{n-j+1}\frac{H_{k} H_{j-1+k}^{2}}{k} -
\sum_{k=1}^{n-j+1}\frac{H_{k}H_{j-1+k}^{(2)}}{k}
+4 \sum_{k=1}^{n-j+1}\frac{H_{k}H_{j-1+k}}{k} -
\sum_{k=1}^{n-j+1}\frac{H_{j-1+k}^{2}}{k^{2}}
+ \sum_{k=1}^{n-j+1}\frac{H_{j-1+k}^{(2)}}{k^{2}} \\
& + 2 \sum_{k=1}^{n-j+1}\frac{H_{j-1+k}^{2}}{k}
- 4 \sum_{k=1}^{n-j+1}\frac{H_{j-1+k}}{k^{2}} - 2
\sum_{k=1}^{n-j+1}\frac{H_{j-1+k}^{(2)}}{k}
+ 8 \sum_{k=1}^{n-j+1}\frac{H_{j-1+k}}{k} \\
& - 2 \sum_{k=1}^{j}\frac{H_{k}^{2}}{k^{2}} - 4
\sum_{k=1}^{j}\frac{H_{k}^{(3)}}{k}
+ 5 \sum_{k=1}^{j}\frac{H_{k}^{(2)}}{k}
- 2 \sum_{k=1}^{n-j+1}\frac{H_{k}^{2}}{k^{2}} - 4
\sum_{k=1}^{n-j+1}\frac{H_{k}^{(3)}}{k}
+ 5 \sum_{k=1}^{n-j+1}\frac{H_{k}^{(2)}}{k} \displaybreak[1] \\
& + \frac{1}{3} H_{j}^{3} - 5 H_{j} H_{j}^{(2)} - \frac{1}{3} H_{j}^{(3)}
- \frac{9}{2} H_{j}^{2} + \frac{12j+9}{2} H_{j}^{(2)} + 2 H_{j}^{(4)} - 22
H_{j} + 4 H_{j} H_{j}^{(3)} \\ & + \frac{1}{3} H_{n-j+1}^{3} - 5 H_{n-j+1}
H_{n-j+1}^{(2)} - \frac{1}{3} H_{n-j+1}^{(3)} - \frac{9}{2} H_{n-j+1}^{2}
+ \frac{12(n-j+1)+9}{2} H_{n-j+1}^{(2)} \\ & + 2 H_{n-j+1}^{(4)} - 22
H_{n-j+1} + 4 H_{n-j+1} H_{n-j+1}^{(3)} \\ & - H_{j}^{2} H_{n} +
H_{j}^{(2)} H_{n} - 4 H_{j} H_{n} - H_{n-j+1}^{2} H_{n} + H_{n-j+1}^{(2)}
H_{n} - 4 H_{n-j+1} H_{n} \\ & + 3 H_{n}^{2} - 6 (n+1) H_{n}^{(2)} + 17
H_{n} + 16. 
\end{align*}

To find an asymptotic equivalent for $b_{n,j}$ for a fixed ratio $j= \rho n$, 
with $0 < \rho < 1$ and $n \to \infty$ we can use the following relations, which
can be obtained by using the asymptotic expansion 
$H_{n-j+k} = \log n + \log(1-\rho) + \gamma + \log\big(1+\frac{k}{(1-\rho)n}\big)
+ \mathcal{O}\left(\frac{1}{n}\right)$ for $1 \le k \le j=\rho n$:
\begin{align*}
   \sum_{k=1}^{j} \frac{H_{k} H_{n-j+k}}{k} & = \frac{1}{2} \log^{3} n + 
   \mathcal{O}\left(\log^{2} n\right) = 
   \sum_{k=1}^{n-j+1} \frac{H_{k} H_{j-1+k}}{k}\, , \\
   \sum_{k=1}^{j} \frac{H_{k} H_{n-j+k}^{2}}{k} & = \frac{1}{2} \log^{4} n + 
   (\log \rho + \log (1-\rho) + 2 \gamma) \log^{3} n + 
   \mathcal{O}\left(\log^{2} n\right) = 
   \sum_{k=1}^{n-j+1} \frac{H_{k} H_{j-1+k}^{2}}{k}\, , \\
   \sum_{k=1}^{j} \frac{H_{n-j+k}^{2}}{k} & = \log^{3} n + 
   \mathcal{O}\left(\log^{2} n\right) = 
   \sum_{k=1}^{n-j+1} \frac{H_{j-1+k}^{2}}{k}\, .
\end{align*}
The assumption that $j= \rho n$
with $0 < \rho < 1$, is a usual one in the analysis of algorithms. 
(For instance in the study of hashing algorithms such a proportionality
factor as $\rho$ is called the {\em load factor}.)
One particular aspect of such an assumption is that one has to consider
only one parameter ($n$) which goes to infinity. Of course, such an
assumption does not necessarily cover instances like $n,j\to\infty$,
as $j\sim\log n$.





The other terms involving sums are bounded by $\mathcal{O}\left(\log^{2} n\right)$.

This leads to the asymptotic expansion for the second factorial moments $b_{n,j}$,
\begin{equation*}
  b_{n,j} = \log^{4} n + \left( 2 \log \rho + 2 \log (1-\rho) + 4 \gamma + \frac{20}{3} 
  \right) \log^{3} n + \mathcal{O}\left(\log^{2} n \right) \, ,
\end{equation*}
for $j = \rho n$ with a fixed $0 < \rho < 1$.
As a consequence we find as the asymptotic equivalent of the variance 
$b_{n,j}+a_{n,j}-a_{n,j}^{2} \sim \frac{2}{3} \log^{3} n$.

\section{Weight of subtrees}

We consider now the total weight of subtrees in binary search trees. 
To do this, we introduce the probability generating function $P_{n,j}(v)$ of
the total weight of the subtree $j$ in a binary search tree with $n$ nodes.
In other words, the coefficient $[v^k]P_{n,j}(v)$ is the probability
that the total weight of the subtree with root $j$ is $k$, assuming that
all permutations of $n$ elements are equally likely. The average 
$a_{n,j}$ is then obtained (as usual) as $a_{n,j}=P'_{n,j}(1)$, 
and the variance as $b_{n,j}+a_{n,j}-\big(a_{n,j}\big)^2$, with
$b_{n,j}=P^{\prime\prime}_{n,j}(1)$.

We get the recurrence
\begin{equation*}
   P_{n,j}(v) = \frac{1}{n}\left(\sum_{k=1}^{j-1}P_{n-k,j-k}(v) +
   \sum_{k=j+1}^{n} P_{k-1,j}(v) + v^{H_{n}} f_{j-1}(v) f_{n-j}(v) \right) \, ,
   \quad P_{1,1}(v) = v \, ,
\end{equation*}
where $f_{n}(v)$ is the probability generating function of the total weight of
a binary search tree, as introduced in section \ref{sect2}. This recursion
is obtained by considering the three cases that the subtree with root $j$
is in the right subtree or left subtree, or that it is the whole tree, in which
case we must add a contribution of $H_n$ to the total count.

For the expectations $a_{n,j} = P_{n,j}'(1)$, we get the following recursion:
\begin{equation*}
   a_{n,j} = \frac{1}{n} \sum_{k=1}^{j-1} a_{n-k,j-k} 
   + \frac{1}{n} \sum_{k=j+1}^{n} a_{k-1,j} + \frac{1}{n}
   (2n + H_{n} - H_{j-1} - H_{n-j} - 2) \, , \quad a_{1,1}=1 \, .
\end{equation*}

An approach with generating functions leads to the result
\begin{equation*}
   a_{n,j} = (n+1) H_{n}^{(2)} - j H_{j}^{(2)} - (n+1-j) H_{n+1-j}^{(2)}
   + 2 H_{j} + 2 H_{n+1-j} + \frac{1}{j} + \frac{1}{n+1-j} - 5 \, ,
\end{equation*}
what means that the weight of the subtree $j$ with $j = \rho n$
and $0 < \rho < 1$ is asymptotically $4 \log n$.

\section{A catalogue of solutions}

In this section we give a list of solutions of
\begin{equation*}
a_{n,j}=\text{\sf{toll-function}}+\frac{1}{n}\bigg[
\sum_{k=1}^{j-1}a_{n-k,j-k}+\sum_{k=j+1}^na_{k-1,j}
\bigg],\quad a_{1,1}=1.
\end{equation*}
 with various toll functions of
the form $n^dH_n$ and (simpler) $n^d$. We write $\j=n+1-j$ for
convenience. (The solutions
 are obtained by an {\it ansatz\/} and matching coefficients.)
Such a catalogue is especially useful for a {\it repertoire approach,} see
\cite{GrKn82}.

For $H_n$ we have


\begin{align*}
a_{n,j}=
\frac{1}{2}\Big(H_{j}^2-H_{j}^{(2)}\Big)
+\frac{1}{2}\Big(H_{\j}^2-H_{\j}^{(2)}\Big)
+2H_j+2H_{\j}
-H_n-2.
\end{align*}


For $nH_n$ we have


\begin{align*}
a_{n,j}&=
(n+1)\Big(H_{n}^2-H_{n}^{(2)}\Big)
-(\tfrac32+j)\Big(H_{j}^2-H_{j}^{(2)}\Big)
-(\tfrac32+\j)\Big(H_{\j}^2-H_{\j}^{(2)}\Big)\\&
+(j-4)H_j+(\j-4)H_{\j}
+(n+4)H_n-n+3.
\end{align*}


For $n^2H_n$ we have


\begin{align*}
a_{n,j}&=
-3(n+1)\Big(H_{n}^2-H_{n}^{(2)}\Big)
+(\tfrac72+3j)\Big(H_{j}^2-H_{j}^{(2)}\Big)
+(\tfrac72+3\j)\Big(H_{\j}^2-H_{\j}^{(2)}\Big)\\&
-(\tfrac{3}{2}j^2+\tfrac{7}{2}j-8)H_j
-(\tfrac{3}{2}\j^2+\tfrac{7}{2}\j-8)H_{\j}
+(n+2)(3n-4)H_n\\&
-\tfrac72j\j-\tfrac{1}{4}n^2+\tfrac{19}{4}n-3.
\end{align*}




For $n^3H_n$ we have


\begin{align*}
a_{n,j}&=
7(n+1)\Big(H_{n}^2-H_{n}^{(2)}\Big)
-(\tfrac{15}2+7j)\Big(H_{j}^2-H_{j}^{(2)}\Big)
-(\tfrac{15}2+7\j)\Big(H_{\j}^2-H_{\j}^{(2)}\Big)\\&
-(\tfrac{2}{3}j^3-5j^2-\tfrac{26}{3}j+16)H_j
-(\tfrac{2}{3}\j^3-5\j^2-\tfrac{26}{3}\j+16)H_{\j}
+(2n^3-6n^2-7n+16)H_n\\&
-(\tfrac76n-\tfrac{61}{6})j\j-\tfrac{1}{9}n^3+\tfrac{3}{2}n^2
-\tfrac{205}{18}n+3.
\end{align*}





For $n^4H_n$ we have


\begin{align*}
a_{n,j}&=
-15(n+1)\Big(H_{n}^2-H_{n}^{(2)}\Big)
+(\tfrac{31}2+15j)\Big(H_{j}^2-H_{j}^{(2)}\Big)
+(\tfrac{31}2+15\j)\Big(H_{\j}^2-H_{\j}^{(2)}\Big)\\&
-(\tfrac{5}{12}j^4-\tfrac{5}{2}j^3+\tfrac{145}{12}j^2+19j-32)H_j
-(\tfrac{5}{12}\j^4-\tfrac{5}{2}\j^3+\tfrac{145}{12}\j^2+19\j-
32)H_{\j}\\&
+(\tfrac{5}{3}n^4-\tfrac{10}{3}n^3+\tfrac{40}{3}n^2
+\tfrac{52}{3}n-32 )H_n\\&
+\tfrac{23}{72}j^2\j^2
-(\tfrac{23}{36}n^2-\tfrac{187}{72}n+\tfrac{214}{9})j\j
-\tfrac{1}{16}n^4+\tfrac {19}{24}n^3
-\tfrac {51}{16}n^2+\tfrac {599}{24}n-3.
\end{align*}



For $n^5H_n$ we have


\begin{align*}
a_{n,j}&=
31(n+1)\Big(H_{n}^2-H_{n}^{(2)}\Big)
-(\tfrac{63}2+31j)\Big(H_{j}^2-H_{j}^{(2)}\Big)
-(\tfrac{63}2+31\j)\Big(H_{\j}^2-H_{\j}^{(2)}\Big)\\&
-(\tfrac{3}{10}j^5-
\tfrac{7}{4}j^4+\tfrac{19}{3}j^3-\tfrac{105}{4}j^2-
\tfrac{1189}{30}j-64)H_j\\&
-(\tfrac{3}{10}\j^5-
\tfrac{7}{4}\j^4+\tfrac{19}{3}\j^3-\tfrac{105}{4}\j^2-
\tfrac{1189}{30}\j-64)H_{\j}\\&
+(\tfrac{3}{2}n^5-\tfrac{5}{2}n^4+\tfrac{15}{2}n^3-
\tfrac{55}{2}n^2
-38n+64)H_n\\&
+(\tfrac{17}{40}n-\tfrac{13}{20})j^2\j^2
-(\tfrac{17}{40}n^3-\tfrac{7}{8}n^2+\tfrac{733}{120}n-
\tfrac{3053}{60})j\j\\&
-\tfrac{1}{25}n^5+\tfrac{41}{80}n^4-\tfrac {497}{360}n^3
+\tfrac {559}{80}n^2-\tfrac {93743}{1800}n+3.
\end{align*}




For $n^6H_n$ we have


\begin{align*}
a_{n,j}&=
-63(n+1)\Big(H_{n}^2-H_{n}^{(2)}\Big)
+(\tfrac{127}2+63j)\Big(H_{j}^2-H_{j}^{(2)}\Big)
+(\tfrac{127}2+63\j)\Big(H_{\j}^2-H_{\j}^{(2)}\Big)\\&
-({\tfrac {7}{30}}{j}^{6}-\tfrac75{j}^{5}+\tfrac{14}{3}{j}^{4}-
14{j}^{3}+{
\tfrac {273}{5}}{j}^{2}+{\tfrac {809}{10}}j-128
)H_j\\&
-({\tfrac {7}{30}}{\j}^{6}-\tfrac75{\j}^{5}+\tfrac{14}{3}{\j}^{4}-
14{\j}^{3}+{
\tfrac {273}{5}}{\j}^{2}+{\tfrac {809}{10}}\j-128
)H_{\j}\\&
+(\tfrac75{n}^{6}-{\tfrac {21}{10}}{n}^{5}+{\tfrac
{35}{6}}{n}^{4}-{
\tfrac {91}{6}}{n}^{3}+{\tfrac {1673}{30}}{n}^{2}+{\tfrac
{1189}{15}}
n-128
)H_n\\&
-\tfrac{47}{450}j^3\j^3
+({\tfrac {47}{100}}{n}^{2}-{\tfrac {301}{600}}n+{\tfrac
{887}{600}})j^2\j^2\\&
-({\tfrac {47}{150}}{n}^{4}-{\tfrac {113}{600}}{n}^{3}+{\tfrac
{491}{
200}}{n}^{2}-{\tfrac {3839}{300}}n+{\tfrac {94553}{900}}
)j\j\\&
-\tfrac1{36}{n}^{6}+{\tfrac {37}{100}}{n}^{5}-{\tfrac
{241}{360}}{n}^{4}+
{\tfrac {31}{10}}{n}^{3}-{\tfrac {5239}{360}}{n}^{2}+{\tfrac
{2657}{
25}}n-3.
\end{align*}


For $n^7H_n$ we have


\begin{align*}
a_{n,j}&=
127(n+1)\Big(H_{n}^2-H_{n}^{(2)}\Big)
-(\tfrac{255}2+127j)\Big(H_{j}^2-H_{j}^{(2)}\Big)
-(\tfrac{255}2+127\j)\Big(H_{\j}^2-H_{\j}^{(2)}\Big)\\&
-({\tfrac {4}{21}}{j}^{7}-\tfrac65{j}^{6}+{\tfrac {59}{15}}{j}^{5}-
\tfrac{21}{2}
{j}^{4}+{\tfrac {88}{3}}{j}^{3}-{\tfrac {1113}{10}}{j}^{2}-{\tfrac
{
5721}{35}}j+256
)H_j\\&
-({\tfrac {4}{21}}{\j}^{7}-\tfrac65{\j}^{6}+{\tfrac
{59}{15}}{\j}^{5}-\tfrac{21}{2}
{j}^{4}+{\tfrac {88}{3}}{\j}^{3}-{\tfrac {1113}{10}}{\j}^{2}-
{\tfrac {
5721}{35}}\j+256
)H_{\j}\\&
+(\tfrac43{n}^{7}-{\tfrac {28}{15}}{n}^{6}+{\tfrac
{77}{15}}{n}^{5}-{
\tfrac {35}{3}}{n}^{4}+{\tfrac {91}{3}}{n}^{3}-{\tfrac
{1687}{15}}{
n}^{2}-{\tfrac {809}{5}}n+256
)H_n\\&
-({\tfrac {31}{126}}n-{\tfrac {347}{3150}}
)j^3\j^3
+({\tfrac {31}{63}}{n}^{3}-{\tfrac {19}{150}}{n}^{2}+{\tfrac
{1237}{900
}}n-{\tfrac {9839}{3150}}
)j^2\j^2\\&
-({\tfrac {31}{126}}{n}^{5}+{\tfrac {509}{3150}}{n}^{4}+{\tfrac
{129}{
100}}{n}^{3}-{\tfrac {30697}{6300}}{n}^{2}+{\tfrac
{164677}{6300}}
n-{\tfrac {112079}{525}}
)j\j\\&
-\tfrac1{49}{n}^{7}+\tfrac27{n}^{6}-{\tfrac
{54}{175}}{n}^{5}+{\tfrac {457}{
280}}{n}^{4}-{\tfrac {389}{60}}{n}^{3}+{\tfrac
{8283}{280}}{n}^{2}
-{\tfrac {3155909}{14700}}n+3.
\end{align*}

And here is the simpler set of solutions:

For $1$,

\begin{align*}
a_{n,j}&=H_j+H_\j-1,
\end{align*}




for $n$,

\begin{align*}
a_{n,j}&=-(3j+2)H_j-(3\j+2)H_\j+2(n+1)H_n+2n+5
\end{align*}





for $n^2$,

\begin{align*}
a_{n,j}&=(6j+7)H_j+(6\j+7)H_\j-6(n+1)H_n+\tfrac32n^2-\tfrac92n-13
-\tfrac{9}{2}n+3j\j-13
\end{align*}




for $n^3$,

\begin{align*}
a_{n,j}&=-(14j+15)H_j-(14\j+15)H_\j+14(n+1)H_n
+(2n-8)j\j+\tfrac43n^3-3n^2+\tfrac{29}3n+29
\end{align*}




for $n^4$,

\begin{align*}
a_{n,j}&=(30j+31)H_j+(30\j+31)H_\j-30(n+1)H_n\\&
+(\tfrac53n^2-\tfrac{25}6n+\tfrac{55}{3}
)j\j+\tfrac54n^4-\tfrac52n^3+\tfrac{25}4n^2-20n-61
-\tfrac56j^2\j^2
\end{align*}






for $n^5$,
\begin{align*}
a_{n,j}&=-(62j+63)H_j-(62\j+63)H_\j+62(n+1)H_n\\&-(\tfrac32n-
2)j^2\j^2
+(\tfrac32n^3-\tfrac52n^2+\tfrac{19}2n-39
)j\j
+\tfrac65n^5-\tfrac94n^4+\tfrac{31}6n^3-
\tfrac{51}4n^2+\tfrac{1219}{30}n+125
\end{align*}

for $n^6$,


\begin{align*}
a_{n,j}&=(126j+127)H_j+(126\j+127)H_\j-126(n+1)H_n\\&
+\tfrac{7}{15}j^3\j^3
-(\tfrac{21}{10}n^2-\tfrac{14}5n+\tfrac{133}{30})j^2\j^2
+(\tfrac75n^4-\tfrac75n^3+\tfrac{91}{15}n^2-
\tfrac{301}{15}n+\tfrac{1204}{15}
)j\j\\&+
\tfrac76{n}^{6}-{\tfrac {21}{10}}{n}^{5}+\tfrac{14}3{n}^{4}-
\tfrac{21}2{n}^{3}+{
\tfrac {77}{3}}{n}^{2}-{\tfrac {819}{10}}n-253,
\end{align*}

for $n^7$,


\begin{align*}
a_{n,j}&=-(254j+255)H_j -(254\j+255)H_\j+254(n+1)H_n\\&
+(\tfrac43n-\tfrac{16}{15})j^3\j^3
-(\tfrac83n^3-\tfrac{14}5n^2+\tfrac{91}{15}n-\tfrac{142}{15}
)j^2\j^2\\&
+(\tfrac43n^5-\tfrac8{15}n^4+\tfrac{21}5n^3-
\tfrac{193}{15}n^2+\tfrac{613}{15}n-\tfrac{814}5
)j\j\\&+
\tfrac87n^7-2n^6+\tfrac{22}5n^5-\tfrac{19}2n^4+21n^3-
\tfrac{103}2n^2+\tfrac{5756}{35}n+509.
\end{align*}
(Several general formul\ae\ for the coefficients can
be spotted from that.)

\section{Results  a la Multiple Quickselect}
Several authors have considered {\it Multiple Quickselect};
see \cite{LeMa96, Prodinger95, PaPr99}. The idea is to use
Quicksort (vs. Quickselect) to find a prescribed set  of
values, e.~g. the third and the tenth; in general we search simultaneously
for $p$ elements with ranks $j_{1}, j_{2}, \dots, 
j_{p}$ where $1 \le j_{1} < j_{2} < \dots < j_{p} \le n$. 
Here we consider analogous recursions with harmonic toll functions:
\begin{align*}
   X[n; j_{1}, \dots, j_{p}] & = H_{n} + \frac{1}{n}\bigg(
   \sum_{t=1}^{p-1} \sum_{k=j_{t}+1}^{j_{t+1}} X[k-1; j_{1}, \dots, j_{t}]
   + \sum_{k=j_{p}+1}^{n} X[k-1; j_{1}, \dots, j_{p}]  \\
   &  +
   \sum_{k=1}^{j_{1}-1} X[n-k; j_{1}-k, \dots, j_{p}-k] +
   \sum_{t=1}^{p-1} \sum_{k=j_{t}}^{j_{t+1}-1} X[n-k; j_{t+1}-k, \dots, j_{p}-k]
   \bigg) \, .
\end{align*}

One finds that the solution of this recursion is given by
\begin{align*}
   X[n; j_{1}, \dots, j_{p}] & = \frac{1}{2}H_{j_{1}}^{2} - \frac{1}{2}H_{j_{1}}^{(2)}
   + 2 H_{j_{1}} + \frac{1}{2} H_{n+1-j_{p}}^{2} - \frac{1}{2} H_{n+1-j_{p}}^{(2)} +
   2 H_{n+1-j_{p}} \\
   & {} + \sum_{t=2}^{p} \left( H_{j_{t}+1-j_{t-1}}^{2}
   -H_{j_{t}+1-j_{t-1}}^{(2)}+2H_{j_{t}+1-j_{t-1}}\right) - H_{n} - 2p \, .
\end{align*}
The proof can be done by induction. The only formul{\ae} which are required, are
\begin{align*}
   & \sum_{k=1}^{n} H_{k}^{2} = (n+1) H_{n}^{2} - (2n+1) H_{n} + 2n \, , \qquad
   \sum_{k=1}^{n} H_{k}^{(2)} = (n+1) H_{n}^{(2)} - H_{n} \, , \\
   & \sum_{k=1}^{n} H_{k} = (n+1) H_{n} - n \, .
\end{align*}
The proof is  omitted here, 
because it is long, but completely routine. It should
be remarked, that for the case $p=2$ the recursion was solved by setting up
generating functions and extracting coefficients of the solution of the resulting
linear differential equation. The solution of the general case for arbitrary $p$ 
was then guessed and verified by computer experiments.

\section{Results  a la Grand Averages}
 For the analysis of Multiple Quickselect, it was suggested that, 
 instead of  a given set of $p$ values
(see \cite{Prodinger95}) it is easier to consider all
$\binom{n}{p}$ subsets of size $p$. This leads to the
{\it Grand Averages\/} of the number of passes and the
number of comparisons. Here again we study such questions together with harmonic
toll functions.

With the notation from the previous section, this means that
we are interested in
\begin{equation*}
   \mathcal{X}_{n,p} = \frac{1}{\binom{n}{p}}\sum_{1 \le j_{1} < j_{2} < \cdots 
                       < j_{p} \le n} X[n;j_{1},\dots,j_{p}] \, .
\end{equation*}

To evaluate this, we use two summation formul{\ae}, which can be found in \cite{Prodinger95}:
\begin{align*}
   \sum_{1 \le j_{1} < j_{2} < \cdots < j_{p} \le n} H_{j_{1}} &= \binom{n+1}{p}
   (H_{n+1}-H_{p}) = \sum_{1 \le j_{1} < j_{2} < \cdots < j_{p} \le n} H_{n+1-j_{p}} \, , \\
   \sum_{t=2}^{p} \sum_{1 \le j_{1} < j_{2} < \cdots < j_{p} \le n} H_{j_{t}+1-j_{t-1}} &=
   (p-1)\binom{n+2}{p} (H_{n+2}-H_{p}) - (p-1) \binom{n}{p-1} \, .
\end{align*}

In a similar manner as shown there, one gets the remaining required formul{\ae}:
\begin{align*}
   \sum_{1 \le j_{1} < j_{2} < \cdots < j_{p} \le n} (H_{j_{1}}^{2}-H_{j_{1}}^{(2)}) 
    &= \binom{n+1}{p} \Big((H_{n+1}-H_{p})^{2}-(H_{n+1}^{(2)}-H_{p}^{(2)})\Big) 
  \\& = \sum_{1 \le j_{1} < j_{2} < \cdots < j_{p} \le n} (H_{n+1-j_{p}}^{2}-H_{n+1-j_{p}}^{(2)}) \, , \\
\end{align*}
 \begin{align*}
   \sum_{t=2}^{p} \sum_{1 \le j_{1} < j_{2} < \cdots < j_{p} \le n} (H_{j_{t}+1-j_{t-1}}^{2}-H_{j_{t}+1-j_{t-1}}^{2}) &=
   (p-1)\binom{n+2}{p} \Big((H_{n+2}-H_{p})^{2}-(H_{n+2}^{(2)}-H_{p}^{(2)})\Big) \, .
\end{align*}

Collecting all contributions and dividing by $\binom{n}{p}$, we get for the grand average:
\begin{align*}
\mathcal{X}_{n,p} &=
\frac{p(n+1)^2}{(n+2-p)(n+1-p)}\Big((H_n-H_p) ^ 2-H^{(2)}_{n}+H^{(2)}_{p}\Big)\\&+\frac{(1+2p)n^2+(1+8p)n+p(7-p)}{(n+2-p)(n+1-p)}H_n
\\&
-\frac{2(n+1)\big((p+1)n+2p+1\big)}{(n+2-p)(n+1-p)}H_p
-\frac{2p(n^2-pn-2)}{(n+2-p)(n+1-p)}.
\end{align*}

\medskip

The formula for $\mathcal{X}_{n,p}$ 
can also be obtained more directly with generating functions.
The starting point for this is the following recursion:

\begin{equation*}
\psi_n(u,v)=v^{H_n}\frac{1+u}n\sum_{j=0}^{n-1}\psi_j(u,v)\psi_{n-1-j}(u,v)
+1-v^{H_n},\quad \psi_0(u,v)=1.
\end{equation*}
Each $[u^p]\psi_n(u,v)\big/\binom np$ is a
probability generating function, namely the coefficient of $v^k$ in it
is the probability that the smallest subtree of a random binary search
tree of size $n$ containing the root and a random sample of $p$ nodes
has total weight $k$. The factor $1+u$ reflects the fact that the root
can be in the sample or not, and the extra terms correspond to the case
of $p=0$; in this case the total weight must be 0. 

If we denote
\begin{equation*}
\sigma_n(u)=\frac{\partial}{\partial v}\psi_n(u,1),
\end{equation*}
then


\begin{align*}
\sigma_n(u)=H_n(1+u)^n+\frac{2(1+u)}n\sum_{j=0}^{n-1}\sigma_j(u)\sigma_
{n-1-j}(u).
\end{align*}
This recursion can again be translated into a differential equation, viz. 

\begin{equation*}
\frac{\partial}{\partial z}G(z,u)=(1+u)S\big(z(1+u)\big)-S(z) +\frac{2(1+u)}{1-z(1+u)}G(z,u),
\qquad G(0,u)=0,
\end{equation*}
again with $S(z)=\frac{1}{(1-z)^2}\log\frac{1}{1-z}+\frac{1}{(1-z)^2}$.
The solution is
\begin{align*}
G(z,u)&=\frac{u(1+u)}{(1-z(1+u))^2}\log^2\frac{1}{1-z}\\
&-\frac{ (1+u )^{2}{z}^{2}-2 (1+u ) (2u+1 )z+2{u}^{2}+4u+1
}{(1-z)(1-z(1+u))^2}\log\frac{1}{1-z}\\&
-\frac{1}{1-z(1+u)}\log\frac{1}{1-z(1+u)}
-\frac{2u(1+u)z}{(1-z(1+u))^2}.
\end{align*}
It is again a routine task to read off the coefficients, and we get for $p\ge1$


\begin{align*}
[u^p]\sigma_n(u)&=p\binom{n+1}{p}\Big[ (H_{n+1}-H_{p})^2-(H^{(2)}_{n+1}-H^{(2)}_{p})\Big]\\
&+(p-1)\binom{n+1}{p-1}\Big[ (H_{n+1}-H_{p-1})^2-(H^{(2)}_{n+1}-H^{(2)}_{p-1})\Big]\\
&+2(p-1)\binom{n+2}{p}\big[ H_{n+2}-H_{p}\big]
+4\binom{n+1}{p}\big[ H_{n+1}-H_{p}\big]\binom npH_n-2\binom{n}{p-1}n.				
\end{align*}
In order to get the desired averages, one must divide that
by $\binom np$, which yields after simplifications

\begin{align*}
\frac{[u^p]\sigma_n(u)}{\binom np}=\mathcal{X}_{n,p}&=
\frac{p(n+1)^2}{(n+2-p)(n+1-p)}\Big((H_n-H_p) ^ 2-H^{(2)}_{n}+H^{(2)}_{p}\Big)\\&+\frac{(1+2p)n^2+(1+8p)n+p(7-p)}{(n+2-p)(n+1-p)}H_n
\\&
-\frac{2(n+1)\big((p+1)n+2p+1\big)}{(n+2-p)(n+1-p)}H_p
-\frac{2p(n^2-pn-2)}{(n+2-p)(n+1-p)}.
\end{align*}
For fixed $p$ (close to the Quickselect instance), this is roughly $p\log^2n$;
for $p\sim n$ (close to the Quicksort instance), this is roughly $2n$, and
for $p\sim \rho n$, with $0<\rho<1$ it is asymptotic to $\big(\frac{\rho}{(1-\rho)^2}\log^2\rho-\frac{2\rho}{(1-\rho)^2}\log\rho-\frac{2\rho}{1-\rho}\big) n$.
The limit of this for $\rho\to1$ is $2n$, which matches the intuition.
\medskip

It is interesting to note the special instance $p=1$:
\begin{align*}
\frac{n+1}{n}(H_n^2-H_n^{(2)})+\frac{n+4}{n}H_n-4.
\end{align*}
Note that this is the same result as we obtained earlier for
\begin{equation*}
\frac1n\sum_{j=1}^na_{n,j},
\end{equation*}
as it should.
The other interesting special case is $p=n$:
\begin{align*}
2n-H_n.
\end{align*}
This is also clear, since there is no random sample involved, and 
only the average total weight of a random binary search tree
has been computed. 

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

\bibitem{Fill96}
J.~Fill.
\newblock On the distribution of binary search trees under the random
  permutation model.
\newblock {\em Random Structures and Algorithms}, 8:1--25, 1996.

\bibitem{GrPr00}
P.~Grabner and H.~Prodinger.
\newblock Sorting algorithms for broadcast communications: {M}athematical
  analysis.
\newblock {\em Theoretical Computer Science, to appear}, 2000.

\bibitem{GrKn82}
D.~H. Greene and D.~E. Knuth.
\newblock {\em Mathematics for the analysis of algorithms}.
\newblock Birkhauser, Boston, second edition, 1982.

\bibitem{KiPr98}
P.~Kirschenhofer and H.~{P}rodinger.
\newblock Comparisons in {H}oare's {F}ind algorithm.
\newblock {\em Combinatorics, Probability, and Computing}, 7:111--120, 1998.

\bibitem{Knuth68}
D.~E. Knuth.
\newblock {\em The Art of Computer Programming}, volume 1: Fundamental
  Algorithms.
\newblock Addison-Wesley, 1968.
\newblock Third edition, 1997.

\bibitem{LeMa96}
J.~Lent and H.~M. Mahmoud.
\newblock Average--case analysis of multiple quickselect: An algorithm for
  finding order statistics.
\newblock {\em Statistics and Probability Letters}, 28:299--310, 1996.

\bibitem{PaPr99}
A.~Panholzer and H.~{P}rodinger.
\newblock A generating functions approach for the analysis of grand averages
  for multiple quickselect.
\newblock {\em Random Structures and Algorithms}, 1999.

\bibitem{Prodinger95}
H.~Prodinger.
\newblock Multiple quickselect --- {H}oare's {F}ind algorithm for several
  elements.
\newblock {\em Information Processing Letters}, 56:123--129, 1995.

\bibitem{SeFl96}
R.~Sedgewick and P.~Flajolet.
\newblock {\em An Introduction to the Analysis of Algorithms}.
\newblock Addison-Wesley, 1996.

\end{thebibliography}



%\bibliography{pro_bib}

\end{document}


