有一些方法可以生成已经排序的样本，但我认为生成部分排序的样本可能会更好。
将输出范围划分为 k 个等宽的桶。每个桶中的样本数将具有等概率的多项分布。对多项分布进行采样的慢速方法是在 [0, k) 中生成 n 个整数。一种更有效的方法是抽取 k 个泊松样本，其速率为 n/k，条件是它们的总和不超过 n，然后使用慢速方式添加另一个 n-sum 样本。完美地对泊松分布进行采样是很棘手的，但是当 n/k 非常大时（就像这里所说的那样），泊松分布可以通过对具有均值和方差 n/k 的正态分布进行四舍五入来很好地逼近。如果这是不可接受的，那么慢速方法确实可以很好地并行化。
给定存储桶计数，计算前缀总和以找到存储桶边界。对于并行的每个桶，在桶范围内生成给定数量的样本并对其进行排序。如果我们选择好 n/k，每个桶几乎肯定会适合 L1 缓存。对于 n = 1e9，我想我会尝试 k = 1e5 或 k = 1e6。
这是一个顺序实现。有点粗糙，因为我们确实需要避免对封闭的桶边界进行 2 倍过采样，但我会把它留给你。我对 OMP 不熟悉，但我认为您可以通过在 SortedUniformSamples 末尾的 for 循环中添加 pragma 来获得相当不错的并行实现。
#include <algorithm>
#include <cmath>
#include <iostream>
#include <numeric>
#include <random>
#include <span>
#include <vector>
template <typename Dist, typename Gen>
void SortedSamples(std::span<double> samples, Dist dist, Gen& gen) {
for (double& sample : samples) {
sample = dist(gen);
}
std::sort(samples.begin(), samples.end());
}
template <typename Gen>
void ApproxMultinomialSample(std::span<std::size_t> samples, std::size_t n,
Gen& gen) {
double lambda = static_cast<double>(n) / samples.size();
std::normal_distribution<double> approx_poisson{lambda, std::sqrt(lambda)};
std::size_t sum;
do {
for (std::size_t& sample : samples) {
sample = std::lrint(approx_poisson(gen));
}
sum = std::accumulate(samples.begin(), samples.end(), std::size_t{0});
} while (sum > n);
std::uniform_int_distribution<std::size_t> uniform{0, samples.size() - 1};
for (; sum < n; sum++) {
samples[uniform(gen)]++;
}
}
template <typename Gen>
void SortedUniformSamples(std::span<double> samples, Gen& gen) {
static constexpr std::size_t kTargetBucketSize = 1024;
if (samples.size() < kTargetBucketSize) {
SortedSamples(samples, std::uniform_real_distribution<double>{0, 1}, gen);
return;
}
std::size_t num_buckets = samples.size() / kTargetBucketSize;
std::vector<std::size_t> bucket_counts(num_buckets);
ApproxMultinomialSample(bucket_counts, samples.size(), gen);
std::vector<std::size_t> prefix_sums(num_buckets + 1);
std::partial_sum(bucket_counts.begin(), bucket_counts.end(),
++prefix_sums.begin());
for (std::size_t i = 0; i < num_buckets; i++) {
SortedSamples(std::span<double>{&samples[prefix_sums[i]],
&samples[prefix_sums[i + 1]]},
std::uniform_real_distribution<double>{
static_cast<double>(i) / num_buckets,
static_cast<double>(i + 1) / num_buckets},
gen);
}
}
int main() {
std::vector<double> samples(100000000);
std::default_random_engine gen;
SortedUniformSamples(samples, gen);
if (std::is_sorted(samples.begin(), samples.end())) {
std::cout << "sorted\n";
}
}
如果您的标准库有 poisson_distribution 的高质量实现，您也可以这样做：
template <typename Gen>
void MultinomialSample(std::span<std::size_t> samples, std::size_t n,
Gen& gen) {
double lambda = static_cast<double>(n) / samples.size();
std::poisson_distribution<std::size_t> poisson{lambda};
std::size_t sum;
do {
for (std::size_t& sample : samples) {
sample = poisson(gen);
}
sum = std::accumulate(samples.begin(), samples.end(), std::size_t{0});
} while (sum > n);
std::uniform_int_distribution<std::size_t> uniform{0, samples.size() - 1};
for (; sum < n; sum++) {
samples[uniform(gen)]++;
}
}