State of affairs
Having a computer sience background and an interest in low-level systems programming, I have come to survey some sortings algorithms. Starting point was Rich Felker's post on musl (here: http://www.etalabs.net/compare_libcs.html). I wondered: Is dietlibc's sorting algorithm really that bad?So I looked it up, and sure enough: It was written exactly as described. And hadn't changed that much in the years since I last looked. Although the code had gained two optional pivot-selection algorithms, it still has largely the same two problems as back then, though one of them is largely mitigated. Maybe I should condense it down first: The algorithm still has a worst case O(n²) time complexity and a worst case O(n) space complexity on the stack. Let me repeat that for people without a background in complexity theory: If an attacker can control length and content of an array to be sorted with qsort(), if an application doing that is linked against dietlibc, the attacker can force a stack smash pretty much all the time. And that code is coming from Fefe, a self-proclaimed security expert. Great.
I really don't know why Quicksort was chosen. Probably the author thought that the function's called qsort(), might as well use that algorithm. But on its face, that was a poor choice: Quicksort is not stable, not adaptive, and not iterative. And not easy to implement: The implementation clocks in at 28 lines of code without the pivot choosing and without the wrapper to make the interface standards compliant. That was the state of the algorithm when I looked at it a few years ago. Now, considering dietlibc's stated goal is to not waste a single byte, Quicksort really was a curious choice. But without the pivot choosing and the two recursions the algorithm has not only a really bad worst case, but those are also known. So an attacker can just hit the algorithm with a known bad case and crash the server.
The code makes reference to some talk slides by Robert Sedgewick (here: http://cs.princeton.edu/~rs/talks/QuicksortIsOptimal.pdf) and was pretty much copied from one of the slides there. Now stop me if this sounds strange to you, but the code was already filling the whole slide. Might it be possible, that Sedgewick only presented a simplified version of the code so as not to have to get out the bigger chalkboard?
Now obviously, that was sarcasm, but I try to let my critique be constructive. This is the code at time of writing:
static void exch(char* base,size_t size,size_t a,size_t b) {
char* x=base+a*size;
char* y=base+b*size;
while (size) {
char z=*x;
*x=*y;
*y=z;
--size; ++x; ++y;
}
}
#define RAND
/* Quicksort with 3-way partitioning, ala Sedgewick */
/* Blame him for the scary variable names */
/* http://www.cs.princeton.edu/~rs/talks/QuicksortIsOptimal.pdf */
static void quicksort(char* base,size_t size,ssize_t l,ssize_t r,
int (*compar)(const void*,const void*)) {
ssize_t i=l-1, j=r, p=l-1, q=r, k;
char* v=base+r*size;
if (r<=l) return;
#ifdef RAND
/*
We chose the rightmost element in the array to be sorted as pivot,
which is OK if the data is random, but which is horrible if the
data is already sorted. Try to improve by exchanging it with a
random other pivot.
*/
exch(base,size,l+(rand()%(r-l)),r);
#elif defined MID
/*
We chose the rightmost element in the array to be sorted as pivot,
which is OK if the data is random, but which is horrible if the
data is already sorted. Try to improve by chosing the middle
element instead.
*/
exch(base,size,l+(r-l)/2,r);
#endif
for (;;) {
while (++i != r && compar(base+i*size,v)<0) ;
while (compar(v,base+(--j)*size)<0) if (j == l) break;
if (i >= j) break;
exch(base,size,i,j);
if (compar(base+i*size,v)==0) exch(base,size,++p,i);
if (compar(v,base+j*size)==0) exch(base,size,j,--q);
}
exch(base,size,i,r); j = i-1; ++i;
for (k=l; k<p; k++, j--) exch(base,size,k,j);
for (k=r-1; k>q; k--, i++) exch(base,size,i,k);
quicksort(base,size,l,j,compar);
quicksort(base,size,i,r,compar);
}
void qsort(void* base,size_t nmemb,size_t size,int (*compar)(const void*,const void*)) {
/* check for integer overflows */
if (nmemb >= (((size_t)-1)>>1) ||
size >= (((size_t)-1)>>1)) return;
#if 0
if (sizeof(size_t) < sizeof(unsigned long long)) {
if ((unsigned long long)size * nmemb > (size_t)-1) return;
} else {
if (size*nmemb/nmemb != size) return;
}
#endif
if (nmemb>1)
quicksort(base,size,0,nmemb-1,compar);
}
Now, with only a bit of tweaking, the devastating O(n) space complexity on the cache can be avoided, no matter the input:
static void quicksort(char* base,size_t size,ssize_t l,ssize_t r,
int (*compar)(const void*,const void*)) {
ssize_t i=l-1, j=r, p=l-1, q=r, k;
char* v=base+r*size;
while (l < r) {
#ifdef RAND
/*
We chose the rightmost element in the array to be sorted as pivot,
which is OK if the data is random, but which is horrible if the
data is already sorted. Try to improve by exchanging it with a
random other pivot.
*/
exch(base,size,l+(rand()%(r-l)),r);
#elif defined MID
/*
We chose the rightmost element in the array to be sorted as pivot,
which is OK if the data is random, but which is horrible if the
data is already sorted. Try to improve by chosing the middle
element instead.
*/
exch(base,size,l+(r-l)/2,r);
#endif
for (;;) {
while (++i != r && compar(base+i*size,v)<0) ;
while (compar(v,base+(--j)*size)<0) if (j == l) break;
if (i >= j) break;
exch(base,size,i,j);
if (compar(base+i*size,v)==0) exch(base,size,++p,i);
if (compar(v,base+j*size)==0) exch(base,size,j,--q);
}
exch(base,size,i,r); j = i-1; ++i;
for (k=l; k<p; k++, j--) exch(base,size,k,j);
for (k=r-1; k>q; k--, i++) exch(base,size,i,k);
if (j - l < r - i) {
quicksort(base, size, l, j, compar);
l = i;
} else {
quicksort(base, size, i, r, compar);
r = j;
}
}
}
The difference is marked in bold. It has the effect that the code will now recurse into the smaller side of the array and tail call into the smaller side. That shifts the worst case for space complexity: The highest number of stack frames can only come from an exactly halfed array, so the number just got cut down to a logarithmic stack usage, which is unproblematic most of the time.
Now, that takes care of the space complexity. The default choice for random pivot selection also means that the code will very likely skirt around the worst case, but we don't need the worst case for pathological behaviour (if the code chooses the pivot so that there are only two elements on one side, that's also pretty bad). Since that code suddenly uses randomness in a sorting algorithm, all we can say about the time complexity is that it's somewhere between O(n log n) and O(n²). And it doesn't get much more precise. Since the code uses the rand()-Interface, it is possible for an attacker to guess the state of the random number generator and calculate a pathological case from that.
Sedgewick himself argued in the very same talk that was linked to, for median-of-three pivot selection. I implemented that today:
[pivot selection]
...
#elif defined MEDIAN_OF_3
/* use median of three algorithm, that is, of left, middle, and right, pick the one that's neither maximum nor minimum. */
ssize_t pivot;
ssize_t mid = (r + l) / 2;
if (compar(base + size * l, base + size * mid) < 0) {
if (compar(base + size * mid, base + size * r) < 0)
pivot = mid;
else if (compar(base + size * l, base + size * r) < 0)
pivot = r;
else
pivot = l;
} else {
if (compar(base + size * l, base + size * r) < 0)
pivot = l;
else if (compar(base + size * mid, base + size * r) < 0)
pivot = r;
else
pivot = mid;
}
if (pivot != r)
exch(base, size, pivot, r);
#endif
That's hardly readable, but correct as far as I can see. But of course there exist pathological cases for this one as well. Randomness is unique in that it has no consistent counter, I guess. The problem is, it adds two or three comparison calls, but will only ever move a single element, and that one to a place where we know it is wrong. Not the best use of that information, but what do you do?
However, I also saw that you can go overboard: In newlib, the sorting algorithm is bubble sort, if there are less than 7 elements to be sorted. I guess they measured that bubble sort is faster below that threshold. Also, the pivot choosing is insane: Not only is it median of three, but if there are more than 40 elements, it will choose the median of three among the medians of three from three cherry picked elements near the bottom, near the top and in the middle. So, hopefully the median of the 9, though you can never know that. And even then: why? Because they could!
musl uses the Smoothsort algorithm, which has the benefit of being adaptive, so it becomes faster, the more sorted the input already is. It's actually the better algorithm here, having a lower worst case complexity, and needing a constant use of stack space.
glibc goes overboard in the other direction: For large inputs, it tries to allocate the memory on the heap, and uses quicksort, if that doesn't work. It also for some reason tries to figure out if such a request can succeed beforehand. I don't know, I always figured that just calling malloc() would tell me if I had the memory or not (overcommit notwithstanding). And if they have the memory, they run the most complicated mergesort I have ever seen! What are they doing?
But still, why didn't the dietlibc people just go with shellsort? Their promise was to create the smallest possible C library, and here they are using an algorithm that's neither the fastest one known, nor the smallest one known. That would be shellsort. So I'm going to try to implement a shellsort here.
Implementing Shellsort
Shellsort is a repeated application of insertion sort, to create a number of interleaved sorted subsequences, where the interval gets shorter each iteration, finally ending in 1. There has been a lot of discussion about the sequence of intervals used, but I'm going to go with Pratt's sequence, because that is very cheap to calculate, and we have both the closed and the iterative form. In closed form the formula is:a(n) = (3^n - 1) / 2
And iterative:
a(0) = 1
a(n + 1) = 3a(n) + 1
The first one lends itself to an analysis of how much space is going to be necessary. Here, we have the problem
(3^n - 1) / 2 > 2^k
We want to figure out, for a given k, how high n must be at least for this to be true. So we see that we can multiply this inequation with 2
3^n - 1 > 2^(k + 1)
Since this is an inequation, I can just change either side in the direction of their side. So, for instance, I can just add 1 to the left side, since it already is bigger than the right side, and adding 1 will only make it bigger still. Then I can push both sides through the binary logarithm, because that way, we can lose a factor:
n lb 3 > k + 1
Now all that's left is dividing by lb 3. Since that's positive, it won't flip the inequality sign.
n > (k + 1) / lb 3
Great, so if I calculate the right side of this with k being the width of a machine word, and round up, I have the number of the first element of the sequence exceeding the maximum number that can be put into a machine word. Thanks to the transdichotomous machine model (I hope I wrote that correctly), that means we can never need more machine words than that to save all elements of Pratt's sequence that's ever going to be smaller than an in-memory object. So that comes out at 21 for for 32 bits and 42 for 64 bits. The easiest is going to be:
size_t gaps[sizeof(size_t) * 6];
This allocates 24 machine words for a 32 bit machine, and 48 for a 64 bit machine. Both should fit easily onto the stack. So, for the interval allocation and calculation, my code is this here:
size_t gaps[sizeof(size_t) * 6];
size_t *gap = gaps;
for (*gap++ = 1; gap[-1] < nmemb / 3; *gap = 3 * gap[-1] + 1, gap++);
Now, insertion sort needs a temporary memory space, which I'm allocating on the stack (A single element should fit onto the stack easily). Then the rest is just:
void qsort(void* base, size_t nmemb, size_t size, int (*cmp)(const void*, const void*)) {
size_t gaps[sizeof(size_t) * 6];
size_t *gap = gaps;
char *temp = alloca(size);
char *cbase = base;
for (*gap++ = 1; gap[-1] < nmemb / 3; *gap = 3 * gap[-1] + 1, gap++);
for (; gap >= gaps; gap--) {
size_t i, j;
for (i = *gap; i < nmemb; i++) {
memcpy(temp, cbase + i * size, size);
for (j = i; j >= *gap && cmp(cbase + (j - *gap) * size, temp) > 0; j -= *gap)
memcpy(cbase + j * size, cbase + (j - *gap) * size, size);
memcpy(cbase + j * size, temp, size);
}
}
}
That's getting pretty concise already. Note that that's the entire algorithm: No wrapper for interface compatibility, no procedure to exchange elements. But one thing is obvious: The variables i, j, and gap are never referred to by themselves, they are always multiplied with size (eventually). Thankfully, we can reduce the number of multiplications heavily, by scaling the gap sequence by size, incrementing i by size instead of 1, and then remove all the multiplications.
Scaling the sequence is actually pretty easy, thanks to the iterative form used for the calculation. If we multiply both formulas with a constant s, we get
sa(0) = s * 1 = s
sa(n + 1) = s(3 a(n) + 1) = 3sa(n) + s
Which means we can just rewrite the above function a bit:
void qsort(void* base, size_t nmemb, size_t size, int (*cmp)(const void*, const void*)) {
size_t high = nmemb * size;
size_t gaps[sizeof(size_t) * 6];
size_t *gap = gaps;
char *temp = alloca(size);
char *cbase = base;
for (*gap++ = size; gap[-1] < high / 3; *gap = 3 * gap[-1] + size, gap++);
for (; gap >= gaps; gap--) {
size_t i, j;
for (i = *gap; i < high; i += size) {
memcpy(temp, cbase + i, size);
for (j = i; j >= *gap && cmp(cbase + j - *gap, temp) > 0; j -= *gap)
memcpy(cbase + j, cbase + j - *gap, size);
memcpy(cbase + j, temp, size);
}
}
}
Changes again in bold. Note how the number of multiplications by size shrunk down to a single one. Oh, and that one cannot overflow: The object has to fit entirely into memory, which means that between base and base + nmemb * size, there must not be an inaccessible location. That would be violated if nmemb * size overflowed. I don't check for that, because qsort() has no way to return failure.
