分段双调排序

今天在"https://blog.csdn.net/u014226072/article/details/56840243"上看到一个关于分段双调排序的问题. 抄录如下:

问题说明：

***********************

给出分成m段的n个浮点数，输入数据已按段号有序，但每段内部无序。用C/C++ 编写一个分段双调排序(Bitonic sort)函数，对每一段内部的浮点数进行排序，但不要改变段间的位置。

接口方式：

***********************

void segmentedBitonicSort(float* data, int*seg_id, int* seg_start, int n, int m);

输入数据中，data包含需要分段排序的n个float值，seg_id给出data中n个元素各 自所在的 段编号。seg_start共有m+1个元素，前m个分别给 出0..m-1共m个段的起 始位置，seg_start[m]保证等于n。

seg_id中的元素保证单调不下降，即对任意的i<j，seg_id[i]<=seg_id[j]。 seg_id所有元 素均在0到m-1范围内。

输出结果覆盖data，保证每一段内排序，但不改变段间元素的顺序。

注意：

***********************

1、必须使用双调排序算法进行排序。

2、可以直接使用从网上下载的双调排序代码，但须注明出处。

样例输入：

***********************

float data[5]={0.8, 0.2, 0.4, 0.6, 0.5};

int seg_id[5]={0,   0,  1,   1,   1}

int seg_start[3]={0,2,5};

int n=5;

int m=2;

样例输出：

***********************

float data[5]={0.2, 0.8, 0.4, 0.5, 0.6};

加分挑战（非必需）：

***********************

1、不递归：segmentedBitonicSort函数及其所调用的任何其他函数都不得直接或间接地进行递归。

2、不调用函数：segmentedBitonicSort不调用除标准库函数外的任何其他函数。

3、内存高效：segmentedBitonicSort及其所调用的任何其他函数都不得进行动态内存分配，包括malloc、new和静态定义的STL容器。

4、可并行：segmentedBitonicSort涉及到的所有时间复杂度O(n)以上的代码都写 在for循 环中，而且每个这样的for循环内部的循环顺序可 以任意改变，不影响程 序结果。注：自己测试时可以用rand()决定循环顺序。

5、不需内存：segmentedBitonicSort不调用任何函数（包括C/C++标准库函数）， 不使用全局变量，所有局部变量都是int、float或指针类 型，C++程序不使用new关键字。

6、绝对鲁棒：在输入数据中包含NaN时（例如sqrt(-1.f)），保证除NaN以外 的数 据正确排序，NaN的个数保持不变。

你的程序每满足以上的一个条件都可以获得额外的加分。

为此, 我们先了解一下双调排序. 下图来自于:

https://www.cs.rutgers.edu/~venugopa/parallel_summer2012/bitonic_overview.html

双调排序是一种基于比较的可并行的排序算法. 

第一点, 如果一个序列前半部分和后半部分都单调, 那么称它为一个双调序列. (图中红色部分是前半部分, 递增. 蓝色部分递减)
第二点, 对一个双调序列, 依次比较两个单调部分的对应位置的元素, 视比较结果交换位置. 可以得到两个序列(s1和s2). 
第三点, 这两个序列都是双调序列, 并且前一个双调序列的每个元素都小于后一个双调序列. 
第四点, 此时, 对于整个数组来说, 前部分小于(同理也可以是大于)后半部分. 那么我们可以递归地使用这种拆分双调序列的方法, 直到序列长度为1, 这时, 整个数组将是有序的. 
------
但是我们需要排序的数组并不一定一开始就是双调序列的, 所以要应用这个排序方法还需要换个思考方式. 
- 注意, 长度为1的序列可以看作是单调的, 那么, 长度为2的序列一定是双调序列. 
- 双调序列可以排序, 所以我们先给这些长度为2的双调序列排序. 

- 排序的时候, 如果排序的方式我们选择递增和递减交替. 排序结束时, 我们将得到长度为4的双调序列. 

好了, 问题已经解决了. 我们终将得到长度为n的双调序列, 也就是把原数组变成了双调序列. 当然, 值得一提的是, 这种方法只能给长度为2的整数次幂的数组排序, 对于一个任意的数组, 需要用一个较大的值填充至长度为2^n之后才能用这种方法. 

不严格的时间复杂度分析: 容易看出, 给一个双调序列排序的时间复杂度是nlogn(一共logn轮比较), 而为了使原数组变成双调序列, 我们要进行logn轮, 每轮n/len个长度为len的排序, 其中在第i轮时, len=logi. 所以有一个很明显的上界O(n*logn*logn). 但是这是串行分析的结果. 事实上在并行条件下可以大大加速. 


具体到这道题, 一个很容易想到的思路是把seg_id当作第一关键字, data当作第二关键字, 把整个数组补齐到2^n后整体排序. 据此, 可以很容易写出满足不递归, 不调用函数, 可并行, 绝对鲁棒的代码. 目前我还想不到如何不动态分配内存. 关于绝对鲁棒, 我是采用手动设置nan大于任何数字来避免因nan而出现排序错乱的问题的. 


但是这里还可以再有一点改进, 我们可以分别给每一段排序. 这样可以把时间复杂度从O(n*log^2n)降为近似O(m * n/m * log^2(n/m)). 所需额外空间也更少. 

最终代码如下:

#include <cstdio>
#include <stdlib.h>
#include <time.h>
#include <assert.h>
#include <cmath>
#include <iostream>

#define SEGN 50
#define SEGLEN 50

#define nan (std::sqrt(-1.f))

//#define debug

void seg3(float* data, int* seg_id, int* seg_start, int n, int m) 
{
    //#ifdef debug
    printf("seg3:\nn = %d, m = %d\n", n, m);
    for (int i = 0; i < n; i++) std::cout << data[i] << ", ";
    printf("\n");
    for (int i = 0; i < n; i++) printf("%d, ", seg_id[i]);
    printf("\n");
    for (int i = 0; i <= m; i++) printf("%d, ", seg_start[i]);
    printf("\n");
    //#endif
    
    for (int seg = 0; seg < m; seg++) {
        int curlen = seg_start[seg + 1] - seg_start[seg];
        int pow = 0;
        int tempn = curlen;
        while (tempn) {
            tempn >>= 1;
            pow++;
        }
        tempn = (1 << pow);
        if (tempn == (curlen << 1)) {
            tempn >>= 1;
        }
        int endid = seg_start[seg] + tempn;
        float* datab;
        int* segb;
        if (endid > n) {
            datab = (float*)malloc(sizeof(float) * (endid - n));
            segb = (int*)malloc(sizeof(int) * (endid - n));
            for (int i = 0; i < endid - n; i++) {
                segb[i] = m;
            }
        }
        else {
            datab = data + seg_start[seg + 1];
            segb = seg_id + seg_start[seg + 1];
        }
        // sort begin
        for (int s = 2; s <= tempn; s *= 2) {
            for (int itr = 0; itr < tempn; itr += s * 2) {
                // mergeup(arr+i, s)
                // float* arr = data + i;
                int offset = itr + seg_start[seg];
                int len = s;
                // begin
                int step = len / 2, i, j, k, temp;
                float tempf;
                while (step > 0) {
                    for (i = 0; i < len; i += step * 2) {
                        for (j = i, k = 0; k < step; j++, k++) {
                            // if arr[j] > arr[j + step], swap
                            int* seg_add_j = offset + j < n ? seg_id + offset + j : segb + offset + j - n;
                            int* seg_add_js = offset + j + step < n ? seg_id + offset + j + step : segb + offset + j + step - n;
                            float* data_add_j = offset + j < n ? data + offset + j : datab + offset + j - n;
                            float* data_add_js = offset + j + step < n ? data + offset + j + step : datab + offset + j + step - n;
                            bool swap = false;
                            
                            if (*seg_add_j > *seg_add_js) {
                                swap = true;
                            }
                            else if (*seg_add_j == *seg_add_js) {  
                                if (*data_add_j > *data_add_js) {
                                    swap = true;
                                }
                                else if (!(*data_add_j > 0) && !(*data_add_j <= 0)) {
                                    // *data_add_j = nan
                                    swap = true;
                                }
                            }
                            
                            if (swap) {
                                temp = *seg_add_j;
                                *seg_add_j = *seg_add_js;
                                *seg_add_js = temp;
                                tempf = *data_add_j;
                                *data_add_j = *data_add_js;
                                *data_add_js = tempf;
                            }
                        }
                    }
                    step >>= 1;
                }
                //mergedown(arr+i+s,s)
                //arr = data + i + s;
                if (s == tempn) continue;
                offset = itr + len + seg_start[seg];
                step = len / 2;
                // begin
                while (step > 0) {
                    for (i = 0; i < len; i += step * 2) {
                        for (j = i, k = 0; k < step; j++, k++) {
                            // if arr[j] < arr[j + step], swap
                            int* seg_add_j = offset + j < n ? seg_id + offset + j : segb + offset + j - n;
                            int* seg_add_js = offset + j + step < n ? seg_id + offset + j + step : segb + offset + j + step - n;
                            float* data_add_j = offset + j < n ? data + offset + j : datab + offset + j - n;
                            float* data_add_js = offset + j + step < n ? data + offset + j + step : datab + offset + j + step - n;
                            
                            bool swap = false;
                            if (*seg_add_j < *seg_add_js) {
                                swap = true;
                            }
                            else if (*seg_add_j == *seg_add_js) {  
                                if (*data_add_j < *data_add_js) {
                                    swap = true;
                                }
                                else if (!(*data_add_js > 0) && !(*data_add_js <= 0)) {
                                    // *data_add_js = nan
                                    swap = true;
                                }
                            }
                            
                            if (swap) {
                                temp = *seg_add_j;
                                *seg_add_j = *seg_add_js;
                                *seg_add_js = temp;
                                tempf = *data_add_j;
                                *data_add_j = *data_add_js;
                                *data_add_js = tempf;
                            }
                        }
                    }
                    step >>= 1;
                }
            }
        }
        if (endid > n) {
            free(datab);
            free(segb); 
        }
    }
}

int main() {
    #ifdef debug
    float data[] = {1.08344, 2.6747, nan, 1.94815, 1.11102, 0.289917, 0.64205,
    6.59573, 0.00369898, 12.7636, 1.56606, 0.0572074, 0.380153, nan, 2.96348,
    3.89681, nan, 2.27465, 1.71466, 3.69074, 0.351484};
  int seg_id[] = {0, 0, 0, 0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4};
  int seg_start[] = {0, 5, 7, 11, 16, 21};
  int n = 21;
  int m = 5;
  seg3(data, seg_id, seg_start, n, m);
        for (int j = 0; j < n; j++) {
            std::cout << data[j] << ' ';
        }
        std::cout << std::endl;
        //
        bool wrong = false;
        for (int j = 1; j < n; j++) {
            if (seg_id[j] < seg_id[j - 1]) {
                wrong = true;
                break;
            }
            
            if (!(seg_id[j] > seg_id[j - 1] || (!(data[j] < data[j - 1])))) {
                std::cout << "j :" << j << std::endl;
                std::cout << data[j] << ' ' << data[j-1] << std::endl;
                std::cout << seg_id[j] << ' ' << seg_id[j-1] << std::endl;
                wrong = true;
                break;
            }
        }
        if (wrong) std::cerr << "wrong^" << std::endl;
    #else
    srand((unsigned)time(NULL));
    // test 
    float* data;
    int* seg_id;
    int* seg_start;
    int n;
    int m;
    for (int i = 0; i < 10; i++) {
        m = rand() % SEGN + 5;
        seg_start = (int*)malloc(sizeof(int) * (m + 1));
        seg_start[0] = 0;
        for (int j = 1; j <= m; j++) {
            seg_start[j] = rand() % SEGLEN + 1 + seg_start[j-1];
        }
        n = seg_start[m];
        seg_id = (int*)malloc(sizeof(int)*n);
        data = (float*)malloc(sizeof(float)*n);
        for (int j = 1; j <= m; j++) {
            for (int k = seg_start[j - 1]; k < seg_start[j]; k++) {
                seg_id[k] = j - 1;
            }
        }
        for (int j = 0; j < n; j++) {
            float tempf = 1.0 * rand() / rand();
            if (tempf < 0.6 && tempf > 0.4) {
                data[j] = std::sqrt(-1.f);
                //data[j] = 0.5;
            }
            else {
                data[j] = tempf;
            }
        }
        // 
        seg3(data, seg_id, seg_start, n, m);
        //
        for (int j = 1; j < n; j++) {
            assert(seg_id[j] >= seg_id[j - 1]);
            assert(seg_id[j] > seg_id[j - 1] || (!(data[j] < data[j - 1])));
        }
        free(seg_start);
        free(seg_id);
        free(data);
    }
    #endif
    return 0;
}

后来我又找到两篇文章, 实现了不用动态分配内存的任意长度数组双调排序, 以后有时间试试把它改为不调用函数的形式. 
https://blog.csdn.net/ljiabin/article/details/8631374
https://blog.csdn.net/lqmiku/article/details/78834178


Ref:
https://blog.csdn.net/u014226072/article/details/56840243
https://www.cs.rutgers.edu/~venugopa/parallel_summer2012/bitonic_overview.html