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(n^3/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.

Posted by programmingpraxis

Filed in Exercises

1 Comment »

One Response to “Shellsort With Three Increments”

Paul said

December 15, 2015 at 1:16 PM

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
    """

Programming Praxis