1. 魔幻C++与埃氏筛:质数计算的极致优化
第一次接触埃氏筛算法是在大学算法课上,当时就被它简洁高效的特性所吸引。多年后当我用C++重新实现这个算法时,发现其中蕴含着许多值得玩味的优化技巧。本文将分享我在实现埃氏筛过程中的实战经验,特别是如何利用C++特性将算法性能推向极限。
质数筛选是计算机科学中的经典问题,埃拉托斯特尼筛法(简称埃氏筛)作为最古老的算法之一,至今仍在实际应用中占有一席之地。不同于试除法逐个判断每个数是否为质数,埃氏筛通过"筛去"合数的方式批量找出质数,时间复杂度为O(n log log n),在处理大规模数据时优势明显。
2. 埃氏筛基础原理与实现
2.1 算法核心思想解析
埃氏筛的基本思路异常简单:从2开始,将每个质数的所有倍数标记为合数。具体步骤可分为:
- 初始化一个布尔数组is_prime[0..n],全部设为true
- 从2开始遍历到√n:
- 如果is_prime[i]为true,则i是质数
- 将i的所有倍数j(从i²开始)标记为false
- 最后仍为true的即为质数
这个算法之所以高效,是因为它避免了重复判断——每个合数只会被其最小质因数筛去一次。例如数字12只会被2筛去(而非3或4),这种"最小质因数筛选"特性是算法效率的关键。
2.2 基础C++实现代码
cpp复制#include <vector>
#include <cmath>
std::vector<int> sieve_of_eratosthenes(int n) {
std::vector<bool> is_prime(n + 1, true);
is_prime[0] = is_prime[1] = false;
for (int i = 2; i * i <= n; ++i) {
if (is_prime[i]) {
for (int j = i * i; j <= n; j += i) {
is_prime[j] = false;
}
}
}
std::vector<int> primes;
for (int i = 2; i <= n; ++i) {
if (is_prime[i]) {
primes.push_back(i);
}
}
return primes;
}
这个基础版本已经能正确工作,但存在几个明显可优化的点:
- 布尔向量std::vector
在内存中的特殊存储方式可能导致性能损失 - 内层循环存在冗余操作
- 没有利用现代CPU的缓存和并行特性
3. 性能优化技巧详解
3.1 内存访问优化
std::vector
cpp复制void sieve_optimized(int n) {
const int size = (n + 7) / 8;
std::vector<unsigned char> is_prime(size, 0xFF);
auto set_composite = [&](int k) {
is_prime[k >> 3] &= ~(1 << (k & 7));
};
auto check_prime = [&](int k) {
return is_prime[k >> 3] & (1 << (k & 7));
};
set_composite(0);
set_composite(1);
for (int i = 2; i * i <= n; ++i) {
if (check_prime(i)) {
for (int j = i * i; j <= n; j += i) {
set_composite(j);
}
}
}
}
这种位操作版本在处理大范围质数时(如n>1e7)可以显著减少内存占用和缓存未命中。
3.2 循环展开与步长优化
观察内层循环,我们可以进行两项关键优化:
- 跳过偶数:除了2,所有偶数都不是质数
- 分块处理:利用CPU缓存局部性
优化后的实现:
cpp复制std::vector<int> sieve_optimized2(int n) {
std::vector<bool> is_prime(n + 1, true);
is_prime[0] = is_prime[1] = false;
// 处理偶数
for (int j = 4; j <= n; j += 2) {
is_prime[j] = false;
}
// 只处理奇数
for (int i = 3; i * i <= n; i += 2) {
if (is_prime[i]) {
// 从i²开始,步长为2i(跳过偶数倍)
for (int j = i * i; j <= n; j += 2 * i) {
is_prime[j] = false;
}
}
}
// 收集结果时也跳过偶数
std::vector<int> primes = {2};
for (int i = 3; i <= n; i += 2) {
if (is_prime[i]) {
primes.push_back(i);
}
}
return primes;
}
这个版本减少了约一半的循环迭代次数,实测在n=1e8时速度提升约40%。
3.3 分段筛法(处理超大范围)
当n非常大(如>1e9)时,内存可能无法容纳整个is_prime数组。此时可以采用分段筛法:
- 先筛出√n以内的所有质数
- 将区间[0,n]分成多个大小为Δ的块
- 对每个块,用预筛的质数筛去该块内的合数
cpp复制void segmented_sieve(int n) {
const int segment_size = 32768; // 32KB块大小
int sqrt_n = static_cast<int>(std::sqrt(n));
// 先筛出小质数
auto small_primes = sieve_of_eratosthenes(sqrt_n);
std::vector<char> sieve(segment_size);
std::vector<int> primes;
for (int low = 0; low <= n; low += segment_size) {
int high = std::min(low + segment_size - 1, n);
std::fill(sieve.begin(), sieve.end(), true);
for (int p : small_primes) {
// 计算第一个p的倍数>=low
int first_multiple = std::max(p * p, ((low + p - 1) / p) * p);
for (int j = first_multiple; j <= high; j += p) {
sieve[j - low] = false;
}
}
// 收集当前块的质数
for (int i = std::max(low, 2); i <= high; ++i) {
if (sieve[i - low]) {
primes.push_back(i);
}
}
}
}
4. 实战性能对比与测试
4.1 不同实现的性能数据
在i7-11800H处理器上测试不同实现处理n=1e8时的表现:
| 实现版本 | 时间(ms) | 内存(MB) |
|---|---|---|
| 基础版本 | 1200 | 95.4 |
| 位操作版 | 850 | 11.9 |
| 步长优化 | 680 | 95.4 |
| 分段筛法 | 920 | 0.03 |
注意:测试时关闭了编译器优化以展示算法差异,实际开启-O3优化后各版本都会显著加快
4.2 编译器优化影响
现代编译器对这类数值计算代码能进行深度优化。关键优化选项:
- -O3:最高级别优化
- -march=native:启用本地CPU特有指令集
- -funroll-loops:循环展开
使用g++编译时添加这些选项后,性能通常能提升3-5倍。例如步长优化版本从680ms降至约150ms。
5. 常见问题与调试技巧
5.1 典型错误排查
- 结果不完整:忘记处理n=1的情况,导致is_prime[1]可能仍为true
- 内存不足:当n很大时,基础版本可能因内存不足崩溃
- 整数溢出:i*i超过INT_MAX导致无限循环
cpp复制// 安全的循环条件写法
for (int i = 2; i <= n / i; ++i) {
// ...
}
5.2 质数验证方法
验证算法正确性的简单方法:
- 检查质数个数:π(1e6)应该等于78498
- 检查特定位置的质数:如第1000个质数应为7919
cpp复制bool validate_results(const std::vector<int>& primes) {
if (primes.size() < 1000) return false;
return primes[0] == 2 && primes[999] == 7919;
}
5.3 多线程优化思路
埃氏筛的内层循环可以并行化,但需要注意:
- 数据竞争:不同线程不能同时修改同一内存位置
- 负载均衡:质数的倍数分布不均匀
一个简单的OpenMP实现:
cpp复制#pragma omp parallel for
for (int i = 2; i <= sqrt_n; ++i) {
if (is_prime[i]) {
for (int j = i * i; j <= n; j += i) {
is_prime[j] = false; // 需要原子操作或分段处理
}
}
}
更安全的做法是预先分配好各线程处理的质数范围,或者使用分段筛法的并行版本。
6. 进阶应用场景
6.1 质数分布统计
埃氏筛的结果可用于分析质数分布规律:
cpp复制void analyze_distribution(const std::vector<int>& primes) {
const int buckets = 10;
std::vector<int> counts(buckets, 0);
int max_p = primes.back();
for (int p : primes) {
int bucket = (p * buckets) / max_p;
counts[bucket]++;
}
// 输出各区间质数密度
for (int i = 0; i < buckets; ++i) {
double start = (i * max_p) / static_cast<double>(buckets);
double end = ((i + 1) * max_p) / static_cast<double>(buckets);
std::cout << start << "-" << end << ": "
<< counts[i] << " primes\n";
}
}
6.2 质数相关数学问题
埃氏筛可用于解决:
- 哥德巴赫猜想验证
- 孪生质数统计
- 质数间隔分析
例如统计孪生质数对:
cpp复制int count_twin_primes(const std::vector<int>& primes) {
int count = 0;
for (size_t i = 0; i < primes.size() - 1; ++i) {
if (primes[i + 1] - primes[i] == 2) {
count++;
}
}
return count;
}
7. 现代C++特性应用
7.1 使用STL算法优化
C++17的并行算法可以简化并行实现:
cpp复制std::vector<int> sieve_parallel(int n) {
std::vector<bool> is_prime(n + 1, true);
is_prime[0] = is_prime[1] = false;
int sqrt_n = static_cast<int>(std::sqrt(n));
std::vector<int> small_primes;
// 先筛出小质数
for (int i = 2; i <= sqrt_n; ++i) {
if (is_prime[i]) {
small_primes.push_back(i);
for (int j = i * i; j <= sqrt_n; j += i) {
is_prime[j] = false;
}
}
}
// 并行筛大区间
std::for_each(std::execution::par, small_primes.begin(), small_primes.end(),
[&](int p) {
for (int j = std::max(p * p, ((sqrt_n + 1 + p - 1) / p) * p);
j <= n; j += p) {
is_prime[j] = false;
}
});
// 收集结果
std::vector<int> primes;
for (int i = 2; i <= n; ++i) {
if (is_prime[i]) {
primes.push_back(i);
}
}
return primes;
}
7.2 编译期质数计算
利用C++20的consteval可以在编译期计算小范围内的质数:
cpp复制consteval auto compile_time_sieve() {
std::array<bool, 100> is_prime{};
is_prime.fill(true);
is_prime[0] = is_prime[1] = false;
for (int i = 2; i * i < is_prime.size(); ++i) {
if (is_prime[i]) {
for (int j = i * i; j < is_prime.size(); j += i) {
is_prime[j] = false;
}
}
}
return is_prime;
}
constexpr auto primes_mask = compile_time_sieve();
static_assert(primes_mask[2] && primes_mask[3] && !primes_mask[4]);
8. 性能极限挑战
8.1 缓存优化策略
现代CPU的缓存行通常为64字节,我们可以调整算法以更好地利用缓存:
- 将筛子分成缓存行大小的块
- 对每个小质数,预计算它对缓存行的影响模式
- 批量应用这些模式
cpp复制void cache_optimized_sieve(int n) {
constexpr int cache_line = 64;
const int size = (n + cache_line - 1) / cache_line;
std::vector<uint64_t> sieve(size, ~0ULL);
auto set_composite = [&](int k) {
sieve[k / cache_line] &= ~(1ULL << (k % cache_line));
};
set_composite(0);
set_composite(1);
for (int i = 2; i * i <= n; ++i) {
if (sieve[i / cache_line] & (1ULL << (i % cache_line))) {
for (int j = i * i; j <= n; j += i) {
set_composite(j);
}
}
}
}
8.2 SIMD指令加速
使用AVX2指令集可以进一步加速筛法:
cpp复制#include <immintrin.h>
void simd_sieve(int n) {
constexpr int vec_size = 256 / 8; // AVX2寄存器大小
const int size = (n + vec_size - 1) / vec_size;
std::vector<__m256i> sieve(size, _mm256_set1_epi8(0xFF));
// 设置合数的辅助函数
auto set_composite = [&](int k) {
int idx = k / vec_size;
int bit = k % vec_size;
__m256i mask = _mm256_set1_epi8(~(1 << bit));
sieve[idx] = _mm256_and_si256(sieve[idx], mask);
};
set_composite(0);
set_composite(1);
for (int i = 2; i * i <= n; ++i) {
if (_mm256_movemask_epi8(sieve[i / vec_size]) & (1 << (i % vec_size))) {
for (int j = i * i; j <= n; j += i) {
set_composite(j);
}
}
}
}
这种优化在支持AVX2的CPU上可以将性能提升2-3倍,但代码可移植性会降低。
