筛子排序（SieveSort) - 5

为什么不只是16，而是还要256或者4096一组？

因为要排序的元素太少的话，启动线程不划算。

要排序的数量越大，用omp的多线程进行排序就越快。

总要有一个标准去分辨快慢，我们可以选择std::sort，这个函数已经被高度优化，可以适合于各种数据类型和长度。

有了16和256这两个基本排序能力。下面的就可以考虑递归了：如果是小于等于16的，用16来排；如果大于16小于等于256的，用256来排。大于256，小于等于4K的，最多分成16分，每一份尽可能达到256，最后一份小于256的，是多少就是多少，也可能等于256。大于4K，小于等于64K的，也最多分成16分，每一份尽可能达到4K，最后一份小于4K的，是多少就多少。然后是64K到1M的，1M到16M的，16M到256M的，256M到4G的……都是这个原则。这显然是一种递归分治的做法。这里有两个难点，一个是这种划分方式怎么实现才能尽可能用最少的代码处理最多的情况，另一个是，这种归并方式超过256之后就无法实现原地排序（不难理解，就不解释了），所以一定会要求更多内存。而栈上的内存肯定是不够的，就必须从堆上分配。然而分配内存又是一个耗时的过程，所以尽可能最先分配完了，后来使用即可，具体在过程中如何处理内存中的数据，也是必须仔细考虑的。

bool get_config(size_t n, size_t& loops, size_t& stride, size_t& reminder, __mmask16& mask, bool& flip, int min_bits = 8) {
	if (n < ((1ULL) << min_bits)) return false;
	int top_bits = get_top_bit_index(n);
	int nb_count = (int)__popcnt64(n);
	int all_bits = nb_count == 1 ? (top_bits - 1) : (top_bits & ~0x3) == top_bits ? top_bits : (top_bits + ((n - ((n >> 4) << 4)) != 0ULL));
	int max_bits = ((all_bits >> 2) + ((all_bits & 0x3) != 0)) << 2;

	stride = (1ULL) << (max_bits - 4);
	reminder = n & (~((~0ULL) << (max_bits - 4)));
	loops = (n - reminder) / stride + (reminder > 0);
	mask = ~((~0U) << (loops));

	flip = (max_bits > 12) && (((max_bits >> 2) & 1) == 0);

	return true;
}

这个get_config函数就给出了分块处理的具体原则：给定数据总量为n，n需要至少大于256（256以内的用别的方法处理），计算n的二进制位数，并把这个位数每四位（对应16倍）一节划分，对于具体的数值，生成其循环的次数（loops），对应的掩码（mask），一个分段的长度（stride）以及最后一段可能的余量（reminder：也可以是0）。还有一个翻转条件（flip）用于处理源和宿之间的关系。把这个配置计算好就可以进行具体的排序处理了。

bool sieve_sort_core(uint32_t* a, size_t n, uint32_t* result, int depth);
bool sieve_sort_omp(uint32_t* a, size_t n, uint32_t* result, int depth) {
	size_t loops = 0, stride = 0, reminder = 0;
	__mmask16 mask = 0;
	bool flip = false;
	if (!get_config(n, loops, stride, reminder, mask, flip)) return false;

	if (depth > 0) {
#pragma omp parallel for
		for (int i = 0; i < loops; i++) {
			sieve_sort_core(a + i * stride,
				(i == loops - 1 && reminder > 0) ? reminder : stride,
				result + i * stride,
				depth - 1);
		}
	}
	else {
		for (int i = 0; i < loops; i++) {
			sieve_sort_core(a + i * stride,
				(i == loops - 1 && reminder > 0) ? reminder : stride,
				result + i * stride,
				depth - 1);
		}
	}
	return sieve_collect(n, loops, stride, reminder, mask, flip, result, a);
}

具体排序过程，分成两种情况，omp并行一种作为可选项，另一种就是普通的顺序执行。

代码是一样的，但执行条件根据depth来判别。如果有一种方法把两种情况都一样处理就好了，就不需要同样的代码写两次了，这是很忌讳的。