算法实现/排序/Smoothsort
此 Smoothsort 的实现与 Dijkstra 的原始实现（在演示方面）有很大不同，经过了一些重大的重构。
为了在算法的固有复杂性下尽可能保持代码整洁，辅助函数被隔离在一个匿名命名空间中。
在实践中，下面的三个部分将按给定的顺序在单独的源文件中连接起来（例如，smoothsort.cpp）。
最后但并非最不重要的一点，代码使用了long long unsigned类型，这是 gcc 的扩展（不，不是）。在安装了gcc不可用的地方，用unsigned long替换（它只影响可排序数组的最大大小）。
/**
** 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;
}
}
#include <stddef.h>
#include <string.h>
/* Begin user-defined types and functions */
/* Type to sort (null-terminated string literals in this case) */
typedef const char *value_t;
/*
* Function that returns qsort-like comparison for parameters. A negative value
* would indicate that a goes before b, a positive value would indicate that a
* goes after b, and zero indicates that the elements are equivalent in order.
* An equivalent function for arithmetic types (e.g. int) would look like this:
static int cmp(value_t a, value_t b) {
return a - b;
}
*/
static int cmp(value_t a, value_t b) {
return strcmp(a, b);
}
/* End user-defined types and functions */
struct state {
value_t *a;
size_t q, r, p, b, c, r1, b1, c1;
};
#if __STDC_VERSION__ < 199901L
# define inline /* Silently prune inline qualifiers away for ANSI C */
#endif
static inline void up(size_t *a, size_t *b) {
size_t tmp;
tmp = *a;
*a += *b + 1;
*b = tmp;
}
static inline void down(size_t *a, size_t *b) {
size_t tmp;
tmp = *b;
*b = *a - *b - 1;
*a = tmp;
}
static void sift(struct state *s) {
size_t r0, r2;
value_t tmp;
r0 = s->r1;
tmp = s->a[r0];
while (s->b1 > 2) {
r2 = s->r1 - s->b1 + s->c1;
if (cmp(s->a[s->r1 - 1], s->a[r2]) >= 0) {
r2 = s->r1 - 1;
down(&s->b1, &s->c1);
}
if (cmp(s->a[r2], tmp) < 0) {
s->b1 = 1;
} else {
s->a[s->r1] = s->a[r2];
s->r1 = r2;
down(&s->b1, &s->c1);
}
}
if (s->r1 - r0 > 0) {
s->a[s->r1] = tmp;
}
}
static void trinkle(struct state *s) {
size_t p1, r0, r2, r3;
value_t tmp;
p1 = s->p;
s->b1 = s->b;
s->c1 = s->c;
r0 = s->r1;
tmp = s->a[r0];
while (p1 > 0) {
while ((p1 & 1) == 0) {
p1 >>= 1;
up(&s->b1, &s->c1);
}
r3 = s->r1 - s->b1;
if (p1 == 1 || cmp(s->a[r3], tmp) < 0) {
p1 = 0;
} else {
p1--;
if (s->b1 == 1) {
s->a[s->r1] = s->a[r3];
s->r1 = r3;
} else if (s->b1 >= 3) {
r2 = s->r1 - s->b1 + s->c1;
if (cmp(s->a[s->r1 - 1], s->a[r2]) >= 0) {
r2 = s->r1 - 1;
down(&s->b1, &s->c1);
p1 <<= 1;
}
if (cmp(s->a[r2], s->a[r3]) < 0) {
s->a[s->r1] = s->a[r3];
s->r1 = r3;
} else {
s->a[s->r1] = s->a[r2];
s->r1 = r2;
down(&s->b1, &s->c1);
p1 = 0;
}
}
}
}
if (s->r1 - r0 != 0) {
s->a[s->r1] = tmp;
}
sift(s);
}
static void semitrinkle(struct state *s) {
value_t tmp;
s->r1 = s->r - s->c;
if (cmp(s->a[s->r1], s->a[s->r]) >= 0) {
tmp = s->a[s->r];
s->a[s->r] = s->a[s->r1];
s->a[s->r1] = tmp;
trinkle(s);
}
}
void smoothsort(value_t *a, size_t n) {
struct state s;
size_t tmp;
s.a = a;
s.r = 0;
s.p = s.b = s.c = 1;
/* Build tree */
for (s.q = 1; s.q < n; s.q++) {
s.r1 = s.r;
if ((s.p & 7) == 3) {
s.b1 = s.b;
s.c1 = s.c;
sift(&s);
s.p = (s.p + 1) >> 2;
/* Two "up"s, saves us a little time */
tmp = s.b + s.c + 1;
s.b += tmp + 1;
s.c = tmp;
} else if ((s.p & 3) == 1) {
if (s.q + s.c < n) {
s.b1 = s.b;
s.c1 = s.c;
sift(&s);
} else {
trinkle(&s);
}
do {
down(&s.b, &s.c);
s.p <<= 1;
} while (s.b > 1);
s.p++;
}
s.r++;
}
s.r1 = s.r;
trinkle(&s);
/* Build sorted array */
while (s.q-- > 1) {
if (s.b == 1) {
s.r--;
s.p--;
while((s.p & 1) == 0) {
s.p >>= 1;
up(&s.b, &s.c);
}
} else if (s.b > 2) {
s.p--;
s.r = s.r - s.b + s.c;
if (s.p > 0) {
semitrinkle(&s);
}
down(&s.b, &s.c);
s.p = (s.p << 1) + 1;
s.r += s.c;
semitrinkle(&s);
down(&s.b, &s.c);
s.p = (s.p << 1) + 1;
}
/* element q processed */
}
/* element 0 processed */
}
/*
* Example main:
#include <stdio.h>
int main(int argc, const char **argv) {
const char *a[] = {
"the",
"quick",
"brown",
"fox",
"jumps",
"over",
"the",
"lazy",
"dog"
};
size_t i;
smoothsort(a, sizeof(a) / sizeof(a[0]));
for (i = 0; i < sizeof(a) / sizeof(a[0]); i++) {
puts(a[i]);
}
return 0;
}
*/
unit USmoothsort;
{ Delphi implementation of Dijkstra's algorithm }
interface
type TItem = Double; { data type }
function IsAscending(v1,v2: TItem): boolean; { comparison function }
{ sorting function }
procedure SmoothSort(var A: array of TItem);
implementation
{ customizable comparison function }
function IsAscending(v1,v2: TItem): boolean;
begin
result := v1<=v2;
end;
{ implementation of Djikstra's algorithm }
procedure SmoothSort(var A: array of TItem);
var q,r,
p,b,c,
r1,b1,c1,
N: integer;
procedure up(var vb,vc: integer);
var temp: integer;
begin
temp := vb;
vb := vb+vc+1;
vc := temp;
end;
procedure down(var vb,vc: integer);
var temp: integer;
begin
temp := vc;
vc := vb-vc-1;
vb := temp;
end;
procedure sift;
var r0, r2: integer;
T: TItem;
begin
r0 := r1;
T := A[r0];
while b1>=3 do
begin
r2 := r1-b1+c1;
if not IsAscending(A[r1-1],A[r2]) then
begin
r2 := r1-1;
down(b1,c1);
end;
if IsAscending(A[r2],T) then b1 := 1 else
begin
A[r1] := A[r2];
r1 := r2;
down(b1,c1);
end;
end;
if r1<>r0 then A[r1] := T;
end;
procedure trinkle;
var p1,r2,r3, r0 : integer;
T: TItem;
begin
p1 := p; b1 := b; c1 := c;
r0 := r1;
T := A[r0];
while p1>0 do
begin
while (p1 and 1)=0 do
begin
p1 := p1 shr 1; up(b1,c1);
end;
r3 := r1-b1;
if (p1=1) or IsAscending(A[r3],T) then p1 := 0 else
{p1>1}
begin
p1 := p1 - 1;
if b1=1 then
begin
A[r1] := A[r3];
r1 := r3;
end else
if b1>=3 then
begin
r2 := r1-b1+c1;
if not IsAscending(A[r1-1],A[r2]) then
begin
r2 := r1-1;
down(b1,c1); p1 := p1 shl 1;
end;
if IsAscending(A[r2],A[r3]) then
begin
A[r1] := A[r3];
r1 := r3;
end else
begin
A[r1] := A[r2];
r1 := r2;
down(b1,c1); p1 := 0;
end;
end;{if b1>=3}
end;{if p1>1}
end;{while}
if r0<>r1 then A[r1] := T;
sift;
end;
procedure semitrinkle;
var T: TItem;
begin
r1 := r-c;
if not IsAscending(A[r1],A[r]) then
begin
T := A[r];
A[r] := A[r1];
A[r1] := T;
trinkle;
end;
end;
begin
N := length(A);
q := 1; r := 0;
p := 1; b := 1; c := 1;
//building tree
while q < N do
begin
r1 := r;
if (p and 7) = 3 then
begin
b1 := b; c1 := c; sift;
p := (p+1) shr 2; up(b,c); up(b,c);
end
else if (p and 3) = 1 then
begin
if q + c < N then
begin
b1 := b; c1 := c; sift;
end else trinkle;
down(b,c); p := p shl 1;
while b > 1 do
begin
down(b,c); p := p shl 1;
end;
p := p+1;
end;
q := q + 1; r := r + 1;
end;
r1 := r; trinkle;
//bulding sorted array
while q>1 do
begin
q := q - 1;
if b = 1 then
begin
r := r - 1;
p := p - 1;
while (p and 1) = 0 do
begin
p := p shr 1; up(b,c);
end;
end else
if b >= 3 then
begin
p := p - 1; r := r-b+c;
if p > 0 then semitrinkle;
down(b,c); p := p shl 1 + 1;
r := r+c; semitrinkle;
down(b,c); p := p shl 1 + 1;
end;
//element q is done
end;
//element 0 is done
end;
end.
// by keeping these constants, we can avoid the tiresome business
// of keeping track of Dijkstra's b and c. Instead of keeping
// b and c, I will keep an index into this array.
static final int LP[] = { 1, 1, 3, 5, 9, 15, 25, 41, 67, 109,
177, 287, 465, 753, 1219, 1973, 3193, 5167, 8361, 13529, 21891,
35421, 57313, 92735, 150049, 242785, 392835, 635621, 1028457,
1664079, 2692537, 4356617, 7049155, 11405773, 18454929, 29860703,
48315633, 78176337, 126491971, 204668309, 331160281, 535828591,
866988873 // the next number is > 31 bits.
};
public static <C extends Comparable<? super C>> void sort(C[] m,
int lo, int hi) {
int head = lo; // the offset of the first element of the prefix into m
// These variables need a little explaining. If our string of heaps
// is of length 38, then the heaps will be of size 25+9+3+1, which are
// Leonardo numbers 6, 4, 2, 1.
// Turning this into a binary number, we get b01010110 = 0x56. We represent
// this number as a pair of numbers by right-shifting all the zeros and
// storing the mantissa and exponent as "p" and "pshift".
// This is handy, because the exponent is the index into L[] giving the
// size of the rightmost heap, and because we can instantly find out if
// the rightmost two heaps are consecutive Leonardo numbers by checking
// (p&3)==3
int p = 1; // the bitmap of the current standard concatenation >> pshift
int pshift = 1;
while (head < hi) {
if ((p & 3) == 3) {
// Add 1 by merging the first two blocks into a larger one.
// The next Leonardo number is one bigger.
sift(m, pshift, head);
p >>>= 2;
pshift += 2;
} else {
// adding a new block of length 1
if (LP[pshift - 1] >= hi - head) {
// this block is its final size.
trinkle(m, p, pshift, head, false);
} else {
// this block will get merged. Just make it trusty.
sift(m, pshift, head);
}
if (pshift == 1) {
// LP[1] is being used, so we add use LP[0]
p <<= 1;
pshift--;
} else {
// shift out to position 1, add LP[1]
p <<= (pshift - 1);
pshift = 1;
}
}
p |= 1;
head++;
}
trinkle(m, p, pshift, head, false);
while (pshift != 1 || p != 1) {
if (pshift <= 1) {
// block of length 1. No fiddling needed
int trail = Integer.numberOfTrailingZeros(p & ~1);
p >>>= trail;
pshift += trail;
} else {
p <<= 2;
p ^= 7;
pshift -= 2;
// This block gets broken into three bits. The rightmost
// bit is a block of length 1. The left hand part is split into
// two, a block of length LP[pshift+1] and one of LP[pshift].
// Both these two are appropriately heapified, but the root
// nodes are not necessarily in order. We therefore semitrinkle
// both of them
trinkle(m, p >>> 1, pshift + 1, head - LP[pshift] - 1, true);
trinkle(m, p, pshift, head - 1, true);
}
head--;
}
}
private static <C extends Comparable<? super C>> void sift(C[] m, int pshift,
int head) {
// we do not use Floyd's improvements to the heapsort sift, because we
// are not doing what heapsort does - always moving nodes from near
// the bottom of the tree to the root.
C val = m[head];
while (pshift > 1) {
int rt = head - 1;
int lf = head - 1 - LP[pshift - 2];
if (val.compareTo(m[lf]) >= 0 && val.compareTo(m[rt]) >= 0)
break;
if (m[lf].compareTo(m[rt]) >= 0) {
m[head] = m[lf];
head = lf;
pshift -= 1;
} else {
m[head] = m[rt];
head = rt;
pshift -= 2;
}
}
m[head] = val;
}
private static <C extends Comparable<? super C>> void trinkle(C[] m, int p,
int pshift, int head, boolean isTrusty) {
C val = m[head];
while (p != 1) {
int stepson = head - LP[pshift];
if (m[stepson].compareTo(val) <= 0)
break; // current node is greater than head. Sift.
// no need to check this if we know the current node is trusty,
// because we just checked the head (which is val, in the first
// iteration)
if (!isTrusty && pshift > 1) {
int rt = head - 1;
int lf = head - 1 - LP[pshift - 2];
if (m[rt].compareTo(m[stepson]) >= 0
|| m[lf].compareTo(m[stepson]) >= 0)
break;
}
m[head] = m[stepson];
head = stepson;
int trail = Integer.numberOfTrailingZeros(p & ~1);
p >>>= trail;
pshift += trail;
isTrusty = false;
}
if (!isTrusty) {
m[head] = val;
sift(m, pshift, head);
}
}