Algorithm Implementation/Sorting/Smoothsort

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 maximum size of the arrays which can be sorted).

Prototypes and forward declarations

edit
 
   /**
    **  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

edit
 
  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

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

 */

Delphi

edit
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.

Java

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