快速排序无限循环数值食谱

Quick Sort infinite loop Numerical Recipies

我正在尝试将一些代码从 MATLAB 导出到 MEX 文件 (a.la C)。对于我要导出的代码段,我需要实现一种可以对二维数组进行排序的算法。打开数字食谱的副本,我找到了基于索引的排序。我已经用零索引实现了这个函数,以一种我稍后将简化的简单方式。目前的问题是,当我 运行 我的函数传递某些随机数集时,它将进入无限循环。例如集合:

0.8147,
0.9058,
0.127,
0.9134,
0.6324,
0.0975,
0.2785,
0.5469,
0.9575,
0.9649,
0.1576,
0.9706,
0.9572,
0.4854,
0.8003,
0.1419,
0.4218,
0.9157,
0.7922,
0.9595

None 这些数字重复。下面是我用 C 代码编写的数字食谱版本。关于什么是错误的任何建议?我在这上面花了 3 个小时。该算法为行为良好的随机数集提供了适当的索引。

#define SWAP(a,b) itemp=(a);(a)=(b);(b)=itemp;
#define M 7
#define NSTACK 50

long* IndexSort(unsigned long int vectorLength, double* column)
{

unsigned long i, indxt, ir = vectorLength, itemp,j, k, l = 1;
int jstack = 0, *istack;
long *indx;
double a;

istack = (int*)malloc(NSTACK*sizeof(int));
//initalize output
indx = (long*)malloc(vectorLength*sizeof(long));

for (j = 0; j < vectorLength; j++)
{
    indx[j] = j;

}

// 
while (true)
{
    if (ir - l < M)
    {   
        for (j = l+1; j <= ir;j++)
        {
            indxt = indx[j-1];
            a = column[indxt];
            for (i = j-1; i >= 1; i--)
            {
                if (column[indx[i - 1]] <= a)
                {
                    break;
                }
                indx[i + 1 - 1] = indx[i - 1];
            }
            indx[i + 1 - 1] = indxt;
        }
        if (jstack == 0)
        {
            break;
        }
        ir = istack[jstack--];
        l = istack[jstack--];
    }
    else
    {   

        k = (l + ir) >> 1;

        SWAP(indx[k - 1], indx[l + 1 - 1])

        if (column[indx[l + 1 - 1]] > column[indx[ir - 1]]){
            SWAP(indx[l + 1 - 1], indx[ir - 1])
        }

        if (column[indx[l - 1]] > column[indx[ir - 1]]){
            SWAP(indx[l - 1], indx[ir - 1])
        }

        if (column[indx[l + 1 - 1]] > column[indx[l - 1]]){
            SWAP(indx[l + 1 - 1], indx[l - 1])
        }

        i = l + 1;
        j = ir;
        indxt = indx[l - 1];
        a = column[indxt];
        while (true)
        {
            do i++; while (column[indx[i - 1]] < a);
            do j--; while (column[indx[j - 1]] > a);
            if (j < i) break;
            SWAP(indx[i - 1], indx[j - 1])
        }
        indx[l - 1] = indx[j - 1];
        indx[j - 1] = indxt;
        jstack += 2;
        if (jstack > NSTACK) error("NSTACK too small for IndexSort");

        if (ir - i + 1 > j - 1){
            istack[jstack] = ir;
            istack[jstack - 1] = l;
            ir = j - 1;
        }

        else{
            istack[jstack] = j - 1;
            istack[jstack - 1] = l;
            l = i;
        }


    }
}
free(istack);
return indx;
}

编辑 1: 我将我的#defines 添加到顶部。 回复第一波评论:

编辑 2:

在此先感谢您的帮助,

如果目标是对给定数组的索引进行排序,则使用 C++ 11、std::sort 和 lambda 的以下代码应该可以完成这项工作:

C++11 示例

#include <algorithm>
#include <iostream>
#include <array>

long* IndexSort(unsigned long int vectorLength, double* column)
{
    long *index = new long[vectorLength](); // I really don't recommend this line
    long n = 0;
    std::generate( index, index + vectorLength, [&] { return n++; });
    std::sort(index, index + vectorLength, [&](long v1, long v2) 
             { return column[v1] < column[v2]; });
    return index;
}

using namespace std;

int main()
{  
  std::array<double,20> testData = 
                        {0.8147,0.9058,0.127,0.9134,0.6324,0.0975,0.2785,
                        0.5469,0.9575,0.9649,0.1576,0.9706,0.9572,0.4854,
                        0.8003,0.1419,0.4218,0.9157,0.7922,0.9595};
  long *indices = IndexSort(testData.size(), &testData[0]);
  for (size_t i = 0; i < testData.size(); ++i )
    cout << testData[indices[i]] << "  has an index of " << indices[i] << "\n";
  delete [] indices;
}

实例:http://ideone.com/EsOFKt

请注意,我们所做的只是构建一个排序的索引数组,并使用传递给 lambda 的两个索引比较原始双精度向量。

请注意,没有必要在开头为索引分配内存,因为 std::vector<long> 是执行此操作的常用方法。但是,我试图模仿您返回指向动态分配内存的指针的原始代码要求。


C++ 0x,98 示例

如果您使用的是 C++ 11 之前的编译器,您可以将代码更改为:

#include <algorithm>
#include <iostream>

struct IndexSorter
{
   double *m_array;
   IndexSorter(double *oArray) : m_array(oArray) {}
   bool operator()(long v1, long v2) const { return m_array[v1] < m_array[v2]; }
};

long* IndexSort(unsigned long int vectorLength, double* column)
{
    long *index = new long[vectorLength](); // I really don't recommend this line
    for (unsigned long i = 0; i < vectorLength; ++i) 
        index[i] = i;

    std::sort(index, index + vectorLength, IndexSorter(column));
    return index;
}

using namespace std;

int main()
{  
  double testData[] =   {0.8147,0.9058,0.127,0.9134,0.6324,0.0975,0.2785,
                        0.5469,0.9575,0.9649,0.1576,0.9706,0.9572,0.4854,
                        0.8003,0.1419,0.4218,0.9157,0.7922,0.9595};
  unsigned long testSize = sizeof(testData) / sizeof(testData[0]);
  long *indices = IndexSort(testSize, &testData[0]);
  for (unsigned long i = 0; i < testSize; ++i )
    cout << testData[indices[i]] << "  has an index of " << indices[i] << "\n";
  delete [] indices;
}

实例:http://ideone.com/nhMwSs


C 例子

由于这个问题也被标记为 C,这里是一个 C 实现:

#include <stdlib.h>
#include <stdio.h>

static double *myArray;

int compareIndices(const void* v1, const void *v2)
{
    long val1 = *(long *)v1;
    long val2 = *(long *)v2;
    if (myArray[val1] < myArray[val2])
        return -1;
    else
    if ( myArray[val1] > myArray[val2])
        return 1;
    return 0;
}

long* IndexSort(unsigned long int vectorLength, double* column)
{
    myArray = column;
    long *index = malloc(vectorLength * sizeof(long));
    if ( index )
    {
        unsigned long i;
        for (i = 0; i < vectorLength; ++i)
            index[i] = i;
        qsort(index, vectorLength, sizeof(long), compareIndices);
        return index;
    }
    return 0;
}

int main()
{
    double testData[] = { 0.8147, 0.9058, 0.127, 0.9134, 0.6324, 0.0975, 0.2785,
        0.5469, 0.9575, 0.9649, 0.1576, 0.9706, 0.9572, 0.4854,
        0.8003, 0.1419, 0.4218, 0.9157, 0.7922, 0.9595 };
    unsigned long numItems = sizeof(testData) / sizeof(testData[0]);
    long *indices = IndexSort(numItems, testData);
    if ( indices )
    {
        unsigned long i;
        for (i = 0; i < numItems; ++i)
            printf("%lf has an index of %ld\n", testData[indices[i]], indices[i]);
        free(indices);
    }
}

实例:http://ideone.com/8SZH1m