std::sort introsort implementation error

• Key: STDCXX-397
• Type: Bug
• Status: Resolved
• Resolution: Fixed
• Priority: Major
• Assignee: Martin Sebor
• Reporter: Joshua Lehrer
• Component/s: 25. Algorithms
• Affects Version/s: 4.1.3
• Fix Version/s: 4.2.0

Description

introsort is designed to detect when an input set would push quicksort into its worst-case scenario N^2 and fall back on a slower, yet still NLogN algorithm.

The implementation in __introsort_loop has a bug, however, and it fails to catch all of the scenarios. While I can not supply an exact input set to demonstrate the bug, I can explain the bug very easily.

First, allow me to paste in the code:

// David R. Musser's Introspective Sorting algorithm
// O(N * log (N)) worst case complexity
_EXPORT
template <class _RandomAccessIter, class _Dist, class _Compare>
void __introsort_loop (_RandomAccessIter __first, _RandomAccessIter __last,
_Dist __max_depth, _Compare __comp)
{
for (; __last - __first > __rw_threshold; __max_depth /= 2) {
if (0 == __max_depth) { __partial_sort (__first, __last, __last, _RWSTD_VALUE_TYPE (_RandomAccessIter), __comp); break; }

_RandomAccessIter __cut =
_unguarded_partition (_first, __last,
_median (*_first,
*(_first + (_last - __first) /2),
*(__last - 1), __comp), __comp);

// limit the depth of the recursion tree to log2 (last - first)
// where first and last are the initial values passed in from sort()
_introsort_loop (_cut, __last, __max_depth, __comp);
__last = __cut;
}
}

the variable '__max_depth' is supposed to be cut in half on each subsequent "recursive" call. Once it reaches zero, LogN recurisve calls have been made, and the algorithm falls back on a different sorting algorithm for the remainder.

The algorithm, as implemented, uses real recursion and tail recursion.

First, the pivot is selected, the pivot is done, and the algorithm has a left and a right half, hopefully balanced.

Consider what happens for the LEFT half, which is done using tail recursion. '_last' gets assigned 'cut', then the code goes to the top of the 'for' loop. The test condition of the loop is run, which divides '_max_depth' by two, bringing it closer to zero.

Now consider what happens for the RIGHT half, which is done using real recursion. The function is called recurisvely on the right. __max_depth is NOT cut in half.

What would happen if a poor pivot was selected causing the right half to be large and the left half to be small? What if that happens again and again? The real-recursion case is failing to decrement __max_depth until it starts working on the left half. You can see how if the algorithm continually built right-halves that were relatively large that __max_depth never gets decremented, and the algorithm never detects that it has made LogN recurisve calls.

I believe the proper fix is as follows:


// David R. Musser's Introspective Sorting algorithm
// O(N * log (N)) worst case complexity
_EXPORT
template <class _RandomAccessIter, class _Dist, class _Compare>
void __introsort_loop (_RandomAccessIter __first, _RandomAccessIter __last,
_Dist __max_depth, _Compare __comp)
{
for (; __last - __first > __rw_threshold; ) {
if (0 == __max_depth) { __partial_sort (__first, __last, __last, _RWSTD_VALUE_TYPE (_RandomAccessIter), __comp); break; } }

_RandomAccessIter __cut =
_unguarded_partition (_first, __last,
_median (*_first,
*(_first + (_last - __first) /2),
*(__last - 1), __comp), __comp);

// limit the depth of the recursion tree to log2 (last - first)
// where first and last are the initial values passed in from sort()
__max_depth /= 2;
_introsort_loop (_cut, __last, __max_depth, __comp);
__last = __cut;
}
}

"__max_depth/=2" is removed from the "for" loop and placed just above the two recursive calls.

This fixes the worst-case sample set that I have generated.

I look forward to your response,

joshua lehrer
http://www.lehrerfamily.com/

Comments

[Joshua Lehrer]

By the way, a killer set can be generated by:

1] change the selection of __cut to be 'first' instead of median of three
2] pass in an already sorted array

The algorithm will explode the stack.

Use my corrected code (still with the selection of __cut being "first") and sort the same already sorted array, and the algorithm will complete in NLogN time.

-josh

[Martin Sebor]

Thanks for the detailed report! We'll need some time to look into it. In the meantime, if you could come up with a test case that would be very very helpful. A user-defined type whose operator<() that keeps track of the number of comparisons should reliably reveal the problem given a suitable input sequence.

Martin,

It is very difficult to come up with a test case given the complexity of the issue. I have been trying and have yet to be able to. I will keep looking.

However, if you change the selection of the __cut to be 'first' instead of the median of three then sort an already sorted array, you will see the number of comparisons explode and O(N) recursive calls.

Looking at the code, you can clearly see that the recursion on the right half never has __max_depth decreased by half, the /2 is only ever applied to the left half. Therefore, the algorithm will detect too many recursive calls to the left, but not to the right.

I will report back if I can come up with an example using the median-of-three pivot selector.

-josh

[Todd Gardner]

Martin,

This code can generate a killer set for the problem:

#include <vector>
#include <fstream>

const int SetSize = 1000;
int main(){
using std::vector;
vector<int> KillerSet;

KillerSet.push_back(9998);
KillerSet.push_back(9997);
KillerSet.push_back(9997);

int CurrNumber = 9996;
while(KillerSet.size() < SetSize) { vector<int>::iterator Middle = KillerSet.end()-((KillerSet.size()/2)); int MidValue = *Middle; *Middle = CurrNumber; KillerSet.push_back(MidValue); KillerSet.push_back(CurrNumber-1); CurrNumber-=2; }
ofstream OutFile("killerset2.txt");
for(vector<int>::reverse_iterator KSCurr = KillerSet.rbegin(), KSEnd = KillerSet.rend(); KSCurr !=KSEnd; ++KSCurr) { OutFile << *KSCurr << "\n"; }
OutFile.close();
}

Modifying the algorithm to keep track of depth on the right side recursion:

// David R. Musser's Introspective Sorting algorithm
// O(N * log (N)) worst case complexity
_EXPORT
template <class _RandomAccessIter, class _Dist, class _Compare>
void __introsort_loop (_RandomAccessIter __first, _RandomAccessIter __last,
_Dist __max_depth, _Compare __comp, int Depth) // Pass in 0 on the initial call
{
std::cout << "Depth: " << Depth << endl;
for (; __last - __first > __rw_threshold; __max_depth /= 2) {
if (0 == __max_depth) { __partial_sort (__first, __last, __last, _RWSTD_VALUE_TYPE (_RandomAccessIter), __comp); break; }

_RandomAccessIter __cut =
_unguarded_partition (_first, __last,
_median (*_first,
*(_first + (_last - __first) /2),
*(__last - 1), __comp), __comp);

// limit the depth of the recursion tree to log2 (last - first)
// where first and last are the initial values passed in from sort()
_introsort_loop (_cut, __last, __max_depth, __comp, Depth+1);
__last = __cut;
}
}

And then running it on the output of the first code gives the output:
Depth: 0
Depth: 1
Depth: 2
Depth: 3
Depth: 4
...
Depth: 489
Depth: 490
Depth: 491
Depth: 492
Depth: 493

Showing that the recursion doesn't have a logarithmic limit on the right hand side.

[Joshua Lehrer]

Here is code which creates a killer set then sorts it, counting the compares.

Without our fix, sorting an array of 99,999 elements requires 2,500,299,945 compares. With our fix it takes only 3,521,701.
Without our fix, sorting an array of 999 elements requires 252,945 compares. With our fix it takes only 21,148.

//begin code
#include <stdlib>
#include <vector>
#include <iostream>
#include <algorithm>

//default comparison with operator< but counts compares
struct MyCompare {
MyCompare(unsigned int &compares) : m_compares(compares) {
}
template <typename T>
inline bool operator()(const T& lhs, const T&rhs) { ++m_compares; return lhs<rhs; }
unsigned int &m_compares;
};

//pass in size of array to create/test
int main(int argc, char **argv) {
using namespace std;

//demonstration works for odd size arrays only
const int mxtemp = atoi(argv[1]);
const int mx = (mxtemp%2) ? mxtemp : mxtemp+1;

cout << "array size=" << mx << endl;

//killer set
vector<int> results;
results.resize(mx);

{
cout << "generating kiler set" << endl;

vector<int> indicies;
indicies.reserve(mx);
for (int i=0; i<mx; ++i) { indicies.push_back(i); }

int i = 0;
while (i<mx) {
results[indicies[i]]=i;
int mid = (mx-i-1) /2 + i;
if (mid != i) { results[indicies[mid]]=i+1; using namespace std; swap(indicies[i+1],indicies[mid]); }
i+=2;
}
}

cout << "sorting" << endl;
unsigned int compares=0;
sort(results.begin(),results.end(),MyCompare(compares));

cout << "total compares=" << compares << endl;
}

[Martin Sebor]

Thanks for the test case!

I've reviewed the Introsort paper (www.cs.rpi.edu/~musser/gp/introsort.ps) and confirmed that your suggested fix is appropriate since it follows the algorithm published in the paper. I'm not sure about the origin of our implementation but it clearly differs in a subtle but important aspect from the reference algorithm. I'll commit a corrected version shortly. Unless you would prefer otherwise I will also add your test case to our regression test suite.

[Joshua Lehrer]

You are welcome to add my test case to your regression test suite. Thanks for your help and confirmation!

-josh

[Martin Sebor]

Resolved. Will close after the test case has been added to the regression test suite.