Quicksort Asymptotics^{1}

James Allen Fill^{2}

Department of Mathematical Sciences

The Johns Hopkins University

            jimfill@jhu.edu and http://www.mts.jhu.edu/~fill/

and Svante Janson

Department of Mathematics

Uppsala University

       svante.janson@math.uu.se and http://www.math.uu.se/~svante/

ABSTRACT

The number of comparisons X_{n} used by Quicksort to sort an array of n distinct numbers has mean μ_{n} of order n log n and standard deviation of order n. Using different methods, Régnier and Rösler each showed that the normalized variate Y_{n} := (X_{n}−μ_{n})/n converges in distribution, say to Y; the distribution of Y can be characterized as the unique fixed point with zero mean of a certain distributional transformation.

We provide the first rates of convergence for the distribution of Y_{n} to that of Y, using various metrics. In particular, we establish the bound 2n^{−1/2} in the d_{2}-metric, and the rate O(n^{ε−(1/2)}) for Kolmogorov–Smirnov distance, for any positive ε.

AMS 2000 subject classifications. Primary 68W40; secondary 68P10, 60F05, 60E05.

Key words and phrases. Quicksort, sorting algorithm, asymptotics, limit distribution, rate of convergence, Kolmogorov–Smirnov distance, density, moment generating function, numerical analysis, d_{p}-metric, coupling.

Date. May 30, 2002.

1This is a preprint of an article accepted for publication in Random Structure & Algoritms \begin{matrix}c \\ ◯ \end{matrix} 2002 John Wiley & Sons, Inc.

2Research supported by NSF grants DMS–9803780 and DMS–0104167, and by The Johns Hopkins University’s Acheson J. Duncan Fund for the Advancement of Research in Statistics; carried out in part while author was Visiting Researcher, Theory Group, Microsoft Research.

1 Introduction and summary

This paper provides the first rates of convergence (as n → ∞) for the distribution of the number of comparisons used by the sorting algorithm Quicksort to sort an array of n distinct numbers. Quicksort is the standard sorting procedure in Unix systems, and has been cited [3] as one of the ten algorithms “with the greatest influence on the development and practice of science and engineering in the 20th century.” We begin with a brief review of what is known about the analysis of Quicksort and a summary of our new results.

The Quicksort algorithm for sorting an array of n distinct numbers is extremely simple to describe. If n = 0 or n = 1, there is nothing to do. If n ≥ 2, pick a number uniformly at random from the given array. Compare the other numbers to it to partition the remaining numbers into two subarrays. Then recursively invoke Quicksort on each of the two subarrays.

Let X_{n} denote the (random) number of comparisons required (so that X_{0} = 0). Then X_{n} satisfies the distributional recurrence relation

X_{n} =^{L} X_{U_{n}−1} + X_{n−U_{n}}^{∗} + n − 1,  n ≥ 1,  (1.1)

where =^{L} denotes equality in law (i.e., in distribution), and where, on the right, U_{n} is distributed uniformly on the set {1, …, n}, X_{j}^{∗} =^{L} X_{j}, and

U_{n}; X_{0}, …, X_{n−1}; X_{0}^{∗}, …, X_{n−1}^{∗}

are all independent.

As is well known and quite easily established, for n ≥ 0 we have

μ_{n} := E X_{n} = 2(n + 1)H_{n} − 4n ∼ 2n ln n,

where H_{n} := ∑_{k=1}^{n} k^{−1} is the nth harmonic number and ∼ denotes asymptotic equivalence. It is also routine to compute explicitly the variance of X_{n} (see Exercise 6.2.2-8 in [15]):

Var X_{n} = 7n^{2} − 4(n + 1)^{2}H_{n}^{(2)} − 2(n + 1)H_{n} + 13n = σ^{2}n^{2} − 2n ln n + O(n) (1.2)

where H_{n}^{(2)} := ∑_{k=1}^{n} k^{−2} are the second-order harmonic numbers and, using the standard notation ≐ for approximate equality,

σ^{2} := 7 − \frac{2}{3}π^{2} ≐ 0.4203.  (1.3)

Consider the normalized variate

Y_{n} := (X_{n} − μ_{n})/n, n ≥ 1.

Then (1.1) implies the recursion

Y_{n} =^{L} \frac{U_{n} − 1}{n}Y_{U_{n}−1} + \frac{n − U_{n}}{n} Y_{n−U_{n}}^{∗} + C_{n}(U_{n}),  n ≥ 1,  (1.4)

with Y_{0} arbitrarily defined (since its coefficient is 0), where on the right, as for X_{n},

we have U_{n} ∼ unif{1, …, n} and Y_{j}^{∗} =^{L} Y_{j}, and U_{n}; Y_{1}, …, Y_{n−1}; Y_{1}^{∗}, …, Y_{n−1}^{∗} are all

independent; further,

C_{n}(i) := \frac{n−1}{n} + \frac{1}{n}(μ_{i−1} + μ_{n−i} − μ_{n}), 1 ≤ i ≤ n.  (1.5)

Note that E Y_{n} = 0 = E C_{n}(U_{n}). We will see below that if n → ∞ and i/n → u ∈ [0, 1], then C_{n}(i) → C(u), where

C(u) := 2u ln u + 2(1 − u) ln(1 − u) + 1, u ∈ [0, 1],

with the natural (continuous) interpretation C(u) := 1 for u = 0, 1.

Moreover, Régnier [18] and Rösler [19] showed, using different methods, that Y_{n} → Y in distribution, and Rösler also showed that Y satisfies the distributional identity

Y =^{L} UY + (1 − U)Y^{∗} + C(U)  (1.6)

obtained by formally taking limits in (1.4), where, on the right, U, Y, and Y^{∗} are independent, with Y^{∗} =^{L} Y and U ∼ unif(0, 1). [Rösler [19] showed further that (1.6) characterizes the limiting law L(Y), subject to E Y = 0 and Var Y < ∞. For a complete characterization of the distributions satisfying (1.6), see [7].]

The purpose of the present paper is to study the rate of convergence of L(Y_{n}) to L(Y), using several different measures of the distance between L(Y_{n}) and L(Y).

First, for real 1 ≤ p < ∞, let ∥X∥_{p} := (E |X|^{p})^{1/p} denote the L^{p}-norm, and let d_{p} denote the metric on the space of all probability distributions with finite pth absolute moment defined by

d_{p}(F, G) := min ∥X − Y ∥_{p},

taking the minimum over all pairs of random variables X and Y (defined on the same probability space) with L(X) = F and L(Y) = G. We will use the fact [1] that the minimum is attained for each 1 ≤ p < ∞ by the same X and Y, viz., X := F^{−1}(u) and Y := G^{−1}(u) defined for u in the probability space (0, 1) (with Lebesgue measure).

We will for simplicity write d_{p}(X, Y) := d_{p}(L(X), L(Y)) for random variables X and Y, but note that this distance depends only on the marginal distributions of X and Y.

Rösler [19] showed that d_{p}(Y_{n}, Y) → 0 as n → ∞ for every 1 ≤ p < ∞. In Sections 2 and 3 we will quantify this and show that

d_{p}(Y_{n}, Y) = O (n^{−1/2})

for every fixed p. In the case p = 2 we will further show the explicit bound

d_{2}(Y_{n}, Y) < 2n^{−1/2}.

The best lower bound we can show (Section 4) is

d_{p}(Y_{n}, Y) ≥ c\frac{ln n}{n}, p ≥ 2,

with c > 0 independent of p. We do not know what the correct rate is, even for p = 2. In an earlier draft of this paper, we conjectured the rate n^{−1/2}. Subsequent to that draft, however, Neininger and Rüschendorf [17] surprisingly showed that for another metric, namely, the Zolotarev metric ζ_{3}, the correct rate is (ln n)/n. We therefore now guess that the rate is also (ln n)/n for the metrics d_{p}, matching our lower bound—but proving this is a challenge.

In Section 5 we use these results to bound the Kolmogorov–Smirnov distance d_{KS}(Y_{n}, Y) between L(Y_{n}) and L(Y). We show that

d_{KS}(Y_{n}, Y) = O (n^{ε−(1/2)})

for every ε > 0. Again we do not know the exact rate, but guess that it is (ln n)/n, as for the ζ_{3} metric considered by Neininger and Rüschendorf [17]. The best lower bound we can prove is c/n with c > 0.

In Section 6 we prove a kind of local limit theorem which enables us to approximate the density function f of Y. (It was proved by Tan and Hadjicostas [20] that Y has a density function; f is bounded and infinitely differentiable by [6].)

Rösler [19] showed that (for fixed λ ∈ R) the moment generating function values E e^{λY_{n}} are bounded and thus converge to E e^{λY}. Again we quantify his bounds and give in Section 7 explicit bounds, based on Rösler’s method.

In several (but not all) bounds we give explicit numerical values to constants. These values are hardly the best possible, but we make some effort to get fairly small values. This includes sometimes the use of extensive numerical verifications by computer for small n. [All numerical calculations have been verified independently by the two authors, the (alphabetically) first using Mathematica and the second using Maple.] Such arguments could be simplified or omitted at the cost of increasing the constants.

1.1 Preliminaries

In order to later estimate C_{n}(i) defined by (1.5) we need some explicit bounds on μ_{n}. First, as mentioned above,

μ_{n} = 2(n + 1)H_{n} − 4n,  (1.7)

which can be rewritten

μ_{n} = 2(n + 1)H_{n+1} − 4n − 2.  (1.8)

Next we use the bounds on the harmonic numbers (see, e.g., Section 1.2.11.2 in [14])

ln n + γ ≤ H_{n} ≤ ln n + γ + \frac{1}{2n},  n ≥ 1,  (1.9)

where γ ≐ 0.5772 is Euler’s constant. Hence, for n ≥ 1, from (1.7)

2(n + 1) ln n + (2γ − 4)n + 2γ ≤ μ_{n} ≤ 2(n + 1) ln n + (2γ − 4)n + 2γ + \frac{n+1}{n} (1.10)

and from (1.8)

2n ln n + (2γ − 4)n + 2 ≤ μ_{n−1} ≤ 2n ln n + (2γ − 4)n + 3.  (1.11)

2 Bounding d_{2}(Y_{n}, Y)

In this section we prove the following explicit estimate for d_{2}(Y_{n}, Y).

Theorem 2.1. For all n ≥ 1, d_{2}(Y_{n}, Y) < 2/\sqrt{n}.

Proof. We basically follow the method of Rösler [19], making all estimates explicit. We study in this paper only properties of the univariate distributions of Y_{n}. We thus take the liberty of letting Y_{n} denote any random variable with the appropriate distribution [Y_{n} =^{L} (X_{n} − μ_{n})/n]. We then may choose Y_{0}, Y_{1}, … defined on the same probability space as Y and such that

∥Y_{n} − Y ∥_{2} = d_{2}(Y_{n}, Y), n ≥ 0.

Further, let (Y^{∗}, Y_{0}^{∗}, Y_{1}^{∗}, …) be an independent copy of (Y, Y_{0}, Y_{1}, …) and let U ∼ unif(0, 1) be independent of everything else. For convenience we write a_{n} := d_{2}(Y_{n}, Y).

Observe, by (1.4), that

Y_{n} =^{L} Ỹ_{n} := \frac{⌈nU⌉ − 1}{n}Y_{⌈nU⌉−1} + \frac{n − ⌈nU⌉}{n}Y_{n−⌈nU⌉}^{∗} + C_{n}(⌈nU⌉)  (2.1)

and recall from (1.6) that

Y =^{L} Ỹ := UY + (1 − U)Y^{∗} + C(U).  (2.2)

Therefore,

a_{n}^{2} = d_{2}^{2}(Y_{n}, Y) ≤ E |Ỹ_{n} − Ỹ |^{2}.  (2.3)

Now

Ỹ_{n} − Ỹ = (\frac{⌈nU⌉ − 1}{n}Y_{⌈nU⌉−1} − UY ) + (\frac{n − ⌈nU⌉}{n}Y_{n−⌈nU⌉}^{∗} − (1 − U)Y^{∗})
+ (C_{n}(⌈nU⌉) − C(U))
=: W_{1} + W_{2} + W_{3},

say. Given U, the random variables W_{1} and W_{2} are independent with zero mean, while W_{3} is a constant. Hence

E (|Ỹ_{n} − Ỹ |^{2} | U) = E ((W_{1} + W_{2} + W_{3})^{2}| U) = E (W_{1}^{2} | U) + E (W_{2}^{2} | U) + W_{3}^{2}

and thus, taking expectations,

E |Ỹ_{n} − Ỹ |^{2} = E W_{1}^{2} + E W_{2}^{2} + E W_{3}^{2}.  (2.4)

By symmetry (replacing U by 1 − U), E W_{2}^{2} = E W_{1}^{2}. We estimate this term by conditioning on U, using the independence of U and Y, Y_{0}, . . .. If U = (k + v)/n, with k ∈ \{0, 1, …, n − 1\} and 0 < v ≤ 1, then ⌈nU⌉ = k + 1 and W_{1} = \frac{k}{n}(Y_{k} − Y) − \frac{v}{n}Y; hence Minkowski’s inequality yields

E (W_{1}^{2} | U = (k + v)/n)^{1/2} ≤ \frac{k}{n}∥Y_{k} − Y ∥_{2} + \frac{v}{n}∥Y ∥_{2}
= \frac{k}{n}a_{k} + \frac{v}{n}σ.

Consequently,

E W_{1}^{2} = \frac{1}{n} ∑_{k=0}^{n−1} ∫_{0}^{1}E (W_{1}^{2} | U = (k + v)/n) dv ≤ \frac{1}{n} ∑_{k=0}^{n−1} ∫_{0}^{1}(\frac{k}{n}a_{k} + \frac{v}{n}σ)^{2} dv

= \frac{1}{n} ∑_{k=0}^{n−1} ∫_{0}^{1}(\frac{k^{2}}{n^{2}} a_{k}^{2} + \frac{2k}{n^{2}} va_{k}σ + \frac{v^{2}}{n^{2}} σ^{2}) dv

= \frac{1}{n} ∑_{k=0}^{n−1} (\frac{k^{2}}{n^{2}} a_{k}^{2} + \frac{k}{n^{2}} a_{k}σ + \frac{σ^{2}}{3n^{2}} ) .  (2.5)

We postpone the estimation of E W_{3}^{2}, and introduce the notation

b_{n} := ∥W_{3}∥_{2} = ∥C_{n}(⌈nU⌉) − C(U)∥_{2}.  (2.6)

Combining (2.3)–(2.6), we obtain our fundamental recursive estimate

a_{n}^{2} ≤ 2E W_{1}^{2} + E W_{3}^{2}
≤ \frac{2}{n^{3}} ∑_{k=1}^{n−1} k^{2}a_{k}^{2} + \frac{2σ}{n^{3}} ∑_{k=1}^{n−1} ka_{k} + \frac{2σ^{2}}{3n^{2}} + b_{n}^{2},  n ≥ 1.  (2.7)

We unwrap this recursion partly, by concentrating on the first sum on the right-hand side and regarding the second as known. Thus, writing

y_{n} := \frac{2σ}{n^{3}} ∑_{k=1}^{n−1} ka_{k} + \frac{2σ^{2}}{3n^{2}} + b_{n}^{2},  (2.8)

we define recursively

x_{n} := \frac{2}{n} ∑_{k=1}^{n−1} x_{k} + n^{2}y_{n},  n ≥ 1,  (2.9)

and find by (2.7) and induction

n^{2}a_{n}^{2} ≤ x_{n}, n ≥ 1.

Now, the recursion (2.9) is easily solved (see, e.g., [5]), giving

a_{n}^{2} ≤ n^{−2}x_{n} = y_{n} + 2\frac{n + 1}{n^{2}} ∑_{j=1}^{n−1} \frac{j^{2}}{(j + 1)(j + 2)}y_{j},  n ≥ 1.  (2.10)

We substitute (2.8), treating the three terms separately, into (2.10). The first term in (2.8) yields the sum

∑_{j=1}^{n−1} \frac{j^{2}}{(j + 1)(j + 2)} \frac{2σ}{j^{3}} ∑_{k=1}^{j−1} ka_{k} = ∑_{k=1}^{n−1} ∑_{k<j<n} σka_{k} \frac{2}{j(j + 1)(j + 2)}
= ∑_{k=1}^{n−1} σka_{k} ∑_{j=k+1}^{n−1} (\frac{1}{j(j + 1)} − \frac{1}{(j + 1)(j + 2)})
= ∑_{k=1}^{n−1} σka_{k} (\frac{1}{(k + 1)(k + 2)} − \frac{1}{n(n + 1)})

and the total contribution

\frac{2σ}{n^{3}} ∑_{k=1}^{n−1} ka_{k} + 2\frac{n + 1}{n^{2}} ∑_{k=1}^{n−1} σka_{k} (\frac{1}{(k + 1)(k + 2)} − \frac{1}{n(n + 1)})

= 2σ \frac{n + 1}{n^{2}} ∑_{k=1}^{n−1} \frac{ka_{k}}{(k + 1)(k + 2)}.  (2.11)

The second term in (2.8) yields the sum

∑_{j=1}^{n−1} \frac{j^{2}}{(j + 1)(j + 2)} \frac{2σ^{2}}{3j^{2}} = \frac{2σ^{2}}{3} ∑_{j=1}^{n−1} (\frac{1}{j + 1} − \frac{1}{j + 2}) = \frac{2σ^{2}}{3} (\frac{1}{2} − \frac{1}{n + 1})

and the total contribution

\frac{2σ^{2}}{3n^{2}} + 2\frac{n + 1}{n^{2}} \frac{2σ^{2}}{3} (\frac{1}{2} − \frac{1}{n + 1}) = \frac{2σ^{2}}{3n^{2}} (1 + n + 1 − 2) = \frac{2σ^{2}}{3n} .  (2.12)

Hence we find from (2.10)

a_{n}^{2} ≤ 2σ \frac{n + 1}{n^{2}} ∑_{k=1}^{n−1} \frac{ka_{k}}{(k + 1)(k + 2)} + \frac{2σ^{2}}{3n} + b_{n}^{2} + 2\frac{n + 1}{n^{2}} ∑_{k=1}^{n−1} \frac{k^{2}b_{k}^{2}}{(k + 1)(k + 2)}. (2.13)

We next use the following estimate of b_{n}, whose proof we postpone. Lemma 2.2. For n ≥ 1,

b_{n} := ∥C_{n}(⌈nU⌉) − C(U)∥_{2} ≤ (3 + \frac{2π}{\sqrt{3}}) \frac{1}{n} < \frac{6.63}{n}.

Using this lemma in (2.13), we find in analogy with (2.12)

b_{n}^{2} + 2\frac{n + 1}{n^{2}} ∑_{k=1}^{n−1} \frac{k^{2}b_{k}^{2}}{(k + 1)(k + 2)} < \frac{(6.63)^{2}}{n} < \frac{44}{n}  (2.14)

and thus

a_{n}^{2} ≤ 2σ \frac{n + 1}{n^{2}} ∑_{k=1}^{n−1} \frac{ka_{k}}{(k + 1)(k + 2)} + (44 + \frac{2σ^{2}}{3} ) \frac{1}{n},  n ≥ 1.  (2.15)

We claim that (2.15) implies the sought estimate a_{n} = O(n^{−1/2}). Indeed, assume that n ≥ 1 and that A > 0 is a number such that

a_{k} ≤ A/\sqrt{k}  (2.16)

for 1 ≤ k ≤ n − 1. Then, using k + 1 ≥ [k(k + 2)]^{1/2},

∑_{k=1}^{n−1} \frac{ka_{k}}{(k + 1)(k + 2)} ≤ A ∑_{k=1}^{n−1} \frac{k^{1/2}}{(k + 1)(k + 2)} ≤ A ∑_{k=1}^{n−1} \frac{1}{(k + 2)^{3/2}}

≤ A ∫_{0}^{n−1}\frac{dx}{(x + 2)^{3/2}} = 2A (2^{−1/2} − (n + 1)^{−1/2}) . (2.17)

In particular, for n ≥ 2,

\frac{1}{n} ∑_{k=1}^{n−1} \frac{ka_{k}}{(k + 1)(k + 2)} ≤ \frac{1}{n}2A ≤ 2A(n + 1)^{−1/2}  (2.18)

and thus (2.17) yields (trivially for n = 1, too)

\frac{n + 1}{n} ∑_{k=1}^{n−1} \frac{ka_{k}}{(k + 1)(k + 2)} ≤ 2^{1/2}A.

Consequently, by (2.15),

na_{n}^{2} ≤ 2^{3/2}σA + 44 + 2\frac{σ^{2}}{3} ≤ 2^{3/2}σA + 45.

If 2^{3/2}σA + 45 ≤ A^{2}, which holds for example for A = 8, then this yields na_{n}^{2} ≤ A^{2}, and thus (2.16) holds for k = n, too. By induction, (2.16) holds for all k ≥ 1, and we have proved the explicit estimate

a_{n} ≤ \frac{8}{\sqrt{n}},  n ≥ 1.  (2.19)

This is the desired estimate, apart from the value of the constant. To improve the constant, we use numerical calculations by computer. Indeed, for (2.6),

b_{n}^{2} = ∑_{i=1}^{n} ∫_{(i−1)/n}^{i/n}(C(u) − C_{n}(i))^{2} du

= ∫_{0}^{1}C(u)^{2} du − 2 ∑_{i=1}^{n} C_{n}(i) ∫_{(i−1)/n}^{i/n}C(u) du + ∑_{i=1}^{n} \frac{1}{n}C_{n}(i)^{2}

= \frac{σ^{2}}{3} − 2 ∑_{i=1}^{n} C_{n}(i) (F (\frac{i}{n}) − F (\frac{i − 1}{n})) + \frac{1}{n} ∑_{i=1}^{n} C_{n}(i)^{2},

where F (u) := u^{2} ln u − (1 − u)^{2} ln(1 − u) and C_{n}(i) is given by (1.5); so, given any integer N, b_{n} can be computed exactly for n ≤ N. Next, for n = 1, …, N, an upper bound ā_{n} to a_{n} can be computed recursively from (2.7) or, equivalently, (2.13), using the already computed ā_{k}, k < n, to bound a_{k} in the right-hand side. (We do not know how to compute a_{n} exactly even for n = 3.) For larger n, we use the estimates (2.16) and Lemma 2.2.

Let

V_{n} := ∑_{k=1}^{n} \frac{ka_{k}}{(k + 1)(k + 2)},

V̄_{n} := ∑_{k=1}^{n} \frac{kā_{k}}{(k + 1)(k + 2)},

W_{n} := ∑_{k=1}^{n} \frac{k^{2}b_{k}^{2}}{(k + 1)(k + 2)}.

Then for n > N, arguing as in (2.17), for any A such that (2.16) holds for all k,

V_{n−1} ≤ V̄_{N} + ∑_{k=N+1}^{n−1} \frac{A}{(k + 2)^{3/2}} ≤ V̄_{N} + 2A ((N + 2)^{−1/2} − (n + 1)^{−1/2})

and thus by (2.18)

\frac{n+1}{n}V_{n−1} = V_{n−1} + \frac{1}{n}V_{n−1} ≤ V̄_{N} + 2A(N + 2)^{−1/2}.  (2.20)

Similarly, with B := (3 + \frac{2π}{\sqrt{3}})^{2} < 44, for n > N, by Lemma 2.2, we have

\frac{n+1}{n}W_{n−1} ≤ W_{N} + ∑_{k=N+1}^{n−1} \frac{B}{(k + 1)(k + 2)} + \frac{1}{n} ∑_{k=1}^{n−1} \frac{B}{(k + 1)(k + 2)}

= W_{N} + B (\frac{1}{N + 2} − \frac{1}{n + 1} + \frac{1}{2n} − \frac{1}{n(n + 1)}) = W_{N} + \frac{B}{N + 2} − \frac{B}{2n}.

Consequently, (2.13) yields, using Lemma 2.2 again and (2.20),

a_{n}^{2} ≤ 2σ \frac{n + 1}{n^{2}}V_{n−1} + \frac{2σ^{2}}{3n} + \frac{B}{n^{2}} + 2\frac{n + 1}{n^{2}}W_{n−1}
≤ \frac{1}{n} (2σV̄_{N} + 4σA(N + 2)^{−1/2} + \frac{2σ^{2}}{3} + 2W_{N} + 2B(N + 2)^{−1}) ,  n > N.

In other words, (2.16) holds for k > N, with A replaced by

A_{N} := (2σV̄_{N} + 4σA(N + 2)^{−1/2} + \frac{2σ^{2}}{3} + 2W_{N} + 2B(N + 2)^{−1})^{1/2} .  (2.21)

For N = 100 we find (using Mathematica or Maple), rounded to four decimal places, V̄_{100} ≐ 1.1995 and W_{100} ≐ 0.3466, and thus, taking A = 8 as in (2.19), A_{100} ≐ 2.3332. Moreover, the computer verifies that n^{1/2}ā_{n} < 1.7 for n ≤ 100; thus (2.16) holds for all k ≥ 1 with A = 2.34. Using this value in (2.21) we find A_{100} ≐ 1.9976, and the theorem is proved.

Remark 2.3. The sequence n^{1/2}ā_{n} seems to increase slowly. For n = 100 the value is (rounded to four decimal places) 1.6018, and hence the bound in Theorem 2.1 cannot be much improved using the present method based on (2.7).

It remains to prove Lemma 2.2 above.

Proof of Lemma 2.2. Let I_{i} := {u : ⌈nu⌉ = i} = ((i − 1)/n, i/n]. Thus I_{1}, …, I_{n} form a partition of (0, 1]. We choose a point t_{i} ∈ Ī_{i} for each i (where the bar here indicates closure) and define

C̃_{n}(u) := C(t_{⌈nu⌉}),

i.e., C̃(u) = C(t_{i}) when u ∈ I_{i}. By Minkowski’s inequality,

b_{n} ≤ ∥C_{n}(⌈nU⌉) − C̃_{n}(U)∥_{2} + ∥C̃_{n}(U) − C(U)∥_{2}.  (2.22)

To estimate the second term in (2.22), note that for u ∈ I_{i},

|C̃_{n}(u) − C(u)| = |C(t_{i}) − C(u)| ≤ ∫_{I_{i}}|C′(x)| dx.

The Cauchy–Schwarz inequality yields

|C̃_{n}(u) − C(u)|^{2} ≤ \frac{1}{n} ∫_{I_{i}}|C′(x)|^{2} dx, u ∈ I_{i},

and thus (for any choice of t_{i} ∈ Ī_{i}),

∥C̃_{n}(U) − C(U)∥_{2}^{2} = ∑_{i=1}^{n} ∫_{I_{i}}|C̃_{n}(u) − C(u)|^{2} du

≤ ∑_{i=1}^{n} \frac{1}{n^{2}} ∫_{I_{i}}|C′(x)|^{2} dx

= \frac{1}{n^{2}} ∫_{0}^{1} |C′(x)|^{2} dx.  (2.23)

We have

C′(x) = 2 ln x − 2 ln(1 − x)

and find

∫_{0}^{1}[ln(1 − x)]^{2} dx = ∫_{0}^{1}(ln x)^{2} dx = ∫_{0}^{∞}y^{2}e^{−y} dy = 2

and

∫_{0}^{1}[ln x][ln(1 − x)] dx = ∑_{k=1}^{∞} \frac{1}{k} ∫_{0}^{1}x^{k} | ln x| dx = ∑_{k=1}^{∞} \frac{1}{k} ∫_{0}^{∞}e^{−ky}ye^{−y} dy

= ∑_{k=1}^{∞} \frac{1}{k} \frac{1}{(k + 1)^{2}} = ∑_{k=1}^{∞} (\frac{1}{k(k + 1)} − \frac{1}{(k + 1)^{2}} )

= 1 − (\frac{π^{2}}{6} − 1) = 2 − \frac{π^{2}}{6} ;

consequently,

∫_{0}^{1}|C′(x)|^{2} dx = 8 ∫_{0}^{1}(ln x)^{2} dx − 8 ∫_{0}^{1}[ln x][ln(1 − x)] dx = \frac{4π^{2}}{3} .

Hence (2.23) yields

∥C̃_{n}(U) − C(U)∥_{2} ≤ \frac{1}{n}∥C′(U)∥_{2} ≤ \frac{2π}{\sqrt{3}n}.  (2.24)

For the first term in (2.22), let us first assume that n ≥ 2. For u ∈ I_{i} we have

C_{n}(⌈nu⌉) − C̃_{n}(u) = C_{n}(i) − C(t_{i})
= −\frac{1}{n} + \frac{1}{n}(μ_{i−1} + μ_{n−i} − μ_{n}) − 2t_{i} ln t_{i} − 2(1 − t_{i}) ln(1 − t_{i}).

For i ≤ ⌈n/2⌉ we choose t_{i} = i/n. This yields, using (1.11) and (1.10),

C_{n}(i) − C(t_{i}) ≤ \frac{1}{n} [−1 + 2i ln i + (2γ − 4)i + 3
+2(n − i + 1) ln(n − i) + (2γ − 4)(n − i) + 2γ + 1 + \frac{1}{n−i}

−2(n + 1) ln n − (2γ − 4)n − 2γ − 2i ln(\frac{i}{n}) − 2(n − i) ln(\frac{n−i}{n})]

= \frac{1}{n} [2 ln(\frac{n−i}{n}) + 3 + \frac{1}{n−i}] ≤ \frac{1}{n} [3 − \frac{2i}{n} + \frac{1}{n−i}]
≤ \frac{3}{n}.  (2.25)

In the opposite direction, by (1.11) and (1.10), still for i ≤ ⌈n/2⌉,

C_{n}(i) − C(t_{i}) ≥ \frac{1}{n} [−1 + 2i ln i + (2γ − 4)i + 2
+2(n − i + 1) ln(n − i) + (2γ − 4)(n − i) + 2γ
−2(n + 1) ln n − (2γ − 4)n − 2γ − 1 − \frac{1}{n}
−2i ln(\frac{i}{n}) − 2(n − i) ln(\frac{n−i}{n})]
= \frac{1}{n} [2 ln(\frac{n−i}{n}) − \frac{1}{n}] ≥ \frac{1}{n} [2 ln(\frac{1}{3}) − \frac{1}{2}] ≥ −\frac{3}{n}.

Consequently, for i ≤ ⌈n/2⌉,

|C_{n}(i) − C(t_{i})| ≤ 3/n.  (2.26)

For i > ⌈n/2⌉ we choose instead t_{i} = (i − 1)/n = 1 − t_{n+1−i}. The symmetries of C_{n} and C then yield C_{n}(i) − C(t_{i}) = C_{n}(n + 1 − i) − C(t_{n+1−i}), and since n + 1 − i ≤ n/2, (2.26) shows that |C_{n}(i)−C(t_{i})| ≤ 3/n for i > ⌈n/2⌉, too, i.e., (2.26) holds for all i ≤ n. In other words,

|C_{n}(⌈nu⌉) − C̃_{n}(u)| = |C_{n}(⌈nu⌉) − C(t_{⌈nu⌉})| ≤ 3/n

for all u ∈ (0, 1]; in particular, ∥C_{n}(⌈nU⌉) − C̃_{n}(U)∥_{2} ≤ 3/n for all n ≥ 2. This holds trivially for n = 1, too, for any choice of t_{1}, and together with (2.22) and (2.24) yields the result.

Remark 2.4. Define

c^{∗} := sup{n^{1/2}d_{2}(Y_{n}, Y) : n ≥ 1},

so that, by Theorem 2.1, c^{∗} ≤ 2. Conversely,

c^{∗} ≥ 2^{1/2}d_{2}(Y_{2}, Y) = 2^{1/2}∥Y ∥_{2} = σ\sqrt{2} > 0.9168;

thus the constant 2 in Theorem 2.1 is no more than about twice the optimal value.

Although we do not know the exact value of d_{2}(Y_{n}, Y) for any n > 2, one can in principle for any n and m compute the exact distributions of Y_{n} and Y_{m} and thus d_{2}(Y_{n}, Y_{m}). We have done this for some m, n ≤ 50 using Mathematica and Maple. The results are consistent with a decay of the type d_{2}(Y_{n}, Y) ∼ cn^{−1/2} with c ≈ 1, but our data are too few to be conclusive.

3 Bounding d_{p}(Y_{n}, Y)

In this section we extend Theorem 2.1 and show that d_{p}(Y_{n}, Y) = O(n^{−1/2}) for every p. In contrast to the style of Section 2, we will make no attempt to keep constants small, nor to keep track of them explicitly.

Theorem 3.1. For every p ≥ 1, there exists a constant c_{p} < ∞ such that

d_{p}(Y_{n}, Y) ≤ c_{p}/\sqrt{n}, n ≥ 1.

Proof. Since d_{p} ≤ d_{q} when p ≤ q, it suffices to consider integer p ≥ 2. The case p = 2 is Theorem 2.1 (with c_{2} = 2), so we assume further that p ≥ 3. We use induction on p and assume that the result holds for smaller positive integer values of p.

Let Y, Y_{n}, Y^{∗}, Y_{n}^{∗}, U be as in Section 2, and note that for every p ≥ 1,

∥Y_{n} − Y ∥_{p} = ∥Y_{n}^{∗} − Y^{∗}∥_{p} = d_{p}(Y_{n}, Y),  n ≥ 0,  (3.1)

by the fact [1] that there is an optimal coupling for d_{2} that is optimal for every d_{p}. Using the notation of Section 2, we have, for n ≥ 1,

d_{p}(Y_{n}, Y) ≤ ∥Ỹ_{n} − Ỹ ∥_{p} = ∥W_{1} + W_{2} + W_{3}∥_{p}.  (3.2)

We use a simple lemma to estimate this.

Lemma 3.2. Let Z_{1}, Z_{2}, and Z_{3} be three independent random variables, and let p ≥ 2 be an integer. Then

E |Z_{1} + Z_{2} + Z_{3}|^{p} ≤ E |Z_{1}|^{p} + E |Z_{2}|^{p} + (∥Z_{1}∥_{p−1} + ∥Z_{2}∥_{p−1} + ∥Z_{3}∥_{p})^{p} .

Proof. By the binomial theorem and independence,

E |Z_{1} +Z_{2} +Z_{3}|^{p} ≤ E (|Z_{1}| + |Z_{2}| + |Z_{3}|)^{p} = ∑_{j,k,l} \binom{p}{j, k, l} (E |Z_{1}|^{j}) (E |Z_{2}|^{k}) (E |Z_{3}|^{l}) .

If j ≤ p − 1 and k ≤ p − 1 we estimate E |Z_{1}|^{j} = ∥Z_{1}∥_{j}^{j} ≤ ∥Z_{1}∥_{p−1}^{j} (which holds also for j = 0, disregarding the central expression) and similarly E |Z_{2}|^{k} ≤ ∥Z_{2}∥_{p−1}^{k} and E |Z_{3}|^{l} ≤ ∥Z_{3}∥_{p}^{l}. Hence all terms in the sum, except E |Z_{1}|^{p} and E |Z_{2}|^{p}, are bounded by the corresponding terms in the trinomial expansion of (∥Z_{1}∥_{p−1} + ∥Z_{2}∥_{p−1} + ∥Z_{3}∥_{p})^{p}.

Conditional on U = u, the three variables W_{1}, W_{2}, and W_{3} are independent, so the lemma is applicable. Fix u ∈ (0, 1) and let i = ⌈nu⌉, so 1 ≤ i ≤ n. Then, given U = u, W_{1} = \frac{i−1}{n}Y_{i−1} − uY and thus, for any q ≥ 1,

E (|W_{1}|^{q} | U = u)^{1/q} = ∥\frac{i−1}{n}Y_{i−1} − uY ∥_{q}
≤ ∥\frac{i−1}{n}(Y_{i−1} − Y)∥_{q} + |\frac{i−1}{n} − u| ∥Y ∥_{q}
≤ \frac{i−1}{n}d_{q}(Y_{i−1}, Y) + \frac{1}{n}∥Y ∥_{q}.  (3.3)

Similarly,

E (|W_{2}|^{q} | U = u)^{1/q} ≤ \frac{n−i}{n}d_{q}(Y_{n−i}, Y) + \frac{1}{n}∥Y ∥_{q}.  (3.4)

Further, given U = u, W_{3} = C_{n}(i) − C(u) is a constant, for which we use the simple estimate (from Proposition 3.2 in [19])

|W_{3}| = |C_{n}(⌈nu⌉) − C(u)| ≤ \frac{6}{n} ln n + O (n^{−1}) = O (n^{−1/2}) .  (3.5)

We first use (3.3) with q = p−1 together with the induction hypothesis d_{p−1}(Y_{i−1}, Y) ≤ c_{p−1}(i − 1)^{−1/2}, i ≥ 2, to obtain (also for i = 1)

E (|W_{1}|^{p−1} | U = u)^{1/(p−1)} ≤ c_{p−1} \frac{(i − 1)^{1/2}}{n} + \frac{1}{n}∥Y ∥_{p−1} ≤ b_{1}n^{−1/2},

where b_{1}, like b_{2}, b_{3}, b_{4} below, denotes some constant depending on p only. By similar argument using (3.4) and (3.5), we obtain

E (|W_{1}|^{p−1} | U = u)^{1/(p−1)}+E (|W_{2}|^{p−1} | U = u)^{1/(p−1)}+E (|W_{3}|^{p} | U = u)^{1/p} ≤ b_{2}n^{−1/2}.

Hence, using (3.3) and (3.4) for q = p, too, Lemma 3.2 yields

E (|W_{1} + W_{2} + W_{3}|^{p} | U = u) ≤ (\frac{i − 1}{n}d_{p}(Y_{i−1}, Y) + b_{3} \frac{1}{n})^{p}

+ (\frac{n − i}{n}d_{p}(Y_{n−i}, Y) + b_{3} \frac{1}{n})^{p} + b_{2}^{p}n^{−p/2}.

Taking the average over all u ∈ (0, 1) we finally find the recursive estimate

d_{p}(Y_{n}, Y)^{p} ≤ E |W_{1} + W_{2} + W_{3}|^{p} = E E (|W_{1} + W_{2} + W_{3}|^{p} | U)
≤ \frac{2}{n} ∑_{j=0}^{n−1} (\frac{j}{n}d_{p}(Y_{j}, Y) + b_{3} \frac{1}{n})^{p} + b_{2}^{p}n^{−p/2}.  (3.6)

The proof is now completed by another induction, this one on n. Suppose that d_{p}(Y_{j}, Y) ≤ cj^{−1/2} for 1 ≤ j ≤ n − 1. Then (3.6) yields

d_{p}(Y_{n}, Y)^{p} ≤ \frac{2}{n} ∑_{j=1}^{n−1} (cj^{1/2}n^{−1} + b_{3}n^{−1})^{p} + \frac{2}{n}b_{3}^{p}n^{−p} + b_{2}^{p}n^{−p/2}

≤ \frac{2}{n}(c + b_{3})^{p} ∑_{j=1}^{n−1} j^{p/2}n^{−p} + b_{4}n^{−p/2}

≤ 2(c + b_{3})^{p} ∫_{0}^{1}x^{p/2}n^{−p/2} dx + b_{4}n^{−p/2}
= [2(c + b_{3})^{p} \frac{1}{(p/2) + 1} + b_{4}] n^{−p/2}.  (3.7)

Since p ≥ 3, we have \frac{2}{(p/2)+1} = \frac{4}{p+2} < 1, and thus, if c is sufficiently large,

\frac{4}{p + 2}(c + b_{3})^{p} + b_{4} ≤ c^{p}.

For such c, (3.7) yields d_{p}(Y_{n}, Y)^{p} ≤ (cn^{−1/2})^{p}, which completes both inductions and the proof.

Note that the arguments used above for p ≥ 3 do not work for p = 2, so we need both the proof here and the proof in Section 2.

4 Lower bounds for d_{p}(Y_{n}, Y)

We do not know whether the upper bounds O(n^{−1/2}) proved in the preceding two sections are sharp. We give in this section two simple lower bounds.

First, d_{p}(Y_{n}, Y) = Ω(n^{−1}) for every p by the following general result.

Proposition 4.1. Let W, W_{1}, W_{2}, … be random variables such that W has an absolutely continuous distribution while, for each n ≥ 1 and some constant b_{n}, n(W_{n} − b_{n}) is integer-valued. Then, for each 1 ≤ p < ∞, d_{p}(W_{n}, W) = Ω(1/n). More precisely,

lim inf_{n→∞} n d_{p}(W_{n}, W) ≥ \frac{1}{2}(p + 1)^{−1/p}.  (4.1)

Proof. Let V_{n} := {n(W − b_{n})}, where {x} := x − ⌊x⌋ denotes the fractional part of x. For any coupling of W and W_{n},

|W − W_{n}| = \frac{1}{n}|n(W − b_{n}) − n(W_{n} − b_{n})| ≥ \frac{1}{n}h(V_{n}),

where h(x) := min(x, 1 − x), 0 ≤ x ≤ 1, and thus

d_{p}(W, W_{n}) ≥ \frac{1}{n}∥h(V_{n})∥_{p}.  (4.2)

We regard V_{n} as a random variable taking values in R/Z ≅ T, and find that its distribution, ν_{n} say, has Fourier coefficients

ν̂_{n}(k) = E (e^{−2πikV_{n}}) = e^{2πiknb_{n}}ϕ(−2πkn),

where ϕ is the characteristic function of W. In particular, |ν̂_{n}(k)| = |ϕ(2πkn)|. By our hypothesis on W and the Riemann–Lebesgue lemma, ϕ(x) → 0 as x → ±∞. Thus, for each fixed k ≠ 0, ν̂_{n}(k) → 0 as n → ∞. This implies that ν_{n} converges weakly (as measures on T) to the uniform distribution, i.e., V_{n} →^{L} U where U ∼ unif(0, 1). Consequently, as n → ∞,

∥h(V_{n})∥_{p}^{p} = E h(V_{n})^{p} → E h(U)^{p} = 2 ∫_{0}^{1/2}x^{p} dx = 2^{−p}/(p + 1),

which together with (4.2) leads to (4.1). The proof of the proposition is completed by

observing d_{p}(W_{n}, W) > 0 for every n, because W_{n} ≠^{L} W.

Note that, in contrast to the asymptotic result (4.1), there is no positive lower bound to d_{p}(W_{n}, W) for a fixed n without further assumptions. Hence the implicit constant in Ω(1/n) in the theorem depends on the variables W, W_{1}, . . ..

For p ≥ 2 we can improve this lower bound by a logarithmic factor by using the known variance of Y_{n}.

Theorem 4.2. If 2 ≤ p < ∞, then

d_{p}(Y_{n}, Y) ≥ d_{2}(Y_{n}, Y) = Ω(\frac{ln n}{n}).

Proof. Recall that Y and Y_{n} have mean 0 and that Var Y = σ^{2} while by (1.2)

Var Y_{n} = σ^{2} − 2\frac{ln n}{n} + O(n^{−1})

and thus

∥Y_{n}∥_{2} = (Var Y_{n})^{1/2} = σ − \frac{1}{σ} \frac{ln n}{n} + O(n^{−1}).

Consequently, for the d_{2}-optimal coupling of Y and Y_{n}, by Minkowski’s inequality,

d_{2}(Y_{n}, Y) = ∥Y_{n} − Y ∥_{2} ≥ ∥Y ∥_{2} − ∥Y_{n}∥_{2} = σ^{−1} \frac{ln n}{n} + O(n^{−1}).

We still have a gap between (ln n)/n and n^{−1/2}.

Remark 4.3. It can be shown that E Y_{n}^{m} = E Y^{m} + O (\frac{ln n}{n}), n ≥ 2, holds also for m = 3, 4, …; cf. the formulas for moments and cumulants by Hennequin [11]. Hence we do not get better lower bounds for d_{p} by considering higher moments.

5 The Kolmogorov–Smirnov distance

Recall that the Kolmogorov–Smirnov distance d_{KS}(F, G) between two distributions is defined as sup_{x∈R} |P (X ≤ x) − P (Y ≤ x)|, when X ∼ F and Y ∼ G. We will in this case also write d_{KS}(X, Y).

To obtain upper bounds for d_{KS}(Y_{n}, Y), we combine the bounds above for d_{p}(Y_{n}, Y) with the following simple general result and the fact [6] that Y has a bounded density function.

Lemma 5.1. Suppose that X and Y are two random variables such that Y is absolutely continuous with a bounded density function f. If M := sup_{y∈R} |f(y)| and 1 ≤ p < ∞, then

d_{KS}(X, Y) ≤ (p + 1)^{1/(p+1)} (Md_{p}(X, Y))^{p/(p+1)} .

Proof. Consider an optimal d_{p}-coupling of X and Y. Then, for x ∈ R and ε > 0, denoting the distribution functions of X and Y by F_{X} and F_{Y},

F_{X}(x) = P (X ≤ x) ≤ P (Y ≤ x) + P (x < Y ≤ x + ε) + P (Y − X > ε)
≤ F_{Y} (x) + Mε + P (Y − X > ε).

Similarly,

F_{Y} (x) ≤ P (X ≤ x) + P (x − ε < Y ≤ x) + P (X − Y > ε)
≤ F_{X}(x) + Mε + P (X − Y > ε).

Consequently,

Δ(x) := |F_{X}(x) − F_{Y} (x)| ≤ Mε + P (|X − Y | > ε)

and thus

d_{p}(X, Y)^{p} = E |X − Y |^{p} = ∫_{0}^{∞}pε^{p−1}P (|X − Y | > ε) dε
≥ ∫_{0}^{Δ(x)/M}pε^{p−1}(Δ(x) − Mε) dε
= \frac{1}{p+1}Δ(x)^{p+1}M^{−p}.

Theorem 5.2. For every ε > 0,

d_{KS}(Y_{n}, Y) = O (n^{ε−(1/2)}) .

Proof. By [6], Y has a bounded density function, so Lemma 5.1 and Theorem 3.1 yield, for every fixed 1 ≤ p < ∞,

d_{KS}(Y_{n}, Y) = O (d_{p}(Y_{n}, Y)^{p/(p+1)}) = O (n^{−p/[2(p+1)]}) .

The result follows by choosing p so large that \frac{p}{2(p+1)} > \frac{1}{2} − ε.

To get an explicit bound we take p = 2 in Lemma 5.1 and use Theorem 2.1. This yields the bound 3^{1/3} (2Mn^{−1/2})^{2/3}, and we know M < 16 from Theorem 3.3 of [6]. Hence,

Theorem 5.3. For n ≥ 1,

d_{KS}(Y_{n}, Y) ≤ (12M^{2})^{1/3}n^{−1/3} < (3072/n)^{1/3} < 15 n^{−1/3}.

Numerical evidence [20] suggests that M < 1, which would give a bound 2.3 n^{−1/3}.

As stated in the introduction, we do not know the right order of decay. The rate O(n^{ε−1/2}) in Theorem 5.2 can be marginally improved to exp(C\sqrt{ln n})n^{−1/2} by checking that the proof of Theorem 3.1 yields c_{p} = O(c^{p}) for some c and then choosing p = (ln n)^{1/2} in Lemma 5.1. We omit the details, since it is likely that this still is far from the truth.

The best lower bound we can prove is Ω(n^{−1}).

Theorem 5.4.

d_{KS}(Y_{n}, Y) > \frac{1}{8(n + 1)}, n ≥ 1.

Again, the lower bound follows from quite general considerations. In this case we use the following lemma.

Lemma 5.5. Suppose that Y and Z are two random variables such that Y has a continuous distribution while a(Z − b) is integer-valued for some real numbers a > 0 and b. If σ_{Z}^{2} := Var Z < ∞, then

d_{KS}(Y, Z) ≥ 1/(12aσ_{Z} + 8).

Proof. For any x ∈ R and δ > 0,

F_{Z}(x + δ) − F_{Y} (x + δ) + F_{Y} (x − δ) − F_{Z}(x − δ) ≤ 2d_{KS}(Y, Z).

Letting δ → 0 we find, since Y is continuous,

P (Z = x) ≤ 2d_{KS}(Y, Z).

The result now follows from the following lemma applied to a(Z − b).

Lemma 5.6. If Z is an integer-valued random variable with finite variance σ_{Z}^{2}, then

sup_{n} P (Z = n) ≥ 1/(6σ_{Z} + 4).

Proof. Let μ := E Z and m := ⌈\frac{3}{2}σ_{Z}⌉. By Chebyshev’s inequality,

P (|Z − μ| ≥ m) ≤ \frac{σ_{Z}^{2}}{m^{2}} ≤ \frac{4}{9} < \frac{1}{2}

and thus

P (μ − m < Z < μ + m) > 1/2.

The interval (μ − m, μ + m) contains at most 2m integers, and thus it must contain an integer n such that

P (Z = n) ≥ \frac{1}{2m}P (μ − m < Z < μ + m) > \frac{1}{4m} > \frac{1}{6σ_{Z} + 4}.

Proof of Theorem 5.4. We apply Lemma 5.5 with a = n and observe that

σ_{Y_{n}} := (Var Y_{n})^{1/2} < σ = (Var Y)^{1/2}  (5.1)

and that 12σ ≐ 7.8 < 8. Indeed, (5.1) is trivial for n = 1 or 2 and easily verified for 3 ≤ n ≤ 6, while for n ≥ 7 it holds because then, by (1.2) and (1.3),

σ^{2} − Var Y_{n} = −4\frac{π^{2}}{6} + 4 (1 + \frac{1}{n})^{2} H_{n}^{(2)} + 2\frac{n + 1}{n^{2}}H_{n} − \frac{13}{n}

> −4 ∑_{k=n+1}^{∞} k^{−2} + \frac{8}{n}H_{n}^{(2)} + \frac{2}{n}H_{n} − \frac{13}{n}

> \frac{1}{n} (−4 + 8H_{n}^{(2)} + 2H_{n} − 13) > 0.

6 Approximating the density of Y

It was shown in [6] that the density f of Y is infinitely differentiable, with all derivatives rapidly decaying. In particular, the derivative f′ is bounded; Theorem 3.3 of [6] gives the explicit bound

M′ := sup_{x∈R} |f′(x)| < 2466.

(This is not very sharp; the true value seems to be less than 2.) The bounds above on the Kolmogorov–Smirnov distance then imply the following local result.

Theorem 6.1. For any x ∈ R and δ > 0,

|\frac{F_{n}(x + \frac{δ}{2}) − F_{n}(x − \frac{δ}{2})}{δ} − f(x)| ≤ \frac{(96M^{2})^{1/3}}{δn^{1/3}} + \frac{M′}{4} δ.

In particular, for any M̄ ≥ M and M̄′ ≥ M′, choosing δ = δ_{n} := 2 (96M̄^{2}(M̄′)^{−3})^{1/6} n^{−1/6} yields

|\frac{F_{n}(x + \frac{δ_{n}}{2} ) − F_{n}(x − \frac{δ_{n}}{2} )}{δ_{n}} − f(x)| ≤ (96M̄^{2}(M̄′)^{3})^{1/6} n^{−1/6}.  (6.1)

The choices M̄ = 16 and M̄′ = 2466 provided by [6] yield the bound 268 n^{−1/6} in (6.1). If M̄ = 1 and M̄′ = 2 could be proven to be legitimate, we could reduce the bound to 3.03 n^{−1/6}.

Proof. By Theorem 5.3,

|F_{n} (x + \frac{δ}{2}) − F_{n} (x − \frac{δ}{2}) − (F (x + \frac{δ}{2}) − F (x − \frac{δ}{2}))| ≤ 2d_{KS}(Y_{n}, Y) ≤ 2(12M^{2})^{1/3}n^{−1/3},

while

|F (x + \frac{δ}{2}) − F (x − \frac{δ}{2}) − δf(x)| = |∫_{−δ/2}^{δ/2}(f(x + y) − f(x)) dy|

≤ ∫_{−δ/2}^{δ/2}M′|y| dy = \frac{M′}{4} δ^{2}.

The first estimate follows, and the second is an immediate consequence.

Theorem 6.1 yields a simple method to numerically calculate the unknown density f up to any given accuracy. For an application, see [2]. (In [2], a preliminary version of Theorem 6.1 with larger constants is used.) Note, however, that the convergence is slow and that it seems impractical to obtain high precision by this method. Other, potentially more powerful, methods to calculate f numerically are discussed in [8].

Open Problem 6.2. Does a local limit theorem hold in the form that

|nP (X_{n} = k) − f (\frac{k − μ_{n}}{n} )| = |nP (Y_{n} = \frac{k − μ_{n}}{n} ) − f (\frac{k − μ_{n}}{n} )| → 0,

perhaps uniformly in k ∈ Z, as n → ∞?

7 Bounds on moment generating functions

Rösler [19] proved that the moment generating functions E e^{λY_{n}} are bounded for fixed λ, and thus E e^{λY_{n}} → E e^{λY} as n → ∞. Rösler did not make his estimates explicit, but his method can be used to obtain explicit bounds. For the limit variable Y, this was done in [8], where we obtained by Rösler’s method (with some refinements) the following explicit estimates for the moment generating function of Y: Let L_{0} ≐ 5.018 be the largest root of e^{L} = 6L^{2}; then

ψ_{Y} (λ) := E e^{λY} ≤ \begin{cases}e^{1.25λ^{2}}, & λ ≤ −0.62, \\ e^{0.5λ^{2}}, & −0.62 ≤ λ ≤ 0, \\ e^{λ^{2}}, & 0 ≤ λ ≤ 0.42, \\ e^{12λ^{2}}, & 0.42 ≤ λ ≤ L_{0}, \\ e^{2e^{λ}}, & L_{0} ≤ λ. \end{cases}  (7.1)

In particular, E e^{λY} ≤ exp (max (12λ^{2}, 2e^{λ})) for all λ ∈ R.

The constants in (7.1) are not sharp, but the doubly exponential growth as λ → +∞ is correct: it was also shown in [8] that ψ_{Y} (λ) ≥ exp (γλ^{−1}e^{λ}) for all large λ whenever γ < 2/e.

In this section we will establish similar bounds for E e^{λY_{n}}. For simplicity we first consider the slight shrinkage

Ŷ_{n} := \frac{n}{n + 1}Y_{n} = \frac{X_{n} − μ_{n}}{n + 1}

of Y_{n}; in particular, Ŷ_{0} := X_{0} − μ_{0} = 0. We then have the following simple result.

Theorem 7.1. E e^{λŶ_{n}} ↑ E e^{λY} as n ↑ ∞. Hence, for any n ≥ 0, E e^{λŶ_{n}} ≤ E e^{λY}, and in particular the upper bounds on E e^{λY} in (7.1) above apply also to E e^{λŶ_{n}}.

Proof. It is well known that the number X_{n} of Quicksort comparisons has the same distribution as the internal path length of a random binary search tree (under the random permutation model) with n internal nodes—see, e.g., [15, Section 6.2.2]. Moreover, it was shown by Régnier [18] that when X_{n} is reinterpreted as the internal path length of an evolving random binary search tree after n keys have been inserted, the process (Ŷ_{n})_{n≥0} is a martingale, which is L^{2}-bounded and thus converges a.s. and in L^{2} to some limit Y. It follows that also Y_{n} → Y a.s., and thus in distribution; hence this random variable Y is (a realization of) the same Y as above.

The martingale property can be written Ŷ_{n} = E (Ŷ_{n+1}|F_{n}), for the appropriate σ-field F_{n}. Since x ↦ e^{λx} is convex, it now follows by Jensen’s inequality for conditional expectations that e^{λŶ_{n}} ≤ E (e^{λŶ_{n+1}}|F_{n}); and thus, taking expectations, E e^{λŶ_{n}} ≤ E e^{λŶ_{n+1}}.

By the same argument, E e^{λŶ_{n}} ≤ E e^{λY} for each n ≥ 0, which together with Fatou’s lemma yields E e^{λŶ_{n}} → E e^{λY} as n → ∞.

Corollary 7.2. For every n ≥ 1, we have

E e^{λY_{n}} ≤ \begin{cases}e^{1.25[1+(1/n)]^{2}λ^{2}}, & λ ≤ 0, \\ e^{0.5[1+(1/n)]^{2}λ^{2}}, & −0.62 n/(n + 1) ≤ λ ≤ 0, \\ e^{[1+(1/n)]^{2}λ^{2}}, & 0 ≤ λ ≤ 0.42 n/(n + 1), \\ e^{12[1+(1/n)]^{2}λ^{2}}, & 0 ≤ λ ≤ L_{0} n/(n + 1), \\ e^{2e^{[1+(1/n)]λ}}, & L_{0} n/(n + 1) ≤ λ. \end{cases}

In particular, E e^{λY_{n}} ≤ exp (max (12[1 + (1/n)]^{2}λ^{2}, 2e^{[1+(1/n)]λ})) for all λ ∈ R.

Proof. λY_{n} = λ_{n}Ŷ_{n} with λ_{n} := [1 + (1/n)]λ.

Remark 7.3. The factors [1+(1/n)] in Corollary 7.2 are annoying but hardly important in applications. With some effort, we have been able to modify the proof in [8] and obtain for λ ≥ −0.58 the same estimates for E e^{λY_{n}} as obtained there for E e^{λY}; for λ < −0.58 we only obtain a slightly weaker bound, which for large n is inferior to the bound in Corollary 7.2. More precisely, we have shown

E e^{λY_{n}} ≤ \begin{cases}e^{1.34λ^{2}}, & λ ≤ −0.58, \\ e^{0.5λ^{2}}, & −0.58 ≤ λ ≤ 0, \\ e^{λ^{2}}, & 0 ≤ λ ≤ 0.42, \\ e^{12λ^{2}}, & 0.42 ≤ λ ≤ L_{0}, \\ e^{2e^{λ}}, & L_{0} ≤ λ. \end{cases}  (7.2)

In particular, E e^{λY_{n}} ≤ exp (max (12λ^{2}, 2e^{λ})) for all λ ∈ R. In other words, we can eliminate the factors [1 + (1/n)] in Corollary 7.2 for λ ≥ −0.58 (and in particular for all positive λ). Since the proof is quite long and the result only marginally improves Corollary 7.2, we give the proof not here but rather in a separate appendix [9].

It seems likely that with further effort one could remove the factor [1 + (1/n)] for λ < −0.58 too, so that all the bounds in (7.1) also would bound e^{λY_{n}}. Moreover, it seems quite likely that E e^{λY_{n}} ≤ E e^{λY} holds for all λ and n, and perhaps even that E e^{λY_{n}} ↑ E e^{λY}, as was proved for Ŷ_{n} in Theorem 7.1.

Theorem 7.1 enables us to get an explicit constant in Rösler’s [19] large deviation bound.

Corollary 7.4. For any ε > 0 and λ > 0,

P (|X_{n} − μ_{n}| ≥ εμ_{n}) ≤ 2 exp [3ελ + max (12λ^{2}, 2e^{λ})] n^{−2ελ}.

Proof. By Markov’s inequality,

P (|X_{n} − μ_{n}| ≥ εμ_{n}) = P (|Ŷ_{n}| ≥ εμ_{n}/(n + 1))
≤ exp(−ελμ_{n}/(n + 1))E e^{λ|Ŷ_{n}|}
≤ exp(−ελμ_{n}/(n + 1)) (E e^{λŶ_{n}} + E e^{−λŶ_{n}}) .

The result follows from Theorem 7.1, since μ_{n}/(n + 1) ≥ 2H_{n} − 4 ≥ 2 ln n − 3 by (1.7) and (1.9).

Corollary 7.5. For any fixed ε > 0,

P (|X_{n} − μ_{n}| ≥ εμ_{n}) ≤ n^{−2ε ln ln n+O(1)}, n ≥ 2.

Proof. Take (for n ≥ 3) λ = ln ln n in Corollary 7.4.

The bound in Corollary 7.5 is essentially the same as the one obtained by McDiarmid and Hayward [16] by different methods; the only difference is a slight improvement in the error term. More generally, McDiarmid and Hayward [16] considered ε varying with n such that 1/ ln n < ε ≤ 1; if we take λ = ln ln n + ln ε in Corollary 7.4, we obtain the bound in their Theorem 1.1 with O(ln ln ln n) improved to O(1). Compare also the related large deviation estimates for the limit distribution Y in [8] (by similar arguments) and Knessl and Szpankowski [13] (much more precise, but assuming an as yet unverified regularity hypothesis).

Finally we show that the rate of convergence of the moment generating functions E e^{λY_{n}} to E e^{λY} also is O(n^{−1/2}). (The same holds for E e^{λŶ_{n}}.)

Theorem 7.6. For any fixed complex λ,

E e^{λY_{n}} = E e^{λY} + O(n^{−1/2}).

Explicitly, with λ_{1} := Re(λ),

|E e^{λY_{n}} − E e^{λY} | ≤ 3|λ| exp [max (24[1 + (1/n)]^{2}λ_{1}^{2}, e^{2[1+(1/n)]λ_{1}})] n^{−1/2}.

Proof. Consider a d_{2}-optimal coupling of Y_{n} and Y. Then, using the mean value theorem, the Cauchy–Schwarz inequality, Corollary 7.2, and (7.1),

|E e^{λY_{n}} − E e^{λY} | ≤ E |e^{λY_{n}} − e^{λY} |
≤ E (|λ||Y_{n} − Y | e^{max(λ_{1}Y_{n},λ_{1}Y)})
≤ |λ| (E |Y_{n} − Y |^{2})^{1/2} (E e^{2 max(λ_{1}Y_{n},λ_{1}Y)})^{1/2}
≤ |λ|d_{2}(Y_{n}, Y) (E e^{2λ_{1}Y_{n}} + E e^{2λ_{1}Y} )^{1/2}
≤ \sqrt{2} |λ| exp [max (24[1 + (1/n)]^{2}λ_{1}^{2}, e^{2[1+(1/n)]λ_{1}})] d_{2}(Y_{n}, Y).

The result follws by Theorem 2.1.

Remark 7.7. By Remark 7.3, the factors [1+(1/n)] can be eliminated in the statement of Theorem 7.6.

Acknowledgment. We thank Anhua Lin and Ludger Rüschendorf for helpful discussions. We also thank two anonymous referees for useful suggestions.

References

[1] Cambanis, S., Simons, G., and Stout, W. Inequalities for E k(X, Y) when the marginals are fixed. Z. Wahrscheinlichkeitstheorie verw. Gebiete 36 (1976), 285– 294.

[2] Devroye, L., Fill, J. A., and Neininger, R. Perfect simulation from the Quicksort limit distribution. Electron. Comm. Probab. 5 (2000), 95–99 (electronic).

[3] Dongarra, J. and Sullivan, F. Guest editors’ introduction: the top 10 algorithms. Computing in Science & Engineering 2:1 (2000), 22–23.

[4] Eddy, W. F. and Schervish, M. J. How many comparisons does Quicksort use? J. Algorithms 19 (1995), 402–431.

[5] Fill, J. A. On the distribution for binary search trees under the random permutation model. Random Structures and Algorithms 8 (1996), 1–25.

[6] Fill, J. A. and Janson, S. Smoothness and decay properties of the limiting Quicksort density function. In D. Gardy and A. Mokkadem, editors, Mathematics and Computer Science: Algorithms, Trees, Combinatorics and Probabilities, Trends in Mathematics, pages 53–64. Birkhäuser Verlag, 2000.

[7] Fill, J. A. and Janson, S. A characterization of the set of fixed points of the Quicksort transformation. Electron. Comm. Probab. 5 (2000), 77–84 (electronic).

[8] Fill, J. A. and Janson, S. Approximating the limiting Quicksort distribution. Random Struct. Alg. 19 (2001), 376–406.

[9] Fill, J. A. and Janson, S. Appendix to Quicksort asymptotics. Unpublished appendix to this paper, available from http://www.mts.jhu.edu/~fill/ or

http://www.math.uu.se/~svante/.

[10] Hennequin, P. Combinatorial analysis of quicksort algorithm. RAIRO Inform. Théor. Appl. 23 (1989), 317–333.

[11] Hennequin, P. Analyse en moyenne d’algorithmes, tri rapide et arbres de recherche. Ph.D. dissertation, L’École Polytechnique Palaiseau.

[12] Hoare, C. A. R. Quicksort. Comput. J. 5 (1962), 10–15.

[13] Knessl, C. and Szpankowski, W. Quicksort algorithm again revisited. Discrete Math. Theor. Comput. Sci. 3 (1999), 43–64.

[14] Knuth, D. E. The Art of Computer Programming. Volume 1. Fundamental Algorithms. Third edition. Addison–Wesley, Reading, Mass., 1997.

[15] Knuth, D. E. The Art of Computer Programming. Volume 3. Sorting and Searching. Second edition. Addison–Wesley, Reading, Mass., 1998.

[16] McDiarmid, C. J. H. and Hayward, R. B. Large deviations for Quicksort. J. Algorithms 21 (1996), 476–507.

[17] Neininger, R. and Rüschendorf, L. Rates of convergence for Quicksort. J. Algorithms, to appear.

[18] Régnier, M. A limiting distribution for quicksort. RAIRO Inform. Théor. Appl. 23 (1989), 335–343.

[19] Rösler, U. A limit theorem for ‘Quicksort’. RAIRO Inform. Théor. Appl. 25 (1991), 85–100.

[20] Tan, K. H. and Hadjicostas, P. Some properties of a limiting distribution in Quicksort. Statist. Probab. Lett. 25 (1995), 87–94.