Quicksort is broken. Is it worth fixing?

Introduction

1. We are told quicksort is good. We are told it usually takes time proportional to N*logN (where N is the number of input records), and indeed it does, for randomly ordered data.
2. Quicksort's worst case takes time proportional to N*N, though that doesn't happen at all often in practice. It's also reputed to do badly on already sorted data. This isn't fair. For some of the really early versions of Quicksort the worst case was an array in sorted order. And the mud stuck.
   My own experiments have shown that it in fact does better for mostly sorted data (we shall see what kinds of ordering upset it later; uniformly ascending or descending data speed it up, if we always choose the middle element of a partition to divide it into smaller partitions!). It's just that other sort algorithms do even better.
3. Duplicated records are the real worry (at least... as Quicksort is usually written - but there are two separate fixes for that given later in this article, one a minor change, the other a near-total rewrite that gears Quicksort for handling duplicated data... at a price).
4. No matter how we choose the pivot element, some input orders will make Quicksort perform very badly indeed. I'll provide an "adversarial" algorithm that will shuffle pre-sorted data to upset most modern Quicksort implementations later in the article (at first glance it's more of curiosity value than of practical use, but there's method to my madness, as you shall see...).
5. Most quicksort implementations use recursion to sort both partitions. This is incorrect. Recursion should only be used for the smaller partition, if at all. If it is used for both partitions, the maximum stack depth that might be required is N-1. In other words, if you call a Quicksort implementation that calls itself recursively to sort both partitions, and supply it a large input dataset, it may run out of stack and crash. [X] [link to an example which illustrates the stack problem for a 16-element array, with "middle element pivot" selection.] [second link to a corrected quicksort algorithm that does NOT have this fault]
6. Procedural recursion is not necessary at all, in Quicksort. If we are prepared to set aside a little memory to store indices of the partition elements used thus far, we can do without recursion altogether. But it turns out this slows the sort down.[X] [link to algorithm]
7. Let's get lateral. Is it worth randomizing the input array order?
   Possibly. Sometimes. But the overhead required to randomize it may outweigh the benefits.
8. But Quicksort has another problem. It is not a stable sort. An explanation is in order. When, say, a user has a sorted array, and sorts by another column, does the sort algorithm always preserve the relative order of the elements that differ for the previous sort key but have the same value in the new sort key? Most users would expect this! A sort algorithm which always preserves the order of elements (which do not differ in the new sort key) is called a "stable sort". One which does not is called an "unstable sort". This can be illustrated with a table containing four records!
   
   

• With an unstable sort this might happen (pink shows the shuffling)	Sorted by name number name 1 Cat 1 Fox 3 Rabbit 2 Rabbit	User sorts by number number name 1 Fox 1 Cat 2 Rabbit 3 Rabbit	User sorts by name again number name 1 Cat 1 Fox 3 Rabbit 2 Rabbit
  With a stable sort the user will always see this	Sorted by name number name 1 Cat 1 Fox 3 Rabbit 2 Rabbit	User sorts by number number name 1 Cat 1 Fox 2 Rabbit 3 Rabbit	User sorts by name again number name 1 Cat 1 Fox 2 Rabbit 3 Rabbit

  
  
  

• For those of you much better at math than me (most of you!), stable sorts are transitive (that is SortonA(SortonB(Data)) is equivalent to SortOnAAndB(Data).
  It turns out that most implementations of Quicksort are unstable sorts. I'll give an example of a stable implementation of Quicksort later in this article. There's also some code illustrating how you can "stabilize" the result of an unstable sort, with an extra pass over the data.

Sorting small partitions

Partition Selection

Reducing Recursion

1. If current partition size=2 then swap if necessary: exit (Note: this can be dispensed with if we have M=3 or more, see above).
2. Otherwise, choose a partition element, p
3. For each *other* element in the partition, check it against the value of element p.
4. As we find pairs of elements (on each side of the partition) that are on the wrong side, swap them. When we have finished this process we will have two sub-partitions.
5. Recursively sort the left-hand sub-partition
6. Recursively sort the right-hand sub-partition

• Recursively sort the smaller sub-partition
• Set the current partition to the larger partition
• Loop back to the first step

Const minSubFile = 9

Sub ISortAlgorithm_Sort(d As Variant, ByVal l As Long, ByVal r As Long, a As Variant)
    If l + minSubFile <= r Then 
        Quicksort d, l, r, a                'Don't bother quicksorting small arrays
    End If
    InsertionSort d, l, r, a                'Straight insertion sort
End Sub

Sub Quicksort(d As Variant, ByVal l As Long, ByVal r As Long, a As Variant)

  'Description:	Two-sided recursion verison (may blow stack limit!)

Dim p As Long     'index of index of partition element
Dim v As Variant  'partition value
Dim rr As Long    'original right boundary
Dim temp As Long  'for swapping indices
    
rr = r                            'remember the index of end of partition
p = Int(l + r) / 2                'pick partition element
v = d(a(p))                       'remember partitioning value
temp = a(p): a(p) = a(l): a(l) = temp 'swap index of partition element to bottom
p = l                             'remember where we put it
Do
  Do
    For l = l + 1 To r - 1
      If v <= d(a(l)) Then Exit For
    Next
    For r = r To l + 1 Step -1
      If d(a(r)) <= v Then Exit Do
    Next
    If rr < l Then  'if we tacked an infinity on the right we could do without this
      l = l - 1
    ElseIf v < d(a(l)) Then
      l = l - 1
    End If
    temp = a(p): a(p) = a(l): a(l) = temp
    If p <= l - minSubFile Then Quicksort d, p, l - 1, a
    If rr - l >= minSubFile Then Quicksort d, l + 1, rr, a
    Exit Sub
  Loop
  temp = a(l): a(l) = a(r): a(r) = temp
  r = r - 1
Loop
End Sub

Private Sub InsertionSort(d As Variant, ByVal l As Long, ByVal r As Long, a As Variant)
    Dim j As Long       'index when searching for inversions
    Dim i As Long       'index when correcting inversions
    Dim v As Variant    'value of element on the move
    Dim oldI As Long    'index of element being moved
    For j = l + 1 To r
        If d(a(j)) < d(a(j - 1)) Then
            i = j - 1
            oldI = a(j)
            v = d(oldI)
            Do
                a(i + 1) = a(i)
                i = i - 1
                If i < l Then Exit Do   'because we don't have a magic low value at the bottom
            Loop Until d(a(i)) <= v
            a(i + 1) = oldI
        End If
    Next
End Sub

  sub randomizeArrayOrder(byref a As Variant)
    dim i    As Long  'index of element to swap
    dim j    As Long  'index of element to swap it with
    dim m    As Long  'number of elements remaining that might be swapped
    dim temp As Long  'for swapping elements (declare as Variant if a isn't an array of longs)
    
    m = ubound(a) - lbound(a) + 1
    for i=lbound(a) to ubound(a)-1
      j = i + int(m*rnd())
      temp=a(i)
      a(i)=a(j)
      a(j)=temp
      m = m - 1
    next    
  end sub

  sub randomizeArrayOrder(byref aIn as Variant, byref aOut as Variant) 'aIn gets trashed
    dim i As Long     'index of element to replace
    dim j AS Long     'index of element to replace it with
    dim m    As Long  'number of elements remaining that might be swapped
    m = ubound(aIn) - lbound(aIn) + 1
    for i=lbound(aIn) to ubound(aIn)
      j = i + int(m*rnd())
      aOut(i)=aIn(j)
      aIn(j)=aIn(i)		'<--two copies, not three, but aIn is trashed
      m=m-1
    next
  end sub

  'Note that for this to work aIn must not be the same array as aOut!  
  'You can trick VB into getting that right by calling it like this: randomizeArrayOrder (a),a
  'The brackets force VB to make a temporary copy of A.  VB can copy an entire array
  'in one hit, saving a lot of time.

  while l<r and element(l)<=pivot
    l = l + 1
  wend
  while l<r and pivot<=element(r)           'technically we don't need to check l<r here
    r = r - 1                               'in most implementations, as the pivot element itself
  wend                                      'will have been swapped with the leftmost element in the
                                            'current partition.
  'if we get here we may have two records to swap                                           

  sub rekeySortOrder(byref d as Variant, byref a as Variant)
    ' a=array of indices resulting from unstable sort pass
    '   containing one copy of each of the values between lbound(a)
    '   and ubound(a).
    ' d=the column input array supplied to the sorting pass.
    '
    dim class       'indicates the equivalence class for each element
    dim classCount  'count of elements in each equivalence class
    dim thisClass   'the current equivalence class number
    dim i           'index into array of indices
    dim j           'index into array of equivalence counts
    
    redim class(lbound(a) to ubound(a))
    redim classCount(0 to ubound(a)-lbound(a))
    
    class(lbound(a)) = 0
    classCount(0) = 1
    thisClass = 0
    for i=lbound(a)+1 to ubound(a)
      if d(a(i-1)) < d(a(i)) then
        thisClass = thisClass + 1
      end if
      class(a(i)) = thisClass
      classCount(thisClass) = classCount(thisClass) + 1
    next
    
    '
    'We now know which equivalence class each element belongs in
    'And we know how many elements there were in each class.
    'So now we just need to know where the rows in each equivalence
    'class will start in the output array, and we can re-order.
    '
    dim offset
    dim classOffset    
    redim classOffset(0 to thisClass)
    
    offset = lbound(a)
    for j=0 to thisClass
      classOffset(j) = offset
      offset = offset + classCount(j)
    next
    
    '
    'At this point we know (a) where each equivalence class should
    'start in the output array (classOffset gives us this),
    'and (b) which class each element is in (class gives us this).
    'So, we can interate through the class array, writing into the
    'part of the output order array for the appropriate equivalence
    'class.
    '
    for i = lbound(a) to ubound(a)
      j = class(i)
      a(classoffset(j)) = i
      classOffset(j) = classOffset(j) + 1
    next    
  end sub

Const minSubFile As Long= 9
Private partition2 As Variant 'Indices from the right that belong on the left   (in order found)
Private partition3 As Variant 'Indices from the left that belong in the centre  (ascending)

Sub ISortAlgorithm_Sort(d As Variant, ByVal l As Long, ByVal r As Long, a As Variant)
    If Not initialized Then
        ReDim partition2(l To r) As Long
        ReDim partition3(l To r) As Long
        initialized = True
    End If
    CheckBoundsEnough 0, r - l, partition2
    CheckBoundsEnough 0, r - l, partition3
    If l + minSubFile <= r Then Quicksort d, l, r, a 'Don't bother quicksorting small arrays
    InsertionSort d, l, r, a                'Straight insertion sort
End Sub

Private Sub InsertionSort(d As Variant, ByVal l As Long, ByVal r As Long, a As Variant)
    Dim j As Long       'index when searching for inversions
    Dim i As Long       'index when correcting inversions
    Dim v As Variant    'value of element on the move
    Dim oldI As Long    'index of element being moved
    For j = l + 1 To r
        If d(a(j)) < d(a(j - 1)) Then
            i = j - 1
            oldI = a(j)
            v = d(oldI)
            Do
                a(i + 1) = a(i)
                i = i - 1
                If i < l Then Exit Do   'because we don't have a magic low value at the bottom
            Loop Until d(a(i)) <= v
            a(i + 1) = oldI
        End If
    Next
End Sub

Private Sub Quicksort(d As Variant, ByVal l As Long, ByVal r As Long, a As Variant)
  Dim n As Long         'Number of elements in area to sort
  Dim m As Long         'Used to calculate stack size required
  Dim p As Long         'index of index of partition element
  Dim v As Variant      'partition value
  Dim i As Long         'Read point in initial index array
  Dim x As Long         'Original index at that point
  Dim vv As Variant     'Value of array element being compared to partition
  Dim p2Count As Long   'Number of elements from left belong in centre
  Dim p3Count As Long   'Number of elements from left belonging at right
  Dim w As Long         'Write point (when writing back new order for the partition)
  Dim lStack() As Long  'Stack of left boundaries
  Dim rStack() As Long  'Stack of right boundaries
  Dim stackSize As Long 'Initially, needed stack size.  Later, current stack size.
  Dim cl As Long        'Destination of first element matching current partition value
  Dim cr As Long        'Destination of last element matching current partition value
                        'cl and cr are used to determine the next partitions to sort.
  n = r - l + 1
  m = 1
  While m < n * 2
    m = m * 2
    stackSize = stackSize + 1
  Wend
  ReDim lStack(stackSize)   'We use a stack so we can avoid the overhead for
  ReDim rStack(stackSize)   'recursive procedure calls.  We only need to stack l and r!
  stackSize = -1            'Nothing in the stack
  
  Do
    p = Int(l + r) / 2                'pick partition element (could be random, doesn't really matter)
    v = d(a(p))                       'remember partitioning value
    p2Count = 0
    p3Count = 0
    w = l
    For i = l To r
        x = a(i)
        vv = d(x)
        If vv < v Then
          a(w) = x
          w = w + 1
        ElseIf v < vv Then
          partition3(p3Count) = x
          p3Count = p3Count + 1
        Else
          partition2(p2Count) = x
          p2Count = p2Count + 1
        End If
    Next
    cl = w   'first point in centre
    For i = 0 To p2Count - 1
        a(w) = partition2(i)
        w = w + 1
    Next
    cr = w - 1 'last point in centre
    For i = 0 To p3Count - 1
        a(w) = partition3(i)
        w = w + 1
    Next
    If cl - l < r - cr Then
        If minSubFile < r - cr Then
            stackSize = stackSize + 1
            lStack(stackSize) = cr + 1
            rStack(stackSize) = r
        End If
        r = cl - 1
    Else
        If minSubFile < cl - l Then
            stackSize = stackSize + 1
            lStack(stackSize) = l
            rStack(stackSize) = cl - 1
        End If
        l = cr + 1
    End If
    If r - l < minSubFile Then
      If stackSize = -1 Then Exit Do
      l = lStack(stackSize)
      r = rStack(stackSize)
      stackSize = stackSize - 1
    End If
  Loop
End Sub