基数排序优化

Question

我正在尝试优化 Radix Sort 代码，因为我觉得它还有空间，因为书籍和网络中的传统代码似乎是彼此的直接副本，而且它们的工作速度非常慢，因为它们采用任意数字，例如10 用于模运算。我已经尽我所能优化了代码，也许我可能错过了一些优化技术。既然如此，请赐教。

优化动机：
http://codercorner.com/RadixSortRevisited.htm
http://stereopsis.com/radix.html
我无法实现文章中的所有优化，主要是超出了我的技能和理解范围，并且没有足够的时间，如果您可以随意实现它们。

编辑 4
这个 Java 版本的基数排序计算 1 次读取中的所有直方图，并且不需要在每次 LSB 排序后用零填充数组 Z，并且如果所有先前的 LSB 都相同，则通常可以跳过排序并跳转到下一个 LSB 排序。像往常一样，这仅适用于 32 位整数，但可以从中创建 64 位版本。

protected static int[] DSC(int A[])// Sorts in descending order
{
    int tmp[] = new int[A.length] ;
    int Z[] = new int[1024] ;
    int i, Jump, Jump2, Jump3, Jump4, swap[] ;

    Jump = A[0] & 255 ;
    Z[Jump] = 1 ;
    Jump2 = ((A[0] >> 8) & 255) + 256 ;
    Z[Jump2] = 1 ;
    Jump3 = ((A[0] >> 16) & 255) + 512 ;
    Z[Jump3] = 1 ;
    Jump4 = (A[0] >> 24) + 768 ;
    Z[Jump4] = 1 ;

    // Histograms creation
    for (i = 1 ; i < A.length; ++i)
    {
        ++Z[A[i] & 255] ;
        ++Z[((A[i] >> 8) & 255) + 256] ;
        ++Z[((A[i] >> 16) & 255) + 512] ;
        ++Z[(A[i] >> 24) + 768] ;
    }

    // 1st LSB Byte Sort
    if( Z[Jump] != A.length )
    {
        Z[0] = A.length - Z[0];
        for (i = 1; i < 256; ++i)
        {
            Z[i] = Z[i - 1] - Z[i];
        }
        for (i = 0; i < A.length; ++i)
        {
            tmp[Z[A[i] & 255]++] = A[i];
        }
        swap = A ; A = tmp ; tmp = swap ;
    }

    // 2nd LSB Byte Sort
    if( Z[Jump2] != A.length )
    {
        Z[256] = A.length - Z[256];
        for (i = 257; i < 512; ++i)
        {
            Z[i] = Z[i - 1] - Z[i];
        }
        for (i = 0; i < A.length; ++i)
        {
            tmp[Z[((A[i] >> 8) & 255) + 256]++] = A[i];
        }
        swap = A ; A = tmp ; tmp = swap ;
    }

    // 3rd LSB Byte Sort
    if( Z[Jump3] != A.length )
    {
        Z[512] = A.length - Z[512];
        for (i = 513; i < 768; ++i)
        {
            Z[i] = Z[i - 1] - Z[i];
        }
        for (i = 0; i < A.length; ++i)
        {
            tmp[Z[((A[i] >> 16) & 255) + 512]++] = A[i];
        }
        swap = A ; A = tmp ; tmp = swap ;
    }

    // 4th LSB Byte Sort
    if( Z[Jump4] != A.length )
    {
        Z[768] = A.length - Z[768];
        for (i = 769; i < Z.length; ++i)
        {
            Z[i] = Z[i - 1] - Z[i];
        }
        for (i = 0; i < A.length; ++i)
        {
            tmp[Z[(A[i] >> 24) + 768]++] = A[i];
        }
        return tmp ;
    }
    return A ;
}

Java 版本使用 != 符号比 == 符号运行得更快

if( Z[Jump] != A.length )
{
    // lines of code
}...

但在 C 中，以下版本平均比使用 != 符号的版本快 25%（使用等号）。您的硬件可能会有不同的反应。

if( Z[Jump] == A.length );
else
{
    // lines of code
}...

下面是 C 代码（我机器上的“long”是 32 位）

long* Radix_2_ac_long(long *A, size_t N, long *Temp)// Sorts in ascending order
{
    size_t Z[1024] = {0};
    long *swp;
    size_t i, Jump, Jump2, Jump3, Jump4;

    // Sort-circuit set-up
    Jump = *A & 255;
    Z[Jump] = 1;
    Jump2 = ((*A >> 8) & 255) + 256;
    Z[Jump2] = 1;
    Jump3 = ((*A >> 16) & 255) + 512;
    Z[Jump3] = 1;
    Jump4 = (*A >> 24) + 768;
    Z[Jump4] = 1;

    // Histograms creation
    for(i = 1 ; i < N ; ++i)
    {
        ++Z[*(A+i) & 255];
        ++Z[((*(A+i) >> 8) & 255) + 256];
        ++Z[((*(A+i) >> 16) & 255) + 512];
        ++Z[(*(A+i) >> 24) + 768];
    }

    // 1st LSB byte sort
    if( Z[Jump] == N );
    else
    {
        for( i = 1 ; i < 256 ; ++i )
        {
            Z[i] = Z[i-1] + Z[i];
        }
        for( i = N-1 ; i < N ; --i )
        {
            *(--Z[*(A+i) & 255] + Temp) = *(A+i);
        }
        swp = A;
        A = Temp;
        Temp = swp;
    }

    // 2nd LSB byte sort
    if( Z[Jump2] == N );
    else
    {
        for( i = 257 ; i < 512 ; ++i )
        {
            Z[i] = Z[i-1] + Z[i];
        }
        for( i = N-1 ; i < N ; --i )
        {
            *(--Z[((*(A+i) >> 8) & 255) + 256] + Temp) = *(A+i);
        }
        swp = A;
        A = Temp;
        Temp = swp;
    }

    // 3rd LSB byte sort
    if( Z[Jump3] == N );
    else
    {
        for( i = 513 ; i < 768 ; ++i )
        {
            Z[i] = Z[i-1] + Z[i];
        }
        for( i = N-1 ; i < N ; --i )
        {
            *(--Z[((*(A+i) >> 16) & 255) + 512] + Temp) = *(A+i);
        }
        swp = A;
        A = Temp;
        Temp = swp;
    }

    // 4th LSB byte sort
    if( Z[Jump4] == N );
    else
    {
        for( i = 769 ; i < 1024 ; ++i )
        {
            Z[i] = Z[i-1] + Z[i];
        }
        for( i = N-1 ; i < N ; --i )
        {
            *(--Z[(*(A+i) >> 24) + 768] + Temp) = *(A+i);
        }
        return Temp;
    }
    return A;
}

编辑 5
排序现在也处理负数。只是对代码进行了一些小的/可忽略的调整。结果它运行得慢了一点，但效果并不显着。用 C 编码，如下（我的系统上的“long”是 32 位）

long* Radix_Sort(long *A, size_t N, long *Temp)
{
    size_t Z[1024] = {0};
    long *swp;
    size_t Jump, Jump2, Jump3, Jump4;
    long i;

    // Sort-circuit set-up
    Jump = *A & 255;
    Z[Jump] = 1;
    Jump2 = ((*A >> 8) & 255) + 256;
    Z[Jump2] = 1;
    Jump3 = ((*A >> 16) & 255) + 512;
    Z[Jump3] = 1;
    Jump4 = ((*A >> 24) & 255) + 768;
    Z[Jump4] = 1;

    // Histograms creation
    for(i = 1 ; i < N ; ++i)
    {
        ++Z[*(A+i) & 255];
        ++Z[((*(A+i) >> 8) & 255) + 256];
        ++Z[((*(A+i) >> 16) & 255) + 512];
        ++Z[((*(A+i) >> 24) & 255) + 768];
    }

    // 1st LSB byte sort
    if( Z[Jump] == N );
    else
    {
        for( i = 1 ; i < 256 ; ++i )
        {
            Z[i] = Z[i-1] + Z[i];
        }
        for( i = N-1 ; i >= 0 ; --i )
        {
            *(--Z[*(A+i) & 255] + Temp) = *(A+i);
        }
        swp = A;
        A = Temp;
        Temp = swp;
    }

    // 2nd LSB byte sort
    if( Z[Jump2] == N );
    else
    {
        for( i = 257 ; i < 512 ; ++i )
        {
            Z[i] = Z[i-1] + Z[i];
        }
        for( i = N-1 ; i >= 0 ; --i )
        {
            *(--Z[((*(A+i) >> 8) & 255) + 256] + Temp) = *(A+i);
        }
        swp = A;
        A = Temp;
        Temp = swp;
    }

    // 3rd LSB byte sort
    if( Z[Jump3] == N );
    else
    {
        for( i = 513 ; i < 768 ; ++i )
        {
            Z[i] = Z[i-1] + Z[i];
        }
        for( i = N-1 ; i >= 0 ; --i )
        {
            *(--Z[((*(A+i) >> 16) & 255) + 512] + Temp) = *(A+i);
        }
        swp = A;
        A = Temp;
        Temp = swp;
    }

    // 4th LSB byte sort and negative numbers sort
    if( Z[Jump4] == N );
    else
    {
        for( i = 897 ; i < 1024 ; ++i )// -ve values frequency starts after index 895, i.e at 896 ( 896 = 768 + 128 ), goes upto 1023
        {
            Z[i] = Z[i-1] + Z[i];
        }
        Z[768] = Z[768] + Z[1023];
        for( i = 769 ; i < 896 ; ++i )
        {
            Z[i] = Z[i-1] + Z[i];
        }
        for( i = N-1 ; i >= 0 ; --i )
        {
            *(--Z[((*(A+i) >> 24) & 255) + 768] + Temp) = *(A+i);
        }
        return Temp;
    }
    return A;
}

编辑 6
下面是指针优化版本（通过指针访问数组位置），它的平均排序时间比上面的版本减少了大约 20%。它还使用 4 个单独的数组来更快地计算地址（我的系统上的“long”是 32 位）。

long* Radix_Sort(long *A, size_t N, long *Temp)
{
    long Z1[256] ;
    long Z2[256] ;
    long Z3[256] ;
    long Z4[256] ;
    long T = 0 ;
    while(T != 256)
    {
        *(Z1+T) = 0 ;
        *(Z2+T) = 0 ;
        *(Z3+T) = 0 ;
        *(Z4+T) = 0 ;
        ++T;
    }
    size_t Jump, Jump2, Jump3, Jump4;

    // Sort-circuit set-up
    Jump = *A & 255 ;
    Z1[Jump] = 1;
    Jump2 = (*A >> 8) & 255 ;
    Z2[Jump2] = 1;
    Jump3 = (*A >> 16) & 255 ;
    Z3[Jump3] = 1;
    Jump4 = (*A >> 24) & 255 ;
    Z4[Jump4] = 1;

    // Histograms creation
    long *swp = A + N;
    long *i = A + 1;
    for( ; i != swp ; ++i)
    {
        ++Z1[*i & 255];
        ++Z2[(*i >> 8) & 255];
        ++Z3[(*i >> 16) & 255];
        ++Z4[(*i >> 24) & 255];
    }

    // 1st LSB byte sort
    if( Z1[Jump] == N );
    else
    {
        swp = Z1+256 ;
        for( i = Z1+1 ; i != swp ; ++i )
        {
            *i = *(i-1) + *i;
        }
        swp = A-1;
        for( i = A+N-1 ; i != swp ; --i )
        {
            *(--Z1[*i & 255] + Temp) = *i;
        }
        swp = A;
        A = Temp;
        Temp = swp;
    }

    // 2nd LSB byte sort
    if( Z2[Jump2] == N );
    else
    {
        swp = Z2+256 ;
        for( i = Z2+1 ; i != swp ; ++i )
        {
            *i = *(i-1) + *i;
        }
        swp = A-1;
        for( i = A+N-1 ; i != swp ; --i )
        {
            *(--Z2[(*i >> 8) & 255] + Temp) = *i;
        }
        swp = A;
        A = Temp;
        Temp = swp;
    }

    // 3rd LSB byte sort
    if( Z3[Jump3] == N );
    else
    {
        swp = Z3 + 256 ;
        for( i = Z3+1 ; i != swp ; ++i )
        {
            *i = *(i-1) + *i;
        }
        swp = A-1;
        for( i = A+N-1 ; i != swp ; --i )
        {
            *(--Z3[(*i >> 16) & 255] + Temp) = *i;
        }
        swp = A;
        A = Temp;
        Temp = swp;
    }

    // 4th LSB byte sort and negative numbers sort
    if( Z4[Jump4] == N );
    else
    {
        swp = Z4 + 256 ;
        for( i = Z4+129 ; i != swp ; ++i )
        {
            *i = *(i-1) + *i;
        }
        *Z4 = *Z4 + *(Z4+255) ;
        swp = Z4 + 128 ;
        for( i = Z4+1 ; i != swp ; ++i )
        {
            *i = *(i-1) + *i;
        }
        swp = A - 1;
        for( i = A+N-1 ; i != swp ; --i )
        {
            *(--Z4[(*i >> 24) & 255] + Temp) = *i;
        }
        return Temp;
    }
    return A;
}

Answer 1

N & 15 , N & 31 , N & 63 .... 等等，这些按位 操作时间最短？

他们是一样的。不要把它弄坏，但是在不知道事情持续多长时间的情况下优化速度可能会很糟糕。即使你知道时间，现在的硬件也非常复杂，而且非常不可预测。你用java编程，那是另一层极其复杂的系统。相同的代码今天可能更快，明天可能更慢。你说approximately 2.232891909840167 times faster。实际上，您使用一组数据对一种硬件和软件配置进行测量，您只能希望测量具有足够的代表性。不幸的是，情况并非总是如此。

我重写了你的函数。它更短更简单，但似乎并不慢。编译器倾向于喜欢不太聪明的代码，因为对于简单的情况有很多优化。对负数的修正不是特别好，不喜欢的可以删掉。它似乎对 8 位和 11 位效果最好，可能是由于缓存大小的原因，看看 rcgldr 的 cmets。

编辑

@ytoamn 你是对的，如果所有内容都在第一个桶中，则循环应该继续，而不是中断。那是一个错误。对于其他更改，我宁愿避免您现在所做的合同。我认为排序功能有三种自然契约。第一个是对原始数组进行排序并返回 null。其次是对原始数组进行排序并返回。第三是返回新的排序数组并保持原始数组不变。我喜欢第一个，因为它的行为是明确的。您现在拥有它的方式应该在文档中添加大警告，原始数组已更改并从函数返回是某些情况，而在其他情况下则不是。我要避免的第二件事是旧的 C 代码风格。如果你只需要它，你应该在循环中定义循环变量。全局定义它会注入可能导致错误的依赖项。而且它在这里没有任何优势，因为正确定义的循环变量最终会共享空间。编译器很清楚作用域，你应该使用你需要的最小作用域。

EDIT2

请直接在我的帖子下发表评论 :-) 局部变量只是堆栈上的地址。您在构造对象时分配内存，这不是这里的情况。至于数组，想想这段代码：

public static void Tst(int[] A) {
    int[] tmp = new int[A.length];
    A[0] = 6;
    A = tmp; // changes what parameter A contains
    A[0] = 7;
}

public static void main(String[] args) {
    int[] A = new int[1];
    A[0] = 5;
    Tst(A);
    System.out.println(A[0]); //prints 6
}

它打印 6。数字 7 仅写入 tmp 数组。 main 中的数组 A 不受影响。

protected static void ASC2(int A[], int bits) {
    int[] origA = A;
    int[] tmp = new int[A.length];
    int[] Z = new int[1 << bits];
    int mask = (1 << bits) - 1;

    for (int shift = 0; shift < 32; shift += bits) {

        if (shift > 0) {
            Arrays.fill(Z, 0);
        }
        for (int i = 0; i < A.length; ++i) {
            Z[(A[i] >> shift) & mask]++;
        }

        if (Z[0] == A.length) {
            continue; // all in first bucket
        }

        Z[Z.length - 1] = A.length - Z[Z.length - 1];
        for (int i = Z.length - 2; i >= 0; --i) {
            Z[i] = Z[i + 1] - Z[i];
        }

        if (shift + bits > 31) { // negative numbers correction
            int halfLength = Z.length / 2;
            int positSum = Z[halfLength];
            int negSum = A.length - positSum;
            if (negSum > 0) {
                for (int i = 0; i < halfLength; ++i) {
                    Z[i] += negSum;
                }
                for (int i = halfLength; i < Z.length; ++i) {
                    Z[i] -= positSum;
                }
            }
        }

        for (int i = 0; i < A.length; ++i) {
            tmp[Z[(A[i] >> shift) & mask]++] = A[i];
        }

        int[] swap = A;
        A = tmp;
        tmp = swap;
    }

    if (A != origA) {
        System.arraycopy(A, 0, origA, 0, A.length);
    }
}

EDIT3

循环展开是一种有效的技术，改善短路真的很好。但是使用数组长度作为常量，你肯定开始太聪明了。如果你硬编码了基本尺寸，为什么不硬编码它：

protected static int[] DSC2(int A[])// sorts in descending order
{
    int tmp[] = new int[A.length];
    int Z[] = new int[256];
    int sample, swap[];

    // 1st LSB byte extraction
    sample = A[0] & 255;
    for (int i = 0; i < A.length; ++i) {
        Z[A[i] & 255]++;
    }

    if (Z[sample] != A.length) {
        Z[0] = A.length - Z[0];
        for (int i = 1; i < Z.length; ++i) {
            Z[i] = Z[i - 1] - Z[i];
        }

        for (int i = 0; i < A.length; ++i) {
            tmp[Z[A[i] & 255]++] = A[i];
        }

        swap = A;
        A = tmp;
        tmp = swap;
        Arrays.fill(Z, 0);
    } else {
        Z[sample] = 0;
    }

    // 2nd LSB byte extraction
    sample = (A[0] >> 8) & 255;
    for (int i = 0; i < A.length; ++i) {
        Z[(A[i] >> 8) & 255]++;
    }

    if (Z[sample] != A.length) {
        Z[0] = A.length - Z[0];
        for (int i = 1; i < Z.length; ++i) {
            Z[i] = Z[i - 1] - Z[i];
        }

        for (int i = 0; i < A.length; ++i) {
            tmp[Z[(A[i] >> 8) & 255]++] = A[i];
        }

        swap = A;
        A = tmp;
        tmp = swap;
        Arrays.fill(Z, 0);
    } else {
        Z[sample] = 0;
    }

    // 3rd LSB byte extraction
    sample = (A[0] >> 16) & 255;
    for (int i = 0; i < A.length; ++i) {
        Z[(A[i] >> 16) & 255]++;
    }

    if (Z[sample] != A.length) {
        Z[0] = A.length - Z[0];
        for (int i = 1; i < Z.length; ++i) {
            Z[i] = Z[i - 1] - Z[i];
        }

        for (int i = 0; i < A.length; ++i) {
            tmp[Z[(A[i] >> 16) & 255]++] = A[i];
        }

        swap = A;
        A = tmp;
        tmp = swap;
        Arrays.fill(Z, 0);
    } else {
        Z[sample] = 0;
    }

    // 4th LSB byte extraction
    sample = (A[0] >> 24) & 255;
    for (int i = 0; i < A.length; ++i) {
        Z[(A[i] >> 24) & 255]++;
    }

    if (Z[sample] != A.length) {
        Z[0] = A.length - Z[0];
        for (int i = 1; i < Z.length; ++i) {
            Z[i] = Z[i - 1] - Z[i];
        }

        for (int i = 0; i < A.length; ++i) {
            tmp[Z[(A[i] >> 24) & 255]++] = A[i];
        }

        A = tmp;
    }

    return A;
}

Answer 2

如果原始数组和临时数组适合缓存，edit 4 版本就足够了。如果数组大小远大于缓存大小，则大部分开销是由于对数组的随机顺序写入。混合 msb/lsb 基数排序可以避免这个问题。例如，根据最高有效字节将数组拆分为 256 个 bin，然后对 256 个 bin 中的每一个进行 lsb 基数排序。这里的想法是一对（原始和临时）bin 将适合缓存，其中随机顺序写入不是问题（对于大多数缓存实现）。

对于 8MB 缓存，目标是每个 bin 的大小

// split array into 256 bins according to most significant byte
void RadixSort(uint32_t * a, size_t count)
{
size_t aIndex[260] = {0};               // count / array
uint32_t * b = new uint32_t [count];    // allocate temp array
size_t i;
    for(i = 0; i < count; i++)          // generate histogram
        aIndex[1+((size_t)(a[i] >> 24))]++;
    for(i = 2; i < 257; i++)            // convert to indices
        aIndex[i] += aIndex[i-1];
    for(i = 0; i < count; i++)          //  sort by msb
        b[aIndex[a[i]>>24]++] = a[i];
    for(i = 256; i; i--)                // restore aIndex
        aIndex[i] = aIndex[i-1];
    aIndex[0] = 0;
    for(i = 0; i < 256; i++)            // radix sort the 256 bins
        RadixSort3(&b[aIndex[i]], &a[aIndex[i]], aIndex[i+1]-aIndex[i]);
    delete[] b;
}

// sort a bin by 3 least significant bytes
void RadixSort3(uint32_t * a, uint32_t *b, size_t count)
{
size_t mIndex[3][256] = {0};            // count / matrix
size_t i,j,m,n;
uint32_t u;
    if(count == 0)
        return;
    for(i = 0; i < count; i++){         // generate histograms
        u = a[i];
        for(j = 0; j < 3; j++){
            mIndex[j][(size_t)(u & 0xff)]++;
            u >>= 8;
        }       
    }
    for(j = 0; j < 3; j++){             // convert to indices
        m = 0;
        for(i = 0; i < 256; i++){
            n = mIndex[j][i];
            mIndex[j][i] = m;
            m += n;
        }       
    }
    for(j = 0; j < 3; j++){             // radix sort
        for(i = 0; i < count; i++){     //  sort by current lsb
            u = a[i];
            m = (size_t)(u>>(j<<3))&0xff;
            b[mIndex[j][m]++] = u;
        }
        std::swap(a, b);                //  swap ptrs
    }
}

---

经典 lsb 基数排序的示例代码：

使用 8、8、8、8 位字段的 C++ lsb 基数排序示例：

typedef unsigned int uint32_t;

void RadixSort(uint32_t * a, size_t count)
{
size_t mIndex[4][256] = {0};            // count / index matrix
uint32_t * b = new uint32_t [count];    // allocate temp array
size_t i,j,m,n;
uint32_t u;
    for(i = 0; i < count; i++){         // generate histograms
        u = a[i];
        for(j = 0; j < 4; j++){
            mIndex[j][(size_t)(u & 0xff)]++;
            u >>= 8;
        }       
    }
    for(j = 0; j < 4; j++){             // convert to indices
        m = 0;
        for(i = 0; i < 256; i++){
            n = mIndex[j][i];
            mIndex[j][i] = m;
            m += n;
        }       
    }
    for(j = 0; j < 4; j++){             // radix sort
        for(i = 0; i < count; i++){     //  sort by current lsb
            u = a[i];
            m = (size_t)(u>>(j<<3))&0xff;
            b[mIndex[j][m]++] = u;
        }
        std::swap(a, b);                //  swap ptrs
    }
    delete[] b;
}

使用 16,16 位字段的示例 C++ 代码：

typedef unsigned int uint32_t;

uint32_t * RadixSort(uint32_t * a, size_t count)
{
size_t mIndex[2][65536] = {0};          // count / index matrix
uint32_t * b = new uint32_t [count];    // allocate temp array
size_t i,j,m,n;
uint32_t u;
    for(i = 0; i < count; i++){         // generate histograms
        u = a[i];
        for(j = 0; j < 2; j++){
            mIndex[j][(size_t)(u & 0xffff)]++;
            u >>= 16;
        }
    }
    for(j = 0; j < 2; j++){             // convert to indices
        m = 0;
        for(i = 0; i < 65536; i++){
            n = mIndex[j][i];
            mIndex[j][i] = m;
            m += n;
        }
    }       
    for(j = 0; j < 2; j++){             // radix sort
        for(i = 0; i < count; i++){     //  sort by current lsb
            u = a[i];
            m = (size_t)(u>>(j<<4))&0xffff;
            b[mIndex[j][m]++] = u;
        }
        std::swap(a, b);                //  swap ptrs
    }
    delete[] b;
    return(a);
}