Shellsort With Three Increments by Svante Janson and Donald E. Knuth

Abstract. A perturbation technique can be used to simplify and sharpen A. C. Yao’s theorems about the behavior of shellsort with increments (h, g, 1). In particular, when h = Θ(n^{7/15}) and g = Θ(h^{1/5}), the average running time is O(n^{23/15}). The proof involves interesting properties of the inversions in random permutations that have been h-sorted and g-sorted.

Shellsort, also known as the “diminishing increment sort” [7, Algorithm 5.2.1D], puts the elements of an array (X_{0}, …, X_{n−1}) into order by successively performing a straight insertion sort on larger and larger subarrays of equally spaced elements. The algorithm consists of t passes defined by increments (h_{t−1}, …, h_{1}, h_{0}), where h_{0} = 1; the jth pass makes X_{k} ≤ X_{l} whenever l − k = h_{t−j}.

A. C. Yao [11] has analyzed the average behavior of shellsort in the general three-pass case when the increments are (h, g, 1). The most interesting part of his analysis dealt with the third pass, where the running time is O(n) plus a term proportional to the average number of inversions that remain after a random permutation has been h-sorted and g-sorted. Yao proved that if g and h are relatively prime, the average number of inversions remaining is

ψ(h, g)n + Ô(n^{2/3}) ,  (0.1)

where the constant implied by Ô depends on g and h. He gave a complicated triple sum for ψ(h, g), which is too difficult to explain here; we will show that

ψ(h, g) = \frac{1}{2} ∑_{d=1}^{g−1} ∑_{r} \binom{h − 1}{r} (\frac{d}{g})^{r} (1 − \frac{d}{g})^{h−1−r} | r − ⌊\frac{hd}{g}⌋| .  (0.2)

Moreover, we will prove that the average number of inversions after such h-sorting and g-sorting is

ψ(h, g)n + O(g^{3}h^{2}) ,  (0.3)

where the constant implied by O is independent of g, h, and n.

The main technique used in proving (0.3) is to consider a stochastic algorithm A whose output has the same distribution as the inversions of the third pass of shellsort. Then by slightly perturbing the probabilities that define A, we will obtain an algorithm A^{∗} whose output has the expected value ψ(h, g)n exactly. Finally we will prove that the perturbations cause the expected value to change by at most O(g^{3}h^{2}).

Section 1 introduces basic techniques for inversion counting, and section 2 adapts those techniques to a random input model. Section 3 proves that the crucial random variables needed for inversion counting are nearly uniform; then section 4 shows that the leading term ψ(h, g)n in (0.3) would be exact if those variables were perfectly uniform. Section 5 explains how to perturb them so that they are indeed uniform, and section 6 shows how this perturbation yields the error term O(g^{3}h^{2}) of (0.3).

The asymptotic value of ψ(h, g) is shown to be (πh/128)^{1/2}g in section 7. The cost of the third pass in (ch, cg, 1)-shellsort for c > 1 is analyzed in section 8. This makes it possible to bound the total running time for all three passes, as shown in section 9, leading to an O(n^{23/15}) average running time when h and g are suitably chosen.

The bound O(g^{3}h^{2}) in (0.3) may not be best possible. Section 10 discusses a conjectured improvement, consistent with computational experiments, which would reduce the average cost to O(n^{3/2}), if it could be proved.

The tantalizing prospect of extending the techniques of this paper to more than three increments is explored briefly in section 11.

1. Counting inversions. We shall assume throughout this paper that g and h are relatively prime. To fix the ideas, suppose h = 5, g = 3, n = 20, and suppose we are sorting the 2-digit numbers

(X_{0}, X_{1}, …, X_{n−1}) = (03, 14, 15, 92, 65, 35, 89, 79, 32, 38, 46, 26, 43, 37, 31, 78, 50, 28, 84, 19) .

(Cf. [6, Eq. 3.3–(1)].) The first pass of shellsort, h-sorting, replaces this array by

(X′_{0}, X′_{1}, …, X′_{n−1}) = (03, 14, 15, 32, 19, 35, 26, 28, 37, 31, 46, 50, 43, 84, 38, 78, 89, 79, 92, 65) .

The second pass, g-sorting, replaces it by

(X′′_{0}, X′′_{1}, …, X′′_{n−1}) = (03, 14, 15, 26, 19, 35, 31, 28, 37, 32, 46, 38, 43, 65, 50, 78, 84, 79, 92, 89) .

Our task is to study the inversions of this list, namely the pairs k, l for which k < l and X′′_{k} > X′′_{l}.

The result of g-sorting is the creation of g ordered bits X′′_{j} < X′′_{j+g} < X′′_{j+2g} < ⋯ for 0 ≤ j < g, each of which contains no inversions within itself. So the inversions remaining are inversions between different sublists. For example, the 20 numbers sorted above lead to

list 0 = (03, 26, 31, 32, 43, 78, 92) ,
list 1 = (14, 19, 28, 46, 65, 84, 89) ,
list 2 = (15, 35, 37, 38, 50, 79) ;

the inversions between list 0 and list 1 are the inversions of

(03, 14, 26, 19, 31, 28, 32, 46, 43, 65, 78, 84, 92, 89) .

It is well known [7, §5.21] that two interleaved ordered lists of lengths m have ∑_{r=0}^{m−1} |r − s_{r}| inversions, where s_{r} of the elements of the second list are less than the (r + 1)st element of the first list; for example, (03, 14, 26, …, 89) has

|0 − 0| + |1 − 2| + |2 − 3| + |3 − 3| + |4 − 3| + |5 − 5| + |6 − 7| = 4

inversions. If r ≥ s_{r}, the (r+1)st element of the first list is inverted by r−s_{r} elements of the second; otherwise it inverts s_{r} − r of those elements. (We assume that the list elements are distinct.) The same formula holds for interleaved ordered lists of lengths m and m − 1, because we can imagine an infinite element at the end of the second list.

Let Y_{kl} be the number of elements X_{k′} such that k′ ≡ k (mod h) and X_{k′} < X_{l}. The n numbers Y_{ll} for 0 ≤ l < n clearly characterize the permutation performed by h-sorting; and it is not hard to see that the full set of hn numbers Y_{kl} for 0 ≤ k < h and 0 ≤ l < n is enough to determine the relative order of all the X’s.

There is a convenient way to enumerate the inversions that remain after g-sorting, using the numbers Y_{kl}. Indeed, let

J_{kl} = (k mod h + hY_{kl}) mod g .  (1.1)

Then X_{l} will appear in list j = J_{ll} after g-sorting. Let S_{jl} be the number of elements X_{k′} such that X_{k′} < X_{l} and X_{k′} is in list j. The inversions between lists j and j′ depend on the difference |S_{jl} − S_{j′l}| when X_{l} goes into list j.

Given any values of j and j′ with 0 ≤ j < j′ < g, let j_{s} = (j + hs) mod g, and let d be minimum with j_{d} = j′. Thus, d is the distance from j to j′ if we count by steps of h modulo g. Let

H = {j_{1}, j_{2}, …, j_{d}}  (1.2)

be the h numbers between j and j′ in this counting process, and let Q_{l} be the number of indices k such that 0 ≤ k < h and J_{kl} ∈ H. Then we can prove the following basic fact:

Lemma 1. Using the notation above, we have

S_{jl} − S_{j′l} = Q_{l} − ⌊hd/g⌋  (1.3)

for all j, j′, and l with 0 ≤ j < j′ < g and 0 ≤ l < n.

Proof. Since the X’s are distinct, there is a permutation (l_{0}, l_{1}, …, l_{n−1}) of {0, 1, …, n − 1} such that X_{l_{0}} < X_{l_{1}} < ⋯ < X_{l_{n−1}}. We will prove (1.3) for l = l_{t} by induction on t.

Suppose first that l = l_{0}, so that X_{l} is the smallest element being sorted. Then Y_{kl} = 0 for all k, hence J_{kl} = k mod g for 0 ≤ k < h. Also S_{jl} = S_{j′l} = 0. Therefore (1.3) is equivalent in this case to the assertion that precisely ⌊hd/g⌋ elements of the multiset

{0 mod g, 1 mod g, …, (h − 1) mod g}

belong to H.

A clever proof of that assertion surely exists, but what is it? We can at any rate use brute force by assuming first that j = 0. Then the number of solutions to x ≡ hd (mod g) and 0 ≤ x < h is the number of integers in the interval [−hd/g . . −h(d − 1)/g), namely ⌈−h(d − 1)/g⌉ − ⌈−hd/g⌉ = ⌊hd/g⌋ − ⌊h(d − 1)/g⌋. Therefore the assertion for j = 0 follows by induction on d. And once we’ve proved it for some pair j < j′, we can prove it for j + 1 < j′ + 1, assuming that j′ + 1 < g: The value of d stays the same, and the values of j_{1}, j_{2}, …, j_{d} increase by 1 (mod g). So we lose one solution if j_{s} ≡ h − 1 (mod g) for some s with 1 ≤ s ≤ d; we gain one solution if j_{s} ≡ −1 (mod g) for some s. Since j_{s} ≡ h − 1 ⟺ j_{s−1} ≡ −1, the net change is zero unless j_{1} ≡ h − 1 (but then j = g − 1) or j_{d} ≡ −1 (but then j′ = g − 1). This completes the proof by brute force when l = l_{0}.

Suppose (1.3) holds for l = l_{t}; we want to show that it also holds for when l is replaced by l′ = l_{t+1}. The numbers Y_{kl} and Y_{kl′} are identical for all but one value of k, since

Y_{kl′} = Y_{kl} + [l ≡ k (mod h)] .

Thus, the values of J_{kl} and J_{kl′} are the same except that J_{kl} increases by h (mod g) when k ≡ l (mod h). It follows that

Q_{l′} = Q_{l} + [J_{ll} = j] − [J_{ll} = j′] .

This completes the proof by induction on t, since S_{jl′} = S_{jl} + [J_{ll} = j] for all j.

Corollary. Using the notations above, the total number of inversions between lists j and j′ is

∑_{l=0}^{n−1} | Q_{l} − ⌊hd/g⌋ | [J_{ll} = j] .  (1.4)

Proof. This is |S_{jl} − S_{j′l}| = |r − s_{r}| summed over all r such that X_{l} is the (r + 1)st element of list j.

In the example of n = 20 two-digit numbers given earlier, with h = 5, g = 3, j = 0, and j′ = 1, we have d = 2, H = {2, 1},

\begin{matrix}l = & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 & 13 & 14 & 15 & 16 & 17 & 18 & 19 \\ X_{l} = & 03 & 14 & 15 & 92 & 65 & 35 & 89 & 79 & 32 & 38 & 46 & 26 & 43 & 37 & 31 & 78 & 50 & 28 & 84 & 19 \\ Y_{0l} = & 0 & 1 & 1 & 4 & 3 & 1 & 4 & 4 & 1 & 2 & 2 & 1 & 2 & 2 & 1 & 3 & 3 & 1 & 4 & 1 \\ Y_{1l} = & 0 & 0 & 1 & 4 & 3 & 2 & 3 & 3 & 2 & 2 & 2 & 1 & 2 & 2 & 2 & 3 & 2 & 2 & 3 & 1 \\ Y_{2l} = & 0 & 0 & 0 & 4 & 3 & 2 & 4 & 3 & 2 & 2 & 3 & 1 & 2 & 2 & 2 & 3 & 3 & 1 & 4 & 1 \\ Y_{3l} = & 0 & 0 & 0 & 3 & 2 & 1 & 3 & 2 & 0 & 2 & 2 & 0 & 2 & 1 & 0 & 2 & 2 & 0 & 2 & 0 \\ Y_{4l} = & 0 & 0 & 0 & 4 & 3 & 2 & 4 & 4 & 2 & 2 & 3 & 1 & 3 & 2 & 1 & 4 & 3 & 1 & 4 & 0 \\ J_{0l} = & \underline{0} & 2 & 2 & 2 & 0 & \underline{2} & 2 & 2 & 2 & 1 & \underline{1} & 2 & 1 & 1 & 2 & \underline{0} & 0 & 2 & 2 & 2 \\ J_{1l} = & 1 & \underline{1} & 0 & 0 & 1 & 2 & \underline{1} & 1 & 2 & 2 & 2 & \underline{0} & 2 & 2 & 2 & 1 & \underline{2} & 2 & 1 & 0 \\ J_{2l} = & 2 & 2 & \underline{2} & 1 & 2 & 0 & 1 & \underline{2} & 0 & 0 & 2 & 1 & \underline{0} & 0 & 0 & 2 & 2 & \underline{1} & 1 & 1 \\ J_{3l} = & 0 & 0 & 0 & \underline{0} & 1 & 2 & 0 & 1 & \underline{0} & 1 & 1 & 0 & 1 & \underline{2} & 0 & 1 & 1 & 0 & \underline{1} & 0 \\ J_{4l} = & 1 & 1 & 1 & 0 & \underline{1} & 2 & 0 & 0 & 2 & \underline{2} & 1 & 0 & 1 & 2 & \underline{0} & 0 & 1 & 0 & 0 & \underline{1} \\ Q_{l} = & 3 & 4 & 3 & 2 & 4 & 4 & 3 & 4 & 3 & 4 & 5 & 2 & 4 & 4 & 2 & 3 & 4 & 3 & 4 & 3 \end{matrix}

and the underlined values J_{ll} are 0 for l = 0, 3, 8, 11, 12, 14, 15 (accounting for the seven elements in list 0). The inversions between lists 0 and 1 are therefore

|3 − 3| + |2 − 3| + |3 − 3| + |2 − 3| + |4 − 3| + |2 − 3| + |3 − 3| = 4

according to (1.4).

2. Random structures. We obtain a random run of shellsort if we assume that the input array (X_{0}, X_{1}, …, X_{n−1}) is a random point in the n-dimensional unit cube. For each integer l in the range 0 ≤ l < n and for each “time” t in the range 0 ≤ t ≤ 1 we will consider the contribution made by X_{l} to the total number of inversions if X_{l} = t.

Thus, instead of the quantities Y_{kl} and J_{kl} defined in the previous section, we define

Y_{kl}(t) = ∑_{{k′≡k ( mod h) \atop 0≤k′<n}} [X_{k′} < t] ,  (2.1)

J_{kl}(t) = (k mod h + hY_{kl}(t)) mod g .  (2.2)

These equations are almost, but not quite, independent of l, because we assume that X_{l} = t while all other X’s are uniformly and independently random.

For each pair of indices j and j′ with 0 ≤ j < j′ < g, we define H as in (1.2), and we let

Q_{l}(t) = ∑_{k=0}^{h−1} [J_{kl}(t) ∈ H][k ≠ l mod h] .  (2.3)

This definition is slightly different from our original definition of Q_{l}, because we have excluded the term for k = l mod h. However, formula (1.4) remains valid because j ∉ H; when J_{ll} = j, the excluded term is therefore zero.

Notice that, for fixed l, the random variables Y_{kl}(t) for 0 ≤ k < h are independent. Therefore the random variables J_{kl}(t) are independent; and Q_{l}(t) is independent of J_{ll}(t). The average contribution of X_{l} to the inversions between lists j and j′ when X_{l} = t is therefore

W_{jj′l}(t) = Pr [J_{ll}(t) = j] E |Q_{l}(t) − ⌊hd/g⌋|  (2.4)

by (1.4), where probabilities and expectations are computed with respect to (X_{0}, …, X_{l−1}, X_{l+1}, …, X_{n−1}). The average total contribution of X_{l} is obtained by integrating over all values of t:

Lemma 2. Let

W_{jj′l} = ∫_{0}^{1} W_{jj′l}(t) dt .  (2.5)

Then the average grand total number of inversions in the third pass of shellsort is

∑_{{0≤j<j′<g \atop 0≤l<n}} W_{jj′l} .  (2.6)

Our goal is to find the asymptotic value of this sum, by proving that it agrees with the estimate (0.3) stated in the introduction.

3. Near uniformity. The complicated formulas of the previous section become vastly simpler when we notice that each random variable J_{kl}(t) is almost uniformly distributed: The probability that J_{kl}(t) = j is very close to 1/g, for each j, as long as t is not too close to 0 or 1. To prove this statement, it suffices to show that Y_{kl}(t) mod g is approximately uniform, because h is relatively prime to g. Notice that Y_{kl}(t) has a binomial distribution, because it is the sum of approximately n/h independent random 0–1 variables that take the value 1 with probability t.

Lemma 3. If Y has the binomial distribution with parameters (m, t), then

|Pr [Y mod g = j] − \frac{1}{g}| < \frac{1}{g} ϕ_{gm}(t)  (3.1)

for 0 ≤ j < g, where

ϕ_{gm}(t) = 2 ∑_{k=1}^{∞} e^{−8t(1−t)k^{2}m/g^{2}} .  (3.2)

Proof. Let y_{j} = Pr [Y mod g = j], and consider the discrete Fourier transform

ŷ_{k} = ∑_{j=0}^{g−1} ω^{kj}y_{j} = E ω^{kY}

where ω = e^{2πi/g}. We have

ŷ_{k} = ∑_{l=0}^{m} \binom{m}{l}t^{l}(1 − t)^{m−l}ω^{kl} = (ω^{k}t + 1 − t)^{m} ,  (3.3)

and

|ω^{k}t + 1 − t|^{2} = t^{2} + (1 − t)^{2} + t(1 − t)(ω^{k} + ω^{−k})
= 1 − 2t(1 − t)(1 − cos 2πk/g)
= 1 − 4t(1 − t) sin^{2} πk/g .  (3.4)

If 0 ≤ x ≤ π/2 we have sin x ≥ 2x/π; hence, if 0 ≤ k ≤ \frac{1}{2}g,

|ω^{k}t + 1 − t|^{2} ≤ 1 − 16t(1 − t)k^{2}/g^{2} < e^{−16t(1−t)k^{2}/g^{2}} .

And if \frac{1}{2}g < k < g we have |ŷ_{k}| = |ŷ_{g−k}|. Therefore

∑_{k=1}^{g−1} |ŷ_{k}| ≤ 2 ∑_{k=1}^{g/2} e^{−8t(1−t)k^{2}m/g^{2}} < ϕ_{gm}(t) .  (3.5)

The desired result follows since

y_{j} = \frac{1}{g} ∑_{k=0}^{g−1} ω^{−kj}ŷ_{k}

and thus

|y_{j} − \frac{1}{g}| = |\frac{1}{g} ∑_{k=1}^{g−1} ω^{−kj}ŷ_{k} | ≤ \frac{1}{g} ∑_{k=1}^{g−1} |ŷ_{k}| .

Corollary. We have

|Pr [J_{kl}(t) = j] − \frac{1}{g}| < \frac{1}{g} ϕ(t)  (3.6)

for 0 ≤ k < h, where

ϕ(t) = \{ \begin{matrix}2 ∑_{k=1}^{∞} e^{−4t(1−t)k^{2}n/g^{2}h} , & if n ≥ 4h; \\ g, & if n < 4h. \end{matrix}  (3.7)

Proof. Each variable Y_{kl}(t) in (2.1) for 0 ≤ k < h has the binomial distribution with parameters (m, t), where if n ≥ 4h

m = ⌈(n − k)/h⌉ − [k = l mod h] ≥ \frac{n}{h} − 2 ≥ \frac{n}{2h} .

Now J_{kl}(t) = j if and only if Y_{kl}(t) has a certain value mod g. The case n < 4h is trivial.

4. Uniformity. Let’s assume now that, for given l and t, the random variables J_{kl}(t) have a perfectly uniform distribution. Since the variables J_{kl}(t) are independent for 0 ≤ k < h, this means that

Pr [J_{0l}(t) = j_{0}, J_{1l}(t) = j_{1}, …, J_{(h−1)l}(t) = j_{h−1}] = \frac{1}{g^{h}}  (4.1)

for all h-tuples (j_{0}, j_{1}, …, j_{h−1}).

In such a case the random variable Q_{l}(t) defined in (2.3) is the sum of h − 1 independent indicator variables, each equal to 1 with probability d/g because H has d elements. Hence Q_{l}(t) has the binomial distribution with parameters (h − 1, d/g), and it is equal to r with probability

\binom{h − 1}{r} (\frac{d}{g})^{r} (1 − \frac{d}{g})^{h−1−r} .  (4.2)

Let W_{jj′l}^{∗}(t) be the value of W_{jj′l}(t) under the assumption of uniformity (see (2.4)). Thus W_{jj′l}^{∗}(t) is independent of t, and we let W_{jj′l}^{∗} = W_{jj′l}^{∗}(t) in accordance with (2.5). Then

W_{jj′l}^{∗} = \frac{1}{g} ∑_{r=0}^{h−1} \binom{h − 1}{r} (\frac{d}{g})^{r} (1 − \frac{d}{g})^{h−1−r} | r − ⌊\frac{hd}{g}⌋| .  (4.3)

For given values of d and j, the index j′ = (j + hd) mod g is at distance d from j. Suppose a(d) of these pairs (j, j′) have j < j′. Then g − a(d) of them have j > j′, and a(g − d) = g − a(d) since j is at distance g −d from j′. The sum of (4.3) over all j < j′ is therefore independent of a(d):

∑_{0≤j<j′<g} W_{jj′l}^{∗} = ∑_{d=1}^{g−1} \frac{a(d)}{g} ∑_{r=0}^{h−1} \binom{h − 1}{r} (\frac{d}{g})^{r} (1 − \frac{d}{g})^{h−1−r} | r − ⌊\frac{hd}{g}⌋|

= ∑_{d=1}^{g−1} \frac{a(g − d)}{g} ∑_{r=0}^{h−1} \binom{h − 1}{h − 1 − r} (\frac{g − d}{g})^{h−1−r} (1 − \frac{g − d}{g})^{r}

× | h − 1 − r − ⌊\frac{h(g − d)}{g}⌋|

= \frac{1}{2} ∑_{d=1}^{g−1} ∑_{r=0}^{h−1} \binom{h − 1}{r} (\frac{d}{g})^{r} (1 − \frac{d}{g})^{h−1−r} | r − ⌊\frac{hd}{g}⌋| .

(We have used the fact that ⌊h(g − d)/g⌋ = h − 1 − ⌊hd/g⌋ when hd/g is not an integer.) But this is just the quantity ψ(h, g) in (0.2), for each value of l. We have proved

Lemma 4. If we assume that the variables J_{kl}(t) have exactly the uniform distribution, the quantity (2.6) is exactly ψ(h, g)n.

5. Perturbation. To complete the proof of (0.3), we use a general technique applicable to the analysis of many algorithms: If a given complicated algorithm A almost always has the same performance characteristics as a simpler algorithm A^{∗}, then the expected performance of A is the same as the performance of A^{∗} plus an error term based on the cases where A and A^{∗} differ. (See, for example, the analysis in [8], where this “principle of negligible perturbation” is applied to a nontrivial branching process.)

In the present situation we retain the (n − 1)-dimensional probability space (X_{0}, …, X_{l−1}, t, X_{l+1}, …, X_{n−1}) on which the random variables J_{kl}(t) were defined in (2.2), and we define a new set of random variables J_{kl}^{∗}(t) on the same space, where J_{kl}^{∗}(t) has exactly a uniform distribution on {0, 1, …, g − 1}. This can be done in such a way that J_{kl}(t) = J_{kl}^{∗}(t) with high probability.

More precisely, when l and t are given, J_{kl}(t) depends only on the variables X_{k′} with k′ ≡ k (mod h) and k′ ≠ l. The unit cube on these variables is partitioned into g parts P_{0}, P_{1}, …, P_{g−1} such that J_{kl}(t) = j when the variables lie in P_{j}; the volume of P_{j} is Pr [J_{kl}(t) = j]. We will divide each P_{j} into g sets P′_{j0}, P′_{j1}, …, P′_{j(g−1)}, and define J_{kl}^{∗}(t) = i on P′_{ji}. This subdivision, performed separately for each k, will yield independent random variables J_{0l}^{∗}(t), J_{1l}^{∗}(t), …, J_{(h−1)l}^{∗}(t). We will show that the subdivision can be done in such a way that

Pr [J_{kl}^{∗}(t) = j] = 1/g ,  (5.1)
Pr [J_{kl}^{∗}(t) ≠ J_{kl}(t)] < ϕ(t) ,  (5.2)

for 0 ≤ j < g and 0 ≤ k < h. Thus, we will have perturbed the values of J_{kl}(t) with low probability when ϕ(t) is small.

The following construction does what we need, and more:

Lemma 5. Let p_{1}, …, p_{m} and p_{1}^{∗}, …, p_{m}^{∗} be nonnegative real numbers with p_{1} + ⋯ + p_{m} = p_{1}^{∗} + ⋯ + p_{m}^{∗} = 1. Then there are nonnegative reals p′_{ij} for 1 ≤ i, j ≤ m such that

p_{i} = ∑_{j=1}^{m} p′_{ij} ,  (5.3)

p_{j}^{∗} = ∑_{i=1}^{m} p′_{ij} ,  (5.4)

and

∑_{i≠j} p′_{ij} = 1 − ∑_{j} p′_{jj} = \frac{1}{2} ∑_{j} |p_{j} − p_{j}^{∗}|.  (5.5)

Proof. This is a special case of “maximal coupling” in probability theory [5; 9, §III.14]; it can be proved as follows.

Let p′_{jj} = min(p_{j}, p_{j}^{∗}), and observe that
∑_{j} p′_{jj} = ∑_{j} min(p_{j}, p_{j}^{∗}) = ∑_{j} \frac{1}{2}(p_{j} + p_{j}^{∗} − |p_{j} − p_{j}^{∗}|) = 1 − \frac{1}{2} ∑_{j} |p_{j} − p_{j}^{∗}| .  (5.6)

The existence of nonnegative p′_{ij}, i ≠ j, such that (5.3) and (5.4) hold follows from the max flow–min cut theorem [4]: Consider a network with a source s, a sink t, and 2m nodes v_{1}, …, v_{m}, v_{1}^{∗}, …, v_{m}^{∗}; the edges are sv_{j} with capacity p_{j}−p′_{jj}, v_{j}^{∗}t with capacity p_{j}^{∗}−p′_{jj}, and v_{i}v_{j}^{∗} with infinite capacity.

6. The effect of perturbation. When independent random variables J_{kl}^{∗}(t) have been defined satisfying (5.1) and (5.2), we can use them to define Q_{l}^{∗}(t) as in (2.3) and W_{jj′l}^{∗}(t) as in (2.4). This value W_{jj′l}^{∗}(t) has already been evaluated in (4.3); we want now to use the idea of perturbation to see how much W_{jj′l}(t) can differ from W_{jj′l}^{∗}(t).

Since Q_{l}(t) = O(h) and
|Q_{l}(t) − Q_{l}^{∗}(t)| ≤ ∑_{k=0}^{h−1} [J_{kl}(t) ≠ J_{kl}^{∗}(t)] ,  (6.1)

we have

|W_{jj′l}(t) − W_{jj′l}^{∗}(t)| = |(Pr [J_{ll}(t) = j] − Pr [J_{ll}^{∗}(t) = j]) E |Q_{l}(t) − ⌊hd/g⌋|
+ Pr [J_{ll}^{∗}(t) = j](E |Q_{l}(t) − ⌊hd/g⌋| − E |Q_{kl}^{∗}(t) − ⌊hd/g⌋|)|

< \frac{1}{g} ϕ(t)O(h) + \frac{1}{g} ∑_{k=0}^{h−1} Pr [J_{kl}(t) ≠ J_{kl}^{∗}(t)]

= O (\frac{h}{g}) ϕ(t) .  (6.2)

(We assume that J_{kl}^{∗}(t) = J_{(k mod h)l}^{∗}(t) when k ≥ h.)

To complete our estimate we need to integrate this difference over all t.

Lemma 6. ∫_{0}^{1} ϕ(t) dt = O(g^{2}h/n).

Proof. The case n < 4h is trivial. Otherwise we have

∫_{0}^{1} ϕ(t) dt = 2 ∫_{0}^{1/2} ϕ(t) dt

< 4 ∫_{0}^{1/2} ∑_{k=1}^{∞} e^{−2tk^{2}n/g^{2}h} dt

< 4 ∫_{0}^{∞} ∑_{k=1}^{∞} e^{−2tk^{2}n/g^{2}h} dt

= 4 ∑_{k=1}^{∞} \frac{g^{2}h}{2k^{2}n} = \frac{π^{2}}{3} \frac{g^{2}h}{n} .

Theorem 1. The average number of inversions remaining after h-sorting and then g-sorting a random permutation of n elements, when h is relatively prime to g, is ψ(h, g)n + O(g^{3}h^{2}), where ψ(h, g) is given by (0.2).

Proof. By (6.2) and Lemmas 2, 4, and 6, the average is ψ(h, g)n plus

\begin{aligned}∑_{{0≤j<j′<g \atop 0≤l<n}} ∫_{0}^{1}(W_{jj′l}(t) − W_{jj′l}^{∗}(t)) dt & = O(g^{2}n)O(h/g) ∫_{0}^{1} ϕ(t) dt \\ & = O(g^{3}h^{2}) . \end{aligned}

Notice that the proof of this theorem implicitly uses Lemma 5 for each choice of l and t, without requiring any sort of continuity between the values of J_{kl}^{∗}(t) as t varies. We could have defined J_{kl}^{∗}(t) in a continuous fashion; indeed, the random variables [X_{k} < t] partition the (n − 1)- cube into 2^{n−1} subrectangles in each of which J_{kl}(t) has a constant value, so we could define J_{kl}^{∗}(t) over (n − 1)-dimensional rectangular prisms with smooth transitions as a function of t. But such complicated refinements are not necessary for the validity of the perturbation argument.

7. Asymptotics. Our next goal is to estimate ψ(h, g) when h and g are large. Notice that

ψ(h, g) = \frac{1}{2} ∑_{d=1}^{g−1} E |Z(h − 1, d/g) − ⌊\frac{hd}{g}⌋|  (7.1)

where Z(m, p) has the binomial distribution with parameters m and p. The mean of Z(h − 1, d/g) is (h−1)d/g = ⌊hd/g⌋+O(1), and the variance is (h−1)d(g −d)/g^{2}. If we replace Z by a normally distributed random variable with this same mean and variance, the expected value of |Z − ⌊hd/g⌋| is approximately (2π)^{−1/2} ∫_{−∞}^{∞} |t|e^{−t^{2}/2} dt = 2/\sqrt{2π} times the standard deviation, so (7.1) will be approximately

\frac{1}{g} \sqrt{\frac{h}{2π}} ∑_{d=1}^{g−1} \sqrt{d(g − d)} .  (7.2)

The detailed calculations in the remainder of this section justify this approximation and provide a rigorous error bound.

Lemma 7. If Z has the binomial distribution with parameters (m, p), and ⌊mp⌋ ≤ a ≤ ⌈mp⌉, then

E |Z − a| = \sqrt{\frac{2p(1 − p)m}{π}} + O(\frac{1}{\sqrt{mp(1 − p)}}) .  (7.3)

Proof. Consider first the case a = mp. By a formula of De Moivre [1, page 101] and Poincaré [10, pages 56–60], see Diaconis and Zabell [2],

E |Z − mp| = 2⌈mp⌉\binom{m}{⌈mp⌉}p^{⌈mp⌉}(1 − p)^{m+1−⌈mp⌉}.  (7.4)

In order to prove (7.3) in this case we may assume that p ≤ 1/2, since |Z −mp| = |m−Z −m(1−p)|. Moreover, we may assume that mp > 1 since (7.3) otherwise is trivial. Then, a routine application of Stirling’s approximation shows that

E |Z − mp| = \sqrt{\frac{2p(1 − p)m}{π}} exp (O(\frac{1}{mp})) .  (7.5)

Next observe that if ⌊mp⌋ ≤ a ≤ ⌈mp⌉, we have

E |Z − a| = E |Z − mp| + (mp − a)(1 − 2 Pr [Z ≤ mp]) .  (7.6)

Since Pr [Z ≤ mp] = \frac{1}{2} + O((mp(1 − p))^{−1/2}), for example by the Berry–Esseen estimate of the error in the central limit theorem [3, §XVI.5], the result follows.

Corollary. The asymptotic value of ψ(h, g) is

ψ(h, g) = \sqrt{\frac{πh}{128}} g + O(g^{−1/2}h^{1/2}) + O(gh^{−1/2}) .  (7.7)

Proof. Since ⌊hd/g⌋ ≤ ⌊(h + 1)d/g⌋ ≤ ⌊(hd + g − 1)/g⌋ = ⌈hd/g⌉, Lemma 7 yields

ψ(h + 1, g) = \frac{1}{2} ∑_{d=1}^{g−1} E |Z(h, d/g) − ⌊\frac{(h + 1)d}{g}⌋|

= ∑_{d=1}^{g−1} (\sqrt{\frac{h}{2π} \frac{d}{g}(1 − \frac{d}{g})} + O ((h\frac{d}{g}(1 − \frac{d}{g}))^{−1/2}))

= \sqrt{\frac{h}{2π}} ∑_{d=1}^{g−1} \sqrt{\frac{d}{g}(1 − \frac{d}{g})} + O(gh^{−1/2}) .

And Euler’s summation formula with f(x) = \sqrt{(x/g)(1 − x/g)} tells us that

\begin{aligned}∑_{d=1}^{g−1} f(d) & = ∫_{1}^{g−1} f(x) dx + \frac{1}{2}f(1) + \frac{1}{2}f(g − 1) + \frac{1}{12}f′(g − 1) − \frac{1}{12}f′(1) − R \\ & = g ∫_{0}^{1} \sqrt{t(1 − t)} dt + O(g^{−1/2}) = \frac{πg}{8} + O(g^{−1/2}) \end{aligned}

because

|R| = |∫_{1}^{g−1} \frac{B_{2}(x mod 1)}{2}f′′(x) dx| ≤ \frac{1}{12} ∫_{1}^{g−1}|f′′(x)| dx = \frac{1}{12}f′(1) − \frac{1}{12}f′(g − 1) .

The error term is thus O(g^{1/2}) when h = g^{2} + 1; for example, we have

h  g  ψ(h, g)  \sqrt{πh/128} g  difference/\sqrt{g}
901  30  140.018  141.076  0.1933
1601  40  249.539  250.741  0.1900
2501  50  390.412  391.739  0.1877

8. Common factors. Now let’s consider the behavior of shellsort with increments (ch, cg, 1), where c is an integer > 1. It is easy to see that the first two passes are equivalent to the first two passes of (h, g, 1) shellsort on c independent subarrays (X_{a}, X_{a+c}, X_{a+2c}, …) of size ⌈(n − a)/c⌉ for 0 ≤ a < c. The inversions that remain are the ψ(h, g)n + O(g^{3}h^{2}c) inversions within these subarrays, plus “cross-inversions” between \binom{c}{2} pairs of subarrays.

Yao [11, Theorem 2] proved that the average number of cross-inversions is \frac{1}{8} \sqrt{πc} (1−c^{−1})n^{3/2}+ O(cghn). The following lemma improves his error term slightly.

Lemma 8. The average number of cross-inversions after ch-sorting and cg-sorting is

\frac{1}{8} \sqrt{πc} (1 − c^{−1})n^{3/2} + O(cgh^{1/2}n) + O(c^{2}g^{3}h^{2}) .  (8.1)

Proof. Let’s consider first the process of h-sorting and g-sorting two independent arrays (X_{0}, X_{1}, …, X_{n−1}) and (X̂_{0}, X̂_{1}, …, X̂_{n−1}), then interleaving the results to obtain (X′′_{0}, X̂′′_{0}, X′′_{1}, X̂′′_{1}, … X′′_{n−1}, X̂′′_{n−1}). The cross inversions are then the pairs {X′′_{l}, X̂′′_{l′}} where either X′′_{l} > X̂′′_{l′} and l ≤ l′ or X′′_{l} < X̂′′_{l′} and l > l′.

Recasting this process in the model of section 2 above, we assume that X_{l} = t, while the other 2n−1 variables (X_{0}, …, X_{l−1}, …, X_{n−1}, X̂_{0}, …, X̂_{n−1}) are independent and uniformly distributed between 0 and 1. We define

Y_{kl}(t) = ∑_{{k′≡k ( mod h) \atop 0≤k′<n}} [X_{k′} < t] , Ŷ_{kl}(t) = ∑_{{k′≡k ( mod h) \atop 0≤k′<n}} [X̂_{k′} < t]  (8.2)

as in (2.1). The elements of each array are divided into h subarrays by h-sorting, and the elements < t have Y_{kl}(t) and Ŷ_{kl}(t) elements in the kth subarrays. Then g-sorting will form g lists, with

L_{jl}(t) = ∑_{k=0}^{h−1} ⌈\frac{Y_{kl}(t) − a_{kj}}{g} ⌉  (8.3)

elements < t in the jth list of the first array, where a_{kj} ∈ {0, 1, …, g − 1} is given by k + a_{kj}h ≡ j (mod g). Similarly, there will be

L̂_{jl}(t) = ∑_{k=0}^{h−1} ⌈\frac{Ŷ_{kl}(t) − a_{kj}}{g} ⌉  (8.4)

elements < t in the jth list of the second. Element X_{l} = t of the first array will go into list j = J_{ll}(t) as before, where J_{kl}(t) is defined in (2.2). The number of cross-inversions between this element and the elements of the second array will then be

V_{l}(t) = ∑_{j′=0}^{g−1} |L̂_{j′l}(t) − L_{jl}(t) − [j′ < j]| .  (8.5)

The average total number of cross-inversions is the sum of E V_{l}(t) over all l, integrated for 0 ≤ t ≤ 1.

We know from Lemma 3 that the numbers Y_{kl}(t) mod g have approximately a uniform distribution. Therefore

⌈\frac{Y_{kl}(t) − a_{kj}}{g} ⌉ = \frac{Y_{kl}(t) − a_{kj} + R_{jkl}(t)}{g}

where R_{jkl}(t) is approximately uniform on {0, 1, …, g − 1}. It follows that

L_{jl}(t) = \frac{Z_{l}(t)}{g} + ∑_{k=0}^{h−1} (\frac{R_{jkl}(t) − a_{kj}}{g} ) ,  (8.6)

where

Z_{l}(t) = ∑_{k=0}^{h−1} Y_{kl}(t)

is the total number of elements in the first array that are < t.

Since R_{jkl}(t) depends on Y_{kl}(t) mod g only, or equivalently on J_{kl}(t), we may use the perturbed truly uniform random variables J_{kl}^{∗}(t) in section 5 (or repeat the argument there with R_{jkl}(t)) and construct random variables R_{jkl}^{∗}(t) that are uniform on {0, 1, …, g − 1} and satisfy Pr [R_{jkl}^{∗}(t) ≠ R_{jkl}(t)] < ϕ(t); moreover, the variables R_{jkl}^{∗}(t) are independent for 0 ≤ k < h and fixed j and l. Consequently

E |R_{jkl}^{∗}(t) − R_{jkl}(t)| ≤ g Pr [R_{jkl}^{∗}(t) ≠ R_{jkl}(t)] < gϕ(t) .  (8.7)

By independence and the fact that E R_{jkl}^{∗}(t) = (g − 1)/2,

E (∑_{k=0}^{h−1} R_{jkl}^{∗}(t) − h(g − 1)/2)^{2} = ∑_{k=0}^{h−1} E (R_{jkl}^{∗}(t) − (g − 1)/2)^{2} < hg^{2},

which by the Cauchy–Schwarz inequality yields

E |∑_{k=0}^{h−1} R_{jkl}^{∗}(t) − h(g − 1)/2| < \sqrt{h}g.  (8.8)

Let W_{jl} = \frac{1}{g}(∑_{k=0}^{h−1} R_{jkl}(t) − h(g − 1)/2) and b_{j} = \frac{1}{g}(h(g − 1)/2 − ∑_{k=0}^{h−1} a_{kj}); then
L_{jl}(t) = \frac{Z_{l}(t)}{g} + W_{jl} + b_{j},  (8.9)

where by (8.7) and (8.8)

E |W_{jl}(t)| < \sqrt{h} + hϕ(t).

A similar argument shows that

L̂_{jl}(t) = \frac{Ẑ_{l}(t)}{g} + Ŵ_{jl} + b_{j}.

Hence

V_{l}(t) = ∑_{j′=0}^{g−1} (\frac{|Ẑ_{l}(t) − Z_{l}(t)|}{g} + O(|W_{j′l}| + |Ŵ_{j′l}| + 1))

and

E V_{l}(t) = E |Ẑ_{l}(t) − Z_{l}(t)| + O(g\sqrt{h}) + O(gh)ϕ(t) .  (8.10)

The quantity |Ẑ_{l}(t) − Z_{l}(t)| is just what we would get if we were counting the cross-inversions between two fully sorted arrays that have been interleaved. Therefore

∫_{0}^{1} ∑_{l=0}^{n−1} E |Ẑ_{l}(t) − Z_{l}(t)| dt

must be the average number of inversions of a random 2-ordered permutation of 2n elements; this, according to Douglas H. Hunt in 1967, is exactly n2^{2n−2}/\binom{2n}{n} [7, exercise 5.2.1–14]. Since \binom{2n}{n} = (1 + O(1/n))4^{n}/\sqrt{πn}, we obtain the desired total

∫_{0}^{1} E ∑_{l=0}^{n−1} V_{l}(t) dt = \frac{\sqrt{π} n^{3/2}}{4} + O(gh^{1/2}n) + O(g^{3}h^{2})  (8.11)

by Lemma 6. Similarly, the same result holds for two arrays of different sizes n + O(1).

Lemma 8 follows if we replace n by n/c + O(1) in (8.11) and multiply by \binom{c}{2}.

9. The total cost. So far we have been considering only the number of inversions removed during the third pass of a three-pass shellsort. But the first two passes can be analyzed as in Yao’s paper [11]:

Theorem 2. Let g and h be relatively prime and let c be a positive integer. The average number of inversions removed when (ch, cg, 1)-shellsort is applied to a random n-element array is

\frac{n^{2}}{4ch} + O(n)  (9.1)

on the first pass,

\frac{1}{8g} \sqrt{\frac{π}{ch}} (h − 1)n^{3/2} + O(hn)  (9.2)

on the second, and

ψ(h, g)n + \frac{1}{8} \sqrt{\frac{π}{c}} (c − 1)n^{3/2} + O((c − 1)gh^{1/2}n) + O(c^{2}g^{3}h^{2})  (9.3)

on the third.

Proof. The first pass removes an average of \frac{1}{4}(n/ch + O(1))^{2} inversions from ch subarrays of size ⌊n/ch⌋ or ⌈n/ch⌉; this proves (9.1). The second pass is equivalent to the second pass of (h, g, 1)- shellsort on c independent subarrays of sizes ⌊n/c⌋ or ⌈n/c⌉. Equation (9.3) is Lemma 8. So the theorem will follow if we can prove (9.2) in the case c = 1. And that case follows from [11, equation (32)], with the O(n) term replaced by O(n/kh) in the notation of that paper. (See also [7, second edition, exercise 5.2.1–40.)

Corollary. If h = Θ(n^{7/15}), g = Θ(n^{1/5}), and gcd(g, h) = 1, the running time of (h, g, 1)-shellsort is O(n^{23/15}).

Proof. The first pass takes time O(n^{2−7/15}), by (9.1); the second takes O(n^{3/2+7/30−1/5}) + O(n^{1+7/15}), by (9.2); and the third takes O(n^{1+1/5+7/30}) + O(n^{3/5+14/15}) by (7.6) and (9.3).

10. Two conjectures. Our estimate O(g^{3}h^{2}) for the difference between ψ(h, g)n and the average number of third-pass inversions may not be the best possible. In fact, the authors conjecture that the difference is at most O(g^{3}h^{3/2}). This sharper bound may perhaps follow from methods analogous to those in the proof of Lemma 8.

If such a conjecture is valid, the running time of (h, g, 1)-shellsort will be O(n^{3/2}) when h ≈ n^{1/2} and g ≈ n^{1/4}. A computer program was written to test this hypothesis by applying (h, g, 1)-shellsort to random arrays of n elements with h = g^{2} + 1 and n = g^{2}h = g^{4} + g^{2}. The following empirical results were obtained, to three significant figures:

g  inversions  ψ(h, g)n  g inversions/10^{5} ψ(h, g)n/10^{5}
1  0 ± 0  0  17 36.6 ± 2.36/32  37.3
2  7.12 ± 2.09/100  7.5  18 51.7 ± 3.35/32  52.6
3  94.4 ± 13.6/100  98.3  19 71.5 ± 4.81/32  72.9
4  563 ± 59.1/100  581  20 97.3 ± 6.14/10  99.2
5  2210 ± 195/100  2280  21 130 ± 8.93/10  133
6  6740 ± 560/100  6910  22 174 ± 12.3/10  176
7 17200 ± 1300/100  17600  23 226 ± 14.0/10  230
8 38600 ± 2820/100  39500  24 291 ± 16.8/10  297
9 78900 ± 5670/100  80600  25 368 ± 23.7/10  380
10 149000 ± 10600/100 152000  26 475 ± 29.1/10  480
11 265000 ± 17200/32  271000  27 595 ± 39.0/10  603
12 447000 ± 30300/32  458000  28 735 ± 44.9/10  750
13 727000 ± 49300/32  742000  29 922 ± 52.1/10  926
14 1140000 ± 75400/32 1160000  30 1110 ± 74.0/10  1140
15 1730000 ± 116000/32 1760000  31 1370 ± 97.9/10  1380
16 2530000 ± 166000/32 2590000  32 1650 ± 101/10  1670

(The inversion counts are given here in the form μ ± σ/\sqrt{r}, where μ and σ are the empirical mean and standard derivation in r independent trials. For example, 10000 trials were made when g ≤ 10, but only 100 trials were made when g ≥ 20.) Both mean and standard derivation seem to be growing proportionately to g^{6} ≈ n^{3/2}, with σ ≈ μ/15 for g ≥ 10.

These data suggest also another conjecture, that the average number of inversions is ≤ ψ(h, g)n when h and g are relatively prime. Indeed, the deviations from uniformity between A and A^{∗} should tend to cause fewer inversions, because A forces the balance condition Y_{kl}(1) = n/h + O(1) for all k and l. This second conjecture obviously implies running time Θ(n^{3/2}) when h = Θ(n^{1/2}) and g = Θ(n^{1/4}).

11. More than three increments? It may be possible to extend this analysis to (h, g, f, 1)- shellsort, by analyzing the following stochastic algorithm. “Initialize two sets of counters (I_{0}, I_{1}, …, I_{g−1}) and (J_{0}, J_{1}, …, J_{h−1}) by setting I_{j} ← j mod f and J_{k} = k mod g for all j and k. Then execute the following procedure n times: Choose a random k in the range 0 ≤ k < h. Set j ← J_{k} and i ← I_{j}; then set J_{k} ← (J_{k} + h) mod g and I_{j} ← (I_{j} + g) mod f.”

Consider the transition from l = l_{t} to l′ = l_{t+1} in the proof of Lemma 1. When elements enter the array in increasing order, the choice of k represents the subarray that will contain a new element X during the h-sort; then X goes into list j during the g-sort, and into list i during the f-sort. We can therefore obtain the contribution of X to the inversions between lists i and i′ for i < i′ < f, by considering a state P_{l} obtained from the I table just as Q_{l} was obtained from the J table in Lemma 1.

References

[1] Abraham De Moivre, Miscellanea Analytica de Seriebus et Quadraturis, (London: J. Tonson and J. Watts, 1730).

[2] Persi Diaconis and Sandy Zabell, “Closed form summation for classical distributions: Variations on a theme of De Moivre,” Statistical Science 6 (1991), 284–302.

[3] William Feller, An Introduction to Probability Theory and Its Applications 2 (New York: Wiley, 1966).

[4] L. R. Ford, Jr., and D. R. Fulkerson, “Maximal flow through a network,” Canadian Journal of Mathematics 8 (1956), 399–404.

[5] Sheldon Goldstein, “Maximal coupling,” Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete 46 (1979), 193–204.

[6] Donald E. Knuth, Seminumerical Algorithms, Volume 2 of The Art of Computer Programming (Reading, Massachusetts: Addison–Wesley, 1969). Second edition, 1981.

[7] Donald E. Knuth, Sorting and Searching, Volume 3 of The Art of Computer Programming (Reading, Massachusetts: Addison–Wesley, 1973). Second edition, 1997.

[8] Donald E. Knuth, Rajeev Motwani, and Boris Pittel, “Stable husbands,” Random Structures and Algorithms 1 (1990), 1–14.

[9] Torgny Lindvall, Lectures on the Coupling Method (New York: Wiley, 1992).

[10] Henri Poincaré, Calcul des Probabilités (Paris: Georges Carré, 1896).

[11] Andrew Chi-Chih Yao, “An analysis of (h, k, 1)-Shellsort,” Journal of Algorithms 1 (1980), 14–50.

Authors’ addresses:

Svante Janson, Department of Mathematics, Uppsala University, P.O.Box 480, 75106 Uppsala,

Sweden; svante.janson@math.uu.se

Donald E. Knuth, Computer Science Department, Gates Building 4B, Stanford University, Stanford

CA 94305–9045 USA; http://www-cs-faculty.stanford.edu/~knuth