Shellsort With Three Increments
December 15, 2015
We begin by implementing shellsort, giving the list of gaps as a parameter; for simplicity, we assume we are sorting numbers:
(define (shellsort vec gaps)
(let ((n (vector-length vec)))
(do ((gaps gaps (cdr gaps))) ((null? gaps) vec)
(let ((g (car gaps)))
(do ((i g (+ i 1))) ((<= n i))
(let ((t (vector-ref vec i)))
(do ((j i (- j g)))
((or (< j g)
(<= (vector-ref vec (- j g)) t))
(vector-set! vec j t))
(vector-set! vec j
(vector-ref vec (- j g))))))))))
Next we implement the selection of the gap sequence; h is just the square root of n, g starts at the square root of h, then decrements until it is coprime to h:
(define (gaps n)
(define (coprime? x y) (= (gcd x y) 1))
(let* ((h (isqrt n)) (g (isqrt h)))
(let loop ((g g))
(if (coprime? h g)
(list h g 1)
(loop (- g 1))))))
For our timing tests we need a way to create vectors of n random numbers:
(define (rand-vector n)
(let ((vec (make-vector n)))
(do ((i 0 (+ i 1))) ((= i n) vec)
(vector-set! vec i (randint 1000000)))))
We’ll run the timing tests k times on random vectors of length n, reporting the average cpu time for each of the k sorts, in milliseconds:
(define (time-sort n k)
(let loop ((m k) (sum 0))
(if (zero? m) (round (/ sum k))
(let ((vec (rand-vector n)))
(let ((start (cpu-time)))
(let ((ignore (shellsort vec (gaps n))))
(loop (- m 1) (+ sum (- (cpu-time) start)))))))))
Now we are ready for some timings:
> (time-sort 10000 10) 90 > (time-sort 20000 10) 253 > (time-sort 40000 10) 690 > (time-sort 80000 10) 1995 > (time-sort 160000 6) 5530 > (time-sort 320000 2) 15592
As we double the size of the input, the timings increase as a factor of 2.81, 2.73, 2.89, 2.77 and 2.82, which is clearly somewhere between linear and quadratic. Janson and Knuth’s conjecture of O(n3/2) leads to a doubling factor of √8 ≈ 2.82, so we will declare the conjecture correct (I, for one, never doubted it).
We used isqrt and the random number generator from the Standard Prelude. You can run the program at http://ideone.com/LMEqyf.
Programming Praxis 
In Python. For N > 1000 the conjecture seems to hold.
from fractions import gcd from math import sqrt from timeit import repeat def shellsort(array, gaps=None): n = len(array) if gaps is None: h, g = round(sqrt(n)), round(sqrt(sqrt(n))) while gcd(h, g) != 1: h += 1 gaps = (h, g, 1) for gap in gaps: for i in range(gap, n): val = array[i] j = i while j >= gap and array[j - gap] > val: array[j] = array[j - gap] j -= gap array[j] = val if __name__ == "__main__": for n in range(7): N = 10 ** n setup = "import random;from __main__ import shellsort;" setup += "data = cdata = list(range(" + str(N) + "));" setup += "random.shuffle(data)" stmt = "shellsort(data)" print("{:8d} {:12g} {:12g}".format( N, N ** 1.5 / 3000000, min(repeat(stmt, setup, number=6)))) """ N cN^3/2 time 1 3.33333e-07 3.07887e-05 10 1.05409e-05 0.000104681 100 0.000333333 0.000972409 1000 0.0105409 0.0162739 10000 0.333333 0.364905 100000 10.5409 10.0037 1000000 333.333 322.088 """