# Algorithm Implementation/Sorting/Smoothsort

## C++

This implementation of Smoothsort is substantially different (in presentation) from Dijkstra's original one, having undergone some serious refactoring.

In order to keep the code as tidy as possible given the inherent complexity of the algorithm, the helper functions are isolated in an anonymous namespace.

In practice, the three sections below would be concatenated in the given order in a separate source file (say, smoothsort.cpp).

Last but not least, the code uses the long long unsigned type, which is a gcc extension (no, it isn't). On installations where gcc is unavailable, replace with unsigned long (it only affects the maximun size of the arrays which can be sorted).

### Prototypes and forward declarations

/**
**  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 ();
}

### Main body

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;
}

### Helper functions

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
**
**    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;
}
}

## C

/* Adapted from Delphi implementation of Dijkstra's algorithm.
Second argument to smoothsort added. Functions IsAscending, UP, DOWN
replaced by macros.
*/

#include <string.h>

/* Type of items to be sorted; for numeric types, replace by int, float or double */

typedef char* SItem;

/* Comparison function; for numeric types, use (A <= B) */

#define IsAscending(A,B) (strcmp(A,B) <= 0)

#define UP(IA,IB)   temp=IA; IA+=IB+1;   IB=temp;
#define DOWN(IA,IB) temp=IB; IB=IA-IB-1; IA=temp;

static int q,r,
p,b,c,
r1,b1,c1;
static SItem *A;

void _inline sift(void){
int r0, r2, temp;
SItem T;
r0 = r1;
T = A[r0];
while (b1 >= 3) {
r2 = r1-b1+c1;
if (! IsAscending(A[r1-1],A[r2])) {
r2 = r1-1;
DOWN(b1,c1)
}
if (IsAscending(A[r2],T)) b1 = 1;
else {
A[r1] = A[r2];
r1 = r2;
DOWN(b1,c1)
}
}
if (r1 - r0) A[r1] = T;
}

void _inline trinkle(void){
int p1,r2,r3, r0, temp;
SItem T;
p1 = p; b1 = b; c1 = c;
r0 = r1; T = A[r0];
while (p1 > 0) {
while ((p1 & 1)==0) {
p1 >>= 1;
UP(b1,c1)
}
r3 = r1-b1;
if ((p1==1) || IsAscending(A[r3], T)) p1 = 0;
else{
p1--;
if (b1==1) {
A[r1] = A[r3];
r1 = r3;
}
else
if (b1 >= 3) {
r2 = r1-b1+c1;
if (! IsAscending(A[r1-1],A[r2])) {
r2 = r1-1;
DOWN(b1,c1)
p1 <<= 1;
}
if (IsAscending(A[r2],A[r3])) {
A[r1] = A[r3]; r1 = r3;
}
else {
A[r1] = A[r2];
r1 = r2;
DOWN(b1,c1)
p1 = 0;
}
}
}
}
if (r0-r1) A[r1] = T;
sift();
}

void _inline semitrinkle(void){
SItem T;
r1 = r-c;
if (! IsAscending(A[r1],A[r])) {
T = A[r]; A[r] = A[r1]; A[r1] = T;
trinkle();
}
}

void smoothsort(SItem Aarg[], const int N){
int temp;
A=Aarg; /* 0-base array; warning: A is shared by other functions */
q = 1; r = 0; p = 1; b = 1; c = 1;

/* building tree */
while (q < N) {
r1 = r;
if ((p & 7)==3) {
b1 = b; c1 = c; sift();
p = (p+1) >> 2;
UP(b,c) UP(b,c)
}
else if ((p & 3)==1) {
if (q + c < N) {
b1 = b; c1 = c; sift();
}
else trinkle();
DOWN(b,c);
p <<= 1;
while (b > 1) {
DOWN(b,c)
p <<= 1;
}
p++;
}
q++; r++;
}
r1 = r; trinkle();

/* building sorted array */
while (q > 1) {
q--;
if (b==1) {
r--; p--;
while ((p & 1)==0) {
p >>= 1;
UP(b,c)
}
}
else
if (b >= 3) {
p--; r = r-b+c;
if (p > 0) semitrinkle();
DOWN(b,c)
p = (p << 1) + 1;
r = r+c;  semitrinkle();
DOWN(b,c)
p = (p << 1) + 1;
}
/* element q processed */
}
/* element 0 processed */
}

## Delphi

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.