Le tri à paniers

#pragma	once

#ifndef __RADIXSORT_H__
#define __RADIXSORT_H__

#define LITTLE_ENDIAN
namespace radix {

    // utility meta function to know type attributes
    template <typename T> inline bool IsFloat();
    template <typename T> inline bool IsSigned();

    // two floating types.
    template <>           inline bool IsFloat<float>()     { return true;}
    template <>           inline bool IsFloat<double>()    { return true;}

    // others are integers.
    template <typename T> inline bool IsFloat()            { return false;}


    // floating types are signed.
    template <> inline bool IsSigned<float>()              { return true; };
    template <> inline bool IsSigned<double>()             { return true; };

    // unsigned types.
    template <> inline bool IsSigned<unsigned short>()     { return false; };
    template <> inline bool IsSigned<unsigned int>()       { return false; };
    template <> inline bool IsSigned<unsigned long>()      { return false; };
    template <> inline bool IsSigned<unsigned long long>() { return false; };
    
    // signed types.
    template <> inline bool IsSigned<signed short>()       { return true; };
    template <> inline bool IsSigned<signed int>()         { return true; };
    template <> inline bool IsSigned<signed long>()        { return true; };
    template <> inline bool IsSigned<signed long long>()   { return true; };
    
    // we need to fix endianess ordering.
    // LE: LSB to MSB, BE: MSB to LSB
    template< typename T>
    inline unsigned int FixMsbEndianess( unsigned int b)
    {
        #ifdef LITTLE_ENDIAN
        return b;
        #else
        return sizeof(T)-1-b;
        #endif
    }

    ///////////////////////////////////////////////////////////////////////////////////
    // Radix sort
    // in: input is the array data to sort
    // in: nb is the array data size
    // out: output, will contains sorted indices
    template < typename T, typename IndexT>
    void Sort(const T* input, unsigned int nb, IndexT* output)
    {
        // Checkings
        if(!input || !output || !nb)	return;

        IndexT* output2 = (IndexT*)alloca(sizeof(IndexT)*nb);

        // Allocate histograms & offsets on the stack
        unsigned int mHistogram[256*sizeof(T)];

        // Create histograms (counters). Counters for all passes are created in one run.
        // Pros:	read input buffer once instead of four times
        // Cons:	mHistogram is sizeof(T)*256 instead of 1Kb

        /* Clear counters */
        memset(mHistogram, 0, 256*sizeof(T)*sizeof(unsigned int));												
        
        unsigned int* h0,*h1,*h2,*h3,*h4,*h5,*h6,*h7;        
        /* Prepare to count */		
        h0= &mHistogram[256*0];		    
        h1= &mHistogram[256*1];	        
        if (sizeof(T) > 2)
        {
            h2= &mHistogram[256*2];	    
            h3= &mHistogram[256*3];	    
        }
        if (sizeof(T) == 8)
        {
            h4= &mHistogram[256*4];	    
            h5= &mHistogram[256*5];	
            h6= &mHistogram[256*6];	
            h7= &mHistogram[256*7];	
        }

        // count each byte in the associate histogram
        for (unsigned int idx = 0;idx<nb;++idx)
        {																							
            /* Read input input in previous sorted order */										
            const T& Val = (T)input[idx];			
            unsigned char* p = (unsigned char*)&Val;	

            /* Create histograms */																	
            h0[*p++]++;	h1[*p++]++;	
            if (sizeof(T) > 2)
            {
                h2[*p++]++;	h3[*p++]++;											
            }
            if (sizeof(T) == 8)
            {
                h4[*p++]++;	h5[*p++]++;	h6[*p++]++;	h7[*p++]++;											
            }
        }				

        // Compute #negative values involved if needed
        unsigned int NbNegativeValues = 0;
        if (IsSigned<T>())
        {
            // An efficient way to compute the number of negatives values we'll have to deal with is simply to sum the 128
            // last values of the last histogram. Last histogram because that's the one for the Most Significant Byte,
            // responsible for the sign. 128 last values because the 128 first ones are related to positive numbers.
            unsigned int* h_msb;
            h_msb = &mHistogram[256*FixMsbEndianess<T>(sizeof(T)-1)];

            for(unsigned int i=128;i<256;i++)	NbNegativeValues += h_msb[i];	// 768 for last histogram, 128 for negative part
        }

        // we'll swap buffer for each byte done.
        IndexT* outputDst = output2;
        IndexT* inputSrc = output;

        // Radix sort, j is the pass number (0=LSB, sizeof(T)-1=MSB)
        for(unsigned int j=0;j<sizeof(T);j++)
        {
            unsigned int* CurCount = &mHistogram[ FixMsbEndianess<T>(j) * 256 ];													

            /* Check pass validity */																	

            /* If all values have the same byte, sorting is useless. */									
            /* It may happen when sorting bytes or words instead of dwords. */							
            /* This routine actually sorts words faster than dwords, and bytes */						
            /* faster than words. Standard running time (O(4*n))is reduced to O(2*n) */					
            /* for words and O(n) for bytes. Running time for floats depends on actual values... */		

            unsigned char UniqueVal;
            /* Get first byte */	
            UniqueVal = *(((unsigned char*)input)+FixMsbEndianess<T>(j));		

            /* Reset flag. The sorting pass is supposed to be performed. (default) */					
            /* Check that byte's counter, j == 0 must be done to initialize output */																
            unsigned int PerformPass = (j == 0) || CurCount[UniqueVal]!=nb;

            // Sometimes the fourth (negative) pass is skipped because all numbers are negative and the MSB is 0xFF (for example). This is
            // not a problem, numbers are correctly sorted anyway.
            if(PerformPass)
            {
                unsigned int mOffset[256];
                
                // Should we care about negative values? only for the msb
                if ( !IsSigned<T>()  || j!=sizeof(T)-1 )
                {
                    // use a register, fix memory load hit store issue
                    unsigned int lastOffset = 0;

                    // Here we deal with positive values only
                    mOffset[0] = 0;

                    // Create offsets
                    for(edU32 i=1;i<256;i++)
                        lastOffset = (mOffset[i] = lastOffset + CurCount[i-1]);
                }
                else
                {
                    // This is a special case to correctly handle negative integers. They're sorted in the right order but at the wrong place.

                    // use a register, fix memory load hit store issue
                    unsigned int lastOffset = NbNegativeValues;

                    // Create biased offsets, in order for negative numbers to be sorted as well
                    mOffset[0] = NbNegativeValues;	
                    // First positive number takes place after the negative ones
                    for(unsigned int i=1;i<128;i++)		
                        lastOffset = (mOffset[i] = lastOffset + CurCount[i-1]);	// 1 to 128 for positive numbers

                    // floating type case differ from integer
                    if (!IsFloat<T>())                    
                    {
                        // use a register, fix memory load hit store issue
                        unsigned int lastOffset = 0;

                        // Fixing the wrong place for negative values
                        mOffset[128] = 0;

                        for(unsigned int i=129;i<256;i++)			
                            lastOffset = (mOffset[i] = lastOffset + CurCount[i-1]);
                    }
                    else
                    {
                        // use a register, fix memory load hit store issue
                        unsigned int lastOffset = 0;

                        // We must reverse the sorting order for negative numbers!
                        mOffset[255] = 0;

                        for(unsigned int i=0;i<127;i++)		
                            lastOffset = (mOffset[254-i] = lastOffset + CurCount[255-i]);	// Fixing the wrong order for negative values
                        for(unsigned int i=128;i<256;i++)	
                            mOffset[i] += CurCount[i];				// Fixing the wrong place for negative values
                    }
                }

                // Perform Radix Sort
                unsigned char* InputBytes	= (unsigned char*)input;

                InputBytes += FixMsbEndianess<T>(j);

                // LSB is special, we initialize fake previous ordering
                if (j==0)
                {
                    for (unsigned int idx = 0;idx<nb;++idx)
                        outputDst[mOffset[ InputBytes[idx*sizeof(T)] ]++] = idx;
                }
                else
                {
                    // common case, only MSB of floating type need something else
                    if (!IsFloat<T>() || j!=sizeof(T)-1)
                    {
                        for (unsigned int idx = 0;idx<nb;++idx)
                        {
                            unsigned int real_idx = inputSrc[idx];
                            outputDst[mOffset[ InputBytes[real_idx*sizeof(T) ]]++] = real_idx;
                        }
                    }
                    else
                    {
                        // Perform Radix Sort
                        for(unsigned int i=0;i<nb;i++)
                        {
                            unsigned int real_idx = inputSrc[i];

                            unsigned int Radix = (*((unsigned int*)input+real_idx))>>24;								// Radix byte, same as above. AND is useless here (edU32).
                            // ### cmp to be killed. Not good. Later.
                            if(Radix<128)		outputDst[mOffset[Radix]++] = real_idx;		// Number is positive, same as above
                            else				outputDst[--mOffset[Radix]] = real_idx;		// Number is negative, flip the sorting order
                        }
                    }
                }

                // Swap pointers for next pass. Valid indices - the most recent ones - are in mIndices after the swap.
                std::swap(outputDst,inputSrc);
            }
            else
            {
                // The pass is useless, yet we still have to reverse the order of current list if all values are negative.
                if(IsFloat<T>() && (j==sizeof(T)-1) && UniqueVal>=128)
                {
                    for(unsigned int i=0;i<nb;i++)	outputDst[i] = inputSrc[nb-i-1];

                    // Swap pointers for next pass. Valid indices - the most recent ones - are in mIndices after the swap.
                    std::swap(outputDst,inputSrc);
                }
            }
        }

        // may happen if not all byte need sorting.
        if (inputSrc != output)
        {
            memcpy(output,inputSrc,sizeof(IndexT)*nb);
        }
    }
}
#endif