 | Wikipedia has related information atSmoothsort |
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 thelong long unsigned type,which is agcc extension (no, it isn't). On installations wheregcc is unavailable, replace withunsigned long (it only affects the maximum size of the arrays which can be sorted).
/** ** 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<typenameT>voidsmoothsort(T*,unsigned)throw();namespace/** ** Helper function's local namespace (declarations). ** **/{classLeonardoNumber/** ** Helper class for manipulation of Leonardo numbers ** **/{public:/** Default ctor. **/LeonardoNumber(void)throw():b(1),c(1){return;}/** Copy ctor. **/LeonardoNumber(constLeonardoNumber&_l)throw():b(_l.b),c(_l.c){return;}/** ** Return the "gap" between the actual Leonardo number and the ** preceding one. **/unsignedgap(void)constthrow(){returnb-c;}/** Perform an "up" operation on the actual number. **/LeonardoNumber&operator++(void)throw(){unsigneds=b;b=b+c+1;c=s;return*this;}/** Perform a "down" operation on the actual number. **/LeonardoNumber&operator--(void)throw(){unsigneds=c;c=b-c-1;b=s;return*this;}/** Return "companion" value. **/unsignedoperator~(void)constthrow(){returnc;}/** Return "actual" value. **/operatorunsigned(void)constthrow(){returnb;}private:unsignedb;/** Actual number. **/unsignedc;/** Companion number. **/};/** Perform a "sift up" operation. **/template<typenameT>inlinevoidsift(T*,unsigned,LeonardoNumber)throw();/** Perform a "semi-trinkle" operation. **/template<typenameT>inlinevoidsemitrinkle(T*,unsigned,unsignedlonglong,LeonardoNumber)throw();/** Perform a "trinkle" operation. **/template<typenameT>inlinevoidtrinkle(T*,unsigned,unsignedlonglong,LeonardoNumber)throw();}
template<typenameT>voidsmoothsort(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;unsignedlonglongp=1;LeonardoNumberb;for(unsignedq=0;++q<_n;++p)if(p%8==3){sift<T>(_m,q-1,b);++++b;p>>=2;}elseif(p%4==1){if(q+~b<_n)sift<T>(_m,q-1,b);elsetrinkle<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;elseif(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<typenameT>inlinevoidsift(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. ** ** **/{unsignedr2;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<typenameT>inlinevoidsemitrinkle(T*_m,unsigned_r,unsignedlonglong_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<typenameT>inlinevoidtrinkle(T*_m,unsigned_r,unsignedlonglong_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;elseif(_b==1){swap(_m[_r],_m[_r-_b]);_r-=_b;}elseif(_b>=3){unsignedr2=_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) */typedefconstchar*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;} */staticintcmp(value_ta,value_tb){returnstrcmp(a,b);}/* End user-defined types and functions */structstate{value_t*a;size_tq,r,p,b,c,r1,b1,c1;};#if __STDC_VERSION__ < 199901L#define inline/* Silently prune inline qualifiers away for ANSI C */#endifstaticinlinevoidup(size_t*a,size_t*b){size_ttmp;tmp=*a;*a+=*b+1;*b=tmp;}staticinlinevoiddown(size_t*a,size_t*b){size_ttmp;tmp=*b;*b=*a-*b-1;*a=tmp;}staticvoidsift(structstate*s){size_tr0,r2;value_ttmp;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;}}staticvoidtrinkle(structstate*s){size_tp1,r0,r2,r3;value_ttmp;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;}elseif(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);}staticvoidsemitrinkle(structstate*s){value_ttmp;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);}}voidsmoothsort(value_t*a,size_tn){structstates;size_ttmp;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;}elseif((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);}}elseif(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;} */
unitUSmoothsort;{ Delphi implementation of Dijkstra's algorithm }interfacetypeTItem=Double;{ data type }functionIsAscending(v1,v2:TItem):boolean;{ comparison function }{ sorting function }procedureSmoothSort(varA:arrayofTItem);implementation{ customizable comparison function }functionIsAscending(v1,v2:TItem):boolean;beginresult:=v1<=v2;end;{ implementation of Djikstra's algorithm }procedureSmoothSort(varA:arrayofTItem);varq,r,p,b,c,r1,b1,c1,N:integer;procedureup(varvb,vc:integer);vartemp:integer;begintemp:=vb;vb:=vb+vc+1;vc:=temp;end;proceduredown(varvb,vc:integer);vartemp:integer;begintemp:=vc;vc:=vb-vc-1;vb:=temp;end;proceduresift;varr0,r2:integer;T:TItem;beginr0:=r1;T:=A[r0];whileb1>=3dobeginr2:=r1-b1+c1;ifnotIsAscending(A[r1-1],A[r2])thenbeginr2:=r1-1;down(b1,c1);end;ifIsAscending(A[r2],T)thenb1:=1elsebeginA[r1]:=A[r2];r1:=r2;down(b1,c1);end;end;ifr1<>r0thenA[r1]:=T;end;proceduretrinkle;varp1,r2,r3,r0:integer;T:TItem;beginp1:=p;b1:=b;c1:=c;r0:=r1;T:=A[r0];whilep1>0dobeginwhile(p1and1)=0dobeginp1:=p1shr1;up(b1,c1);end;r3:=r1-b1;if(p1=1)orIsAscending(A[r3],T)thenp1:=0else{p1>1}beginp1:=p1-1;ifb1=1thenbeginA[r1]:=A[r3];r1:=r3;endelseifb1>=3thenbeginr2:=r1-b1+c1;ifnotIsAscending(A[r1-1],A[r2])thenbeginr2:=r1-1;down(b1,c1);p1:=p1shl1;end;ifIsAscending(A[r2],A[r3])thenbeginA[r1]:=A[r3];r1:=r3;endelsebeginA[r1]:=A[r2];r1:=r2;down(b1,c1);p1:=0;end;end;{if b1>=3}end;{if p1>1}end;{while}ifr0<>r1thenA[r1]:=T;sift;end;proceduresemitrinkle;varT:TItem;beginr1:=r-c;ifnotIsAscending(A[r1],A[r])thenbeginT:=A[r];A[r]:=A[r1];A[r1]:=T;trinkle;end;end;beginN:=length(A);q:=1;r:=0;p:=1;b:=1;c:=1;//building treewhileq<Ndobeginr1:=r;if(pand7)=3thenbeginb1:=b;c1:=c;sift;p:=(p+1)shr2;up(b,c);up(b,c);endelseif(pand3)=1thenbeginifq+c<Nthenbeginb1:=b;c1:=c;sift;endelsetrinkle;down(b,c);p:=pshl1;whileb>1dobegindown(b,c);p:=pshl1;end;p:=p+1;end;q:=q+1;r:=r+1;end;r1:=r;trinkle;//bulding sorted arraywhileq>1dobeginq:=q-1;ifb=1thenbeginr:=r-1;p:=p-1;while(pand1)=0dobeginp:=pshr1;up(b,c);end;endelseifb>=3thenbeginp:=p-1;r:=r-b+c;ifp>0thensemitrinkle;down(b,c);p:=pshl1+1;r:=r+c;semitrinkle;down(b,c);p:=pshl1+1;end;//element q is doneend;//element 0 is doneend;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.staticfinalintLP[]={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.};publicstatic<CextendsComparable<?superC>>voidsort(C[]m,intlo,inthi){inthead=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)==3intp=1;// the bitmap of the current standard concatenation >> pshiftintpshift=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 1if(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 neededinttrail=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 themtrinkle(m,p>>>1,pshift+1,head-LP[pshift]-1,true);trinkle(m,p,pshift,head-1,true);}head--;}}privatestatic<CextendsComparable<?superC>>voidsift(C[]m,intpshift,inthead){// 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.Cval=m[head];while(pshift>1){intrt=head-1;intlf=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;}privatestatic<CextendsComparable<?superC>>voidtrinkle(C[]m,intp,intpshift,inthead,booleanisTrusty){Cval=m[head];while(p!=1){intstepson=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){intrt=head-1;intlf=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;inttrail=Integer.numberOfTrailingZeros(p&~1);p>>>=trail;pshift+=trail;isTrusty=false;}if(!isTrusty){m[head]=val;sift(m,pshift,head);}}