Hi there, long time no see.
I had a question about smoothsort, I can't figure it out, I know it's like heapsort and I'm familiar with the heapsort implementation, but I just can't get my head around smoothsort.
For those unfamiliar with smoothsort, it's a variation of heapsort designed to run close to O(n) if the list is initially mostly sorted.
Here is a sample implementation I found at wikibooks, note, I do understand what each step of the code is doing (syntax, bitshifts, templates, etc), what I don't understand is the logic behind what's happening, I need to know the why to the code, not the what. And it's the what I'm having trouble with. :)
As for a "why are you using this?" I'm doing it to make my life a bit more interesting in a boring C++ class.
/**
** SmoothSort function template + helper functions.
**
** Formal type T should have a comparison operator >= with prototype:
**
** bool T::operator >= (const T &) const throw ();
**
** which should compare its arguments by the given relation
** (possibly taking advantage of the type itself).
**
**
**/
/** Sort an array in ascending order. **/
template <typename T>
void smoothsort (T *, unsigned) throw ();
namespace
/**
** Helper function's local namespace (declarations).
**
**/
{
class LeonardoNumber
/**
** Helper class for manipulation of Leonardo numbers
**
**/
{
public:
/** Default ctor. **/
LeonardoNumber (void) throw () : b (1), c (1)
{ return; }
/** Copy ctor. **/
LeonardoNumber (const LeonardoNumber & _l) throw () : b (_l.b), c (_l.c)
{ return; }
/**
** Return the "gap" between the actual Leonardo number and the
** preceding one.
**/
unsigned gap (void) const throw ()
{ return b - c; }
/** Perform an "up" operation on the actual number. **/
LeonardoNumber & operator ++ (void) throw ()
{ unsigned s = b; b = b + c + 1; c = s; return * this; }
/** Perform a "down" operation on the actual number. **/
LeonardoNumber & operator -- (void) throw ()
{ unsigned s = c; c = b - c - 1; b = s; return * this; }
/** Return "companion" value. **/
unsigned operator ~ (void) const throw ()
{ return c; }
/** Return "actual" value. **/
operator unsigned (void) const throw ()
{ return b; }
private:
unsigned b; /** Actual number. **/
unsigned c; /** Companion number. **/
};
/** Perform a "sift up" operation. **/
template <typename T>
inline void sift (T *, unsigned, LeonardoNumber) throw ();
/** Perform a "semi-trinkle" operation. **/
template <typename T>
inline void semitrinkle (T *, unsigned, unsigned long long, LeonardoNumber) throw ();
/** Perform a "trinkle" operation. **/
template <typename T>
inline void trinkle (T *, unsigned, unsigned long long, LeonardoNumber) throw ();
}
template <typename T>
void smoothsort (T * _m, unsigned _n) throw ()
/**
** Sorts the given array in ascending order.
**
** Usage: smoothsort (<array>, <size>)
**
** Where: <array> pointer to the first element of the array in question.
** <size> length of the array to be sorted.
**
**
**/
{
if (!(_m && _n)) return;
unsigned long long p = 1;
LeonardoNumber b;
for (unsigned q = 0; ++q < _n ; ++p)
if (p % 8 == 3)
{
sift<T> (_m, q - 1, b);
++++b; p >>= 2;
}
else if (p % 4 == 1)
{
if (q + ~b < _n) sift<T> (_m, q - 1, b);
else trinkle<T> (_m, q - 1, p, b);
for (p <<= 1; --b > 1; p <<= 1) ;
}
trinkle<T> (_m, _n - 1, p, b);
for (--p; _n-- > 1; --p)
if (b == 1)
for ( ; !(p % 2); p >>= 1) ++b;
else if (b >= 3)
{
if (p) semitrinkle<T> (_m, _n - b.gap (), p, b);
--b; p <<= 1; ++p;
semitrinkle<T> (_m, _n - 1, p, b);
--b; p <<= 1; ++p;
}
return;
}
namespace
/**
** Helper function's local namespace (definitions).
**
**/
{
template <typename T>
inline void sift (T * _m, unsigned _r, LeonardoNumber _b) throw ()
/**
** Sifts up the root of the stretch in question.
**
** Usage: sift (<array>, <root>, <number>)
**
** Where: <array> Pointer to the first element of the array in
** question.
** <root> Index of the root of the array in question.
** <number> Current Leonardo number.
**
**
**/
{
unsigned r2;
while (_b >= 3)
{
if (_m [_r - _b.gap ()] >= _m [_r - 1])
r2 = _r - _b.gap ();
else
{ r2 = _r - 1; --_b; }
if (_m [_r] >= _m [r2]) break;
else
{ swap(_m [_r], _m [r2]); _r = r2; --_b; }
}
return;
}
template <typename T>
inline void semitrinkle (T * _m, unsigned _r, unsigned long long _p,
LeonardoNumber _b) throw ()
/**
** Trinkles the roots of the stretches of a given array and root when the
** adjacent stretches are trusty.
**
** Usage: semitrinkle (<array>, <root>, <standard_concat>, <number>)
**
** Where: <array> Pointer to the first element of the array in
** question.
** <root> Index of the root of the array in question.
** <standard_concat> Standard concatenation's codification.
** <number> Current Leonardo number.
**
**
**/
{
if (_m [_r - ~_b] >= _m [_r])
{
swap(_m [_r], _m [_r - ~_b]);
trinkle<T> (_m, _r - ~_b, _p, _b);
}
return;
}
template <typename T>
inline void trinkle (T * _m, unsigned _r, unsigned long long _p,
LeonardoNumber _b) throw ()
/**
** Trinkles the roots of the stretches of a given array and root.
**
** Usage: trinkle (<array>, <root>, <standard_concat>, <number>)
**
** Where: <array> Pointer to the first element of the array in
** question.
** <root> Index of the root of the array in question.
** <standard_concat> Standard concatenation's codification.
** <number> Current Leonardo number.
**
**
**/
{
while (_p)
{
for ( ; !(_p % 2); _p >>= 1) ++_b;
if (!--_p || (_m [_r] >= _m [_r - _b])) break;
else
if (_b == 1)
{ swap(_m [_r], _m [_r - _b]); _r -= _b; }
else if (_b >= 3)
{
unsigned r2 = _r - _b.gap (), r3 = _r - _b;
if (_m [_r - 1] >= _m [r2])
{ r2 = _r - 1; _p <<= 1; --_b; }
if (_m [r3] >= _m [r2])
{ swap(_m [_r], _m [r3]); _r = r3; }
else
{ swap(_m [_r], _m [r2]); _r = r2; --_b; break; }
}
}
sift<T> (_m, _r, _b);
return;
}
}
Build a man a fire, warm him for a day,
Set a man on fire, warm him for the rest of his life.
Short answer: I've seen that paper, and I actually don't understand it at all. :/
Long answer: I've tried trudging through those; however, I can't seem to determine how the tree is set up, and how the algorithm actually sorts.
Heapsort uses a binary heap for it's tree and the largest values are at the base of the tree.
Smoothsort uses (not a binary?) tree based on "Leonardo Numbers", and I don't know how the values are arranged and pruned.
I don't understand the difference between "trusty" and "dubious" "stretches." It does explain that in the paper, and I've found the explanation but it still befuddles me.
I've tried working the algorithm out on paper to see if I can grok it by doing, which is always a good method, but I screwed something up because it didn't sort correctly, and I couldn't always figure out ahead of time where the logic would go.
EDIT:
The code you quoted is pretty obfuscated despite the airy and well-commented appearance.
Eh, I can read the code and I understand what each lines does, the code in the paper you linked is just as bad, probably worse.
There was also a Delphi implementation on the wikibooks site that was a bit easier on the eyes. But I don't know Delphi like I do C.
Build a man a fire, warm him for a day,
Set a man on fire, warm him for the rest of his life.
The code you quoted is pretty obfuscated despite the airy and well-commented appearance.
Eh, I can read the code and I understand what each lines does, the code in the paper you linked is just as bad, probably worse.
Well, sure, I can too, but that doesn't mean I understand its intent any more than you do. That's what I refer to by obfuscated.
For example, I can tell just fine what each statement in this code (by Peter Klausler, from IOCCC 2005) does:
typedef unsigned char B;
char *x[] = {
#include "dict.h"
0
};
typedef struct L
{
B *s;
struct L *n;
} L;
L *h[128], *l[128], *s[128], Z[sizeof x / sizeof *x], *F = Z;
int c[256], m, a = 1;
int k (B * q)
{
int g = 0;
B *p = q;
while (*p)
g |= !c[*p++]--;
return g - 1 & p - q;
}
void u (B * p)
{
while (*p)
c[*p++]++;
}
void S (int N, int r, int t, L * W)
{
L *w;
int i, n;
for (n = r < N ? r : N; n > 0; n--)
for (w = n == N ? W : h[n];
s[t] = w;
u (w->s), w = w->n)
if (k (w->s))
if (n == r)
{
if (t == m - 1)
for (i = a = 0; i <= t; i++)
printf ("%s%c", s[i]->s, i < t ? ' ' : '\n');
}
else if (t < m - 1)
S (n, r - n, t + 1, s[t] = w);
}
int main (int C, B ** A)
{
int i = 0, g, n = 0;
B *p;
while (--C)
for (p = *++A; n < 127 && *p;)
c[*p++]++, n++;
for (; p = x[i++]; u (p))
if (g = k (p))
(l[g] = *(l[g] ? &l[g]->n : &h[g]) = F++)->s = p;
while (++m < 128)
S (127, n, 0, h[127]);
return a;
}
However, catching the gist of its purpose, or the algorithms employed, is another thing. (For the record, it "finds best anagrams for word or phrase given as program arguments".)
Although I'd argue that understanding this IOCCC program is easier than understanding that shuffle sort algorithm.
Joined: 4/25/2004
Posts: 615
Location: The Netherlands
From hits like http://stochastix.wordpress.com/2007/11/10/smooth-sorting/ I conclude that it uses some hardcore math (explaining the obfuscated code, I guess there's no way out of that) and I literally dont have the time to get into it now (eg. trying to understand it).
The code you quoted is pretty obfuscated despite the airy and well-commented appearance.
<nitpicking> Not to talk about the C++ standard saying that symbol names starting with an underscore are reserved for the use of the compiler and should not be used by the user... </nitpicking>
Still compiles and works, has less obnoxious operator overloads and less _ characters.
And is shorter and therefore less demanding on human memory, and thus easier to understand.
Though its working still defies me.
I can't even find the word "trinkle" in a dictionary.
I did some quick test on the smoothsort. Instead of measuring run time, I am just counting comparisons. "Almost sorted" isn't really very meaningful so instead let us look at making small changes to a sorted array and how many more comparisons it requires. Some possible changes are swap two adjacent items, swap two distant items, and move an item to some other location, sliding items in between.
Comparison counts, sorting a list with smoothsort of size 100,000
199971 in order already
4289498 random shuffled
3253547 reversed
201517 random adjacent swaps (1000 of them), avg +1.5 per swap
309452 random distance swaps (1000 of them), avg +109 per swap
1411142 random slides (1000 of them), avg +1211 per slide
It seems to handle slides really poorly, even though I would consider those as nearly sorted, for example 0 2 3 4 5 6 7 1 8 9 is pretty close to sorted.... Also alot of those 1000 slides, were canceling each out to some degree so the situation is even worse, most elements were less than 10 positions from their sorted locations, with the 1% of points that were slid being about 40,000 away.
I created a new algorithm to see how it does. The algorithm is very simple, it basically repeatedly does a binary search to find where an element goes, then places it there. To make it work well for nearly sorted data, it starts the search where the previous element was placed. Source is here
template<typename RI, typename IntType>
void insert(IntType val, RI a, RI b) {
for(; a<=b; a++) swap(val, *a);
}
template<typename RI, typename IntType>
RI bsxfind(IntType v,RI a, RI b, RI x, int width, RI lb, RI rb) {
if(width >= (rb - lb) * 2) return lower_bound(lb,rb,v);
if(v<*x) return bsxfind(v,a,b,x-width,width*2, lb, x);
else return bsxfind(v,a,b,x+width,width*2, x+1, rb);
}
template<typename RandomIterator>
void mysort(RandomIterator begin, RandomIterator end) {
RandomIterator i,last;
for(i=begin; i<end; i++) {
last=bsxfind(*i,begin,i,last,1,begin,i);
insert(*i,last,i);
}
}
Since I know what the algorithm does I will list its asymptotic behavior. Note that I am only counting comparisons, my algorithm is actually n^2 for random data because the insert function takes linear time due to my data structure just being an array.
already sorted: O(n)
reversed: O(n)
randomly shuffled: O(n log(n))
k adjacent swaps: O(n+k)
k distance swaps: O(n+k*log(n))
k slides: O(n+k*log(n))
Now repeating same measurements.
99999 in order already
2648698 random shuffled
99999 3253547 reversed
101991 random adjacent swaps (1000 of them)
275229 random distance swaps (1000 of them)
236662 random slides (1000 of them)
So it beats smoothsort in every category. I did not list change avg increase per swap because for already sorted data it needs only 1 comparison per item. But if the first element is actually the largest it requires 2 for per item. This has no effect on asymptotic behavior, but would make the results less meaningful. So here are results again, but first and last element swapped before all other changes.
200028 in order already
200027 reversed
201064 random adjacent swaps (1000 of them), avg +1 per swap
275513 random distance swaps (1000 of them), avg +75 per swap
236956 random slides (1000 of them), avg +37 per slide
So it is not hard to make an algorithm that beats smoothsort in comparisons. I would guess that smoothsort is similar, but overcomes the linear insertion time using fibonacci heaps (I don't know what those are yet) but at the cost of more comparisons (but not asymptotically more except slides/reverse). However I could be completely wrong and it uses some other method to seed nearly sorted things well, that could explain its bad slide handling.
Omni, that was a fun experiment for me but I also hope it was of some use for you in understanding how smoothsort works, although maybe you were already beyond that part and now focussed on how fibonacci heaps work.
template<typename RI, typename IntType>
void insert(IntType val, RI a, RI b) {
for(; a<=b; a++) swap(val, *a);
}
That's a really obfuscated way of performing an insertion. It took me something like 5 minutes to understand what it's doing.
Anyways, the most efficient way of performing an insertion into an array is not to swap each element with something else. By performing swaps you are performing 2 useless assignments per iteration. This is a more efficient way of doing it:
template<typename RI, typename IntType>
void insert(IntType val, RI where, RI last)
{
for(; last != where; --last) *last = *(last-1);
*where = val;
}
It only performs (last-where+1) assignments rather than 3*(last-where+1).
(Btw, with the next C++ standard the elements should be *moved* to their places rather than assigned. That makes it more efficient for elements which are aware of move semantics. Your swapping solution already takes automatically care of the move semantics in the next standard, but it nevertheless performs too many of them.)
