/ray/src/lib/tl/sort.h

Go to the documentation of this file.

/*
 * lib/tl/sort.h 
 * 
 * Array sorting templates. 
 * 
 * Copyright (c) 2004 by Wolfgang Wieser ] wwieser (a) gmx <*> de [ 
 * 
 * This file may be distributed and/or modified under the terms of the 
 * GNU General Public License version 2 as published by the Free Software 
 * Foundation. (See COPYING.GPL for details.)
 * 
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 * 
 */

#ifndef _TemplateLibrary_SORT_H_
#define _TemplateLibrary_SORT_H_ 1

#include <lib/sconfig.h>    /* MUST be first */

#include <lib/tl/defop.h>


template<typename T,typename _OP>inline void _SortSwap(
    T *a,size_t i0,size_t i1,T &tmp,const _OP &op)
{
    op.ass(tmp,op.idx(a,i0));
    op.ass(op.idx(a,i0),op.idx(a,i1));
    op.ass(op.idx(a,i1),tmp);
}

template<typename T,typename _OP>void InsertionSortBS(
    T *a,size_t nelem,
    const _OP &op=TLDefaultOperators_CDT<T>())
{
    T tmp;
    if(nelem>=2 && op.lt(op.idx(a,1),op.idx(a,0)))  // a[1]<a[0]
    {  _SortSwap(a,0,1,tmp,op);  }
    for(size_t i=2; i<nelem; i++)
    {
        if(op.le(op.idx(a,i-1),op.idx(a,i)))  continue;  // if(a[i-1]<=a[i])
        // Need to move a[i]; find insertion index. 
        T &key=op.idx(a,i);
        size_t aa=0,bb=i-1;  // bb>aa here as we start with i=2. 
        while(bb-aa>1)
        {  size_t m=(aa+bb)/2;  op.le(op.idx(a,m),key) ? aa=m : bb=m;  }
        size_t ii = (!aa && op.lt(key,op.idx(a,0))) ? 0 : bb;
        // Move the elements (ii to i). 
        op.ass(tmp,op.idx(a,i));   // tmp=a[i];
        for(size_t j=i; j>ii; j--)  op.ass(op.idx(a,j),op.idx(a,j-1));
        op.ass(op.idx(a,ii),tmp);
    }
}

template<typename T,typename _OP>void InsertionSort(
    T *a,size_t nelem,
    const _OP &op=TLDefaultOperators_CDT<T>())
{
    T tmp;
    for(size_t i=1; i<nelem; i++)
    {
        if(op.le(op.idx(a,i-1),op.idx(a,i)))  continue;  // if(a[i-1]<=a[i])
        // Need to move a[i]. 
        op.ass(tmp,op.idx(a,i));   // tmp=a[i];
        size_t j=i;
        // See if it is smaller than a[0] (to get around index check in loop). 
        if(op.lt(op.idx(a,i),op.idx(a,0)))   // if(a[i]<a[0])
            for(; j; j--)  op.ass(op.idx(a,j),op.idx(a,j-1));
        else
            do { op.ass(op.idx(a,j),op.idx(a,j-1));  --j; }
            while(op.lt(tmp,op.idx(a,j-1)));
        // a[i] now moved. 
        op.ass(op.idx(a,j),tmp);
    }
}


template<typename T,typename _OP>void HeapSort(
    T *_a,size_t nelem,
    const _OP &op=TLDefaultOperators_CDT<T>())
{
    // Heap addressing is shifted by 1 element so that address calculation 
    // gets easier. 
    T *a=(T*)(((char*)_a)-op.size());
    
    T tmp;
    // Build up heap: Largest element at the root. Use upheap operation. 
    for(size_t i=2; i<=nelem; i++)
    {
        // upheap(i) with heap size i. 
        // For simple types, these 3 possibilities are equally fast. 
        #if 1
            if(op.le(op.idx(a,i),op.idx(a,i/2)))  continue;
            size_t j=i,k=j/2;
            op.ass(tmp,op.idx(a,j));
            do
            {  op.ass(op.idx(a,j),op.idx(a,k)); j=k; k/=2;  }
            while(j>1 && op.lt(op.idx(a,k),tmp));
            op.ass(op.idx(a,j),tmp);
        #elif 0
            size_t j=i,k=j/2;
            while(j>1 && op.lt(op.idx(a,k),op.idx(a,j)))
            {  _SortSwap(a,j,k,tmp,op); j=k; k/=2;  }
        #else
            size_t j=i,k=j/2;
            op.ass(tmp,op.idx(a,j));
            while(j>1 && op.lt(op.idx(a,k),tmp))
            {  op.ass(op.idx(a,j),op.idx(a,k)); j=k; k/=2;  }
            op.ass(op.idx(a,j),tmp);
        #endif
    }
    
    // Remove elements from the heap: 
    // This is the part which should be improved. 
    for(size_t n=nelem; n>1; )
    {
        op.ass(tmp,op.idx(a,n));
        op.ass(op.idx(a,n--),op.idx(a,1));
        // downheap(1) with heap size n and element 1 already in tmp. 
        size_t k=1,e=n/2;
        while(k<=e)
        {
            #if 1  /* The second version is not faster, either. */
                size_t j=k+k;
                if(j<n && op.lt(op.idx(a,j),op.idx(a,j+1)))  ++j;
                if(op.le(op.idx(a,j),tmp))  break;
                op.ass(op.idx(a,k),op.idx(a,j));  k=j;
            #else
                size_t j=k+k;
                char *t=(char*)&op.idx(a,j),*t1=t+op.size();
                if(j<n && op.lt(*(T*)t,*(T*)t1))
                {  ++j;  t=t1;  }
                if(op.le(*(T*)t,tmp))  break;
                op.ass(op.idx(a,k),*(T*)t);  k=j;
            #endif
        }
        op.ass(op.idx(a,k),tmp);
    }
}


template<typename T,typename _OP>inline void _QuickSort_Final(
    T *a,size_t nelem,T &tmp,
    const _OP &op)
{
    if(nelem<=2)
    {
        if(nelem<2)  return;
        if(op.lt(op.idx(a,1),op.idx(a,0)))  // a[1]<a[0]
            _SortSwap(a,0,1,tmp,op);
        return;
    }
    //if(nelem==3)  {  Sort3() }
    
    // One could either use one insertion sort step at the end for the 
    // complete array or calling it directly here for each sub-array. 
    // Tests showed that it is slightly faster here (in amount of compares) 
    // and also CPU caching should be better. 
    InsertionSort(a,nelem,op);
}


template<typename T,typename _OP>void _DumbQuickSort_R(
    T *a,size_t nelem,size_t insertion_sort_thresh,
    const _OP &op)
{
    T tmp;
    
    // Sort small subarrays using insertion sort. 
    if(nelem<insertion_sort_thresh)  // insertion_sort_thresh>=2 (!!)
    {  _QuickSort_Final(a,nelem,tmp,op);  return;  }
    
    // MUST USE nelem>=2 here. 
    // Choose element to place: always the first element. 
    
    // The famous quicksort array walk: 
    T &key=op.idx(a,0);
    size_t aa=1,bb=nelem-1;
    while(aa<bb)
    {
        #if 1  /* Seems slightly faster for simple types. */
            while(aa<bb && op.le(op.idx(a,aa),key))  ++aa;
            while(aa<bb && op.le(key,op.idx(a,bb)))  --bb;
            if(aa>=bb)  break;
        #else
            for(;;++aa)
            {   if(aa>=bb)  goto breakout;
                if(op.lt(key,op.idx(a,aa)))  break;  }
            for(;;--bb)
            {   if(aa>=bb)  goto breakout;
                if(op.lt(op.idx(a,bb),key))  break;  }
            //Assert(aa<bb);
        #endif
        _SortSwap(a,aa++,bb--,tmp,op);
    } // breakout:;    // <-- Un-comment when using second variant above. 
    
    // Place key element: 
    if(op.gt(op.idx(a,bb),key))
    { Assert(aa==bb); --bb;  }
    _SortSwap(a,0,bb,tmp,op);
    
    // Recurse: 
    _DumbQuickSort_R(a,bb,insertion_sort_thresh,op);
    _DumbQuickSort_R(&op.idx(a,bb+1),nelem-(bb+1),insertion_sort_thresh,op);
}

template<typename T,typename _OP>void DumbQuickSort(
    T *a,size_t nelem,
    size_t insertion_sort_thresh=16,
    const _OP &op=TLDefaultOperators_CDT<T>())
{
    _DumbQuickSort_R(a,nelem,
        insertion_sort_thresh<2 ? 2 : insertion_sort_thresh,op);
    // The insertion sort step is done at end of recursion for small 
    // subarray and not for all the array. 
    //InsertionSort(a,nelem,op);
}


#define TLIB_QUICKSORT_TAIL_CALL_OPTIMIZE 1
template<typename T,typename _OP>void _CleverQuickSort_R(
    T *a,size_t nelem,size_t insertion_sort_thresh,
    const _OP &op)
{
    T tmp,key;
    
#if TLIB_QUICKSORT_TAIL_CALL_OPTIMIZE
for(;;) {
#endif
    // Sort small subarrays using insertion sort. 
    if(nelem<insertion_sort_thresh)  // insertion_sort_thresh>=3 (!!)
    {  _QuickSort_Final(a,nelem,tmp,op);  return;  }
    
    // MUST USE nelem>=3 here. 
    // Choose element to place: median of first, middle and last element. 
    // Sort these 3 elements here to provide marks so that we do not have 
    // to check the indices in the array walk. 
    // In the end, the largest of the 3 elements is at index bb, the smallest 
    // at index 1 (not 0) and the median ist saved in key. 
    size_t aa=0,bb=nelem-1,m=(aa+bb)/2;  // Start layout:  aa .. m .. bb
                                         // Target layout: m aa ... bb
    if(op.le(op.idx(a,aa),op.idx(a,bb)))  // aa<=bb
    {
        if(op.lt(op.idx(a,m),op.idx(a,aa)))  // m<aa  -> order: m aa bb
        {
            // largest (bb) at end: OK
            // median (aa) at beginning: OK
            op.ass(key,op.idx(a,aa));
            op.ass(op.idx(a,aa),op.idx(a,m));   // from here: swap m,1
            op.ass(op.idx(a,m),op.idx(a,1));
            op.ass(op.idx(a,1),op.idx(a,aa));
        }
        else if(op.le(op.idx(a,m),op.idx(a,bb)))  // m<=bb  -> order: aa m bb
        {
            // largest (bb) at end: OK
            op.ass(key,op.idx(a,m));
            op.ass(op.idx(a,m),op.idx(a,1));
            op.ass(op.idx(a,1),op.idx(a,aa));
            //op.ass(op.idx(a,aa),op.idx(a,m));  // (not needed)
        }
        else  // order: aa bb m
        {
            op.ass(key,op.idx(a,bb));
            op.ass(op.idx(a,bb),op.idx(a,m));  // largest to end
            op.ass(op.idx(a,m),op.idx(a,1));
            op.ass(op.idx(a,1),op.idx(a,aa));
        }
    }
    else if(op.lt(op.idx(a,bb),op.idx(a,m)))  // bb<m
    {
        if(op.le(op.idx(a,m),op.idx(a,aa)))  // m<=aa  -> order bb m aa
        {
            op.ass(key,op.idx(a,m));
            op.ass(op.idx(a,m),op.idx(a,1));
            op.ass(op.idx(a,1),op.idx(a,bb));
            op.ass(op.idx(a,bb),op.idx(a,aa));
        }
        else  // order: bb aa m
        {
            op.ass(key,op.idx(a,aa));
            op.ass(op.idx(a,aa),op.idx(a,m));
            op.ass(op.idx(a,m),op.idx(a,1));
            op.ass(op.idx(a,1),op.idx(a,bb));
            op.ass(op.idx(a,bb),op.idx(a,aa));
        }
    }
    else  // order: m bb aa
    {
        op.ass(key,op.idx(a,bb));
        op.ass(op.idx(a,bb),op.idx(a,aa));
        _SortSwap(a,m,1,tmp,op);
        /*op.ass(op.idx(a,aa),op.idx(a,m));
        op.ass(op.idx(a,m),op.idx(a,1));
        op.ass(op.idx(a,1),op.idx(a,aa));*/
    }
    
    // The famous quicksort array walk: 
    // We must use lt() and not le() because of the marks (to avoid the 
    // index check). 
    ++aa;
    while(aa<bb)
    {
        while(op.lt(op.idx(a,aa),key))  ++aa;
        while(op.lt(key,op.idx(a,bb)))  --bb;
        if(aa>=bb)  break;
        _SortSwap(a,aa++,bb--,tmp,op);
    }
    
    // Place key element: 
    if(op.gt(op.idx(a,bb),key))
    { Assert(aa==bb); --bb;  }
    op.ass(op.idx(a,0),op.idx(a,bb));
    op.ass(op.idx(a,bb),key);
    
#if TLIB_QUICKSORT_TAIL_CALL_OPTIMIZE
    // Recurse on the smaller subarray and do the larger one by tail 
    // recursion elimination. This guarantees O(log(n)) memory 
    // consumption. 
    aa=nelem-(bb+1);
    if(bb<aa)
    {
        _CleverQuickSort_R(a,bb,insertion_sort_thresh,op);
        a=&op.idx(a,bb+1);  nelem=aa;
    } else {
        _CleverQuickSort_R(&op.idx(a,bb+1),aa,insertion_sort_thresh,op);
        nelem=bb;
    }
} // end for(;;) {
#else
    // Recurse: (simple version) 
    _CleverQuickSort_R(a,bb,insertion_sort_thresh,op);
    _CleverQuickSort_R(&op.idx(a,bb+1),nelem-(bb+1),insertion_sort_thresh,op);
#endif
}

template<typename T,typename _OP>void CleverQuickSort(
    T *a,size_t nelem,
    size_t insertion_sort_thresh=16,
    const _OP &op=TLDefaultOperators_CDT<T>())
{
    _CleverQuickSort_R(a,nelem,
        insertion_sort_thresh<3 ? 3 : insertion_sort_thresh,op);
    // The insertion sort step is done at end of recursion for small 
    // subarray and not for all the array. 
    //InsertionSort(a,nelem,op);
}

#endif  /* _TemplateLibrary_SORT_H_ */

---