Smoothness and Decay Properties of the Limiting Quicksort Density Function

James Allen Fill^{1} 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

Using Fourier analysis, we prove that the limiting distribution of the standardized random number of comparisons used by Quicksort to sort an array of n numbers has an everywhere positive and infinitely differentiable density f, and that each derivative f^{(k)} enjoys superpolynomial decay at ±∞. In particular, each f^{(k)} is bounded. Our method is sufficiently computational to prove, for example, that f is bounded by 16. AMS 2000 subject classifications. Primary 68W40; secondary 68P10, 60E05, 60E10. Key words and phrases. Quicksort, density, characteristic function, sorting algorithm, Fourier analysis, rapidly decreasing C^{∞} function, tempered distribution, integral equation. Date. March 31, 2000.

1Research supported by NSF grant DMS–9803780, and by the Acheson J. Duncan Fund for the Advancement of Research in Statistics.

1 Introduction and summary

The Quicksort algorithm of Hoare [7] is “one of the fastest, the best-known, the most generalized, the most completely analyzed, and the most widely used algorithms for sorting an array of numbers” [2]. Quicksort is the standard sorting procedure in Unix systems, and Philippe Flajolet, a leader in the field of analysis of algorithms, has noted that it is among “some of the most basic algorithms—the ones that do deserve deep investigation” [4]. Our goal in this introductory section is to review briefly some of what is known about the analysis of Quicksort and to summarize how this paper advances that analysis.

The Quicksort algorithm for sorting an array of n 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,

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} := EX_{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 standard deviation of X_{n} (see Exercise 6.2.2-8 in [9]), which turns out to be ∼ n\sqrt{7 − \frac{2}{3}π^{2}}.

Consider the standardized variate

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

Régnier [11] showed using martingale arguments that Y_{n} → Y in distribution, with Y satisfying the distributional identity

Y =^{L} UY + (1 − U)Z + g(U) =: h_{Y,Z}(U), (1.1)

where

g(u) := 2u ln u + 2(1 − u) ln(1 − u) + 1, (1.2) and where, on the right of =^{L} in (1.1), U, Y, and Z are independent, with Z =^{L} Y and U ∼ unif(0, 1). Rösler [12] showed that (1.1) characterizes the limiting law L(Y), in the precise sense that F := L(Y) is the unique fixed point of the operator

G = L(V) ↦ SG := L(UV + (1 − U)V^{∗} + g(U))

(in what should now be obvious notation) subject to

EV = 0, VarV < ∞.

Thus it is clear that fundamental (asymptotic) probabilistic understanding of Quicksort’s behavior relies on fundamental understanding of the limiting distribution F. In this regard, Rösler [12] showed that

(1.3) the moment generating function (mgf) of Y is everywhere finite,

and Hennequin [5] [6] and Rösler showed how all the moments of Y can be pumped out one at a time, though there is no known expression for the mgf nor for the general pth moment in terms of p. Tan and Hadjicostas [15] proved that F has a density f which is almost everywhere positive, but their proof does not even show whether f is continuous.

The main goal of this paper is to prove that F has a density f which is infinitely differentiable, and that each derivative f^{(k)}(y) decays as y → ±∞ more rapidly than any power of |y|^{−1}: this is our main Theorem 3.1. In particular, it follows that each f^{(k)} is bounded (cf. Theorem 3.3).

Our main tool will be Fourier analysis. We begin in Section 2 by showing (see Theorem 2.9) that the characteristic function ϕ for F has rapidly decaying derivatives of every order. Standard arguments reviewed briefly at the outset of Section 3 then immediately carry this result over from ϕ to f. Finally, in Section 4 we will use the boundedness and continuity of f to establish an integral equation for f (Theorem 4.1). As a corollary, f is everywhere positive (Corollary 4.2).

Remark 1.1. (a) Our method is sufficiently computational that we will prove, for example, that f is bounded by 16. This is not sharp numerically, as Figure 4 of [15] strongly suggests that the maximum value of f is about 2/3. However, in future work we will rigorously justify (and discuss how to obtain bounds on the error in) the numerical computations used to obtain that figure, and the rather crude bounds on f and its derivatives obtained in the present paper are needed as a starting point for that more refined work.

(b) Very little is known rigorously about f. For example, the figure discussed in (a) indicates that f is unimodal. Can this be proved? Is f in fact strongly unimodal (i.e., log-concave)? What can one say about changes of signs for the derivatives of f?

(c) Knessl and Szpankowski [8] purport to prove very sharp estimates of the rates of decay of f(y) as y → −∞ and as y → ∞. Roughly put, they assert that the left tail of f decays doubly exponentially (like the tail of an extreme-value density) and that the right tail decays exponentially. But their results rely on several unproven assumptions (as noted in their paper). Among these, for example, is their assumption (59) that

Ee^{−λY} ∼ exp(αλ ln λ + βλ + γ ln λ + δ) as λ → ∞

for some constants α(> 0), β, γ, δ. (Having assumed this, they derive the values of α, γ, and δ exactly, and the value of β numerically.)

2 Bounds on the limiting Quicksort characteristic function

We will in this section prove the following result on superpolynomial decay of the characteristic function of the limit variable Y.

Theorem 2.1. For every real p ≥ 0 there is a smallest constant 0 < c_{p} < ∞ such that the characteristic function ϕ(t) :≡ Ee^{itY} satisfies

|ϕ(t)| ≤ c_{p}|t|^{−p} for all t ∈ R. (2.1)

These best possible constants c_{p} satisfy c_{0} = 1, c_{1/2} ≤ 2, c_{3/4} ≤ \sqrt{8π}, c_{1} ≤ 4π, c_{3/2} < 187, c_{5/2} < 103215, c_{7/2} < 197102280, and the relations

(2.2)  c_{p_{1}}^{1/p_{1}} ≤ c_{p_{2}}^{1/p_{2}},  0 < p_{1} ≤ p_{2};
(2.3)  c_{p+1} ≤ 2^{p+1}c_{p}^{1+(1/p)}p/(p − 1),  p > 1;
(2.4)  c_{p} ≤ 2^{p^{2}+6p},  p > 0.

[The numerical bounds are not sharp (except in the trivial case of c_{0}); they are the best that we can get without too much work, but we expect that substantial improvements are possible.]

Proof. The basic approach is to use the fundamental relation (1.1). We will first show, using a method of van der Corput [1, 10], that the characteristic function of h_{y,z}(U) is bounded by 2|t|^{−1/2} for each y, z. Mixing, this yields Theorem 2.1 for p = 1/2. Then we will use another consequence of (1.1), namely, the functional equation

ϕ(t) = ∫_{u=0}^{1}ϕ(ut) ϕ((1 − u)t) e^{itg(u)} du, t ∈ R, (2.5)

or rather its consequence

|ϕ(t)| ≤ ∫_{u=0}^{1}|ϕ(ut)| |ϕ((1 − u)t)| du, (2.6)

and obtain successive improvements in the exponent p.

We give the details as a series of lemmas, beginning with a standard calculus estimate [10]. Note that it suffices to consider t > 0 in the proofs because ϕ(−t) = \overline{ϕ(t)} and thus |ϕ(−t)| = |ϕ(t)|. Note also that the best constants satisfy c_{p} = sup_{t>0} t^{p}|ϕ(t)| (although we do not know in advance of proving Theorem 2.1 that these are finite), and thus c_{p}^{1/p} = sup_{t>0} t|ϕ(t)|^{1/p}, which clearly satisfies (2.2) because |ϕ(t)| ≤ 1.

Lemma 2.2. Suppose that a function h is twice continuously differentiable on an open interval (a, b) with

h′(x) ≥ c > 0 and h′′(x) ≥ 0 for x ∈ (a, b).

Then

|∫_{x=a}^{b}e^{ith(x)} dx| ≤ \frac{2}{ct} for all t > 0.

Proof. By considering subintervals (a + ε, b − ε) and letting ε → 0, we may without loss of generality assume that h is defined and twice differentiable at the endpoints, too. Then, using integration by parts, we calculate

∫_{x=a}^{b}e^{ith(x)} dx = \frac{1}{it} ∫_{x=a}^{b}[\frac{d}{dx}e^{ith(x)}] \frac{dx}{h′(x)}

= \frac{1}{it} \{ \frac{e^{ith(x)}}{h′(x)} |_{x=a}^{b} − ∫_{x=a}^{b}e^{ith(x)} d (\frac{1}{h′(x)})\} .

So

|∫_{x=a}^{b}e^{ith(x)} dx| ≤ \frac{1}{t} \{(\frac{1}{h′(b)} + \frac{1}{h′(a)}) + ∫_{x=a}^{b}|d (\frac{1}{h′(x)})| dx\}

= \frac{1}{t} \{(\frac{1}{h′(b)} + \frac{1}{h′(a)}) + ∫_{x=a}^{b}[−d (\frac{1}{h′(x)})] dx\}

= \frac{1}{t} \{(\frac{1}{h′(b)} + \frac{1}{h′(a)}) + (\frac{1}{h′(a)} − \frac{1}{h′(b)})\}

= \frac{2}{th′(a)} ≤ \frac{2}{ct}.

Lemma 2.3. For any real numbers y and z, the random variable h_{y,z}(U) defined by (1.1) satisfies

|Ee^{ith_{y,z}(U)}| ≤ 2|t|^{−1/2}.

Proof. We will apply Lemma 2.2, taking h to be h_{y,z}. Observe that

h′′_{y,z}(u) = 2 (\frac{1}{u} + \frac{1}{1 − u}) = \frac{2}{u(1 − u)} ≥ 8 for u ∈ (0, 1)

and that

h′_{y,z}(u) = 0 if and only if u = α_{y,z} := \frac{1}{1 + exp (\frac{1}{2}(y − z))} ∈ (0, 1).

Let t > 0 and γ > 0. If in Lemma 2.2 we take a := α_{y,z} + γt^{−1/2} and b := 1, and assume that a < b, then note

h′(u) = h′_{y,z}(u) = ∫_{x=α_{y,z}}^{u} h′′_{y,z}(x) dx ≥ 8(u − α_{y,z}) ≥ 8γt^{−1/2} for all u ∈ (a, b).

So, by Lemma 2.2,

|∫_{u=α_{y,z}+γt^{−1/2}}^{1}e^{ith_{y,z}(u)} du| ≤ \frac{2}{t}[8γt^{−1/2}]^{−1} = \frac{1}{4γ} t^{−1/2}.

Trivially,

|∫_{u=α_{y,z}}^{α_{y,z}+γt^{−1/2}}e^{ith_{y,z}(u)} du| ≤ γt^{−1/2},

so we can conclude

|∫_{u=α_{y,z}}^{1}e^{ith_{y,z}(u)} du| ≤ [(4γ)^{−1} + γ]t^{−1/2}.

This result is trivially also true when a = α_{y,z} + γt^{−1/2} ≥ b = 1, so it holds for all t, γ > 0. The optimal choice of γ here is γ_{opt} = 1/2, which yields

|∫_{u=α_{y,z}}^{1}e^{ith_{y,z}(u)} du| ≤ t^{−1/2} for all t > 0.

Similarly, for example by considering u ↦ h(1 − u),

|∫_{0}^{α_{y,z}}e^{ith_{y,z}(u)} du| ≤ t^{−1/2} for all t > 0,

and we conclude that the lemma holds for all t > 0, and thus for all real t.

Lemma 2.4. For any real t, |ϕ(t)| ≤ 2|t|^{−1/2}.

Proof. Lemma 2.3 shows that

|E(e^{ith_{Y,Z}(U)} | Y, Z)| ≤ 2|t|^{−1/2}

and thus

|ϕ(t)| = |Ee^{ith_{Y,Z}(U)}| ≤ E|E(e^{ith_{Y,Z}(U)} | Y, Z)| ≤ 2|t|^{−1/2}.

The preceding lemma is the case p = 1/2 of Theorem 2.1. We now improve the exponent.

Lemma 2.5. Let 0 < p < 1. Then

c_{2p} ≤ \frac{[Γ(1 − p)]^{2}}{Γ(2 − 2p)} c_{p}^{2}.

Proof. By (2.6) and the definition of c_{p},

|ϕ(t)| ≤ ∫_{u=0}^{1} c_{p}^{2}|ut|^{−p}|(1 − u)t|^{−p} du = c_{p}^{2}|t|^{−2p} ∫_{u=0}^{1} u^{−p}(1 − u)^{−p} du,

and the result follows by evaluating the beta integral.

In particular, recalling Γ(1/2) = \sqrt{π}, Lemmas 2.4 and 2.5 yield

|ϕ(t)| ≤ \frac{4π}{|t|} . (2.7)

This proves (2.1) for p = 1, with c_{1} ≤ 4π, and thus by (2.2) for every p ≤ 1 with c_{p} ≤ (4π)^{p}; applying Lemma 2.5 again, we obtain the finiteness of c_{p} in (2.1) for all p < 2. Somewhat better numerical bounds are obtained for 1/2 < p < 1 by taking a geometric average between the cases p = 1/2 and p = 1: the inequality

|ϕ(t)| ≤ (2t^{−1/2})^{2−2p}(4πt^{−1})^{2p−1} = 2^{2p}π^{2p−1}t^{−p}, t > 0,

shows that c_{p} ≤ 2^{2p}π^{2p−1}, 1/2 ≤ p ≤ 1. In particular, we have c_{3/4} ≤ \sqrt{8π}, and thus, by Lemma 2.5, c_{3/2} ≤ 8π^{1/2}[Γ(1/4)]^{2} < 186.4 < 187.

Lemma 2.6. Let p > 1. Then

c_{p+1} ≤ 2^{p+1}c_{p}^{1+(1/p)}p/(p − 1).

Proof. Assume that t ≥ 2c_{p}^{1/p}. Then, again using (2.6),

|ϕ(t)| ≤ ∫_{u=0}^{1}min (\frac{c_{p}}{(ut)^{p}} , 1) min (\frac{c_{p}}{[(1 − u)t]^{p}} , 1) du

= 2 ∫_{u=0}^{c_{p}^{1/p}t^{−1}} \frac{c_{p}}{[(1 − u)t]^{p}} du + ∫_{u=c_{p}^{1/p}t^{−1}}^{1−c_{p}^{1/p}t^{−1}} \frac{c_{p}^{2}}{[u(1 − u)t^{2}]^{p}} du

≤ \frac{2}{[1 − c_{p}^{1/p}t^{−1}]^{p}} \frac{c_{p}^{1+(1/p)}}{t^{p+1}} + 2\frac{c_{p}^{2}}{t^{2p}} ∫_{u=c_{p}^{1/p}t^{−1}}^{1/2} \frac{du}{[u(1 − u)]^{p}}

≤ \frac{2}{(1/2)^{p}} c_{p}^{1+(1/p)}t^{−(p+1)} + \frac{2}{(1/2)^{p}} \frac{c_{p}^{2}}{t^{2p}} ∫_{u=c_{p}^{1/p}t^{−1}}^{1/2}u^{−p} du

≤ 2^{p+1} \{c_{p}^{1+(1/p)}t^{−(p+1)} + \frac{1}{p − 1}c_{p}^{2}t^{−2p} [c_{p}^{1/p}t^{−1}]^{−(p−1)}\}

= 2^{p+1}c_{p}^{1+(1/p)} \frac{p}{p − 1}t^{−(p+1)}.

We have derived the desired bound for all t ≥ 2c_{p}^{1/p}. But also, for all 0 < t < 2c_{p}^{1/p}, we have

2^{p+1}c_{p}^{1+(1/p)} \frac{p}{p − 1}t^{−(p+1)} ≥ \frac{p}{p − 1} ≥ 1 ≥ |ϕ(t)|,

so the estimate holds for all t > 0.

Lemma 2.6 completes the proof of finiteness of every c_{p} in (2.1) (by induction), and of the estimate (2.3). The bound for c_{3/2} obtained above now shows (using Maple) that c_{5/2} < 103215, which then gives c_{7/2} < 197102280.

We can rewrite (2.3) as

c_{p+1}^{1/(p+1)} ≤ 2c_{p}^{1/p}(1 + \frac{1}{p − 1})^{1/(p+1)} ≤ 2c_{p}^{1/p} exp(\frac{1}{(p − 1)(p + 1)})

= 2c_{p}^{1/p} exp(\frac{1}{2(p − 1)} − \frac{1}{2(p + 1)}).

Hence, by induction, if p = n + \frac{5}{2} for a nonnegative integer n, then

c_{p}^{1/p} ≤ 2^{n}c_{5/2}^{2/5}e^{(1/3)+(1/5)} = C2^{p},

where C := 2^{−5/2}e^{8/15}c_{5/2}^{2/5} < 30.6 < 2^{5}, using the above estimate of c_{5/2}. Consequently,

c_{p}^{1/p} < 2^{p+5} when p = n + \frac{5}{2}. For general p > 3/2 we now use (2.2) with p_{1} = p and p_{2} = ⌈p − \frac{5}{2}⌉ + \frac{5}{2}, obtaining c_{p}^{1/p} < 2^{p_{2}+5} < 2^{p+6}; the case p ≤ 3/2 follows from (2.2) and the estimate c_{3/2}^{2/3} < 33 < 2^{6}. This completes the proof of (2.4) and hence of Theorem 2.1.

Remark 2.7. We used (1.1) in two different ways. In the first step we conditioned on the values of Y and Z, while in the inductive steps we conditioned on U.

Remark 2.8. A variety of other bounds are possible. For example, if we begin with the inequality (2.7), use (2.6), and proceed just as in the proof of Lemma 2.6, we can easily derive the following result in the case t ≥ 8π:

(2.8) |ϕ(t)| ≤ \frac{32π^{2}}{t^{2}} (ln(\frac{t}{4π} ) + 2) ≤ \frac{32π^{2} ln t}{t^{2}} for all t ≥ 1.72.

The result is trivial for 1.72 ≤ t < 8π, since then the bounds exceed unity.

Since Y has finite moments of all orders [recall (1.3)], the characteristic function ϕ is infinitely differentiable. Theorem 2.1 implies a rapid decrease of all derivatives, too.

Theorem 2.9. For each real p ≥ 0 and integer k ≥ 0, there is a constant c_{p,k} such that

|ϕ^{(k)}(t)| ≤ c_{p,k}|t|^{−p} for all t ∈ R.

Proof. The case k = 0 is Theorem 2.1, and the case p = 0 follows by |ϕ^{(k)}(t)| ≤ E|Y |^{k}. The remaining cases follows from these cases by induction on k and the following calculus lemma.

Lemma 2.10. Suppose that g is a complex-valued function on (0, ∞) and that A, B, p > 0 are such that |g(t)| ≤ At^{−p} and |g′′(t)| ≤ B for all t > 0. Then |g′(t)| ≤ 2\sqrt{AB}t^{−p/2}.

Proof. Fix t > 0 and let θ = arg(g′(t)). For s > t,

Re(e^{−iθ}g′(s)) ≥ Re(e^{−iθ}g′(t)) − |g′(s) − g′(t)| ≥ |g′(t)| − B(s − t)

and thus, integrating from t to t_{1} := t + (|g′(t)|/B),

\begin{aligned}Re(e^{−iθ}(g(t_{1}) − g(t))) & ≥ ∫_{t}^{t_{1}}(|g′(t)| − B(s − t)) ds \\ & = (t_{1} − t)|g′(t)| − \frac{1}{2}B(t_{1} − t)^{2} = |g′(t)|^{2}/(2B). \end{aligned}

Consequently,

|g′(t)|^{2}/(2B) ≤ |g(t)| + |g(t_{1})| ≤ 2At^{−p},

and the result follows.

In other words, the characteristic function ϕ belongs to the class S of infinitely differentiable functions that, together with all derivatives, decrease more rapidly than any power. (This is the important class of test functions for tempered distributions, introduced by Schwartz [14]; it is often called the class of rapidly decreasing C^{∞} functions.)

3 The limiting Quicksort density f and its derivatives

We can now improve the result by Tan and Hadjicostas [15] on existence of a density f for Y. It is an immediate consequence of Theorem 2.1, with p = 0 and p = 2, say, that the characteristic function ϕ is integrable over the real line. It is well-known—see, e.g., [3, Theorem XV.3.3]—that this implies that Y has a bounded continuous density f given by the Fourier inversion formula

f(x) = \frac{1}{2π} ∫_{t=−∞}^{∞}e^{−itx} ϕ(t) dt, x ∈ R. (3.1)

Moreover, using Theorem 2.1 with p = k + 2, we see that t^{k}ϕ(t) is also integrable for each integer k ≥ 0, which by a standard argument (cf. [3, Section XV.4]) shows that f is infinitely smooth, with a kth derivative (k ≥ 0) given by

f^{(k)}(x) = \frac{1}{2π} ∫_{t=−∞}^{∞}(−it)^{k} e^{−itx} ϕ(t) dt, x ∈ R. (3.2)

It follows further that the derivatives are bounded, with

sup_{x} |f^{(k)}(x)| ≤ \frac{1}{2π} ∫_{t=−∞}^{∞} |t|^{k} |ϕ(t)| dt (k ≥ 0), (3.3)

and these bounds in turn can be estimated using Theorem 2.1. Moreover, as is well known [14], [13, Theorem 7.4], an extension of this argument shows that the class S discussed at the end of Section 2 is preserved by the Fourier transform, and thus Theorem 2.9 implies that f ∈ S:

Theorem 3.1. The Quicksort limiting distribution has an infinitely differentiable density function f. For each real p ≥ 0 and integer k ≥ 0, there is a constant C_{p,k} such that

|f^{(k)}(x)| ≤ C_{p,k}|x|^{−p} for all x ∈ R.

For numerical bounds on f, we can use (3.3) with k = 0 and Theorem 2.1 for several different p (in different intervals); for example, using p = 0, 1/2, 1, 3/2, 5/2, 7/2, and taking t_{1} = 4, t_{2} = 4π^{2}, t_{3} = (187/(4π))^{2}, t_{4} = 103215/187, t_{5} = 197102280/103215,

\begin{aligned}f(x) & ≤ \frac{1}{2π} ∫_{t=−∞}^{∞} |ϕ(t)| dt = \frac{1}{π} ∫_{t=0}^{∞} |ϕ(t)| dt \\ & ≤ \frac{1}{π} ∫_{t=0}^{∞} min(1, 2t^{−1/2}, 4πt^{−1}, 187t^{−3/2}, 103215 t^{−5/2}, 197102280 t^{−7/2}) dt \\ & = \frac{1}{π} (∫_{t=0}^{t_{1}} dt + ∫_{t=t_{1}}^{t_{2}} 2t^{−1/2} dt + ∫_{t=t_{2}}^{t_{3}} 4πt^{−1} dt + ∫_{t=t_{3}}^{t_{4}} 187 t^{−3/2} dt \\ & + ∫_{t=t_{4}}^{t_{5}} 103215 t^{−5/2} dt + ∫_{t=t_{5}}^{∞} 197102280 t^{−7/2} dt) \\ & ≤ 18.2. \end{aligned} (3.4)

Remark 3.2. We can do somewhat better by using the first bound in (2.8) over the interval (103.18, 1984) instead of (as above) Theorem 2.1 with p = 1, 3/2, 5/2, 7/2 over (103.18, t_{3}), (t_{3}, t_{4}), (t_{4}, t_{5}), (t_{5}, 1984), respectively. This gives

f(x) < 15.3.

Similarly, (3.3) with k = 1 and the same estimates of |ϕ(t)| as in (3.4) yield

|f′(x)| ≤ \frac{1}{2π} ∫_{t=−∞}^{∞} |t||ϕ(t)| dt = \frac{1}{π} ∫_{t=0}^{∞} t|ϕ(t)| dt < 3652.1,

which can be reduced to 2492.1 by proceeding as in Remark 3.2. The bound can be further improved to 2465.9 by using also p = 9/2.

Somewhat better bounds are obtained by using more values of p in the estimates of the integrals, but the improvements obtained in this way seem to be slight. We summarize the bounds we have obtained.

Theorem 3.3. The limiting Quicksort density function f satisfies max_{x} f(x) < 16 and max_{x} |f′(x)| < 2466.

The numerical bounds obtained here are far from sharp; examination of Figure 4 of [15] suggests that max f < 1 and max |f′| < 2. Our present technique cannot hope to produce a better bound on f than 4/π > 1.27, since neither Lemma 2.3 nor (2.6) can improve on the bound |ϕ(t)| ≤ 1 for |t| ≤ 4. Further, no technique based on (3.3) can hope to do better than the actual value of (2π)^{−1} ∫_{t=−∞}^{∞}|ϕ(t)| dt, which from cursory examination of Figure 6 of [15] appears to be about 2.

4 An integral equation for the density f

Our estimates are readily used to justify rigorously the following functional equation.

Theorem 4.1. The continuous limiting Quicksort density f satisfies (pointwise) the integral equation

f(x) = ∫_{u=0}^{1} ∫_{y∈R} f(y) f (\frac{x − g(u) − (1 − u)y}{u} ) \frac{1}{u} dy du, x ∈ R,

where g(·) is as in (1.2).

Proof. For each u with 0 < u < 1, the random variable

uY + (1 − u)Z + g(u) (4.1)

[with notation as in (1.1)] has the density function

f_{u}(x) := ∫_{z∈R} f(z) f (\frac{x − g(u) − (1 − u)z}{u} ) \frac{1}{u} dz, (4.2)

where the integral converges for each x since, using Theorem 3.3, the integrand is bounded by f(z)(max f)/u ≤ 16f(z)/u; dominated convergence using the continuity of f and the same bound shows further that f_{u} is continuous.

This argument yields the bound f_{u}(x) ≤ 16/u, and since f_{u} = f_{1−u} by symmetry in (4.1), we have f_{u}(x) ≤ 16/ max(u, 1 − u) ≤ 32. This uniform bound, (1.1), and dominated convergence again imply that ∫_{0}^{1} f_{u}(x) du is a continuous density for Y, and thus equals f(x) for every x.

It was shown in [15] that f is positive almost everywhere; we now can improve this by removing the qualifier “almost.”

Corollary 4.2. The continuous limiting Quicksort density function is everywhere positive.

Proof. We again use the notation (4.2) from the proof of Theorem 4.1. Fix x ∈ R and u ∈ (0, 1). Since f is almost everywhere positive [15], the integrand in (4.2) is positive almost everywhere. Therefore f_{u}(x) > 0. Now we integrate over u ∈ (0, 1) to conclude that f(x) > 0.

Alternatively, Corollary 4.2 can be derived directly from Theorem 4.1, without recourse to [15]. Indeed, if f(y_{0}) > 0 and u_{0} ∈ (0, 1), set x = y_{0} + g(y_{0}); then the integrand in the double integral for f(x) in Theorem 4.1 is postive for (u, y) equal to (u_{0}, y_{0}), and therefore, by continuity, also in some small neighborhood thereof. It follows that f(y_{0} + g(u_{0})) > 0. Since u_{0} is arbitrary and the image of (0, 1) under g is (−(2 ln 2 − 1), 1), an open interval containing the origin, Corollary 4.2 follows readily.

Remark 4.3. In future work, we will use arguments similar to those of this paper, together with other arguments, to show that when one applies the method of successive substitutions to the integral equation in Theorem 4.1, the iterates enjoy exponentialrate uniform convergence to f. This will settle an issue raised in the third paragraph of Section 3 in [15].

References

[1] van der Corput, J. G. Zahlentheoretische Abschätzungen. Math. Ann. 84 (1921), 53–79.

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

[3] Feller, W. An Introduction to Probability Theory and its Applications. Vol. II. Second edition. Wiley, New York, 1971.

[4] Flajolet, Ph. See http://algo.inria.fr/AofA/Research/05-99.html .

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

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

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

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

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

[10] Montgomery, H. L. Ten Lectures on the Interface Between Analytic Number Theory and Harmonic Analysis. CBMS Reg. Conf. Ser. Math. 84, AMS, Providence, R.I., 1994.

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

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

[13] Rudin, W. Functional Analysis. Second edition. McGraw–Hill, New York, 1991.

[14] Schwartz, L. Théorie des Distributions. Second edition. Hermann, Paris, 1966.

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